MOM_diagnose_MLD.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!> Provides functions for some diabatic processes such as fraxil, brine rejection,
6!! tendency due to surface flux divergence.
8
9use mom_diag_mediator, only : post_data
11use mom_eos, only : calculate_density, calculate_tfreeze, eos_domain
12use mom_eos, only : calculate_specific_vol_derivs, calculate_density_derivs
13use mom_error_handler, only : mom_error, fatal, warning
14use mom_grid, only : ocean_grid_type
15use mom_interface_heights, only : thickness_to_dz
19
20implicit none ; private
21
22#include <MOM_memory.h>
23
25
26! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
27! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
28! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
29! vary with the Boussinesq approximation, the Boussinesq variant is given first.
30
31contains
32!> Diagnose a mixed layer depth (MLD) determined by a given density difference with the surface.
33!> This routine is appropriate in MOM_diabatic_aux due to its position within the time stepping.
34subroutine diagnosemldbydensitydifference(id_MLD, h, tv, densityDiff, G, GV, US, diagPtr, &
35 ref_h_mld, id_ref_z, id_ref_rho, id_N2subML, id_MLDsq, &
36 dz_subML, MLD_out)
37 type(ocean_grid_type), intent(in) :: g !< Grid type
38 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
39 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
40 integer, intent(in) :: id_mld !< Handle (ID) of MLD diagnostic
41 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
42 intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
43 type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any
44 !! available thermodynamic fields.
45 real, intent(in) :: densitydiff !< Density difference to determine MLD [R ~> kg m-3]
46 type(diag_ctrl), pointer :: diagptr !< Diagnostics structure
47 real, intent(in) :: ref_h_mld !< Depth of the calculated "surface" densisty [Z ~> m]
48 integer, intent(in) :: id_ref_z !< Handle (ID) of reference depth diagnostic
49 integer, intent(in) :: id_ref_rho !< Handle (ID) of reference density diagnostic
50 integer, optional, intent(in) :: id_n2subml !< Optional handle (ID) of subML stratification
51 integer, optional, intent(in) :: id_mldsq !< Optional handle (ID) of squared MLD
52 real, optional, intent(in) :: dz_subml !< The distance over which to calculate N2subML
53 !! or 50 m if missing [Z ~> m]
54 real, dimension(SZI_(G),SZJ_(G)), &
55 optional, intent(out) :: mld_out !< Send MLD to other routines [Z ~> m]
56
57 ! Local variables
58 real, dimension(SZI_(G)) :: deltarhoatkm1, deltarhoatk ! Density differences [R ~> kg m-3].
59 real, dimension(SZI_(G)) :: pref_mld, pref_n2 ! Reference pressures [R L2 T-2 ~> Pa].
60 real, dimension(SZI_(G)) :: href_mld ! Reference depth [Z ~> m].
61 real, dimension(SZI_(G)) :: h_subml, dh_n2 ! Summed thicknesses used in N2 calculation [H ~> m or kg m-2]
62 real, dimension(SZI_(G)) :: dz_n2 ! Summed vertical distance used in N2 calculation [Z ~> m]
63 real, dimension(SZI_(G)) :: t_subml, t_deeper ! Temperatures used in the N2 calculation [C ~> degC].
64 real, dimension(SZI_(G)) :: s_subml, s_deeper ! Salinities used in the N2 calculation [S ~> ppt].
65 real, dimension(SZI_(G)) :: rho_subml, rho_deeper ! Densities used in the N2 calculation [R ~> kg m-3].
66 real, dimension(SZI_(G),SZK_(GV)) :: dz_2d ! Layer thicknesses in depth units [Z ~> m]
67 real, dimension(SZI_(G)) :: dz, dzm1 ! Layer thicknesses associated with interfaces [Z ~> m]
68 real, dimension(SZI_(G)) :: rhosurf ! Density used in finding the mixed layer depth [R ~> kg m-3].
69 real, dimension(SZI_(G), SZJ_(G)) :: z_ref_diag ! The actual depth of the reference density [Z ~> m].
70 real, dimension(SZI_(G), SZJ_(G)) :: mld ! Diagnosed mixed layer depth [Z ~> m].
71 real, dimension(SZI_(G), SZJ_(G)) :: submln2 ! Diagnosed stratification below ML [T-2 ~> s-2].
72 real, dimension(SZI_(G), SZJ_(G)) :: mld2 ! Diagnosed MLD^2 [Z2 ~> m2].
73 logical, dimension(SZI_(G)) :: n2_region_set ! If true, all necessary values for calculating N2
74 ! have been stored already.
75 real :: ge_rho0 ! The gravitational acceleration, sometimes divided by the Boussinesq
76 ! reference density [H T-2 R-1 ~> m4 s-2 kg-1 or m s-2].
77 real :: dz_sub_ml ! Depth below ML over which to diagnose stratification [Z ~> m]
78 real :: afac ! A nondimensional factor [nondim]
79 real :: ddrho ! A density difference [R ~> kg m-3]
80 real :: dddpth ! A depth difference [Z ~> m]
81 real :: rhosurf_k, rhosurf_km1 ! Desisty in the layers below and above the target reference depth [R ~> kg m-3].
82 real, dimension(SZI_(G), SZJ_(G)) :: rhosurf_2d ! The density that is considered the "surface" when calculating
83 ! the MLD. It can be saved as a diagnostic [R ~> kg m-3].
84 integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
85 integer :: i, j, is, ie, js, je, k, nz, id_n2, id_sq
86
87 id_sq = -1 ; if (PRESENT(id_mldsq)) id_sq = id_mldsq
88
89 id_n2 = -1
90 if (present(id_n2subml)) then
91 if (present(dz_subml)) then
92 id_n2 = id_n2subml
93 dz_sub_ml = dz_subml
94 else
95 call mom_error(fatal, "When the diagnostic of the subML stratification is "//&
96 "requested by providing id_N2_subML to diagnoseMLDbyDensityDifference, "//&
97 "the distance over which to calculate that distance must also be provided.")
98 endif
99 endif
100
101 ge_rho0 = gv%g_Earth_Z_T2 / gv%H_to_RZ
102
103 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
104
105 href_mld(:) = ref_h_mld
106 pref_mld(:) = gv%H_to_RZ*gv%g_Earth*ref_h_mld
107 z_ref_diag(:,:) = 0.
108
109 eosdom(:) = eos_domain(g%HI)
110 do j=js,je
111 ! Find the vertical distances across layers.
112 call thickness_to_dz(h, tv, dz_2d, j, g, gv)
113
114 if (pref_mld(is) /= 0.0) then
115 rhosurf(:) = 0.0
116 do i=is,ie
117 dz(i) = 0.5 * dz_2d(i,1) ! Depth of center of surface layer
118 if (dz(i) >= href_mld(i)) then
119 call calculate_density(tv%T(i,j,1), tv%S(i,j,1), pref_mld(i), rhosurf_k, tv%eqn_of_state)
120 rhosurf(i) = rhosurf_k
121 endif
122 enddo
123 do k=2,nz
124 do i=is,ie
125 dzm1(i) = dz(i) ! Depth of center of layer K-1
126 dz(i) = dz(i) + 0.5 * ( dz_2d(i,k) + dz_2d(i,k-1) ) ! Depth of center of layer K
127 dddpth = dz(i) - dzm1(i)
128 if ((rhosurf(i) == 0.) .and. &
129 (dzm1(i) < href_mld(i)) .and. (dz(i) >= href_mld(i))) then
130 afac = ( href_mld(i) - dzm1(i) ) / dddpth
131 z_ref_diag(i,j) = (dz(i) * afac + dzm1(i) * (1. - afac))
132 call calculate_density(tv%T(i,j,k) , tv%S(i,j,k) , pref_mld(i), rhosurf_k, tv%eqn_of_state)
133 call calculate_density(tv%T(i,j,k-1), tv%S(i,j,k-1), pref_mld(i), rhosurf_km1, tv%eqn_of_state)
134 rhosurf(i) = (rhosurf_k * afac + rhosurf_km1 * (1. - afac))
135 h_subml(i) = h(i,j,k)
136 elseif ((rhosurf(i) == 0.) .and. (k >= nz)) then
137 call calculate_density(tv%T(i,j,1), tv%S(i,j,1), pref_mld(i), rhosurf_k, tv%eqn_of_state)
138 rhosurf(i) = rhosurf_k
139 endif
140 enddo
141 enddo
142 do i=is,ie
143 dz(i) = 0.5 * dz_2d(i,1) ! reset dZ to surface depth
144 rhosurf_2d(i,j) = rhosurf(i)
145 deltarhoatk(i) = 0.
146 mld(i,j) = 0.
147 if (id_n2>0) then
148 submln2(i,j) = 0.0
149 dh_n2(i) = 0.0 ; dz_n2(i) = 0.0
150 t_subml(i) = 0.0 ; s_subml(i) = 0.0 ; t_deeper(i) = 0.0 ; s_deeper(i) = 0.0
151 n2_region_set(i) = (g%mask2dT(i,j)<0.5) ! Only need to work on ocean points.
152 endif
153 enddo
154 elseif (pref_mld(is) == 0.0) then
155 rhosurf(:) = 0.0
156 do i=is,ie ; dz(i) = 0.5 * dz_2d(i,1) ; enddo ! Depth of center of surface layer
157 call calculate_density(tv%T(:,j,1), tv%S(:,j,1), pref_mld, rhosurf, tv%eqn_of_state, eosdom)
158 do i=is,ie
159 rhosurf_2d(i,j) = rhosurf(i)
160 deltarhoatk(i) = 0.
161 mld(i,j) = 0.
162 if (id_n2>0) then
163 submln2(i,j) = 0.0
164 h_subml(i) = h(i,j,1) ; dh_n2(i) = 0.0 ; dz_n2(i) = 0.0
165 t_subml(i) = 0.0 ; s_subml(i) = 0.0 ; t_deeper(i) = 0.0 ; s_deeper(i) = 0.0
166 n2_region_set(i) = (g%mask2dT(i,j)<0.5) ! Only need to work on ocean points.
167 endif
168 enddo
169 endif
170
171 do k=2,nz
172 do i=is,ie
173 dzm1(i) = dz(i) ! Depth of center of layer K-1
174 dz(i) = dz(i) + 0.5 * ( dz_2d(i,k) + dz_2d(i,k-1) ) ! Depth of center of layer K
175 enddo
176
177 ! Prepare to calculate stratification, N2, immediately below the mixed layer by finding
178 ! the cells that extend over at least dz_subML.
179 if (id_n2>0) then
180 do i=is,ie
181 if (mld(i,j) == 0.0) then ! Still in the mixed layer.
182 h_subml(i) = h_subml(i) + h(i,j,k)
183 elseif (.not.n2_region_set(i)) then ! This block is below the mixed layer, but N2 has not been found yet.
184 if (dz_n2(i) == 0.0) then ! Record the temperature, salinity, pressure, immediately below the ML
185 t_subml(i) = tv%T(i,j,k) ; s_subml(i) = tv%S(i,j,k)
186 h_subml(i) = h_subml(i) + 0.5 * h(i,j,k) ! Start midway through this layer.
187 dh_n2(i) = 0.5 * h(i,j,k)
188 dz_n2(i) = 0.5 * dz_2d(i,k)
189 elseif (dz_n2(i) + dz_2d(i,k) < dz_sub_ml) then
190 dh_n2(i) = dh_n2(i) + h(i,j,k)
191 dz_n2(i) = dz_n2(i) + dz_2d(i,k)
192 else ! This layer includes the base of the region where N2 is calculated.
193 t_deeper(i) = tv%T(i,j,k) ; s_deeper(i) = tv%S(i,j,k)
194 dh_n2(i) = dh_n2(i) + 0.5 * h(i,j,k)
195 dz_n2(i) = dz_n2(i) + 0.5 * dz_2d(i,k)
196 n2_region_set(i) = .true.
197 endif
198 endif
199 enddo ! i-loop
200 endif ! id_N2>0
201
202 ! Mixed-layer depth, using sigma-0 (surface reference pressure)
203 do i=is,ie ; deltarhoatkm1(i) = deltarhoatk(i) ; enddo ! Store value from previous iteration of K
204 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pref_mld, deltarhoatk, tv%eqn_of_state, eosdom)
205 do i = is, ie
206 deltarhoatk(i) = deltarhoatk(i) - rhosurf(i) ! Density difference between layer K and surface
207 ddrho = deltarhoatk(i) - deltarhoatkm1(i)
208 if ((mld(i,j) == 0.) .and. (ddrho > 0.) .and. &
209 (deltarhoatkm1(i) < densitydiff) .and. (deltarhoatk(i) >= densitydiff)) then
210 afac = ( densitydiff - deltarhoatkm1(i) ) / ddrho
211 mld(i,j) = (dz(i) * afac + dzm1(i) * (1. - afac))
212 endif
213 if (id_sq > 0) mld2(i,j) = mld(i,j)**2
214 enddo ! i-loop
215 enddo ! k-loop
216 do i=is,ie
217 if ((mld(i,j) == 0.) .and. (deltarhoatk(i) < densitydiff)) mld(i,j) = dz(i) ! Mixing goes to the bottom
218 enddo
219
220 if (id_n2>0) then ! Now actually calculate stratification, N2, below the mixed layer.
221 do i=is,ie ; pref_n2(i) = (gv%g_Earth * gv%H_to_RZ) * (h_subml(i) + 0.5*dh_n2(i)) ; enddo
222 ! if ((.not.N2_region_set(i)) .and. (dZ_N2(i) > 0.5*dZ_sub_ML)) then
223 ! ! Use whatever stratification we can, measured over whatever distance is available?
224 ! T_deeper(i) = tv%T(i,j,nz) ; S_deeper(i) = tv%S(i,j,nz)
225 ! N2_region_set(i) = .true.
226 ! endif
227 call calculate_density(t_subml, s_subml, pref_n2, rho_subml, tv%eqn_of_state, eosdom)
228 call calculate_density(t_deeper, s_deeper, pref_n2, rho_deeper, tv%eqn_of_state, eosdom)
229 do i=is,ie ; if ((g%mask2dT(i,j) > 0.0) .and. n2_region_set(i)) then
230 submln2(i,j) = ge_rho0 * (rho_deeper(i) - rho_subml(i)) / dh_n2(i)
231 endif ; enddo
232 endif
233 enddo ! j-loop
234
235 if (id_mld > 0) call post_data(id_mld, mld, diagptr)
236 if (id_n2 > 0) call post_data(id_n2, submln2, diagptr)
237 if (id_sq > 0) call post_data(id_sq, mld2, diagptr)
238
239 if ((id_ref_z > 0) .and. (pref_mld(is)/=0.)) call post_data(id_ref_z, z_ref_diag , diagptr)
240 if (id_ref_rho > 0) call post_data(id_ref_rho, rhosurf_2d , diagptr)
241
242 if (present(mld_out)) then
243 mld_out(:,:) = 0.0
244 mld_out(is:ie,js:je) = mld(is:ie,js:je)
245 endif
246
248
249!> Diagnose a mixed layer depth (MLD) determined by the depth a given energy value would mix.
250!> This routine is appropriate in MOM_diabatic_aux due to its position within the time stepping.
251subroutine diagnosemldbyenergy(id_MLD, h, tv, G, GV, US, Mixing_Energy, k_bounds, diagPtr, OM4_iteration, MLD_out)
252 ! Author: Brandon Reichl
253 ! Date: October 2, 2020
254 ! //
255 ! *Note that gravity is assumed constant everywhere and divided out of all calculations.
256 !
257 ! This code has been written to step through the columns layer by layer, summing the PE
258 ! change inferred by mixing the layer with all layers above. When the change exceeds a
259 ! threshold (determined by input array Mixing_Energy), the code needs to solve for how far
260 ! into this layer the threshold PE change occurs (assuming constant density layers).
261 ! This is expressed here via solving the function F(X) = 0 where:
262 ! F(X) = 0.5 * ( Ca*X^3/(D1+X) + Cb*X^2/(D1+X) + Cc*X/(D1+X) + Dc/(D1+X)
263 ! + Ca2*X^2 + Cb2*X + Cc2)
264 ! where all coefficients are determined by the previous mixed layer depth, the
265 ! density of the previous mixed layer, the present layer thickness, and the present
266 ! layer density. This equation is worked out by computing the total PE assuming constant
267 ! density in the mixed layer as well as in the remaining part of the present layer that is
268 ! not mixed.
269 ! To solve for X in this equation a Newton's method iteration is employed, which
270 ! converges extremely quickly (usually 1 guess) since this equation turns out to be rather
271 ! linear for PE change with increasing X.
272 ! Input parameters:
273 integer, dimension(3), intent(in) :: id_mld !< Energy output diagnostic IDs
274 type(ocean_grid_type), intent(in) :: g !< Grid type
275 type(verticalgrid_type), intent(in) :: gv !< ocean vertical grid structure
276 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
277 real, dimension(3), intent(in) :: mixing_energy !< Energy values for up to 3 MLDs [R Z3 T-2 ~> J m-2]
278 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
279 intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
280 type(thermo_var_ptrs), intent(in) :: tv !< Structure containing pointers to any
281 !! available thermodynamic fields.
282 type(diag_ctrl), pointer :: diagptr !< Diagnostics structure
283 integer, dimension(2), intent(in) :: k_bounds !< vertical interface bounds to apply calculations
284 logical, optional, intent(in) :: om4_iteration !< Uses a legacy version of the MLD iteration
285 !! it is kept to reproduce OM4 output
286 real, dimension(SZI_(G),SZJ_(G)), &
287 optional, intent(out) :: mld_out !< Send MLD to other routines [Z ~> m]
288
289 ! Local variables
290 real, dimension(SZI_(G),SZJ_(G),3) :: mld ! Diagnosed mixed layer depth [Z ~> m].
291 real, dimension(SZK_(GV)+1) :: z_int ! Depths of the interfaces from the surface [Z ~> m]
292 real, dimension(SZI_(G),SZK_(GV)) :: dz ! Layer thicknesses in depth units [Z ~> m]
293 real, dimension(SZI_(G),SZK_(GV)) :: rho_c ! Columns of layer densities [R ~> kg m-3]
294 real, dimension(SZI_(G)) :: pref_mld ! The reference pressure for the mixed layer
295 ! depth calculation [R L2 T-2 ~> Pa]
296 real, dimension(3) :: pe_threshold ! The energy threshold divided by g [R Z2 ~> kg m-1]
297
298 real :: pe_threshold_fraction ! The fractional tolerance of the specified energy
299 ! for the energy used to mix to the diagnosed depth [nondim]
300 real :: h_ml ! The accumulated depth of the mixed layer [Z ~> m]
301 real :: pe ! The cumulative potential energy of the unmixed water column to a depth
302 ! of H_ML, divided by the gravitational acceleration [R Z2 ~> kg m-1]
303 real :: pe_mixed ! The potential energy of the completely mixed water column to a depth
304 ! of H_ML, divided by the gravitational acceleration [R Z2 ~> kg m-1]
305 real :: rhodz_ml ! The depth integrated density of the mixed layer [R Z ~> kg m-2]
306 real :: h_ml_tst ! A new test value for the depth of the mixed layer [Z ~> m]
307 real :: pe_mixed_tst ! The potential energy of the completely mixed water column to a depth
308 ! of H_ML_TST, divided by the gravitational acceleration [R Z2 ~> kg m-1]
309 real :: rhodz_ml_tst ! A test value of the new depth integrated density of the mixed layer [R Z ~> kg m-2]
310 real :: rho_ml ! The average density of the mixed layer [R ~> kg m-3]
311
312 ! These are all temporary variables used to shorten the expressions in the iterations.
313 real :: r1, r2, ca, ca2 ! Some densities [R ~> kg m-3]
314 real :: d1, d2, x, x2 ! Some thicknesses [Z ~> m]
315 real :: cb, cb2 ! A depth integrated density [R Z ~> kg m-2]
316 real :: c, d ! A depth squared [Z2 ~> m2]
317 real :: cc, cc2 ! A density times a depth squared [R Z2 ~> kg m-1]
318 real :: cd ! A density times a depth cubed [R Z3 ~> kg]
319 real :: gx ! A triple integral in depth of density [R Z3 ~> kg]
320 real :: gpx ! The derivative of Gx with x [R Z2 ~> kg m-1]
321 real :: hx ! The vertical integral depth [Z ~> m]
322 real :: ihx ! The inverse of Hx [Z-1 ~> m-1]
323 real :: hpx ! The derivative of Hx with x, since H(x) = constant + x, its derivative is 1. [nondim]
324 real :: ix ! A double integral in depth of density [R Z2 ~> kg m-1]
325 real :: ipx ! The derivative of Ix with x [R Z ~> kg m-2]
326 real :: fgx ! The mixing energy difference from the target [R Z2 ~> kg m-1]
327 real :: fpx ! The derivative of Fgx with x [R Z ~> kg m-2]
328 real :: zr ! An upper (lower) bound for the PE integration in surface (bottom) mixed layer mode [Z ~> m]
329 integer :: k_zr ! Sets the index of Zr
330 real :: pe_dir ! A factor that is used to generalize the iteration for upper and lower mixed layers
331 integer :: k_int ! Controls the direction of the loop to be forward or backward
332 logical :: use_om4_iteration ! A logical to use the OM4_iteration if the optional argument is present
333
334 integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
335 integer :: it, im
336 integer :: i, j, is, ie, js, je, k, nz
337
338 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
339
340 if (present(om4_iteration)) then
341 use_om4_iteration = om4_iteration
342 endif
343
344 pref_mld(:) = 0.0
345 mld(:,:,:) = 0.0
346 pe_threshold_fraction = 1.e-4 !Fixed threshold of 0.01%, could be runtime.
347
348 ! The derivative of H(x) is always 1., so it is moved outside the loops.
349 hpx = 1.
350
351 do im=1,3
352 pe_threshold(im) = mixing_energy(im) / gv%g_Earth_Z_T2
353 enddo
354
355 eosdom(:) = eos_domain(g%HI)
356
357 if (k_bounds(1)<k_bounds(2)) then
358 k_int = 1 ! k_interval is forward in k-space
359 pe_dir = -1. ! A "down" factor so calculations forward and backward use the same algorithm
360 k_zr = k_bounds(1) !top of cell indicated by k_bounds(1)
361 else
362 k_int = -1 ! k_interval is backward in k-space
363 pe_dir = 1. ! An "up" factor so calculations forward and backward use the same algorithm
364 k_zr = k_bounds(1)+1 !bottom of cell indicated by k_bounds(1)
365 endif
366
367 do j=js,je
368 ! Find the vertical distances across layers.
369 call thickness_to_dz(h, tv, dz, j, g, gv)
370
371 if (pe_dir>0) then
372 ! We want to reference pressure to bottom for upward calculation
373 pref_mld(:) = 0.0
374 do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
375 do k=1,nz
376 pref_mld(i) = pref_mld(i) + h(i,j,k)*gv%H_to_RZ*gv%g_Earth
377 enddo
378 endif ; enddo
379 endif
380
381 do k=1,nz
382 call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pref_mld(:), rho_c(:,k), tv%eqn_of_state, eosdom)
383 enddo
384
385 do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
386
387 !We reference everything to the SSH, so that Z_int(1) is defined where Z=0.
388 ! All presently implemented calculations are not sensitive to this choice.
389 ! If "use_OM4_iteration = .true." setting this non-zero would break the iteration
390 z_int(1) = 0.0
391 do k=1,nz
392 z_int(k+1) = z_int(k) - dz(i,k)
393 enddo
394
395 ! Set the reference for the upper (lower) bound of the mixing integral as the surface
396 ! or the bottom depending on the direction of the calculation (as determined by
397 ! the interface bounds k_bounds)
398 zr = z_int(k_zr)
399
400 do im=1,3
401
402 ! Initialize these for each column-wise calculation
403 pe = 0.0
404 rhodz_ml = 0.0
405 h_ml = 0.0
406 rhodz_ml_tst = 0.0
407 h_ml_tst = 0.0
408 pe_mixed = 0.0
409
410 do k=k_bounds(1),k_bounds(2),k_int
411
412 ! This is the unmixed PE cumulative sum in the direction k_int
413 ! The first expression preserves OM4 diagnostic answers, the second is more robust
414 if (use_om4_iteration) then
415 pe = pe + 0.5 * rho_c(i,k) * (z_int(k)**2 - z_int(k+1)**2)
416 else
417 pe = pe + 0.5 * (rho_c(i,k) * dz(i,k)) * (z_int(k) + z_int(k+1))
418 endif
419
420 ! This is the depth and integral of density
421 h_ml_tst = h_ml + dz(i,k)
422 rhodz_ml_tst = rhodz_ml + rho_c(i,k) * dz(i,k)
423
424 ! The average density assuming all layers including this were mixed
425 rho_ml = rhodz_ml_tst/h_ml_tst
426
427 ! The PE assuming all layers including this were mixed
428 ! Zr is the upper (lower) bound of the integral when operating in surface (bottom)
429 ! mixed layer calculation mode.
430 !These are mathematically equivalent, the latter is numerically well-behaved, but the
431 ! former is kept as a comment as it may be more intuitive how it is derived.
432 !PE_Mixed_TST = (0.5 * (Rho_ML*pe_dir)) * ( (Zr + pe_dir*H_ML_TST)**2 - Zr**2.)
433 pe_mixed_tst = (0.5 * (rho_ml*pe_dir)) * (h_ml_tst * (h_ml_tst + 2.0*pe_dir*zr))
434
435 ! Check if we supplied enough energy to mix to this layer
436 if (pe_mixed_tst - pe <= pe_threshold(im)) then
437 h_ml = h_ml_tst
438 rhodz_ml = rhodz_ml_tst
439 else ! If not, we need to solve where the energy ran out within the layer
440 ! This will be done with a Newton's method iteration:
441
442 ! First guess for an iteration using Newton's method
443 x = dz(i,k) * 0.5
444
445 ! We are trying to solve the function:
446 ! F(x) = G(x)/H(x)+I(x)
447 ! for where F(x) = PE+PE_threshold, or equivalently for where
448 ! F(x) = G(x)/H(x)+I(x) - (PE+PE_threshold) = 0
449 ! We also need the derivative of this function for the Newton's method iteration
450 ! F'(x) = (G'(x)H(x)-G(x)H'(x))/H(x)^2 + I'(x)
451 !
452 !For the Surface Boundary Layer:
453 ! The total function F(x) adds the PE of the top layer with some entrained distance X
454 ! to the PE of the bottom layer below the entrained distance:
455 ! (Rho1*D1+Rho2*x)
456 ! PE = ---------------- (Zr^2 - (Zr-D1-x)^2) + Rho2 * ((Zr-D1-x)^2 - (Zr-D1-D2)^2)
457 ! (D1 + x)
458 !
459 ! where Rho1 is the mixed density, D1 is the mixed thickness, Rho2 is the unmixed density,
460 ! D2 is the unmixed thickness, Zr is the top surface height, and x is the fraction of the
461 ! unmixed region that becomes mixed.
462 !
463 !//
464 !G(x) = (Rho1*D1+Rho2*x)*(Zr^2 - (Zr-(D1+x))^2)
465 !
466 ! = -Rho2 * x^3 + (-Rho1*D1-2*Rho2*D1+2*Rho2*Zr)*x^2
467 ! \-Ca-/ \--------Cb----------------/
468 !
469 ! + (-2*Rho1*D1^2+2*Rho1*D1*Zr-Rho2*D1^2+Rho2*2*D1*Zr)*X + Rho1*(-D1^3+2*D1^2*Zr)
470 ! \----------------------Cc----------------------/ \-------Cd----------/
471 !
472 !//
473 !H(x) = D1 + x
474 !
475 !//
476 !I(x) = Rho2 * ((Zr-(D1+x))^2-(Zr-(D1+D2))^2)
477 ! = Rho2 * x^2 + Rho2*(2*D1-2*Zr) * X + Rho2*(D1^2-2*D1*Zr-D2^2+D1^2-2*D1*Zr-2*D2*Zr+2*D1*D2)
478 ! \Ca2/ \-----Cb2-----/ \-------------------Cc2----------------------------/
479 !
480 !
481 !For the Bottom Boundary Layer:
482 ! The total function is relative to Zr as the bottom interface height, so slightly different:
483 ! (Rho1*D1+Rho2*X)
484 ! PE = ---------------- ((Zr+D1+X)^2 - Zr^2) + Rho2 * ((Zr+D1+D2)^2 - (Zr+D1+X)^2)
485 ! (D1 + X)
486 ! These differences propagate through and are accounted for via the factor pe_dir
487 !
488 ! Set these coefficients before the iteration
489 r1 = rhodz_ml / h_ml ! The density of the mixed layer (not including this layer)
490 d1 = h_ml ! The thickness of the mixed layer (not including this layer)
491 r2 = rho_c(i,k) ! The density of this layer to be mixed
492 d2 = dz(i,k) ! The thickness of this layer to be mixed
493
494 ! This sets Zr to "0", which only works for the downward surface mixed layer calculation.
495 ! it should give the same answer at roundoff as the more general expressions below.
496 if (k_int>0 .and. use_om4_iteration) then
497 ca = -(r2)
498 cb = -(r1 * d1 + r2 * (2. * d1))
499 d = d1**2
500 cc = -(r1 * d1 * (2. * d1) + (r2 * d))
501 cd = -r1 * (d1 * d)
502 ca2 = r2
503 cb2 = r2 * (2. * d1)
504 c = d2**2 + d1**2 + 2. * (d1 * d2)
505 d = d1**2
506 cc2 = r2 * (d - c)
507 else
508 ! recall pe_dir = -1 for down, pe_dir = 1 for up.
509 !down Ca = -R2
510 !up Ca = R2
511 ca = pe_dir * r2 ! Density of layer to be mixed
512 !down Cb = -(R1*D1) - 2.*R2*D1 + 2.*Zr*R2
513 !up Cb = (R1*D1) + 2.*R2*D1 + 2.*Zr*R2
514 cb = pe_dir * ( (r1 * d1) + (2. * r2) * ( d1 + zr ) )
515 !down Cc = -2.*R1*D1**2 - R2*D1**2 + 2.*R2*D1*Zr + 2.*Zr*R1*D1
516 !up Cc = 2.*R1*D1**2 + R2*D1**2 + 2.*R2*D1*Zr + 2.*Zr*R1*D1
517 cc = ( pe_dir * d1**2 ) * ( r2 + 2.*r1 ) + ( 2. * ( zr * d1 ) ) * ( r2 + r1 )
518 !down Cd = R1*(-D1**3+2.*D1**2*Zr)
519 !up Cd = R1*( D1**3+2.*D1**2*Zr)
520 cd = ( r1 * d1**2 ) * ( pe_dir * d1 + 2. * zr )
521 !down Ca2 = R2
522 !up Ca2 = -R2
523 ca2 = ( -1. * pe_dir ) * r2
524 !down Cb2 = R2*(2*D1-2*Zr)
525 !up Cb2 = R2*(-2*D1-2*Zr)
526 cb2 = ( 2. * r2 ) * ( (-1.*pe_dir)*d1 - zr )
527 !down Cc2 = R2*(2.*Zr*D2-2.*D1*D2-D2**2)
528 !up Cc2 = R2*(2.*Zr*D2+2.*D1*D2+D2**2)
529 cc2 = ( r2 * d2 ) * ( 2.* zr + pe_dir * ( 2. * d1 + d2 ) )
530 endif
531
532 it=0
533 do while(it<10)!We can iterate up to 10 times
534
535 ! G and its derivative
536 gx = 0.5 * (ca * (x*x*x) + cb * x**2 + cc * x + cd)
537 gpx = 0.5 * (3. * (ca * x**2) + 2. * (cb * x) + cc)
538 ! H, its inverse, and its derivative
539 hx = d1 + x
540 ihx = 1. / hx
541 !Hpx = 1. ! The derivative is always 1 so it was moved outside the loop
542 ! I and its derivative
543 ix = 0.5 * (ca2 * x**2 + cb2 * x + cc2)
544 ipx = 0.5 * (2. * ca2 * x + cb2)
545
546 ! The Function and its derivative:
547 pe_mixed = gx * ihx + ix
548 fgx = pe_mixed - (pe + pe_threshold(im))
549 fpx = (gpx * hx - hpx * gx) * ihx**2 + ipx
550
551 ! Check if our solution is within the threshold bounds, if not update
552 ! using Newton's method. This appears to converge almost always in
553 ! one step because the function is very close to linear in most applications.
554 if (abs(fgx) > pe_threshold(im) * pe_threshold_fraction) then
555 x2 = x - fgx / fpx
556 it = it + 1
557 if (x2 < 0. .or. x2 > dz(i,k)) then
558 ! The iteration seems to be robust, but we need to do something *if*
559 ! things go wrong... How should we treat failed iteration?
560 ! Present solution: Stop trying to compute and just say we can't mix this layer.
561 x=0
562 exit
563 else
564 x = x2
565 endif
566 else
567 exit! Quit the iteration
568 endif
569 enddo
570 h_ml = h_ml + x
571 exit! Quit looping through the column
572 endif
573 enddo
574 mld(i,j,im) = h_ml
575 enddo
576 endif ; enddo
577 enddo
578
579 if (id_mld(1) > 0) call post_data(id_mld(1), mld(:,:,1), diagptr)
580 if (id_mld(2) > 0) call post_data(id_mld(2), mld(:,:,2), diagptr)
581 if (id_mld(3) > 0) call post_data(id_mld(3), mld(:,:,3), diagptr)
582
583 if (present(mld_out)) then
584 mld_out(:,:) = 0.0
585 mld_out(is:ie,js:je) = mld(is:ie,js:je,1)
586 endif
587
588end subroutine diagnosemldbyenergy
589
590!> \namespace mom_diagnose_mld
591!!
592!! This module contains subroutines that apply various diabatic processes. Usually these
593!! subroutines are called from the MOM_diabatic module. All of these routines use appropriate
594!! limiters or logic to work properly with arbitrary layer thicknesses (including massless layers)
595!! and an arbitrarily large timestep.
596!!
597!! The subroutine diagnoseMLDbyDensityDifference diagnoses a mixed layer depth based on a
598!! density difference criterion, and may also estimate the stratification of the water below
599!! this diagnosed mixed layer.
600!!
601!! The subroutine diagnoseMLDbyEnergy diagnoses a mixed layer depth based on a mixing-energy
602!! criterion, as described by Reichl et al., 2022, JGR: Oceans, doi:10.1029/2021JC018140.
603!!
604
605end module mom_diagnose_mld