MOM6
MOM_regularize_layers.F90
1 !> Provides regularization of layers in isopycnal mode
2 module mom_regularize_layers
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
7 use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
8 use mom_diag_mediator, only : time_type, diag_ctrl
9 use mom_domains, only : pass_var
10 use mom_error_handler, only : mom_error, fatal, warning
11 use mom_file_parser, only : get_param, log_version, param_file_type
12 use mom_grid, only : ocean_grid_type
13 use mom_variables, only : thermo_var_ptrs
14 use mom_verticalgrid, only : verticalgrid_type
15 use mom_eos, only : calculate_density, calculate_density_derivs
16 
17 implicit none ; private
18 
19 #include <MOM_memory.h>
20 #undef DEBUG_CODE
21 
22 public regularize_layers, regularize_layers_init
23 
24 !> This control structure holds parameters used by the MOM_regularize_layers module
25 type, public :: regularize_layers_cs ; private
26  logical :: regularize_surface_layers !< If true, vertically restructure the
27  !! near-surface layers when they have too much
28  !! lateral variations to allow for sensible lateral
29  !! barotropic transports.
30  logical :: reg_sfc_detrain !< If true, allow the buffer layers to detrain into the
31  !! interior as a part of the restructuring when
32  !! regularize_surface_layers is true
33  real :: h_def_tol1 !< The value of the relative thickness deficit at
34  !! which to start modifying the structure, 0.5 by
35  !! default (or a thickness ratio of 5.83).
36  real :: h_def_tol2 !< The value of the relative thickness deficit at
37  !! which to the structure modification is in full
38  !! force, now 20% of the way from h_def_tol1 to 1.
39  real :: h_def_tol3 !< The value of the relative thickness deficit at which to start
40  !! detrainment from the buffer layers to the interior, now 30% of
41  !! the way from h_def_tol1 to 1.
42  real :: h_def_tol4 !< The value of the relative thickness deficit at which to do
43  !! detrainment from the buffer layers to the interior at full
44  !! force, now 50% of the way from h_def_tol1 to 1.
45  real :: hmix_min !< The minimum mixed layer thickness [H ~> m or kg m-2].
46  type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
47  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
48  !! regulate the timing of diagnostic output.
49  logical :: debug !< If true, do more thorough checks for debugging purposes.
50 
51  integer :: id_def_rat = -1 !< A diagnostic ID
52  logical :: allow_clocks_in_omp_loops !< If true, clocks can be called from inside loops that
53  !! can be threaded. To run with multiple threads, set to False.
54 #ifdef DEBUG_CODE
55  !>@{ Diagnostic IDs
56  integer :: id_def_rat_2 = -1, id_def_rat_3 = -1
57  integer :: id_def_rat_u = -1, id_def_rat_v = -1
58  integer :: id_e1 = -1, id_e2 = -1, id_e3 = -1
59  integer :: id_def_rat_u_1b = -1, id_def_rat_v_1b = -1
60  integer :: id_def_rat_u_2 = -1, id_def_rat_u_2b = -1
61  integer :: id_def_rat_v_2 = -1, id_def_rat_v_2b = -1
62  integer :: id_def_rat_u_3 = -1, id_def_rat_u_3b = -1
63  integer :: id_def_rat_v_3 = -1, id_def_rat_v_3b = -1
64  !!@}
65 #endif
66 end type regularize_layers_cs
67 
68 !>@{ Clock IDs
69 !! \todo Should these be global?
70 integer :: id_clock_pass, id_clock_eos
71 !!@}
72 
73 contains
74 
75 !> This subroutine partially steps the bulk mixed layer model.
76 !! The following processes are executed, in the order listed.
77 subroutine regularize_layers(h, tv, dt, ea, eb, G, GV, 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_(G)), &
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 ptrs.
85  real, intent(in) :: dt !< Time increment [s].
86  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
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_(G)), &
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(regularize_layers_cs), pointer :: cs !< The control structure returned by a previous
95  !! call to regularize_layers_init.
96  ! Local variables
97  integer :: i, j, k, is, ie, js, je, nz
98 
99  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
100 
101  if (.not. associated(cs)) call mom_error(fatal, "MOM_regularize_layers: "//&
102  "Module must be initialized before it is used.")
103 
104  if (cs%regularize_surface_layers) &
105  call pass_var(h, g%Domain, clock=id_clock_pass)
106 
107  if (cs%regularize_surface_layers) then
108  call regularize_surface(h, tv, dt, ea, eb, g, gv, cs)
109  endif
110 
111 end subroutine regularize_layers
112 
113 !> This subroutine ensures that there is a degree of horizontal smoothness
114 !! in the depths of the near-surface interfaces.
115 subroutine regularize_surface(h, tv, dt, ea, eb, G, GV, CS)
116  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
117  type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
118  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
119  intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2].
120  type(thermo_var_ptrs), intent(inout) :: tv !< A structure containing pointers to any
121  !! available thermodynamic fields. Absent fields
122  !! have NULL ptrs.
123  real, intent(in) :: dt !< Time increment [s].
124  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
125  intent(inout) :: ea !< The amount of fluid moved downward into a
126  !! layer; this should be increased due to mixed
127  !! layer detrainment [H ~> m or kg m-2].
128  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
129  intent(inout) :: eb !< The amount of fluid moved upward into a layer
130  !! this should be increased due to mixed layer
131  !! entrainment [H ~> m or kg m-2].
132  type(regularize_layers_cs), pointer :: CS !< The control structure returned by a previous
133  !! call to regularize_layers_init.
134  ! Local variables
135  real, dimension(SZIB_(G),SZJ_(G)) :: &
136  def_rat_u ! The ratio of the thickness deficit to the minimum depth [nondim].
137  real, dimension(SZI_(G),SZJB_(G)) :: &
138  def_rat_v ! The ratio of the thickness deficit to the minimum depth [nondim].
139  real, dimension(SZI_(G),SZJ_(G)) :: &
140  def_rat_h ! The ratio of the thickness deficit to the minimum depth [nondim].
141  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
142  e ! The interface depths [H ~> m or kg m-2], positive upward.
143 
144 #ifdef DEBUG_CODE
145  real, dimension(SZIB_(G),SZJ_(G)) :: &
146  def_rat_u_1b, def_rat_u_2, def_rat_u_2b, def_rat_u_3, def_rat_u_3b
147  real, dimension(SZI_(G),SZJB_(G)) :: &
148  def_rat_v_1b, def_rat_v_2, def_rat_v_2b, def_rat_v_3, def_rat_v_3b
149  real, dimension(SZI_(G),SZJB_(G)) :: &
150  def_rat_h2, def_rat_h3
151  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
152  ef ! The filtered interface depths [H ~> m or kg m-2], positive upward.
153 #endif
154 
155  real, dimension(SZI_(G),SZK_(G)+1) :: &
156  e_filt, e_2d ! The interface depths [H ~> m or kg m-2], positive upward.
157  real, dimension(SZI_(G),SZK_(G)) :: &
158  h_2d, & ! A 2-d version of h [H ~> m or kg m-2].
159  T_2d, & ! A 2-d version of tv%T [degC].
160  S_2d, & ! A 2-d version of tv%S [ppt].
161  Rcv, & ! A 2-d version of the coordinate density [kg m-3].
162  h_2d_init, & ! The initial value of h_2d [H ~> m or kg m-2].
163  T_2d_init, & ! THe initial value of T_2d [degC].
164  S_2d_init, & ! The initial value of S_2d [ppt].
165  d_eb, & ! The downward increase across a layer in the entrainment from
166  ! below [H ~> m or kg m-2]. The sign convention is that positive values of
167  ! d_eb correspond to a gain in mass by a layer by upward motion.
168  d_ea ! The upward increase across a layer in the entrainment from
169  ! above [H ~> m or kg m-2]. The sign convention is that positive values of
170  ! d_ea mean a net gain in mass by a layer from downward motion.
171  real, dimension(SZI_(G)) :: &
172  p_ref_cv, & ! Reference pressure for the potential density which defines
173  ! the coordinate variable, set to P_Ref [Pa].
174  rcv_tol, & ! A tolerence, relative to the target density differences
175  ! between layers, for detraining into the interior [nondim].
176  h_add_tgt, h_add_tot, &
177  h_tot1, th_tot1, sh_tot1, &
178  h_tot3, th_tot3, sh_tot3, &
179  h_tot2, th_tot2, sh_tot2
180  real, dimension(SZK_(G)) :: &
181  h_prev_1d ! The previous thicknesses [H ~> m or kg m-2].
182  real :: I_dtol ! The inverse of the tolerance changes [nondim].
183  real :: I_dtol34 ! The inverse of the tolerance changes [nondim].
184  real :: h1, h2 ! Temporary thicknesses [H ~> m or kg m-2].
185  real :: e_e, e_w, e_n, e_s ! Temporary interface heights [H ~> m or kg m-2].
186  real :: wt ! The weight of the filted interfaces in setting the targets [nondim].
187  real :: scale ! A scaling factor [nondim].
188  real :: h_neglect ! A thickness that is so small it is usually lost
189  ! in roundoff and can be neglected [H ~> m or kg m-2].
190  real, dimension(SZK_(G)+1) :: &
191  int_flux, int_Tflux, int_Sflux, int_Rflux
192  real :: h_add
193  real :: h_det_tot
194  real :: max_def_rat
195  real :: Rcv_min_det ! The lightest (min) and densest (max) coordinate density
196  real :: Rcv_max_det ! that can detrain into a layer [kg m-3].
197 
198  real :: int_top, int_bot
199  real :: h_predicted
200  real :: h_prev
201  real :: h_deficit
202 
203  logical :: cols_left, ent_any, more_ent_i(SZI_(G)), ent_i(SZI_(G))
204  logical :: det_any, det_i(SZI_(G))
205  logical :: do_j(SZJ_(G)), do_i(SZI_(G)), find_i(SZI_(G))
206  logical :: debug = .false.
207  logical :: fatal_error
208  character(len=256) :: mesg ! Message for error messages.
209  integer :: i, j, k, is, ie, js, je, nz, nkmb, nkml, k1, k2, k3, ks, nz_filt
210 
211  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
212 
213  if (.not. associated(cs)) call mom_error(fatal, "MOM_regularize_layers: "//&
214  "Module must be initialized before it is used.")
215 
216  if (gv%nkml<1) return
217  nkmb = gv%nk_rho_varies ; nkml = gv%nkml
218  if (.not.associated(tv%eqn_of_state)) call mom_error(fatal, &
219  "MOM_regularize_layers: This module now requires the use of temperature and "//&
220  "an equation of state.")
221 
222  h_neglect = gv%H_subroundoff
223  debug = (debug .or. cs%debug)
224 #ifdef DEBUG_CODE
225  debug = .true.
226  if (cs%id_def_rat_2 > 0) then ! Calculate over a slightly larger domain.
227  is = g%isc-1 ; ie = g%iec+1 ; js = g%jsc-1 ; je = g%jec+1
228  endif
229 #endif
230 
231  i_dtol = 1.0 / max(cs%h_def_tol2 - cs%h_def_tol1, 1e-40)
232  i_dtol34 = 1.0 / max(cs%h_def_tol4 - cs%h_def_tol3, 1e-40)
233 
234  p_ref_cv(:) = tv%P_Ref
235 
236  do j=js-1,je+1 ; do i=is-1,ie+1
237  e(i,j,1) = 0.0
238  enddo ; enddo
239  do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
240  e(i,j,k+1) = e(i,j,k) - h(i,j,k)
241  enddo ; enddo ; enddo
242 
243 #ifdef DEBUG_CODE
244  call find_deficit_ratios(e, def_rat_u, def_rat_v, g, gv, cs, def_rat_u_1b, def_rat_v_1b, 1, h)
245 #else
246  call find_deficit_ratios(e, def_rat_u, def_rat_v, g, gv, cs, h=h)
247 #endif
248  ! Determine which columns are problematic
249  do j=js,je ; do_j(j) = .false. ; enddo
250  do j=js,je ; do i=is,ie
251  def_rat_h(i,j) = max(def_rat_u(i-1,j), def_rat_u(i,j), &
252  def_rat_v(i,j-1), def_rat_v(i,j))
253  if (def_rat_h(i,j) > cs%h_def_tol1) do_j(j) = .true.
254  enddo ; enddo
255 
256 #ifdef DEBUG_CODE
257  if ((cs%id_def_rat_3 > 0) .or. (cs%id_e3 > 0) .or. &
258  (cs%id_def_rat_u_3 > 0) .or. (cs%id_def_rat_u_3b > 0) .or. &
259  (cs%id_def_rat_v_3 > 0) .or. (cs%id_def_rat_v_3b > 0) ) then
260  do j=js-1,je+1 ; do i=is-1,ie+1
261  ef(i,j,1) = 0.0
262  enddo ; enddo
263  do k=2,nz+1 ; do j=js,je ; do i=is,ie
264  if (g%mask2dCu(i,j) <= 0.0) then ; e_e = e(i,j,k) ; else
265  e_e = max(e(i+1,j,k) + min(e(i,j,k) - e(i+1,j,nz+1), 0.0), e(i,j,nz+1))
266  endif
267  if (g%mask2dCu(i-1,j) <= 0.0) then ; e_w = e(i,j,k) ; else
268  e_w = max(e(i-1,j,k) + min(e(i,j,k) - e(i-1,j,nz+1), 0.0), e(i,j,nz+1))
269  endif
270  if (g%mask2dCv(i,j) <= 0.0) then ; e_n = e(i,j,k) ; else
271  e_n = max(e(i,j+1,k) + min(e(i,j,k) - e(i,j+1,nz+1), 0.0), e(i,j,nz+1))
272  endif
273  if (g%mask2dCv(i,j-1) <= 0.0) then ; e_s = e(i,j,k) ; else
274  e_s = max(e(i,j-1,k) + min(e(i,j,k) - e(i,j-1,nz+1), 0.0), e(i,j,nz+1))
275  endif
276 
277  wt = 1.0
278  ef(i,j,k) = (1.0 - 0.5*wt) * e(i,j,k) + &
279  wt * 0.125 * ((e_e + e_w) + (e_n + e_s))
280  enddo ; enddo ; enddo
281  call find_deficit_ratios(ef, def_rat_u_3, def_rat_v_3, g, gv, cs, def_rat_u_3b, def_rat_v_3b)
282 
283  ! Determine which columns are problematic
284  do j=js,je ; do i=is,ie
285  def_rat_h3(i,j) = max(def_rat_u_3(i-1,j), def_rat_u_3(i,j), &
286  def_rat_v_3(i,j-1), def_rat_v_3(i,j))
287  enddo ; enddo
288 
289  if (cs%id_e3 > 0) call post_data(cs%id_e3, ef, cs%diag)
290  if (cs%id_def_rat_3 > 0) call post_data(cs%id_def_rat_3, def_rat_h3, cs%diag)
291  if (cs%id_def_rat_u_3 > 0) call post_data(cs%id_def_rat_u_3, def_rat_u_3, cs%diag)
292  if (cs%id_def_rat_u_3b > 0) call post_data(cs%id_def_rat_u_3b, def_rat_u_3b, cs%diag)
293  if (cs%id_def_rat_v_3 > 0) call post_data(cs%id_def_rat_v_3, def_rat_v_3, cs%diag)
294  if (cs%id_def_rat_v_3b > 0) call post_data(cs%id_def_rat_v_3b, def_rat_v_3b, cs%diag)
295  endif
296 #endif
297 
298 
299  ! Now restructure the layers.
300 !$OMP parallel do default(none) shared(is,ie,js,je,nz,do_j,def_rat_h,CS,nkmb,G,GV,&
301 !$OMP e,I_dtol,h,tv,debug,h_neglect,p_ref_cv,ea, &
302 !$OMP eb,id_clock_EOS,nkml) &
303 !$OMP private(d_ea,d_eb,max_def_rat,do_i,nz_filt,e_e,e_w,&
304 !$OMP e_n,e_s,wt,e_filt,e_2d,h_2d,T_2d,S_2d, &
305 !$OMP h_2d_init,T_2d_init,S_2d_init,ent_any, &
306 !$OMP more_ent_i,ent_i,h_add_tgt,h_add_tot, &
307 !$OMP cols_left,h_add,h_prev,ks,det_any,det_i, &
308 !$OMP Rcv_tol,Rcv,k1,k2,h_det_tot,Rcv_min_det, &
309 !$OMP Rcv_max_det,h_deficit,h_tot3,Th_tot3, &
310 !$OMP Sh_tot3,scale,int_top,int_flux,int_Rflux, &
311 !$OMP int_Tflux,int_Sflux,int_bot,h_prev_1d, &
312 !$OMP h_tot1,Th_tot1,Sh_tot1,h_tot2,Th_tot2, &
313 !$OMP Sh_tot2,h_predicted,fatal_error,mesg )
314  do j=js,je ; if (do_j(j)) then
315 
316 ! call cpu_clock_begin(id_clock_EOS)
317 ! call calculate_density_derivs(T(:,1), S(:,1), p_ref_cv, dRcv_dT, dRcv_dS, &
318 ! is, ie-is+1, tv%eqn_of_state)
319 ! call cpu_clock_end(id_clock_EOS)
320 
321  do k=1,nz ; do i=is,ie ; d_ea(i,k) = 0.0 ; d_eb(i,k) = 0.0 ; enddo ; enddo
322 
323  max_def_rat = 0.0
324  do i=is,ie
325  do_i(i) = def_rat_h(i,j) > cs%h_def_tol1
326  if (def_rat_h(i,j) > max_def_rat) max_def_rat = def_rat_h(i,j)
327  enddo
328  nz_filt = nkmb+1 ; if (max_def_rat > cs%h_def_tol3) nz_filt = nz+1
329 
330  ! Find a 2-D 1-2-1 filtered version of e to target. Area weights are
331  ! deliberately omitted here. This is slightly more complicated than a
332  ! simple filter so that the effects of topography are eliminated.
333  do k=1,nz_filt ; do i=is,ie ; if (do_i(i)) then
334  if (g%mask2dCu(i,j) <= 0.0) then ; e_e = e(i,j,k) ; else
335  e_e = max(e(i+1,j,k) + min(e(i,j,k) - e(i+1,j,nz+1), 0.0), &
336  e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
337 
338  endif
339  if (g%mask2dCu(i-1,j) <= 0.0) then ; e_w = e(i,j,k) ; else
340  e_w = max(e(i-1,j,k) + min(e(i,j,k) - e(i-1,j,nz+1), 0.0), &
341  e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
342  endif
343  if (g%mask2dCv(i,j) <= 0.0) then ; e_n = e(i,j,k) ; else
344  e_n = max(e(i,j+1,k) + min(e(i,j,k) - e(i,j+1,nz+1), 0.0), &
345  e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
346  endif
347  if (g%mask2dCv(i,j-1) <= 0.0) then ; e_s = e(i,j,k) ; else
348  e_s = max(e(i,j-1,k) + min(e(i,j,k) - e(i,j-1,nz+1), 0.0), &
349  e(i,j,nz+1) + (nz+1-k)*gv%Angstrom_H)
350  endif
351 
352  wt = max(0.0, min(1.0, i_dtol*(def_rat_h(i,j)-cs%h_def_tol1)))
353 
354  e_filt(i,k) = (1.0 - 0.5*wt) * e(i,j,k) + &
355  wt * 0.125 * ((e_e + e_w) + (e_n + e_s))
356  e_2d(i,k) = e(i,j,k)
357  endif ; enddo ; enddo
358  do k=1,nz ; do i=is,ie
359  h_2d(i,k) = h(i,j,k)
360  t_2d(i,k) = tv%T(i,j,k) ; s_2d(i,k) = tv%S(i,j,k)
361  enddo ; enddo
362 
363  if (debug) then
364  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
365  h_2d_init(i,k) = h(i,j,k)
366  t_2d_init(i,k) = tv%T(i,j,k) ; s_2d_init(i,k) = tv%S(i,j,k)
367  endif ; enddo ; enddo
368  endif
369 
370  ! First, try to entrain from the interior.
371  ent_any = .false.
372  do i=is,ie
373  more_ent_i(i) = .false. ; ent_i(i) = .false.
374  h_add_tgt(i) = 0.0 ; h_add_tot(i) = 0.0
375  if (do_i(i) .and. (e_2d(i,nkmb+1) > e_filt(i,nkmb+1))) then
376  more_ent_i(i) = .true. ; ent_i(i) = .true. ; ent_any = .true.
377  h_add_tgt(i) = e_2d(i,nkmb+1) - e_filt(i,nkmb+1)
378  endif
379  enddo
380 
381  if (ent_any) then
382  do k=nkmb+1,nz
383  cols_left = .false.
384  do i=is,ie ; if (more_ent_i(i)) then
385  if (h_2d(i,k) - gv%Angstrom_H > h_neglect) then
386  if (e_2d(i,nkmb+1)-e_filt(i,nkmb+1) > h_2d(i,k) - gv%Angstrom_H) then
387  h_add = h_2d(i,k) - gv%Angstrom_H
388  h_2d(i,k) = gv%Angstrom_H
389  else
390  h_add = e_2d(i,nkmb+1)-e_filt(i,nkmb+1)
391  h_2d(i,k) = h_2d(i,k) - h_add
392  endif
393  d_eb(i,k-1) = d_eb(i,k-1) + h_add
394  h_add_tot(i) = h_add_tot(i) + h_add
395  e_2d(i,nkmb+1) = e_2d(i,nkmb+1) - h_add
396  h_prev = h_2d(i,nkmb)
397  h_2d(i,nkmb) = h_2d(i,nkmb) + h_add
398 
399  t_2d(i,nkmb) = (h_prev*t_2d(i,nkmb) + h_add*t_2d(i,k)) / h_2d(i,nkmb)
400  s_2d(i,nkmb) = (h_prev*s_2d(i,nkmb) + h_add*s_2d(i,k)) / h_2d(i,nkmb)
401 
402  if ((e_2d(i,nkmb+1) <= e_filt(i,nkmb+1)) .or. &
403  (h_add_tot(i) > 0.6*h_add_tgt(i))) then !### 0.6 is adjustable?.
404  more_ent_i(i) = .false.
405  else
406  cols_left = .true.
407  endif
408  else
409  cols_left = .true.
410  endif
411  endif ; enddo
412  if (.not.cols_left) exit
413  enddo
414 
415  ks = min(k-1,nz-1)
416  do k=ks,nkmb,-1 ; do i=is,ie ; if (ent_i(i)) then
417  d_eb(i,k) = d_eb(i,k) + d_eb(i,k+1)
418  endif ; enddo ; enddo
419  endif ! ent_any
420 
421  ! This is where code to detrain to the interior will go.
422  ! The buffer layers can only detrain water into layers when the buffer
423  ! layer potential density is between (c*Rlay(k-1) + (1-c)*Rlay(k)) and
424  ! (c*Rlay(k+1) + (1-c)*Rlay(k)), where 0.5 <= c < 1.0.
425  ! Do not detrain if the 2-layer deficit ratio is not significant.
426  ! Detrainment must be able to come from all mixed and buffer layers.
427  ! All water is moved out of the buffer layers below before moving from
428  ! a shallower layer (characteristics do not cross).
429  det_any = .false.
430  if ((max_def_rat > cs%h_def_tol3) .and. (cs%reg_sfc_detrain)) then
431  do i=is,ie
432  det_i(i) = .false. ; rcv_tol(i) = 0.0
433  if (do_i(i) .and. (e_2d(i,nkmb+1) < e_filt(i,nkmb+1)) .and. &
434  (def_rat_h(i,j) > cs%h_def_tol3)) then
435  det_i(i) = .true. ; det_any = .true.
436  rcv_tol(i) = min((def_rat_h(i,j) - cs%h_def_tol3), 1.0)
437  endif
438  enddo
439  endif
440  if (det_any) then
441  call cpu_clock_begin(id_clock_eos)
442  do k=1,nkmb
443  call calculate_density(t_2d(:,k),s_2d(:,k),p_ref_cv,rcv(:,k), &
444  is,ie-is+1,tv%eqn_of_state)
445  enddo
446  call cpu_clock_end(id_clock_eos)
447 
448  do i=is,ie ; if (det_i(i)) then
449  k1 = nkmb ; k2 = nz
450  h_det_tot = 0.0
451  do ! This loop is terminated by exits.
452  if (k1 <= 1) exit
453  if (k2 <= nkmb) exit
454  ! ### The 0.6 here should be adjustable? It gives 20% overlap for now.
455  rcv_min_det = gv%Rlay(k2) + 0.6*rcv_tol(i)*(gv%Rlay(k2-1)-gv%Rlay(k2))
456  if (k2 < nz) then
457  rcv_max_det = gv%Rlay(k2) + 0.6*rcv_tol(i)*(gv%Rlay(k2+1)-gv%Rlay(k2))
458  else
459  rcv_max_det = gv%Rlay(nz) + 0.6*rcv_tol(i)*(gv%Rlay(nz)-gv%Rlay(nz-1))
460  endif
461  if (rcv(i,k1) > rcv_max_det) &
462  exit ! All shallower interior layers are too light for detrainment.
463 
464  h_deficit = (e_filt(i,k2)-e_filt(i,k2+1)) - h_2d(i,k2)
465  if ((e_filt(i,k2) > e_2d(i,k1+1)) .and. (h_deficit > 0.0) .and. &
466  (rcv(i,k1) < rcv_max_det) .and. (rcv(i,k1) > rcv_min_det)) then
467  ! Detrainment will occur.
468  h_add = min(e_filt(i,k2) - e_2d(i,k2), h_deficit )
469  if (h_add < h_2d(i,k1)) then
470  ! Only part of layer k1 detrains.
471  if (h_add > 0.0) then
472  h_prev = h_2d(i,k2)
473  h_2d(i,k2) = h_2d(i,k2) + h_add
474  e_2d(i,k2) = e_2d(i,k2+1) + h_2d(i,k2)
475  d_ea(i,k2) = d_ea(i,k2) + h_add
476  ! ### THIS IS UPWIND. IT SHOULD BE HIGHER ORDER...
477  t_2d(i,k2) = (h_prev*t_2d(i,k2) + h_add*t_2d(i,k1)) / h_2d(i,k2)
478  s_2d(i,k2) = (h_prev*s_2d(i,k2) + h_add*s_2d(i,k1)) / h_2d(i,k2)
479  h_det_tot = h_det_tot + h_add
480 
481  h_2d(i,k1) = h_2d(i,k1) - h_add
482  do k3=k1,nkmb ; e_2d(i,k3+1) = e_2d(i,k3) - h_2d(i,k3) ; enddo
483  do k3=k1+1,nkmb ; d_ea(i,k3) = d_ea(i,k3) + h_add ; enddo
484  else
485  if (h_add < 0.0) &
486  call mom_error(fatal, "h_add is negative. Some logic is wrong.")
487  h_add = 0.0 ! This usually should not happen...
488  endif
489 
490  ! Move up to the next target layer.
491  k2 = k2-1
492  if (k2>nkmb+1) e_2d(i,k2) = e_2d(i,k2) + h_det_tot
493  else
494  h_add = h_2d(i,k1)
495  h_prev = h_2d(i,k2)
496  h_2d(i,k2) = h_2d(i,k2) + h_add
497  e_2d(i,k2) = e_2d(i,k2+1) + h_2d(i,k2)
498  d_ea(i,k2) = d_ea(i,k2) + h_add
499  t_2d(i,k2) = (h_prev*t_2d(i,k2) + h_add*t_2d(i,k1)) / h_2d(i,k2)
500  s_2d(i,k2) = (h_prev*s_2d(i,k2) + h_add*s_2d(i,k1)) / h_2d(i,k2)
501  h_det_tot = h_det_tot + h_add
502 
503  h_2d(i,k1) = 0.0
504  do k3=k1,nkmb ; e_2d(i,k3+1) = e_2d(i,k3) - h_2d(i,k3) ; enddo
505  do k3=k1+1,nkmb ; d_ea(i,k3) = d_ea(i,k3) + h_add ; enddo
506 
507  ! Move up to the next source layer.
508  k1 = k1-1
509  endif
510 
511  else
512  ! Move up to the next target layer.
513  k2 = k2-1
514  if (k2>nkmb+1) e_2d(i,k2) = e_2d(i,k2) + h_det_tot
515  endif
516 
517  enddo ! exit terminated loop.
518  endif ; enddo
519  ! ### This could be faster if the deepest k with nonzero d_ea were kept.
520  do k=nz-1,nkmb+1,-1 ; do i=is,ie ; if (det_i(i)) then
521  d_ea(i,k) = d_ea(i,k) + d_ea(i,k+1)
522  endif ; enddo ; enddo
523  endif ! Detrainment to the interior.
524  if (debug) then
525  do i=is,ie ; h_tot3(i) = 0.0 ; th_tot3(i) = 0.0 ; sh_tot3(i) = 0.0 ; enddo
526  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
527  h_tot3(i) = h_tot3(i) + h_2d(i,k)
528  th_tot3(i) = th_tot3(i) + h_2d(i,k) * t_2d(i,k)
529  sh_tot3(i) = sh_tot3(i) + h_2d(i,k) * s_2d(i,k)
530  endif ; enddo ; enddo
531  endif
532 
533  do i=is,ie ; if (do_i(i)) then
534  ! Rescale the interface targets so the depth at the bottom of the deepest
535  ! buffer layer matches.
536  scale = e_2d(i,nkmb+1) / e_filt(i,nkmb+1)
537  do k=2,nkmb+1 ; e_filt(i,k) = e_filt(i,k) * scale ; enddo
538 
539  ! Ensure that layer 1 only has water from layers 1 to nkml and rescale
540  ! the remaining layer thicknesses if necessary.
541  if (e_filt(i,2) < e_2d(i,nkml)) then
542  scale = (e_2d(i,nkml) - e_filt(i,nkmb+1)) / &
543  ((e_filt(i,2) - e_filt(i,nkmb+1)) + h_neglect)
544  do k=3,nkmb
545  e_filt(i,k) = e_filt(i,nkmb+1) + scale * (e_filt(i,k) - e_filt(i,nkmb+1))
546  enddo
547  e_filt(i,2) = e_2d(i,nkml)
548  endif
549 
550  ! Map the water back into the layers.
551  k1 = 1 ; k2 = 1
552  int_top = 0.0
553  do k=1,nkmb+1
554  int_flux(k) = 0.0 ; int_rflux(k) = 0.0
555  int_tflux(k) = 0.0 ; int_sflux(k) = 0.0
556  enddo
557  do k=1,2*nkmb
558  int_bot = max(e_2d(i,k1+1),e_filt(i,k2+1))
559  h_add = int_top - int_bot
560 
561  if (k2 > k1) then
562  do k3=k1+1,k2
563  d_ea(i,k3) = d_ea(i,k3) + h_add
564  int_flux(k3) = int_flux(k3) + h_add
565  int_tflux(k3) = int_tflux(k3) + h_add*t_2d(i,k1)
566  int_sflux(k3) = int_sflux(k3) + h_add*s_2d(i,k1)
567  enddo
568  elseif (k1 > k2) then
569  do k3=k2,k1-1
570  d_eb(i,k3) = d_eb(i,k3) + h_add
571  int_flux(k3+1) = int_flux(k3+1) - h_add
572  int_tflux(k3+1) = int_tflux(k3+1) - h_add*t_2d(i,k1)
573  int_sflux(k3+1) = int_sflux(k3+1) - h_add*s_2d(i,k1)
574  enddo
575  endif
576 
577  if (int_bot <= e_filt(i,k2+1)) then
578  ! Increment the target layer.
579  k2 = k2 + 1
580  elseif (int_bot <= e_2d(i,k1+1)) then
581  ! Increment the source layer.
582  k1 = k1 + 1
583  else
584  call mom_error(fatal, &
585  "Regularize_surface: Could not increment target or source.")
586  endif
587  if ((k1 > nkmb) .or. (k2 > nkmb)) exit
588  int_top = int_bot
589  enddo
590  if (k2 < nkmb) &
591  call mom_error(fatal, "Regularize_surface: Did not assign fluid to layer nkmb.")
592 
593  ! Note that movement of water across the base of the bottommost buffer
594  ! layer has already been dealt with separately.
595  do k=1,nkmb ; h_prev_1d(k) = h_2d(i,k) ; enddo
596  h_2d(i,1) = h_2d(i,1) - int_flux(2)
597  do k=2,nkmb-1
598  h_2d(i,k) = h_2d(i,k) + (int_flux(k) - int_flux(k+1))
599  enddo
600  ! Note that movement of water across the base of the bottommost buffer
601  ! layer has already been dealt with separately.
602  h_2d(i,nkmb) = h_2d(i,nkmb) + int_flux(nkmb)
603 
604  t_2d(i,1) = (t_2d(i,1)*h_prev_1d(1) - int_tflux(2)) / h_2d(i,1)
605  s_2d(i,1) = (s_2d(i,1)*h_prev_1d(1) - int_sflux(2)) / h_2d(i,1)
606  do k=2,nkmb-1
607  t_2d(i,k) = (t_2d(i,k)*h_prev_1d(k) + (int_tflux(k) - int_tflux(k+1))) / h_2d(i,k)
608  s_2d(i,k) = (s_2d(i,k)*h_prev_1d(k) + (int_sflux(k) - int_sflux(k+1))) / h_2d(i,k)
609  enddo
610  t_2d(i,nkmb) = (t_2d(i,nkmb)*h_prev_1d(nkmb) + int_tflux(nkmb) ) / h_2d(i,nkmb)
611  s_2d(i,nkmb) = (s_2d(i,nkmb)*h_prev_1d(nkmb) + int_sflux(nkmb) ) / h_2d(i,nkmb)
612 
613  endif ; enddo ! i-loop
614 
615  ! Copy the interior thicknesses and other fields back to the 3-d arrays.
616  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
617  h(i,j,k) = h_2d(i,k)
618  tv%T(i,j,k) = t_2d(i,k) ; tv%S(i,j,k) = s_2d(i,k)
619  ea(i,j,k) = ea(i,j,k) + d_ea(i,k)
620  eb(i,j,k) = eb(i,j,k) + d_eb(i,k)
621  endif ; enddo ; enddo
622 
623  if (debug) then
624  do i=is,ie ; h_tot1(i) = 0.0 ; th_tot1(i) = 0.0 ; sh_tot1(i) = 0.0 ; enddo
625  do i=is,ie ; h_tot2(i) = 0.0 ; th_tot2(i) = 0.0 ; sh_tot2(i) = 0.0 ; enddo
626 
627  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
628  h_tot1(i) = h_tot1(i) + h_2d_init(i,k)
629  h_tot2(i) = h_tot2(i) + h(i,j,k)
630 
631  th_tot1(i) = th_tot1(i) + h_2d_init(i,k) * t_2d_init(i,k)
632  th_tot2(i) = th_tot2(i) + h(i,j,k) * tv%T(i,j,k)
633  sh_tot1(i) = sh_tot1(i) + h_2d_init(i,k) * s_2d_init(i,k)
634  sh_tot2(i) = sh_tot2(i) + h(i,j,k) * tv%S(i,j,k)
635  if (h(i,j,k) < 0.0) &
636  call mom_error(fatal,"regularize_surface: Negative thicknesses.")
637  if (k==1) then ; h_predicted = h_2d_init(i,k) + (d_eb(i,k) - d_ea(i,k+1))
638  elseif (k==nz) then ; h_predicted = h_2d_init(i,k) + (d_ea(i,k) - d_eb(i,k-1))
639  else
640  h_predicted = h_2d_init(i,k) + ((d_ea(i,k) - d_eb(i,k-1)) + &
641  (d_eb(i,k) - d_ea(i,k+1)))
642  endif
643  if (abs(h(i,j,k) - h_predicted) > max(1e-9*abs(h_predicted),gv%Angstrom_H)) &
644  call mom_error(fatal, "regularize_surface: d_ea mismatch.")
645  endif ; enddo ; enddo
646  do i=is,ie ; if (do_i(i)) then
647  fatal_error = .false.
648  if (abs(h_tot1(i) - h_tot2(i)) > 1e-12*h_tot1(i)) then
649  write(mesg,'(ES11.4," became ",ES11.4," diff ",ES11.4)') &
650  h_tot1(i), h_tot2(i), (h_tot1(i) - h_tot2(i))
651  call mom_error(warning, "regularize_surface: Mass non-conservation."//&
652  trim(mesg), .true.)
653  fatal_error = .true.
654  endif
655  if (abs(th_tot1(i) - th_tot2(i)) > 1e-12*(th_tot1(i)+10.0*h_tot1(i))) then
656  write(mesg,'(ES11.4," became ",ES11.4," diff ",ES11.4," int diff ",ES11.4)') &
657  th_tot1(i), th_tot2(i), (th_tot1(i) - th_tot2(i)), (th_tot1(i) - th_tot3(i))
658  call mom_error(warning, "regularize_surface: Heat non-conservation."//&
659  trim(mesg), .true.)
660  fatal_error = .true.
661  endif
662  if (abs(sh_tot1(i) - sh_tot2(i)) > 1e-12*(sh_tot1(i)+10.0*h_tot1(i))) then
663  write(mesg,'(ES11.4," became ",ES11.4," diff ",ES11.4," int diff ",ES11.4)') &
664  sh_tot1(i), sh_tot2(i), (sh_tot1(i) - sh_tot2(i)), (sh_tot1(i) - sh_tot3(i))
665  call mom_error(warning, "regularize_surface: Salinity non-conservation."//&
666  trim(mesg), .true.)
667  fatal_error = .true.
668  endif
669  if (fatal_error) then
670  write(mesg,'("Error at lat/lon ",2(ES11.4))') g%geoLatT(i,j), g%geoLonT(i,j)
671  call mom_error(fatal, "regularize_surface: Terminating with fatal error. "//&
672  trim(mesg))
673  endif
674  endif ; enddo
675  endif
676 
677  endif ; enddo ! j-loop.
678 
679  if (cs%id_def_rat > 0) call post_data(cs%id_def_rat, def_rat_h, cs%diag)
680 
681 #ifdef DEBUG_CODE
682  if (cs%id_e1 > 0) call post_data(cs%id_e1, e, cs%diag)
683  if (cs%id_def_rat_u > 0) call post_data(cs%id_def_rat_u, def_rat_u, cs%diag)
684  if (cs%id_def_rat_u_1b > 0) call post_data(cs%id_def_rat_u_1b, def_rat_u_1b, cs%diag)
685  if (cs%id_def_rat_v > 0) call post_data(cs%id_def_rat_v, def_rat_v, cs%diag)
686  if (cs%id_def_rat_v_1b > 0) call post_data(cs%id_def_rat_v_1b, def_rat_v_1b, cs%diag)
687 
688  if ((cs%id_def_rat_2 > 0) .or. &
689  (cs%id_def_rat_u_2 > 0) .or. (cs%id_def_rat_u_2b > 0) .or. &
690  (cs%id_def_rat_v_2 > 0) .or. (cs%id_def_rat_v_2b > 0) ) then
691  do j=js-1,je+1 ; do i=is-1,ie+1
692  e(i,j,1) = 0.0
693  enddo ; enddo
694  do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
695  e(i,j,k+1) = e(i,j,k) - h(i,j,k)
696  enddo ; enddo ; enddo
697 
698  call find_deficit_ratios(e, def_rat_u_2, def_rat_v_2, g, gv, cs, def_rat_u_2b, def_rat_v_2b, h=h)
699 
700  ! Determine which columns are problematic
701  do j=js,je ; do i=is,ie
702  def_rat_h2(i,j) = max(def_rat_u_2(i-1,j), def_rat_u_2(i,j), &
703  def_rat_v_2(i,j-1), def_rat_v_2(i,j))
704  enddo ; enddo
705 
706  if (cs%id_def_rat_2 > 0) call post_data(cs%id_def_rat_2, def_rat_h2, cs%diag)
707  if (cs%id_e2 > 0) call post_data(cs%id_e2, e, cs%diag)
708  if (cs%id_def_rat_u_2 > 0) call post_data(cs%id_def_rat_u_2, def_rat_u_2, cs%diag)
709  if (cs%id_def_rat_u_2b > 0) call post_data(cs%id_def_rat_u_2b, def_rat_u_2b, cs%diag)
710  if (cs%id_def_rat_v_2 > 0) call post_data(cs%id_def_rat_v_2, def_rat_v_2, cs%diag)
711  if (cs%id_def_rat_v_2b > 0) call post_data(cs%id_def_rat_v_2b, def_rat_v_2b, cs%diag)
712  endif
713 #endif
714 
715 end subroutine regularize_surface
716 
717 !> This subroutine determines the amount by which the harmonic mean
718 !! thickness at velocity points differ from the arithmetic means, relative to
719 !! the the arithmetic means, after eliminating thickness variations that are
720 !! solely due to topography and aggregating all interior layers into one.
721 subroutine find_deficit_ratios(e, def_rat_u, def_rat_v, G, GV, CS, &
722  def_rat_u_2lay, def_rat_v_2lay, halo, h)
723  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
724  type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
725  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
726  intent(in) :: e !< Interface depths [H ~> m or kg m-2]
727  real, dimension(SZIB_(G),SZJ_(G)), &
728  intent(out) :: def_rat_u !< The thickness deficit ratio at u points,
729  !! [nondim].
730  real, dimension(SZI_(G),SZJB_(G)), &
731  intent(out) :: def_rat_v !< The thickness deficit ratio at v points,
732  !! [nondim].
733  type(regularize_layers_cs), pointer :: CS !< The control structure returned by a
734  !! previous call to regularize_layers_init.
735  real, dimension(SZIB_(G),SZJ_(G)), &
736  optional, intent(out) :: def_rat_u_2lay !< The thickness deficit ratio at u
737  !! points when the mixed and buffer layers
738  !! are aggregated into 1 layer [nondim].
739  real, dimension(SZI_(G),SZJB_(G)), &
740  optional, intent(out) :: def_rat_v_2lay !< The thickness deficit ratio at v
741  !! pointswhen the mixed and buffer layers
742  !! are aggregated into 1 layer [nondim].
743  integer, optional, intent(in) :: halo !< An extra-wide halo size, 0 by default.
744  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
745  optional, intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
746  !! If h is not present, vertical differences
747  !! in interface heights are used instead.
748  ! Local variables
749  real, dimension(SZIB_(G),SZJ_(G)) :: &
750  h_def_u, & ! The vertically summed thickness deficits at u-points [H ~> m or kg m-2].
751  h_norm_u, & ! The vertically summed arithmetic mean thickness by which
752  ! h_def_u is normalized [H ~> m or kg m-2].
753  h_def2_u
754  real, dimension(SZI_(G),SZJB_(G)) :: &
755  h_def_v, & ! The vertically summed thickness deficits at v-points [H ~> m or kg m-2].
756  h_norm_v, & ! The vertically summed arithmetic mean thickness by which
757  ! h_def_v is normalized [H ~> m or kg m-2].
758  h_def2_v
759  real :: h_neglect ! A thickness that is so small it is usually lost
760  ! in roundoff and can be neglected [H ~> m or kg m-2].
761  real :: Hmix_min ! A local copy of CS%Hmix_min [H ~> m or kg m-2].
762  real :: h1, h2 ! Temporary thicknesses [H ~> m or kg m-2].
763  integer :: i, j, k, is, ie, js, je, nz, nkmb
764 
765  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
766  if (present(halo)) then
767  is = g%isc-halo ; ie = g%iec+halo ; js = g%jsc-halo ; je = g%jec+halo
768  endif
769  nkmb = gv%nk_rho_varies
770  h_neglect = gv%H_subroundoff
771  hmix_min = cs%Hmix_min
772 
773  ! Determine which zonal faces are problematic.
774  do j=js,je ; do i=is-1,ie
775  ! Aggregate all water below the mixed and buffer layers for the purposes of
776  ! this diagnostic.
777  h1 = e(i,j,nkmb+1)-e(i,j,nz+1) ; h2 = e(i+1,j,nkmb+1)-e(i+1,j,nz+1)
778  if (e(i,j,nz+1) < e(i+1,j,nz+1)) then
779  if (h1 > h2) h1 = max(e(i,j,nkmb+1)-e(i+1,j,nz+1), h2)
780  elseif (e(i+1,j,nz+1) < e(i,j,nz+1)) then
781  if (h2 > h1) h2 = max(e(i+1,j,nkmb+1)-e(i,j,nz+1), h1)
782  endif
783  h_def_u(i,j) = 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
784  h_norm_u(i,j) = 0.5*(h1+h2)
785  enddo ; enddo
786  if (present(def_rat_u_2lay)) then ; do j=js,je ; do i=is-1,ie
787  ! This is a particular metric of the aggregation into two layers.
788  h1 = e(i,j,1)-e(i,j,nkmb+1) ; h2 = e(i+1,j,1)-e(i+1,j,nkmb+1)
789  if (e(i,j,nkmb+1) < e(i+1,j,nz+1)) then
790  if (h1 > h2) h1 = max(e(i,j,1)-e(i+1,j,nz+1), h2)
791  elseif (e(i+1,j,nkmb+1) < e(i,j,nz+1)) then
792  if (h2 > h1) h2 = max(e(i+1,j,1)-e(i,j,nz+1), h1)
793  endif
794  h_def2_u(i,j) = h_def_u(i,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
795  enddo ; enddo ; endif
796  do k=1,nkmb ; do j=js,je ; do i=is-1,ie
797  if (present(h)) then
798  h1 = h(i,j,k) ; h2 = h(i+1,j,k)
799  else
800  h1 = e(i,j,k)-e(i,j,k+1) ; h2 = e(i+1,j,k)-e(i+1,j,k+1)
801  endif
802  ! Thickness deficits can not arise simply because a layer's bottom is bounded
803  ! by the bathymetry.
804  if (e(i,j,k+1) < e(i+1,j,nz+1)) then
805  if (h1 > h2) h1 = max(e(i,j,k)-e(i+1,j,nz+1), h2)
806  elseif (e(i+1,j,k+1) < e(i,j,nz+1)) then
807  if (h2 > h1) h2 = max(e(i+1,j,k)-e(i,j,nz+1), h1)
808  endif
809  h_def_u(i,j) = h_def_u(i,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
810  h_norm_u(i,j) = h_norm_u(i,j) + 0.5*(h1+h2)
811  enddo ; enddo ; enddo
812  if (present(def_rat_u_2lay)) then ; do j=js,je ; do i=is-1,ie
813  def_rat_u(i,j) = g%mask2dCu(i,j) * h_def_u(i,j) / &
814  (max(hmix_min, h_norm_u(i,j)) + h_neglect)
815  def_rat_u_2lay(i,j) = g%mask2dCu(i,j) * h_def2_u(i,j) / &
816  (max(hmix_min, h_norm_u(i,j)) + h_neglect)
817  enddo ; enddo ; else ; do j=js,je ; do i=is-1,ie
818  def_rat_u(i,j) = g%mask2dCu(i,j) * h_def_u(i,j) / &
819  (max(hmix_min, h_norm_u(i,j)) + h_neglect)
820  enddo ; enddo ; endif
821 
822  ! Determine which meridional faces are problematic.
823  do j=js-1,je ; do i=is,ie
824  ! Aggregate all water below the mixed and buffer layers for the purposes of
825  ! this diagnostic.
826  h1 = e(i,j,nkmb+1)-e(i,j,nz+1) ; h2 = e(i,j+1,nkmb+1)-e(i,j+1,nz+1)
827  if (e(i,j,nz+1) < e(i,j+1,nz+1)) then
828  if (h1 > h2) h1 = max(e(i,j,nkmb+1)-e(i,j+1,nz+1), h2)
829  elseif (e(i,j+1,nz+1) < e(i,j,nz+1)) then
830  if (h2 > h1) h2 = max(e(i,j+1,nkmb+1)-e(i,j,nz+1), h1)
831  endif
832  h_def_v(i,j) = 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
833  h_norm_v(i,j) = 0.5*(h1+h2)
834  enddo ; enddo
835  if (present(def_rat_v_2lay)) then ; do j=js-1,je ; do i=is,ie
836  ! This is a particular metric of the aggregation into two layers.
837  h1 = e(i,j,1)-e(i,j,nkmb+1) ; h2 = e(i,j+1,1)-e(i,j+1,nkmb+1)
838  if (e(i,j,nkmb+1) < e(i,j+1,nz+1)) then
839  if (h1 > h2) h1 = max(e(i,j,1)-e(i,j+1,nz+1), h2)
840  elseif (e(i,j+1,nkmb+1) < e(i,j,nz+1)) then
841  if (h2 > h1) h2 = max(e(i,j+1,1)-e(i,j,nz+1), h1)
842  endif
843  h_def2_v(i,j) = h_def_v(i,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
844  enddo ; enddo ; endif
845  do k=1,nkmb ; do j=js-1,je ; do i=is,ie
846  if (present(h)) then
847  h1 = h(i,j,k) ; h2 = h(i,j+1,k)
848  else
849  h1 = e(i,j,k)-e(i,j,k+1) ; h2 = e(i,j+1,k)-e(i,j+1,k+1)
850  endif
851  ! Thickness deficits can not arise simply because a layer's bottom is bounded
852  ! by the bathymetry.
853  if (e(i,j,k+1) < e(i,j+1,nz+1)) then
854  if (h1 > h2) h1 = max(e(i,j,k)-e(i,j+1,nz+1), h2)
855  elseif (e(i,j+1,k+1) < e(i,j,nz+1)) then
856  if (h2 > h1) h2 = max(e(i,j+1,k)-e(i,j,nz+1), h1)
857  endif
858  h_def_v(i,j) = h_def_v(i,j) + 0.5*(h1-h2)**2 / ((h1 + h2) + h_neglect)
859  h_norm_v(i,j) = h_norm_v(i,j) + 0.5*(h1+h2)
860  enddo ; enddo ; enddo
861  if (present(def_rat_v_2lay)) then ; do j=js-1,je ; do i=is,ie
862  def_rat_v(i,j) = g%mask2dCv(i,j) * h_def_v(i,j) / &
863  (max(hmix_min, h_norm_v(i,j)) + h_neglect)
864  def_rat_v_2lay(i,j) = g%mask2dCv(i,j) * h_def2_v(i,j) / &
865  (max(hmix_min, h_norm_v(i,j)) + h_neglect)
866  enddo ; enddo ; else ; do j=js-1,je ; do i=is,ie
867  def_rat_v(i,j) = g%mask2dCv(i,j) * h_def_v(i,j) / &
868  (max(hmix_min, h_norm_v(i,j)) + h_neglect)
869  enddo ; enddo ; endif
870 
871 end subroutine find_deficit_ratios
872 
873 !> Initializes the regularize_layers control structure
874 subroutine regularize_layers_init(Time, G, GV, param_file, diag, CS)
875  type(time_type), target, intent(in) :: time !< The current model time.
876  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure.
877  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
878  type(param_file_type), intent(in) :: param_file !< A structure to parse for
879  !! run-time parameters.
880  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate
881  !! diagnostic output.
882  type(regularize_layers_cs), pointer :: cs !< A pointer that is set to point to the
883  !! control structure for this module.
884 #include "version_variable.h"
885  character(len=40) :: mdl = "MOM_regularize_layers" ! This module's name.
886  logical :: use_temperature
887  integer :: isd, ied, jsd, jed
888  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
889 
890  if (associated(cs)) then
891  call mom_error(warning, "regularize_layers_init called with an associated"// &
892  "associated control structure.")
893  return
894  else ; allocate(cs) ; endif
895 
896  cs%diag => diag
897  cs%Time => time
898 
899 ! Set default, read and log parameters
900  call log_version(param_file, mdl, version, "")
901  call get_param(param_file, mdl, "REGULARIZE_SURFACE_LAYERS", cs%regularize_surface_layers, &
902  "If defined, vertically restructure the near-surface "//&
903  "layers when they have too much lateral variations to "//&
904  "allow for sensible lateral barotropic transports.", &
905  default=.false.)
906  if (cs%regularize_surface_layers) then
907  call get_param(param_file, mdl, "REGULARIZE_SURFACE_DETRAIN", cs%reg_sfc_detrain, &
908  "If true, allow the buffer layers to detrain into the "//&
909  "interior as a part of the restructuring when "//&
910  "REGULARIZE_SURFACE_LAYERS is true.", default=.true.)
911  endif
912 
913  call get_param(param_file, mdl, "HMIX_MIN", cs%Hmix_min, &
914  "The minimum mixed layer depth if the mixed layer depth "//&
915  "is determined dynamically.", units="m", default=0.0, scale=gv%m_to_H)
916  call get_param(param_file, mdl, "REG_SFC_DEFICIT_TOLERANCE", cs%h_def_tol1, &
917  "The value of the relative thickness deficit at which "//&
918  "to start modifying the layer structure when "//&
919  "REGULARIZE_SURFACE_LAYERS is true.", units="nondim", &
920  default=0.5)
921  cs%h_def_tol2 = 0.2 + 0.8*cs%h_def_tol1
922  cs%h_def_tol3 = 0.3 + 0.7*cs%h_def_tol1
923  cs%h_def_tol4 = 0.5 + 0.5*cs%h_def_tol1
924 
925  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
926 ! if (.not. CS%debug) &
927 ! call get_param(param_file, mdl, "DEBUG_CONSERVATION", CS%debug, &
928 ! "If true, monitor conservation and extrema.", default=.false.)
929 
930  call get_param(param_file, mdl, "ALLOW_CLOCKS_IN_OMP_LOOPS", cs%allow_clocks_in_omp_loops, &
931  "If true, clocks can be called from inside loops that can "//&
932  "be threaded. To run with multiple threads, set to False.", &
933  default=.true.)
934 
935  cs%id_def_rat = register_diag_field('ocean_model', 'deficit_ratio', diag%axesT1, &
936  time, 'Max face thickness deficit ratio', 'nondim')
937 
938 #ifdef DEBUG_CODE
939  cs%id_def_rat_2 = register_diag_field('ocean_model', 'deficit_rat2', diag%axesT1, &
940  time, 'Corrected thickness deficit ratio', 'nondim')
941  cs%id_def_rat_3 = register_diag_field('ocean_model', 'deficit_rat3', diag%axesT1, &
942  time, 'Filtered thickness deficit ratio', 'nondim')
943  cs%id_e1 = register_diag_field('ocean_model', 'er_1', diag%axesTi, &
944  time, 'Intial interface depths before remapping', 'm')
945  cs%id_e2 = register_diag_field('ocean_model', 'er_2', diag%axesTi, &
946  time, 'Intial interface depths after remapping', 'm')
947  cs%id_e3 = register_diag_field('ocean_model', 'er_3', diag%axesTi, &
948  time, 'Intial interface depths filtered', 'm')
949 
950  cs%id_def_rat_u = register_diag_field('ocean_model', 'defrat_u', diag%axesCu1, &
951  time, 'U-point thickness deficit ratio', 'nondim')
952  cs%id_def_rat_u_1b = register_diag_field('ocean_model', 'defrat_u_1b', diag%axesCu1, &
953  time, 'U-point 2-layer thickness deficit ratio', 'nondim')
954  cs%id_def_rat_u_2 = register_diag_field('ocean_model', 'defrat_u_2', diag%axesCu1, &
955  time, 'U-point corrected thickness deficit ratio', 'nondim')
956  cs%id_def_rat_u_2b = register_diag_field('ocean_model', 'defrat_u_2b', diag%axesCu1, &
957  time, 'U-point corrected 2-layer thickness deficit ratio', 'nondim')
958  cs%id_def_rat_u_3 = register_diag_field('ocean_model', 'defrat_u_3', diag%axesCu1, &
959  time, 'U-point filtered thickness deficit ratio', 'nondim')
960  cs%id_def_rat_u_3b = register_diag_field('ocean_model', 'defrat_u_3b', diag%axesCu1, &
961  time, 'U-point filtered 2-layer thickness deficit ratio', 'nondim')
962 
963  cs%id_def_rat_v = register_diag_field('ocean_model', 'defrat_v', diag%axesCv1, &
964  time, 'V-point thickness deficit ratio', 'nondim')
965  cs%id_def_rat_v_1b = register_diag_field('ocean_model', 'defrat_v_1b', diag%axesCv1, &
966  time, 'V-point 2-layer thickness deficit ratio', 'nondim')
967  cs%id_def_rat_v_2 = register_diag_field('ocean_model', 'defrat_v_2', diag%axesCv1, &
968  time, 'V-point corrected thickness deficit ratio', 'nondim')
969  cs%id_def_rat_v_2b = register_diag_field('ocean_model', 'defrat_v_2b', diag%axesCv1, &
970  time, 'V-point corrected 2-layer thickness deficit ratio', 'nondim')
971  cs%id_def_rat_v_3 = register_diag_field('ocean_model', 'defrat_v_3', diag%axesCv1, &
972  time, 'V-point filtered thickness deficit ratio', 'nondim')
973  cs%id_def_rat_v_3b = register_diag_field('ocean_model', 'defrat_v_3b', diag%axesCv1, &
974  time, 'V-point filtered 2-layer thickness deficit ratio', 'nondim')
975 #endif
976 
977  if (cs%allow_clocks_in_omp_loops) then
978  id_clock_eos = cpu_clock_id('(Ocean regularize_layers EOS)', grain=clock_routine)
979  endif
980  id_clock_pass = cpu_clock_id('(Ocean regularize_layers halo updates)', grain=clock_routine)
981 
982 end subroutine regularize_layers_init
983 
984 end module mom_regularize_layers
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:82
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_domains::pass_var
Do a halo update on an array.
Definition: MOM_domains.F90:49
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_regularize_layers::regularize_layers_cs
This control structure holds parameters used by the MOM_regularize_layers module.
Definition: MOM_regularize_layers.F90:25
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_regularize_layers
Provides regularization of layers in isopycnal mode.
Definition: MOM_regularize_layers.F90:2
mom_cpu_clock
Wraps the MPP cpu clock functions.
Definition: MOM_cpu_clock.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_eos::calculate_density_derivs
Calculate the derivatives of density with temperature and salinity from T, S, and P.
Definition: MOM_EOS.F90:70
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:25
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60