MOM_mixed_layer_restrat.F90

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.
7
8use mom_debugging, only : hchksum
10use mom_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type
12use mom_domains, only : pass_var, to_west, to_south, omit_corners
13use mom_error_handler, only : mom_error, fatal, warning
14use mom_file_parser, only : get_param, log_param, log_version, param_file_type
16use mom_forcing_type, only : mech_forcing, find_ustar
17use mom_grid, only : ocean_grid_type
19use mom_interpolate, only : init_external_field, time_interp_external, time_interp_external_init
20use mom_interpolate, only : external_field
22use mom_io, only : slasher, mom_read_data
24use mom_restart, only : register_restart_field, query_initialized, mom_restart_cs
28use mom_eos, only : calculate_density, calculate_spec_vol, eos_domain
29
30implicit none ; private
31
32#include <MOM_memory.h>
33
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
45type, 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
143
144character(len=40) :: mdl = "MOM_mixed_layer_restrat" !< This module's name.
145
146contains
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.
151subroutine 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 if (.not. cs%initialized) call mom_error(fatal, "mixedlayer_restrat: "// &
175 "Module must be initialized before it is used.")
176
177 if (gv%nkml>0) then
178 ! Original form, written for the isopycnal model with a bulk mixed layer
179 call mixedlayer_restrat_bml(h, uhtr, vhtr, tv, forces, dt, g, gv, us, cs)
180 elseif (cs%use_Bodner) then
181 ! Implementation of Bodner et al., 2023
182 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 call mixedlayer_restrat_om4(h, uhtr, vhtr, tv, forces, dt, h_mld, varmix, g, gv, us, cs)
186 endif
187
188end subroutine mixedlayer_restrat
189
190!> Calculates a restratifying flow in the mixed layer, following the formulation used in OM4
191subroutine 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 real :: uhml(SZIB_(G),SZJ_(G),SZK_(GV)) ! Restratifying zonal thickness transports [H L2 T-1 ~> m3 s-1 or kg s-1]
212 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 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 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 mld_fast, & ! Mixed layer depth actually used in MLE restratification parameterization [H ~> m or kg m-2]
221 htot_fast, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2]
222 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 mld_slow, & ! Mixed layer depth actually used in MLE restratification parameterization [H ~> m or kg m-2]
225 mle_fl_2d, & ! MLE frontal length-scale [L ~> m]
226 htot_slow, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2]
227 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 real :: rho_ml(SZI_(G)) ! Potential density relative to the surface [R ~> kg m-3]
232 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 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 real :: SpV_ml(SZI_(G)) ! Specific volume evaluated at the surface pressure [R-1 ~> m3 kg-1]
235 real :: SpV_int_fast(SZI_(G)) ! Specific volume integrated through the mixed layer [H R-1 ~> m4 kg-1 or m]
236 real :: SpV_int_slow(SZI_(G)) ! Specific volume integrated through the mixed layer [H R-1 ~> m4 kg-1 or m]
237 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 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 real :: b(SZK_(GV)) ! As for a(k) but for the slow-filtered MLD [nondim]
254 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 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 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 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 real :: utimescale_diag(SZIB_(G),SZJ_(G)) ! Zonal restratification timescale [T ~> s], stored for diagnostics.
261 real :: vtimescale_diag(SZI_(G),SZJB_(G)) ! Meridional restratification timescale [T ~> s], stored for diagnostics.
262 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 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 real, dimension(SZI_(G)) :: covTS, & ! SGS TS covariance in Stanley param; currently 0 [C S ~> degC ppt]
265 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 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
282 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
283
284 h_min = 0.5*gv%Angstrom_H ! This should be GV%Angstrom_H, but that value would change answers.
285 covts(:) = 0.0 !!Functionality not implemented yet; in future, should be passed in tv
286 vars(:) = 0.0
287 mle_fl_2d(:,:) = 0.0
288 vonkar_x_pi2 = cs%vonKar * 9.8696
289
290 if (.not.associated(tv%eqn_of_state)) call mom_error(fatal, "mixedlayer_restrat_OM4: "// &
291 "An equation of state must be used with this module.")
292 if (.not. allocated(varmix%Rd_dx_h) .and. cs%front_length > 0.) &
293 call mom_error(fatal, "mixedlayer_restrat_OM4: "// &
294 "The resolution argument, Rd/dx, was not associated.")
295 if (cs%use_Stanley_ML .and. .not.gv%Boussinesq) call mom_error(fatal, &
296 "MOM_mixedlayer_restrat: The Stanley parameterization is not "//&
297 "available without the Boussinesq approximation.")
298
299 ! Extract the friction velocity from the forcing type.
300 call find_ustar(forces, tv, u_star_2d, g, gv, us, halo=1, h_t_units=.true.)
301
302 if (cs%MLE_density_diff > 0.) then ! We need to calculate a mixed layer depth, MLD.
303 call detect_mld(h, tv, mld_fast, g, gv, cs)
304 elseif (cs%MLE_use_PBL_MLD) then
305 do j=js-1,je+1 ; do i=is-1,ie+1
306 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 "No MLD to use for MLE parameterization.")
311 endif
312
313 ! Apply time filter (to remove diurnal cycle)
314 if (cs%MLE_MLD_decay_time>0.) then
315 if (cs%debug) then
316 call hchksum(cs%MLD_filtered, 'mixed_layer_restrat: MLD_filtered', g%HI, haloshift=1, unscale=gv%H_to_mks)
317 call hchksum(h_mld, 'mixed_layer_restrat: MLD in', g%HI, haloshift=1, unscale=gv%H_to_mks)
318 endif
319 afac = cs%MLE_MLD_decay_time / ( dt + cs%MLE_MLD_decay_time )
320 bfac = dt / ( dt + cs%MLE_MLD_decay_time )
321 do j=js-1,je+1 ; do i=is-1,ie+1
322 ! Expression bFac*MLD_fast(i,j) + aFac*CS%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 cs%MLD_filtered(i,j) = max( mld_fast(i,j), bfac*mld_fast(i,j) + afac*cs%MLD_filtered(i,j) )
326 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 if (cs%MLE_MLD_decay_time2>0.) then
332 if (cs%debug) then
333 call hchksum(cs%MLD_filtered_slow, 'mixed_layer_restrat: MLD_filtered_slow', g%HI, &
334 haloshift=1, unscale=gv%H_to_mks)
335 call hchksum(mld_fast, 'mixed_layer_restrat: MLD fast', g%HI, haloshift=1, unscale=gv%H_to_mks)
336 endif
337 afac = cs%MLE_MLD_decay_time2 / ( dt + cs%MLE_MLD_decay_time2 )
338 bfac = dt / ( dt + cs%MLE_MLD_decay_time2 )
339 do j=js-1,je+1 ; do i=is-1,ie+1
340 ! Expression bFac*MLD_fast(i,j) + aFac*CS%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 cs%MLD_filtered_slow(i,j) = max( mld_fast(i,j), bfac*mld_fast(i,j) + afac*cs%MLD_filtered_slow(i,j) )
344 mld_slow(i,j) = cs%MLD_filtered_slow(i,j)
345 enddo ; enddo
346 else
347 do j=js-1,je+1 ; do i=is-1,ie+1
348 mld_slow(i,j) = mld_fast(i,j)
349 enddo ; enddo
350 endif
351
352 udml(:) = 0.0 ; vdml(:) = 0.0
353 udml_slow(:) = 0.0 ; vdml_slow(:) = 0.0
354 i4dt = 0.25 / dt
355 g_rho0 = gv%H_to_Z * gv%g_Earth / gv%Rho0
356 h_neglect = gv%H_subroundoff
357 if (cs%front_length>0.) then
358 res_upscale = .true.
359 do j=js-1,je+1 ; do i=is-1,ie+1
360 mle_fl_2d(i,j) = cs%front_length
361 enddo ; enddo
362 elseif (cs%front_length == 0. .and. cs%fl_from_file) then
363 res_upscale = .true.
364 call time_interp_external(cs%sbc_fl, cs%Time, mle_fl_2d, turns=g%HI%turns, scale=us%m_to_L)
365 call pass_var(mle_fl_2d, g%domain, halo=1)
366 do j=js,je ; do i=is,ie
367 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 mle_fl_2d(i,j), i, j, g%geoLonT(i,j), g%geoLatT(i,j)
371 call mom_error(fatal, "MOM_mixed_layer_restrat mixedlayer_restrat_OM4: "//trim(mesg))
372 endif
373 enddo ; enddo
374 else
375 res_upscale = .false.
376 endif
377
378 p0(:) = 0.0
379 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 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
387 !$OMP do
388 do j=js-1,je+1
389 do i=is-1,ie+1
390 htot_fast(i,j) = 0.0 ; rml_int_fast(i) = 0.0
391 htot_slow(i,j) = 0.0 ; rml_int_slow(i) = 0.0
392 enddo
393 keep_going = .true.
394 do k=1,nz
395 do i=is-1,ie+1
396 h_avail(i,j,k) = max(i4dt*g%areaT(i,j)*(h(i,j,k)-gv%Angstrom_H),0.0)
397 enddo
398 if (keep_going) then
399 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 rho_ml(:), tv%eqn_of_state, eosdom)
402 else
403 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, rho_ml(:), tv%eqn_of_state, eosdom)
404 endif
405 line_is_empty = .true.
406 do i=is-1,ie+1
407 if (htot_fast(i,j) < mld_fast(i,j)) then
408 dh = min( h(i,j,k), mld_fast(i,j)-htot_fast(i,j) )
409 rml_int_fast(i) = rml_int_fast(i) + dh*rho_ml(i)
410 htot_fast(i,j) = htot_fast(i,j) + dh
411 line_is_empty = .false.
412 endif
413 if (htot_slow(i,j) < mld_slow(i,j)) then
414 dh = min( h(i,j,k), mld_slow(i,j)-htot_slow(i,j) )
415 rml_int_slow(i) = rml_int_slow(i) + dh*rho_ml(i)
416 htot_slow(i,j) = htot_slow(i,j) + dh
417 line_is_empty = .false.
418 endif
419 enddo
420 if (line_is_empty) keep_going=.false.
421 endif
422 enddo
423
424 do i=is-1,ie+1
425 rml_av_fast(i,j) = -(g_rho0*rml_int_fast(i)) / (htot_fast(i,j) + h_neglect)
426 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 do j=js-1,je+1
432 do i=is-1,ie+1
433 htot_fast(i,j) = 0.0 ; spv_int_fast(i) = 0.0
434 htot_slow(i,j) = 0.0 ; spv_int_slow(i) = 0.0
435 enddo
436 keep_going = .true.
437 do k=1,nz
438 do i=is-1,ie+1
439 h_avail(i,j,k) = max(i4dt*g%areaT(i,j)*(h(i,j,k)-gv%Angstrom_H),0.0)
440 enddo
441 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 call calculate_spec_vol(tv%T(:,j,k), tv%S(:,j,k), p0, spv_ml, tv%eqn_of_state, eosdom)
447 ! endif
448 line_is_empty = .true.
449 do i=is-1,ie+1
450 if (htot_fast(i,j) < mld_fast(i,j)) then
451 dh = min( h(i,j,k), mld_fast(i,j)-htot_fast(i,j) )
452 spv_int_fast(i) = spv_int_fast(i) + dh*spv_ml(i)
453 htot_fast(i,j) = htot_fast(i,j) + dh
454 line_is_empty = .false.
455 endif
456 if (htot_slow(i,j) < mld_slow(i,j)) then
457 dh = min( h(i,j,k), mld_slow(i,j)-htot_slow(i,j) )
458 spv_int_slow(i) = spv_int_slow(i) + dh*spv_ml(i)
459 htot_slow(i,j) = htot_slow(i,j) + dh
460 line_is_empty = .false.
461 endif
462 enddo
463 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 do i=is-1,ie+1
469 rml_av_fast(i,j) = (gv%H_to_RZ*gv%g_Earth * spv_int_fast(i)) / (htot_fast(i,j) + h_neglect)
470 rml_av_slow(i,j) = (gv%H_to_RZ*gv%g_Earth * spv_int_slow(i)) / (htot_slow(i,j) + h_neglect)
471 enddo
472 enddo
473 endif
474
475 if (cs%debug) then
476 call hchksum(h, 'mixed_layer_restrat: h', g%HI, haloshift=1, unscale=gv%H_to_mks)
477 call hchksum(u_star_2d, 'mixed_layer_restrat: u*', g%HI, haloshift=1, unscale=gv%H_to_m*us%s_to_T)
478 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 unscale=gv%m_to_H*us%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 do j=js,je ; do i=is-1,ie
490 u_star = max(cs%ustar_min, 0.5*(u_star_2d(i,j) + u_star_2d(i+1,j)))
491
492 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 lfront = 0.5 * (mle_fl_2d(i,j) + mle_fl_2d(i+1,j))
495 ! Adcroft reciprocal
496 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 if (res_upscale) res_scaling_fac = &
499 ( sqrt( 0.5 * ( (g%dxCu(i,j)**2) + (g%dyCu(i,j)**2) ) ) * i_lfront ) &
500 * 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_ml*u_star)/(absf*h_ml + 4.0*u_star)
503 ! momentum mixing rate: pi^2*visc/h_ml^2
504 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_pi2*u_star**2 / &
510 (absf*h_vel**2 + 4.0*(h_vel+h_neglect)*u_star)
511 timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
512 timescale = timescale * cs%ml_restrat_coef
513
514 if (res_upscale) timescale = timescale * res_scaling_fac
515 udml(i) = timescale * g%dyCu(i,j)*g%IdxCu_OBCmask(i,j) * &
516 (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 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_pi2*u_star**2 / &
525 (absf*h_vel**2 + 4.0*(h_vel+h_neglect)*u_star)
526 timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
527 timescale = timescale * cs%ml_restrat_coef2
528
529 if (res_upscale) timescale = timescale * res_scaling_fac
530 udml_slow(i) = timescale * g%dyCu(i,j)*g%IdxCu_OBCmask(i,j) * &
531 (rml_av_slow(i+1,j)-rml_av_slow(i,j)) * (h_vel**2)
532
533 if (udml(i) + udml_slow(i) == 0.) then
534 do k=1,nz ; uhml(i,j,k) = 0.0 ; enddo
535 else
536 ihtot = 2.0 / ((htot_fast(i,j) + htot_fast(i+1,j)) + h_neglect)
537 ihtot_slow = 2.0 / ((htot_slow(i,j) + htot_slow(i+1,j)) + h_neglect)
538 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 do k=1,nz
542 hatvel = 0.5*(h(i,j,k) + h(i+1,j,k))
543 a(k) = mu(zpa, cs%MLE_tail_dh) ! mu(z/MLD) for upper interface
544 zpa = zpa - (hatvel * ihtot) ! z/H for lower interface
545 a(k) = a(k) - mu(zpa, cs%MLE_tail_dh) ! Transport profile
546 ! Limit magnitude (uDml) if it would violate CFL
547 if (a(k)*udml(i) > 0.0) then
548 if (a(k)*udml(i) > h_avail(i,j,k)) udml(i) = h_avail(i,j,k) / a(k)
549 elseif (a(k)*udml(i) < 0.0) then
550 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 do k=1,nz
554 ! Transport for slow-filtered MLD
555 hatvel = 0.5*(h(i,j,k) + h(i+1,j,k))
556 b(k) = mu(zpb, cs%MLE_tail_dh) ! mu(z/MLD) for upper interface
557 zpb = zpb - (hatvel * ihtot_slow) ! z/H for lower interface
558 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 if (b(k)*udml_slow(i) > 0.0) then
561 if (b(k)*udml_slow(i) > h_avail(i,j,k) - a(k)*udml(i)) &
562 udml_slow(i) = max( 0., h_avail(i,j,k) - a(k)*udml(i) ) / b(k)
563 elseif (b(k)*udml_slow(i) < 0.0) then
564 if (-b(k)*udml_slow(i) > h_avail(i+1,j,k) + a(k)*udml(i)) &
565 udml_slow(i) = -max( 0., h_avail(i+1,j,k) + a(k)*udml(i) ) / b(k)
566 endif
567 enddo
568 do k=1,nz
569 uhml(i,j,k) = a(k)*udml(i) + b(k)*udml_slow(i)
570 uhtr(i,j,k) = uhtr(i,j,k) + uhml(i,j,k)*dt
571 enddo
572 endif
573
574 utimescale_diag(i,j) = timescale
575 udml_diag(i,j) = udml(i)
576 enddo ; enddo
577
578! V- component
579 !$OMP do
580 do j=js-1,je ; do i=is,ie
581 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 lfront = 0.5 * (mle_fl_2d(i,j) + mle_fl_2d(i,j+1))
584 ! Adcroft reciprocal
585 i_lfront = 0.0 ; if (lfront /= 0.0) i_lfront = 1.0/lfront
586 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 if (res_upscale) res_scaling_fac = &
589 ( sqrt( 0.5 * ( (g%dxCv(i,j)**2) + (g%dyCv(i,j)**2) ) ) * i_lfront ) &
590 * 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_ml*u_star)/(absf*h_ml + 4.0*u_star)
593 ! momentum mixing rate: pi^2*visc/h_ml^2
594 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_pi2*u_star**2 / &
600 (absf*h_vel**2 + 4.0*(h_vel+h_neglect)*u_star)
601 timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
602 timescale = timescale * cs%ml_restrat_coef
603
604 if (res_upscale) timescale = timescale * res_scaling_fac
605 vdml(i) = timescale * g%dxCv(i,j)*g%IdyCv_OBCmask(i,j) * &
606 (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 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_pi2*u_star**2 / &
615 (absf*h_vel**2 + 4.0*(h_vel+h_neglect)*u_star)
616 timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
617 timescale = timescale * cs%ml_restrat_coef2
618
619 if (res_upscale) timescale = timescale * res_scaling_fac
620 vdml_slow(i) = timescale * g%dxCv(i,j)*g%IdyCv_OBCmask(i,j) * &
621 (rml_av_slow(i,j+1)-rml_av_slow(i,j)) * (h_vel**2)
622
623 if (vdml(i) + vdml_slow(i) == 0.) then
624 do k=1,nz ; vhml(i,j,k) = 0.0 ; enddo
625 else
626 ihtot = 2.0 / ((htot_fast(i,j) + htot_fast(i,j+1)) + h_neglect)
627 ihtot_slow = 2.0 / ((htot_slow(i,j) + htot_slow(i,j+1)) + h_neglect)
628 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 do k=1,nz
632 hatvel = 0.5*(h(i,j,k) + h(i,j+1,k))
633 a(k) = mu(zpa, cs%MLE_tail_dh) ! mu(z/MLD) for upper interface
634 zpa = zpa - (hatvel * ihtot) ! z/H for lower interface
635 a(k) = a(k) - mu(zpa, cs%MLE_tail_dh) ! Transport profile
636 ! Limit magnitude (vDml) if it would violate CFL
637 if (a(k)*vdml(i) > 0.0) then
638 if (a(k)*vdml(i) > h_avail(i,j,k)) vdml(i) = h_avail(i,j,k) / a(k)
639 elseif (a(k)*vdml(i) < 0.0) then
640 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 do k=1,nz
644 ! Transport for slow-filtered MLD
645 hatvel = 0.5*(h(i,j,k) + h(i,j+1,k))
646 b(k) = mu(zpb, cs%MLE_tail_dh) ! mu(z/MLD) for upper interface
647 zpb = zpb - (hatvel * ihtot_slow) ! z/H for lower interface
648 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 if (b(k)*vdml_slow(i) > 0.0) then
651 if (b(k)*vdml_slow(i) > h_avail(i,j,k) - a(k)*vdml(i)) &
652 vdml_slow(i) = max( 0., h_avail(i,j,k) - a(k)*vdml(i) ) / b(k)
653 elseif (b(k)*vdml_slow(i) < 0.0) then
654 if (-b(k)*vdml_slow(i) > h_avail(i,j+1,k) + a(k)*vdml(i)) &
655 vdml_slow(i) = -max( 0., h_avail(i,j+1,k) + a(k)*vdml(i) ) / b(k)
656 endif
657 enddo
658 do k=1,nz
659 vhml(i,j,k) = a(k)*vdml(i) + b(k)*vdml_slow(i)
660 vhtr(i,j,k) = vhtr(i,j,k) + vhml(i,j,k)*dt
661 enddo
662 endif
663
664 vtimescale_diag(i,j) = timescale
665 vdml_diag(i,j) = vdml(i)
666 enddo ; enddo
667
668 !$OMP do
669 do j=js,je ; do k=1,nz ; do i=is,ie
670 h(i,j,k) = h(i,j,k) - dt*g%IareaT(i,j) * &
671 ((uhml(i,j,k) - uhml(i-1,j,k)) + (vhml(i,j,k) - vhml(i,j-1,k)))
672 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 if (cs%id_uhml > 0 .or. cs%id_vhml > 0) &
679 ! Remapped uhml and vhml require east/north halo updates of h
680 call pass_var(h, g%domain, to_west+to_south+omit_corners, halo=1)
681 call diag_update_remap_grids(cs%diag)
682
683 ! Offer diagnostic fields for averaging.
684 if (query_averaging_enabled(cs%diag)) then
685 if (cs%id_urestrat_time > 0) call post_data(cs%id_urestrat_time, utimescale_diag, cs%diag)
686 if (cs%id_vrestrat_time > 0) call post_data(cs%id_vrestrat_time, vtimescale_diag, cs%diag)
687 if (cs%id_uhml > 0) call post_data(cs%id_uhml, uhml, cs%diag)
688 if (cs%id_vhml > 0) call post_data(cs%id_vhml, vhml, cs%diag)
689 if (cs%id_BLD > 0) call post_data(cs%id_BLD, mld_fast, cs%diag)
690 if (cs%id_MLD > 0) call post_data(cs%id_MLD, mld_slow, cs%diag)
691 if (cs%id_Rml > 0) call post_data(cs%id_Rml, rml_av_fast, cs%diag)
692 if (cs%id_uDml > 0) call post_data(cs%id_uDml, udml_diag, cs%diag)
693 if (cs%id_vDml > 0) call post_data(cs%id_vDml, vdml_diag, cs%diag)
694 if (cs%id_mle_fl > 0) call post_data(cs%id_mle_fl, mle_fl_2d, cs%diag)
695
696 if (cs%id_uml > 0) then
697 do j=js,je ; do i=is-1,ie
698 h_vel = 0.5*((htot_fast(i,j) + htot_fast(i+1,j)) + h_neglect)
699 udml_diag(i,j) = udml_diag(i,j) / (0.01*h_vel) * g%IdyCu(i,j) * (mu(0.,0.)-mu(-.01,0.))
700 enddo ; enddo
701 call post_data(cs%id_uml, udml_diag, cs%diag)
702 endif
703 if (cs%id_vml > 0) then
704 do j=js-1,je ; do i=is,ie
705 h_vel = 0.5*((htot_fast(i,j) + htot_fast(i,j+1)) + h_neglect)
706 vdml_diag(i,j) = vdml_diag(i,j) / (0.01*h_vel) * g%IdxCv(i,j) * (mu(0.,0.)-mu(-.01,0.))
707 enddo ; enddo
708 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 call diag_update_remap_grids(cs%diag)
715
716end subroutine mixedlayer_restrat_om4
717
718!> Stream function shape as a function of non-dimensional position within mixed-layer [nondim]
719real 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 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 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 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 bottop = 0.5*(1. - sign(1., sigma + 0.5)) ! =0 for sigma>-0.5, =1 for sigma<-0.5
751
752 mu = max(mu, dd*bottop) ! Combines original psi1 with tail
753end function mu
754
755!> Calculates a restratifying flow in the mixed layer, following the formulation
756!! used in Bodner et al., 2023 (B22)
757subroutine 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 real :: uhml(SZIB_(G),SZJ_(G),SZK_(GV)) ! zonal mixed layer transport [H L2 T-1 ~> m3 s-1 or kg s-1]
780 real :: vhml(SZI_(G),SZJB_(G),SZK_(GV)) ! merid mixed layer transport [H L2 T-1 ~> m3 s-1 or kg s-1]
781 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 little_h, & ! "Little h" representing active mixing layer depth [H ~> m or kg m-2]
785 big_H, & ! "Big H" representing the mixed layer depth [H ~> m or kg m-2]
786 mld, & ! The mixed layer depth returned by detect_mld [H ~> m or kg m-2]
787 htot, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2]
788 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 wpup ! Turbulent vertical momentum [L H T-2 ~> m2 s-2 or kg m-1 s-2]
791 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 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 real :: lf_bodner_diag(SZI_(G),SZJ_(G)) ! Front width as in Bodner et al., 2023 (B22), eq 24 [L ~> m]
794 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 real :: covTS(SZI_(G)) ! SGS TS covariance in Stanley param; currently 0 [C S ~> degC ppt]
798 real :: varS(SZI_(G)) ! SGS S variance in Stanley param; currently 0 [S2 ~> ppt2]
799 real :: dmu(SZK_(GV)) ! Change in mu(z) across layer k [nondim]
800 real :: Rml_int(SZI_(G)) ! Potential density integrated through the mixed layer [R H ~> kg m-2 or kg2 m-5]
801 real :: SpV_ml(SZI_(G)) ! Specific volume evaluated at the surface pressure [R-1 ~> m3 kg-1]
802 real :: SpV_int(SZI_(G)) ! Specific volume integrated through the mixed layer [H R-1 ~> m4 kg-1 or m]
803 real :: rho_ml(SZI_(G)) ! Potential density relative to the surface [R ~> kg m-3]
804 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 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
838 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
839
840 i4dt = 0.25 / dt
841 g_rho0 = gv%H_to_Z * gv%g_Earth / gv%Rho0
842 h_neglect = gv%H_subroundoff
843
844 covts(:) = 0.0 ! Might be in tv% in the future. Not implemented for the time being.
845 vars(:) = 0.0 ! Ditto.
846
847 ! This value is roughly (pi / (the age of the universe) )^2.
848 absurdly_small_freq2 = 1e-34*us%T_to_s**2
849
850 if (.not.associated(tv%eqn_of_state)) call mom_error(fatal, "mixedlayer_restrat_Bodner: "// &
851 "An equation of state must be used with this module.")
852 if (cs%MLE_use_PBL_MLD) then
853 if (cs%MLE_density_diff > 0.) call mom_error(fatal, "mixedlayer_restrat_Bodner: "// &
854 "MLE_density_diff is +ve and should not be in mixedlayer_restrat_Bodner.")
855 if (.not.associated(bflux)) call mom_error(fatal, "mixedlayer_restrat_Bodner: "// &
856 "Surface buoyancy flux was not associated.")
857 else
858 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 "Bodner_detect_MLD must be True.")
861 endif
862
863 if (associated(bflux)) &
864 call pass_var(bflux, g%domain, halo=1)
865
866 ! Extract the friction velocity from the forcing type.
867 call find_ustar(forces, tv, u_star_2d, g, gv, us, halo=1)
868
869 if (cs%debug) then
870 call hchksum(h,'mixed_Bodner: h', g%HI, haloshift=1, unscale=gv%H_to_mks)
871 call hchksum(bld, 'mle_Bodner: BLD', g%HI, haloshift=1, unscale=us%Z_to_m)
872 call hchksum(h_mld, 'mle_Bodner: h_MLD', g%HI, haloshift=1, unscale=gv%H_to_mks)
873 if (associated(bflux)) &
874 call hchksum(bflux, 'mle_Bodner: bflux', g%HI, haloshift=1, unscale=us%Z_to_m**2*us%s_to_T**3)
875 call hchksum(u_star_2d, 'mle_Bodner: u*', g%HI, haloshift=1, unscale=us%Z_to_m*us%s_to_T)
876 call hchksum(cs%MLD_filtered, 'mle_Bodner: MLD_filtered 1', &
877 g%HI, haloshift=1, unscale=gv%H_to_mks)
878 call hchksum(cs%MLD_filtered_slow,'mle_Bodner: MLD_filtered_slow 1', &
879 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 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 cs%BLD_growing_Tfilt, cs%BLD_decaying_Tfilt, dt)
887 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 if (cs%MLD_grid) then
892 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 cs%MLD_growing_Tfilt, cs%MLD_Tfilt_space(i,j), dt)
895 enddo ; enddo
896 elseif (cs%Bodner_detect_MLD) then
897 call detect_mld(h, tv, mld, g, gv, cs)
898 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 cs%MLD_growing_Tfilt, cs%MLD_decaying_Tfilt, dt)
901 enddo ; enddo
902 else
903 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 cs%MLD_growing_Tfilt, cs%MLD_decaying_Tfilt, dt)
906 enddo ; enddo
907 endif
908 do j=js-1,je+1 ; do i=is-1,ie+1
909 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 if (allocated(tv%SpV_avg) .and. .not.(gv%Boussinesq .or. gv%semi_Boussinesq)) then
914 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 * (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 elseif (cs%answer_date < 20240201) then
930 z3_t3_to_m3_s3 = (us%Z_to_m * us%s_to_T)**3
931 m2_s2_to_z2_t2 = (us%m_to_Z * us%T_to_s)**2
932 do j=js-1,je+1 ; do i=is-1,ie+1
933 w_star3 = max(0., -bflux(i,j)) * bld(i,j) ! In [Z3 T-3 ~> m3 s-3]
934 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 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 do j=js-1,je+1 ; do i=is-1,ie+1
940 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 * 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 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 cs%BLD_growing_Tfilt, cs%BLD_decaying_Tfilt, dt)
950 cs%wpup_filtered(i,j) = wpup(i,j)
951 enddo ; enddo
952
953 if (cs%id_lfbod > 0) then
954 do j=js-1,je+1 ; do i=is-1,ie+1
955 ! Calculate front length used in B22 formula (eq 24).
956 w_star3 = max(0., -bflux(i,j)) * bld(i,j)
957 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 + (g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j-1)))
962 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 / (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 if (allocated(tv%SpV_avg) .and. .not.(gv%Boussinesq .or. gv%semi_Boussinesq)) then
971 do j=js-1,je+1 ; do i=is-1,ie+1
972 lf_bodner_diag(i,j) = lf_bodner_diag(i,j) &
973 * (us%Z_to_L * gv%RZ_to_H / tv%SpV_avg(i,j,1))
974 enddo ; enddo
975 else
976 do j=js-1,je+1 ; do i=is-1,ie+1
977 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 if (cs%debug) then
983 call hchksum(little_h,'mle_Bodner: little_h', g%HI, haloshift=1, unscale=gv%H_to_mks)
984 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 g%HI, haloshift=1, unscale=gv%H_to_mks)
987 call hchksum(cs%MLD_filtered_slow,'mle_Bodner: MLD_filtered_slow 2', &
988 g%HI, haloshift=1, unscale=gv%H_to_mks)
989 call hchksum(wpup,'mle_Bodner: wpup', g%HI, haloshift=1, unscale=us%L_to_m*gv%H_to_mks*us%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 p0(:) = 0.0
996 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 do j=js-1,je+1
1005 rho_ml(:) = 0.0 ; spv_ml(:) = 0.0
1006 do i=is-1,ie+1
1007 htot(i,j) = 0.0 ; rml_int(i) = 0.0 ; spv_int(i) = 0.0
1008 enddo
1009 keep_going = .true.
1010 do k=1,nz
1011 do i=is-1,ie+1
1012 vol_dt_avail(i,j,k) = max(i4dt*g%areaT(i,j)*(h(i,j,k)-gv%Angstrom_H),0.0)
1013 enddo
1014 if (keep_going) then
1015 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
1016 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 rho_ml, tv%eqn_of_state, eosdom)
1019 else
1020 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, rho_ml, tv%eqn_of_state, eosdom)
1021 endif
1022 else
1023 call calculate_spec_vol(tv%T(:,j,k), tv%S(:,j,k), p0, spv_ml, tv%eqn_of_state, eosdom)
1024 endif
1025 line_is_empty = .true.
1026 do i=is-1,ie+1
1027 if (htot(i,j) < big_h(i,j)) then
1028 dh = min( h(i,j,k), big_h(i,j) - htot(i,j) )
1029 rml_int(i) = rml_int(i) + dh*rho_ml(i) ! Rml_int has units of [R H ~> kg m-2]
1030 spv_int(i) = spv_int(i) + dh*spv_ml(i) ! SpV_int has units of [H R-1 ~> m4 kg-1 or m]
1031 htot(i,j) = htot(i,j) + dh
1032 line_is_empty = .false.
1033 endif
1034 enddo
1035 if (line_is_empty) keep_going=.false.
1036 endif
1037 enddo
1038
1039 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
1040 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 buoy_av(i,j) = -( g_rho0 * rml_int(i) ) / (htot(i,j) + h_neglect)
1043 enddo
1044 else
1045 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 buoy_av(i,j) = (gv%H_to_RZ*gv%g_Earth * spv_int(i)) / (htot(i,j) + h_neglect)
1048 enddo
1049 endif
1050 enddo
1051
1052 if (cs%debug) then
1053 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 unscale=us%L_to_m**2*gv%H_to_mks*us%s_to_T)
1056 call hchksum(buoy_av,'mle_Bodner: buoy_av', g%HI, haloshift=1, unscale=gv%m_to_H*us%L_T_to_m_s**2)
1057 endif
1058
1059 ! U - Component
1060 !$OMP do
1061 do j=js,je ; do i=is-1,ie
1062 if (g%OBCmaskCu(i,j) > 0.) then
1063 grid_dsd = sqrt(0.5*( g%dxCu(i,j)**2 + g%dyCu(i,j)**2 )) * g%dyCu(i,j) ! [L2 ~> m2]
1064 absf = 0.5*(abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i,j))) ! [T-1 ~> s-1]
1065 h_sml = 0.5*( little_h(i,j) + little_h(i+1,j) ) ! [H ~> m or kg m-2]
1066 h_big = 0.5*( big_h(i,j) + big_h(i+1,j) ) ! [H ~> m or kg m-2]
1067 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 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 * ( 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 psi_mag = 0.0
1073 endif
1074
1075 ihtot = 2.0 / ((htot(i,j) + htot(i+1,j)) + h_neglect) ! [H-1 ~> m-1 or m2 kg-1]
1076 sigint = 0.0
1077 muzb = 0.0 ! This will be the first value of muza = mu(z=0)
1078 do k=1,nz
1079 muza = muzb ! mu(z/MLD) for upper interface [nondim]
1080 hatvel = 0.5*(h(i,j,k) + h(i+1,j,k)) ! Thickness at velocity point [H ~> m or kg m-2]
1081 sigint = sigint - (hatvel * ihtot) ! z/H for lower interface [nondim]
1082 muzb = mu(sigint, cs%MLE_tail_dh) ! mu(z/MLD) for lower interface [nondim]
1083 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 if (dmu(k)*psi_mag > 0.0) then
1087 if (dmu(k)*psi_mag > vol_dt_avail(i,j,k)) psi_mag = vol_dt_avail(i,j,k) / dmu(k)
1088 elseif (dmu(k)*psi_mag < 0.0) then
1089 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 do k=1,nz
1093 uhml(i,j,k) = dmu(k) * psi_mag ! [L2 H T-1 ~> m3 s-1 or kg s-1]
1094 uhtr(i,j,k) = uhtr(i,j,k) + uhml(i,j,k) * dt ! [L2 H ~> m3 or kg]
1095 enddo
1096
1097 udml_diag(i,j) = psi_mag
1098 enddo ; enddo
1099
1100 ! V- component
1101 !$OMP do
1102 do j=js-1,je ; do i=is,ie
1103 if (g%OBCmaskCv(i,j) > 0.) then
1104 grid_dsd = sqrt(0.5*( g%dxCv(i,j)**2 + g%dyCv(i,j)**2 )) * g%dxCv(i,j) ! [L2 ~> m2]
1105 absf = 0.5*(abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j))) ! [T-1 ~> s-1]
1106 h_sml = 0.5*( little_h(i,j) + little_h(i,j+1) ) ! [H ~> m or kg m-2]
1107 h_big = 0.5*( big_h(i,j) + big_h(i,j+1) ) ! [H ~> m or kg m-2]
1108 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 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 * ( 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 psi_mag = 0.0
1114 endif
1115
1116 ihtot = 2.0 / ((htot(i,j) + htot(i,j+1)) + h_neglect) ! [H-1 ~> m-1 or m2 kg-1]
1117 sigint = 0.0
1118 muzb = 0.0 ! This will be the first value of muza = mu(z=0)
1119 do k=1,nz
1120 muza = muzb ! mu(z/MLD) for upper interface [nondim]
1121 hatvel = 0.5*(h(i,j,k) + h(i,j+1,k)) ! Thickness at velocity point [H ~> m or kg m-2]
1122 sigint = sigint - (hatvel * ihtot) ! z/H for lower interface [nondim]
1123 muzb = mu(sigint, cs%MLE_tail_dh) ! mu(z/MLD) for lower interface [nondim]
1124 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 if (dmu(k)*psi_mag > 0.0) then
1128 if (dmu(k)*psi_mag > vol_dt_avail(i,j,k)) psi_mag = vol_dt_avail(i,j,k) / dmu(k)
1129 elseif (dmu(k)*psi_mag < 0.0) then
1130 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 do k=1,nz
1134 vhml(i,j,k) = dmu(k) * psi_mag ! [L2 H T-1 ~> m3 s-1 or kg s-1]
1135 vhtr(i,j,k) = vhtr(i,j,k) + vhml(i,j,k) * dt ! [L2 H ~> m3 or kg]
1136 enddo
1137
1138 vdml_diag(i,j) = psi_mag
1139 enddo ; enddo
1140
1141 !$OMP do
1142 do j=js,je ; do k=1,nz ; do i=is,ie
1143 h(i,j,k) = h(i,j,k) - dt*g%IareaT(i,j) * &
1144 ((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 if (cs%id_uhml > 0 .or. cs%id_vhml > 0) &
1149 ! Remapped uhml and vhml require east/north halo updates of h
1150 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 call diag_update_remap_grids(cs%diag)
1154
1155 ! Offer diagnostic fields for averaging.
1156 if (query_averaging_enabled(cs%diag)) then
1157 if (cs%id_ustar > 0) call post_data(cs%id_ustar, u_star_2d, cs%diag)
1158 if (cs%id_bflux > 0) call post_data(cs%id_bflux, bflux, cs%diag)
1159 if (cs%id_wpup > 0) call post_data(cs%id_wpup, wpup, cs%diag)
1160 if (cs%id_Rml > 0) call post_data(cs%id_Rml, buoy_av, cs%diag)
1161 if (cs%id_BLD > 0) call post_data(cs%id_BLD, little_h, cs%diag)
1162 if (cs%id_MLD > 0) call post_data(cs%id_MLD, big_h, cs%diag)
1163 if (cs%id_uhml > 0) call post_data(cs%id_uhml, uhml, cs%diag)
1164 if (cs%id_vhml > 0) call post_data(cs%id_vhml, vhml, cs%diag)
1165 if (cs%id_uDml > 0) call post_data(cs%id_uDml, udml_diag, cs%diag)
1166 if (cs%id_vDml > 0) call post_data(cs%id_vDml, vdml_diag, cs%diag)
1167 if (cs%id_lfbod > 0) call post_data(cs%id_lfbod, lf_bodner_diag, cs%diag)
1168
1169 if (cs%id_uml > 0) then
1170 do j=js,je ; do i=is-1,ie
1171 h_vel = 0.5*((htot(i,j) + htot(i+1,j)) + h_neglect)
1172 udml_diag(i,j) = udml_diag(i,j) / (0.01*h_vel) * g%IdyCu(i,j) * (mu(0.,0.)-mu(-.01,0.))
1173 enddo ; enddo
1174 call post_data(cs%id_uml, udml_diag, cs%diag)
1175 endif
1176 if (cs%id_vml > 0) then
1177 do j=js-1,je ; do i=is,ie
1178 h_vel = 0.5*((htot(i,j) + htot(i,j+1)) + h_neglect)
1179 vdml_diag(i,j) = vdml_diag(i,j) / (0.01*h_vel) * g%IdxCv(i,j) * (mu(0.,0.)-mu(-.01,0.))
1180 enddo ; enddo
1181 call post_data(cs%id_vml, vdml_diag, cs%diag)
1182 endif
1183 endif
1184
1185end 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.
1199real 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 if (signal>=filtered) then
1211 rt = 1.0 / ( dt + tau_growing )
1212 afac = tau_growing * rt
1213 bfac = 1. - afac
1214 else
1215 rt = 1.0 / ( dt + tau_decaying )
1216 afac = tau_decaying * rt
1217 bfac = 1. - afac
1218 endif
1219
1220 rmean2ts = afac * filtered + bfac * signal
1221
1222end function rmean2ts
1223
1224!> Calculates a restratifying flow assuming a 2-layer bulk mixed layer.
1225subroutine 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 real :: uhml(SZIB_(G),SZJ_(G),SZK_(GV)) ! Restratifying zonal thickness transports [H L2 T-1 ~> m3 s-1 or kg s-1]
1241 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 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 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 htot, & ! The sum of the thicknesses of layers in the mixed layer [H ~> m or kg m-2]
1250 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 real :: Rho_ml(SZI_(G)) ! Potential density relative to the surface [R ~> kg m-3]
1255 real :: rho_int(SZI_(G)) ! The integral of density over the mixed layer depth [R H ~> kg m-2 or kg2 m-5]
1256 real :: SpV_ml(SZI_(G)) ! Specific volume evaluated at the surface pressure [R-1 ~> m3 kg-1]
1257 real :: SpV_int(SZI_(G)) ! Specific volume integrated through the surface layer [H R-1 ~> m4 kg-1 or m]
1258 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 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 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 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 real :: utimescale_diag(SZIB_(G),SZJ_(G)) ! Zonal restratification timescale [T ~> s], stored for diagnostics.
1280 real :: vtimescale_diag(SZI_(G),SZJB_(G)) ! Meridional restratification timescale [T ~> s], stored for diagnostics.
1281 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 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 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1288 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nkml = gv%nkml
1289
1290 if (.not. cs%initialized) call mom_error(fatal, "mixedlayer_restrat_BML: "// &
1291 "Module must be initialized before it is used.")
1292
1293 if ((nkml<2) .or. (cs%ml_restrat_coef<=0.0)) return
1294
1295
1296 h_min = 0.5*gv%Angstrom_H ! This should be GV%Angstrom_H, but that value would change answers.
1297 udml(:) = 0.0 ; vdml(:) = 0.0
1298 i4dt = 0.25 / dt
1299 g_rho0 = gv%H_to_Z * gv%g_Earth / gv%Rho0
1300 vonkar_x_pi2 = cs%vonKar * 9.8696
1301 use_eos = associated(tv%eqn_of_state)
1302 h_neglect = gv%H_subroundoff
1303
1304 if (.not.use_eos) call mom_error(fatal, "mixedlayer_restrat_BML: "// &
1305 "An equation of state must be used with this module.")
1306
1307 if (cs%use_Stanley_ML) call mom_error(fatal, "mixedlayer_restrat_BML: "// &
1308 "The Stanley parameterization is not available with the BML.")
1309
1310 ! Extract the friction velocity from the forcing type.
1311 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 p0(:) = 0.0
1316 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 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
1322 !$OMP do
1323 do j=js-1,je+1
1324 do i=is-1,ie+1
1325 htot(i,j) = 0.0 ; rho_int(i) = 0.0
1326 enddo
1327 do k=1,nkml
1328 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p0, rho_ml(:), tv%eqn_of_state, eosdom)
1329 do i=is-1,ie+1
1330 rho_int(i) = rho_int(i) + h(i,j,k)*rho_ml(i)
1331 htot(i,j) = htot(i,j) + h(i,j,k)
1332 h_avail(i,j,k) = max(i4dt*g%areaT(i,j)*(h(i,j,k)-gv%Angstrom_H),0.0)
1333 enddo
1334 enddo
1335
1336 do i=is-1,ie+1
1337 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 do j=js-1,je+1
1343 do i=is-1,ie+1
1344 htot(i,j) = 0.0 ; spv_int(i) = 0.0
1345 enddo
1346 do k=1,nkml
1347 call calculate_spec_vol(tv%T(:,j,k), tv%S(:,j,k), p0, spv_ml, tv%eqn_of_state, eosdom)
1348 do i=is-1,ie+1
1349 spv_int(i) = spv_int(i) + h(i,j,k)*spv_ml(i) ! [H R-1 ~> m4 kg-1 or m]
1350 htot(i,j) = htot(i,j) + h(i,j,k)
1351 h_avail(i,j,k) = max(i4dt*g%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 do i=is-1,ie+1
1357 rml_av(i,j) = (-gv%H_to_RZ*gv%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 do j=js,je ; do i=is-1,ie
1369 h_vel = 0.5*(htot(i,j) + htot(i+1,j))
1370
1371 u_star = max(cs%ustar_min, 0.5*(u_star_2d(i,j) + u_star_2d(i+1,j)))
1372
1373 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_ml*u_star)/(absf*h_ml + 4.0*u_star)
1379 ! momentum mixing rate: pi^2*visc/h_ml^2
1380 mom_mixrate = vonkar_x_pi2*u_star**2 / &
1381 (absf*h_vel**2 + 4.0*(h_vel+h_neglect)*u_star)
1382 timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
1383
1384 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 (rml_av(i+1,j)-rml_av(i,j)) * (h_vel**2)
1389
1390 if (udml(i) == 0) then
1391 do k=1,nkml ; uhml(i,j,k) = 0.0 ; enddo
1392 else
1393 i2htot = 1.0 / (htot(i,j) + htot(i+1,j) + h_neglect)
1394 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 do k=1,nkml
1398 hx2 = (h(i,j,k) + h(i+1,j,k) + h_neglect)
1399 a(k) = (hx2 * i2htot) * (2.0 - 4.0*(z_topx2+0.5*hx2)*i2htot)
1400 z_topx2 = z_topx2 + hx2
1401 if (a(k)*udml(i) > 0.0) then
1402 if (a(k)*udml(i) > h_avail(i,j,k)) udml(i) = h_avail(i,j,k) / a(k)
1403 else
1404 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 do k=1,nkml
1408 uhml(i,j,k) = a(k)*udml(i)
1409 uhtr(i,j,k) = uhtr(i,j,k) + uhml(i,j,k)*dt
1410 enddo
1411 endif
1412
1413 udml_diag(i,j) = udml(i)
1414 utimescale_diag(i,j) = timescale
1415 enddo ; enddo
1416
1417! V- component
1418 !$OMP do
1419 do j=js-1,je ; do i=is,ie
1420 h_vel = 0.5*(htot(i,j) + htot(i,j+1))
1421
1422 u_star = max(cs%ustar_min, 0.5*(u_star_2d(i,j) + u_star_2d(i,j+1)))
1423
1424 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_ml*u_star)/(absf*h_ml + 4.0*u_star)
1430 ! momentum mixing rate: pi^2*visc/h_ml^2
1431 mom_mixrate = vonkar_x_pi2*u_star**2 / &
1432 (absf*h_vel**2 + 4.0*(h_vel+h_neglect)*u_star)
1433 timescale = 0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + mom_mixrate**2)
1434
1435 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 (rml_av(i,j+1)-rml_av(i,j)) * (h_vel**2)
1440 if (vdml(i) == 0) then
1441 do k=1,nkml ; vhml(i,j,k) = 0.0 ; enddo
1442 else
1443 i2htot = 1.0 / (htot(i,j) + htot(i,j+1) + h_neglect)
1444 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 do k=1,nkml
1448 hx2 = (h(i,j,k) + h(i,j+1,k) + h_neglect)
1449 a(k) = (hx2 * i2htot) * (2.0 - 4.0*(z_topx2+0.5*hx2)*i2htot)
1450 z_topx2 = z_topx2 + hx2
1451 if (a(k)*vdml(i) > 0.0) then
1452 if (a(k)*vdml(i) > h_avail(i,j,k)) vdml(i) = h_avail(i,j,k) / a(k)
1453 else
1454 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 do k=1,nkml
1458 vhml(i,j,k) = a(k)*vdml(i)
1459 vhtr(i,j,k) = vhtr(i,j,k) + vhml(i,j,k)*dt
1460 enddo
1461 endif
1462
1463 vtimescale_diag(i,j) = timescale
1464 vdml_diag(i,j) = vdml(i)
1465 enddo ; enddo
1466
1467 !$OMP do
1468 do j=js,je ; do k=1,nkml ; do i=is,ie
1469 h(i,j,k) = h(i,j,k) - dt*g%IareaT(i,j) * &
1470 ((uhml(i,j,k) - uhml(i-1,j,k)) + (vhml(i,j,k) - vhml(i,j-1,k)))
1471 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 if (cs%id_uhml > 0 .or. cs%id_vhml > 0) &
1478 ! Remapped uhml and vhml require east/north halo updates of h
1479 call pass_var(h, g%domain, to_west+to_south+omit_corners, halo=1)
1480 call diag_update_remap_grids(cs%diag)
1481
1482 ! Offer diagnostic fields for averaging.
1483 if (query_averaging_enabled(cs%diag) .and. &
1484 ((cs%id_urestrat_time > 0) .or. (cs%id_vrestrat_time > 0))) then
1485 call post_data(cs%id_urestrat_time, utimescale_diag, cs%diag)
1486 call post_data(cs%id_vrestrat_time, vtimescale_diag, cs%diag)
1487 endif
1488 if (query_averaging_enabled(cs%diag) .and. &
1489 ((cs%id_uhml>0) .or. (cs%id_vhml>0))) then
1490 do k=nkml+1,nz
1491 do j=js,je ; do i=isq,ieq ; uhml(i,j,k) = 0.0 ; enddo ; enddo
1492 do j=jsq,jeq ; do i=is,ie ; vhml(i,j,k) = 0.0 ; enddo ; enddo
1493 enddo
1494 if (cs%id_uhml > 0) call post_data(cs%id_uhml, uhml, cs%diag)
1495 if (cs%id_vhml > 0) call post_data(cs%id_vhml, vhml, cs%diag)
1496 if (cs%id_MLD > 0) call post_data(cs%id_MLD, htot, cs%diag)
1497 if (cs%id_Rml > 0) call post_data(cs%id_Rml, rml_av, cs%diag)
1498 if (cs%id_uDml > 0) call post_data(cs%id_uDml, udml_diag, cs%diag)
1499 if (cs%id_vDml > 0) call post_data(cs%id_vDml, vdml_diag, cs%diag)
1500 endif
1501
1502end subroutine mixedlayer_restrat_bml
1503
1504!> Detects the mixed layer depth using a density difference criterion (MLE_DENSITY_DIFF)
1505subroutine 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 real, dimension(SZI_(G)) :: pRef_MLD ! A reference pressure for calculating the mixed layer
1515 ! densities [R L2 T-2 ~> Pa].
1516 real, dimension(SZI_(G)) :: rhoSurf, deltaRhoAtKm1, deltaRhoAtK ! Densities and density differences [R ~> kg m-3]
1517 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 real :: covTS(SZI_(G)) ! SGS TS covariance in Stanley param; currently 0 [C S ~> degC ppt]
1521 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 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1526
1527 covts(:) = 0.0 ! Might be in tv% in the future. Not implemented for the time being.
1528 vars(:) = 0.0 ! Ditto.
1529
1530 !! TODO: use derivatives and mid-MLD pressure. Currently this is sigma-0. -AJA
1531 pref_mld(:) = 0.
1532 eosdom(:) = eos_domain(g%HI, halo=1)
1533 do j=js-1,je+1
1534 dk(:) = 0.5 * h(:,j,1) ! Depth of center of surface layer
1535 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 rhosurf, tv%eqn_of_state, eosdom)
1538 else
1539 call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pref_mld, rhosurf, tv%eqn_of_state, eosdom)
1540 endif
1541 deltarhoatk(:) = 0.
1542 mld_fast(:,j) = 0.
1543 do k=2,nz
1544 dkm1(:) = dk(:) ! Depth of center of layer K-1
1545 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 deltarhoatkm1(:) = deltarhoatk(:) ! Store value from previous iteration of K
1548 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 deltarhoatk, tv%eqn_of_state, eosdom)
1551 else
1552 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pref_mld, deltarhoatk, tv%eqn_of_state, eosdom)
1553 endif
1554 do i=is-1,ie+1
1555 deltarhoatk(i) = deltarhoatk(i) - rhosurf(i) ! Density difference between layer K and surface
1556 enddo
1557 do i=is-1,ie+1
1558 ddrho = deltarhoatk(i) - deltarhoatkm1(i)
1559 if ((mld_fast(i,j)==0.) .and. (ddrho>0.) .and. &
1560 (deltarhoatkm1(i)<cs%MLE_density_diff) .and. (deltarhoatk(i)>=cs%MLE_density_diff)) then
1561 afac = ( cs%MLE_density_diff - deltarhoatkm1(i) ) / ddrho
1562 mld_fast(i,j) = dk(i) * afac + dkm1(i) * (1. - afac)
1563 endif
1564 enddo ! i-loop
1565 enddo ! k-loop
1566 do i=is-1,ie+1
1567 mld_fast(i,j) = cs%MLE_MLD_stretch * mld_fast(i,j)
1568 if ((mld_fast(i,j)==0.) .and. (deltarhoatk(i)<cs%MLE_density_diff)) &
1569 mld_fast(i,j) = dk(i) ! Assume mixing to the bottom
1570 enddo
1571 enddo ! j-loop
1572end 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]
1578real 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_ml*u_star)/(absf*h_ml + 4.0*u_star) + Kv_water
1597 ! momentum mixing rate: pi^2*visc/h_ml^2
1598 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 if (kv_rest <= 0.0) then
1601 ! This case reproduces the previous answers, but the extra h_neg is otherwise unnecessary.
1602 mom_mixrate = (pi2*vonkar)*u_star**2 / (absf*hbl**2 + 4.0*(hbl + h_neg)*u_star)
1603 growth_time = restrat_coef * (0.0625 * (absf + 2.0*mom_mixrate) / (absf**2 + 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 = pi2*vonKar*u_star**2 / (absf*hBL**2 + 4.0*hBL*u_star) + pi2*Kv_rest / hBL**2
1607 if (absf*hbl <= 4.0e-16*u_star) then
1608 kv_eff = pi2 * (kv_rest + 0.25*vonkar*hbl*u_star)
1609 else
1610 kv_eff = pi2 * (kv_rest + vonkar*u_star**2*hbl / (absf*hbl + 4.0*u_star))
1611 endif
1612 growth_time = (restrat_coef*0.0625) * ((hbl**2*(hbl**2*absf + 2.0*kv_eff)) / ((hbl**2*absf)**2 + kv_eff**2))
1613 endif
1614
1615end function growth_time
1616
1617!> Initialize the mixed layer restratification module
1618logical 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 default=.false., do_not_log=.true.)
1648 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 "BULKMIXEDLAYER is true.", default=.false.)
1654 if (.not. mixedlayer_restrat_init) return
1655
1656 cs%initialized = .true.
1657 cs%Time => time
1658
1659 ! Nonsense values to cause problems when these parameters are not used
1660 cs%MLE_MLD_decay_time = -9.e9*us%s_to_T
1661 cs%MLE_density_diff = -9.e9*us%kg_m3_to_R
1662 cs%MLE_tail_dh = -9.e9
1663 cs%MLE_use_PBL_MLD = .false.
1664 cs%MLE_MLD_stretch = -9.e9
1665 cs%use_Stanley_ML = .false.
1666 cs%use_Bodner = .false.
1667 cs%fl_from_file = .false.
1668 cs%MLD_grid = .false.
1669 cs%Cr_grid = .false.
1670
1671 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 default=99991231, do_not_log=.true.)
1675 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1676 call openparameterblock(param_file,'MLE') ! Prepend MLE% to all parameters
1677 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 default=.false.)
1682 endif
1683 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 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 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 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 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 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 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 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 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 units="m2 s-2", default=1.0e-24, scale=us%m_to_Z**2*us%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 "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 "in ML restrat code.", default=.false.)
1738 call get_param(param_file, mdl, "USE_CR_GRID", cs%Cr_grid, &
1739 "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 "If true, read in a spatially varying MLD_decaying_Tfilt field.", default=.false.)
1742 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 default="")
1746 call get_param(param_file, mdl, "MLD_TFILT_VAR", varname, &
1747 "The variable name for MLD_decaying_Tfilt field.", &
1748 default="MLD_tfilt")
1749 filename = trim(inputdir) // "/" // trim(filename)
1750 allocate(cs%MLD_Tfilt_space(g%isd:g%ied,g%jsd:g%jed), source=0.0)
1751 call mom_read_data(filename, varname, cs%MLD_Tfilt_space, g%domain, scale=us%s_to_T)
1752 call pass_var(cs%MLD_Tfilt_space, g%domain)
1753 endif
1754 allocate(cs%Cr_space(g%isd:g%ied,g%jsd:g%jed), source=cs%Cr)
1755 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 default="")
1759 call get_param(param_file, mdl, "CR_VAR", varname, &
1760 "The variable name for Cr field.", &
1761 default="Cr")
1762 filename = trim(inputdir) // "/" // trim(filename)
1763 call mom_read_data(filename, varname, cs%Cr_space, g%domain)
1764 call pass_var(cs%Cr_space, g%domain)
1765 endif
1766 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 "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 "detected via the density difference criterion MLE_DENSITY_DIFF.", default=.false.)
1776 if (.not.(cs%MLE_use_PBL_MLD.or.cs%Bodner_detect_MLD)) call mom_error(fatal, "mixedlayer_restrat_init: "// &
1777 "To use MLE%USE_BODNER23=True then MLE_USE_PBL_MLD or BODNER_DETECT_MLD must be true.")
1778 if (cs%MLE_use_PBL_MLD.and.cs%Bodner_detect_MLD) call mom_error(fatal, "mixedlayer_restrat_init: "// &
1779 "MLE_USE_PBL_MLD and BODNER_DETECT_MLD cannot both be true.")
1780 else
1781 call closeparameterblock(param_file) ! The remaining parameters do not have MLE% prepended
1782 endif
1783
1784 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 "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 "in ML restrat code.", default=.false.)
1799 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 units="nondim", default=-1.0, do_not_log=.true.)
1803 if (stanley_coeff < 0.0) call mom_error(fatal, &
1804 "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 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 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 "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 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 "If true, the MLE front-length scale is read from a file.", default=.false.)
1823 if (cs%fl_from_file) then
1824 call time_interp_external_init()
1825
1826 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 "It is used when MLE_FRONT_LENGTH_FROM_FILE is true.", fail_if_missing=.true.)
1830 mle_fl_filename = trim(slasher(inputdir))//trim(mle_fl_file)
1831 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 "Name of MLE front-length scale variable in MLE_FL_FILE.", default='mle_fl')
1834 if (modulo(g%Domain%turns, 4) /= 0) then
1835 cs%sbc_fl = init_external_field(mle_fl_filename, trim(fl_varname), mom_domain=g%Domain%domain_in)
1836 else
1837 cs%sbc_fl = init_external_field(mle_fl_filename, trim(fl_varname), mom_domain=g%Domain)
1838 endif
1839 endif
1840 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 "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 "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 "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 "is instantaneously set to the current MLD.", units="s", default=0., scale=us%s_to_T)
1859 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 "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 "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 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 units="m2 s-1", default=0.0, scale=gv%m2_s_to_HZ_T*(us%Z_to_m*gv%m_to_H))
1879 call get_param(param_file, mdl, "OMEGA", omega, &
1880 "The rotation rate of the earth.", &
1881 units="s-1", default=7.2921e-5, scale=us%T_to_s)
1882 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 units="m s-1", default=us%Z_to_m*us%s_to_T*ustar_min_dflt, scale=gv%m_to_H*us%T_to_s)
1888 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 "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 units="nondim", default=1.0)
1897 endif
1898
1899 cs%diag => diag
1900
1901 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 '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 '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 '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 '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 '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 '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 '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 'm3 s-1', conversion=gv%H_to_m*(us%L_to_m**2)*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 'm3 s-1', conversion=gv%H_to_m*(us%L_to_m**2)*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 '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 'm s-1', conversion=us%L_T_to_m_s)
1934 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 'm2 s-2', conversion=us%L_to_m*gv%H_to_m*us%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 '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 'm2 s-3', conversion=(us%Z_to_m**2*us%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 '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 '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 if (allocated(cs%MLD_filtered)) call pass_var(cs%MLD_filtered, g%domain)
1955 if (allocated(cs%MLD_filtered_slow)) call pass_var(cs%MLD_filtered_slow, g%domain)
1956 if (allocated(cs%wpup_filtered)) call pass_var(cs%wpup_filtered, g%domain)
1957
1958end function mixedlayer_restrat_init
1959
1960!> Allocate and register fields in the mixed layer restratification structure for restarts
1961subroutine 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 default=.false., do_not_log=.true.)
1977 if (.not. mixedlayer_restrat_init) return
1978
1979 call get_param(param_file, mdl, "MLE_MLD_DECAY_TIME", cs%MLE_MLD_decay_time, &
1980 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 units="s", default=0., scale=us%s_to_T, do_not_log=.true.)
1983 call openparameterblock(param_file, 'MLE', do_not_log=.true.)
1984 call get_param(param_file, mdl, "USE_BODNER23", use_bodner, &
1985 default=.false., do_not_log=.true.)
1986 call closeparameterblock(param_file)
1987 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 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 units=get_thickness_units(gv), conversion=gv%H_to_MKS)
1993 endif
1994 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 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 units=get_thickness_units(gv), conversion=gv%H_to_MKS)
2000 endif
2001 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 mom_flux_units = "m2 s-2" ; if (.not.gv%Boussinesq) mom_flux_units = "kg m-1 s-2"
2004 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 units=mom_flux_units, conversion=us%L_to_m*gv%H_to_mks*us%s_to_T**2 )
2008 endif
2009
2011
2012!> Returns true if a unit test of functions in MOM_mixedlayer_restrat fail.
2013!! Returns false otherwise.
2014logical 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 print *,'===== mixedlayer_restrat: mixedlayer_restrat_unit_tests =================='
2021
2022 ! Tests of the shape function mu(z)
2023 this_test = &
2024 test_answer(verbose, mu(3.,0.), 0., 'mu(3)=0')
2025 this_test = this_test .or. &
2026 test_answer(verbose, mu(0.,0.), 0., 'mu(0)=0')
2027 this_test = this_test .or. &
2028 test_answer(verbose, mu(-0.25,0.), 0.7946428571428572, 'mu(-0.25)=0.7946...', tol=epsilon(1.))
2029 this_test = this_test .or. &
2030 test_answer(verbose, mu(-0.5,0.), 1., 'mu(-0.5)=1')
2031 this_test = this_test .or. &
2032 test_answer(verbose, mu(-0.75,0.), 0.7946428571428572, 'mu(-0.75)=0.7946...', tol=epsilon(1.))
2033 this_test = this_test .or. &
2034 test_answer(verbose, mu(-1.,0.), 0., 'mu(-1)=0')
2035 this_test = this_test .or. &
2036 test_answer(verbose, mu(-3.,0.), 0., 'mu(-3)=0')
2037 this_test = this_test .or. &
2038 test_answer(verbose, mu(-0.5,0.5), 1., 'mu(-0.5,0.5)=1')
2039 this_test = this_test .or. &
2040 test_answer(verbose, mu(-1.,0.5), 0.25, 'mu(-1,0.5)=0.25')
2041 this_test = this_test .or. &
2042 test_answer(verbose, mu(-1.5,0.5), 0., 'mu(-1.5,0.5)=0')
2043 if (.not. this_test) print '(a)',' Passed tests of mu(z)'
2045
2046 ! Tests of the two time-scale running mean function
2047 this_test = &
2048 test_answer(verbose, rmean2ts(3.,2.,0.,0.,3.), 3., 'rmean2ts(3,2,0,0,3)=3')
2049 this_test = this_test .or. &
2050 test_answer(verbose, rmean2ts(1.,2.,0.,0.,3.), 1., 'rmean2ts(1,2,0,0,3)=1')
2051 this_test = this_test .or. &
2052 test_answer(verbose, rmean2ts(4.,0.,3.,0.,1.), 1., 'rmean2ts(4,0,3,0,1)=1')
2053 this_test = this_test .or. &
2054 test_answer(verbose, rmean2ts(0.,4.,0.,3.,1.), 3., 'rmean2ts(0,4,0,3,1)=3')
2055 if (.not. this_test) print '(a)',' Passed tests of rmean2ts(s,f,g,d,dt)'
2057
2059
2060!> Returns true if any cell of u and u_true are not identical. Returns false otherwise.
2061logical 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 tolerance = 0.0 ; if (present(tol)) tolerance = tol
2071 test_answer = .false.
2072
2073 if (abs(u - u_true) > tolerance) test_answer = .true.
2074 if (test_answer .or. verbose) then
2075 if (test_answer) then
2076 print '(3(a,1pe24.16),1x,a,1x,a)','computed =',u,' correct =',u_true, &
2077 ' err=',u-u_true,' < wrong',label
2078 else
2079 print '(2(a,1pe24.16),1x,a)','computed =',u,' correct =',u_true,label
2080 endif
2081 endif
2082
2083end 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
2252end module mom_mixed_layer_restrat