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