diff --git a/src/core/MOM.F90 b/src/core/MOM.F90 index 4b7e0fdad1..900271bdf5 100644 --- a/src/core/MOM.F90 +++ b/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] diff --git a/src/diagnostics/MOM_diagnostics.F90 b/src/diagnostics/MOM_diagnostics.F90 index 5bfafd83b7..0f437d0307 100644 --- a/src/diagnostics/MOM_diagnostics.F90 +++ b/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 diff --git a/src/parameterizations/vertical/MOM_diabatic_driver.F90 b/src/parameterizations/vertical/MOM_diabatic_driver.F90 index 295e782632..9e7b1fad72 100644 --- a/src/parameterizations/vertical/MOM_diabatic_driver.F90 +++ b/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 diff --git a/src/tracer/MOM_tracer_Z_init.F90 b/src/tracer/MOM_tracer_Z_init.F90 index 9730c7ca30..7817c0d022 100644 --- a/src/tracer/MOM_tracer_Z_init.F90 +++ b/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)