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