Fix bugs introduced during merge

src/core/MOM.F90

@@ -241,7 +241,6 @@ module MOM
                                      !! number of dynamics steps in nstep_tot
   logical :: debug                   !< If true, write verbose checksums for debugging purposes.
   integer :: ntrunc                  !< number u,v truncations since last call to write_energy
-
   integer :: cont_stencil            !< The stencil for thickness from the continuity solver.
   ! These elements are used to control the dynamics updates.
   logical :: do_dynamics             !< If false, does not call step_MOM_dyn_*. This is an
@@ -1186,7 +1185,7 @@ subroutine step_MOM_tracer_dyn(CS, G, GV, US, h, Time_local)
   call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo)
 
   call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics)
-  call post_transport_diagnostics(G, GV, CS%uhtr, CS%vhtr, h, CS%transport_IDs, &
+  call post_transport_diagnostics(G, GV, US, CS%uhtr, CS%vhtr, h, CS%transport_IDs, &
            CS%diag_pre_dyn, CS%diag, CS%t_dyn_rel_adv, CS%diag_to_Z_CSp, CS%tracer_reg)
   ! Rebuild the remap grids now that we've posted the fields which rely on thicknesses
   ! from before the dynamics calls
@@ -1643,6 +1642,7 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
   integer :: IsdB, IedB, JsdB, JedB
   real    :: dtbt        ! The barotropic timestep [s]
+  real    :: Z_diag_int  ! minimum interval between calc depth-space diagnostics (sec)
 
   real, allocatable, dimension(:,:)   :: eta ! free surface height or column mass [H ~> m or kg m-2]
   real, allocatable, dimension(:,:)   :: area_shelf_h ! area occupied by ice shelf [L2 ~> m2]

src/diagnostics/MOM_diagnostics.F90

@@ -1340,7 +1340,7 @@ end subroutine post_surface_thermo_diags
 
 !> This routine posts diagnostics of the transports, including the subgridscale
 !! contributions.
-subroutine post_transport_diagnostics(G, GV, uhtr, vhtr, h, IDs, diag_pre_dyn, diag, dt_trans, &
+subroutine post_transport_diagnostics(G, GV, US, uhtr, vhtr, h, IDs, diag_pre_dyn, diag, dt_trans, &
                                       diag_to_Z_CSp, Reg)
   type(ocean_grid_type),    intent(inout) :: G   !< ocean grid structure
   type(verticalGrid_type),  intent(in)    :: GV  !< ocean vertical grid structure

src/parameterizations/vertical/MOM_diabatic_driver.F90

@@ -229,6 +229,7 @@ module MOM_diabatic_driver
   type(tracer_flow_control_CS), pointer :: tracer_flow_CSp       => NULL() !< Control structure for a child module
   type(optics_type),            pointer :: optics                => NULL() !< Control structure for a child module
   type(KPP_CS),                 pointer :: KPP_CSp               => NULL() !< Control structure for a child module
+  type(diag_to_Z_CS),           pointer :: diag_to_Z_CSp         => NULL()
   type(tidal_mixing_cs),        pointer :: tidal_mixing_csp      => NULL() !< Control structure for a child module
   type(CVMix_conv_cs),          pointer :: CVMix_conv_csp        => NULL() !< Control structure for a child module
   type(diapyc_energy_req_CS),   pointer :: diapyc_en_rec_CSp     => NULL() !< Control structure for a child module
@@ -1926,7 +1927,7 @@ subroutine diabatic_ALE(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_end,
   if (num_z_diags > 0) &
     call calc_Zint_diags(h, z_ptrs, z_ids, num_z_diags, G, GV, CS%diag_to_Z_CSp)
 
-  if (CS%debugConservation) call MOM_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, G)
+  if (CS%debugConservation) call MOM_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, G, GV, US)
   if (showCallTree) call callTree_leave("diabatic()")
 
 end subroutine diabatic_ALE
@@ -2914,7 +2915,7 @@ subroutine layered_diabatic(u, v, h, tv, Hml, fluxes, visc, ADp, CDp, dt, Time_e
   if (num_z_diags > 0) &
     call calc_Zint_diags(h, z_ptrs, z_ids, num_z_diags, G, GV, CS%diag_to_Z_CSp)
 
-  if (CS%debugConservation) call MOM_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, G)
+  if (CS%debugConservation) call MOM_state_stats('leaving diabatic', u, v, h, tv%T, tv%S, G, GV, US)
   if (showCallTree) call callTree_leave("diabatic()")
 
 end subroutine layered_diabatic

src/tracer/MOM_tracer_Z_init.F90

@@ -3,7 +3,7 @@ module MOM_tracer_Z_init
 
 ! This file is part of MOM6. See LICENSE.md for the license.
 
-use MOM_diag_to_Z, only : find_overlap, find_limited_slope
+use MOM_diag_to_Z, only : find_overlap
 use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe
 ! use MOM_file_parser, only : get_param, log_version, param_file_type
 use MOM_grid, only : ocean_grid_type
@@ -509,75 +509,6 @@ subroutine read_Z_edges(filename, tr_name, z_edges, nz_out, has_edges, &
 
 end subroutine read_Z_edges
 
-!### `find_overlap` and `find_limited_slope` were previously part of
-!    MOM_diag_to_Z.F90, and are nearly identical to `find_overlap` in
-!    `midas_vertmap.F90` with some slight differences.  We keep it here for
-!    reproducibility, but the two should be merged at some point
-
-!> Determines the layers bounded by interfaces e that overlap
-!! with the depth range between Z_top and Z_bot, and the fractional weights
-!! of each layer. It also calculates the normalized relative depths of the range
-!! of each layer that overlaps that depth range.
-subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z2)
-  real, dimension(:), intent(in)  :: e      !< Column interface heights, [Z ~> m] or other units.
-  real,               intent(in)  :: Z_top  !< Top of range being mapped to, in the units of e [Z ~> m].
-  real,               intent(in)  :: Z_bot  !< Bottom of range being mapped to, in the units of e [Z ~> m].
-  integer,            intent(in)  :: k_max  !< Number of valid layers.
-  integer,            intent(in)  :: k_start !< Layer at which to start searching.
-  integer,            intent(out) :: k_top  !< Indices of top layers that overlap with the depth range.
-  integer,            intent(out) :: k_bot  !< Indices of bottom layers that overlap with the depth range.
-  real, dimension(:), intent(out) :: wt     !< Relative weights of each layer from k_top to k_bot [nondim].
-  real, dimension(:), intent(out) :: z1     !< Depth of the top limits of the part of
-       !! a layer that contributes to a depth level, relative to the cell center and normalized
-       !! by the cell thickness [nondim].  Note that -1/2 <= z1 < z2 <= 1/2.
-  real, dimension(:), intent(out)   :: z2     !< Depths of the bottom limit of the part of
-       !! a layer that contributes to a depth level, relative to the cell center and normalized
-       !! by the cell thickness [nondim].  Note that -1/2 <= z1 < z2 <= 1/2.
-  ! Local variables
-  real    :: Ih, e_c, tot_wt, I_totwt
-  integer :: k
-
-  wt(:) = 0.0 ; z1(:) = 0.0 ; z2(:) = 0.0 ; k_bot = k_max
-
-  do k=k_start,k_max ; if (e(K+1) < Z_top) exit ; enddo
-  k_top = k
-  if (k_top > k_max) return
-
-  ! Determine the fractional weights of each layer.
-  ! Note that by convention, e and Z_int decrease with increasing k.
-  if (e(K+1) <= Z_bot) then
-    wt(k) = 1.0 ; k_bot = k
-    Ih = 0.0 ; if (e(K) /= e(K+1)) Ih = 1.0 / (e(K)-e(K+1))
-    e_c = 0.5*(e(K)+e(K+1))
-    z1(k) = (e_c - MIN(e(K), Z_top)) * Ih
-    z2(k) = (e_c - Z_bot) * Ih
-  else
-    wt(k) = MIN(e(K),Z_top) - e(K+1) ; tot_wt = wt(k) ! These are always > 0.
-    if (e(K) /= e(K+1)) then
-      z1(k) = (0.5*(e(K)+e(K+1)) - MIN(e(K), Z_top)) / (e(K)-e(K+1))
-    else ; z1(k) = -0.5 ; endif
-    z2(k) = 0.5
-    k_bot = k_max
-    do k=k_top+1,k_max
-      if (e(K+1) <= Z_bot) then
-        k_bot = k
-        wt(k) = e(K) - Z_bot ; z1(k) = -0.5
-        if (e(K) /= e(K+1)) then
-          z2(k) = (0.5*(e(K)+e(K+1)) - Z_bot) / (e(K)-e(K+1))
-        else ; z2(k) = 0.5 ; endif
-      else
-        wt(k) = e(K) - e(K+1) ; z1(k) = -0.5 ; z2(k) = 0.5
-      endif
-      tot_wt = tot_wt + wt(k) ! wt(k) is always > 0.
-      if (k>=k_bot) exit
-    enddo
-
-    I_totwt = 0.0 ; if (tot_wt > 0.0) I_totwt = 1.0 / tot_wt
-    do k=k_top,k_bot ; wt(k) = I_totwt*wt(k) ; enddo
-  endif
-
-end subroutine find_overlap
-
 !> This subroutine determines a limited slope for val to be advected with
 !! a piecewise limited scheme.
 function find_limited_slope(val, e, k) result(slope)