src/parameterizations/lateral/MOM_mixed_layer_restrat.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

1		! This file is part of MOM6, the Modular Ocean Model version 6.
2		! See the LICENSE file for licensing information.
3		! SPDX-License-Identifier: Apache-2.0
4
5		!> \brief Parameterization of mixed layer restratification by unresolved mixed-layer eddies.
6		module MOM_mixed_layer_restrat
7
8		use MOM_debugging, only : hchksum
9		use MOM_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl
10		use MOM_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type
11		use MOM_diag_mediator, only : diag_update_remap_grids
12		use MOM_domains, only : pass_var, To_West, To_South, Omit_Corners
13		use MOM_error_handler, only : MOM_error, FATAL, WARNING
14		use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
15		use MOM_file_parser, only : openParameterBlock, closeParameterBlock
16		use MOM_forcing_type, only : mech_forcing, find_ustar
17		use MOM_grid, only : ocean_grid_type
18		use MOM_hor_index, only : hor_index_type
19		use MOM_interpolate, only : init_external_field, time_interp_external, time_interp_external_init
20		use MOM_interpolate, only : external_field
21		use MOM_intrinsic_functions, only : cuberoot
22		use MOM_io, only : slasher, MOM_read_data
23		use MOM_lateral_mixing_coeffs, only : VarMix_CS
24		use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS
25		use MOM_unit_scaling, only : unit_scale_type
26		use MOM_variables, only : thermo_var_ptrs
27		use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units
28		use MOM_EOS, only : calculate_density, calculate_spec_vol, EOS_domain
29
30		implicit none ; private
31
32		#include <MOM_memory.h>
33		#include "do_concurrent_compat.h"
34
35		#ifdef __NVCOMPILER_OPENMP_GPU
36		integer, parameter :: default_nkblock = 0 !< Default k block size for the mixed layer density integral [nondim]
37		integer, parameter :: default_njblock = 0 !< Default j block size for the mixed layer density integral [nondim]
38		#else
39		integer, parameter :: default_nkblock = 1 !< Default k block size for the mixed layer density integral [nondim]
40		integer, parameter :: default_njblock = 1 !< Default j block size for the mixed layer density integral [nondim]
41		#endif
42
43		public mixedlayer_restrat
44		public mixedlayer_restrat_init
45		public mixedlayer_restrat_register_restarts
46		public mixedlayer_restrat_unit_tests
47
48		! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
49		! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
50		! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
51		! vary with the Boussinesq approximation, the Boussinesq variant is given first.
52
53		!> Control structure for mom_mixed_layer_restrat
54		type, public :: mixedlayer_restrat_CS ; private
55		logical :: initialized = .false. !< True if this control structure has been initialized.
56		real :: ml_restrat_coef !< A non-dimensional factor by which the instability is enhanced
57		!! over what would be predicted based on the resolved gradients
58		!! [nondim]. This increases with grid spacing^2, up to something
59		!! of order 500.
60		real :: ml_restrat_coef2 !< As for ml_restrat_coef but using the slow filtered MLD [nondim].
61		real :: front_length !< If non-zero, is the frontal-length scale [L ~> m] used to calculate the
62		!! upscaling of buoyancy gradients that is otherwise represented
63		!! by the parameter FOX_KEMPER_ML_RESTRAT_COEF. If MLE_FRONT_LENGTH is
64		!! non-zero, it is recommended to set FOX_KEMPER_ML_RESTRAT_COEF=1.0.
65		logical :: fl_from_file !< If true, read the MLE front-length scale from a netCDF file.
66		logical :: MLE_use_PBL_MLD !< If true, use the MLD provided by the PBL parameterization.
67		!! if false, MLE will calculate a MLD based on a density difference
68		!! based on the parameter MLE_DENSITY_DIFF.
69		logical :: Bodner_detect_MLD !< If true, detect the MLD based on given density difference criterion
70		!! (MLE_DENSITY_DIFF) in the Bodner et al. parameterization.
71		real :: vonKar !< The von Karman constant as used for mixed layer viscosity [nondim]
72		real :: MLE_MLD_decay_time !< Time-scale to use in a running-mean when MLD is retreating [T ~> s].
73		real :: MLE_MLD_decay_time2 !< Time-scale to use in a running-mean when filtered MLD is retreating [T ~> s].
74		real :: MLE_density_diff !< Density difference used in detecting mixed-layer depth [R ~> kg m-3].
75		real :: MLE_tail_dh !< Fraction by which to extend the mixed-layer restratification
76		!! depth used for a smoother stream function at the base of
77		!! the mixed-layer [nondim].
78		real :: MLE_MLD_stretch !< A scaling coefficient for stretching/shrinking the MLD used in
79		!! the MLE scheme [nondim]. This simply multiplies MLD wherever used.
80
81		! The following parameters are used in the Bodner et al., 2023, parameterization
82		logical :: use_Bodner = .false. !< If true, use the Bodner et al., 2023, parameterization.
83		real :: Cr !< Efficiency coefficient from Bodner et al., 2023 [nondim]
84		real :: mstar !< The m* value used to estimate the turbulent vertical momentum flux [nondim]
85		real :: nstar !< The n* value used to estimate the turbulent vertical momentum flux [nondim]
86		real :: min_wstar2 !< The minimum lower bound to apply to the vertical momentum flux,
87		!! w'u', in the Bodner et al., restratification parameterization
88		!! [Z2 T-2 ~> m2 s-2]. This avoids a division-by-zero in the limit when u*
89		!! and the buoyancy flux are zero.
90		real :: BLD_growing_Tfilt !< The time-scale for a running-mean filter applied to the boundary layer
91		!! depth (BLD) when the BLD is deeper than the running mean [T ~> s].
92		!! A value of 0 instantaneously sets the running mean to the current value of BLD.
93		real :: BLD_decaying_Tfilt !< The time-scale for a running-mean filter applied to the boundary layer
94		!! depth (BLD) when the BLD is shallower than the running mean [T ~> s].
95		!! A value of 0 instantaneously sets the running mean to the current value of BLD.
96		real :: MLD_decaying_Tfilt !< The time-scale for a running-mean filter applied to the time-filtered
97		!! MLD, when the latter is shallower than the running mean [T ~> s].
98		!! A value of 0 instantaneously sets the running mean to the current value of MLD.
99		real :: MLD_growing_Tfilt !< The time-scale for a running-mean filter applied to the time-filtered
100		!! MLD, when the latter is deeper than the running mean [T ~> s].
101		!! A value of 0 instantaneously sets the running mean to the current value of MLD.
102		integer :: answer_date !< The vintage of the order of arithmetic and expressions in the
103		!! mixed layer restrat calculations. Values below 20240201 recover
104		!! the answers from the end of 2023, while higher values use the new
105		!! cuberoot function in the Bodner code to avoid needing to undo
106		!! dimensional rescaling.
107
108		logical :: debug = .false. !< If true, calculate checksums of fields for debugging.
109
110		type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
111		!! timing of diagnostic output.
112		type(external_field) :: sbc_fl !< A handle used in time interpolation of
113		!! front-length scales read from a file.
114		type(time_type), pointer :: Time => NULL() !< A pointer to the ocean model's clock.
115		logical :: use_Stanley_ML !< If true, use the Stanley parameterization of SGS T variance
116		real :: ustar_min !< A minimum value of ustar in thickness units to avoid numerical
117		!! problems [H T-1 ~> m s-1 or kg m-2 s-1]
118		real :: Kv_restrat !< A viscosity that sets a floor on the momentum mixing rate
119		!! during restratification, rescaled into thickness-based
120		!! units [H2 T-1 ~> m2 s-1 or kg2 m-4 s-1]
121		logical :: MLD_grid !< If true, read a spacially varying field for MLD_decaying_Tfilt
122		logical :: Cr_grid !< If true, read a spacially varying field for Cr
123		integer :: nkblock = default_nkblock !< The number of layers whose density is evaluated in one call to the
124		!! equation of state in the Bodner mixed-layer buoyancy integral, or 0
125		!! to do the whole column in a single block [nondim].
126		integer :: njblock = default_njblock !< The number of rows whose density is evaluated in one call to the
127		!! equation of state in the Bodner mixed-layer buoyancy integral, or 0
128		!! to do the whole j compute domain in a single block [nondim].
129
130		real, dimension(:,:), allocatable :: &
131		MLD_filtered, & !< Time-filtered MLD [H ~> m or kg m-2]
132		MLD_filtered_slow, & !< Slower time-filtered MLD [H ~> m or kg m-2]
133		wpup_filtered, & !< Time-filtered vertical momentum flux [H L T-2 ~> m2 s-2 or kg m-1 s-2]
134		MLD_Tfilt_space, & !< Spatially varying time scale for MLD filter [T ~> s]
135		Cr_space !< Spatially varying Cr coefficient [nondim]
136
137		!>@{
138		!! Diagnostic identifier
139		integer :: id_urestrat_time = -1
140		integer :: id_vrestrat_time = -1
141		integer :: id_uhml = -1
142		integer :: id_vhml = -1
143		integer :: id_MLD = -1
144		integer :: id_BLD = -1
145		integer :: id_Rml = -1
146		integer :: id_uDml = -1
147		integer :: id_vDml = -1
148		integer :: id_uml = -1
149		integer :: id_vml = -1
150		integer :: id_wpup = -1
151		integer :: id_ustar = -1
152		integer :: id_bflux = -1
153		integer :: id_lfbod = -1
154		integer :: id_mle_fl = -1
155		!>@}
156
157		end type mixedlayer_restrat_CS
158
159		character(len=40) :: mdl = "MOM_mixed_layer_restrat" !< This module's name.
160
161		contains
162
163		!> Driver for the mixed-layer restratification parameterization.
164		!! The code branches between two different implementations depending
165		!! on whether the bulk-mixed layer or a general coordinate are in use.
166	24	subroutine mixedlayer_restrat(h, uhtr, vhtr, tv, forces, dt, MLD, h_MLD, bflux, VarMix, G, GV, US, CS)
167		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
168		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
169		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
170		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
171		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: uhtr !< Accumulated zonal mass flux
172		!! [H L2 ~> m3 or kg]
173		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: vhtr !< Accumulated meridional mass flux
174		!! [H L2 ~> m3 or kg]
175		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables structure
176		type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
177		real, intent(in) :: dt !< Time increment [T ~> s]
178		real, dimension(:,:), pointer :: MLD !< Mixed layer depth provided by the
179		!! planetary boundary layer scheme [Z ~> m]
180		real, dimension(:,:), pointer :: h_MLD !< Mixed layer thickness provided
181		!! by the planetary boundary layer
182		!! scheme [H ~> m or kg m-2]
183		real, dimension(:,:), pointer :: bflux !< Surface buoyancy flux provided by the
184		!! PBL scheme [Z2 T-3 ~> m2 s-3]
185		type(VarMix_CS), intent(in) :: VarMix !< Variable mixing control structure
186		type(mixedlayer_restrat_CS), intent(inout) :: CS !< Module control structure
187
188
189	24	if (.not. CS%initialized) call MOM_error(FATAL, "mixedlayer_restrat: "// &
190	0	"Module must be initialized before it is used.")
191
192	24	if (GV%nkml>0) then
193		! Original form, written for the isopycnal model with a bulk mixed layer
194	0	call mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
195	24	elseif (CS%use_Bodner) then
196		! Implementation of Bodner et al., 2023
197	24	call mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, dt, MLD, h_MLD, bflux)
198		else
199		! Implementation of Fox-Kemper et al., 2008, to work in general coordinates
200	0	call mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, h_MLD, VarMix, G, GV, US, CS)
201		endif
202
203	24	end subroutine mixedlayer_restrat
204
205		!> Calculates a restratifying flow in the mixed layer, following the formulation used in OM4
206	0	subroutine mixedlayer_restrat_OM4(h, uhtr, vhtr, tv, forces, dt, h_MLD, VarMix, G, GV, US, CS)
207		! Arguments
208		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
209		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
210		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
211		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
212		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: uhtr !< Accumulated zonal mass flux
213		!! [H L2 ~> m3 or kg]
214		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: vhtr !< Accumulated meridional mass flux
215		!! [H L2 ~> m3 or kg]
216		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables structure
217		type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
218		real, intent(in) :: dt !< Time increment [T ~> s]
219		real, dimension(:,:), pointer :: h_MLD !< Thickness of water within the
220		!! mixed layer depth provided by
221		!! the PBL scheme [H ~> m or kg m-2]
222		type(VarMix_CS), intent(in) :: VarMix !< Variable mixing control structure
223		type(mixedlayer_restrat_CS), intent(inout) :: CS !< Module control structure
224
225		! Local variables
226	0	real :: uhml(SZIB_(G),SZJ_(G),SZK_(GV)) ! Restratifying zonal thickness transports [H L2 T-1 ~> m3 s-1 or kg s-1]
227	0	real :: vhml(SZI_(G),SZJB_(G),SZK_(GV)) ! Restratifying meridional thickness transports [H L2 T-1 ~> m3 s-1 or kg s-1]
228		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
229	0	h_avail ! The volume available for diffusion out of each face of each
230		! sublayer of the mixed layer, divided by dt [H L2 T-1 ~> m3 s-1 or kg s-1].
231		real, dimension(SZI_(G),SZJ_(G)) :: &
232	0	U_star_2d, & ! The wind friction velocity in thickness-based units, calculated using
233		! the Boussinesq reference density or the time-evolving surface density
234		! in non-Boussinesq mode [H T-1 ~> m s-1 or kg m-2 s-1]
235	0	MLD_fast, & ! Mixed layer depth actually used in MLE restratification parameterization [H ~> m or kg m-2]
236	0	htot_fast, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2]
237	0	Rml_av_fast, & ! Negative g_Rho0 times the average mixed layer density or G_Earth
238		! times the average specific volume [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2]
239	0	MLD_slow, & ! Mixed layer depth actually used in MLE restratification parameterization [H ~> m or kg m-2]
240	0	mle_fl_2d, & ! MLE frontal length-scale [L ~> m]
241	0	htot_slow, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2]
242	0	Rml_av_slow ! Negative g_Rho0 times the average mixed layer density or G_Earth
243		! times the average specific volume [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2]
244		real :: g_Rho0 ! G_Earth/Rho0 times a thickness conversion factor
245		! [L2 H-1 T-2 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2]
246	0	real :: rho_ml(SZI_(G)) ! Potential density relative to the surface [R ~> kg m-3]
247	0	real :: rml_int_fast(SZI_(G)) ! The integral of density over the mixed layer depth [R H ~> kg m-2 or kg2 m-5]
248	0	real :: rml_int_slow(SZI_(G)) ! The integral of density over the mixed layer depth [R H ~> kg m-2 or kg2 m-5]
249	0	real :: SpV_ml(SZI_(G)) ! Specific volume evaluated at the surface pressure [R-1 ~> m3 kg-1]
250	0	real :: SpV_int_fast(SZI_(G)) ! Specific volume integrated through the mixed layer [H R-1 ~> m4 kg-1 or m]
251	0	real :: SpV_int_slow(SZI_(G)) ! Specific volume integrated through the mixed layer [H R-1 ~> m4 kg-1 or m]
252	0	real :: p0(SZI_(G)) ! A pressure of 0 [R L2 T-2 ~> Pa]
253
254		real :: h_vel ! htot interpolated onto velocity points [H ~> m or kg m-2]
255		real :: absf ! absolute value of f, interpolated to velocity points [T-1 ~> s-1]
256		real :: u_star ! surface friction velocity, interpolated to velocity points and recast into
257		! thickness-based units [H T-1 ~> m s-1 or kg m-2 s-1].
258		real :: mom_mixrate ! rate at which momentum is homogenized within mixed layer [T-1 ~> s-1]
259		real :: timescale ! mixing growth timescale [T ~> s]
260		real :: h_min ! The minimum layer thickness [H ~> m or kg m-2]. h_min could be 0.
261		real :: h_neglect ! tiny thickness usually lost in roundoff so can be neglected [H ~> m or kg m-2]
262		real :: I4dt ! 1/(4 dt) [T-1 ~> s-1]
263		real :: Ihtot, Ihtot_slow ! Inverses of the total mixed layer thickness [H-1 ~> m-1 or m2 kg-1]
264	0	real :: a(SZK_(GV)) ! A non-dimensional value relating the overall flux
265		! magnitudes (uDml & vDml) to the realized flux in a
266		! layer [nondim]. The vertical sum of a() through the pieces of
267		! the mixed layer must be 0.
268	0	real :: b(SZK_(GV)) ! As for a(k) but for the slow-filtered MLD [nondim]
269	0	real :: uDml(SZIB_(G)) ! Zonal volume fluxes in the upper half of the mixed layer [H L2 T-1 ~> m3 s-1 or kg s-1]
270	0	real :: vDml(SZI_(G)) ! Meridional volume fluxes in the upper half of the mixed layer [H L2 T-1 ~> m3 s-1 or kg s-1]
271	0	real :: uDml_slow(SZIB_(G)) ! Zonal volume fluxes in the upper half of the boundary layer to
272		! restratify the time-filtered boundary layer depth [H L2 T-1 ~> m3 s-1 or kg s-1]
273	0	real :: vDml_slow(SZI_(G)) ! Meridional volume fluxes in the upper half of the boundary layer to
274		! restratify the time-filtered boundary layer depth [H L2 T-1 ~> m3 s-1 or kg s-1]
275	0	real :: utimescale_diag(SZIB_(G),SZJ_(G)) ! Zonal restratification timescale [T ~> s], stored for diagnostics.
276	0	real :: vtimescale_diag(SZI_(G),SZJB_(G)) ! Meridional restratification timescale [T ~> s], stored for diagnostics.
277	0	real :: uDml_diag(SZIB_(G),SZJ_(G)) ! A 2D copy of uDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1]
278	0	real :: vDml_diag(SZI_(G),SZJB_(G)) ! A 2D copy of vDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1]
279	0	real, dimension(SZI_(G)) :: covTS, & ! SGS TS covariance in Stanley param; currently 0 [C S ~> degC ppt]
280	0	varS ! SGS S variance in Stanley param; currently 0 [S2 ~> ppt2]
281		real :: aFac, bFac ! Nondimensional ratios [nondim]
282		real :: hAtVel ! Thickness at the velocity points [H ~> m or kg m-2]
283		real :: zpa ! Fractional position within the mixed layer of the interface above a layer [nondim]
284		real :: zpb ! Fractional position within the mixed layer of the interface below a layer [nondim]
285		real :: dh ! Portion of the layer thickness that is in the mixed layer [H ~> m or kg m-2]
286		real :: res_scaling_fac ! The resolution-dependent scaling factor [nondim]
287		real :: lfront ! Frontal length scale at velocity points [L ~> m]
288		real :: I_LFront ! The inverse of the frontal length scale [L-1 ~> m-1]
289		real :: vonKar_x_pi2 ! A scaling constant that is approximately the von Karman constant times
290		! pi squared [nondim]
291		character(len=128) :: mesg
292		logical :: line_is_empty, keep_going, res_upscale
293		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
294		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
295
296	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
297	0	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
298
299	0	h_min = 0.5*GV%Angstrom_H ! This should be GV%Angstrom_H, but that value would change answers.
300	0	covTS(:) = 0.0 !!Functionality not implemented yet; in future, should be passed in tv
301	0	varS(:) = 0.0
302	0	mle_fl_2d(:,:) = 0.0
303	0	vonKar_x_pi2 = CS%vonKar * 9.8696
304
305	0	if (.not.associated(tv%eqn_of_state)) call MOM_error(FATAL, "mixedlayer_restrat_OM4: "// &
306	0	"An equation of state must be used with this module.")
307	0	if (.not. allocated(VarMix%Rd_dx_h) .and. CS%front_length > 0.) &
308		call MOM_error(FATAL, "mixedlayer_restrat_OM4: "// &
309	0	"The resolution argument, Rd/dx, was not associated.")
310	0	if (CS%use_Stanley_ML .and. .not.GV%Boussinesq) call MOM_error(FATAL, &
311		"MOM_mixedlayer_restrat: The Stanley parameterization is not "//&
312	0	"available without the Boussinesq approximation.")
313
314		! Extract the friction velocity from the forcing type.
315	0	call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1, H_T_units=.true.)
316
317	0	if (CS%MLE_density_diff > 0.) then ! We need to calculate a mixed layer depth, MLD.
318	0	call detect_mld(h, tv, MLD_fast, G, GV, CS)
319	0	elseif (CS%MLE_use_PBL_MLD) then
320	0	do j=js-1,je+1 ; do i=is-1,ie+1
321	0	MLD_fast(i,j) = CS%MLE_MLD_stretch * h_MLD(i,j)
322		enddo ; enddo
323		else
324		call MOM_error(FATAL, "mixedlayer_restrat_OM4: "// &
325	0	"No MLD to use for MLE parameterization.")
326		endif
327
328		! Apply time filter (to remove diurnal cycle)
329	0	if (CS%MLE_MLD_decay_time>0.) then
330	0	if (CS%debug) then
331	0	call hchksum(CS%MLD_filtered, 'mixed_layer_restrat: MLD_filtered', G%HI, haloshift=1, unscale=GV%H_to_mks)
332	0	call hchksum(h_MLD, 'mixed_layer_restrat: MLD in', G%HI, haloshift=1, unscale=GV%H_to_mks)
333		endif
334	0	aFac = CS%MLE_MLD_decay_time / ( dt + CS%MLE_MLD_decay_time )
335	0	bFac = dt / ( dt + CS%MLE_MLD_decay_time )
336	0	do j=js-1,je+1 ; do i=is-1,ie+1
337		! Expression bFacMLD_fast(i,j) + aFacCS%MLD_filtered(i,j) is the time-filtered
338		! (running mean) of MLD. The max() allows the "running mean" to be reset
339		! instantly to a deeper MLD.
340	0	CS%MLD_filtered(i,j) = max( MLD_fast(i,j), bFacMLD_fast(i,j) + aFacCS%MLD_filtered(i,j) )
341	0	MLD_fast(i,j) = CS%MLD_filtered(i,j)
342		enddo ; enddo
343		endif
344
345		! Apply slower time filter (to remove seasonal cycle) on already filtered MLD_fast
346	0	if (CS%MLE_MLD_decay_time2>0.) then
347	0	if (CS%debug) then
348		call hchksum(CS%MLD_filtered_slow, 'mixed_layer_restrat: MLD_filtered_slow', G%HI, &
349	0	haloshift=1, unscale=GV%H_to_mks)
350	0	call hchksum(MLD_fast, 'mixed_layer_restrat: MLD fast', G%HI, haloshift=1, unscale=GV%H_to_mks)
351		endif
352	0	aFac = CS%MLE_MLD_decay_time2 / ( dt + CS%MLE_MLD_decay_time2 )
353	0	bFac = dt / ( dt + CS%MLE_MLD_decay_time2 )
354	0	do j=js-1,je+1 ; do i=is-1,ie+1
355		! Expression bFacMLD_fast(i,j) + aFacCS%MLD_filtered(i,j) is the time-filtered
356		! (running mean) of MLD. The max() allows the "running mean" to be reset
357		! instantly to a deeper MLD.
358	0	CS%MLD_filtered_slow(i,j) = max( MLD_fast(i,j), bFacMLD_fast(i,j) + aFacCS%MLD_filtered_slow(i,j) )
359	0	MLD_slow(i,j) = CS%MLD_filtered_slow(i,j)
360		enddo ; enddo
361		else
362	0	do j=js-1,je+1 ; do i=is-1,ie+1
363	0	MLD_slow(i,j) = MLD_fast(i,j)
364		enddo ; enddo
365		endif
366
367	0	uDml(:) = 0.0 ; vDml(:) = 0.0
368	0	uDml_slow(:) = 0.0 ; vDml_slow(:) = 0.0
369	0	I4dt = 0.25 / dt
370	0	g_Rho0 = GV%H_to_Z * GV%g_Earth / GV%Rho0
371	0	h_neglect = GV%H_subroundoff
372	0	if (CS%front_length>0.) then
373	0	res_upscale = .true.
374	0	do j=js-1,je+1 ; do i=is-1,ie+1
375	0	mle_fl_2d(i,j) = CS%front_length
376		enddo ; enddo
377	0	elseif (CS%front_length == 0. .and. CS%fl_from_file) then
378	0	res_upscale = .true.
379	0	call time_interp_external(CS%sbc_fl, CS%Time, mle_fl_2d, turns=G%HI%turns, scale=US%m_to_L)
380	0	call pass_var(mle_fl_2d, G%domain, halo=1)
381	0	do j=js,je ; do i=is,ie
382	0	if ((G%mask2dT(i,j) > 0.0) .and. (mle_fl_2d(i,j) < 0.0)) then
383		write(mesg,'(" Time_interp negative MLE frontal-length scale of ",(1pe12.4)," at i,j = ",&
384		& I0,", ",I0," lon/lat = ",(1pe12.4)," E ", (1pe12.4), " N.")') &
385	0	mle_fl_2d(i,j), i, j, G%geoLonT(i,j), G%geoLatT(i,j)
386	0	call MOM_error(FATAL, "MOM_mixed_layer_restrat mixedlayer_restrat_OM4: "//trim(mesg))
387		endif
388		enddo ; enddo
389		else
390	0	res_upscale = .false.
391		endif
392
393	0	p0(:) = 0.0
394	0	EOSdom(:) = EOS_domain(G%HI, halo=1)
395		!$OMP parallel default(shared) private(rho_ml,h_vel,u_star,absf,mom_mixrate,timescale, &
396		!$OMP SpV_ml,SpV_int_fast,SpV_int_slow,Rml_int_fast,Rml_int_slow, &
397		!$OMP line_is_empty,keep_going,res_scaling_fac, &
398		!$OMP a,IhTot,b,Ihtot_slow,zpb,hAtVel,zpa,dh,lfront,I_LFront) &
399		!$OMP firstprivate(uDml,vDml,uDml_slow,vDml_slow)
400
401	0	if (GV%Boussinesq .or. GV%semi_Boussinesq) then
402		!$OMP do
403	0	do j=js-1,je+1
404	0	do i=is-1,ie+1
405	0	htot_fast(i,j) = 0.0 ; Rml_int_fast(i) = 0.0
406	0	htot_slow(i,j) = 0.0 ; Rml_int_slow(i) = 0.0
407		enddo
408	0	keep_going = .true.
409	0	do k=1,nz
410	0	do i=is-1,ie+1
411	0	h_avail(i,j,k) = max(I4dtG%areaT(i,j)(h(i,j,k)-GV%Angstrom_H),0.0)
412		enddo
413	0	if (keep_going) then
414	0	if (CS%use_Stanley_ML) then
415		call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, tv%varT(:,j,k), covTS, varS, &
416	0	rho_ml(:), tv%eqn_of_state, EOSdom)
417		else
418	0	call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, rho_ml(:), tv%eqn_of_state, EOSdom)
419		endif
420	0	line_is_empty = .true.
421	0	do i=is-1,ie+1
422	0	if (htot_fast(i,j) < MLD_fast(i,j)) then
423	0	dh = min( h(i,j,k), MLD_fast(i,j)-htot_fast(i,j) )
424	0	Rml_int_fast(i) = Rml_int_fast(i) + dh*rho_ml(i)
425	0	htot_fast(i,j) = htot_fast(i,j) + dh
426	0	line_is_empty = .false.
427		endif
428	0	if (htot_slow(i,j) < MLD_slow(i,j)) then
429	0	dh = min( h(i,j,k), MLD_slow(i,j)-htot_slow(i,j) )
430	0	Rml_int_slow(i) = Rml_int_slow(i) + dh*rho_ml(i)
431	0	htot_slow(i,j) = htot_slow(i,j) + dh
432	0	line_is_empty = .false.
433		endif
434		enddo
435	0	if (line_is_empty) keep_going=.false.
436		endif
437		enddo
438
439	0	do i=is-1,ie+1
440	0	Rml_av_fast(i,j) = -(g_Rho0*Rml_int_fast(i)) / (htot_fast(i,j) + h_neglect)
441	0	Rml_av_slow(i,j) = -(g_Rho0*Rml_int_slow(i)) / (htot_slow(i,j) + h_neglect)
442		enddo
443		enddo
444		else ! This is only used in non-Boussinesq mode.
445		!$OMP do
446	0	do j=js-1,je+1
447	0	do i=is-1,ie+1
448	0	htot_fast(i,j) = 0.0 ; SpV_int_fast(i) = 0.0
449	0	htot_slow(i,j) = 0.0 ; SpV_int_slow(i) = 0.0
450		enddo
451	0	keep_going = .true.
452	0	do k=1,nz
453	0	do i=is-1,ie+1
454	0	h_avail(i,j,k) = max(I4dtG%areaT(i,j)(h(i,j,k)-GV%Angstrom_H),0.0)
455		enddo
456	0	if (keep_going) then
457		! if (CS%use_Stanley_ML) then ! This is not implemented yet in the EoS code.
458		! call calculate_spec_vol(tv%T(:,j,k), tv%S(:,j,k), p0, tv%varT(:,j,k), covTS, varS, &
459		! rho_ml(:), tv%eqn_of_state, EOSdom)
460		! else
461	0	call calculate_spec_vol(tv%T(:,j,k), tv%S(:,j,k), p0, SpV_ml, tv%eqn_of_state, EOSdom)
462		! endif
463	0	line_is_empty = .true.
464	0	do i=is-1,ie+1
465	0	if (htot_fast(i,j) < MLD_fast(i,j)) then
466	0	dh = min( h(i,j,k), MLD_fast(i,j)-htot_fast(i,j) )
467	0	SpV_int_fast(i) = SpV_int_fast(i) + dh*SpV_ml(i)
468	0	htot_fast(i,j) = htot_fast(i,j) + dh
469	0	line_is_empty = .false.
470		endif
471	0	if (htot_slow(i,j) < MLD_slow(i,j)) then
472	0	dh = min( h(i,j,k), MLD_slow(i,j)-htot_slow(i,j) )
473	0	SpV_int_slow(i) = SpV_int_slow(i) + dh*SpV_ml(i)
474	0	htot_slow(i,j) = htot_slow(i,j) + dh
475	0	line_is_empty = .false.
476		endif
477		enddo
478	0	if (line_is_empty) keep_going=.false.
479		endif
480		enddo
481
482		! Convert the vertically integrated specific volume into a positive variable with units of density.
483	0	do i=is-1,ie+1
484	0	Rml_av_fast(i,j) = (GV%H_to_RZGV%g_Earth SpV_int_fast(i)) / (htot_fast(i,j) + h_neglect)
485	0	Rml_av_slow(i,j) = (GV%H_to_RZGV%g_Earth SpV_int_slow(i)) / (htot_slow(i,j) + h_neglect)
486		enddo
487		enddo
488		endif
489
490	0	if (CS%debug) then
491	0	call hchksum(h, 'mixed_layer_restrat: h', G%HI, haloshift=1, unscale=GV%H_to_mks)
492	0	call hchksum(U_star_2d, 'mixed_layer_restrat: u', G%HI, haloshift=1, unscale=GV%H_to_mUS%s_to_T)
493	0	call hchksum(MLD_fast, 'mixed_layer_restrat: MLD', G%HI, haloshift=1, unscale=GV%H_to_mks)
494		call hchksum(Rml_av_fast, 'mixed_layer_restrat: rml', G%HI, haloshift=1, &
495	0	unscale=GV%m_to_HUS%L_T_to_m_s*2)
496		endif
497
498		! TO DO:
499		! 1. Mixing extends below the mixing layer to the mixed layer. Find it!
500		! 2. Add exponential tail to stream-function?
501
502		! U - Component
503		!$OMP do
504	0	do j=js,je ; do I=is-1,ie
505	0	u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i+1,j)))
506
507	0	absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J)))
508		! Compute I_LFront = 1 / (frontal length scale) [L-1 ~> m-1]
509	0	lfront = 0.5 * (mle_fl_2d(i,j) + mle_fl_2d(i+1,j))
510		! Adcroft reciprocal
511	0	I_LFront = 0.0 ; if (lfront /= 0.0) I_LFront = 1.0/lfront
512		! If needed, res_scaling_fac = min( ds, L_d ) / l_f
513	0	if (res_upscale) res_scaling_fac = &
514		( sqrt( 0.5 * ( (G%dxCu(I,j)2) + (G%dyCu(I,j)2) ) ) * I_LFront ) &
515	0	* min( 1., 0.5*( VarMix%Rd_dx_h(i,j) + VarMix%Rd_dx_h(i+1,j) ) )
516
517		! peak ML visc: u_star * von_Karman * (h_mlu_star)/(absfh_ml + 4.0*u_star)
518		! momentum mixing rate: pi^2*visc/h_ml^2
519	0	h_vel = 0.5*((htot_fast(i,j) + htot_fast(i+1,j)) + h_neglect)
520
521		! NOTE: growth_time changes answers on some systems, see below.
522		! timescale = growth_time(u_star, h_vel, absf, h_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef)
523
524		mom_mixrate = vonKar_x_pi2u_star*2 / &
525	0	(absfh_vel2 + 4.0(h_vel+h_neglect)*u_star)
526	0	timescale = 0.0625 * (absf + 2.0mom_mixrate) / (absf2 + mom_mixrate*2)
527	0	timescale = timescale * CS%ml_restrat_coef
528
529	0	if (res_upscale) timescale = timescale * res_scaling_fac
530		uDml(I) = timescale * G%dyCu(I,j)G%IdxCu_OBCmask(I,j) &
531	0	(Rml_av_fast(i+1,j)-Rml_av_fast(i,j)) * (h_vel**2)
532
533		! As above but using the slow filtered MLD
534	0	h_vel = 0.5*((htot_slow(i,j) + htot_slow(i+1,j)) + h_neglect)
535
536		! NOTE: growth_time changes answers on some systems, see below.
537		! timescale = growth_time(u_star, h_vel, absf, h_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef2)
538
539		mom_mixrate = vonKar_x_pi2u_star*2 / &
540	0	(absfh_vel2 + 4.0(h_vel+h_neglect)*u_star)
541	0	timescale = 0.0625 * (absf + 2.0mom_mixrate) / (absf2 + mom_mixrate*2)
542	0	timescale = timescale * CS%ml_restrat_coef2
543
544	0	if (res_upscale) timescale = timescale * res_scaling_fac
545		uDml_slow(I) = timescale * G%dyCu(I,j)G%IdxCu_OBCmask(I,j) &
546	0	(Rml_av_slow(i+1,j)-Rml_av_slow(i,j)) * (h_vel**2)
547
548	0	if (uDml(I) + uDml_slow(I) == 0.) then
549	0	do k=1,nz ; uhml(I,j,k) = 0.0 ; enddo
550		else
551	0	IhTot = 2.0 / ((htot_fast(i,j) + htot_fast(i+1,j)) + h_neglect)
552	0	IhTot_slow = 2.0 / ((htot_slow(i,j) + htot_slow(i+1,j)) + h_neglect)
553	0	zpa = 0.0 ; zpb = 0.0
554		! a(k) relates the sublayer transport to uDml with a linear profile.
555		! The sum of a(k) through the mixed layers must be 0.
556	0	do k=1,nz
557	0	hAtVel = 0.5*(h(i,j,k) + h(i+1,j,k))
558	0	a(k) = mu(zpa, CS%MLE_tail_dh) ! mu(z/MLD) for upper interface
559	0	zpa = zpa - (hAtVel * IhTot) ! z/H for lower interface
560	0	a(k) = a(k) - mu(zpa, CS%MLE_tail_dh) ! Transport profile
561		! Limit magnitude (uDml) if it would violate CFL
562	0	if (a(k)*uDml(I) > 0.0) then
563	0	if (a(k)*uDml(I) > h_avail(i,j,k)) uDml(I) = h_avail(i,j,k) / a(k)
564	0	elseif (a(k)*uDml(I) < 0.0) then
565	0	if (-a(k)*uDml(I) > h_avail(i+1,j,k)) uDml(I) = -h_avail(i+1,j,k) / a(k)
566		endif
567		enddo
568	0	do k=1,nz
569		! Transport for slow-filtered MLD
570	0	hAtVel = 0.5*(h(i,j,k) + h(i+1,j,k))
571	0	b(k) = mu(zpb, CS%MLE_tail_dh) ! mu(z/MLD) for upper interface
572	0	zpb = zpb - (hAtVel * IhTot_slow) ! z/H for lower interface
573	0	b(k) = b(k) - mu(zpb, CS%MLE_tail_dh) ! Transport profile
574		! Limit magnitude (uDml_slow) if it would violate CFL when added to uDml
575	0	if (b(k)*uDml_slow(I) > 0.0) then
576	0	if (b(k)uDml_slow(I) > h_avail(i,j,k) - a(k)uDml(I)) &
577	0	uDml_slow(I) = max( 0., h_avail(i,j,k) - a(k)*uDml(I) ) / b(k)
578	0	elseif (b(k)*uDml_slow(I) < 0.0) then
579	0	if (-b(k)uDml_slow(I) > h_avail(i+1,j,k) + a(k)uDml(I)) &
580	0	uDml_slow(I) = -max( 0., h_avail(i+1,j,k) + a(k)*uDml(I) ) / b(k)
581		endif
582		enddo
583	0	do k=1,nz
584	0	uhml(I,j,k) = a(k)uDml(I) + b(k)uDml_slow(I)
585	0	uhtr(I,j,k) = uhtr(I,j,k) + uhml(I,j,k)*dt
586		enddo
587		endif
588
589	0	utimescale_diag(I,j) = timescale
590	0	uDml_diag(I,j) = uDml(I)
591		enddo ; enddo
592
593		! V- component
594		!$OMP do
595	0	do J=js-1,je ; do i=is,ie
596	0	u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i,j+1)))
597		! Compute I_LFront = 1 / (frontal length scale) [L-1 ~> m-1]
598	0	lfront = 0.5 * (mle_fl_2d(i,j) + mle_fl_2d(i,j+1))
599		! Adcroft reciprocal
600	0	I_LFront = 0.0 ; if (lfront /= 0.0) I_LFront = 1.0/lfront
601	0	absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))
602		! If needed, res_scaling_fac = min( ds, L_d ) / l_f
603	0	if (res_upscale) res_scaling_fac = &
604		( sqrt( 0.5 * ( (G%dxCv(i,J)2) + (G%dyCv(i,J)2) ) ) * I_LFront ) &
605	0	* min( 1., 0.5*( VarMix%Rd_dx_h(i,j) + VarMix%Rd_dx_h(i,j+1) ) )
606
607		! peak ML visc: u_star * von_Karman * (h_mlu_star)/(absfh_ml + 4.0*u_star)
608		! momentum mixing rate: pi^2*visc/h_ml^2
609	0	h_vel = 0.5*((htot_fast(i,j) + htot_fast(i,j+1)) + h_neglect)
610
611		! NOTE: growth_time changes answers on some systems, see below.
612		! timescale = growth_time(u_star, h_vel, absf, h_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef)
613
614		mom_mixrate = vonKar_x_pi2u_star*2 / &
615	0	(absfh_vel2 + 4.0(h_vel+h_neglect)*u_star)
616	0	timescale = 0.0625 * (absf + 2.0mom_mixrate) / (absf2 + mom_mixrate*2)
617	0	timescale = timescale * CS%ml_restrat_coef
618
619	0	if (res_upscale) timescale = timescale * res_scaling_fac
620		vDml(i) = timescale * G%dxCv(i,J)G%IdyCv_OBCmask(i,J) &
621	0	(Rml_av_fast(i,j+1)-Rml_av_fast(i,j)) * (h_vel**2)
622
623		! As above but using the slow filtered MLD
624	0	h_vel = 0.5*((htot_slow(i,j) + htot_slow(i,j+1)) + h_neglect)
625
626		! NOTE: growth_time changes answers on some systems, see below.
627		! timescale = growth_time(u_star, h_vel, absf, h_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef2)
628
629		mom_mixrate = vonKar_x_pi2u_star*2 / &
630	0	(absfh_vel2 + 4.0(h_vel+h_neglect)*u_star)
631	0	timescale = 0.0625 * (absf + 2.0mom_mixrate) / (absf2 + mom_mixrate*2)
632	0	timescale = timescale * CS%ml_restrat_coef2
633
634	0	if (res_upscale) timescale = timescale * res_scaling_fac
635		vDml_slow(i) = timescale * G%dxCv(i,J)G%IdyCv_OBCmask(i,J) &
636	0	(Rml_av_slow(i,j+1)-Rml_av_slow(i,j)) * (h_vel**2)
637
638	0	if (vDml(i) + vDml_slow(i) == 0.) then
639	0	do k=1,nz ; vhml(i,J,k) = 0.0 ; enddo
640		else
641	0	IhTot = 2.0 / ((htot_fast(i,j) + htot_fast(i,j+1)) + h_neglect)
642	0	IhTot_slow = 2.0 / ((htot_slow(i,j) + htot_slow(i,j+1)) + h_neglect)
643	0	zpa = 0.0 ; zpb = 0.0
644		! a(k) relates the sublayer transport to vDml with a linear profile.
645		! The sum of a(k) through the mixed layers must be 0.
646	0	do k=1,nz
647	0	hAtVel = 0.5*(h(i,j,k) + h(i,j+1,k))
648	0	a(k) = mu(zpa, CS%MLE_tail_dh) ! mu(z/MLD) for upper interface
649	0	zpa = zpa - (hAtVel * IhTot) ! z/H for lower interface
650	0	a(k) = a(k) - mu(zpa, CS%MLE_tail_dh) ! Transport profile
651		! Limit magnitude (vDml) if it would violate CFL
652	0	if (a(k)*vDml(i) > 0.0) then
653	0	if (a(k)*vDml(i) > h_avail(i,j,k)) vDml(i) = h_avail(i,j,k) / a(k)
654	0	elseif (a(k)*vDml(i) < 0.0) then
655	0	if (-a(k)*vDml(i) > h_avail(i,j+1,k)) vDml(i) = -h_avail(i,j+1,k) / a(k)
656		endif
657		enddo
658	0	do k=1,nz
659		! Transport for slow-filtered MLD
660	0	hAtVel = 0.5*(h(i,j,k) + h(i,j+1,k))
661	0	b(k) = mu(zpb, CS%MLE_tail_dh) ! mu(z/MLD) for upper interface
662	0	zpb = zpb - (hAtVel * IhTot_slow) ! z/H for lower interface
663	0	b(k) = b(k) - mu(zpb, CS%MLE_tail_dh) ! Transport profile
664		! Limit magnitude (vDml_slow) if it would violate CFL when added to vDml
665	0	if (b(k)*vDml_slow(i) > 0.0) then
666	0	if (b(k)vDml_slow(i) > h_avail(i,j,k) - a(k)vDml(i)) &
667	0	vDml_slow(i) = max( 0., h_avail(i,j,k) - a(k)*vDml(i) ) / b(k)
668	0	elseif (b(k)*vDml_slow(i) < 0.0) then
669	0	if (-b(k)vDml_slow(i) > h_avail(i,j+1,k) + a(k)vDml(i)) &
670	0	vDml_slow(i) = -max( 0., h_avail(i,j+1,k) + a(k)*vDml(i) ) / b(k)
671		endif
672		enddo
673	0	do k=1,nz
674	0	vhml(i,J,k) = a(k)vDml(i) + b(k)vDml_slow(i)
675	0	vhtr(i,J,k) = vhtr(i,J,k) + vhml(i,J,k)*dt
676		enddo
677		endif
678
679	0	vtimescale_diag(i,J) = timescale
680	0	vDml_diag(i,J) = vDml(i)
681		enddo ; enddo
682
683		!$OMP do
684	0	do j=js,je ; do k=1,nz ; do i=is,ie
685		h(i,j,k) = h(i,j,k) - dtG%IareaT(i,j) &
686	0	((uhml(I,j,k) - uhml(I-1,j,k)) + (vhml(i,J,k) - vhml(i,J-1,k)))
687	0	if (h(i,j,k) < h_min) h(i,j,k) = h_min
688		enddo ; enddo ; enddo
689		!$OMP end parallel
690
691		! Whenever thickness changes let the diag manager know, target grids
692		! for vertical remapping may need to be regenerated.
693	0	if (CS%id_uhml > 0 .or. CS%id_vhml > 0) &
694		! Remapped uhml and vhml require east/north halo updates of h
695	0	call pass_var(h, G%domain, To_West+To_South+Omit_Corners, halo=1)
696	0	call diag_update_remap_grids(CS%diag)
697
698		! Offer diagnostic fields for averaging.
699	0	if (query_averaging_enabled(CS%diag)) then
700	0	if (CS%id_urestrat_time > 0) call post_data(CS%id_urestrat_time, utimescale_diag, CS%diag)
701	0	if (CS%id_vrestrat_time > 0) call post_data(CS%id_vrestrat_time, vtimescale_diag, CS%diag)
702	0	if (CS%id_uhml > 0) call post_data(CS%id_uhml, uhml, CS%diag)
703	0	if (CS%id_vhml > 0) call post_data(CS%id_vhml, vhml, CS%diag)
704	0	if (CS%id_BLD > 0) call post_data(CS%id_BLD, MLD_fast, CS%diag)
705	0	if (CS%id_MLD > 0) call post_data(CS%id_MLD, MLD_slow, CS%diag)
706	0	if (CS%id_Rml > 0) call post_data(CS%id_Rml, Rml_av_fast, CS%diag)
707	0	if (CS%id_uDml > 0) call post_data(CS%id_uDml, uDml_diag, CS%diag)
708	0	if (CS%id_vDml > 0) call post_data(CS%id_vDml, vDml_diag, CS%diag)
709	0	if (CS%id_mle_fl > 0) call post_data(CS%id_mle_fl, mle_fl_2d, CS%diag)
710
711	0	if (CS%id_uml > 0) then
712	0	do j=js,je ; do I=is-1,ie
713	0	h_vel = 0.5*((htot_fast(i,j) + htot_fast(i+1,j)) + h_neglect)
714	0	uDml_diag(I,j) = uDml_diag(I,j) / (0.01h_vel) G%IdyCu(I,j) * (mu(0.,0.)-mu(-.01,0.))
715		enddo ; enddo
716	0	call post_data(CS%id_uml, uDml_diag, CS%diag)
717		endif
718	0	if (CS%id_vml > 0) then
719	0	do J=js-1,je ; do i=is,ie
720	0	h_vel = 0.5*((htot_fast(i,j) + htot_fast(i,j+1)) + h_neglect)
721	0	vDml_diag(i,J) = vDml_diag(i,J) / (0.01h_vel) G%IdxCv(i,J) * (mu(0.,0.)-mu(-.01,0.))
722		enddo ; enddo
723	0	call post_data(CS%id_vml, vDml_diag, CS%diag)
724		endif
725		endif
726		! Whenever thickness changes let the diag manager know, target grids
727		! for vertical remapping may need to be regenerated.
728		! This needs to happen after the H update and before the next post_data.
729	0	call diag_update_remap_grids(CS%diag)
730
731	0	end subroutine mixedlayer_restrat_OM4
732
733		!> Stream function shape as a function of non-dimensional position within mixed-layer [nondim]
734	26244000	pure real function mu(sigma, dh)
735		real, intent(in) :: sigma !< Fractional position within mixed layer [nondim]
736		!! z=0 is surface, z=-1 is the bottom of the mixed layer
737		real, intent(in) :: dh !< Non-dimensional distance over which to extend stream
738		!! function to smooth transport at base [nondim]
739		! Local variables
740		real :: xp !< A linear function from mid-point of the mixed-layer
741		!! to the extended mixed-layer bottom [nondim]
742		real :: bottop !< A mask, 0 in upper half of mixed layer, 1 otherwise [nondim]
743		real :: dd !< A cubic(-ish) profile in lower half of extended mixed
744		!! layer to smooth out the parameterized transport [nondim]
745
746		! Lower order shape (not used), see eq 10 from FK08b.
747		! Apparently used in CM2G, see eq 14 of FK11.
748		!mu = max(0., (1. - (2.sigma + 1.)*2))
749
750		! Second order, in Rossby number, shape. See eq 21 from FK08a, eq 9 from FK08b, eq 5 FK11
751	26244000	mu = max(0., (1. - (2.sigma + 1.)2) (1. + (5./21.)(2.sigma + 1.)**2))
752
753		! -0.5 < sigma : xp(sigma)=0 (upper half of mixed layer)
754		! -1.0+dh < sigma < -0.5 : xp(sigma)=linear (lower half +dh of mixed layer)
755		! sigma < -1.0+dh : xp(sigma)=1 (below mixed layer + dh)
756	26244000	xp = max(0., min(1., (-sigma - 0.5)2. / (1. + 2.dh)))
757
758		! -0.5 < sigma : dd(sigma)=1 (upper half of mixed layer)
759		! -1.0+dh < sigma < -0.5 : dd(sigma)=cubic (lower half +dh of mixed layer)
760		! sigma < -1.0+dh : dd(sigma)=0 (below mixed layer + dh)
761	26244000	dd = (max(1. - xp*2 (3. - 2.xp), 0.))(1. + 2.dh)
762
763		! -0.5 < sigma : bottop(sigma)=0 (upper half of mixed layer)
764		! sigma < -0.5 : bottop(sigma)=1 (below upper half)
765	26244000	bottop = 0.5*(1. - sign(1., sigma + 0.5)) ! =0 for sigma>-0.5, =1 for sigma<-0.5
766
767	26244000	mu = max(mu, dd*bottop) ! Combines original psi1 with tail
768	26244000	end function mu
769
770		!> Calculates a restratifying flow in the mixed layer, following the formulation
771		!! used in Bodner et al., 2023 (B22)
772	24	subroutine mixedlayer_restrat_Bodner(CS, G, GV, US, h, uhtr, vhtr, tv, forces, dt, BLD, h_MLD, bflux)
773		! Arguments
774		type(mixedlayer_restrat_CS), intent(inout) :: CS !< Module control structure
775		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
776		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
777		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
778		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
779		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: uhtr !< Accumulated zonal mass flux
780		!! [H L2 ~> m3 or kg]
781		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: vhtr !< Accumulated meridional mass flux
782		!! [H L2 ~> m3 or kg]
783		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables structure
784		type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
785		real, intent(in) :: dt !< Time increment [T ~> s]
786		real, dimension(:,:), pointer :: BLD !< Active boundary layer depth provided by the
787		!! PBL scheme [Z ~> m] (not H)
788		real, dimension(:,:), pointer :: h_MLD !< Thickness of water within the
789		!! active boundary layer depth provided by
790		!! the PBL scheme [H ~> m or kg m-2]
791		real, dimension(:,:), pointer :: bflux !< Surface buoyancy flux provided by the
792		!! PBL scheme [Z2 T-3 ~> m2 s-3]
793		! Local variables
794	48	real :: uhml(SZIB_(G),SZJ_(G),SZK_(GV)) ! zonal mixed layer transport [H L2 T-1 ~> m3 s-1 or kg s-1]
795	48	real :: vhml(SZI_(G),SZJB_(G),SZK_(GV)) ! merid mixed layer transport [H L2 T-1 ~> m3 s-1 or kg s-1]
796	48	real :: vol_dt_avail(SZI_(G),SZJ_(G),SZK_(GV)) ! The volume available for exchange out of each face of
797		! each layer, divided by dt [H L2 T-1 ~> m3 s-1 or kg s-1]
798		real, dimension(SZI_(G),SZJ_(G)) :: &
799	48	little_h, & ! "Little h" representing active mixing layer depth [H ~> m or kg m-2]
800	48	big_H, & ! "Big H" representing the mixed layer depth [H ~> m or kg m-2]
801	48	mld, & ! The mixed layer depth returned by detect_mld [H ~> m or kg m-2]
802	48	htot, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2]
803	48	buoy_av, & ! g_Rho0 times the average mixed layer density or G_Earth
804		! times the average specific volume [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2]
805	48	wpup ! Turbulent vertical momentum [L H T-2 ~> m2 s-2 or kg m-1 s-2]
806	48	real :: uDml_diag(SZIB_(G),SZJ_(G)) ! A 2D copy of uDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1]
807	48	real :: vDml_diag(SZI_(G),SZJB_(G)) ! A 2D copy of vDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1]
808	48	real :: lf_bodner_diag(SZI_(G),SZJ_(G)) ! Front width as in Bodner et al., 2023 (B22), eq 24 [L ~> m]
809	48	real :: U_star_2d(SZI_(G),SZJ_(G)) ! The wind friction velocity, calculated using the Boussinesq
810		! reference density or the time-evolving surface density in non-Boussinesq
811		! mode [Z T-1 ~> m s-1]
812	48	real :: covTS(SZI_(G)) ! SGS TS covariance in Stanley param; currently 0 [C S ~> degC ppt]
813	48	real :: varS(SZI_(G)) ! SGS S variance in Stanley param; currently 0 [S2 ~> ppt2]
814	48	real :: Rml_int(SZI_(G)) ! Potential density integrated through the mixed layer [R H ~> kg m-2 or kg2 m-5]
815	48	real :: SpV_ml(SZI_(G)) ! Specific volume evaluated at the surface pressure [R-1 ~> m3 kg-1]
816	48	real :: SpV_int(SZI_(G)) ! Specific volume integrated through the mixed layer [H R-1 ~> m4 kg-1 or m]
817	48	real :: rho_ml(SZI_(G)) ! Potential density relative to the surface [R ~> kg m-3]
818		real :: rho_blk(SZI_(G), merge(G%jed-G%jsd+1, CS%njblock, CS%njblock==0), &
819	24	merge(GV%ke, CS%nkblock, CS%nkblock==0))
820		! Potential density relative to the surface for a tile [R ~> kg m-3]
821		real :: p_blk(SZI_(G), merge(G%jed-G%jsd+1, CS%njblock, CS%njblock==0), &
822	24	merge(GV%ke, CS%nkblock, CS%nkblock==0))
823		! A pressure of 0 for a tile [R L2 T-2 ~> Pa]
824	24	real :: Rml_blk(SZI_(G), merge(G%jed-G%jsd+1, CS%njblock, CS%njblock==0))
825		! Density integrated through the mixed layer, for the rows of a tile,
826		! carried across the k blocks of a column [R H ~> kg m-2 or kg2 m-5]
827	24	real :: p0(SZI_(G)) ! A pressure of 0 [R L2 T-2 ~> Pa]
828		real :: g_Rho0 ! G_Earth/Rho0 times a thickness conversion factor
829		! [L2 H-1 T-2 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2]
830		real :: h_vel ! htot interpolated onto velocity points [H ~> m or kg m-2]
831		real :: w_star3 ! Cube of turbulent convective velocity [Z3 T-3 ~> m3 s-3]
832		real :: u_star3 ! Cube of surface friction velocity [Z3 T-3 ~> m3 s-3]
833		real :: r_wpup ! reciprocal of vertical momentum flux [T2 L-1 H-1 ~> s2 m-2 or m s2 kg-1]
834		real :: absf ! absolute value of f, interpolated to velocity points [T-1 ~> s-1]
835		real :: f_h ! Coriolis parameter at h-points [T-1 ~> s-1]
836		real :: f2_h ! Coriolis parameter at h-points squared [T-2 ~> s-2]
837		real :: absurdly_small_freq2 ! Frequency squared used to avoid division by 0 [T-2 ~> s-2]
838		real :: grid_dsd ! combination of grid scales [L2 ~> m2]
839		real :: h_sml ! "Little h", the active mixing depth with diurnal cycle removed [H ~> m or kg m-2]
840		real :: h_big ! "Big H", the mixed layer depth based on a time filtered "little h" [H ~> m or kg m-2]
841		real :: grd_b ! The vertically average gradient of buoyancy [L H-1 T-2 ~> s-2 or m3 kg-1 s-2]
842		real :: psi_mag ! Magnitude of stream function [L2 H T-1 ~> m3 s-1 or kg s-1]
843		real :: h_neglect ! tiny thickness usually lost in roundoff so can be neglected [H ~> m or kg m-2]
844		real :: I4dt ! 1/(4 dt) [T-1 ~> s-1]
845		real :: Ihtot ! Inverses of the total mixed layer thickness [H-1 ~> m-1 or m2 kg-1]
846		real :: hAtVel ! Thickness at the velocity points [H ~> m or kg m-2]
847		real :: sigint ! Fractional position within the mixed layer of the interface above a layer [nondim]
848		real :: muzb ! mu(z) at bottom of the layer [nondim]
849		real :: muza ! mu(z) at top of the layer [nondim]
850		real :: dh ! Portion of the layer thickness that is in the mixed layer [H ~> m or kg m-2]
851		real :: Z3_T3_to_m3_s3 ! Conversion factors to undo scaling and permit terms to be raised to a
852		! fractional power [T3 m3 Z-3 s-3 ~> 1]
853		real :: m2_s2_to_Z2_T2 ! Conversion factors to restore scaling after a term is raised to a
854		! fractional power [Z2 s2 T-2 m-2 ~> 1]
855		real, parameter :: two_thirds = 2./3. ! [nondim]
856		real :: dmu ! Change in mu(z) across a layer [nondim]
857		logical :: line_is_empty, keep_going
858		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
859		integer :: EOSdom3(3,2) ! The (i,j,k) computational domain for the blocked equation of state calls
860		integer :: nkblock, njblock ! The number of layers and rows in each tile of the density integral [nondim]
861		integer :: kstart, kend ! The first and last layer of the tile being worked on [nondim]
862		integer :: jstart, jend ! The first and last row of the tile being worked on [nondim]
863		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
864
865	24	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
866	24	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
867
868	24	I4dt = 0.25 / dt
869	24	g_Rho0 = GV%H_to_Z * GV%g_Earth / GV%Rho0
870	24	h_neglect = GV%H_subroundoff
871
872	3096	covTS(:) = 0.0 ! Might be in tv% in the future. Not implemented for the time being.
873	3096	varS(:) = 0.0 ! Ditto.
874
875		! This value is roughly (pi / (the age of the universe) )^2.
876	24	absurdly_small_freq2 = 1e-34US%T_to_s*2
877
878	24	if (.not.associated(tv%eqn_of_state)) call MOM_error(FATAL, "mixedlayer_restrat_Bodner: "// &
879	0	"An equation of state must be used with this module.")
880	24	if (CS%MLE_use_PBL_MLD) then
881	24	if (CS%MLE_density_diff > 0.) call MOM_error(FATAL, "mixedlayer_restrat_Bodner: "// &
882	0	"MLE_density_diff is +ve and should not be in mixedlayer_restrat_Bodner.")
883	24	if (.not.associated(bflux)) call MOM_error(FATAL, "mixedlayer_restrat_Bodner: "// &
884	0	"Surface buoyancy flux was not associated.")
885		else
886	0	if (.not.CS%Bodner_detect_MLD) call MOM_error(FATAL, "mixedlayer_restrat_Bodner: "// &
887		"To use the Bodner et al., 2023, MLE parameterization, either MLE_USE_PBL_MLD or "// &
888	0	"Bodner_detect_MLD must be True.")
889		endif
890
891	24	if (associated(bflux)) &
892	24	call pass_var(bflux, G%domain, halo=1)
893
894		! Extract the friction velocity from the forcing type.
895	24	call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1)
896
897	24	if (CS%debug) then
898	0	call hchksum(h,'mixed_Bodner: h', G%HI, haloshift=1, unscale=GV%H_to_mks)
899	0	call hchksum(BLD, 'mle_Bodner: BLD', G%HI, haloshift=1, unscale=US%Z_to_m)
900	0	call hchksum(h_MLD, 'mle_Bodner: h_MLD', G%HI, haloshift=1, unscale=GV%H_to_mks)
901	0	if (associated(bflux)) &
902	0	call hchksum(bflux, 'mle_Bodner: bflux', G%HI, haloshift=1, unscale=US%Z_to_m*2US%s_to_T**3)
903	0	call hchksum(U_star_2d, 'mle_Bodner: u', G%HI, haloshift=1, unscale=US%Z_to_mUS%s_to_T)
904		call hchksum(CS%MLD_filtered, 'mle_Bodner: MLD_filtered 1', &
905	0	G%HI, haloshift=1, unscale=GV%H_to_mks)
906		call hchksum(CS%MLD_filtered_slow,'mle_Bodner: MLD_filtered_slow 1', &
907	0	G%HI, haloshift=1, unscale=GV%H_to_mks)
908		endif
909
910		!$omp target enter data map(to: U_star_2d, h_MLD)
911		!$omp target enter data map(alloc: little_h, big_H, wpup, htot, buoy_av, uDml_diag, vDml_diag)
912		!$omp target enter data map(alloc: vol_dt_avail, uhml, vhml)
913
914		! Apply time filter to h_MLD (to remove diurnal cycle) to obtain "little h".
915		! "little h" is representative of the active mixing layer depth, used in B22 formula (eq 27).
916	2952	do concurrent (j=js-1:je+1, i=is-1:ie+1)
917		little_h(i,j) = rmean2ts(h_MLD(i,j), CS%MLD_filtered(i,j), &
918	181536	CS%BLD_growing_Tfilt, CS%BLD_decaying_Tfilt, dt)
919	184488	CS%MLD_filtered(i,j) = little_h(i,j)
920		enddo
921
922		! Calculate "big H", representative of the mixed layer depth, used in B22 formula (eq 27).
923		! not offloaded, send back to the host
924	24	if (CS%MLD_grid) then
925		!$omp target update from(little_h, CS%MLD_filtered_slow)
926	0	do j=js-1,je+1 ; do i=is-1,ie+1
927		big_H(i,j) = rmean2ts(little_h(i,j), CS%MLD_filtered_slow(i,j), &
928	0	CS%MLD_growing_Tfilt, CS%MLD_Tfilt_space(i,j), dt)
929		enddo ; enddo
930		!$omp target update to(big_H)
931	24	elseif (CS%Bodner_detect_MLD) then
932		!$omp target update from(CS%MLD_filtered_slow)
933	0	call detect_mld(h, tv, MLD, G, GV, CS)
934	0	do j=js-1,je+1 ; do i=is-1,ie+1
935		big_H(i,j) = rmean2ts(MLD(i,j), CS%MLD_filtered_slow(i,j), &
936	0	CS%MLD_growing_Tfilt, CS%MLD_decaying_Tfilt, dt)
937		enddo ; enddo
938		!$omp target update to(big_H)
939		else
940	2952	do concurrent (j=js-1:je+1, i=is-1:ie+1)
941		big_H(i,j) = rmean2ts(little_h(i,j), CS%MLD_filtered_slow(i,j), &
942	184488	CS%MLD_growing_Tfilt, CS%MLD_decaying_Tfilt, dt)
943		enddo
944		endif
945	2952	do concurrent (j=js-1:je+1, i=is-1:ie+1)
946	184488	CS%MLD_filtered_slow(i,j) = big_H(i,j)
947		enddo
948
949		! Estimate w'u' at h-points, with a floor to avoid division by zero later.
950		! needed here, otehrwe death
951		!$omp target enter data map(to: BLD, bflux)
952	24	if (allocated(tv%SpV_avg) .and. .not.(GV%Boussinesq .or. GV%semi_Boussinesq)) then
953	0	do j=js-1,je+1 ; do i=is-1,ie+1
954		! This expression differs by a factor of 1. / (Rho_0 * SpV_avg) compared with the other
955		! expressions below, and it is invariant to the value of Rho_0 in non-Boussinesq mode.
956		wpup(i,j) = max((cuberoot( CS%mstar * U_star_2d(i,j)**3 + &
957		CS%nstar * max(0., -bflux(i,j)) * BLD(i,j) ))**2, CS%min_wstar2) &
958	0	* (US%Z_to_L * GV%RZ_to_H / tv%SpV_avg(i,j,1))
959		! The final line above converts from [Z2 T-2 ~> m2 s-2] to [L H T-2 ~> m2 s-2 or Pa].
960		! Some rescaling factors and the division by specific volume compensating for other
961		! factors that are in find_ustar_mech, and others effectively converting the wind
962		! stresses from [R L Z T-2 ~> Pa] to [L H T-2 ~> m2 s-2 or Pa]. The rescaling factors
963		! and density being applied to the buoyancy flux are not so neatly explained because
964		! fractional powers cancel out or combine with terms in the definitions of BLD and
965		! bflux (such as SpV_avg**-2/3 combining with other terms in bflux to give the thermal
966		! expansion coefficient) and because the specific volume does vary within the mixed layer.
967		enddo ; enddo
968		!$omp target update to(wpup)
969	24	elseif (CS%answer_date < 20240201) then
970	0	Z3_T3_to_m3_s3 = (US%Z_to_m * US%s_to_T)**3
971	0	m2_s2_to_Z2_T2 = (US%m_to_Z * US%T_to_s)**2
972	0	do j=js-1,je+1 ; do i=is-1,ie+1
973	0	w_star3 = max(0., -bflux(i,j)) * BLD(i,j) ! In [Z3 T-3 ~> m3 s-3]
974	0	u_star3 = U_star_2d(i,j)**3 ! In [Z3 T-3 ~> m3 s-3]
975		wpup(i,j) = max(m2_s2_to_Z2_T2 * (Z3_T3_to_m3_s3 * ( CS%mstar * u_star3 + CS%nstar * w_star3 ) )**two_thirds, &
976	0	CS%min_wstar2) * US%Z_to_L * GV%Z_to_H ! In [L H T-2 ~> m2 s-2 or kg m-1 s-2]
977		enddo ; enddo
978		!$omp target update to(wpup)
979		else
980	2952	do concurrent (j=js-1:je+1, i=is-1:ie+1)
981		wpup(i,j) = max( (cuberoot(CS%mstar * U_star_2d(i,j)**3 &
982		+ CS%nstar * (max(0., -bflux(i,j)) * BLD(i,j))))**2, CS%min_wstar2 ) &
983	184488	* US%Z_to_L * GV%Z_to_H ! In [L H T-2 ~> m2 s-2 or kg m-1 s-2]
984		enddo
985		endif
986		!$omp target exit data map(release: BLD, bflux)
987
988
989		! We filter w'u' with the same time scales used for "little h"
990	2952	do concurrent (j=js-1:je+1, i=is-1:ie+1)
991		wpup(i,j) = rmean2ts(wpup(i,j), CS%wpup_filtered(i,j), &
992	181536	CS%BLD_growing_Tfilt, CS%BLD_decaying_Tfilt, dt)
993	184488	CS%wpup_filtered(i,j) = wpup(i,j)
994		enddo
995
996	24	if (CS%id_lfbod > 0) then
997		! This diagnostic is evaluated on the host, so it needs the device's "little h".
998		!$omp target update from(little_h)
999	0	do j=js-1,je+1 ; do i=is-1,ie+1
1000		! Calculate front length used in B22 formula (eq 24).
1001	0	w_star3 = max(0., -bflux(i,j)) * BLD(i,j)
1002	0	u_star3 = U_star_2d(i,j)**3
1003
1004		! Include an absurdly_small_freq2 to prevent division by zero.
1005		f_h = 0.25 * ((G%CoriolisBu(I,J) + G%CoriolisBu(I-1,J-1)) &
1006	0	+ (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J-1)))
1007	0	f2_h = max(f_h**2, absurdly_small_freq2)
1008
1009		lf_bodner_diag(i,j) = &
1010		0.25 * cuberoot(CS%mstar * u_star3 + CS%nstar * w_star3)**2 &
1011	0	/ (f2_h * max(little_h(i,j), GV%Angstrom_H))
1012		enddo ; enddo
1013
1014		! Rescale from [Z2 H-1 ~> m or m4 kg-1] to [L ~> m]
1015	0	if (allocated(tv%SpV_avg) .and. .not.(GV%Boussinesq .or. GV%semi_Boussinesq)) then
1016	0	do j=js-1,je+1 ; do i=is-1,ie+1
1017		lf_bodner_diag(i,j) = lf_bodner_diag(i,j) &
1018	0	* (US%Z_to_L * GV%RZ_to_H / tv%SpV_avg(i,j,1))
1019		enddo ; enddo
1020		else
1021	0	do j=js-1,je+1 ; do i=is-1,ie+1
1022	0	lf_bodner_diag(i,j) = lf_bodner_diag(i,j) * US%Z_to_L * GV%Z_to_H
1023		enddo ; enddo
1024		endif
1025		endif
1026
1027	24	if (CS%debug) then
1028		!$omp target update from(little_h, big_H, wpup, CS%MLD_filtered, CS%MLD_filtered_slow)
1029	0	call hchksum(little_h,'mle_Bodner: little_h', G%HI, haloshift=1, unscale=GV%H_to_mks)
1030	0	call hchksum(big_H,'mle_Bodner: big_H', G%HI, haloshift=1, unscale=GV%H_to_mks)
1031		call hchksum(CS%MLD_filtered,'mle_Bodner: MLD_filtered 2', &
1032	0	G%HI, haloshift=1, unscale=GV%H_to_mks)
1033		call hchksum(CS%MLD_filtered_slow,'mle_Bodner: MLD_filtered_slow 2', &
1034	0	G%HI, haloshift=1, unscale=GV%H_to_mks)
1035	0	call hchksum(wpup,'mle_Bodner: wpup', G%HI, haloshift=1, unscale=US%L_to_mGV%H_to_mksUS%s_to_T**2)
1036		endif
1037
1038		! Calculate the average density in the "mixed layer".
1039		! Notice we use p=0 (sigma_0) since horizontal differences of vertical averages of
1040		! in-situ density would contain the MLD gradient (through the pressure dependence).
1041	3096	p0(:) = 0.0
1042	24	EOSdom(:) = EOS_domain(G%HI, halo=1)
1043
1044	2952	do concurrent (k=1:nz, j=js-1:je+1, i=is-1:ie+1)
1045	13799688	vol_dt_avail(i,j,k) = max(I4dtG%areaT(i,j)(h(i,j,k)-GV%Angstrom_H), 0.0)
1046		enddo
1047
1048	24	if ((GV%Boussinesq .or. GV%semi_Boussinesq) .and. .not.CS%use_Stanley_ML) then
1049		! Walk the mixed layer in tiles of njblock rows by nkblock layers: evaluate the density for a
1050		! tile, then integrate it. The columns are independent of one another, so the integral runs as
1051		! one concurrent loop over (i,j) with the sequential k-dependence kept inside each column, and
1052		! the running totals for the tile's rows carry across its k blocks.
1053	24	nkblock = merge(nz, CS%nkblock, CS%nkblock==0)
1054	24	njblock = merge((je+1) - (js-1) + 1, CS%njblock, CS%njblock==0)
1055	72	EOSdom3(1,:) = EOS_domain(G%HI, halo=1)
1056
1057		! The tile scratch never leaves the device, and no other branch uses it.
1058		!$omp target enter data map(alloc: rho_blk, p_blk, Rml_blk)
1059
1060		if (nkblock < nz) then
1061		! The early-exit test below compares htot against big_H on the host, and big_H has been
1062		! computed on the device, so fetch it once. With a single k block that test never runs.
1063		!$omp target update from(big_H)
1064		endif
1065
1066	24	do concurrent (k=1:nkblock, j=1:njblock, i=is-1:ie+1)
1067	8808	p_blk(i,j,k) = 0.0
1068		enddo
1069
1070	1512	do jstart=js-1,je+1,njblock
1071	1488	jend = min(jstart+njblock-1, je+1)
1072	4464	EOSdom3(2,:) = [1, jend-jstart+1]
1073
1074	1488	do concurrent (j=jstart:jend, i=is-1:ie+1)
1075	364560	htot(i,j) = 0.0 ; Rml_blk(i,j-jstart+1) = 0.0
1076		enddo
1077
1078	1488	keep_going = .true.
1079	156953	do kstart=1,nz,nkblock ; if (keep_going) then
1080	43865	kend = min(kstart+nkblock-1, nz)
1081	131595	EOSdom3(3,:) = [1, kend-kstart+1]
1082		call calculate_density(tv%T(:,jstart:jend,kstart:kend), tv%S(:,jstart:jend,kstart:kend), &
1083	43865	p_blk, rho_blk, tv%eqn_of_state, EOSdom3)
1084
1085	43865	do concurrent (j=jstart:jend, i=is-1:ie+1) DO_LOCALITY(local(k, dh))
1086	16098455	do k=kstart,kend
1087	10703060	if (htot(i,j) < big_H(i,j)) then
1088	1228014	dh = min( h(i,j,k), big_H(i,j) - htot(i,j) )
1089		! Rml_blk is in [R H ~> kg m-2]
1090	1228014	Rml_blk(i,j-jstart+1) = Rml_blk(i,j-jstart+1) + dh*rho_blk(i,j-jstart+1,k-kstart+1)
1091	1228014	htot(i,j) = htot(i,j) + dh
1092		endif
1093		enddo
1094		enddo
1095
1096	43865	if (nkblock < nz) then
1097		! Stop calling the equation of state once every column in the tile has been filled to
1098		! "big H". With a single k block there is nothing left to skip, and this test would drag
1099		! htot back to the host for nothing.
1100		!$omp target update from(htot)
1101	43865	keep_going = .false.
1102	5439260	do j=jstart,jend ; do i=is-1,ie+1
1103	5395395	if (htot(i,j) < big_H(i,j)) keep_going = .true.
1104		enddo ; enddo
1105		endif
1106		endif ; enddo
1107
1108	24	do concurrent (j=jstart:jend, i=is-1:ie+1)
1109		! Buoy_av has units (L2 H-1 T-2 R-1) * (R H) * H-1 = [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2]
1110	364560	buoy_av(i,j) = -( g_Rho0 * Rml_blk(i,j-jstart+1) ) / (htot(i,j) + h_neglect)
1111		enddo
1112		enddo
1113
1114		!$omp target exit data map(release: rho_blk, p_blk, Rml_blk)
1115		else
1116		! The Stanley and non-Boussinesq variants of the equation of state have no array-of-columns
1117		! interface to block over, so they retain the original row-at-a-time form on the host.
1118		!$OMP parallel do default(shared) &
1119		!$OMP private(i, k, keep_going, line_is_empty, dh, rho_ml, SpV_ml, Rml_int, SpV_int)
1120	0	do j=js-1,je+1
1121	0	rho_ml(:) = 0.0 ; SpV_ml(:) = 0.0
1122	0	do i=is-1,ie+1
1123	0	htot(i,j) = 0.0 ; Rml_int(i) = 0.0 ; SpV_int(i) = 0.0
1124		enddo
1125	0	keep_going = .true.
1126	0	do k=1,nz
1127	0	if (keep_going) then
1128	0	if (GV%Boussinesq .or. GV%semi_Boussinesq) then
1129		call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, tv%varT(:,j,k), covTS, varS, &
1130	0	rho_ml, tv%eqn_of_state, EOSdom)
1131		else
1132	0	call calculate_spec_vol(tv%T(:,j,k), tv%S(:,j,k), p0, SpV_ml, tv%eqn_of_state, EOSdom)
1133		endif
1134	0	line_is_empty = .true.
1135	0	do i=is-1,ie+1
1136	0	if (htot(i,j) < big_H(i,j)) then
1137	0	dh = min( h(i,j,k), big_H(i,j) - htot(i,j) )
1138	0	Rml_int(i) = Rml_int(i) + dh*rho_ml(i) ! Rml_int has units of [R H ~> kg m-2]
1139	0	SpV_int(i) = SpV_int(i) + dh*SpV_ml(i) ! SpV_int has units of [H R-1 ~> m4 kg-1 or m]
1140	0	htot(i,j) = htot(i,j) + dh
1141	0	line_is_empty = .false.
1142		endif
1143		enddo
1144	0	if (line_is_empty) keep_going=.false.
1145		endif
1146		enddo
1147
1148	0	if (GV%Boussinesq .or. GV%semi_Boussinesq) then
1149	0	do i=is-1,ie+1
1150		! Buoy_av has units (L2 H-1 T-2 R-1) * (R H) * H-1 = [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2]
1151	0	buoy_av(i,j) = -( g_Rho0 * Rml_int(i) ) / (htot(i,j) + h_neglect)
1152		enddo
1153		else
1154	0	do i=is-1,ie+1
1155		! Buoy_av has units (R L2 H-1 T-2) * (R-1 H) * H-1 = [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2]
1156	0	buoy_av(i,j) = (GV%H_to_RZGV%g_Earth SpV_int(i)) / (htot(i,j) + h_neglect)
1157		enddo
1158		endif
1159		enddo
1160		!$omp target update to(htot, buoy_av)
1161		endif
1162
1163	24	if (CS%debug) then
1164		!$omp target update from(htot, vol_dt_avail, buoy_av)
1165	0	call hchksum(htot,'mle_Bodner: htot', G%HI, haloshift=1, unscale=GV%H_to_mks)
1166		call hchksum(vol_dt_avail,'mle_Bodner: vol_dt_avail', G%HI, haloshift=1, &
1167	0	unscale=US%L_to_m*2GV%H_to_mks*US%s_to_T)
1168	0	call hchksum(buoy_av,'mle_Bodner: buoy_av', G%HI, haloshift=1, unscale=GV%m_to_HUS%L_T_to_m_s*2)
1169		endif
1170
1171		! U - Component
1172		do concurrent (j=js:je, I=is-1:ie) &
1173	24	DO_LOCALITY(local(k,dmu,grid_dsd,absf,h_sml,h_big,grd_b,r_wpup,psi_mag,IhTot,sigint,muzb,muza,hAtVel))
1174	174240	if (G%OBCmaskCu(I,j) > 0.) then
1175	118752	grid_dsd = sqrt(0.5( G%dxCu(I,j)2 + G%dyCu(I,j)2 )) G%dyCu(I,j) ! [L2 ~> m2]
1176	118752	absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J))) ! [T-1 ~> s-1]
1177	118752	h_sml = 0.5*( little_h(i,j) + little_h(i+1,j) ) ! [H ~> m or kg m-2]
1178	118752	h_big = 0.5*( big_H(i,j) + big_H(i+1,j) ) ! [H ~> m or kg m-2]
1179	118752	grd_b = ( buoy_av(i+1,j) - buoy_av(i,j) ) * G%IdxCu(I,j) ! [L H-1 T-2 ~> s-2 or m3 kg-1 s-2]
1180	118752	r_wpup = 2. / ( wpup(i,j) + wpup(i+1,j) ) ! [T2 L-1 H-1 ~> s2 m-2 or m s2 kg-1]
1181		psi_mag = ( ( ( (0.5(CS%Cr_space(i,j) + CS%Cr_space(i+1,j))) grid_dsd ) & ! [L2 H T-1 ~> m3 s-1 or kg s-1]
1182	118752	* ( absf * h_sml ) ) * ( ( h_big*2 ) grd_b ) ) * r_wpup
1183		else ! There is no flux on land and no gradient at open boundary points.
1184	55488	psi_mag = 0.0
1185		endif
1186
1187	174240	IhTot = 2.0 / ((htot(i,j) + htot(i+1,j)) + h_neglect) ! [H-1 ~> m-1 or m2 kg-1]
1188	174240	sigint = 0.0
1189	174240	muzb = 0.0 ! This will be the first value of muza = mu(z=0)
1190	13242240	do k=1,nz
1191	13068000	muza = muzb ! mu(z/MLD) for upper interface [nondim]
1192	13068000	hAtVel = 0.5*(h(i,j,k) + h(i+1,j,k)) ! Thickness at velocity point [H ~> m or kg m-2]
1193	13068000	sigint = sigint - (hAtVel * IhTot) ! z/H for lower interface [nondim]
1194	13068000	muzb = mu(sigint, CS%MLE_tail_dh) ! mu(z/MLD) for lower interface [nondim]
1195	13068000	dmu = muza - muzb ! Change in mu(z) across layer [nondim]
1196	13068000	uhml(I,j,k) = dmu ! Stash dmu in uhml: the columns run concurrently, so there is nowhere
1197		! else to keep a per-column profile. It is scaled by psi_mag below.
1198		! dmu*psi_mag is the transport in this layer [L2 H T-1 ~> m3 s-1]
1199		! Limit magnitude (psi_mag) if it would violate CFL
1200	13242240	if (dmu*psi_mag > 0.0) then
1201	554486	if (dmu*psi_mag > vol_dt_avail(i,j,k)) psi_mag = vol_dt_avail(i,j,k) / dmu
1202	12513514	elseif (dmu*psi_mag < 0.0) then
1203	567746	if (-dmu*psi_mag > vol_dt_avail(i+1,j,k)) psi_mag = -vol_dt_avail(i+1,j,k) / dmu
1204		endif
1205		enddo ! These loops cannot be fused because psi_mag applies to the whole column
1206	13242240	do k=1,nz
1207	13068000	uhml(I,j,k) = uhml(I,j,k) * psi_mag ! [L2 H T-1 ~> m3 s-1 or kg s-1]
1208	13242240	uhtr(I,j,k) = uhtr(I,j,k) + uhml(I,j,k) * dt ! [L2 H ~> m3 or kg]
1209		enddo
1210
1211	177168	uDml_diag(I,j) = psi_mag
1212		enddo
1213
1214		! V- component
1215		do concurrent (J=js-1:je, i=is:ie) &
1216	24	DO_LOCALITY(local(k,dmu,grid_dsd,absf,h_sml,h_big,grd_b,r_wpup,psi_mag,IhTot,sigint,muzb,muza,hAtVel))
1217	175680	if (G%OBCmaskCv(i,J) > 0.) then
1218	117456	grid_dsd = sqrt(0.5( G%dxCv(i,J)2 + G%dyCv(i,J)2 )) G%dxCv(i,J) ! [L2 ~> m2]
1219	117456	absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J))) ! [T-1 ~> s-1]
1220	117456	h_sml = 0.5*( little_h(i,j) + little_h(i,j+1) ) ! [H ~> m or kg m-2]
1221	117456	h_big = 0.5*( big_H(i,j) + big_H(i,j+1) ) ! [H ~> m or kg m-2]
1222	117456	grd_b = ( buoy_av(i,j+1) - buoy_av(i,j) ) * G%IdyCv(I,j) ! [L H-1 T-2 ~> s-2 or m3 kg-1 s-2]
1223	117456	r_wpup = 2. / ( wpup(i,j) + wpup(i,j+1) ) ! [T2 L-1 H-1 ~> s2 m-2 or m s2 kg-1]
1224		psi_mag = ( ( ( (0.5(CS%Cr_space(i,j) + CS%Cr_space(i,j+1))) grid_dsd ) & ! [L2 H T-1 ~> m3 s-1 or kg s-1]
1225	117456	* ( absf * h_sml ) ) * ( ( h_big*2 ) grd_b ) ) * r_wpup
1226		else ! There is no flux on land and no gradient at open boundary points.
1227	58224	psi_mag = 0.0
1228		endif
1229
1230	175680	IhTot = 2.0 / ((htot(i,j) + htot(i,j+1)) + h_neglect) ! [H-1 ~> m-1 or m2 kg-1]
1231	175680	sigint = 0.0
1232	175680	muzb = 0.0 ! This will be the first value of muza = mu(z=0)
1233	13351680	do k=1,nz
1234	13176000	muza = muzb ! mu(z/MLD) for upper interface [nondim]
1235	13176000	hAtVel = 0.5*(h(i,j,k) + h(i,j+1,k)) ! Thickness at velocity point [H ~> m or kg m-2]
1236	13176000	sigint = sigint - (hAtVel * IhTot) ! z/H for lower interface [nondim]
1237	13176000	muzb = mu(sigint, CS%MLE_tail_dh) ! mu(z/MLD) for lower interface [nondim]
1238	13176000	dmu = muza - muzb ! Change in mu(z) across layer [nondim]
1239	13176000	vhml(i,J,k) = dmu ! Stash dmu in vhml, as for uhml above; scaled by psi_mag below.
1240		! dmu*psi_mag is the transport in this layer [L2 H T-1 ~> m3 s-1 or kg s-1]
1241		! Limit magnitude (psi_mag) if it would violate CFL
1242	13351680	if (dmu*psi_mag > 0.0) then
1243	493134	if (dmu*psi_mag > vol_dt_avail(i,j,k)) psi_mag = vol_dt_avail(i,j,k) / dmu
1244	12682866	elseif (dmu*psi_mag < 0.0) then
1245	586332	if (-dmu*psi_mag > vol_dt_avail(i,j+1,k)) psi_mag = -vol_dt_avail(i,j+1,k) / dmu
1246		endif
1247		enddo ! These loops cannot be fused because psi_mag applies to the whole column
1248	13351680	do k=1,nz
1249	13176000	vhml(i,J,k) = vhml(i,J,k) * psi_mag ! [L2 H T-1 ~> m3 s-1 or kg s-1]
1250	13351680	vhtr(i,J,k) = vhtr(i,J,k) + vhml(i,J,k) * dt ! [L2 H ~> m3 or kg]
1251		enddo
1252
1253	178584	vDml_diag(i,J) = psi_mag
1254		enddo
1255
1256	24	do concurrent (j=js:je, k=1:nz, i=is:ie)
1257		h(i,j,k) = h(i,j,k) - dtG%IareaT(i,j) &
1258	13178904	((uhml(I,j,k) - uhml(I-1,j,k)) + (vhml(i,J,k) - vhml(i,J-1,k)))
1259		enddo
1260
1261		! h, uhtr and vhtr stay mapped for the whole run, but step_MOM_dynamics pushes the host copies
1262		! back to the device after this call, so the values just computed here have to reach the host.
1263		!$omp target update from(h, uhtr, vhtr)
1264
1265		! uhml and vhml are otherwise device-only; retrieve them only when they are being diagnosed.
1266		if (CS%id_uhml > 0) then
1267		!$omp target update from(uhml)
1268		endif
1269		if (CS%id_vhml > 0) then
1270		!$omp target update from(vhml)
1271		endif
1272
1273		!$omp target exit data map(from: little_h, big_H, wpup, htot, buoy_av, uDml_diag, vDml_diag)
1274		!$omp target exit data map(release: vol_dt_avail, uhml, vhml, U_star_2d, h_MLD)
1275
1276	24	if (CS%id_uhml > 0 .or. CS%id_vhml > 0) &
1277		! Remapped uhml and vhml require east/north halo updates of h
1278	0	call pass_var(h, G%domain, To_West+To_South+Omit_Corners, halo=1)
1279		! Whenever thickness changes let the diag manager know, target grids
1280		! for vertical remapping may need to be regenerated.
1281	24	call diag_update_remap_grids(CS%diag)
1282
1283		! Offer diagnostic fields for averaging.
1284	24	if (query_averaging_enabled(CS%diag)) then
1285	24	if (CS%id_ustar > 0) call post_data(CS%id_ustar, U_star_2d, CS%diag)
1286	24	if (CS%id_bflux > 0) call post_data(CS%id_bflux, bflux, CS%diag)
1287	24	if (CS%id_wpup > 0) call post_data(CS%id_wpup, wpup, CS%diag)
1288	24	if (CS%id_Rml > 0) call post_data(CS%id_Rml, buoy_av, CS%diag)
1289	24	if (CS%id_BLD > 0) call post_data(CS%id_BLD, little_h, CS%diag)
1290	24	if (CS%id_MLD > 0) call post_data(CS%id_MLD, big_H, CS%diag)
1291	24	if (CS%id_uhml > 0) call post_data(CS%id_uhml, uhml, CS%diag)
1292	24	if (CS%id_vhml > 0) call post_data(CS%id_vhml, vhml, CS%diag)
1293	24	if (CS%id_uDml > 0) call post_data(CS%id_uDml, uDml_diag, CS%diag)
1294	24	if (CS%id_vDml > 0) call post_data(CS%id_vDml, vDml_diag, CS%diag)
1295	24	if (CS%id_lfbod > 0) call post_data(CS%id_lfbod, lf_bodner_diag, CS%diag)
1296
1297	24	if (CS%id_uml > 0) then
1298	0	do j=js,je ; do I=is-1,ie
1299	0	h_vel = 0.5*((htot(i,j) + htot(i+1,j)) + h_neglect)
1300	0	uDml_diag(I,j) = uDml_diag(I,j) / (0.01h_vel) G%IdyCu(I,j) * (mu(0.,0.)-mu(-.01,0.))
1301		enddo ; enddo
1302	0	call post_data(CS%id_uml, uDml_diag, CS%diag)
1303		endif
1304	24	if (CS%id_vml > 0) then
1305	0	do J=js-1,je ; do i=is,ie
1306	0	h_vel = 0.5*((htot(i,j) + htot(i,j+1)) + h_neglect)
1307	0	vDml_diag(i,J) = vDml_diag(i,J) / (0.01h_vel) G%IdxCv(i,J) * (mu(0.,0.)-mu(-.01,0.))
1308		enddo ; enddo
1309	0	call post_data(CS%id_vml, vDml_diag, CS%diag)
1310		endif
1311		endif
1312
1313	24	end subroutine mixedlayer_restrat_Bodner
1314
1315		!> Two time-scale running mean in the same arbitrary units as "signal" and "filtered"
1316		!!
1317		!! If signal > filtered, returns running-mean with time scale "tau_growing".
1318		!! If signal <= filtered, returns running-mean with time scale "tau_decaying".
1319		!!
1320		!! The running mean of \f$ s \f$ with time scale "of \f$ \tau \f$ is:
1321		!! \f[
1322		!! \bar{s} <- ( \Delta t * s + \tau * \bar{s} ) / ( \Delta t + \tau )
1323		!! \f]
1324		!!
1325		!! Note that if \f$ tau=0 \f$, then the running mean equals the signal. Thus,
1326		!! rmean2ts with tau_growing=0 recovers the "resetting running mean" used in OM4.
1327	544608	real elemental function rmean2ts(signal, filtered, tau_growing, tau_decaying, dt)
1328		! Arguments
1329		real, intent(in) :: signal ! Unfiltered signal in arbitrary units [A]
1330		real, intent(in) :: filtered ! Current value of running mean in the same arbitrary units [A]
1331		real, intent(in) :: tau_growing ! Time scale for growing signal [T ~> s]
1332		real, intent(in) :: tau_decaying ! Time scale for decaying signal [T ~> s]
1333		real, intent(in) :: dt ! Time step [T ~> s]
1334		! Local variables
1335		real :: afac, bfac ! Non-dimensional fractional weights [nondim]
1336		real :: rt ! Reciprocal time scale [T-1 ~> s-1]
1337
1338	544608	if (signal>=filtered) then
1339	530126	rt = 1.0 / ( dt + tau_growing )
1340	530126	aFac = tau_growing * rt
1341	530126	bFac = 1. - aFac
1342		else
1343	14482	rt = 1.0 / ( dt + tau_decaying )
1344	14482	aFac = tau_decaying * rt
1345	14482	bFac = 1. - aFac
1346		endif
1347
1348	544608	rmean2ts = aFac * filtered + bFac * signal
1349
1350	544608	end function rmean2ts
1351
1352		!> Calculates a restratifying flow assuming a 2-layer bulk mixed layer.
1353	0	subroutine mixedlayer_restrat_BML(h, uhtr, vhtr, tv, forces, dt, G, GV, US, CS)
1354		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1355		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1356		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1357		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
1358		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: uhtr !< Accumulated zonal mass flux
1359		!! [H L2 ~> m3 or kg]
1360		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: vhtr !< Accumulated meridional mass flux
1361		!! [H L2 ~> m3 or kg]
1362		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables structure
1363		type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
1364		real, intent(in) :: dt !< Time increment [T ~> s]
1365		type(mixedlayer_restrat_CS), intent(inout) :: CS !< Module control structure
1366
1367		! Local variables
1368	0	real :: uhml(SZIB_(G),SZJ_(G),SZK_(GV)) ! Restratifying zonal thickness transports [H L2 T-1 ~> m3 s-1 or kg s-1]
1369	0	real :: vhml(SZI_(G),SZJB_(G),SZK_(GV)) ! Restratifying meridional thickness transports [H L2 T-1 ~> m3 s-1 or kg s-1]
1370		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
1371	0	h_avail ! The volume available for diffusion out of each face of each
1372		! sublayer of the mixed layer, divided by dt [H L2 T-1 ~> m3 s-1 or kg s-1].
1373		real, dimension(SZI_(G),SZJ_(G)) :: &
1374	0	U_star_2d, & ! The wind friction velocity in thickness-based units, calculated using
1375		! the Boussinesq reference density or the time-evolving surface density
1376		! in non-Boussinesq mode [H T-1 ~> m s-1 or kg m-2 s-1]
1377	0	htot, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2]
1378	0	Rml_av ! g_Rho0 times the average mixed layer density or negative G_Earth
1379		! times the average specific volume [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2]
1380		real :: g_Rho0 ! G_Earth/Rho0 times a thickness conversion factor
1381		! [L2 H-1 T-2 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2]
1382	0	real :: Rho_ml(SZI_(G)) ! Potential density relative to the surface [R ~> kg m-3]
1383	0	real :: rho_int(SZI_(G)) ! The integral of density over the mixed layer depth [R H ~> kg m-2 or kg2 m-5]
1384	0	real :: SpV_ml(SZI_(G)) ! Specific volume evaluated at the surface pressure [R-1 ~> m3 kg-1]
1385	0	real :: SpV_int(SZI_(G)) ! Specific volume integrated through the surface layer [H R-1 ~> m4 kg-1 or m]
1386	0	real :: p0(SZI_(G)) ! A pressure of 0 [R L2 T-2 ~> Pa]
1387
1388		real :: h_vel ! htot interpolated onto velocity points [H ~> m or kg m-2]
1389		real :: absf ! absolute value of f, interpolated to velocity points [T-1 ~> s-1]
1390		real :: u_star ! surface friction velocity, interpolated to velocity points and recast into
1391		! thickness-based units [H T-1 ~> m s-1 or kg m-2 s-1].
1392		real :: vonKar_x_pi2 ! A scaling constant that is approximately the von Karman constant times
1393		! pi squared [nondim]
1394		real :: mom_mixrate ! rate at which momentum is homogenized within mixed layer [T-1 ~> s-1]
1395		real :: timescale ! mixing growth timescale [T ~> s]
1396		real :: h_min ! The minimum layer thickness [H ~> m or kg m-2]. h_min could be 0.
1397		real :: h_neglect ! tiny thickness usually lost in roundoff and can be neglected [H ~> m or kg m-2]
1398		real :: I4dt ! 1/(4 dt) [T-1 ~> s-1]
1399		real :: I2htot ! Twice the total mixed layer thickness at velocity points [H ~> m or kg m-2]
1400		real :: z_topx2 ! depth of the top of a layer at velocity points [H ~> m or kg m-2]
1401		real :: hx2 ! layer thickness at velocity points [H ~> m or kg m-2]
1402	0	real :: a(SZK_(GV)) ! A non-dimensional value relating the overall flux magnitudes (uDml & vDml)
1403		! to the realized flux in a layer [nondim]. The vertical sum of a()
1404		! through the pieces of the mixed layer must be 0.
1405	0	real :: uDml(SZIB_(G)) ! Zonal volume fluxes in the upper half of the mixed layer [H L2 T-1 ~> m3 s-1 or kg s-1]
1406	0	real :: vDml(SZI_(G)) ! Meridional volume fluxes in the upper half of the mixed layer [H L2 T-1 ~> m3 s-1 or kg s-1]
1407	0	real :: utimescale_diag(SZIB_(G),SZJ_(G)) ! Zonal restratification timescale [T ~> s], stored for diagnostics.
1408	0	real :: vtimescale_diag(SZI_(G),SZJB_(G)) ! Meridional restratification timescale [T ~> s], stored for diagnostics.
1409	0	real :: uDml_diag(SZIB_(G),SZJ_(G)) ! A 2D copy of uDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1]
1410	0	real :: vDml_diag(SZI_(G),SZJB_(G)) ! A 2D copy of vDml for diagnostics [H L2 T-1 ~> m3 s-1 or kg s-1]
1411		logical :: use_EOS ! If true, density is calculated from T & S using an equation of state.
1412
1413		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
1414		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkml
1415	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1416	0	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nkml = GV%nkml
1417
1418	0	if (.not. CS%initialized) call MOM_error(FATAL, "mixedlayer_restrat_BML: "// &
1419	0	"Module must be initialized before it is used.")
1420
1421	0	if ((nkml<2) .or. (CS%ml_restrat_coef<=0.0)) return
1422
1423
1424	0	h_min = 0.5*GV%Angstrom_H ! This should be GV%Angstrom_H, but that value would change answers.
1425	0	uDml(:) = 0.0 ; vDml(:) = 0.0
1426	0	I4dt = 0.25 / dt
1427	0	g_Rho0 = GV%H_to_Z * GV%g_Earth / GV%Rho0
1428	0	vonKar_x_pi2 = CS%vonKar * 9.8696
1429	0	use_EOS = associated(tv%eqn_of_state)
1430	0	h_neglect = GV%H_subroundoff
1431
1432	0	if (.not.use_EOS) call MOM_error(FATAL, "mixedlayer_restrat_BML: "// &
1433	0	"An equation of state must be used with this module.")
1434
1435	0	if (CS%use_Stanley_ML) call MOM_error(FATAL, "mixedlayer_restrat_BML: "// &
1436	0	"The Stanley parameterization is not available with the BML.")
1437
1438		! Extract the friction velocity from the forcing type.
1439	0	call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1, H_T_units=.true.)
1440
1441		! Fix this later for nkml >= 3.
1442
1443	0	p0(:) = 0.0
1444	0	EOSdom(:) = EOS_domain(G%HI, halo=1)
1445		!$OMP parallel default(shared) private(Rho_ml,rho_int,h_vel,u_star,absf,mom_mixrate,timescale, &
1446		!$OMP SpV_ml,SpV_int,I2htot,z_topx2,hx2,a) &
1447		!$OMP firstprivate(uDml,vDml)
1448
1449	0	if (GV%Boussinesq .or. GV%semi_Boussinesq) then
1450		!$OMP do
1451	0	do j=js-1,je+1
1452	0	do i=is-1,ie+1
1453	0	htot(i,j) = 0.0 ; rho_int(i) = 0.0
1454		enddo
1455	0	do k=1,nkml
1456	0	call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, Rho_ml(:), tv%eqn_of_state, EOSdom)
1457	0	do i=is-1,ie+1
1458	0	rho_int(i) = rho_int(i) + h(i,j,k)*Rho_ml(i)
1459	0	htot(i,j) = htot(i,j) + h(i,j,k)
1460	0	h_avail(i,j,k) = max(I4dtG%areaT(i,j)(h(i,j,k)-GV%Angstrom_H),0.0)
1461		enddo
1462		enddo
1463
1464	0	do i=is-1,ie+1
1465	0	Rml_av(i,j) = (g_Rho0*rho_int(i)) / (htot(i,j) + h_neglect)
1466		enddo
1467		enddo
1468		else ! This is only used in non-Boussinesq mode.
1469		!$OMP do
1470	0	do j=js-1,je+1
1471	0	do i=is-1,ie+1
1472	0	htot(i,j) = 0.0 ; SpV_int(i) = 0.0
1473		enddo
1474	0	do k=1,nkml
1475	0	call calculate_spec_vol(tv%T(:,j,k), tv%S(:,j,k), p0, SpV_ml, tv%eqn_of_state, EOSdom)
1476	0	do i=is-1,ie+1
1477	0	SpV_int(i) = SpV_int(i) + h(i,j,k)*SpV_ml(i) ! [H R-1 ~> m4 kg-1 or m]
1478	0	htot(i,j) = htot(i,j) + h(i,j,k)
1479	0	h_avail(i,j,k) = max(I4dtG%areaT(i,j)(h(i,j,k)-GV%Angstrom_H),0.0)
1480		enddo
1481		enddo
1482
1483		! Convert the vertically integrated specific volume into a negative variable with units of density.
1484	0	do i=is-1,ie+1
1485	0	Rml_av(i,j) = (-GV%H_to_RZGV%g_Earth SpV_int(i)) / (htot(i,j) + h_neglect)
1486		enddo
1487		enddo
1488		endif
1489
1490		! TO DO:
1491		! 1. Mixing extends below the mixing layer to the mixed layer. Find it!
1492		! 2. Add exponential tail to stream-function?
1493
1494		! U - Component
1495		!$OMP do
1496	0	do j=js,je ; do I=is-1,ie
1497	0	h_vel = 0.5*(htot(i,j) + htot(i+1,j))
1498
1499	0	u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i+1,j)))
1500
1501	0	absf = 0.5*(abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J)))
1502
1503		! NOTE: growth_time changes answers on some systems, see below.
1504		! timescale = growth_time(u_star, h_vel, absf, h_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef)
1505
1506		! peak ML visc: u_star * von_Karman * (h_mlu_star)/(absfh_ml + 4.0*u_star)
1507		! momentum mixing rate: pi^2*visc/h_ml^2
1508		mom_mixrate = vonKar_x_pi2u_star*2 / &
1509	0	(absfh_vel2 + 4.0(h_vel+h_neglect)*u_star)
1510	0	timescale = 0.0625 * (absf + 2.0mom_mixrate) / (absf2 + mom_mixrate*2)
1511
1512	0	timescale = timescale * CS%ml_restrat_coef
1513		! timescale = timescale(2?)(L_def/L_MLI) * min(EKE/MKE,1.0 + (G%dyCv(i,j)/L_def)**2)
1514
1515		uDml(I) = timescale * G%dyCu(I,j)G%IdxCu_OBCmask(I,j) &
1516	0	(Rml_av(i+1,j)-Rml_av(i,j)) * (h_vel**2)
1517
1518	0	if (uDml(I) == 0) then
1519	0	do k=1,nkml ; uhml(I,j,k) = 0.0 ; enddo
1520		else
1521	0	I2htot = 1.0 / (htot(i,j) + htot(i+1,j) + h_neglect)
1522	0	z_topx2 = 0.0
1523		! a(k) relates the sublayer transport to uDml with a linear profile.
1524		! The sum of a(k) through the mixed layers must be 0.
1525	0	do k=1,nkml
1526	0	hx2 = (h(i,j,k) + h(i+1,j,k) + h_neglect)
1527	0	a(k) = (hx2 * I2htot) * (2.0 - 4.0(z_topx2+0.5hx2)*I2htot)
1528	0	z_topx2 = z_topx2 + hx2
1529	0	if (a(k)*uDml(I) > 0.0) then
1530	0	if (a(k)*uDml(I) > h_avail(i,j,k)) uDml(I) = h_avail(i,j,k) / a(k)
1531		else
1532	0	if (-a(k)*uDml(I) > h_avail(i+1,j,k)) uDml(I) = -h_avail(i+1,j,k)/a(k)
1533		endif
1534		enddo
1535	0	do k=1,nkml
1536	0	uhml(I,j,k) = a(k)*uDml(I)
1537	0	uhtr(I,j,k) = uhtr(I,j,k) + uhml(I,j,k)*dt
1538		enddo
1539		endif
1540
1541	0	uDml_diag(I,j) = uDml(I)
1542	0	utimescale_diag(I,j) = timescale
1543		enddo ; enddo
1544
1545		! V- component
1546		!$OMP do
1547	0	do J=js-1,je ; do i=is,ie
1548	0	h_vel = 0.5*(htot(i,j) + htot(i,j+1))
1549
1550	0	u_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i,j+1)))
1551
1552	0	absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))
1553
1554		! NOTE: growth_time changes answers on some systems, see below.
1555		! timescale = growth_time(u_star, h_vel, absf, h_neglect, CS%vonKar, CS%Kv_restrat, CS%ml_restrat_coef)
1556
1557		! peak ML visc: u_star * von_Karman * (h_mlu_star)/(absfh_ml + 4.0*u_star)
1558		! momentum mixing rate: pi^2*visc/h_ml^2
1559		mom_mixrate = vonKar_x_pi2u_star*2 / &
1560	0	(absfh_vel2 + 4.0(h_vel+h_neglect)*u_star)
1561	0	timescale = 0.0625 * (absf + 2.0mom_mixrate) / (absf2 + mom_mixrate*2)
1562
1563	0	timescale = timescale * CS%ml_restrat_coef
1564		! timescale = timescale(2?)(L_def/L_MLI) * min(EKE/MKE,1.0 + (G%dyCv(i,j)/L_def)**2)
1565
1566		vDml(i) = timescale * G%dxCv(i,J)G%IdyCv_OBCmask(i,J) &
1567	0	(Rml_av(i,j+1)-Rml_av(i,j)) * (h_vel**2)
1568	0	if (vDml(i) == 0) then
1569	0	do k=1,nkml ; vhml(i,J,k) = 0.0 ; enddo
1570		else
1571	0	I2htot = 1.0 / (htot(i,j) + htot(i,j+1) + h_neglect)
1572	0	z_topx2 = 0.0
1573		! a(k) relates the sublayer transport to vDml with a linear profile.
1574		! The sum of a(k) through the mixed layers must be 0.
1575	0	do k=1,nkml
1576	0	hx2 = (h(i,j,k) + h(i,j+1,k) + h_neglect)
1577	0	a(k) = (hx2 * I2htot) * (2.0 - 4.0(z_topx2+0.5hx2)*I2htot)
1578	0	z_topx2 = z_topx2 + hx2
1579	0	if (a(k)*vDml(i) > 0.0) then
1580	0	if (a(k)*vDml(i) > h_avail(i,j,k)) vDml(i) = h_avail(i,j,k) / a(k)
1581		else
1582	0	if (-a(k)*vDml(i) > h_avail(i,j+1,k)) vDml(i) = -h_avail(i,j+1,k)/a(k)
1583		endif
1584		enddo
1585	0	do k=1,nkml
1586	0	vhml(i,J,k) = a(k)*vDml(i)
1587	0	vhtr(i,J,k) = vhtr(i,J,k) + vhml(i,J,k)*dt
1588		enddo
1589		endif
1590
1591	0	vtimescale_diag(i,J) = timescale
1592	0	vDml_diag(i,J) = vDml(i)
1593		enddo ; enddo
1594
1595		!$OMP do
1596	0	do j=js,je ; do k=1,nkml ; do i=is,ie
1597		h(i,j,k) = h(i,j,k) - dtG%IareaT(i,j) &
1598	0	((uhml(I,j,k) - uhml(I-1,j,k)) + (vhml(i,J,k) - vhml(i,J-1,k)))
1599	0	if (h(i,j,k) < h_min) h(i,j,k) = h_min
1600		enddo ; enddo ; enddo
1601		!$OMP end parallel
1602
1603		! Whenever thickness changes let the diag manager know, target grids
1604		! for vertical remapping may need to be regenerated.
1605	0	if (CS%id_uhml > 0 .or. CS%id_vhml > 0) &
1606		! Remapped uhml and vhml require east/north halo updates of h
1607	0	call pass_var(h, G%domain, To_West+To_South+Omit_Corners, halo=1)
1608	0	call diag_update_remap_grids(CS%diag)
1609
1610		! Offer diagnostic fields for averaging.
1611	0	if (query_averaging_enabled(CS%diag) .and. &
1612		((CS%id_urestrat_time > 0) .or. (CS%id_vrestrat_time > 0))) then
1613	0	call post_data(CS%id_urestrat_time, utimescale_diag, CS%diag)
1614	0	call post_data(CS%id_vrestrat_time, vtimescale_diag, CS%diag)
1615		endif
1616	0	if (query_averaging_enabled(CS%diag) .and. &
1617		((CS%id_uhml>0) .or. (CS%id_vhml>0))) then
1618	0	do k=nkml+1,nz
1619	0	do j=js,je ; do I=Isq,Ieq ; uhml(I,j,k) = 0.0 ; enddo ; enddo
1620	0	do J=Jsq,Jeq ; do i=is,ie ; vhml(i,J,k) = 0.0 ; enddo ; enddo
1621		enddo
1622	0	if (CS%id_uhml > 0) call post_data(CS%id_uhml, uhml, CS%diag)
1623	0	if (CS%id_vhml > 0) call post_data(CS%id_vhml, vhml, CS%diag)
1624	0	if (CS%id_MLD > 0) call post_data(CS%id_MLD, htot, CS%diag)
1625	0	if (CS%id_Rml > 0) call post_data(CS%id_Rml, Rml_av, CS%diag)
1626	0	if (CS%id_uDml > 0) call post_data(CS%id_uDml, uDml_diag, CS%diag)
1627	0	if (CS%id_vDml > 0) call post_data(CS%id_vDml, vDml_diag, CS%diag)
1628		endif
1629
1630		end subroutine mixedlayer_restrat_BML
1631
1632		!> Detects the mixed layer depth using a density difference criterion (MLE_DENSITY_DIFF)
1633	0	subroutine detect_mld(h, tv, MLD_fast, G, GV, CS)
1634		type(mixedlayer_restrat_CS), intent(inout) :: CS !< Module control structure
1635		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
1636		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1637		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
1638		real, dimension(SZI_(G),SZJ_(G)), intent(out) :: MLD_fast !< detected mixed layer depth [H ~> m or kg m-2]
1639		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamic variables structure
1640
1641		! Local variables
1642	0	real, dimension(SZI_(G)) :: pRef_MLD ! A reference pressure for calculating the mixed layer
1643		! densities [R L2 T-2 ~> Pa].
1644	0	real, dimension(SZI_(G)) :: rhoSurf, deltaRhoAtKm1, deltaRhoAtK ! Densities and density differences [R ~> kg m-3]
1645	0	real, dimension(SZI_(G)) :: dK, dKm1 ! Depths of layer centers [H ~> m or kg m-2].
1646		real :: ddRho ! A density difference [R ~> kg m-3]
1647		real :: aFac ! A nondimensional ratio [nondim]
1648	0	real :: covTS(SZI_(G)) ! SGS TS covariance in Stanley param; currently 0 [C S ~> degC ppt]
1649	0	real :: varS(SZI_(G)) ! SGS S variance in Stanley param; currently 0 [S2 ~> ppt2]
1650		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
1651		integer :: i, j, k, is, ie, js, je, nz
1652
1653	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1654
1655	0	covTS(:) = 0.0 ! Might be in tv% in the future. Not implemented for the time being.
1656	0	varS(:) = 0.0 ! Ditto.
1657
1658		!! TODO: use derivatives and mid-MLD pressure. Currently this is sigma-0. -AJA
1659	0	pRef_MLD(:) = 0.
1660	0	EOSdom(:) = EOS_domain(G%HI, halo=1)
1661	0	do j=js-1,je+1
1662	0	dK(:) = 0.5 * h(:,j,1) ! Depth of center of surface layer
1663	0	if (CS%use_Stanley_ML) then
1664		call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, tv%varT(:,j,1), covTS, varS, &
1665	0	rhoSurf, tv%eqn_of_state, EOSdom)
1666		else
1667	0	call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pRef_MLD, rhoSurf, tv%eqn_of_state, EOSdom)
1668		endif
1669	0	deltaRhoAtK(:) = 0.
1670	0	MLD_fast(:,j) = 0.
1671	0	do k=2,nz
1672	0	dKm1(:) = dK(:) ! Depth of center of layer K-1
1673	0	dK(:) = dK(:) + 0.5 * ( h(:,j,k) + h(:,j,k-1) ) ! Depth of center of layer K
1674		! Mixed-layer depth, using sigma-0 (surface reference pressure)
1675	0	deltaRhoAtKm1(:) = deltaRhoAtK(:) ! Store value from previous iteration of K
1676	0	if (CS%use_Stanley_ML) then
1677		call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, tv%varT(:,j,k), covTS, varS, &
1678	0	deltaRhoAtK, tv%eqn_of_state, EOSdom)
1679		else
1680	0	call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pRef_MLD, deltaRhoAtK, tv%eqn_of_state, EOSdom)
1681		endif
1682	0	do i=is-1,ie+1
1683	0	deltaRhoAtK(i) = deltaRhoAtK(i) - rhoSurf(i) ! Density difference between layer K and surface
1684		enddo
1685	0	do i=is-1,ie+1
1686	0	ddRho = deltaRhoAtK(i) - deltaRhoAtKm1(i)
1687		if ((MLD_fast(i,j)==0.) .and. (ddRho>0.) .and. &
1688	0	(deltaRhoAtKm1(i)<CS%MLE_density_diff) .and. (deltaRhoAtK(i)>=CS%MLE_density_diff)) then
1689	0	aFac = ( CS%MLE_density_diff - deltaRhoAtKm1(i) ) / ddRho
1690	0	MLD_fast(i,j) = dK(i) * aFac + dKm1(i) * (1. - aFac)
1691		endif
1692		enddo ! i-loop
1693		enddo ! k-loop
1694	0	do i=is-1,ie+1
1695	0	MLD_fast(i,j) = CS%MLE_MLD_stretch * MLD_fast(i,j)
1696	0	if ((MLD_fast(i,j)==0.) .and. (deltaRhoAtK(i)<CS%MLE_density_diff)) &
1697	0	MLD_fast(i,j) = dK(i) ! Assume mixing to the bottom
1698		enddo
1699		enddo ! j-loop
1700	0	end subroutine detect_mld
1701
1702		! NOTE: This function appears to change answers on some platforms, so it is
1703		! currently unused in the model, but we intend to introduce it in the future.
1704
1705		!> Return the growth timescale for the submesoscale mixed layer eddies in [T ~> s]
1706	0	real function growth_time(u_star, hBL, absf, h_neg, vonKar, Kv_rest, restrat_coef)
1707		real, intent(in) :: u_star !< Surface friction velocity in thickness-based units [H T-1 ~> m s-1 or kg m-2 s-1]
1708		real, intent(in) :: hBL !< Boundary layer thickness including at least a negligible
1709		!! value to keep it positive definite [H ~> m or kg m-2]
1710		real, intent(in) :: absf !< Absolute value of the Coriolis parameter [T-1 ~> s-1]
1711		real, intent(in) :: h_neg !< A tiny thickness that is usually lost in roundoff so can be
1712		!! neglected [H ~> m or kg m-2]
1713		real, intent(in) :: Kv_rest !< The background laminar vertical viscosity used for restratification,
1714		!! rescaled into thickness-based units [H2 T-1 ~> m2 s-1 or kg2 m-4 s-1]
1715		real, intent(in) :: vonKar !< The von Karman constant, used to scale the turbulent limits
1716		!! on the restratification timescales [nondim]
1717		real, intent(in) :: restrat_coef !< An overall scaling factor for the restratification timescale [nondim]
1718
1719		! Local variables
1720		real :: mom_mixrate ! rate at which momentum is homogenized within mixed layer [T-1 ~> s-1]
1721		real :: Kv_eff ! An effective overall viscosity in thickness-based units [H2 T-1 ~> m2 s-1 or kg2 m-4 s-1]
1722		real :: pi2 ! A scaling constant that is approximately pi^2 [nondim]
1723
1724		! peak ML visc: u_star * von_Karman * (h_mlu_star)/(absfh_ml + 4.0*u_star) + Kv_water
1725		! momentum mixing rate: pi^2*visc/h_ml^2
1726	0	pi2 = 9.8696 ! Approximately pi^2. This is more accurate than the overall uncertainty of the
1727		! scheme, with a value that is chosen to reproduce previous answers.
1728	0	if (Kv_rest <= 0.0) then
1729		! This case reproduces the previous answers, but the extra h_neg is otherwise unnecessary.
1730	0	mom_mixrate = (pi2vonKar)u_star*2 / (absfhBL*2 + 4.0(hBL + h_neg)*u_star)
1731	0	growth_time = restrat_coef * (0.0625 * (absf + 2.0mom_mixrate) / (absf2 + mom_mixrate*2))
1732		else
1733		! Set the mixing rate to the sum of a turbulent mixing rate and a laminar viscous rate.
1734		! mom_mixrate = pi2vonKaru_star*2 / (absfhBL*2 + 4.0hBLu_star) + pi2Kv_rest / hBL**2
1735	0	if (absfhBL <= 4.0e-16u_star) then
1736	0	Kv_eff = pi2 * (Kv_rest + 0.25vonKarhBL*u_star)
1737		else
1738	0	Kv_eff = pi2 * (Kv_rest + vonKaru_star2hBL / (absfhBL + 4.0u_star))
1739		endif
1740	0	growth_time = (restrat_coef0.0625) ((hBL*2(hBL*2absf + 2.0Kv_eff)) / ((hBL2absf)2 + Kv_eff2))
1741		endif
1742
1743	0	end function growth_time
1744
1745		!> Initialize the mixed layer restratification module
1746	1	logical function mixedlayer_restrat_init(Time, G, GV, US, param_file, diag, CS, restart_CS)
1747		type(time_type), target, intent(in) :: Time !< Current model time
1748		type(ocean_grid_type), intent(inout) :: G !< Ocean grid structure
1749		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1750		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1751		type(param_file_type), intent(in) :: param_file !< Parameter file to parse
1752		type(diag_ctrl), target, intent(inout) :: diag !< Regulate diagnostics
1753		type(mixedlayer_restrat_CS), intent(inout) :: CS !< Module control structure
1754		type(MOM_restart_CS), intent(in) :: restart_CS !< MOM restart control structure
1755
1756		! Local variables
1757		real :: flux_to_kg_per_s ! A unit conversion factor for fluxes. [kg T s-1 H-1 L-2 ~> kg m-3 or 1]
1758		real :: omega ! The Earth's rotation rate [T-1 ~> s-1].
1759		real :: ustar_min_dflt ! The default value for RESTRAT_USTAR_MIN [Z T-1 ~> m s-1]
1760		real :: Stanley_coeff ! Coefficient relating the temperature gradient and sub-gridscale
1761		! temperature variance [nondim]
1762		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags
1763		! This include declares and sets the variable "version".
1764		character(len=200) :: inputdir ! The directory where NetCDF input files
1765		character(len=240) :: mle_fl_filename ! A file from which chl_a concentrations are to be read.
1766		character(len=128) :: mle_fl_file ! Data containing MLE front-length scale. Used
1767		! when reading from file.
1768		character(len=32) :: fl_varname ! Name of front-length scale variable in mle_fl_file.
1769
1770		# include "version_variable.h"
1771		character(len=200) :: filename, varname
1772
1773		! Read all relevant parameters and write them to the model log.
1774		call get_param(param_file, mdl, "MIXEDLAYER_RESTRAT", mixedlayer_restrat_init, &
1775	1	default=.false., do_not_log=.true.)
1776	1	call log_version(param_file, mdl, version, "", all_default=.not.mixedlayer_restrat_init)
1777		call get_param(param_file, mdl, "MIXEDLAYER_RESTRAT", mixedlayer_restrat_init, &
1778		"If true, a density-gradient dependent re-stratifying "//&
1779		"flow is imposed in the mixed layer. Can be used in ALE mode "//&
1780		"without restriction but in layer mode can only be used if "//&
1781	1	"BULKMIXEDLAYER is true.", default=.false.)
1782	1	if (.not. mixedlayer_restrat_init) return
1783
1784	1	CS%initialized = .true.
1785	1	CS%Time => Time
1786
1787		! Nonsense values to cause problems when these parameters are not used
1788	1	CS%MLE_MLD_decay_time = -9.e9*US%s_to_T
1789	1	CS%MLE_density_diff = -9.e9*US%kg_m3_to_R
1790	1	CS%MLE_tail_dh = -9.e9
1791	1	CS%MLE_use_PBL_MLD = .false.
1792	1	CS%MLE_MLD_stretch = -9.e9
1793	1	CS%use_Stanley_ML = .false.
1794	1	CS%use_Bodner = .false.
1795	1	CS%fl_from_file = .false.
1796	1	CS%MLD_grid = .false.
1797	1	CS%Cr_grid = .false.
1798
1799	1	call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false., do_not_log=.true.)
1800		call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
1801		"This sets the default value for the various _ANSWER_DATE parameters.", &
1802	1	default=99991231, do_not_log=.true.)
1803	1	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1804	1	call openParameterBlock(param_file,'MLE') ! Prepend MLE% to all parameters
1805	1	if (GV%nkml==0) then
1806		call get_param(param_file, mdl, "USE_BODNER23", CS%use_Bodner, &
1807		"If true, use the Bodner et al., 2023, formulation of the re-stratifying "//&
1808		"mixed-layer restratification parameterization. This only works in ALE mode.", &
1809	1	default=.false.)
1810		endif
1811	1	if (CS%use_Bodner) then
1812		call get_param(param_file, mdl, "CR", CS%Cr, &
1813		"The efficiency coefficient in eq 27 of Bodner et al., 2023.", &
1814	1	units="nondim", default=0.0)
1815		call get_param(param_file, mdl, "BODNER_NSTAR", CS%Nstar, &
1816		"The n* value used to estimate the turbulent vertical momentum flux "//&
1817		"in Bodner et al., 2023, eq. 18. This is independent of the value used in "//&
1818		"the PBL scheme but should be set to be the same for consistency.", &
1819	1	units="nondim", default=0.066)
1820		call get_param(param_file, mdl, "BODNER_MSTAR", CS%Mstar, &
1821		"The m* value used to estimate the turbulent vertical momentum flux "//&
1822		"in Bodner et al., 2023, eq. 18. This is independent of the value used in "//&
1823		"the PBL scheme but should be set to be the same for consistency.", &
1824	1	units="nondim", default=0.5)
1825		call get_param(param_file, mdl, "BLD_GROWING_TFILTER", CS%BLD_growing_Tfilt, &
1826		"The time-scale for a running-mean filter applied to the boundary layer "//&
1827		"depth (BLD) when the BLD is deeper than the running mean. A value of 0 "//&
1828		"instantaneously sets the running mean to the current value of BLD.", &
1829	1	units="s", default=0., scale=US%s_to_T)
1830		call get_param(param_file, mdl, "BLD_DECAYING_TFILTER", CS%BLD_decaying_Tfilt, &
1831		"The time-scale for a running-mean filter applied to the boundary layer "//&
1832		"depth (BLD) when the BLD is shallower than the running mean. A value of 0 "//&
1833		"instantaneously sets the running mean to the current value of BLD.", &
1834	1	units="s", default=0., scale=US%s_to_T)
1835		call get_param(param_file, mdl, "MLD_GROWING_TFILTER", CS%MLD_growing_Tfilt, &
1836		"The time-scale for a running-mean filter applied to the time-filtered "//&
1837		"BLD, when the latter is deeper than the running mean. A value of 0 "//&
1838		"instantaneously sets the running mean to the current value filtered BLD.", &
1839	1	units="s", default=0., scale=US%s_to_T)
1840		call get_param(param_file, mdl, "MLD_DECAYING_TFILTER", CS%MLD_decaying_Tfilt, &
1841		"The time-scale for a running-mean filter applied to the time-filtered "//&
1842		"BLD, when the latter is shallower than the running mean. A value of 0 "//&
1843		"instantaneously sets the running mean to the current value filtered BLD.", &
1844	1	units="s", default=0., scale=US%s_to_T)
1845		call get_param(param_file, mdl, "ML_RESTRAT_ANSWER_DATE", CS%answer_date, &
1846		"The vintage of the order of arithmetic and expressions in the mixed layer "//&
1847		"restrat calculations. Values below 20240201 recover the answers from the end "//&
1848		"of 2023, while higher values use the new cuberoot function in the Bodner code "//&
1849		"to avoid needing to undo dimensional rescaling.", &
1850		default=default_answer_date, &
1851	1	do_not_log=.not.(CS%use_Bodner.and.(GV%Boussinesq.or.GV%semi_Boussinesq)))
1852		call get_param(param_file, mdl, "MIN_WSTAR2", CS%min_wstar2, &
1853		"The minimum lower bound to apply to the vertical momentum flux, w'u', "//&
1854		"in the Bodner et al., restratification parameterization. This avoids "//&
1855		"a division-by-zero in the limit when u* and the buoyancy flux are zero. "//&
1856		"The default is less than the molecular viscosity of water times the Coriolis "//&
1857		"parameter a micron away from the equator.", &
1858	1	units="m2 s-2", default=1.0e-24, scale=US%m_to_Z*2US%T_to_s**2)
1859		call get_param(param_file, mdl, "TAIL_DH", CS%MLE_tail_dh, &
1860		"Fraction by which to extend the mixed-layer restratification "//&
1861		"depth used for a smoother stream function at the base of "//&
1862	1	"the mixed-layer.", units="nondim", default=0.0)
1863		call get_param(param_file, mdl, "USE_STANLEY_TVAR", CS%use_Stanley_ML, &
1864		"If true, turn on Stanley SGS T variance parameterization "// &
1865	1	"in ML restrat code.", default=.false.)
1866		call get_param(param_file, mdl, "MLE_NKBLOCK", CS%nkblock, &
1867		"Vertical block size for the mixed layer density integral. "//&
1868		"0 processes all layers in a single block.", &
1869	1	default=default_nkblock)
1870		call get_param(param_file, mdl, "MLE_NJBLOCK", CS%njblock, &
1871		"j-tile size for the mixed layer density integral. "//&
1872		"0 uses the full j compute domain.", &
1873	1	default=default_njblock)
1874		call get_param(param_file, mdl, "USE_CR_GRID", CS%Cr_grid, &
1875	1	"If true, read in a spatially varying Cr field.", default=.false.)
1876		call get_param(param_file, mdl, "USE_MLD_GRID", CS%MLD_grid, &
1877	1	"If true, read in a spatially varying MLD_decaying_Tfilt field.", default=.false.)
1878	1	if (CS%MLD_grid) then
1879		call get_param(param_file, mdl, "MLD_TFILT_FILE", filename, &
1880		"The path to the file containing the MLD_decaying_Tfilt fields.", &
1881	0	default="")
1882		call get_param(param_file, mdl, "MLD_TFILT_VAR", varname, &
1883		"The variable name for MLD_decaying_Tfilt field.", &
1884	0	default="MLD_tfilt")
1885	0	filename = trim(inputdir) // "/" // trim(filename)
1886	0	allocate(CS%MLD_Tfilt_space(G%isd:G%ied,G%jsd:G%jed), source=0.0)
1887	0	call MOM_read_data(filename, varname, CS%MLD_Tfilt_space, G%domain, scale=US%s_to_T)
1888	0	call pass_var(CS%MLD_Tfilt_space, G%domain)
1889		endif
1890	8773	allocate(CS%Cr_space(G%isd:G%ied,G%jsd:G%jed), source=CS%Cr)
1891	1	if (CS%Cr_grid) then
1892		call get_param(param_file, mdl, "CR_FILE", filename, &
1893		"The path to the file containing the Cr fields.", &
1894	0	default="")
1895		call get_param(param_file, mdl, "CR_VAR", varname, &
1896		"The variable name for Cr field.", &
1897	0	default="Cr")
1898	0	filename = trim(inputdir) // "/" // trim(filename)
1899	0	call MOM_read_data(filename, varname, CS%Cr_space, G%domain)
1900	0	call pass_var(CS%Cr_space, G%domain)
1901		endif
1902	1	call closeParameterBlock(param_file) ! The remaining parameters do not have MLE% prepended
1903		call get_param(param_file, mdl, "MLE_USE_PBL_MLD", CS%MLE_use_PBL_MLD, &
1904		"If true, the MLE parameterization will use the mixed-layer "//&
1905		"depth provided by the active PBL parameterization. If false, "//&
1906		"MLE will estimate a MLD based on a density difference with the "//&
1907		"surface using the parameter MLE_DENSITY_DIFF, unless "//&
1908	1	"BODNER_DETECT_MLD is true.", default=.false.)
1909		call get_param(param_file, mdl, "BODNER_DETECT_MLD", CS%Bodner_detect_MLD, &
1910		"If true, the Bodner parameterization will use the mixed-layer depth "//&
1911	1	"detected via the density difference criterion MLE_DENSITY_DIFF.", default=.false.)
1912	1	if (.not.(CS%MLE_use_PBL_MLD.or.CS%Bodner_detect_MLD)) call MOM_error(FATAL, "mixedlayer_restrat_init: "// &
1913	0	"To use MLE%USE_BODNER23=True then MLE_USE_PBL_MLD or BODNER_DETECT_MLD must be true.")
1914	1	if (CS%MLE_use_PBL_MLD.and.CS%Bodner_detect_MLD) call MOM_error(FATAL, "mixedlayer_restrat_init: "// &
1915	0	"MLE_USE_PBL_MLD and BODNER_DETECT_MLD cannot both be true.")
1916		else
1917	0	call closeParameterBlock(param_file) ! The remaining parameters do not have MLE% prepended
1918		endif
1919
1920	1	if (.not.CS%use_Bodner) then
1921		! This coefficient is used in both layered and ALE versions of Fox-Kemper but not Bodner
1922		call get_param(param_file, mdl, "FOX_KEMPER_ML_RESTRAT_COEF", CS%ml_restrat_coef, &
1923		"A nondimensional coefficient that is proportional to "//&
1924		"the ratio of the deformation radius to the dominant "//&
1925		"lengthscale of the submesoscale mixed layer "//&
1926		"instabilities, times the minimum of the ratio of the "//&
1927		"mesoscale eddy kinetic energy to the large-scale "//&
1928		"geostrophic kinetic energy or 1 plus the square of the "//&
1929		"grid spacing over the deformation radius, as detailed "//&
1930	0	"by Fox-Kemper et al. (2011)", units="nondim", default=0.0)
1931		! These parameters are only used in the OM4-era version of Fox-Kemper
1932		call get_param(param_file, mdl, "USE_STANLEY_ML", CS%use_Stanley_ML, &
1933		"If true, turn on Stanley SGS T variance parameterization "// &
1934	0	"in ML restrat code.", default=.false.)
1935	0	if (CS%use_Stanley_ML) then
1936		call get_param(param_file, mdl, "STANLEY_COEFF", Stanley_coeff, &
1937		"Coefficient correlating the temperature gradient and SGS T variance.", &
1938	0	units="nondim", default=-1.0, do_not_log=.true.)
1939	0	if (Stanley_coeff < 0.0) call MOM_error(FATAL, &
1940	0	"STANLEY_COEFF must be set >= 0 if USE_STANLEY_ML is true.")
1941		endif
1942		call get_param(param_file, mdl, 'VON_KARMAN_CONST', CS%vonKar, &
1943		'The value the von Karman constant as used for mixed layer viscosity.', &
1944	0	units='nondim', default=0.41)
1945		! We use GV%nkml to distinguish between the old and new implementation of MLE.
1946		! The old implementation only works for the layer model with nkml>0.
1947	0	if (GV%nkml==0) then
1948		call get_param(param_file, mdl, "FOX_KEMPER_ML_RESTRAT_COEF2", CS%ml_restrat_coef2, &
1949		"As for FOX_KEMPER_ML_RESTRAT_COEF but used in a second application "//&
1950	0	"of the MLE restratification parameterization.", units="nondim", default=0.0)
1951		call get_param(param_file, mdl, "MLE_FRONT_LENGTH", CS%front_length, &
1952		"If non-zero, is the frontal-length scale used to calculate the "//&
1953		"upscaling of buoyancy gradients that is otherwise represented "//&
1954		"by the parameter FOX_KEMPER_ML_RESTRAT_COEF. If MLE_FRONT_LENGTH is "//&
1955		"non-zero, it is recommended to set FOX_KEMPER_ML_RESTRAT_COEF=1.0.",&
1956	0	units="m", default=0.0, scale=US%m_to_L)
1957		call get_param(param_file, mdl, "MLE_FRONT_LENGTH_FROM_FILE", CS%fl_from_file, &
1958	0	"If true, the MLE front-length scale is read from a file.", default=.false.)
1959	0	if (CS%fl_from_file) then
1960	0	call time_interp_external_init()
1961
1962	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1963		call get_param(param_file, mdl, "MLE_FL_FILE", mle_fl_file, &
1964		"MLE_FL_FILE is the file containing MLE front-length scale. "//&
1965	0	"It is used when MLE_FRONT_LENGTH_FROM_FILE is true.", fail_if_missing=.true.)
1966	0	mle_fl_filename = trim(slasher(inputdir))//trim(mle_fl_file)
1967	0	call log_param(param_file, mdl, "INPUTDIR/MLE_FL_FILE", mle_fl_filename)
1968		call get_param(param_file, mdl, "FL_VARNAME", fl_varname, &
1969	0	"Name of MLE front-length scale variable in MLE_FL_FILE.", default='mle_fl')
1970	0	if (modulo(G%Domain%turns, 4) /= 0) then
1971	0	CS%sbc_fl = init_external_field(mle_fl_filename, trim(fl_varname), MOM_domain=G%Domain%domain_in)
1972		else
1973	0	CS%sbc_fl = init_external_field(mle_fl_filename, trim(fl_varname), MOM_domain=G%Domain)
1974		endif
1975		endif
1976	0	if (CS%fl_from_file .and. CS%front_length>0.0) call MOM_error(FATAL, "mixedlayer_restrat_init: "// &
1977		"MLE_FRONT_LENGTH_FROM_FILE cannot be true when MLE_FRONT_LENGTH > 0.0. "// &
1978		"If you want to use MLE_FRONT_LENGTH, set MLE_FRONT_LENGTH_FROM_FILE to false. " // &
1979	0	"If you want to use MLE_FRONT_LENGTH_FROM_FILE, set MLE_FRONT_LENGTH to 0.0.")
1980		call get_param(param_file, mdl, "MLE_USE_PBL_MLD", CS%MLE_use_PBL_MLD, &
1981		"If true, the MLE parameterization will use the mixed-layer "//&
1982		"depth provided by the active PBL parameterization. If false, "//&
1983		"MLE will estimate a MLD based on a density difference with the "//&
1984	0	"surface using the parameter MLE_DENSITY_DIFF.", default=.false.)
1985		call get_param(param_file, mdl, "MLE_MLD_DECAY_TIME", CS%MLE_MLD_decay_time, &
1986		"The time-scale for a running-mean filter applied to the mixed-layer "//&
1987		"depth used in the MLE restratification parameterization. When "//&
1988		"the MLD deepens below the current running-mean the running-mean "//&
1989	0	"is instantaneously set to the current MLD.", units="s", default=0., scale=US%s_to_T)
1990		call get_param(param_file, mdl, "MLE_MLD_DECAY_TIME2", CS%MLE_MLD_decay_time2, &
1991		"The time-scale for a running-mean filter applied to the filtered "//&
1992		"mixed-layer depth used in a second MLE restratification parameterization. "//&
1993		"When the MLD deepens below the current running-mean the running-mean "//&
1994	0	"is instantaneously set to the current MLD.", units="s", default=0., scale=US%s_to_T)
1995	0	if (.not. CS%MLE_use_PBL_MLD) then
1996		call get_param(param_file, mdl, "MLE_DENSITY_DIFF", CS%MLE_density_diff, &
1997		"Density difference used to detect the mixed-layer "//&
1998		"depth used for the mixed-layer eddy parameterization "//&
1999	0	"by Fox-Kemper et al. (2011)", units="kg/m3", default=0.03, scale=US%kg_m3_to_R)
2000		endif
2001		call get_param(param_file, mdl, "MLE_TAIL_DH", CS%MLE_tail_dh, &
2002		"Fraction by which to extend the mixed-layer restratification "//&
2003		"depth used for a smoother stream function at the base of "//&
2004	0	"the mixed-layer.", units="nondim", default=0.0)
2005		call get_param(param_file, mdl, "MLE_MLD_STRETCH", CS%MLE_MLD_stretch, &
2006		"A scaling coefficient for stretching/shrinking the MLD "//&
2007		"used in the MLE scheme. This simply multiplies MLD wherever used.",&
2008	0	units="nondim", default=1.0)
2009		endif
2010		call get_param(param_file, mdl, "KV_RESTRAT", CS%Kv_restrat, &
2011		"A small viscosity that sets a floor on the momentum mixing rate during "//&
2012		"restratification. If this is positive, it will prevent some possible "//&
2013		"divisions by zero even if ustar, RESTRAT_USTAR_MIN, and f are all 0.", &
2014	0	units="m2 s-1", default=0.0, scale=GV%m2_s_to_HZ_T(US%Z_to_mGV%m_to_H))
2015		call get_param(param_file, mdl, "OMEGA", omega, &
2016		"The rotation rate of the earth.", &
2017	0	units="s-1", default=7.2921e-5, scale=US%T_to_s)
2018	0	ustar_min_dflt = 2.0e-4 * omega * (GV%Angstrom_Z + GV%dZ_subroundoff)
2019		call get_param(param_file, mdl, "RESTRAT_USTAR_MIN", CS%ustar_min, &
2020		"The minimum value of ustar that will be used by the mixed layer "//&
2021		"restratification module. This can be tiny, but if this is greater than 0, "//&
2022		"it will prevent divisions by zero when f and KV_RESTRAT are zero.", &
2023	0	units="m s-1", default=US%Z_to_mUS%s_to_Tustar_min_dflt, scale=GV%m_to_H*US%T_to_s)
2024	1	elseif (CS%Bodner_detect_MLD) then
2025		call get_param(param_file, mdl, "MLE_DENSITY_DIFF", CS%MLE_density_diff, &
2026		"Density difference used to detect the mixed-layer "//&
2027		"depth used for the mixed-layer eddy parameterization "//&
2028	0	"by Fox-Kemper et al. (2010)", units="kg/m3", default=0.03, scale=US%kg_m3_to_R)
2029		call get_param(param_file, mdl, "MLE_MLD_STRETCH", CS%MLE_MLD_stretch, &
2030		"A scaling coefficient for stretching/shrinking the MLD "//&
2031		"used in the MLE scheme. This simply multiplies MLD wherever used.",&
2032	0	units="nondim", default=1.0)
2033		endif
2034
2035	1	CS%diag => diag
2036
2037	1	flux_to_kg_per_s = GV%H_to_kg_m2 * US%L_to_m*2 US%s_to_T
2038
2039		CS%id_uhml = register_diag_field('ocean_model', 'uhml', diag%axesCuL, Time, &
2040		'Zonal Thickness Flux to Restratify Mixed Layer', &
2041	1	'kg s-1', conversion=flux_to_kg_per_s, y_cell_method='sum', v_extensive=.true.)
2042		CS%id_vhml = register_diag_field('ocean_model', 'vhml', diag%axesCvL, Time, &
2043		'Meridional Thickness Flux to Restratify Mixed Layer', &
2044	1	'kg s-1', conversion=flux_to_kg_per_s, x_cell_method='sum', v_extensive=.true.)
2045		CS%id_urestrat_time = register_diag_field('ocean_model', 'MLu_restrat_time', diag%axesCu1, Time, &
2046	1	'Mixed Layer Zonal Restratification Timescale', 's', conversion=US%T_to_s)
2047		CS%id_vrestrat_time = register_diag_field('ocean_model', 'MLv_restrat_time', diag%axesCv1, Time, &
2048	1	'Mixed Layer Meridional Restratification Timescale', 's', conversion=US%T_to_s)
2049		CS%id_MLD = register_diag_field('ocean_model', 'MLD_restrat', diag%axesT1, Time, &
2050		'Mixed Layer Depth as used in the mixed-layer restratification parameterization', &
2051	1	'm', conversion=GV%H_to_m)
2052		CS%id_BLD = register_diag_field('ocean_model', 'BLD_restrat', diag%axesT1, Time, &
2053		'Boundary Layer Depth as used in the mixed-layer restratification parameterization', &
2054	1	'm', conversion=GV%H_to_m)
2055		CS%id_Rml = register_diag_field('ocean_model', 'ML_buoy_restrat', diag%axesT1, Time, &
2056		'Mixed Layer Buoyancy as used in the mixed-layer restratification parameterization', &
2057	1	'm s-2', conversion=GV%m_to_H(US%L_T_to_m_s*2))
2058		CS%id_uDml = register_diag_field('ocean_model', 'udml_restrat', diag%axesCu1, Time, &
2059		'Transport stream function amplitude for zonal restratification of mixed layer', &
2060	1	'm3 s-1', conversion=GV%H_to_m(US%L_to_m2)US%s_to_T)
2061		CS%id_vDml = register_diag_field('ocean_model', 'vdml_restrat', diag%axesCv1, Time, &
2062		'Transport stream function amplitude for meridional restratification of mixed layer', &
2063	1	'm3 s-1', conversion=GV%H_to_m(US%L_to_m2)US%s_to_T)
2064		CS%id_uml = register_diag_field('ocean_model', 'uml_restrat', diag%axesCu1, Time, &
2065		'Surface zonal velocity component of mixed layer restratification', &
2066	1	'm s-1', conversion=US%L_T_to_m_s)
2067		CS%id_vml = register_diag_field('ocean_model', 'vml_restrat', diag%axesCv1, Time, &
2068		'Surface meridional velocity component of mixed layer restratification', &
2069	1	'm s-1', conversion=US%L_T_to_m_s)
2070	1	if (CS%use_Bodner) then
2071		CS%id_wpup = register_diag_field('ocean_model', 'MLE_wpup', diag%axesT1, Time, &
2072		'Vertical turbulent momentum flux in Bodner mixed layer restratification parameterization', &
2073	1	'm2 s-2', conversion=US%L_to_mGV%H_to_mUS%s_to_T**2)
2074		CS%id_ustar = register_diag_field('ocean_model', 'MLE_ustar', diag%axesT1, Time, &
2075		'Surface turbulent friction velocity, u*, in Bodner mixed layer restratification parameterization', &
2076	1	'm s-1', conversion=(US%Z_to_m*US%s_to_T))
2077		CS%id_bflux = register_diag_field('ocean_model', 'MLE_bflux', diag%axesT1, Time, &
2078		'Surface buoyancy flux, B0, in Bodner mixed layer restratification parameterization', &
2079	1	'm2 s-3', conversion=(US%Z_to_m*2US%s_to_T**3))
2080		CS%id_lfbod = register_diag_field('ocean_model', 'lf_bodner', diag%axesT1, Time, &
2081		'Front length in Bodner mixed layer restratificiation parameterization', &
2082	1	'm', conversion=US%L_to_m)
2083		else
2084		CS%id_mle_fl = register_diag_field('ocean_model', 'mle_fl', diag%axesT1, Time, &
2085		'Frontal length scale used in the mixed layer restratificiation parameterization', &
2086	0	'm', conversion=US%L_to_m)
2087		endif
2088
2089		! If MLD_filtered is being used, we need to update halo regions after a restart
2090	1	if (allocated(CS%MLD_filtered)) call pass_var(CS%MLD_filtered, G%domain)
2091	1	if (allocated(CS%MLD_filtered_slow)) call pass_var(CS%MLD_filtered_slow, G%domain)
2092	1	if (allocated(CS%wpup_filtered)) call pass_var(CS%wpup_filtered, G%domain)
2093
2094		if (CS%use_Bodner) then
2095		! very important!
2096		!$omp target update to(CS)
2097		!$omp target enter data map(to: CS%Cr_space)
2098		!$omp target enter data map(to: CS%MLD_filtered, CS%MLD_filtered_slow, CS%wpup_filtered)
2099		endif
2100
2101	1	end function mixedlayer_restrat_init
2102
2103		!> Allocate and register fields in the mixed layer restratification structure for restarts
2104	1	subroutine mixedlayer_restrat_register_restarts(HI, GV, US, param_file, CS, restart_CS)
2105		! Arguments
2106		type(hor_index_type), intent(in) :: HI !< Horizontal index structure
2107		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
2108		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2109		type(param_file_type), intent(in) :: param_file !< Parameter file to parse
2110		type(mixedlayer_restrat_CS), intent(inout) :: CS !< Module control structure
2111		type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure
2112
2113		! Local variables
2114		character(len=64) :: mom_flux_units
2115		logical :: mixedlayer_restrat_init, use_Bodner
2116
2117		! Check to see if this module will be used
2118		call get_param(param_file, mdl, "MIXEDLAYER_RESTRAT", mixedlayer_restrat_init, &
2119	1	default=.false., do_not_log=.true.)
2120	1	if (.not. mixedlayer_restrat_init) return
2121
2122		call get_param(param_file, mdl, "MLE_MLD_DECAY_TIME", CS%MLE_MLD_decay_time, &
2123	1	units="s", default=0., scale=US%s_to_T, do_not_log=.true.)
2124		call get_param(param_file, mdl, "MLE_MLD_DECAY_TIME2", CS%MLE_MLD_decay_time2, &
2125	1	units="s", default=0., scale=US%s_to_T, do_not_log=.true.)
2126	1	call openParameterBlock(param_file, 'MLE', do_not_log=.true.)
2127		call get_param(param_file, mdl, "USE_BODNER23", use_Bodner, &
2128	1	default=.false., do_not_log=.true.)
2129	1	call closeParameterBlock(param_file)
2130	1	if (CS%MLE_MLD_decay_time>0. .or. CS%MLE_MLD_decay_time2>0. .or. use_Bodner) then
2131		! CS%MLD_filtered is used to keep a running mean of the PBL's actively mixed MLD.
2132	8773	allocate(CS%MLD_filtered(HI%isd:HI%ied,HI%jsd:HI%jed), source=0.)
2133		call register_restart_field(CS%MLD_filtered, "MLD_MLE_filtered", .false., restart_CS, &
2134		longname="Time-filtered MLD for use in MLE", &
2135	1	units=get_thickness_units(GV), conversion=GV%H_to_MKS)
2136		endif
2137	1	if (CS%MLE_MLD_decay_time2>0. .or. use_Bodner) then
2138		! CS%MLD_filtered_slow is used to keep a running mean of the PBL's seasonal or winter MLD.
2139	8773	allocate(CS%MLD_filtered_slow(HI%isd:HI%ied,HI%jsd:HI%jed), source=0.)
2140		call register_restart_field(CS%MLD_filtered_slow, "MLD_MLE_filtered_slow", .false., restart_CS, &
2141		longname="Slower time-filtered MLD for use in MLE", &
2142	1	units=get_thickness_units(GV), conversion=GV%H_to_MKS)
2143		endif
2144	1	if (use_Bodner) then
2145		! CS%MLD_filtered_slow is used to keep a running mean of the PBL's seasonal or winter MLD.
2146	1	mom_flux_units = "m2 s-2" ; if (.not.GV%Boussinesq) mom_flux_units = "kg m-1 s-2"
2147	8773	allocate(CS%wpup_filtered(HI%isd:HI%ied,HI%jsd:HI%jed), source=0.)
2148		call register_restart_field(CS%wpup_filtered, "MLE_Bflux", .false., restart_CS, &
2149		longname="Time-filtered vertical turbulent momentum flux for use in MLE", &
2150	1	units=mom_flux_units, conversion=US%L_to_mGV%H_to_mksUS%s_to_T**2 )
2151		endif
2152
2153		end subroutine mixedlayer_restrat_register_restarts
2154
2155		!> Returns true if a unit test of functions in MOM_mixedlayer_restrat fail.
2156		!! Returns false otherwise.
2157	0	logical function mixedlayer_restrat_unit_tests(verbose)
2158		logical, intent(in) :: verbose !< If true, write results to stdout
2159
2160		! Local variables
2161		logical :: this_test
2162
2163	0	print *,'===== mixedlayer_restrat: mixedlayer_restrat_unit_tests =================='
2164
2165		! Tests of the shape function mu(z)
2166		this_test = &
2167	0	test_answer(verbose, mu(3.,0.), 0., 'mu(3)=0')
2168		this_test = this_test .or. &
2169	0	test_answer(verbose, mu(0.,0.), 0., 'mu(0)=0')
2170		this_test = this_test .or. &
2171	0	test_answer(verbose, mu(-0.25,0.), 0.7946428571428572, 'mu(-0.25)=0.7946...', tol=epsilon(1.))
2172		this_test = this_test .or. &
2173	0	test_answer(verbose, mu(-0.5,0.), 1., 'mu(-0.5)=1')
2174		this_test = this_test .or. &
2175	0	test_answer(verbose, mu(-0.75,0.), 0.7946428571428572, 'mu(-0.75)=0.7946...', tol=epsilon(1.))
2176		this_test = this_test .or. &
2177	0	test_answer(verbose, mu(-1.,0.), 0., 'mu(-1)=0')
2178		this_test = this_test .or. &
2179	0	test_answer(verbose, mu(-3.,0.), 0., 'mu(-3)=0')
2180		this_test = this_test .or. &
2181	0	test_answer(verbose, mu(-0.5,0.5), 1., 'mu(-0.5,0.5)=1')
2182		this_test = this_test .or. &
2183	0	test_answer(verbose, mu(-1.,0.5), 0.25, 'mu(-1,0.5)=0.25')
2184		this_test = this_test .or. &
2185	0	test_answer(verbose, mu(-1.5,0.5), 0., 'mu(-1.5,0.5)=0')
2186	0	if (.not. this_test) print '(a)',' Passed tests of mu(z)'
2187	0	mixedlayer_restrat_unit_tests = this_test
2188
2189		! Tests of the two time-scale running mean function
2190		this_test = &
2191	0	test_answer(verbose, rmean2ts(3.,2.,0.,0.,3.), 3., 'rmean2ts(3,2,0,0,3)=3')
2192		this_test = this_test .or. &
2193	0	test_answer(verbose, rmean2ts(1.,2.,0.,0.,3.), 1., 'rmean2ts(1,2,0,0,3)=1')
2194		this_test = this_test .or. &
2195	0	test_answer(verbose, rmean2ts(4.,0.,3.,0.,1.), 1., 'rmean2ts(4,0,3,0,1)=1')
2196		this_test = this_test .or. &
2197	0	test_answer(verbose, rmean2ts(0.,4.,0.,3.,1.), 3., 'rmean2ts(0,4,0,3,1)=3')
2198	0	if (.not. this_test) print '(a)',' Passed tests of rmean2ts(s,f,g,d,dt)'
2199	0	mixedlayer_restrat_unit_tests = mixedlayer_restrat_unit_tests .or. this_test
2200
2201	0	end function mixedlayer_restrat_unit_tests
2202
2203		!> Returns true if any cell of u and u_true are not identical. Returns false otherwise.
2204	0	logical function test_answer(verbose, u, u_true, label, tol)
2205		logical, intent(in) :: verbose !< If true, write results to stdout
2206		real, intent(in) :: u !< Values to test in arbitrary units [A]
2207		real, intent(in) :: u_true !< Values to test against (correct answer) [A]
2208		character(len=*), intent(in) :: label !< Message
2209		real, optional, intent(in) :: tol !< The tolerance for differences between u and u_true [A]
2210		! Local variables
2211		real :: tolerance ! The tolerance for differences between u and u_true [A]
2212
2213	0	tolerance = 0.0 ; if (present(tol)) tolerance = tol
2214	0	test_answer = .false.
2215
2216	0	if (abs(u - u_true) > tolerance) test_answer = .true.
2217	0	if (test_answer .or. verbose) then
2218	0	if (test_answer) then
2219	0	print '(3(a,1pe24.16),1x,a,1x,a)','computed =',u,' correct =',u_true, &
2220	0	' err=',u-u_true,' < wrong',label
2221		else
2222	0	print '(2(a,1pe24.16),1x,a)','computed =',u,' correct =',u_true,label
2223		endif
2224		endif
2225
2226	0	end function test_answer
2227
2228		!> \namespace mom_mixed_layer_restrat
2229		!!
2230		!! \section section_mle Mixed-layer eddy parameterization module
2231		!!
2232		!! The subroutines in this module implement a parameterization of unresolved viscous
2233		!! mixed layer restratification of the mixed layer as described in \cite fox-kemper2008,
2234		!! and whose impacts are described in \cite fox-kemper2011.
2235		!! This is derived in part from the older parameterization that is described in
2236		!! \cite Hallberg2003, which this new parameterization surpasses, which
2237		!! in turn is based on the sub-inertial mixed layer theory of \cite Young1994.
2238		!! There is no net horizontal volume transport due to this parameterization, and
2239		!! no direct effect below the mixed layer. A revised version of the parameterization by
2240		!! \cite Bodner2023 is also available as an option.
2241		!!
2242		!! This parameterization sets the restratification timescale to agree with
2243		!! high-resolution studies of mixed layer restratification.
2244		!!
2245		!! The run-time parameter <code>FOX_KEMPER_ML_RESTRAT_COEF</code> is a non-dimensional number of
2246		!! order a few tens, proportional to the ratio of the deformation radius or the
2247		!! grid scale (whichever is smaller) to the dominant horizontal length-scale of the
2248		!! sub-meso-scale mixed layer instabilities.
2249		!!
2250		!! \subsection section_mle_nutshell "Sub-meso" in a nutshell
2251		!!
2252		!! The parameterization is colloquially referred to as "sub-meso".
2253		!!
2254		!! The original \cite fox-kemper2008-2 paper proposed a quasi-Stokes
2255		!! advection described by the stream function (eq. 5 of \cite fox-kemper2011):
2256		!! \f[
2257		!! {\bf \Psi}_o = C_e \frac{ H^2 \nabla \bar{b} \times \hat{\bf z} }{ \|f\| } \mu(z)
2258		!! \f]
2259		!!
2260		!! where the vertical profile function is
2261		!! \f[
2262		!! \mu(z) = \max \left\{ 0, \left[ 1 - \left(\frac{2z}{H}+1\right)^2 \right]
2263		!! \left[ 1 + \frac{5}{21} \left(\frac{2z}{H}+1\right)^2 \right] \right\}
2264		!! \f]
2265		!! and \f$ H \f$ is the mixed-layer depth, \f$ f \f$ is the local Coriolis parameter, \f$ C_e \sim 0.06-0.08 \f$ and
2266		!! \f$ \nabla \bar{b} \f$ is a depth mean buoyancy gradient averaged over the mixed layer.
2267		!!
2268		!! For use in coarse-resolution models, an upscaling of the buoyancy gradients and adaption for the equator
2269		!! leads to the following parameterization (eq. 6 of \cite fox-kemper2011):
2270		!! \f[
2271		!! {\bf \Psi} = C_e \Gamma_\Delta \frac{\Delta s}{l_f} \frac{ H^2 \nabla \bar{b} \times \hat{\bf z} }
2272		!! { \sqrt{ f^2 + \tau^{-2}} } \mu(z)
2273		!! \f]
2274		!! where \f$ \Delta s \f$ is the minimum of grid-scale and deformation radius,
2275		!! \f$ l_f \f$ is the width of the mixed-layer fronts, and \f$ \Gamma_\Delta=1 \f$.
2276		!! \f$ \tau \f$ is a time-scale for mixing momentum across the mixed layer.
2277		!! \f$ l_f \f$ is thought to be of order hundreds of meters.
2278		!!
2279		!! The upscaling factor \f$ \frac{\Delta s}{l_f} \f$ can be a global constant, model parameter
2280		!! <code>FOX_KEMPER_ML_RESTRAT</code>, so that in practice the parameterization is:
2281		!! \f[
2282		!! {\bf \Psi} = C_e \Gamma_\Delta \frac{ H^2 \nabla \bar{b} \times \hat{\bf z} }{ \sqrt{ f^2 + \tau^{-2}} } \mu(z)
2283		!! \f]
2284		!! with non-unity \f$ \Gamma_\Delta \f$.
2285		!!
2286		!! \f$ C_e \f$ is hard-coded as 0.0625. \f$ \tau \f$ is calculated from the surface friction velocity \f$ u^* \f$.
2287		!!
2288		!! \todo Explain expression for momentum mixing time-scale.
2289		!!
2290		!! \| Symbol \| Module parameter \|
2291		!! \| ---------------------------- \| --------------------- \|
2292		!! \| \f$ \Gamma_\Delta \f$ \| FOX_KEMPER_ML_RESTRAT \|
2293		!! \| \f$ l_f \f$ \| MLE_FRONT_LENGTH \|
2294		!! \| \f$ \Delta \rho \f$ \| MLE_DENSITY_DIFF \|
2295		!!
2296		!! \subsection section_mle_filtering Time-filtering of mixed-layer depth
2297		!!
2298		!! Using the instantaneous mixed-layer depth is inconsistent with the finite life-time of
2299		!! mixed-layer instabilities. We provide a one-sided running-mean filter of mixed-layer depth, \f$ H \f$, of the form:
2300		!! \f[
2301		!! \bar{H} \leftarrow \max \left( H, \frac{ \Delta t H + \tau_h \bar{H} }{ \Delta t + \tau_h } \right)
2302		!! \f]
2303		!! which allows the effective mixed-layer depth seen by the parameterization, \f$\bar{H}\f$, to instantaneously deepen
2304		!! but to decay with time-scale \f$ \tau_h \f$.
2305		!! \f$ \bar{H} \f$ is substituted for \f$ H \f$ in the above equations.
2306		!!
2307		!! \| Symbol \| Module parameter \|
2308		!! \| ---------------------------- \| --------------------- \|
2309		!! \| \f$ \tau_h \f$ \| MLE_MLD_DECAY_TIME \|
2310		!!
2311		!! \subsection section_mle_mld Defining the mixed-layer-depth
2312		!!
2313		!! If the parameter MLE_USE_PBL_MLD=True then the mixed-layer depth is defined/diagnosed by the
2314		!! boundary-layer parameterization (e.g. ePBL, KPP, etc.).
2315		!!
2316		!! If the parameter MLE_USE_PBL_MLD=False then the mixed-layer depth is diagnosed in this module
2317		!! as the depth of a given density difference, \f$ \Delta \rho \f$, with the surface where the
2318		!! density difference is the parameter MLE_DENSITY_DIFF.
2319		!!
2320		!! \subsection The Bodner (2023) modification
2321		!!
2322		!! To use this variant of the parameterization, set MLE\%USE_BODNER23=True which then changes the
2323		!! available parameters.
2324		!! MLE_USE_PBL_MLD must be True to use the B23 modification.
2325		!!
2326		!! \cite Bodner2023, (B23) use an expression for the frontal width which changes the scaling from \f$ H^2 \f$
2327		!! to \f$ h H^2 \f$:
2328		!! \f[
2329		!! {\bf \Psi} = C_r \frac{\Delta s \|f\| \bar{h} \bar{H}^2 \nabla \bar{b} \times \hat{\bf z} }
2330		!! { \left( m_u_^3 + n_* w_*^3 \right)^{2/3} } \mu(z)
2331		!! \f]
2332		!! (see eq. 27 of B23).
2333		!! Here, the \f$h\f$ is the activate boundary layer depth, and \f$H\f$ is the mixed layer depth.
2334		!! The denominator is an approximation of the vertical turbulent momentum flux \f$\overline{w'u'}\f$ (see
2335		!! eq. 18 of B23) calculated from the surface friction velocity \f$u_*\f$, and from the surface buoyancy flux,
2336		!! \f$B\f$, using the relation \f$ w_*^3 \sim -B h \f$.
2337		!! An advantage of this form of "sub-meso" is the denominator is well behaved at the equator but we apply a
2338		!! lower bound of \f$w_{min}^2\f$ to avoid division by zero under zero forcing.
2339		!! As for the original Fox-Kemper parameterization, \f$\nabla \bar{b}\f$ is the buoyancy gradient averaged
2340		!! over the mixed-layer.
2341		!!
2342		!! The instantaneous boundary layer depth, \f$h\f$, is time filtered primarily to remove the diurnal cycle:
2343		!! \f[
2344		!! \bar{h} \leftarrow \max \left(
2345		!! \min \left( h, \frac{ \Delta t h + \tau_{h+} \bar{h} }{ \Delta t + \tau_{h+} } \right),
2346		!! \frac{ \Delta t h + \tau_{h-} \bar{h} }{ \Delta t + \tau_{h-} } \right)
2347		!! \f]
2348		!! Setting \f$ \tau_{h+}=0 \f$ means that when \f$ h>\bar{h} \f$ then \f$\bar{h}\leftarrow h\f$, i.e. the
2349		!! effective (filtered) depth, \f$\bar{h}\f$, is instantly deepened. When \f$h<\bar{h}\f$ then the effective
2350		!! depth shoals with time-scale \f$\tau_{h-}\f$.
2351		!!
2352		!! A second filter is applied to \f$\bar{h}\f$ to yield and effective "mixed layer depth", \f$\bar{H}\f$,
2353		!! defined as the deepest the boundary layer over some time-scale \f$\tau_{H-}\f$:
2354		!! \f[
2355		!! \bar{H} \leftarrow \max \left(
2356		!! \min \left( \bar{h}, \frac{ \Delta t \bar{h} + \tau_{H+} \bar{H} }{ \Delta t + \tau_{H+} } \right),
2357		!! \frac{ \Delta t \bar{h} + \tau_{h-} \bar{H} }{ \Delta t + \tau_{H-} } \right)
2358		!! \f]
2359		!! Again, setting \f$ \tau_{H+}=0 \f$ allows the effective mixed layer to instantly deepend to \f$ \bar{h} \f$.
2360		!!
2361		!! \| Symbol \| Module parameter \|
2362		!! \| ---------------------------- \| ------------------------- \|
2363		!! \| \f$ C_r \f$ \| MLE\%CR \|
2364		!! \| \f$ n_* \f$ \| MLE\%BODNER_NSTAR \|
2365		!! \| \f$ m_* \f$ \| MLE\%BODNER_MSTAR \|
2366		!! \| \f$ w_* \f$ \| MLE\%BODNER_MSTAR \|
2367		!! \| \f$ w_{min}^2 \f$ \| MLE\%MIN_WSTAR2 \|
2368		!! \| \f$ \tau_{h+} \f$ \| MLE\%BLD_GROWING_TFILTER \|
2369		!! \| \f$ \tau_{h-} \f$ \| MLE\%BLD_DECAYING_TFILTER \|
2370		!! \| \f$ \tau_{H+} \f$ \| MLE\%MLD_GROWING_TFILTER \|
2371		!! \| \f$ \tau_{H-} \f$ \| MLE\%MLD_DECAYING_TFILTER \|
2372		!!
2373		!! \subsection section_mle_ref References
2374		!!
2375		!! Fox-Kemper, B., Ferrari, R. and Hallberg, R., 2008:
2376		!! Parameterization of Mixed Layer Eddies. Part I: Theory and Diagnosis
2377		!! J. Phys. Oceangraphy, 38 (6), p1145-1165.
2378		!! https://doi.org/10.1175/2007JPO3792.1
2379		!!
2380		!! Fox-Kemper, B. and Ferrari, R. 2008:
2381		!! Parameterization of Mixed Layer Eddies. Part II: Prognosis and Impact
2382		!! J. Phys. Oceangraphy, 38 (6), p1166-1179.
2383		!! https://doi.org/10.1175/2007JPO3788.1
2384		!!
2385		!! B. Fox-Kemper, G. Danabasoglu, R. Ferrari, S.M. Griffies, R.W. Hallberg, M.M. Holland, M.E. Maltrud,
2386		!! S. Peacock, and B.L. Samuels, 2011: Parameterization of mixed layer eddies. III: Implementation and impact
2387		!! in global ocean climate simulations. Ocean Modell., 39(1), p61-78.
2388		!! https://doi.org/10.1016/j.ocemod.2010.09.002
2389		!!
2390		!! A.S. Bodner, B. Fox-Kemper, L. Johnson, L. P. Van Roekel, J. C. McWilliams, P. P. Sullivan, P. S. Hall,
2391		!! and J. Dong, 2023: Modifying the Mixed Layer Eddy Parameterization to Include Frontogenesis Arrest by
2392		!! Boundary Layer Turbulence. J. Phys. Oceanogr., 53(1), p323-339.
2393		!! https://doi.org/10.1175/JPO-D-21-0297.1
2394
2395	0	end module MOM_mixed_layer_restrat