MOM_regularize_layers.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 regularization of layers in isopycnal mode
7
8use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
9use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
10use mom_diag_mediator, only : time_type, diag_ctrl
11use mom_domains, only : pass_var
12use mom_error_handler, only : mom_error, fatal, warning
13use mom_file_parser, only : get_param, log_version, param_file_type
14use mom_grid, only : ocean_grid_type
18use mom_eos, only : calculate_density, eos_domain
19
20implicit none ; private
21
22#include <MOM_memory.h>
23
25
26!> This control structure holds parameters used by the MOM_regularize_layers module
27type, public :: regularize_layers_cs ; private
28 logical :: initialized = .false. !< True if this control structure has been initialized.
29 logical :: regularize_surface_layers !< If true, vertically restructure the
30 !! near-surface layers when they have too much
31 !! lateral variations to allow for sensible lateral
32 !! barotropic transports.
33 logical :: reg_sfc_detrain !< If true, allow the buffer layers to detrain into the
34 !! interior as a part of the restructuring when
35 !! regularize_surface_layers is true
36 real :: density_match_tol !< A relative tolerance for how well the densities must match
37 !! with the target densities during detrainment when regularizing
38 !! the near-surface layers [nondim]
39 real :: sufficient_adjustment !< The fraction of the target entrainment of mass to the mixed
40 !! and buffer layers that is enough for one timestep when regularizing
41 !! the near-surface layers [nondim]. No more mass will be sought from
42 !! deeper layers in the interior after this fraction is exceeded.
43 real :: h_def_tol1 !< The value of the relative thickness deficit at
44 !! which to start modifying the structure, 0.5 by
45 !! default (or a thickness ratio of 5.83) [nondim].
46 real :: h_def_tol2 !< The value of the relative thickness deficit at
47 !! which to the structure modification is in full
48 !! force, now 20% of the way from h_def_tol1 to 1 [nondim].
49 real :: h_def_tol3 !< The value of the relative thickness deficit at which to start
50 !! detrainment from the buffer layers to the interior, now 30% of
51 !! the way from h_def_tol1 to 1 [nondim].
52 real :: h_def_tol4 !< The value of the relative thickness deficit at which to do
53 !! detrainment from the buffer layers to the interior at full
54 !! force, now 50% of the way from h_def_tol1 to 1 [nondim].
55 real :: hmix_min !< The minimum mixed layer thickness [H ~> m or kg m-2].
56 type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
57 type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
58 !! regulate the timing of diagnostic output.
59 integer :: answer_date !< The vintage of the order of arithmetic and expressions in this module's
60 !! calculations. Values below 20190101 recover the answers from the
61 !! end of 2018, while higher values use updated and more robust forms
62 !! of the same expressions.
63 logical :: debug !< If true, do more thorough checks for debugging purposes.
64
65 integer :: id_def_rat = -1 !< A diagnostic ID
67
68!>@{ Clock IDs
69!! \todo Should these be global?
70integer :: id_clock_pass
71!>@}
72
73contains
74
75!> This subroutine partially steps the bulk mixed layer model.
76!! The following processes are executed, in the order listed.
77subroutine regularize_layers(h, tv, dt, ea, eb, G, GV, US, CS)
78 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
79 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
80 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
81 intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2].
82 type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
83 !! available thermodynamic fields. Absent fields
84 !! have NULL pointers.
85 real, intent(in) :: dt !< Time increment [T ~> s].
86 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
87 intent(inout) :: ea !< The amount of fluid moved downward into a
88 !! layer; this should be increased due to mixed
89 !! layer detrainment [H ~> m or kg m-2].
90 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
91 intent(inout) :: eb !< The amount of fluid moved upward into a layer
92 !! this should be increased due to mixed layer
93 !! entrainment [H ~> m or kg m-2].
94 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
95 type(regularize_layers_cs), intent(in) :: cs !< Regularize layer control structure
96
97 if (.not. cs%initialized) call mom_error(fatal, "MOM_regularize_layers: "//&
98 "Module must be initialized before it is used.")
99
100 if (cs%regularize_surface_layers) then
101 call pass_var(h, g%Domain, clock=id_clock_pass)
102 call regularize_surface(h, tv, dt, ea, eb, g, gv, us, cs)
103 endif
104
105end subroutine regularize_layers
106
107!> This subroutine ensures that there is a degree of horizontal smoothness
108!! in the depths of the near-surface interfaces.
109subroutine regularize_surface(h, tv, dt, ea, eb, G, GV, US, CS)
110 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
111 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
112 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
113 intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2].
114 type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
115 !! available thermodynamic fields. Absent fields
116 !! have NULL pointers.
117 real, intent(in) :: dt !< Time increment [T ~> s].
118 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
119 intent(inout) :: ea !< The amount of fluid moved downward into a
120 !! layer; this should be increased due to mixed
121 !! layer detrainment [H ~> m or kg m-2].
122 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
123 intent(inout) :: eb !< The amount of fluid moved upward into a layer
124 !! this should be increased due to mixed layer
125 !! entrainment [H ~> m or kg m-2].
126 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
127 type(regularize_layers_cs), intent(in) :: CS !< Regularize layer control structure
128
129 ! Local variables
130 real, dimension(SZIB_(G),SZJ_(G)) :: &
131 def_rat_u ! The ratio of the thickness deficit to the minimum depth [nondim].
132 real, dimension(SZI_(G),SZJB_(G)) :: &
133 def_rat_v ! The ratio of the thickness deficit to the minimum depth [nondim].
134 real, dimension(SZI_(G),SZJ_(G)) :: &
135 def_rat_h ! The ratio of the thickness deficit to the minimum depth [nondim].
136 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
137 e ! The interface depths [H ~> m or kg m-2], positive upward.
138
139 real, dimension(SZI_(G),SZK_(GV)+1) :: &
140 e_filt, e_2d ! The interface depths [H ~> m or kg m-2], positive upward.
141 real, dimension(SZI_(G),SZK_(GV)) :: &
142 h_2d, & ! A 2-d version of h [H ~> m or kg m-2].
143 T_2d, & ! A 2-d version of tv%T [C ~> degC].
144 S_2d, & ! A 2-d version of tv%S [S ~> ppt].
145 Rcv, & ! A 2-d version of the coordinate density [R ~> kg m-3].
146 h_2d_init, & ! The initial value of h_2d [H ~> m or kg m-2].
147 T_2d_init, & ! The initial value of T_2d [C ~> degC].
148 S_2d_init, & ! The initial value of S_2d [S ~> ppt].
149 d_eb, & ! The downward increase across a layer in the entrainment from
150 ! below [H ~> m or kg m-2]. The sign convention is that positive values of
151 ! d_eb correspond to a gain in mass by a layer by upward motion.
152 d_ea ! The upward increase across a layer in the entrainment from
153 ! above [H ~> m or kg m-2]. The sign convention is that positive values of
154 ! d_ea mean a net gain in mass by a layer from downward motion.
155 real, dimension(SZI_(G)) :: &
156 p_ref_cv, & ! Reference pressure for the potential density which defines
157 ! the coordinate variable, set to P_Ref [R L2 T-2 ~> Pa].
158 rcv_tol, & ! A tolerance, relative to the target density differences
159 ! between layers, for detraining into the interior [nondim].
160 h_add_tgt, & ! The target for the thickness to add to the mixed layers [H ~> m or kg m-2]
161 h_add_tot, & ! The net thickness added to the mixed layers [H ~> m or kg m-2]
162 h_tot1, h_tot2, h_tot3, & ! Debugging diagnostics of total thicknesses [H ~> m or kg m-2]
163 th_tot1, th_tot2, th_tot3, & ! Debugging diagnostics of integrated temperatures [C H ~> degC m or degC kg m-2]
164 sh_tot1, sh_tot2, sh_tot3 ! Debugging diagnostics of integrated salinities [S H ~> ppt m or ppt kg m-2]
165 real, dimension(SZK_(GV)) :: &
166 h_prev_1d ! The previous thicknesses [H ~> m or kg m-2].
167 real :: I_dtol ! The inverse of the tolerance changes [nondim].
168 real :: I_dtol34 ! The inverse of the tolerance changes [nondim].
169 real :: e_e, e_w, e_n, e_s ! Temporary interface heights [H ~> m or kg m-2].
170 real :: wt ! The weight of the filtered interfaces in setting the targets [nondim].
171 real :: scale ! A scaling factor [nondim].
172 real :: h_neglect ! A thickness that is so small it is usually lost
173 ! in roundoff and can be neglected [H ~> m or kg m-2].
174 real, dimension(SZK_(GV)+1) :: &
175 int_flux, & ! Mass flux across the interfaces [H ~> m or kg m-2]
176 int_Tflux, & ! Temperature flux across the interfaces [C H ~> degC m or degC kg m-2]
177 int_Sflux ! Salinity flux across the interfaces [S H ~> ppt m or ppt kg m-2]
178 real :: h_add ! The thickness to add to the layers above an interface [H ~> m or kg m-2]
179 real :: h_det_tot ! The total thickness detrained by the mixed layers [H ~> m or kg m-2]
180 real :: max_def_rat ! The maximum value of the ratio of the thickness deficit to the minimum depth [nondim]
181 real :: Rcv_min_det ! The lightest coordinate density that can detrain into a layer [R ~> kg m-3]
182 real :: Rcv_max_det ! The densest coordinate density that can detrain into a layer [R ~> kg m-3]
183
184 real :: int_top, int_bot ! The interface depths above and below a layer [H ~> m or kg m-2], positive upward.
185 real :: h_predicted ! An updated thickness [H ~> m or kg m-2]
186 real :: h_prev ! The previous thickness [H ~> m or kg m-2]
187 real :: h_deficit ! The difference between the layer thickness and the value estimated from the
188 ! filtered interface depths [H ~> m or kg m-2]
189
190 logical :: cols_left, ent_any, more_ent_i(SZI_(G)), ent_i(SZI_(G))
191 logical :: det_any, det_i(SZI_(G))
192 logical :: do_j(SZJ_(G)), do_i(SZI_(G))
193 logical :: debug = .false.
194 logical :: fatal_error
195 character(len=256) :: mesg ! Message for error messages.
196 integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
197 integer :: i, j, k, is, ie, js, je, nz, nkmb, nkml, k1, k2, k3, ks, nz_filt, kmax_d_ea
198
199 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
200
201 if (.not. cs%initialized) call mom_error(fatal, "MOM_regularize_layers: "//&
202 "Module must be initialized before it is used.")
203
204 if (gv%nkml<1) return
205 nkmb = gv%nk_rho_varies ; nkml = gv%nkml
206 if (.not.associated(tv%eqn_of_state)) call mom_error(fatal, &
207 "MOM_regularize_layers: This module now requires the use of temperature and "//&
208 "an equation of state.")
209
210 h_neglect = gv%H_subroundoff
211 debug = (debug .or. cs%debug)
212
213 i_dtol = 1.0 / max(cs%h_def_tol2 - cs%h_def_tol1, 1e-40)
214 i_dtol34 = 1.0 / max(cs%h_def_tol4 - cs%h_def_tol3, 1e-40)
215
216 p_ref_cv(:) = tv%P_Ref
217 eosdom(:) = eos_domain(g%HI)
218
219 do j=js-1,je+1 ; do i=is-1,ie+1
220 e(i,j,1) = 0.0
221 enddo ; enddo
222 do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
223 e(i,j,k+1) = e(i,j,k) - h(i,j,k)
224 enddo ; enddo ; enddo
225
226 call find_deficit_ratios(e, def_rat_u, def_rat_v, g, gv, cs, h)
227
228 ! Determine which columns are problematic
229 do j=js,je ; do_j(j) = .false. ; enddo
230 do j=js,je ; do i=is,ie
231 def_rat_h(i,j) = max(def_rat_u(i-1,j), def_rat_u(i,j), &
232 def_rat_v(i,j-1), def_rat_v(i,j))
233 if (def_rat_h(i,j) > cs%h_def_tol1) do_j(j) = .true.
234 enddo ; enddo
235
236 ! Now restructure the layers.
237 !$OMP parallel do default(private) shared(is,ie,js,je,nz,do_j,def_rat_h,CS,nkmb,G,GV,US, &
238 !$OMP e,I_dtol,h,tv,debug,h_neglect,p_ref_cv,ea, &
239 !$OMP eb,nkml,EOSdom)
240 do j=js,je ; if (do_j(j)) then
241
242 do k=1,nz ; do i=is,ie ; d_ea(i,k) = 0.0 ; d_eb(i,k) = 0.0 ; enddo ; enddo
243 kmax_d_ea = 0
244
245 max_def_rat = 0.0
246 do i=is,ie
247 do_i(i) = def_rat_h(i,j) > cs%h_def_tol1
248 if (def_rat_h(i,j) > max_def_rat) max_def_rat = def_rat_h(i,j)
249 enddo
250 nz_filt = nkmb+1 ; if (max_def_rat > cs%h_def_tol3) nz_filt = nz+1
251
252 ! Find a 2-D 1-2-1 filtered version of e to target. Area weights are
253 ! deliberately omitted here. This is slightly more complicated than a
254 ! simple filter so that the effects of topography are eliminated.
255 do k=1,nz_filt ; do i=is,ie ; if (do_i(i)) then
256 if (g%mask2dCu(i,j) <= 0.0) then ; e_e = e(i,j,k) ; else
257 e_e = max(e(i+1,j,k) + min(e(i,j,k) - e(i+1,j,nz+1), 0.0), &
258 e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
259
260 endif
261 if (g%mask2dCu(i-1,j) <= 0.0) then ; e_w = e(i,j,k) ; else
262 e_w = max(e(i-1,j,k) + min(e(i,j,k) - e(i-1,j,nz+1), 0.0), &
263 e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
264 endif
265 if (g%mask2dCv(i,j) <= 0.0) then ; e_n = e(i,j,k) ; else
266 e_n = max(e(i,j+1,k) + min(e(i,j,k) - e(i,j+1,nz+1), 0.0), &
267 e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
268 endif
269 if (g%mask2dCv(i,j-1) <= 0.0) then ; e_s = e(i,j,k) ; else
270 e_s = max(e(i,j-1,k) + min(e(i,j,k) - e(i,j-1,nz+1), 0.0), &
271 e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
272 endif
273
274 wt = max(0.0, min(1.0, i_dtol*(def_rat_h(i,j)-cs%h_def_tol1)))
275
276 e_filt(i,k) = (1.0 - 0.5*wt) * e(i,j,k) + &
277 wt * 0.125 * ((e_e + e_w) + (e_n + e_s))
278 e_2d(i,k) = e(i,j,k)
279 endif ; enddo ; enddo
280 do k=1,nz ; do i=is,ie
281 h_2d(i,k) = h(i,j,k)
282 t_2d(i,k) = tv%T(i,j,k) ; s_2d(i,k) = tv%S(i,j,k)
283 enddo ; enddo
284
285 if (debug) then
286 do k=1,nz ; do i=is,ie ; if (do_i(i)) then
287 h_2d_init(i,k) = h(i,j,k)
288 t_2d_init(i,k) = tv%T(i,j,k) ; s_2d_init(i,k) = tv%S(i,j,k)
289 endif ; enddo ; enddo
290 endif
291
292 ! First, try to entrain from the interior.
293 ent_any = .false.
294 do i=is,ie
295 more_ent_i(i) = .false. ; ent_i(i) = .false.
296 h_add_tgt(i) = 0.0 ; h_add_tot(i) = 0.0
297 if (do_i(i) .and. (e_2d(i,nkmb+1) > e_filt(i,nkmb+1))) then
298 more_ent_i(i) = .true. ; ent_i(i) = .true. ; ent_any = .true.
299 h_add_tgt(i) = e_2d(i,nkmb+1) - e_filt(i,nkmb+1)
300 endif
301 enddo
302
303 if (ent_any) then
304 do k=nkmb+1,nz
305 cols_left = .false.
306 do i=is,ie ; if (more_ent_i(i)) then
307 if (h_2d(i,k) - gv%Angstrom_H > h_neglect) then
308 if (e_2d(i,nkmb+1)-e_filt(i,nkmb+1) > h_2d(i,k) - gv%Angstrom_H) then
309 h_add = h_2d(i,k) - gv%Angstrom_H
310 h_2d(i,k) = gv%Angstrom_H
311 e_2d(i,nkmb+1) = e_2d(i,nkmb+1) - h_add
312 else
313 h_add = e_2d(i,nkmb+1) - e_filt(i,nkmb+1)
314 h_2d(i,k) = h_2d(i,k) - h_add
315 if (cs%answer_date < 20190101) then
316 e_2d(i,nkmb+1) = e_2d(i,nkmb+1) - h_add
317 else
318 e_2d(i,nkmb+1) = e_filt(i,nkmb+1)
319 endif
320 endif
321 d_eb(i,k-1) = d_eb(i,k-1) + h_add
322 h_add_tot(i) = h_add_tot(i) + h_add
323 h_prev = h_2d(i,nkmb)
324 h_2d(i,nkmb) = h_2d(i,nkmb) + h_add
325
326 t_2d(i,nkmb) = (h_prev*t_2d(i,nkmb) + h_add*t_2d(i,k)) / h_2d(i,nkmb)
327 s_2d(i,nkmb) = (h_prev*s_2d(i,nkmb) + h_add*s_2d(i,k)) / h_2d(i,nkmb)
328
329 if ((e_2d(i,nkmb+1) <= e_filt(i,nkmb+1)) .or. &
330 (h_add_tot(i) > cs%sufficient_adjustment*h_add_tgt(i))) then
331 more_ent_i(i) = .false.
332 else
333 cols_left = .true.
334 endif
335 else
336 cols_left = .true.
337 endif
338 endif ; enddo
339 if (.not.cols_left) exit
340 enddo
341
342 ks = min(k-1,nz-1)
343 do k=ks,nkmb,-1 ; do i=is,ie ; if (ent_i(i)) then
344 d_eb(i,k) = d_eb(i,k) + d_eb(i,k+1)
345 endif ; enddo ; enddo
346 endif ! ent_any
347
348 ! This is where code to detrain to the interior will go.
349 ! The buffer layers can only detrain water into layers when the buffer
350 ! layer potential density is between (c*Rlay(k-1) + (1-c)*Rlay(k)) and
351 ! (c*Rlay(k+1) + (1-c)*Rlay(k)), where 0.5 <= c < 1.0.
352 ! Do not detrain if the 2-layer deficit ratio is not significant.
353 ! Detrainment must be able to come from all mixed and buffer layers.
354 ! All water is moved out of the buffer layers below before moving from
355 ! a shallower layer (characteristics do not cross).
356 det_any = .false.
357 if ((max_def_rat > cs%h_def_tol3) .and. (cs%reg_sfc_detrain)) then
358 do i=is,ie
359 det_i(i) = .false. ; rcv_tol(i) = 0.0
360 if (do_i(i) .and. (e_2d(i,nkmb+1) < e_filt(i,nkmb+1)) .and. &
361 (def_rat_h(i,j) > cs%h_def_tol3)) then
362 det_i(i) = .true. ; det_any = .true.
363 ! The CS%density_match_tol default value of 0.6 gives 20% overlap in acceptable densities.
364 rcv_tol(i) = cs%density_match_tol * min((def_rat_h(i,j) - cs%h_def_tol3), 1.0)
365 endif
366 enddo
367 endif
368 if (det_any) then
369 do k=1,nkmb
370 call calculate_density(t_2d(:,k), s_2d(:,k), p_ref_cv, rcv(:,k), tv%eqn_of_state, eosdom)
371 enddo
372
373 do i=is,ie ; if (det_i(i)) then
374 k1 = nkmb ; k2 = nz
375 h_det_tot = 0.0
376 do ! This loop is terminated by exits.
377 if (k1 <= 1) exit
378 if (k2 <= nkmb) exit
379 rcv_min_det = (gv%Rlay(k2) + rcv_tol(i)*(gv%Rlay(k2-1)-gv%Rlay(k2)))
380 if (k2 < nz) then
381 rcv_max_det = (gv%Rlay(k2) + rcv_tol(i)*(gv%Rlay(k2+1)-gv%Rlay(k2)))
382 else
383 rcv_max_det = (gv%Rlay(nz) + rcv_tol(i)*(gv%Rlay(nz)-gv%Rlay(nz-1)))
384 endif
385 if (rcv(i,k1) > rcv_max_det) &
386 exit ! All shallower interior layers are too light for detrainment.
387
388 h_deficit = (e_filt(i,k2)-e_filt(i,k2+1)) - h_2d(i,k2)
389 if ((e_filt(i,k2) > e_2d(i,k1+1)) .and. (h_deficit > 0.0) .and. &
390 (rcv(i,k1) < rcv_max_det) .and. (rcv(i,k1) > rcv_min_det)) then
391 ! Detrainment will occur.
392 h_add = min(e_filt(i,k2) - e_2d(i,k2), h_deficit )
393 if (h_add < h_2d(i,k1)) then
394 ! Only part of layer k1 detrains.
395 if (h_add > 0.0) then
396 h_prev = h_2d(i,k2)
397 h_2d(i,k2) = h_2d(i,k2) + h_add
398 e_2d(i,k2) = e_2d(i,k2+1) + h_2d(i,k2)
399 d_ea(i,k2) = d_ea(i,k2) + h_add
400 kmax_d_ea = max(kmax_d_ea, k2)
401 ! This is upwind. It should perhaps be higher order...
402 t_2d(i,k2) = (h_prev*t_2d(i,k2) + h_add*t_2d(i,k1)) / h_2d(i,k2)
403 s_2d(i,k2) = (h_prev*s_2d(i,k2) + h_add*s_2d(i,k1)) / h_2d(i,k2)
404 h_det_tot = h_det_tot + h_add
405
406 h_2d(i,k1) = h_2d(i,k1) - h_add
407 do k3=k1,nkmb ; e_2d(i,k3+1) = e_2d(i,k3) - h_2d(i,k3) ; enddo
408 do k3=k1+1,nkmb ; d_ea(i,k3) = d_ea(i,k3) + h_add ; enddo
409 else
410 if (h_add < 0.0) &
411 call mom_error(fatal, "h_add is negative. Some logic is wrong.")
412 h_add = 0.0 ! This usually should not happen...
413 endif
414
415 ! Move up to the next target layer.
416 k2 = k2-1
417 if (k2>nkmb+1) e_2d(i,k2) = e_2d(i,k2) + h_det_tot
418 else
419 h_add = h_2d(i,k1)
420 h_prev = h_2d(i,k2)
421 h_2d(i,k2) = h_2d(i,k2) + h_add
422 e_2d(i,k2) = e_2d(i,k2+1) + h_2d(i,k2)
423 d_ea(i,k2) = d_ea(i,k2) + h_add
424 kmax_d_ea = max(kmax_d_ea, k2)
425 t_2d(i,k2) = (h_prev*t_2d(i,k2) + h_add*t_2d(i,k1)) / h_2d(i,k2)
426 s_2d(i,k2) = (h_prev*s_2d(i,k2) + h_add*s_2d(i,k1)) / h_2d(i,k2)
427 h_det_tot = h_det_tot + h_add
428
429 h_2d(i,k1) = 0.0
430 do k3=k1,nkmb ; e_2d(i,k3+1) = e_2d(i,k3) - h_2d(i,k3) ; enddo
431 do k3=k1+1,nkmb ; d_ea(i,k3) = d_ea(i,k3) + h_add ; enddo
432
433 ! Move up to the next source layer.
434 k1 = k1-1
435 endif
436
437 else
438 ! Move up to the next target layer.
439 k2 = k2-1
440 if (k2>nkmb+1) e_2d(i,k2) = e_2d(i,k2) + h_det_tot
441 endif
442
443 enddo ! exit terminated loop.
444 endif ; enddo
445 do k=kmax_d_ea-1,nkmb+1,-1 ; do i=is,ie ; if (det_i(i)) then
446 d_ea(i,k) = d_ea(i,k) + d_ea(i,k+1)
447 endif ; enddo ; enddo
448 endif ! Detrainment to the interior.
449 if (debug) then
450 do i=is,ie ; h_tot3(i) = 0.0 ; th_tot3(i) = 0.0 ; sh_tot3(i) = 0.0 ; enddo
451 do k=1,nz ; do i=is,ie ; if (do_i(i)) then
452 h_tot3(i) = h_tot3(i) + h_2d(i,k)
453 th_tot3(i) = th_tot3(i) + h_2d(i,k) * t_2d(i,k)
454 sh_tot3(i) = sh_tot3(i) + h_2d(i,k) * s_2d(i,k)
455 endif ; enddo ; enddo
456 endif
457
458 do i=is,ie ; if (do_i(i)) then
459 ! Rescale the interface targets so the depth at the bottom of the deepest
460 ! buffer layer matches.
461 scale = e_2d(i,nkmb+1) / e_filt(i,nkmb+1)
462 do k=2,nkmb+1 ; e_filt(i,k) = e_filt(i,k) * scale ; enddo
463
464 ! Ensure that layer 1 only has water from layers 1 to nkml and rescale
465 ! the remaining layer thicknesses if necessary.
466 if (e_filt(i,2) < e_2d(i,nkml)) then
467 scale = (e_2d(i,nkml) - e_filt(i,nkmb+1)) / &
468 ((e_filt(i,2) - e_filt(i,nkmb+1)) + h_neglect)
469 do k=3,nkmb
470 e_filt(i,k) = e_filt(i,nkmb+1) + scale * (e_filt(i,k) - e_filt(i,nkmb+1))
471 enddo
472 e_filt(i,2) = e_2d(i,nkml)
473 endif
474
475 ! Map the water back into the layers. There are not mixed or buffer layers that are exceedingly
476 ! small compared to the others, so the code here is less prone to roundoff than elsewhere in MOM6.
477 k1 = 1 ; k2 = 1
478 int_top = 0.0
479 do k=1,nkmb+1
480 int_flux(k) = 0.0
481 int_tflux(k) = 0.0 ; int_sflux(k) = 0.0
482 enddo
483 do k=1,2*nkmb
484 int_bot = max(e_2d(i,k1+1),e_filt(i,k2+1))
485 h_add = int_top - int_bot
486
487 if (k2 > k1) then
488 do k3=k1+1,k2
489 d_ea(i,k3) = d_ea(i,k3) + h_add
490 int_flux(k3) = int_flux(k3) + h_add
491 int_tflux(k3) = int_tflux(k3) + h_add*t_2d(i,k1)
492 int_sflux(k3) = int_sflux(k3) + h_add*s_2d(i,k1)
493 enddo
494 elseif (k1 > k2) then
495 do k3=k2,k1-1
496 d_eb(i,k3) = d_eb(i,k3) + h_add
497 int_flux(k3+1) = int_flux(k3+1) - h_add
498 int_tflux(k3+1) = int_tflux(k3+1) - h_add*t_2d(i,k1)
499 int_sflux(k3+1) = int_sflux(k3+1) - h_add*s_2d(i,k1)
500 enddo
501 endif
502
503 if (int_bot <= e_filt(i,k2+1)) then
504 ! Increment the target layer.
505 k2 = k2 + 1
506 elseif (int_bot <= e_2d(i,k1+1)) then
507 ! Increment the source layer.
508 k1 = k1 + 1
509 else
510 call mom_error(fatal, &
511 "Regularize_surface: Could not increment target or source.")
512 endif
513 if ((k1 > nkmb) .or. (k2 > nkmb)) exit
514 int_top = int_bot
515 enddo
516 if (k2 < nkmb) &
517 call mom_error(fatal, "Regularize_surface: Did not assign fluid to layer nkmb.")
518
519 ! Note that movement of water across the base of the bottommost buffer
520 ! layer has already been dealt with separately.
521 do k=1,nkmb ; h_prev_1d(k) = h_2d(i,k) ; enddo
522 h_2d(i,1) = h_2d(i,1) - int_flux(2)
523 do k=2,nkmb-1
524 h_2d(i,k) = h_2d(i,k) + (int_flux(k) - int_flux(k+1))
525 enddo
526 ! Note that movement of water across the base of the bottommost buffer
527 ! layer has already been dealt with separately.
528 h_2d(i,nkmb) = h_2d(i,nkmb) + int_flux(nkmb)
529
530 t_2d(i,1) = (t_2d(i,1)*h_prev_1d(1) - int_tflux(2)) / h_2d(i,1)
531 s_2d(i,1) = (s_2d(i,1)*h_prev_1d(1) - int_sflux(2)) / h_2d(i,1)
532 do k=2,nkmb-1
533 t_2d(i,k) = (t_2d(i,k)*h_prev_1d(k) + (int_tflux(k) - int_tflux(k+1))) / h_2d(i,k)
534 s_2d(i,k) = (s_2d(i,k)*h_prev_1d(k) + (int_sflux(k) - int_sflux(k+1))) / h_2d(i,k)
535 enddo
536 t_2d(i,nkmb) = (t_2d(i,nkmb)*h_prev_1d(nkmb) + int_tflux(nkmb) ) / h_2d(i,nkmb)
537 s_2d(i,nkmb) = (s_2d(i,nkmb)*h_prev_1d(nkmb) + int_sflux(nkmb) ) / h_2d(i,nkmb)
538
539 endif ; enddo ! i-loop
540
541 ! Copy the interior thicknesses and other fields back to the 3-d arrays.
542 do k=1,nz ; do i=is,ie ; if (do_i(i)) then
543 h(i,j,k) = h_2d(i,k)
544 tv%T(i,j,k) = t_2d(i,k) ; tv%S(i,j,k) = s_2d(i,k)
545 ea(i,j,k) = ea(i,j,k) + d_ea(i,k)
546 eb(i,j,k) = eb(i,j,k) + d_eb(i,k)
547 endif ; enddo ; enddo
548
549 if (debug) then
550 do i=is,ie ; h_tot1(i) = 0.0 ; th_tot1(i) = 0.0 ; sh_tot1(i) = 0.0 ; enddo
551 do i=is,ie ; h_tot2(i) = 0.0 ; th_tot2(i) = 0.0 ; sh_tot2(i) = 0.0 ; enddo
552
553 do k=1,nz ; do i=is,ie ; if (do_i(i)) then
554 h_tot1(i) = h_tot1(i) + h_2d_init(i,k)
555 h_tot2(i) = h_tot2(i) + h(i,j,k)
556
557 th_tot1(i) = th_tot1(i) + h_2d_init(i,k) * t_2d_init(i,k)
558 th_tot2(i) = th_tot2(i) + h(i,j,k) * tv%T(i,j,k)
559 sh_tot1(i) = sh_tot1(i) + h_2d_init(i,k) * s_2d_init(i,k)
560 sh_tot2(i) = sh_tot2(i) + h(i,j,k) * tv%S(i,j,k)
561 if (h(i,j,k) < 0.0) &
562 call mom_error(fatal,"regularize_surface: Negative thicknesses.")
563 if (k==1) then ; h_predicted = h_2d_init(i,k) + (d_eb(i,k) - d_ea(i,k+1))
564 elseif (k==nz) then ; h_predicted = h_2d_init(i,k) + (d_ea(i,k) - d_eb(i,k-1))
565 else
566 h_predicted = h_2d_init(i,k) + ((d_ea(i,k) - d_eb(i,k-1)) + &
567 (d_eb(i,k) - d_ea(i,k+1)))
568 endif
569 if (abs(h(i,j,k) - h_predicted) > max(1e-9*abs(h_predicted),gv%Angstrom_H)) &
570 call mom_error(fatal, "regularize_surface: d_ea mismatch.")
571 endif ; enddo ; enddo
572 do i=is,ie ; if (do_i(i)) then
573 fatal_error = .false.
574 if (abs(h_tot1(i) - h_tot2(i)) > 1e-12*h_tot1(i)) then
575 write(mesg,'(ES11.4," became ",ES11.4," diff ",ES11.4)') &
576 h_tot1(i), h_tot2(i), (h_tot1(i) - h_tot2(i))
577 call mom_error(warning, "regularize_surface: Mass non-conservation. "//&
578 trim(mesg), .true.)
579 fatal_error = .true.
580 endif
581 if (abs(th_tot1(i) - th_tot2(i)) > 1e-12*abs(th_tot1(i) + 10.0*us%degC_to_C*h_tot1(i))) then
582 write(mesg,'(ES11.4," became ",ES11.4," diff ",ES11.4," int diff ",ES11.4)') &
583 th_tot1(i), th_tot2(i), (th_tot1(i) - th_tot2(i)), (th_tot1(i) - th_tot3(i))
584 call mom_error(warning, "regularize_surface: Heat non-conservation. "//&
585 trim(mesg), .true.)
586 fatal_error = .true.
587 endif
588 if (abs(sh_tot1(i) - sh_tot2(i)) > 1e-12*abs(sh_tot1(i) + 10.0*us%ppt_to_S*h_tot1(i))) then
589 write(mesg,'(ES11.4," became ",ES11.4," diff ",ES11.4," int diff ",ES11.4)') &
590 sh_tot1(i), sh_tot2(i), (sh_tot1(i) - sh_tot2(i)), (sh_tot1(i) - sh_tot3(i))
591 call mom_error(warning, "regularize_surface: Salinity non-conservation. "//&
592 trim(mesg), .true.)
593 fatal_error = .true.
594 endif
595 if (fatal_error) then
596 write(mesg,'("Error at lat/lon ",2(ES11.4))') g%geoLatT(i,j), g%geoLonT(i,j)
597 call mom_error(fatal, "regularize_surface: Terminating with fatal error. "//&
598 trim(mesg))
599 endif
600 endif ; enddo
601 endif
602
603 endif ; enddo ! j-loop.
604
605 if (cs%id_def_rat > 0) call post_data(cs%id_def_rat, def_rat_h, cs%diag)
606
607end subroutine regularize_surface
608
609!> This subroutine determines the amount by which the harmonic mean
610!! thickness at velocity points differ from the arithmetic means, relative to
611!! the arithmetic means, after eliminating thickness variations that are
612!! solely due to topography and aggregating all interior layers into one.
613subroutine find_deficit_ratios(e, def_rat_u, def_rat_v, G, GV, CS, h)
614 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
615 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
616 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
617 intent(in) :: e !< Interface depths [H ~> m or kg m-2]
618 real, dimension(SZIB_(G),SZJ_(G)), &
619 intent(out) :: def_rat_u !< The thickness deficit ratio at u points,
620 !! [nondim].
621 real, dimension(SZI_(G),SZJB_(G)), &
622 intent(out) :: def_rat_v !< The thickness deficit ratio at v points,
623 !! [nondim].
624 type(regularize_layers_cs), intent(in) :: CS !< Regularize layer control structure
625 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
626 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
627
628 ! Local variables
629 real, dimension(SZIB_(G),SZJ_(G)) :: &
630 h_def_u, & ! The vertically summed thickness deficits at u-points [H ~> m or kg m-2].
631 h_norm_u ! The vertically summed arithmetic mean thickness by which
632 ! h_def_u is normalized [H ~> m or kg m-2].
633 real, dimension(SZI_(G),SZJB_(G)) :: &
634 h_def_v, & ! The vertically summed thickness deficits at v-points [H ~> m or kg m-2].
635 h_norm_v ! The vertically summed arithmetic mean thickness by which
636 ! h_def_v is normalized [H ~> m or kg m-2].
637 real :: h_neglect ! A thickness that is so small it is usually lost
638 ! in roundoff and can be neglected [H ~> m or kg m-2].
639 real :: Hmix_min ! A local copy of CS%Hmix_min [H ~> m or kg m-2].
640 real :: h1, h2 ! Temporary thicknesses [H ~> m or kg m-2].
641 integer :: i, j, k, is, ie, js, je, nz, nkmb
642
643 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
644 nkmb = gv%nk_rho_varies
645 h_neglect = gv%H_subroundoff
646 hmix_min = cs%Hmix_min
647
648 ! Determine which zonal faces are problematic.
649 do j=js,je ; do i=is-1,ie
650 ! Aggregate all water below the mixed and buffer layers for the purposes of
651 ! this diagnostic.
652 h1 = e(i,j,nkmb+1)-e(i,j,nz+1) ; h2 = e(i+1,j,nkmb+1)-e(i+1,j,nz+1)
653 if (e(i,j,nz+1) < e(i+1,j,nz+1)) then
654 if (h1 > h2) h1 = max(e(i,j,nkmb+1)-e(i+1,j,nz+1), h2)
655 elseif (e(i+1,j,nz+1) < e(i,j,nz+1)) then
656 if (h2 > h1) h2 = max(e(i+1,j,nkmb+1)-e(i,j,nz+1), h1)
657 endif
658 h_def_u(i,j) = 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
659 h_norm_u(i,j) = 0.5*(h1+h2)
660 enddo ; enddo
661 do k=1,nkmb ; do j=js,je ; do i=is-1,ie
662 h1 = h(i,j,k) ; h2 = h(i+1,j,k)
663 ! Thickness deficits can not arise simply because a layer's bottom is bounded
664 ! by the bathymetry.
665 if (e(i,j,k+1) < e(i+1,j,nz+1)) then
666 if (h1 > h2) h1 = max(e(i,j,k)-e(i+1,j,nz+1), h2)
667 elseif (e(i+1,j,k+1) < e(i,j,nz+1)) then
668 if (h2 > h1) h2 = max(e(i+1,j,k)-e(i,j,nz+1), h1)
669 endif
670 h_def_u(i,j) = h_def_u(i,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
671 h_norm_u(i,j) = h_norm_u(i,j) + 0.5*(h1+h2)
672 enddo ; enddo ; enddo
673 do j=js,je ; do i=is-1,ie
674 def_rat_u(i,j) = g%mask2dCu(i,j) * h_def_u(i,j) / &
675 (max(hmix_min, h_norm_u(i,j)) + h_neglect)
676 enddo ; enddo
677
678 ! Determine which meridional faces are problematic.
679 do j=js-1,je ; do i=is,ie
680 ! Aggregate all water below the mixed and buffer layers for the purposes of
681 ! this diagnostic.
682 h1 = e(i,j,nkmb+1)-e(i,j,nz+1) ; h2 = e(i,j+1,nkmb+1)-e(i,j+1,nz+1)
683 if (e(i,j,nz+1) < e(i,j+1,nz+1)) then
684 if (h1 > h2) h1 = max(e(i,j,nkmb+1)-e(i,j+1,nz+1), h2)
685 elseif (e(i,j+1,nz+1) < e(i,j,nz+1)) then
686 if (h2 > h1) h2 = max(e(i,j+1,nkmb+1)-e(i,j,nz+1), h1)
687 endif
688 h_def_v(i,j) = 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
689 h_norm_v(i,j) = 0.5*(h1+h2)
690 enddo ; enddo
691 do k=1,nkmb ; do j=js-1,je ; do i=is,ie
692 h1 = h(i,j,k) ; h2 = h(i,j+1,k)
693 ! Thickness deficits can not arise simply because a layer's bottom is bounded
694 ! by the bathymetry.
695 if (e(i,j,k+1) < e(i,j+1,nz+1)) then
696 if (h1 > h2) h1 = max(e(i,j,k)-e(i,j+1,nz+1), h2)
697 elseif (e(i,j+1,k+1) < e(i,j,nz+1)) then
698 if (h2 > h1) h2 = max(e(i,j+1,k)-e(i,j,nz+1), h1)
699 endif
700 h_def_v(i,j) = h_def_v(i,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
701 h_norm_v(i,j) = h_norm_v(i,j) + 0.5*(h1+h2)
702 enddo ; enddo ; enddo
703 do j=js-1,je ; do i=is,ie
704 def_rat_v(i,j) = g%mask2dCv(i,j) * h_def_v(i,j) / &
705 (max(hmix_min, h_norm_v(i,j)) + h_neglect)
706 enddo ; enddo
707
708end subroutine find_deficit_ratios
709
710!> Initializes the regularize_layers control structure
711subroutine regularize_layers_init(Time, G, GV, param_file, diag, CS)
712 type(time_type), target, intent(in) :: time !< The current model time.
713 type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
714 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
715 type(param_file_type), intent(in) :: param_file !< A structure to parse for
716 !! run-time parameters.
717 type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate
718 !! diagnostic output.
719 type(regularize_layers_cs), intent(inout) :: cs !< Regularize layer control structure
720
721# include "version_variable.h"
722 character(len=40) :: mdl = "MOM_regularize_layers" ! This module's name.
723 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags
724 logical :: just_read
725 integer :: isd, ied, jsd, jed
726 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
727
728 cs%initialized = .true.
729
730 cs%diag => diag
731 cs%Time => time
732
733! Set default, read and log parameters
734 call get_param(param_file, mdl, "REGULARIZE_SURFACE_LAYERS", cs%regularize_surface_layers, &
735 default=.false., do_not_log=.true.)
736 call log_version(param_file, mdl, version, "", all_default=.not.cs%regularize_surface_layers)
737 call get_param(param_file, mdl, "REGULARIZE_SURFACE_LAYERS", cs%regularize_surface_layers, &
738 "If defined, vertically restructure the near-surface "//&
739 "layers when they have too much lateral variations to "//&
740 "allow for sensible lateral barotropic transports.", &
741 default=.false.)
742 just_read = .not.cs%regularize_surface_layers
743 if (cs%regularize_surface_layers) then
744 call get_param(param_file, mdl, "REGULARIZE_SURFACE_DETRAIN", cs%reg_sfc_detrain, &
745 "If true, allow the buffer layers to detrain into the "//&
746 "interior as a part of the restructuring when "//&
747 "REGULARIZE_SURFACE_LAYERS is true.", default=.true., do_not_log=just_read)
748 call get_param(param_file, mdl, "REG_SFC_DENSE_MATCH_TOLERANCE", cs%density_match_tol, &
749 "A relative tolerance for how well the densities must match with the target "//&
750 "densities during detrainment when regularizing the near-surface layers. The "//&
751 "default of 0.6 gives 20% overlaps in density", &
752 units="nondim", default=0.6, do_not_log=just_read)
753 call get_param(param_file, mdl, "REG_SFC_SUFFICIENT_ADJ", cs%sufficient_adjustment, &
754 "The fraction of the target entrainment of mass to the mixed and buffer layers "//&
755 "that is enough for one timestep when regularizing the near-surface layers. "//&
756 "No more mass will be sought from deeper layers in the interior after this "//&
757 "fraction is exceeded.", units="nondim", default=0.6, do_not_log=just_read)
758 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
759 "This sets the default value for the various _ANSWER_DATE parameters.", &
760 default=99991231, do_not_log=just_read)
761 call get_param(param_file, mdl, "REGULARIZE_LAYERS_ANSWER_DATE", cs%answer_date, &
762 "The vintage of the order of arithmetic and expressions in the regularize "//&
763 "layers calculations. Values below 20190101 recover the answers from the "//&
764 "end of 2018, while higher values use updated and more robust forms of the "//&
765 "same expressions.", &
766 default=default_answer_date, do_not_log=.not.gv%Boussinesq)
767 if (.not.gv%Boussinesq) cs%answer_date = max(cs%answer_date, 20230701)
768 endif
769
770 call get_param(param_file, mdl, "HMIX_MIN", cs%Hmix_min, &
771 "The minimum mixed layer depth if the mixed layer depth is determined "//&
772 "dynamically.", units="m", default=0.0, scale=gv%m_to_H, do_not_log=just_read)
773 call get_param(param_file, mdl, "REG_SFC_DEFICIT_TOLERANCE", cs%h_def_tol1, &
774 "The value of the relative thickness deficit at which "//&
775 "to start modifying the layer structure when "//&
776 "REGULARIZE_SURFACE_LAYERS is true.", units="nondim", &
777 default=0.5, do_not_log=just_read)
778 cs%h_def_tol2 = 0.2 + 0.8*cs%h_def_tol1
779 cs%h_def_tol3 = 0.3 + 0.7*cs%h_def_tol1
780 cs%h_def_tol4 = 0.5 + 0.5*cs%h_def_tol1
781
782 call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
783! if (.not. CS%debug) &
784! call get_param(param_file, mdl, "DEBUG_CONSERVATION", CS%debug, &
785! "If true, monitor conservation and extrema.", default=.false., do_not_log=just_read)
786
787 if (.not.cs%regularize_surface_layers) return
788
789 cs%id_def_rat = register_diag_field('ocean_model', 'deficit_ratio', diag%axesT1, &
790 time, 'Max face thickness deficit ratio', 'nondim')
791
792 id_clock_pass = cpu_clock_id('(Ocean regularize_layers halo updates)', grain=clock_routine)
793
794end subroutine regularize_layers_init
795
796end module mom_regularize_layers