Obc tracer fix (#507)

* +Fix for issue #506, tracer OBC bug

 - it only happens in the advection for certain flow directions,
   inflow from OBC plus along-boundary flow.

* Tracer OBCs need more of an h halo update.

- This one should finally fix the processor count issues with OBCs.

* Correct the "if" statement.

* +Adding more halo points to an exchange

 - This will change answers if you start with a non-zero velocity.
   You need three halo points (or maybe cont_stencil) for the
   continuity solver.
 - Also trying to put in some initial DEBUG_OBC code.

* Fix some DEBUG_OBC logic

* Writing to temporary arrays for tres_xy

 - Way to trick some compiler.

* Another fix for DEBUG_OBC

* Fixing the whalo troubles for grids that are 2 wide/long.

* Exchange all the h_new points.

 - without this, because we're remapping all the tres points, it
 dies in debug mode on bad h_new values.

* Trying a different exchange

 - as it was, it was passing my tests, failing the auto-tests.

* Fixing the DEBUG_OBC logging

* Putting the logging statement back.

 - Making an error more verbose too.

src/core/MOM.F90

@@ -1627,9 +1627,11 @@ subroutine step_MOM_thermo(CS, G, GV, US, u, v, h, tv, fluxes, dtdia, &
         if (CS%split) &
           call remap_dyn_split_RK2_aux_vars(G, GV, CS%dyn_split_RK2_CSp, h, h_new, CS%ALE_CSp, CS%OBC, dzRegrid)
 
-        if (associated(CS%OBC)) &
-          call pass_var(h_new, G%Domain)
+        if (associated(CS%OBC)) then
+          call pass_var(h, G%Domain, complete=.false.)
+          call pass_var(h_new, G%Domain, complete=.true.)
           call remap_OBC_fields(G, GV, h, h_new, CS%OBC, PCM_cell=PCM_cell)
+        endif
 
         call remap_vertvisc_aux_vars(G, GV, CS%visc, h, h_new, CS%ALE_CSp, CS%OBC)
         if (associated(CS%visc%Kv_shear)) &
@@ -3016,10 +3018,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, &
     call cpu_clock_begin(id_clock_pass_init)
     call create_group_pass(tmp_pass_uv_T_S_h, CS%u, CS%v, G%Domain)
     if (use_temperature) then
-      call create_group_pass(tmp_pass_uv_T_S_h, CS%tv%T, G%Domain, halo=1)
-      call create_group_pass(tmp_pass_uv_T_S_h, CS%tv%S, G%Domain, halo=1)
+      call create_group_pass(tmp_pass_uv_T_S_h, CS%tv%T, G%Domain)
+      call create_group_pass(tmp_pass_uv_T_S_h, CS%tv%S, G%Domain)
     endif
-    call create_group_pass(tmp_pass_uv_T_S_h, CS%h, G%Domain, halo=1)
+    call create_group_pass(tmp_pass_uv_T_S_h, CS%h, G%Domain)
     call do_group_pass(tmp_pass_uv_T_S_h, G%Domain)
     call cpu_clock_end(id_clock_pass_init)
 

src/core/MOM_open_boundary.F90

@@ -7,33 +7,33 @@ module MOM_open_boundary
 use MOM_array_transform,      only : allocate_rotated_array
 use MOM_coms,                 only : sum_across_PEs, Set_PElist, Get_PElist, PE_here, num_PEs
 use MOM_cpu_clock,            only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE
+use MOM_debugging,            only : hchksum, uvchksum
 use MOM_diag_mediator,        only : diag_ctrl, time_type
 use MOM_domains,              only : pass_var, pass_vector
 use MOM_domains,              only : create_group_pass, do_group_pass, group_pass_type
 use MOM_domains,              only : To_All, EAST_FACE, NORTH_FACE, SCALAR_PAIR, CGRID_NE, CORNER
+use MOM_dyn_horgrid,          only : dyn_horgrid_type
 use MOM_error_handler,        only : MOM_mesg, MOM_error, FATAL, WARNING, NOTE, is_root_pe
 use MOM_file_parser,          only : get_param, log_version, param_file_type, log_param
 use MOM_grid,                 only : ocean_grid_type, hor_index_type
-use MOM_dyn_horgrid,          only : dyn_horgrid_type
 use MOM_interface_heights,    only : thickness_to_dz
+use MOM_interpolate,          only : init_external_field, time_interp_external, time_interp_external_init
+use MOM_interpolate,          only : external_field
 use MOM_io,                   only : slasher, field_size, SINGLE_FILE
 use MOM_io,                   only : vardesc, query_vardesc, var_desc
+use MOM_obsolete_params,      only : obsolete_logical, obsolete_int, obsolete_real, obsolete_char
+use MOM_regridding,           only : regridding_CS
+use MOM_remapping,            only : remappingSchemesDoc, remappingDefaultScheme, remapping_CS
+use MOM_remapping,            only : initialize_remapping, remapping_core_h, end_remapping
 use MOM_restart,              only : register_restart_field, register_restart_pair
 use MOM_restart,              only : query_initialized, MOM_restart_CS
-use MOM_obsolete_params,      only : obsolete_logical, obsolete_int, obsolete_real, obsolete_char
-use MOM_string_functions,     only : extract_word, remove_spaces, uppercase
+use MOM_string_functions,     only : extract_word, remove_spaces, uppercase, lowercase
 use MOM_tidal_forcing,        only : astro_longitudes, astro_longitudes_init, eq_phase, nodal_fu, tidal_frequency
 use MOM_time_manager,         only : set_date, time_type, time_type_to_real, operator(-)
 use MOM_tracer_registry,      only : tracer_type, tracer_registry_type, tracer_name_lookup
-use MOM_interpolate,          only : init_external_field, time_interp_external, time_interp_external_init
-use MOM_interpolate,          only : external_field
-use MOM_remapping,            only : remappingSchemesDoc, remappingDefaultScheme, remapping_CS
-use MOM_remapping,            only : initialize_remapping, remapping_core_h, end_remapping
-use MOM_regridding,           only : regridding_CS
 use MOM_unit_scaling,         only : unit_scale_type
 use MOM_variables,            only : thermo_var_ptrs
 use MOM_verticalGrid,         only : verticalGrid_type
-use MOM_string_functions,     only : lowercase
 
 implicit none ; private
 
@@ -528,13 +528,16 @@ subroutine open_boundary_config(G, US, param_file, OBC)
       OBC%add_tide_constituents = .false.
     endif
 
-    call get_param(param_file, mdl, "DEBUG", OBC%debug, default=.false.)
-    call get_param(param_file, mdl, "DEBUG_OBC", debug_OBC, default=.false.)
-    if (debug_OBC .or. OBC%debug) &
+    call get_param(param_file, mdl, "DEBUG", debug_OBC, default=.false.)
+    call get_param(param_file, mdl, "DEBUG_OBC", debug_OBC, default=debug_OBC, &
+                   do_not_log=.not.debug_OBC)
+    if (debug_OBC) then
       call log_param(param_file, mdl, "DEBUG_OBC", debug_OBC, &
                  "If true, do additional calls to help debug the performance "//&
                  "of the open boundary condition code.", default=.false., &
                  debuggingParam=.true.)
+      OBC%debug = debug_OBC
+    endif
 
     call get_param(param_file, mdl, "OBC_SILLY_THICK", OBC%silly_h, &
                  "A silly value of thicknesses used outside of open boundary "//&
@@ -854,6 +857,8 @@ subroutine initialize_segment_data(G, GV, US, OBC, PF)
 !       if (siz(4) == 1) segment%values_needed = .false.
         if (segment%on_pe) then
           if (OBC%brushcutter_mode .and. (modulo(siz(1),2) == 0 .or. modulo(siz(2),2) == 0)) then
+            write(mesg,'("Brushcutter mode sizes ", I6, I6))') siz(1), siz(2)
+            call MOM_error(WARNING, mesg // " " // trim(filename) // " " // trim(fieldname))
             call MOM_error(FATAL,'segment data are not on the supergrid')
           endif
           siz2(1)=1
@@ -2224,6 +2229,8 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US,
   type(OBC_segment_type), pointer :: segment => NULL()
   integer :: i, j, k, is, ie, js, je, m, nz, n
   integer :: is_obc, ie_obc, js_obc, je_obc
+  logical :: sym
+  character(len=3) :: var_num
 
   is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
 
@@ -3298,6 +3305,29 @@ subroutine radiation_open_bdry_conds(OBC, u_new, u_old, v_new, v_old, G, GV, US,
 
   call pass_vector(u_new, v_new, G%Domain, clock=id_clock_pass)
 
+  if (OBC%debug) then
+    sym = G%Domain%symmetric
+    if (OBC%radiation_BCs_exist_globally) then
+      call uvchksum("radiation_OBCs: OBC%r[xy]_normal", OBC%rx_normal, OBC%ry_normal, G%HI, &
+                  haloshift=0, symmetric=sym, scale=1.0)
+    endif
+    if (OBC%oblique_BCs_exist_globally) then
+      call uvchksum("radiation_OBCs: OBC%r[xy]_oblique_[uv]", OBC%rx_oblique_u, OBC%ry_oblique_v, G%HI, &
+                  haloshift=0, symmetric=sym, scale=1.0/US%L_T_to_m_s**2)
+      call uvchksum("radiation_OBCs: OBC%r[yx]_oblique_[uv]", OBC%ry_oblique_u, OBC%rx_oblique_v, G%HI, &
+                  haloshift=0, symmetric=sym, scale=1.0/US%L_T_to_m_s**2)
+      call uvchksum("radiation_OBCs: OBC%cff_normal_[uv]", OBC%cff_normal_u, OBC%cff_normal_v, G%HI, &
+                  haloshift=0, symmetric=sym, scale=1.0/US%L_T_to_m_s**2)
+    endif
+    if (OBC%ntr == 0) return
+    if (.not. allocated (OBC%tres_x) .or. .not. allocated (OBC%tres_y)) return
+    do m=1,OBC%ntr
+      write(var_num,'(I3.3)') m
+      call uvchksum("radiation_OBCs: OBC%tres_[xy]_"//var_num, OBC%tres_x(:,:,:,m), OBC%tres_y(:,:,:,m), G%HI, &
+                    haloshift=0, symmetric=sym, scale=1.0)
+    enddo
+  endif
+
 end subroutine radiation_open_bdry_conds
 
 !> Applies OBC values stored in segments to 3d u,v fields
@@ -5638,9 +5668,6 @@ subroutine remap_OBC_fields(G, GV, h_old, h_new, OBC, PCM_cell)
                      To_All+Scalar_Pair)
   if (OBC%oblique_BCs_exist_globally) then
     call do_group_pass(OBC%pass_oblique, G%Domain)
-!   call pass_vector(OBC%rx_oblique_u, OBC%ry_oblique_v, G%Domain, To_All+Scalar_Pair)
-!   call pass_vector(OBC%ry_oblique_u, OBC%rx_oblique_v, G%Domain, To_All+Scalar_Pair)
-!   call pass_vector(OBC%cff_normal_u, OBC%cff_normal_v, G%Domain, To_All+Scalar_Pair)
   endif
 
 end subroutine remap_OBC_fields

src/framework/MOM_restart.F90

@@ -660,7 +660,7 @@ end subroutine register_restart_field_0d
 
 
 !> query_initialized_name determines whether a named field has been successfully
-!! read from a restart file or has otherwise been recored as being initialzed.
+!! read from a restart file or has otherwise been recorded as being initialized.
 function query_initialized_name(name, CS) result(query_initialized)
   character(len=*),     intent(in) :: name  !< The name of the field that is being queried
   type(MOM_restart_CS), intent(in) :: CS    !< MOM restart control struct
@@ -1271,7 +1271,7 @@ subroutine only_read_restart_pair_3d(a_ptr, b_ptr, a_name, b_name, G, CS, &
 
 end subroutine only_read_restart_pair_3d
 
-!> Return an indicationof whether the named variable is the restart files, and provie the full path
+!> Return an indication of whether the named variable is in the restart files, and provide the full path
 !! to the restart file in which a variable is found.
 function find_var_in_restart_files(varname, G, CS, file_path, filename, directory, is_global) result (found)
   character(len=*),                intent(in)    :: varname   !< The variable name to be used in the restart file

src/parameterizations/lateral/MOM_hor_visc.F90

@@ -414,12 +414,13 @@ subroutine horizontal_viscosity(u, v, h, diffu, diffv, MEKE, VarMix, G, GV, US,
   Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
 
   h_neglect  = GV%H_subroundoff
-  h_neglect3 = h_neglect**3
+  !h_neglect3 = h_neglect**3
+  h_neglect3 = h_neglect*h_neglect*h_neglect
   inv_PI3 = 1.0/((4.0*atan(1.0))**3)
   inv_PI2 = 1.0/((4.0*atan(1.0))**2)
   inv_PI6 = inv_PI3 * inv_PI3
 
-  m_leithy(:,:) = 0. ! Initialize
+  m_leithy(:,:) = 0.0 ! Initialize
 
   if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then
     apply_OBC = OBC%Flather_u_BCs_exist_globally .or. OBC%Flather_v_BCs_exist_globally