MOM6
MOM_thickness_diffuse.F90
1 !> Thickness diffusion (or Gent McWilliams)
2 module mom_thickness_diffuse
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_debugging, only : hchksum, uvchksum
7 use mom_diag_mediator, only : post_data, query_averaging_enabled, diag_ctrl
8 use mom_diag_mediator, only : register_diag_field, safe_alloc_ptr, time_type
9 use mom_diag_mediator, only : diag_update_remap_grids
10 use mom_domains, only : pass_var, corner, pass_vector
11 use mom_error_handler, only : mom_error, fatal, warning, is_root_pe
12 use mom_eos, only : calculate_density, calculate_density_derivs, eos_domain
13 use mom_eos, only : calculate_density_second_derivs
14 use mom_file_parser, only : get_param, log_version, param_file_type
15 use mom_grid, only : ocean_grid_type
16 use mom_interface_heights, only : find_eta
17 use mom_isopycnal_slopes, only : vert_fill_ts
18 use mom_lateral_mixing_coeffs, only : varmix_cs
19 use mom_meke_types, only : meke_type
20 use mom_unit_scaling, only : unit_scale_type
21 use mom_variables, only : thermo_var_ptrs, cont_diag_ptrs
22 use mom_verticalgrid, only : verticalgrid_type
23 implicit none ; private
24 
25 #include <MOM_memory.h>
26 
27 public thickness_diffuse, thickness_diffuse_init, thickness_diffuse_end
28 ! public vert_fill_TS
29 public thickness_diffuse_get_kh
30 
31 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
32 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
33 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
34 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
35 
36 !> Control structure for thickness diffusion
37 type, public :: thickness_diffuse_cs ; private
38  real :: khth !< Background interface depth diffusivity [L2 T-1 ~> m2 s-1]
39  real :: khth_slope_cff !< Slope dependence coefficient of Khth [nondim]
40  real :: max_khth_cfl !< Maximum value of the diffusive CFL for thickness diffusion
41  real :: khth_min !< Minimum value of Khth [L2 T-1 ~> m2 s-1]
42  real :: khth_max !< Maximum value of Khth [L2 T-1 ~> m2 s-1], or 0 for no max
43  real :: slope_max !< Slopes steeper than slope_max are limited in some way [nondim].
44  real :: kappa_smooth !< Vertical diffusivity used to interpolate more
45  !! sensible values of T & S into thin layers [Z2 T-1 ~> m2 s-1].
46  logical :: thickness_diffuse !< If true, interfaces heights are diffused.
47  logical :: use_fgnv_streamfn !< If true, use the streamfunction formulation of
48  !! Ferrari et al., 2010, which effectively emphasizes
49  !! graver vertical modes by smoothing in the vertical.
50  real :: fgnv_scale !< A coefficient scaling the vertical smoothing term in the
51  !! Ferrari et al., 2010, streamfunction formulation [nondim].
52  real :: fgnv_c_min !< A minimum wave speed used in the Ferrari et al., 2010,
53  !! streamfunction formulation [L T-1 ~> m s-1].
54  real :: n2_floor !< A floor for Brunt-Vasaila frequency in the Ferrari et al., 2010,
55  !! streamfunction formulation [T-2 ~> s-2].
56  logical :: detangle_interfaces !< If true, add 3-d structured interface height
57  !! diffusivities to horizontally smooth jagged layers.
58  real :: detangle_time !< If detangle_interfaces is true, this is the
59  !! timescale over which maximally jagged grid-scale
60  !! thickness variations are suppressed [T ~> s]. This must be
61  !! longer than DT, or 0 (the default) to use DT.
62  integer :: nkml !< number of layers within mixed layer
63  logical :: debug !< write verbose checksums for debugging purposes
64  logical :: use_gme_thickness_diffuse !< If true, passes GM coefficients to MOM_hor_visc for use
65  !! with GME closure.
66  logical :: meke_geometric !< If true, uses the GM coefficient formulation from the GEOMETRIC
67  !! framework (Marshall et al., 2012)
68  real :: meke_geometric_alpha!< The nondimensional coefficient governing the efficiency of
69  !! the GEOMETRIC thickness difussion [nondim]
70  real :: meke_geometric_epsilon !< Minimum Eady growth rate for the GEOMETRIC thickness
71  !! diffusivity [T-1 ~> s-1].
72  logical :: meke_geom_answers_2018 !< If true, use expressions in the MEKE_GEOMETRIC calculation
73  !! that recover the answers from the original implementation.
74  !! Otherwise, use expressions that satisfy rotational symmetry.
75  logical :: use_kh_in_meke !< If true, uses the thickness diffusivity calculated here to diffuse MEKE.
76  logical :: gm_src_alt !< If true, use the GM energy conversion form S^2*N^2*kappa rather
77  !! than the streamfunction for the GM source term.
78  logical :: use_gm_work_bug !< If true, use the incorrect sign for the
79  !! top-level work tendency on the top layer.
80  real :: stanley_det_coeff !< The coefficient correlating SGS temperature variance with the mean
81  !! temperature gradient in the deterministic part of the Stanley parameterization.
82  !! Negative values disable the scheme." [nondim]
83 
84  type(diag_ctrl), pointer :: diag => null() !< structure used to regulate timing of diagnostics
85  real, pointer :: gmwork(:,:) => null() !< Work by thickness diffusivity [R Z L2 T-3 ~> W m-2]
86  real, pointer :: diagslopex(:,:,:) => null() !< Diagnostic: zonal neutral slope [nondim]
87  real, pointer :: diagslopey(:,:,:) => null() !< Diagnostic: zonal neutral slope [nondim]
88 
89  real, dimension(:,:,:), pointer :: &
90  kh_u_gme => null(), & !< interface height diffusivities in u-columns [L2 T-1 ~> m2 s-1]
91  kh_v_gme => null() !< interface height diffusivities in v-columns [L2 T-1 ~> m2 s-1]
92 
93  !>@{
94  !! Diagnostic identifier
95  integer :: id_uhgm = -1, id_vhgm = -1, id_gmwork = -1
96  integer :: id_kh_u = -1, id_kh_v = -1, id_kh_t = -1
97  integer :: id_kh_u1 = -1, id_kh_v1 = -1, id_kh_t1 = -1
98  integer :: id_slope_x = -1, id_slope_y = -1
99  integer :: id_sfn_unlim_x = -1, id_sfn_unlim_y = -1, id_sfn_x = -1, id_sfn_y = -1
100  !>@}
101 end type thickness_diffuse_cs
102 
103 contains
104 
105 !> Calculates thickness diffusion coefficients and applies thickness diffusion to layer
106 !! thicknesses, h. Diffusivities are limited to ensure stability.
107 !! Also returns along-layer mass fluxes used in the continuity equation.
108 subroutine thickness_diffuse(h, uhtr, vhtr, tv, dt, G, GV, US, MEKE, VarMix, CDp, CS)
109  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
110  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
111  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
112  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
113  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uhtr !< Accumulated zonal mass flux
114  !! [L2 H ~> m3 or kg]
115  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vhtr !< Accumulated meridional mass flux
116  !! [L2 H ~> m3 or kg]
117  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
118  real, intent(in) :: dt !< Time increment [T ~> s]
119  type(meke_type), pointer :: meke !< MEKE control structure
120  type(varmix_cs), pointer :: varmix !< Variable mixing coefficients
121  type(cont_diag_ptrs), intent(inout) :: cdp !< Diagnostics for the continuity equation
122  type(thickness_diffuse_cs), pointer :: cs !< Control structure for thickness diffusion
123  ! Local variables
124  real :: e(szi_(g), szj_(g), szk_(g)+1) ! heights of interfaces, relative to mean
125  ! sea level [Z ~> m], positive up.
126  real :: uhd(szib_(g), szj_(g), szk_(g)) ! Diffusive u*h fluxes [L2 H T-1 ~> m3 s-1 or kg s-1]
127  real :: vhd(szi_(g), szjb_(g), szk_(g)) ! Diffusive v*h fluxes [L2 H T-1 ~> m3 s-1 or kg s-1]
128 
129  real, dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: &
130  kh_u, & ! interface height diffusivities in u-columns [L2 T-1 ~> m2 s-1]
131  int_slope_u ! A nondimensional ratio from 0 to 1 that gives the relative
132  ! weighting of the interface slopes to that calculated also
133  ! using density gradients at u points. The physically correct
134  ! slopes occur at 0, while 1 is used for numerical closures [nondim].
135  real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: &
136  kh_v, & ! interface height diffusivities in v-columns [L2 T-1 ~> m2 s-1]
137  int_slope_v ! A nondimensional ratio from 0 to 1 that gives the relative
138  ! weighting of the interface slopes to that calculated also
139  ! using density gradients at v points. The physically correct
140  ! slopes occur at 0, while 1 is used for numerical closures [nondim].
141  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
142  kh_t ! diagnosed diffusivity at tracer points [L2 T-1 ~> m2 s-1]
143 
144  real, dimension(SZIB_(G), SZJ_(G)) :: &
145  kh_u_cfl ! The maximum stable interface height diffusivity at u grid points [L2 T-1 ~> m2 s-1]
146  real, dimension(SZI_(G), SZJB_(G)) :: &
147  kh_v_cfl ! The maximum stable interface height diffusivity at v grid points [L2 T-1 ~> m2 s-1]
148  real :: khth_loc_u(szib_(g), szj_(g))
149  real :: khth_loc_v(szi_(g), szjb_(g))
150  real :: khth_loc(szib_(g), szjb_(g)) ! locally calculated thickness diffusivity [L2 T-1 ~> m2 s-1]
151  real :: h_neglect ! A thickness that is so small it is usually lost
152  ! in roundoff and can be neglected [H ~> m or kg m-2].
153  real, dimension(:,:), pointer :: cg1 => null() !< Wave speed [L T-1 ~> m s-1]
154  logical :: use_varmix, resoln_scaled, depth_scaled, use_stored_slopes, khth_use_ebt_struct, use_visbeck
155  logical :: use_qg_leith
156  integer :: i, j, k, is, ie, js, je, nz
157  real :: hu(szi_(g), szj_(g)) ! u-thickness [H ~> m or kg m-2]
158  real :: hv(szi_(g), szj_(g)) ! v-thickness [H ~> m or kg m-2]
159  real :: kh_u_lay(szi_(g), szj_(g)) ! layer ave thickness diffusivities [L2 T-1 ~> m2 s-1]
160  real :: kh_v_lay(szi_(g), szj_(g)) ! layer ave thickness diffusivities [L2 T-1 ~> m2 s-1]
161 
162  if (.not. associated(cs)) call mom_error(fatal, "MOM_thickness_diffuse: "//&
163  "Module must be initialized before it is used.")
164 
165  if ((.not.cs%thickness_diffuse) .or. &
166  .not.( cs%Khth > 0.0 .or. associated(varmix) .or. associated(meke) ) ) return
167 
168  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
169  h_neglect = gv%H_subroundoff
170 
171  if (associated(meke)) then
172  if (associated(meke%GM_src)) then
173  do j=js,je ; do i=is,ie ; meke%GM_src(i,j) = 0. ; enddo ; enddo
174  endif
175  endif
176 
177  use_varmix = .false. ; resoln_scaled = .false. ; use_stored_slopes = .false.
178  khth_use_ebt_struct = .false. ; use_visbeck = .false. ; use_qg_leith = .false.
179  depth_scaled = .false.
180 
181  if (associated(varmix)) then
182  use_varmix = varmix%use_variable_mixing .and. (cs%KHTH_Slope_Cff > 0.)
183  resoln_scaled = varmix%Resoln_scaled_KhTh
184  depth_scaled = varmix%Depth_scaled_KhTh
185  use_stored_slopes = varmix%use_stored_slopes
186  khth_use_ebt_struct = varmix%khth_use_ebt_struct
187  use_visbeck = varmix%use_Visbeck
188  use_qg_leith = varmix%use_QG_Leith_GM
189  if (associated(varmix%cg1)) cg1 => varmix%cg1
190  else
191  cg1 => null()
192  endif
193 
194 
195 !$OMP parallel do default(none) shared(is,ie,js,je,KH_u_CFL,dt,G,CS)
196  do j=js,je ; do i=is-1,ie
197  kh_u_cfl(i,j) = (0.25*cs%max_Khth_CFL) / &
198  (dt * (g%IdxCu(i,j)*g%IdxCu(i,j) + g%IdyCu(i,j)*g%IdyCu(i,j)))
199  enddo ; enddo
200 !$OMP parallel do default(none) shared(is,ie,js,je,KH_v_CFL,dt,G,CS)
201  do j=js-1,je ; do i=is,ie
202  kh_v_cfl(i,j) = (0.25*cs%max_Khth_CFL) / &
203  (dt * (g%IdxCv(i,j)*g%IdxCv(i,j) + g%IdyCv(i,j)*g%IdyCv(i,j)))
204  enddo ; enddo
205 
206  ! Calculates interface heights, e, in [Z ~> m].
207  call find_eta(h, tv, g, gv, us, e, halo_size=1)
208 
209  ! Set the diffusivities.
210 !$OMP parallel default(none) shared(is,ie,js,je,Khth_Loc_u,CS,use_VarMix,VarMix, &
211 !$OMP MEKE,Resoln_scaled,KH_u,G,use_QG_Leith,use_Visbeck,&
212 !$OMP KH_u_CFL,nz,Khth_Loc,KH_v,KH_v_CFL,int_slope_u, &
213 !$OMP int_slope_v,khth_use_ebt_struct, Depth_scaled, &
214 !$OMP Khth_loc_v)
215 !$OMP do
216  do j=js,je; do i=is-1,ie
217  khth_loc_u(i,j) = cs%Khth
218  enddo ; enddo
219 
220  if (use_varmix) then
221  if (use_visbeck) then
222 !$OMP do
223  do j=js,je ; do i=is-1,ie
224  khth_loc_u(i,j) = khth_loc_u(i,j) + &
225  cs%KHTH_Slope_Cff*varmix%L2u(i,j) * varmix%SN_u(i,j)
226  enddo ; enddo
227  endif
228  endif
229 
230  if (associated(meke)) then ; if (associated(meke%Kh)) then
231  if (cs%MEKE_GEOMETRIC) then
232 !$OMP do
233  do j=js,je ; do i=is-1,ie
234  khth_loc_u(i,j) = khth_loc_u(i,j) + g%mask2dCu(i,j) * cs%MEKE_GEOMETRIC_alpha * &
235  0.5*(meke%MEKE(i,j)+meke%MEKE(i+1,j)) / &
236  (varmix%SN_u(i,j) + cs%MEKE_GEOMETRIC_epsilon)
237  enddo ; enddo
238  else
239  do j=js,je ; do i=is-1,ie
240  khth_loc_u(i,j) = khth_loc_u(i,j) + meke%KhTh_fac*sqrt(meke%Kh(i,j)*meke%Kh(i+1,j))
241  enddo ; enddo
242  endif
243  endif ; endif
244 
245  if (resoln_scaled) then
246 !$OMP do
247  do j=js,je; do i=is-1,ie
248  khth_loc_u(i,j) = khth_loc_u(i,j) * varmix%Res_fn_u(i,j)
249  enddo ; enddo
250  endif
251 
252  if (depth_scaled) then
253 !$OMP do
254  do j=js,je; do i=is-1,ie
255  khth_loc_u(i,j) = khth_loc_u(i,j) * varmix%Depth_fn_u(i,j)
256  enddo ; enddo
257  endif
258 
259  if (cs%Khth_Max > 0) then
260 !$OMP do
261  do j=js,je; do i=is-1,ie
262  khth_loc_u(i,j) = max(cs%Khth_Min, min(khth_loc_u(i,j), cs%Khth_Max))
263  enddo ; enddo
264  else
265 !$OMP do
266  do j=js,je; do i=is-1,ie
267  khth_loc_u(i,j) = max(cs%Khth_Min, khth_loc_u(i,j))
268  enddo ; enddo
269  endif
270 !$OMP do
271  do j=js,je; do i=is-1,ie
272  kh_u(i,j,1) = min(kh_u_cfl(i,j), khth_loc_u(i,j))
273  enddo ; enddo
274 
275  if (khth_use_ebt_struct) then
276 !$OMP do
277  do k=2,nz+1 ; do j=js,je ; do i=is-1,ie
278  kh_u(i,j,k) = kh_u(i,j,1) * 0.5 * ( varmix%ebt_struct(i,j,k-1) + varmix%ebt_struct(i+1,j,k-1) )
279  enddo ; enddo ; enddo
280  else
281 !$OMP do
282  do k=2,nz+1 ; do j=js,je ; do i=is-1,ie
283  kh_u(i,j,k) = kh_u(i,j,1)
284  enddo ; enddo ; enddo
285  endif
286 
287  if (use_varmix) then
288  if (use_qg_leith) then
289 !$OMP do
290  do k=1,nz ; do j=js,je ; do i=is-1,ie
291  kh_u(i,j,k) = varmix%KH_u_QG(i,j,k)
292  enddo ; enddo ; enddo
293  endif
294  endif
295 
296  if (cs%use_GME_thickness_diffuse) then
297 !$OMP do
298  do k=1,nz+1 ; do j=js,je ; do i=is-1,ie
299  cs%KH_u_GME(i,j,k) = kh_u(i,j,k)
300  enddo ; enddo ; enddo
301  endif
302 
303 !$OMP do
304  do j=js-1,je ; do i=is,ie
305  khth_loc_v(i,j) = cs%Khth
306  enddo ; enddo
307 
308  if (use_varmix) then
309  if (use_visbeck) then
310 !$OMP do
311  do j=js-1,je ; do i=is,ie
312  khth_loc_v(i,j) = khth_loc_v(i,j) + cs%KHTH_Slope_Cff*varmix%L2v(i,j)*varmix%SN_v(i,j)
313  enddo ; enddo
314  endif
315  endif
316  if (associated(meke)) then ; if (associated(meke%Kh)) then
317  if (cs%MEKE_GEOMETRIC) then
318 !$OMP do
319  do j=js-1,je ; do i=is,ie
320  khth_loc_v(i,j) = khth_loc_v(i,j) + g%mask2dCv(i,j) * cs%MEKE_GEOMETRIC_alpha * &
321  0.5*(meke%MEKE(i,j)+meke%MEKE(i,j+1)) / &
322  (varmix%SN_v(i,j) + cs%MEKE_GEOMETRIC_epsilon)
323  enddo ; enddo
324  else
325  do j=js-1,je ; do i=is,ie
326  khth_loc_v(i,j) = khth_loc_v(i,j) + meke%KhTh_fac*sqrt(meke%Kh(i,j)*meke%Kh(i,j+1))
327  enddo ; enddo
328  endif
329  endif ; endif
330 
331  if (resoln_scaled) then
332 !$OMP do
333  do j=js-1,je; do i=is,ie
334  khth_loc_v(i,j) = khth_loc_v(i,j) * varmix%Res_fn_v(i,j)
335  enddo ; enddo
336  endif
337 
338  if (depth_scaled) then
339 !$OMP do
340  do j=js-1,je ; do i=is,ie
341  khth_loc_v(i,j) = khth_loc_v(i,j) * varmix%Depth_fn_v(i,j)
342  enddo ; enddo
343  endif
344 
345  if (cs%Khth_Max > 0) then
346 !$OMP do
347  do j=js-1,je ; do i=is,ie
348  khth_loc_v(i,j) = max(cs%Khth_Min, min(khth_loc_v(i,j), cs%Khth_Max))
349  enddo ; enddo
350  else
351 !$OMP do
352  do j=js-1,je ; do i=is,ie
353  khth_loc_v(i,j) = max(cs%Khth_Min, khth_loc_v(i,j))
354  enddo ; enddo
355  endif
356 
357  if (cs%max_Khth_CFL > 0.0) then
358 !$OMP do
359  do j=js-1,je ; do i=is,ie
360  kh_v(i,j,1) = min(kh_v_cfl(i,j), khth_loc_v(i,j))
361  enddo ; enddo
362  endif
363 
364  if (khth_use_ebt_struct) then
365 !$OMP do
366  do k=2,nz+1 ; do j=js-1,je ; do i=is,ie
367  kh_v(i,j,k) = kh_v(i,j,1) * 0.5 * ( varmix%ebt_struct(i,j,k-1) + varmix%ebt_struct(i,j+1,k-1) )
368  enddo ; enddo ; enddo
369  else
370 !$OMP do
371  do k=2,nz+1 ; do j=js-1,je ; do i=is,ie
372  kh_v(i,j,k) = kh_v(i,j,1)
373  enddo ; enddo ; enddo
374  endif
375 
376  if (use_varmix) then
377  if (use_qg_leith) then
378 !$OMP do
379  do k=1,nz ; do j=js-1,je ; do i=is,ie
380  kh_v(i,j,k) = varmix%KH_v_QG(i,j,k)
381  enddo ; enddo ; enddo
382  endif
383  endif
384 
385  if (cs%use_GME_thickness_diffuse) then
386 !$OMP do
387  do k=1,nz+1 ; do j=js-1,je ; do i=is,ie
388  cs%KH_v_GME(i,j,k) = kh_v(i,j,k)
389  enddo ; enddo ; enddo
390  endif
391 
392  if (associated(meke)) then ; if (associated(meke%Kh)) then
393  if (cs%MEKE_GEOMETRIC) then
394  if (cs%MEKE_GEOM_answers_2018) then
395  !$OMP do
396  do j=js,je ; do i=is,ie
397  ! This does not give bitwise rotational symmetry.
398  meke%Kh(i,j) = cs%MEKE_GEOMETRIC_alpha * meke%MEKE(i,j) / &
399  (0.25*(varmix%SN_u(i,j)+varmix%SN_u(i-1,j) + &
400  varmix%SN_v(i,j)+varmix%SN_v(i,j-1)) + &
401  cs%MEKE_GEOMETRIC_epsilon)
402  enddo ; enddo
403  else
404  !$OMP do
405  do j=js,je ; do i=is,ie
406  ! With the additional parentheses this gives bitwise rotational symmetry.
407  meke%Kh(i,j) = cs%MEKE_GEOMETRIC_alpha * meke%MEKE(i,j) / &
408  (0.25*((varmix%SN_u(i,j)+varmix%SN_u(i-1,j)) + &
409  (varmix%SN_v(i,j)+varmix%SN_v(i,j-1))) + &
410  cs%MEKE_GEOMETRIC_epsilon)
411  enddo ; enddo
412  endif
413  endif
414  endif ; endif
415 
416 
417 !$OMP do
418  do k=1,nz+1 ; do j=js,je ; do i=is-1,ie ; int_slope_u(i,j,k) = 0.0 ; enddo ; enddo ; enddo
419 !$OMP do
420  do k=1,nz+1 ; do j=js-1,je ; do i=is,ie ; int_slope_v(i,j,k) = 0.0 ; enddo ; enddo ; enddo
421 !$OMP end parallel
422 
423  if (cs%detangle_interfaces) then
424  call add_detangling_kh(h, e, kh_u, kh_v, kh_u_cfl, kh_v_cfl, tv, dt, g, gv, us, &
425  cs, int_slope_u, int_slope_v)
426  endif
427 
428  if (cs%debug) then
429  call uvchksum("Kh_[uv]", kh_u, kh_v, g%HI, haloshift=0, &
430  scale=(us%L_to_m**2)*us%s_to_T, scalar_pair=.true.)
431  call uvchksum("int_slope_[uv]", int_slope_u, int_slope_v, g%HI, haloshift=0)
432  call hchksum(h, "thickness_diffuse_1 h", g%HI, haloshift=1, scale=gv%H_to_m)
433  call hchksum(e, "thickness_diffuse_1 e", g%HI, haloshift=1, scale=us%Z_to_m)
434  if (use_stored_slopes) then
435  call uvchksum("VarMix%slope_[xy]", varmix%slope_x, varmix%slope_y, &
436  g%HI, haloshift=0)
437  endif
438  if (associated(tv%eqn_of_state)) then
439  call hchksum(tv%T, "thickness_diffuse T", g%HI, haloshift=1)
440  call hchksum(tv%S, "thickness_diffuse S", g%HI, haloshift=1)
441  endif
442  endif
443 
444  ! Calculate uhD, vhD from h, e, KH_u, KH_v, tv%T/S
445  if (use_stored_slopes) then
446  call thickness_diffuse_full(h, e, kh_u, kh_v, tv, uhd, vhd, cg1, dt, g, gv, us, meke, cs, &
447  int_slope_u, int_slope_v, varmix%slope_x, varmix%slope_y)
448  else
449  call thickness_diffuse_full(h, e, kh_u, kh_v, tv, uhd, vhd, cg1, dt, g, gv, us, meke, cs, &
450  int_slope_u, int_slope_v)
451  endif
452 
453  if (associated(meke) .AND. associated(varmix)) then
454  if (associated(meke%Rd_dx_h) .and. associated(varmix%Rd_dx_h)) then
455 !$OMP parallel do default(none) shared(is,ie,js,je,MEKE,VarMix)
456  do j=js,je ; do i=is,ie
457  meke%Rd_dx_h(i,j) = varmix%Rd_dx_h(i,j)
458  enddo ; enddo
459  endif
460  endif
461 
462  ! offer diagnostic fields for averaging
463  if (query_averaging_enabled(cs%diag)) then
464  if (cs%id_uhGM > 0) call post_data(cs%id_uhGM, uhd, cs%diag)
465  if (cs%id_vhGM > 0) call post_data(cs%id_vhGM, vhd, cs%diag)
466  if (cs%id_GMwork > 0) call post_data(cs%id_GMwork, cs%GMwork, cs%diag)
467  if (cs%id_KH_u > 0) call post_data(cs%id_KH_u, kh_u, cs%diag)
468  if (cs%id_KH_v > 0) call post_data(cs%id_KH_v, kh_v, cs%diag)
469  if (cs%id_KH_u1 > 0) call post_data(cs%id_KH_u1, kh_u(:,:,1), cs%diag)
470  if (cs%id_KH_v1 > 0) call post_data(cs%id_KH_v1, kh_v(:,:,1), cs%diag)
471 
472  ! Diagnose diffusivity at T-cell point. Do simple average, rather than
473  ! thickness-weighted average, in order that KH_t is depth-independent
474  ! in the case where KH_u and KH_v are depth independent. Otherwise,
475  ! if use thickness weighted average, the variations of thickness with
476  ! depth will place a spurious depth dependence to the diagnosed KH_t.
477  if (cs%id_KH_t > 0 .or. cs%id_KH_t1 > 0 .or. cs%Use_KH_in_MEKE) then
478  do k=1,nz
479  ! thicknesses across u and v faces, converted to 0/1 mask
480  ! layer average of the interface diffusivities KH_u and KH_v
481  do j=js,je ; do i=is-1,ie
482  hu(i,j) = 2.0*h(i,j,k)*h(i+1,j,k)/(h(i,j,k)+h(i+1,j,k)+h_neglect)
483  if (hu(i,j) /= 0.0) hu(i,j) = 1.0
484  kh_u_lay(i,j) = 0.5*(kh_u(i,j,k)+kh_u(i,j,k+1))
485  enddo ; enddo
486  do j=js-1,je ; do i=is,ie
487  hv(i,j) = 2.0*h(i,j,k)*h(i,j+1,k)/(h(i,j,k)+h(i,j+1,k)+h_neglect)
488  if (hv(i,j) /= 0.0) hv(i,j) = 1.0
489  kh_v_lay(i,j) = 0.5*(kh_v(i,j,k)+kh_v(i,j,k+1))
490  enddo ; enddo
491  ! diagnose diffusivity at T-point
492  do j=js,je ; do i=is,ie
493  kh_t(i,j,k) = ((hu(i-1,j)*kh_u_lay(i-1,j)+hu(i,j)*kh_u_lay(i,j)) &
494  +(hv(i,j-1)*kh_v_lay(i,j-1)+hv(i,j)*kh_v_lay(i,j))) &
495  / (hu(i-1,j)+hu(i,j)+hv(i,j-1)+hv(i,j)+h_neglect)
496  enddo ; enddo
497  enddo
498 
499  if (cs%Use_KH_in_MEKE) then
500  meke%Kh_diff(:,:) = 0.0
501  do k=1,nz
502  do j=js,je ; do i=is,ie
503  meke%Kh_diff(i,j) = meke%Kh_diff(i,j) + kh_t(i,j,k) * h(i,j,k)
504  enddo; enddo
505  enddo
506 
507  do j=js,je ; do i=is,ie
508  meke%Kh_diff(i,j) = meke%Kh_diff(i,j) / max(1.0,g%bathyT(i,j))
509  enddo ; enddo
510  endif
511 
512  if (cs%id_KH_t > 0) call post_data(cs%id_KH_t, kh_t, cs%diag)
513  if (cs%id_KH_t1 > 0) call post_data(cs%id_KH_t1, kh_t(:,:,1), cs%diag)
514  endif
515 
516  endif
517 
518  !$OMP parallel do default(none) shared(is,ie,js,je,nz,uhtr,uhD,dt,vhtr,CDp,vhD,h,G,GV)
519  do k=1,nz
520  do j=js,je ; do i=is-1,ie
521  uhtr(i,j,k) = uhtr(i,j,k) + uhd(i,j,k) * dt
522  if (associated(cdp%uhGM)) cdp%uhGM(i,j,k) = uhd(i,j,k)
523  enddo ; enddo
524  do j=js-1,je ; do i=is,ie
525  vhtr(i,j,k) = vhtr(i,j,k) + vhd(i,j,k) * dt
526  if (associated(cdp%vhGM)) cdp%vhGM(i,j,k) = vhd(i,j,k)
527  enddo ; enddo
528  do j=js,je ; do i=is,ie
529  h(i,j,k) = h(i,j,k) - dt * g%IareaT(i,j) * &
530  ((uhd(i,j,k) - uhd(i-1,j,k)) + (vhd(i,j,k) - vhd(i,j-1,k)))
531  if (h(i,j,k) < gv%Angstrom_H) h(i,j,k) = gv%Angstrom_H
532  enddo ; enddo
533  enddo
534 
535  ! Whenever thickness changes let the diag manager know, target grids
536  ! for vertical remapping may need to be regenerated.
537  ! This needs to happen after the H update and before the next post_data.
538  call diag_update_remap_grids(cs%diag)
539 
540  if (cs%debug) then
541  call uvchksum("thickness_diffuse [uv]hD", uhd, vhd, &
542  g%HI, haloshift=0, scale=gv%H_to_m*us%L_to_m**2*us%s_to_T)
543  call uvchksum("thickness_diffuse [uv]htr", uhtr, vhtr, &
544  g%HI, haloshift=0, scale=us%L_to_m**2*gv%H_to_m)
545  call hchksum(h, "thickness_diffuse h", g%HI, haloshift=0, scale=gv%H_to_m)
546  endif
547 
548 end subroutine thickness_diffuse
549 
550 !> Calculates parameterized layer transports for use in the continuity equation.
551 !! Fluxes are limited to give positive definite thicknesses.
552 !! Called by thickness_diffuse().
553 subroutine thickness_diffuse_full(h, e, Kh_u, Kh_v, tv, uhD, vhD, cg1, dt, G, GV, US, MEKE, &
554  CS, int_slope_u, int_slope_v, slope_x, slope_y)
555  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
556  type(verticalgrid_type), intent(in) :: GV !< Vertical grid structure
557  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
558  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
559  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: e !< Interface positions [Z ~> m]
560  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(in) :: Kh_u !< Thickness diffusivity on interfaces
561  !! at u points [L2 T-1 ~> m2 s-1]
562  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(in) :: Kh_v !< Thickness diffusivity on interfaces
563  !! at v points [L2 T-1 ~> m2 s-1]
564  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
565  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(out) :: uhD !< Zonal mass fluxes
566  !! [H L2 T-1 ~> m3 s-1 or kg s-1]
567  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(out) :: vhD !< Meridional mass fluxes
568  !! [H L2 T-1 ~> m3 s-1 or kg s-1]
569  real, dimension(:,:), pointer :: cg1 !< Wave speed [L T-1 ~> m s-1]
570  real, intent(in) :: dt !< Time increment [T ~> s]
571  type(meke_type), pointer :: MEKE !< MEKE control structure
572  type(thickness_diffuse_cs), pointer :: CS !< Control structure for thickness diffusion
573  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), optional, intent(in) :: int_slope_u !< Ratio that determine how much of
574  !! the isopycnal slopes are taken directly from
575  !! the interface slopes without consideration of
576  !! density gradients [nondim].
577  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), optional, intent(in) :: int_slope_v !< Ratio that determine how much of
578  !! the isopycnal slopes are taken directly from
579  !! the interface slopes without consideration of
580  !! density gradients [nondim].
581  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), optional, intent(in) :: slope_x !< Isopycnal slope at u-points
582  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), optional, intent(in) :: slope_y !< Isopycnal slope at v-points
583  ! Local variables
584  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: &
585  T, & ! The temperature (or density) [degC], with the values in
586  ! in massless layers filled vertically by diffusion.
587  s, & ! The filled salinity [ppt], with the values in
588  ! in massless layers filled vertically by diffusion.
589  h_avail, & ! The mass available for diffusion out of each face, divided
590  ! by dt [H L2 T-1 ~> m3 s-1 or kg s-1].
591  h_frac ! The fraction of the mass in the column above the bottom
592  ! interface of a layer that is within a layer [nondim]. 0<h_frac<=1
593  real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: &
594  Slope_y_PE, & ! 3D array of neutral slopes at v-points, set equal to Slope (below, nondim)
595  hN2_y_PE ! thickness in m times Brunt-Vaisala freqeuncy at v-points [L2 Z-1 T-2 ~> m s-2],
596  ! used for calculating PE release
597  real, dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: &
598  Slope_x_PE, & ! 3D array of neutral slopes at u-points, set equal to Slope (below, nondim)
599  hN2_x_PE ! thickness in m times Brunt-Vaisala freqeuncy at u-points [L2 Z-1 T-2 ~> m s-2],
600  ! used for calculating PE release
601  real, dimension(SZI_(G), SZJ_(G), SZK_(G)+1) :: &
602  pres, & ! The pressure at an interface [R L2 T-2 ~> Pa].
603  h_avail_rsum ! The running sum of h_avail above an interface [H L2 T-1 ~> m3 s-1 or kg s-1].
604  real, dimension(SZIB_(G)) :: &
605  drho_dT_u, & ! The derivative of density with temperature at u points [R degC-1 ~> kg m-3 degC-1]
606  drho_dS_u, & ! The derivative of density with salinity at u points [R ppt-1 ~> kg m-3 ppt-1].
607  drho_dT_dT_u ! The second derivative of density with temperature at u points [R degC-2 ~> kg m-3 degC-2]
608  real, dimension(SZIB_(G)) :: scrap ! An array to pass to calculate_density_second_derivs() that will be ingored.
609  real, dimension(SZI_(G)) :: &
610  drho_dT_v, & ! The derivative of density with temperature at v points [R degC-1 ~> kg m-3 degC-1]
611  drho_dS_v, & ! The derivative of density with salinity at v points [R ppt-1 ~> kg m-3 ppt-1].
612  drho_dT_dT_v ! The second derivative of density with temperature at v points [R degC-2 ~> kg m-3 degC-2]
613  real :: uhtot(SZIB_(G), SZJ_(G)) ! The vertical sum of uhD [H L2 T-1 ~> m3 s-1 or kg s-1].
614  real :: vhtot(SZI_(G), SZJB_(G)) ! The vertical sum of vhD [H L2 T-1 ~> m3 s-1 or kg s-1].
615  real, dimension(SZIB_(G)) :: &
616  T_u, & ! Temperature on the interface at the u-point [degC].
617  S_u, & ! Salinity on the interface at the u-point [ppt].
618  pres_u ! Pressure on the interface at the u-point [R L2 T-2 ~> Pa].
619  real, dimension(SZI_(G)) :: &
620  T_v, & ! Temperature on the interface at the v-point [degC].
621  S_v, & ! Salinity on the interface at the v-point [ppt].
622  pres_v ! Pressure on the interface at the v-point [R L2 T-2 ~> Pa].
623  real :: Work_u(SZIB_(G), SZJ_(G)) ! The work being done by the thickness
624  real :: Work_v(SZI_(G), SZJB_(G)) ! diffusion integrated over a cell [R Z L4 T-3 ~> W ]
625  real :: Work_h ! The work averaged over an h-cell [R Z L2 T-3 ~> W m-2].
626  real :: PE_release_h ! The amount of potential energy released by GM averaged over an h-cell [L4 Z-1 T-3 ~> m3 s-3]
627  ! The calculation is equal to h * S^2 * N^2 * kappa_GM.
628  real :: I4dt ! 1 / 4 dt [T-1 ~> s-1].
629  real :: drdiA, drdiB ! Along layer zonal- and meridional- potential density
630  real :: drdjA, drdjB ! gradients in the layers above (A) and below(B) the
631  ! interface times the grid spacing [R ~> kg m-3].
632  real :: drdkL, drdkR ! Vertical density differences across an interface [R ~> kg m-3].
633  real :: drdi_u(SZIB_(G), SZK_(G)) ! Copy of drdi at u-points [R ~> kg m-3].
634  real :: drdj_v(SZI_(G), SZK_(G)) ! Copy of drdj at v-points [R ~> kg m-3].
635  real :: drdkDe_u(SZIB_(G),SZK_(G)+1) ! Lateral difference of product of drdk and e at u-points
636  ! [Z R ~> kg m-2].
637  real :: drdkDe_v(SZI_(G),SZK_(G)+1) ! Lateral difference of product of drdk and e at v-points
638  ! [Z R ~> kg m-2].
639  real :: hg2A, hg2B, hg2L, hg2R ! Squares of geometric mean thicknesses [H2 ~> m2 or kg2 m-4].
640  real :: haA, haB, haL, haR ! Arithmetic mean thicknesses [H ~> m or kg m-2].
641  real :: dzaL, dzaR ! Temporary thicknesses [Z ~> m].
642  real :: wtA, wtB, wtL, wtR ! Unscaled weights, with various units.
643  real :: drdx, drdy ! Zonal and meridional density gradients [R L-1 ~> kg m-4].
644  real :: drdz ! Vertical density gradient [R Z-1 ~> kg m-4].
645  real :: h_harm ! Harmonic mean layer thickness [H ~> m or kg m-2].
646  real :: c2_h_u(SZIB_(G), SZK_(G)+1) ! Wave speed squared divided by h at u-points [L2 Z-1 T-2 ~> m s-2].
647  real :: c2_h_v(SZI_(G), SZK_(G)+1) ! Wave speed squared divided by h at v-points [L2 Z-1 T-2 ~> m s-2].
648  real :: hN2_u(SZIB_(G), SZK_(G)+1) ! Thickness in m times N2 at interfaces above u-points [L2 Z-1 T-2 ~> m s-2].
649  real :: hN2_v(SZI_(G), SZK_(G)+1) ! Thickness in m times N2 at interfaces above v-points [L2 Z-1 T-2 ~> m s-2].
650  real :: Sfn_est ! A preliminary estimate (before limiting) of the overturning
651  ! streamfunction [Z L2 T-1 ~> m3 s-1].
652  real :: Sfn_unlim_u(SZIB_(G), SZK_(G)+1) ! Streamfunction for u-points [Z L2 T-1 ~> m3 s-1].
653  real :: Sfn_unlim_v(SZI_(G), SZK_(G)+1) ! Streamfunction for v-points [Z L2 T-1 ~> m3 s-1].
654  real :: slope2_Ratio_u(SZIB_(G), SZK_(G)+1) ! The ratio of the slope squared to slope_max squared.
655  real :: slope2_Ratio_v(SZI_(G), SZK_(G)+1) ! The ratio of the slope squared to slope_max squared.
656  real :: Sfn_in_h ! The overturning streamfunction [H L2 T-1 ~> m3 s-1 or kg s-1] (note that
657  ! the units are different from other Sfn vars).
658  real :: Sfn_safe ! The streamfunction that goes linearly back to 0 at the surface. This is a
659  ! good thing to use when the slope is so large as to be meaningless [Z L2 T-1 ~> m3 s-1].
660  real :: Slope ! The slope of density surfaces, calculated in a way
661  ! that is always between -1 and 1, nondimensional.
662  real :: mag_grad2 ! The squared magnitude of the 3-d density gradient [R2 L-2 ~> kg2 m-8].
663  real :: I_slope_max2 ! The inverse of slope_max squared, nondimensional.
664  real :: h_neglect ! A thickness that is so small it is usually lost
665  ! in roundoff and can be neglected [H ~> m or kg m-2].
666  real :: h_neglect2 ! h_neglect^2 [H2 ~> m2 or kg2 m-4].
667  real :: dz_neglect ! A thickness [Z ~> m], that is so small it is usually lost
668  ! in roundoff and can be neglected [Z ~> m].
669  real :: G_scale ! The gravitational acceleration times a unit conversion
670  ! factor [L2 H-1 T-2 ~> m s-2 or m4 kg-1 s-2].
671  logical :: use_EOS ! If true, density is calculated from T & S using an
672  ! equation of state.
673  logical :: find_work ! If true, find the change in energy due to the fluxes.
674  integer :: nk_linear ! The number of layers over which the streamfunction goes to 0.
675  real :: G_rho0 ! g/Rho0 [L2 R-1 Z-1 T-2 ~> m4 kg-1 s-2].
676  real :: N2_floor ! A floor for N2 to avoid degeneracy in the elliptic solver
677  ! times unit conversion factors [T-2 L2 Z-2 ~> s-2]
678  real :: Tl(5) ! copy and T in local stencil [degC]
679  real :: mn_T ! mean of T in local stencil [degC]
680  real :: mn_T2 ! mean of T**2 in local stencil [degC]
681  real :: hl(5) ! Copy of local stencil of H [H ~> m]
682  real :: r_sm_H ! Reciprocal of sum of H in local stencil [H-1 ~> m-1]
683  real, dimension(SZI_(G), SZJ_(G), SZK_(G)) :: Tsgs2 ! Sub-grid temperature variance [degC2]
684 
685  real, dimension(SZIB_(G), SZJ_(G), SZK_(G)+1) :: diag_sfn_x, diag_sfn_unlim_x ! Diagnostics
686  real, dimension(SZI_(G), SZJB_(G), SZK_(G)+1) :: diag_sfn_y, diag_sfn_unlim_y ! Diagnostics
687  logical :: present_int_slope_u, present_int_slope_v
688  logical :: present_slope_x, present_slope_y, calc_derivatives
689  integer, dimension(2) :: EOSdom_u ! The shifted i-computational domain to use for equation of
690  ! state calculations at u-points.
691  integer, dimension(2) :: EOSdom_v ! The shifted I-computational domain to use for equation of
692  ! state calculations at v-points.
693  logical :: use_Stanley
694  integer :: is, ie, js, je, nz, IsdB, halo
695  integer :: i, j, k
696  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke ; isdb = g%IsdB
697 
698  i4dt = 0.25 / dt
699  i_slope_max2 = 1.0 / (cs%slope_max**2)
700  g_scale = gv%g_Earth * gv%H_to_Z
701 
702  h_neglect = gv%H_subroundoff ; h_neglect2 = h_neglect**2
703  dz_neglect = gv%H_subroundoff*gv%H_to_Z
704  g_rho0 = gv%g_Earth / gv%Rho0
705  n2_floor = cs%N2_floor*us%Z_to_L**2
706 
707  use_eos = associated(tv%eqn_of_state)
708  present_int_slope_u = PRESENT(int_slope_u)
709  present_int_slope_v = PRESENT(int_slope_v)
710  present_slope_x = PRESENT(slope_x)
711  present_slope_y = PRESENT(slope_y)
712  use_stanley = cs%Stanley_det_coeff >= 0.
713 
714  nk_linear = max(gv%nkml, 1)
715 
716  slope_x_pe(:,:,:) = 0.0
717  slope_y_pe(:,:,:) = 0.0
718  hn2_x_pe(:,:,:) = 0.0
719  hn2_y_pe(:,:,:) = 0.0
720 
721  find_work = .false.
722  if (associated(meke)) find_work = associated(meke%GM_src)
723  find_work = (associated(cs%GMwork) .or. find_work)
724 
725  if (use_eos) then
726  halo = 1 ! Default halo to fill is 1
727  if (use_stanley) halo = 2 ! Need wider valid halo for gradients of T
728  call vert_fill_ts(h, tv%T, tv%S, cs%kappa_smooth*dt, t, s, g, gv, halo, larger_h_denom=.true.)
729  endif
730 
731  if (cs%use_FGNV_streamfn .and. .not. associated(cg1)) call mom_error(fatal, &
732  "cg1 must be associated when using FGNV streamfunction.")
733 
734 !$OMP parallel default(none) shared(is,ie,js,je,h_avail_rsum,pres,h_avail,I4dt, use_Stanley, &
735 !$OMP CS,G,GV,tv,h,h_frac,nz,uhtot,Work_u,vhtot,Work_v,Tsgs2,T, &
736 !$OMP diag_sfn_x, diag_sfn_y, diag_sfn_unlim_x, diag_sfn_unlim_y ) &
737 !$OMP private(hl,r_sm_H,Tl,mn_T,mn_T2)
738  ! Find the maximum and minimum permitted streamfunction.
739 !$OMP do
740  do j=js-1,je+1 ; do i=is-1,ie+1
741  h_avail_rsum(i,j,1) = 0.0
742  pres(i,j,1) = 0.0
743  if (associated(tv%p_surf)) then ; pres(i,j,1) = tv%p_surf(i,j) ; endif
744 
745  h_avail(i,j,1) = max(i4dt*g%areaT(i,j)*(h(i,j,1)-gv%Angstrom_H),0.0)
746  h_avail_rsum(i,j,2) = h_avail(i,j,1)
747  h_frac(i,j,1) = 1.0
748  pres(i,j,2) = pres(i,j,1) + (gv%g_Earth*gv%H_to_RZ) * h(i,j,1)
749  enddo ; enddo
750  if (use_stanley) then
751 !$OMP do
752  do k=1, nz ; do j=js-1,je+1 ; do i=is-1,ie+1
753  !! SGS variance in i-direction [degC2]
754  !dTdi2 = ( ( G%mask2dCu(I ,j) * G%IdxCu(I ,j) * ( T(i+1,j,k) - T(i,j,k) ) &
755  ! + G%mask2dCu(I-1,j) * G%IdxCu(I-1,j) * ( T(i,j,k) - T(i-1,j,k) ) &
756  ! ) * G%dxT(i,j) * 0.5 )**2
757  !! SGS variance in j-direction [degC2]
758  !dTdj2 = ( ( G%mask2dCv(i,J ) * G%IdyCv(i,J ) * ( T(i,j+1,k) - T(i,j,k) ) &
759  ! + G%mask2dCv(i,J-1) * G%IdyCv(i,J-1) * ( T(i,j,k) - T(i,j-1,k) ) &
760  ! ) * G%dyT(i,j) * 0.5 )**2
761  !Tsgs2(i,j,k) = CS%Stanley_det_coeff * 0.5 * ( dTdi2 + dTdj2 )
762  ! This block does a thickness weighted variance calculation and helps control for
763  ! extreme gradients along layers which are vanished against topography. It is
764  ! still a poor approximation in the interior when coordinates are strongly tilted.
765  hl(1) = h(i,j,k) * g%mask2dT(i,j)
766  hl(2) = h(i-1,j,k) * g%mask2dCu(i-1,j)
767  hl(3) = h(i+1,j,k) * g%mask2dCu(i,j)
768  hl(4) = h(i,j-1,k) * g%mask2dCv(i,j-1)
769  hl(5) = h(i,j+1,k) * g%mask2dCv(i,j)
770  r_sm_h = 1. / ( ( hl(1) + ( ( hl(2) + hl(3) ) + ( hl(4) + hl(5) ) ) ) + gv%H_subroundoff )
771  ! Mean of T
772  tl(1) = t(i,j,k) ; tl(2) = t(i-1,j,k) ; tl(3) = t(i+1,j,k)
773  tl(4) = t(i,j-1,k) ; tl(5) = t(i,j+1,k)
774  mn_t = ( hl(1)*tl(1) + ( ( hl(2)*tl(2) + hl(3)*tl(3) ) + ( hl(4)*tl(4) + hl(5)*tl(5) ) ) ) * r_sm_h
775  ! Adjust T vectors to have zero mean
776  tl(:) = tl(:) - mn_t ; mn_t = 0.
777  ! Variance of T
778  mn_t2 = ( hl(1)*tl(1)*tl(1) + ( ( hl(2)*tl(2)*tl(2) + hl(3)*tl(3)*tl(3) ) &
779  + ( hl(4)*tl(4)*tl(4) + hl(5)*tl(5)*tl(5) ) ) ) * r_sm_h
780  ! Variance should be positive but round-off can violate this. Calculating
781  ! variance directly would fix this but requires more operations.
782  tsgs2(i,j,k) = cs%Stanley_det_coeff * max(0., mn_t2)
783  enddo ; enddo ; enddo
784  endif
785 !$OMP do
786  do j=js-1,je+1
787  do k=2,nz ; do i=is-1,ie+1
788  h_avail(i,j,k) = max(i4dt*g%areaT(i,j)*(h(i,j,k)-gv%Angstrom_H),0.0)
789  h_avail_rsum(i,j,k+1) = h_avail_rsum(i,j,k) + h_avail(i,j,k)
790  h_frac(i,j,k) = 0.0 ; if (h_avail(i,j,k) > 0.0) &
791  h_frac(i,j,k) = h_avail(i,j,k) / h_avail_rsum(i,j,k+1)
792  pres(i,j,k+1) = pres(i,j,k) + (gv%g_Earth*gv%H_to_RZ) * h(i,j,k)
793  enddo ; enddo
794  enddo
795 !$OMP do
796  do j=js,je ; do i=is-1,ie
797  uhtot(i,j) = 0.0 ; work_u(i,j) = 0.0
798  diag_sfn_x(i,j,1) = 0.0 ; diag_sfn_unlim_x(i,j,1) = 0.0
799  diag_sfn_x(i,j,nz+1) = 0.0 ; diag_sfn_unlim_x(i,j,nz+1) = 0.0
800  enddo ; enddo
801 !$OMP do
802  do j=js-1,je ; do i=is,ie
803  vhtot(i,j) = 0.0 ; work_v(i,j) = 0.0
804  diag_sfn_y(i,j,1) = 0.0 ; diag_sfn_unlim_y(i,j,1) = 0.0
805  diag_sfn_y(i,j,nz+1) = 0.0 ; diag_sfn_unlim_y(i,j,nz+1) = 0.0
806  enddo ; enddo
807 !$OMP end parallel
808 
809  eosdom_u(1) = (is-1) - (g%IsdB-1) ; eosdom_u(2) = ie - (g%IsdB-1)
810 !$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S, &
811 !$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect, &
812 !$OMP I_slope_max2,h_neglect2,present_int_slope_u, &
813 !$OMP int_slope_u,KH_u,uhtot,h_frac,h_avail_rsum, &
814 !$OMP uhD,h_avail,G_scale,Work_u,CS,slope_x,cg1, &
815 !$OMP diag_sfn_x, diag_sfn_unlim_x,N2_floor,EOSdom_u, &
816 !$OMP use_stanley, Tsgs2, &
817 !$OMP present_slope_x,G_rho0,Slope_x_PE,hN2_x_PE) &
818 !$OMP private(drdiA,drdiB,drdkL,drdkR,pres_u,T_u,S_u, &
819 !$OMP drho_dT_u,drho_dS_u,hg2A,hg2B,hg2L,hg2R,haA, &
820 !$OMP drho_dT_dT_u,scrap, &
821 !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, &
822 !$OMP drdx,mag_grad2,Slope,slope2_Ratio_u,hN2_u, &
823 !$OMP Sfn_unlim_u,drdi_u,drdkDe_u,h_harm,c2_h_u, &
824 !$OMP Sfn_safe,Sfn_est,Sfn_in_h,calc_derivatives)
825  do j=js,je
826  do i=is-1,ie ; hn2_u(i,1) = 0. ; hn2_u(i,nz+1) = 0. ; enddo
827  do k=nz,2,-1
828  if (find_work .and. .not.(use_eos)) then
829  drdia = 0.0 ; drdib = 0.0
830  drdkl = gv%Rlay(k) - gv%Rlay(k-1) ; drdkr = drdkl
831  endif
832 
833  calc_derivatives = use_eos .and. (k >= nk_linear) .and. &
834  (find_work .or. .not. present_slope_x .or. cs%use_FGNV_streamfn .or. use_stanley)
835 
836  ! Calculate the zonal fluxes and gradients.
837  if (calc_derivatives) then
838  do i=is-1,ie
839  pres_u(i) = 0.5*(pres(i,j,k) + pres(i+1,j,k))
840  t_u(i) = 0.25*((t(i,j,k) + t(i+1,j,k)) + (t(i,j,k-1) + t(i+1,j,k-1)))
841  s_u(i) = 0.25*((s(i,j,k) + s(i+1,j,k)) + (s(i,j,k-1) + s(i+1,j,k-1)))
842  enddo
843  call calculate_density_derivs(t_u, s_u, pres_u, drho_dt_u, drho_ds_u, &
844  tv%eqn_of_state, eosdom_u)
845  endif
846  if (use_stanley) then
847  ! The second line below would correspond to arguments
848  ! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, &
849  call calculate_density_second_derivs(t_u, s_u, pres_u, &
850  scrap, scrap, drho_dt_dt_u, scrap, scrap, &
851  (is-isdb+1)-1, ie-is+2, tv%eqn_of_state)
852  endif
853 
854  do i=is-1,ie
855  if (calc_derivatives) then
856  ! Estimate the horizontal density gradients along layers.
857  drdia = drho_dt_u(i) * (t(i+1,j,k-1)-t(i,j,k-1)) + &
858  drho_ds_u(i) * (s(i+1,j,k-1)-s(i,j,k-1))
859  drdib = drho_dt_u(i) * (t(i+1,j,k)-t(i,j,k)) + &
860  drho_ds_u(i) * (s(i+1,j,k)-s(i,j,k))
861 
862  ! Estimate the vertical density gradients times the grid spacing.
863  drdkl = (drho_dt_u(i) * (t(i,j,k)-t(i,j,k-1)) + &
864  drho_ds_u(i) * (s(i,j,k)-s(i,j,k-1)))
865  drdkr = (drho_dt_u(i) * (t(i+1,j,k)-t(i+1,j,k-1)) + &
866  drho_ds_u(i) * (s(i+1,j,k)-s(i+1,j,k-1)))
867  drdkde_u(i,k) = drdkr * e(i+1,j,k) - drdkl * e(i,j,k)
868  elseif (find_work) then ! This is used in pure stacked SW mode
869  drdkde_u(i,k) = drdkr * e(i+1,j,k) - drdkl * e(i,j,k)
870  endif
871  if (use_stanley) then
872  ! Correction to the horizontal density gradient due to nonlinearity in
873  ! the EOS rectifying SGS temperature anomalies
874  drdia = drdia + drho_dt_dt_u(i) * 0.5 * ( tsgs2(i+1,j,k-1)-tsgs2(i,j,k-1) )
875  drdib = drdib + drho_dt_dt_u(i) * 0.5 * ( tsgs2(i+1,j,k)-tsgs2(i,j,k) )
876  endif
877  if (find_work) drdi_u(i,k) = drdib
878 
879  if (k > nk_linear) then
880  if (use_eos) then
881  if (cs%use_FGNV_streamfn .or. find_work .or. .not.present_slope_x) then
882  hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
883  hg2r = h(i+1,j,k-1)*h(i+1,j,k) + h_neglect2
884  hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
885  har = 0.5*(h(i+1,j,k-1) + h(i+1,j,k)) + h_neglect
886  if (gv%Boussinesq) then
887  dzal = hal * gv%H_to_Z ; dzar = har * gv%H_to_Z
888  else
889  dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
890  dzar = 0.5*(e(i+1,j,k-1) - e(i+1,j,k+1)) + dz_neglect
891  endif
892  ! Use the harmonic mean thicknesses to weight the horizontal gradients.
893  ! These unnormalized weights have been rearranged to minimize divisions.
894  wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)
895 
896  drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
897  ! The expression for drdz above is mathematically equivalent to:
898  ! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
899  ! ((hg2L/haL) + (hg2R/haR))
900  hg2a = h(i,j,k-1)*h(i+1,j,k-1) + h_neglect2
901  hg2b = h(i,j,k)*h(i+1,j,k) + h_neglect2
902  haa = 0.5*(h(i,j,k-1) + h(i+1,j,k-1)) + h_neglect
903  hab = 0.5*(h(i,j,k) + h(i+1,j,k)) + h_neglect
904 
905  ! hN2_u is used with the FGNV streamfunction formulation
906  hn2_u(i,k) = (0.5 * gv%H_to_Z * ( hg2a / haa + hg2b / hab )) * &
907  max(drdz*g_rho0, n2_floor)
908  endif
909  if (present_slope_x) then
910  slope = slope_x(i,j,k)
911  slope2_ratio_u(i,k) = slope**2 * i_slope_max2
912  else
913  ! Use the harmonic mean thicknesses to weight the horizontal gradients.
914  ! These unnormalized weights have been rearranged to minimize divisions.
915  wta = hg2a*hab ; wtb = hg2b*haa
916  ! This is the gradient of density along geopotentials.
917  drdx = ((wta * drdia + wtb * drdib) / (wta + wtb) - &
918  drdz * (e(i,j,k)-e(i+1,j,k))) * g%IdxCu(i,j)
919 
920  ! This estimate of slope is accurate for small slopes, but bounded
921  ! to be between -1 and 1.
922  mag_grad2 = drdx**2 + (us%L_to_Z*drdz)**2
923  if (mag_grad2 > 0.0) then
924  slope = drdx / sqrt(mag_grad2)
925  slope2_ratio_u(i,k) = slope**2 * i_slope_max2
926  else ! Just in case mag_grad2 = 0 ever.
927  slope = 0.0
928  slope2_ratio_u(i,k) = 1.0e20 ! Force the use of the safe streamfunction.
929  endif
930  endif
931 
932  ! Adjust real slope by weights that bias towards slope of interfaces
933  ! that ignore density gradients along layers.
934  if (present_int_slope_u) then
935  slope = (1.0 - int_slope_u(i,j,k)) * slope + &
936  int_slope_u(i,j,k) * us%Z_to_L*((e(i+1,j,k)-e(i,j,k)) * g%IdxCu(i,j))
937  slope2_ratio_u(i,k) = (1.0 - int_slope_u(i,j,k)) * slope2_ratio_u(i,k)
938  endif
939 
940  slope_x_pe(i,j,k) = min(slope,cs%slope_max)
941  hn2_x_pe(i,j,k) = hn2_u(i,k)
942  if (cs%id_slope_x > 0) cs%diagSlopeX(i,j,k) = slope
943 
944  ! Estimate the streamfunction at each interface [Z L2 T-1 ~> m3 s-1].
945  sfn_unlim_u(i,k) = -((kh_u(i,j,k)*g%dy_Cu(i,j))*us%L_to_Z*slope)
946 
947  ! Avoid moving dense water upslope from below the level of
948  ! the bottom on the receiving side.
949  if (sfn_unlim_u(i,k) > 0.0) then ! The flow below this interface is positive.
950  if (e(i,j,k) < e(i+1,j,nz+1)) then
951  sfn_unlim_u(i,k) = 0.0 ! This is not uhtot, because it may compensate for
952  ! deeper flow in very unusual cases.
953  elseif (e(i+1,j,nz+1) > e(i,j,k+1)) then
954  ! Scale the transport with the fraction of the donor layer above
955  ! the bottom on the receiving side.
956  sfn_unlim_u(i,k) = sfn_unlim_u(i,k) * ((e(i,j,k) - e(i+1,j,nz+1)) / &
957  ((e(i,j,k) - e(i,j,k+1)) + dz_neglect))
958  endif
959  else
960  if (e(i+1,j,k) < e(i,j,nz+1)) then ; sfn_unlim_u(i,k) = 0.0
961  elseif (e(i,j,nz+1) > e(i+1,j,k+1)) then
962  sfn_unlim_u(i,k) = sfn_unlim_u(i,k) * ((e(i+1,j,k) - e(i,j,nz+1)) / &
963  ((e(i+1,j,k) - e(i+1,j,k+1)) + dz_neglect))
964  endif
965  endif
966 
967  else ! .not. use_EOS
968  if (present_slope_x) then
969  slope = slope_x(i,j,k)
970  else
971  slope = us%Z_to_L*((e(i,j,k)-e(i+1,j,k))*g%IdxCu(i,j)) * g%mask2dCu(i,j)
972  endif
973  if (cs%id_slope_x > 0) cs%diagSlopeX(i,j,k) = slope
974  sfn_unlim_u(i,k) = ((kh_u(i,j,k)*g%dy_Cu(i,j))*us%L_to_Z*slope)
975  hn2_u(i,k) = gv%g_prime(k)
976  endif ! if (use_EOS)
977  else ! if (k > nk_linear)
978  hn2_u(i,k) = n2_floor * dz_neglect
979  sfn_unlim_u(i,k) = 0.
980  endif ! if (k > nk_linear)
981  if (cs%id_sfn_unlim_x>0) diag_sfn_unlim_x(i,j,k) = sfn_unlim_u(i,k)
982  enddo ! i-loop
983  enddo ! k-loop
984 
985  if (cs%use_FGNV_streamfn) then
986  do k=1,nz ; do i=is-1,ie ; if (g%mask2dCu(i,j)>0.) then
987  h_harm = max( h_neglect, &
988  2. * h(i,j,k) * h(i+1,j,k) / ( ( h(i,j,k) + h(i+1,j,k) ) + h_neglect ) )
989  c2_h_u(i,k) = cs%FGNV_scale * &
990  ( 0.5*( cg1(i,j) + cg1(i+1,j) ) )**2 / (gv%H_to_Z*h_harm)
991  endif ; enddo ; enddo
992 
993  ! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010.
994  do i=is-1,ie
995  if (g%mask2dCu(i,j)>0.) then
996  sfn_unlim_u(i,:) = ( 1. + cs%FGNV_scale ) * sfn_unlim_u(i,:)
997  call streamfn_solver(nz, c2_h_u(i,:), hn2_u(i,:), sfn_unlim_u(i,:))
998  else
999  sfn_unlim_u(i,:) = 0.
1000  endif
1001  enddo
1002  endif
1003 
1004  do k=nz,2,-1
1005  do i=is-1,ie
1006  if (k > nk_linear) then
1007  if (use_eos) then
1008 
1009  if (uhtot(i,j) <= 0.0) then
1010  ! The transport that must balance the transport below is positive.
1011  sfn_safe = uhtot(i,j) * (1.0 - h_frac(i,j,k)) * gv%H_to_Z
1012  else ! (uhtot(I,j) > 0.0)
1013  sfn_safe = uhtot(i,j) * (1.0 - h_frac(i+1,j,k)) * gv%H_to_Z
1014  endif
1015 
1016  ! The actual streamfunction at each interface.
1017  sfn_est = (sfn_unlim_u(i,k) + slope2_ratio_u(i,k)*sfn_safe) / (1.0 + slope2_ratio_u(i,k))
1018  else ! With .not.use_EOS, the layers are constant density.
1019  sfn_est = sfn_unlim_u(i,k)
1020  endif
1021 
1022  ! Make sure that there is enough mass above to allow the streamfunction
1023  ! to satisfy the boundary condition of 0 at the surface.
1024  sfn_in_h = min(max(sfn_est * gv%Z_to_H, -h_avail_rsum(i,j,k)), h_avail_rsum(i+1,j,k))
1025 
1026  ! The actual transport is limited by the mass available in the two
1027  ! neighboring grid cells.
1028  uhd(i,j,k) = max(min((sfn_in_h - uhtot(i,j)), h_avail(i,j,k)), &
1029  -h_avail(i+1,j,k))
1030 
1031  if (cs%id_sfn_x>0) diag_sfn_x(i,j,k) = diag_sfn_x(i,j,k+1) + uhd(i,j,k)
1032 ! sfn_x(I,j,K) = max(min(Sfn_in_h, uhtot(I,j)+h_avail(i,j,k)), &
1033 ! uhtot(I,j)-h_avail(i+1,j,K))
1034 ! sfn_slope_x(I,j,K) = max(uhtot(I,j)-h_avail(i+1,j,k), &
1035 ! min(uhtot(I,j)+h_avail(i,j,k), &
1036 ! min(h_avail_rsum(i+1,j,K), max(-h_avail_rsum(i,j,K), &
1037 ! (KH_u(I,j,K)*G%dy_Cu(I,j)) * ((e(i,j,K)-e(i+1,j,K))*G%IdxCu(I,j)) )) ))
1038  else ! k <= nk_linear
1039  ! Balance the deeper flow with a return flow uniformly distributed
1040  ! though the remaining near-surface layers. This is the same as
1041  ! using Sfn_safe above. There is no need to apply the limiters in
1042  ! this case.
1043  if (uhtot(i,j) <= 0.0) then
1044  uhd(i,j,k) = -uhtot(i,j) * h_frac(i,j,k)
1045  else ! (uhtot(I,j) > 0.0)
1046  uhd(i,j,k) = -uhtot(i,j) * h_frac(i+1,j,k)
1047  endif
1048 
1049  if (cs%id_sfn_x>0) diag_sfn_x(i,j,k) = diag_sfn_x(i,j,k+1) + uhd(i,j,k)
1050 ! if (sfn_slope_x(I,j,K+1) <= 0.0) then
1051 ! sfn_slope_x(I,j,K) = sfn_slope_x(I,j,K+1) * (1.0 - h_frac(i,j,k))
1052 ! else
1053 ! sfn_slope_x(I,j,K) = sfn_slope_x(I,j,K+1) * (1.0 - h_frac(i+1,j,k))
1054 ! endif
1055  endif
1056 
1057  uhtot(i,j) = uhtot(i,j) + uhd(i,j,k)
1058 
1059  if (find_work) then
1060  ! This is the energy tendency based on the original profiles, and does
1061  ! not include any nonlinear terms due to a finite time step (which would
1062  ! involve interactions between the fluxes through the different faces.
1063  ! A second order centered estimate is used for the density transfered
1064  ! between water columns.
1065 
1066  work_u(i,j) = work_u(i,j) + g_scale * &
1067  ( uhtot(i,j) * drdkde_u(i,k) - &
1068  (uhd(i,j,k) * drdi_u(i,k)) * 0.25 * &
1069  ((e(i,j,k) + e(i,j,k+1)) + (e(i+1,j,k) + e(i+1,j,k+1))) )
1070  endif
1071 
1072  enddo
1073  enddo ! end of k-loop
1074  enddo ! end of j-loop
1075 
1076  ! Calculate the meridional fluxes and gradients.
1077  eosdom_v(:) = eos_domain(g%HI)
1078 !$OMP parallel do default(none) shared(nz,is,ie,js,je,find_work,use_EOS,G,GV,US,pres,T,S, &
1079 !$OMP nk_linear,IsdB,tv,h,h_neglect,e,dz_neglect, &
1080 !$OMP I_slope_max2,h_neglect2,present_int_slope_v, &
1081 !$OMP int_slope_v,KH_v,vhtot,h_frac,h_avail_rsum, &
1082 !$OMP vhD,h_avail,G_scale,Work_v,CS,slope_y,cg1, &
1083 !$OMP diag_sfn_y,diag_sfn_unlim_y,N2_floor,EOSdom_v,&
1084 !$OMP use_stanley, Tsgs2, &
1085 !$OMP present_slope_y,G_rho0,Slope_y_PE,hN2_y_PE) &
1086 !$OMP private(drdjA,drdjB,drdkL,drdkR,pres_v,T_v,S_v, &
1087 !$OMP drho_dT_v,drho_dS_v,hg2A,hg2B,hg2L,hg2R,haA, &
1088 !$OMP drho_dT_dT_v,scrap, &
1089 !$OMP haB,haL,haR,dzaL,dzaR,wtA,wtB,wtL,wtR,drdz, &
1090 !$OMP drdy,mag_grad2,Slope,slope2_Ratio_v,hN2_v, &
1091 !$OMP Sfn_unlim_v,drdj_v,drdkDe_v,h_harm,c2_h_v, &
1092 !$OMP Sfn_safe,Sfn_est,Sfn_in_h,calc_derivatives)
1093  do j=js-1,je
1094  do k=nz,2,-1
1095  if (find_work .and. .not.(use_eos)) then
1096  drdja = 0.0 ; drdjb = 0.0
1097  drdkl = gv%Rlay(k) - gv%Rlay(k-1) ; drdkr = drdkl
1098  endif
1099 
1100  calc_derivatives = use_eos .and. (k >= nk_linear) .and. &
1101  (find_work .or. .not. present_slope_y .or. cs%use_FGNV_streamfn .or. use_stanley)
1102 
1103  if (calc_derivatives) then
1104  do i=is,ie
1105  pres_v(i) = 0.5*(pres(i,j,k) + pres(i,j+1,k))
1106  t_v(i) = 0.25*((t(i,j,k) + t(i,j+1,k)) + (t(i,j,k-1) + t(i,j+1,k-1)))
1107  s_v(i) = 0.25*((s(i,j,k) + s(i,j+1,k)) + (s(i,j,k-1) + s(i,j+1,k-1)))
1108  enddo
1109  call calculate_density_derivs(t_v, s_v, pres_v, drho_dt_v, drho_ds_v, &
1110  tv%eqn_of_state, eosdom_v)
1111  endif
1112  if (use_stanley) then
1113  ! The second line below would correspond to arguments
1114  ! drho_dS_dS, drho_dS_dT, drho_dT_dT, drho_dS_dP, drho_dT_dP, &
1115  call calculate_density_second_derivs(t_v, s_v, pres_v, &
1116  scrap, scrap, drho_dt_dt_v, scrap, scrap, &
1117  is, ie-is+1, tv%eqn_of_state)
1118  endif
1119  do i=is,ie
1120  if (calc_derivatives) then
1121  ! Estimate the horizontal density gradients along layers.
1122  drdja = drho_dt_v(i) * (t(i,j+1,k-1)-t(i,j,k-1)) + &
1123  drho_ds_v(i) * (s(i,j+1,k-1)-s(i,j,k-1))
1124  drdjb = drho_dt_v(i) * (t(i,j+1,k)-t(i,j,k)) + &
1125  drho_ds_v(i) * (s(i,j+1,k)-s(i,j,k))
1126 
1127  ! Estimate the vertical density gradients times the grid spacing.
1128  drdkl = (drho_dt_v(i) * (t(i,j,k)-t(i,j,k-1)) + &
1129  drho_ds_v(i) * (s(i,j,k)-s(i,j,k-1)))
1130  drdkr = (drho_dt_v(i) * (t(i,j+1,k)-t(i,j+1,k-1)) + &
1131  drho_ds_v(i) * (s(i,j+1,k)-s(i,j+1,k-1)))
1132  drdkde_v(i,k) = drdkr * e(i,j+1,k) - drdkl * e(i,j,k)
1133  elseif (find_work) then ! This is used in pure stacked SW mode
1134  drdkde_v(i,k) = drdkr * e(i,j+1,k) - drdkl * e(i,j,k)
1135  endif
1136  if (use_stanley) then
1137  ! Correction to the horizontal density gradient due to nonlinearity in
1138  ! the EOS rectifying SGS temperature anomalies
1139  drdja = drdja + drho_dt_dt_v(i) * 0.5 * ( tsgs2(i,j+1,k-1)-tsgs2(i,j,k-1) )
1140  drdjb = drdjb + drho_dt_dt_v(i) * 0.5 * ( tsgs2(i,j+1,k)-tsgs2(i,j,k) )
1141  endif
1142 
1143  if (find_work) drdj_v(i,k) = drdjb
1144 
1145  if (k > nk_linear) then
1146  if (use_eos) then
1147  if (cs%use_FGNV_streamfn .or. find_work .or. .not. present_slope_y) then
1148  hg2l = h(i,j,k-1)*h(i,j,k) + h_neglect2
1149  hg2r = h(i,j+1,k-1)*h(i,j+1,k) + h_neglect2
1150  hal = 0.5*(h(i,j,k-1) + h(i,j,k)) + h_neglect
1151  har = 0.5*(h(i,j+1,k-1) + h(i,j+1,k)) + h_neglect
1152  if (gv%Boussinesq) then
1153  dzal = hal * gv%H_to_Z ; dzar = har * gv%H_to_Z
1154  else
1155  dzal = 0.5*(e(i,j,k-1) - e(i,j,k+1)) + dz_neglect
1156  dzar = 0.5*(e(i,j+1,k-1) - e(i,j+1,k+1)) + dz_neglect
1157  endif
1158  ! Use the harmonic mean thicknesses to weight the horizontal gradients.
1159  ! These unnormalized weights have been rearranged to minimize divisions.
1160  wtl = hg2l*(har*dzar) ; wtr = hg2r*(hal*dzal)
1161 
1162  drdz = (wtl * drdkl + wtr * drdkr) / (dzal*wtl + dzar*wtr)
1163  ! The expression for drdz above is mathematically equivalent to:
1164  ! drdz = ((hg2L/haL) * drdkL/dzaL + (hg2R/haR) * drdkR/dzaR) / &
1165  ! ((hg2L/haL) + (hg2R/haR))
1166  hg2a = h(i,j,k-1)*h(i,j+1,k-1) + h_neglect2
1167  hg2b = h(i,j,k)*h(i,j+1,k) + h_neglect2
1168  haa = 0.5*(h(i,j,k-1) + h(i,j+1,k-1)) + h_neglect
1169  hab = 0.5*(h(i,j,k) + h(i,j+1,k)) + h_neglect
1170 
1171  ! hN2_v is used with the FGNV streamfunction formulation
1172  hn2_v(i,k) = (0.5 * gv%H_to_Z * ( hg2a / haa + hg2b / hab )) * &
1173  max(drdz*g_rho0, n2_floor)
1174  endif
1175  if (present_slope_y) then
1176  slope = slope_y(i,j,k)
1177  slope2_ratio_v(i,k) = slope**2 * i_slope_max2
1178  else
1179  ! Use the harmonic mean thicknesses to weight the horizontal gradients.
1180  ! These unnormalized weights have been rearranged to minimize divisions.
1181  wta = hg2a*hab ; wtb = hg2b*haa
1182  ! This is the gradient of density along geopotentials.
1183  drdy = ((wta * drdja + wtb * drdjb) / (wta + wtb) - &
1184  drdz * (e(i,j,k)-e(i,j+1,k))) * g%IdyCv(i,j)
1185 
1186  ! This estimate of slope is accurate for small slopes, but bounded
1187  ! to be between -1 and 1.
1188  mag_grad2 = drdy**2 + (us%L_to_Z*drdz)**2
1189  if (mag_grad2 > 0.0) then
1190  slope = drdy / sqrt(mag_grad2)
1191  slope2_ratio_v(i,k) = slope**2 * i_slope_max2
1192  else ! Just in case mag_grad2 = 0 ever.
1193  slope = 0.0
1194  slope2_ratio_v(i,k) = 1.0e20 ! Force the use of the safe streamfunction.
1195  endif
1196  endif
1197 
1198  ! Adjust real slope by weights that bias towards slope of interfaces
1199  ! that ignore density gradients along layers.
1200  if (present_int_slope_v) then
1201  slope = (1.0 - int_slope_v(i,j,k)) * slope + &
1202  int_slope_v(i,j,k) * us%Z_to_L*((e(i,j+1,k)-e(i,j,k)) * g%IdyCv(i,j))
1203  slope2_ratio_v(i,k) = (1.0 - int_slope_v(i,j,k)) * slope2_ratio_v(i,k)
1204  endif
1205 
1206  slope_y_pe(i,j,k) = min(slope,cs%slope_max)
1207  hn2_y_pe(i,j,k) = hn2_v(i,k)
1208  if (cs%id_slope_y > 0) cs%diagSlopeY(i,j,k) = slope
1209 
1210  ! Estimate the streamfunction at each interface [Z L2 T-1 ~> m3 s-1].
1211  sfn_unlim_v(i,k) = -((kh_v(i,j,k)*g%dx_Cv(i,j))*us%L_to_Z*slope)
1212 
1213  ! Avoid moving dense water upslope from below the level of
1214  ! the bottom on the receiving side.
1215  if (sfn_unlim_v(i,k) > 0.0) then ! The flow below this interface is positive.
1216  if (e(i,j,k) < e(i,j+1,nz+1)) then
1217  sfn_unlim_v(i,k) = 0.0 ! This is not vhtot, because it may compensate for
1218  ! deeper flow in very unusual cases.
1219  elseif (e(i,j+1,nz+1) > e(i,j,k+1)) then
1220  ! Scale the transport with the fraction of the donor layer above
1221  ! the bottom on the receiving side.
1222  sfn_unlim_v(i,k) = sfn_unlim_v(i,k) * ((e(i,j,k) - e(i,j+1,nz+1)) / &
1223  ((e(i,j,k) - e(i,j,k+1)) + dz_neglect))
1224  endif
1225  else
1226  if (e(i,j+1,k) < e(i,j,nz+1)) then ; sfn_unlim_v(i,k) = 0.0
1227  elseif (e(i,j,nz+1) > e(i,j+1,k+1)) then
1228  sfn_unlim_v(i,k) = sfn_unlim_v(i,k) * ((e(i,j+1,k) - e(i,j,nz+1)) / &
1229  ((e(i,j+1,k) - e(i,j+1,k+1)) + dz_neglect))
1230  endif
1231  endif
1232 
1233  else ! .not. use_EOS
1234  if (present_slope_y) then
1235  slope = slope_y(i,j,k)
1236  else
1237  slope = us%Z_to_L*((e(i,j,k)-e(i,j+1,k))*g%IdyCv(i,j)) * g%mask2dCv(i,j)
1238  endif
1239  if (cs%id_slope_y > 0) cs%diagSlopeY(i,j,k) = slope
1240  sfn_unlim_v(i,k) = ((kh_v(i,j,k)*g%dx_Cv(i,j))*us%L_to_Z*slope)
1241  hn2_v(i,k) = gv%g_prime(k)
1242  endif ! if (use_EOS)
1243  else ! if (k > nk_linear)
1244  hn2_v(i,k) = n2_floor * dz_neglect
1245  sfn_unlim_v(i,k) = 0.
1246  endif ! if (k > nk_linear)
1247  if (cs%id_sfn_unlim_y>0) diag_sfn_unlim_y(i,j,k) = sfn_unlim_v(i,k)
1248  enddo ! i-loop
1249  enddo ! k-loop
1250 
1251  if (cs%use_FGNV_streamfn) then
1252  do k=1,nz ; do i=is,ie ; if (g%mask2dCv(i,j)>0.) then
1253  h_harm = max( h_neglect, &
1254  2. * h(i,j,k) * h(i,j+1,k) / ( ( h(i,j,k) + h(i,j+1,k) ) + h_neglect ) )
1255  c2_h_v(i,k) = cs%FGNV_scale * &
1256  ( 0.5*( cg1(i,j) + cg1(i,j+1) ) )**2 / (gv%H_to_Z*h_harm)
1257  endif ; enddo ; enddo
1258 
1259  ! Solve an elliptic equation for the streamfunction following Ferrari et al., 2010.
1260  do i=is,ie
1261  if (g%mask2dCv(i,j)>0.) then
1262  sfn_unlim_v(i,:) = ( 1. + cs%FGNV_scale ) * sfn_unlim_v(i,:)
1263  call streamfn_solver(nz, c2_h_v(i,:), hn2_v(i,:), sfn_unlim_v(i,:))
1264  else
1265  sfn_unlim_v(i,:) = 0.
1266  endif
1267  enddo
1268  endif
1269 
1270  do k=nz,2,-1
1271  do i=is,ie
1272  if (k > nk_linear) then
1273  if (use_eos) then
1274 
1275  if (vhtot(i,j) <= 0.0) then
1276  ! The transport that must balance the transport below is positive.
1277  sfn_safe = vhtot(i,j) * (1.0 - h_frac(i,j,k)) * gv%H_to_Z
1278  else ! (vhtot(I,j) > 0.0)
1279  sfn_safe = vhtot(i,j) * (1.0 - h_frac(i,j+1,k)) * gv%H_to_Z
1280  endif
1281 
1282  ! The actual streamfunction at each interface.
1283  sfn_est = (sfn_unlim_v(i,k) + slope2_ratio_v(i,k)*sfn_safe) / (1.0 + slope2_ratio_v(i,k))
1284  else ! With .not.use_EOS, the layers are constant density.
1285  sfn_est = sfn_unlim_v(i,k)
1286  endif
1287 
1288  ! Make sure that there is enough mass above to allow the streamfunction
1289  ! to satisfy the boundary condition of 0 at the surface.
1290  sfn_in_h = min(max(sfn_est * gv%Z_to_H, -h_avail_rsum(i,j,k)), h_avail_rsum(i,j+1,k))
1291 
1292  ! The actual transport is limited by the mass available in the two
1293  ! neighboring grid cells.
1294  vhd(i,j,k) = max(min((sfn_in_h - vhtot(i,j)), h_avail(i,j,k)), -h_avail(i,j+1,k))
1295 
1296  if (cs%id_sfn_y>0) diag_sfn_y(i,j,k) = diag_sfn_y(i,j,k+1) + vhd(i,j,k)
1297 ! sfn_y(i,J,K) = max(min(Sfn_in_h, vhtot(i,J)+h_avail(i,j,k)), &
1298 ! vhtot(i,J)-h_avail(i,j+1,k))
1299 ! sfn_slope_y(i,J,K) = max(vhtot(i,J)-h_avail(i,j+1,k), &
1300 ! min(vhtot(i,J)+h_avail(i,j,k), &
1301 ! min(h_avail_rsum(i,j+1,K), max(-h_avail_rsum(i,j,K), &
1302 ! (KH_v(i,J,K)*G%dx_Cv(i,J)) * ((e(i,j,K)-e(i,j+1,K))*G%IdyCv(i,J)) )) ))
1303  else ! k <= nk_linear
1304  ! Balance the deeper flow with a return flow uniformly distributed
1305  ! though the remaining near-surface layers. This is the same as
1306  ! using Sfn_safe above. There is no need to apply the limiters in
1307  ! this case.
1308  if (vhtot(i,j) <= 0.0) then
1309  vhd(i,j,k) = -vhtot(i,j) * h_frac(i,j,k)
1310  else ! (vhtot(i,J) > 0.0)
1311  vhd(i,j,k) = -vhtot(i,j) * h_frac(i,j+1,k)
1312  endif
1313 
1314  if (cs%id_sfn_y>0) diag_sfn_y(i,j,k) = diag_sfn_y(i,j,k+1) + vhd(i,j,k)
1315 ! if (sfn_slope_y(i,J,K+1) <= 0.0) then
1316 ! sfn_slope_y(i,J,K) = sfn_slope_y(i,J,K+1) * (1.0 - h_frac(i,j,k))
1317 ! else
1318 ! sfn_slope_y(i,J,K) = sfn_slope_y(i,J,K+1) * (1.0 - h_frac(i,j+1,k))
1319 ! endif
1320  endif
1321 
1322  vhtot(i,j) = vhtot(i,j) + vhd(i,j,k)
1323 
1324  if (find_work) then
1325  ! This is the energy tendency based on the original profiles, and does
1326  ! not include any nonlinear terms due to a finite time step (which would
1327  ! involve interactions between the fluxes through the different faces.
1328  ! A second order centered estimate is used for the density transfered
1329  ! between water columns.
1330 
1331  work_v(i,j) = work_v(i,j) + g_scale * &
1332  ( vhtot(i,j) * drdkde_v(i,k) - &
1333  (vhd(i,j,k) * drdj_v(i,k)) * 0.25 * &
1334  ((e(i,j,k) + e(i,j,k+1)) + (e(i,j+1,k) + e(i,j+1,k+1))) )
1335  endif
1336 
1337  enddo
1338  enddo ! end of k-loop
1339  enddo ! end of j-loop
1340 
1341  ! In layer 1, enforce the boundary conditions that Sfn(z=0) = 0.0
1342  if (.not.find_work .or. .not.(use_eos)) then
1343  do j=js,je ; do i=is-1,ie ; uhd(i,j,1) = -uhtot(i,j) ; enddo ; enddo
1344  do j=js-1,je ; do i=is,ie ; vhd(i,j,1) = -vhtot(i,j) ; enddo ; enddo
1345  else
1346  eosdom_u(1) = (is-1) - (g%IsdB-1) ; eosdom_u(2) = ie - (g%IsdB-1)
1347  !$OMP parallel do default(shared) private(pres_u,T_u,S_u,drho_dT_u,drho_dS_u,drdiB)
1348  do j=js,je
1349  if (use_eos) then
1350  do i=is-1,ie
1351  pres_u(i) = 0.5*(pres(i,j,1) + pres(i+1,j,1))
1352  t_u(i) = 0.5*(t(i,j,1) + t(i+1,j,1))
1353  s_u(i) = 0.5*(s(i,j,1) + s(i+1,j,1))
1354  enddo
1355  call calculate_density_derivs(t_u, s_u, pres_u, drho_dt_u, drho_ds_u, &
1356  tv%eqn_of_state, eosdom_u )
1357  endif
1358  do i=is-1,ie
1359  uhd(i,j,1) = -uhtot(i,j)
1360 
1361  if (use_eos) then
1362  drdib = drho_dt_u(i) * (t(i+1,j,1)-t(i,j,1)) + &
1363  drho_ds_u(i) * (s(i+1,j,1)-s(i,j,1))
1364  endif
1365  if (cs%use_GM_work_bug) then
1366  work_u(i,j) = work_u(i,j) + g_scale * &
1367  ( (uhd(i,j,1) * drdib) * 0.25 * &
1368  ((e(i,j,1) + e(i,j,2)) + (e(i+1,j,1) + e(i+1,j,2))) )
1369  else
1370  work_u(i,j) = work_u(i,j) - g_scale * &
1371  ( (uhd(i,j,1) * drdib) * 0.25 * &
1372  ((e(i,j,1) + e(i,j,2)) + (e(i+1,j,1) + e(i+1,j,2))) )
1373  endif
1374  enddo
1375  enddo
1376 
1377  eosdom_v(:) = eos_domain(g%HI)
1378  !$OMP parallel do default(shared) private(pres_v,T_v,S_v,drho_dT_v,drho_dS_v,drdjB)
1379  do j=js-1,je
1380  if (use_eos) then
1381  do i=is,ie
1382  pres_v(i) = 0.5*(pres(i,j,1) + pres(i,j+1,1))
1383  t_v(i) = 0.5*(t(i,j,1) + t(i,j+1,1))
1384  s_v(i) = 0.5*(s(i,j,1) + s(i,j+1,1))
1385  enddo
1386  call calculate_density_derivs(t_v, s_v, pres_v, drho_dt_v, drho_ds_v, &
1387  tv%eqn_of_state, eosdom_v)
1388  endif
1389  do i=is,ie
1390  vhd(i,j,1) = -vhtot(i,j)
1391 
1392  if (use_eos) then
1393  drdjb = drho_dt_v(i) * (t(i,j+1,1)-t(i,j,1)) + &
1394  drho_ds_v(i) * (s(i,j+1,1)-s(i,j,1))
1395  endif
1396  work_v(i,j) = work_v(i,j) - g_scale * &
1397  ( (vhd(i,j,1) * drdjb) * 0.25 * &
1398  ((e(i,j,1) + e(i,j,2)) + (e(i,j+1,1) + e(i,j+1,2))) )
1399  enddo
1400  enddo
1401  endif
1402 
1403  if (find_work) then ; do j=js,je ; do i=is,ie
1404  ! Note that the units of Work_v and Work_u are W, while Work_h is W m-2.
1405  work_h = 0.5 * g%IareaT(i,j) * &
1406  ((work_u(i-1,j) + work_u(i,j)) + (work_v(i,j-1) + work_v(i,j)))
1407  if (associated(cs%GMwork)) cs%GMwork(i,j) = work_h
1408  if (associated(meke) .and. .not.cs%GM_src_alt) then ; if (associated(meke%GM_src)) then
1409  meke%GM_src(i,j) = meke%GM_src(i,j) + work_h
1410  endif ; endif
1411  enddo ; enddo ; endif
1412 
1413  if (find_work .and. cs%GM_src_alt .and. associated(meke)) then ; if (associated(meke%GM_src)) then
1414  do j=js,je ; do i=is,ie ; do k=nz,1,-1
1415  pe_release_h = -0.25*(kh_u(i,j,k)*(slope_x_pe(i,j,k)**2) * hn2_x_pe(i,j,k) + &
1416  kh_u(i-1,j,k)*(slope_x_pe(i-1,j,k)**2) * hn2_x_pe(i-1,j,k) + &
1417  kh_v(i,j,k)*(slope_y_pe(i,j,k)**2) * hn2_y_pe(i,j,k) + &
1418  kh_v(i,j-1,k)*(slope_y_pe(i,j-1,k)**2) * hn2_y_pe(i,j-1,k))
1419  meke%GM_src(i,j) = meke%GM_src(i,j) + us%L_to_Z**2 * gv%Rho0 * pe_release_h
1420  enddo ; enddo ; enddo
1421  endif ; endif
1422 
1423  if (cs%id_slope_x > 0) call post_data(cs%id_slope_x, cs%diagSlopeX, cs%diag)
1424  if (cs%id_slope_y > 0) call post_data(cs%id_slope_y, cs%diagSlopeY, cs%diag)
1425  if (cs%id_sfn_x > 0) call post_data(cs%id_sfn_x, diag_sfn_x, cs%diag)
1426  if (cs%id_sfn_y > 0) call post_data(cs%id_sfn_y, diag_sfn_y, cs%diag)
1427  if (cs%id_sfn_unlim_x > 0) call post_data(cs%id_sfn_unlim_x, diag_sfn_unlim_x, cs%diag)
1428  if (cs%id_sfn_unlim_y > 0) call post_data(cs%id_sfn_unlim_y, diag_sfn_unlim_y, cs%diag)
1429 
1430 end subroutine thickness_diffuse_full
1431 
1432 !> Tridiagonal solver for streamfunction at interfaces
1433 subroutine streamfn_solver(nk, c2_h, hN2, sfn)
1434  integer, intent(in) :: nk !< Number of layers
1435  real, dimension(nk), intent(in) :: c2_h !< Wave speed squared over thickness in layers [L2 Z-1 T-2 ~> m s-2]
1436  real, dimension(nk+1), intent(in) :: hN2 !< Thickness times N2 at interfaces [L2 Z-1 T-2 ~> m s-2]
1437  real, dimension(nk+1), intent(inout) :: sfn !< Streamfunction [Z L2 T-1 ~> m3 s-1] or arbitrary units
1438  !! On entry, equals diffusivity times slope.
1439  !! On exit, equals the streamfunction.
1440  ! Local variables
1441  integer :: k
1442 
1443  real :: b_denom, beta, d1, c1(nk)
1444 
1445  sfn(1) = 0.
1446  b_denom = hn2(2) + c2_h(1)
1447  beta = 1.0 / ( b_denom + c2_h(2) )
1448  d1 = beta * b_denom
1449  sfn(2) = ( beta * hn2(2) )*sfn(2)
1450  do k=3,nk
1451  c1(k-1) = beta * c2_h(k-1)
1452  b_denom = hn2(k) + d1*c2_h(k-1)
1453  beta = 1.0 / (b_denom + c2_h(k))
1454  d1 = beta * b_denom
1455  sfn(k) = beta * (hn2(k)*sfn(k) + c2_h(k-1)*sfn(k-1))
1456  enddo
1457  c1(nk) = beta * c2_h(nk)
1458  sfn(nk+1) = 0.
1459  do k=nk,2,-1
1460  sfn(k) = sfn(k) + c1(k)*sfn(k+1)
1461  enddo
1462 
1463 end subroutine streamfn_solver
1464 
1465 !> Modifies thickness diffusivities to untangle layer structures
1466 subroutine add_detangling_kh(h, e, Kh_u, Kh_v, KH_u_CFL, KH_v_CFL, tv, dt, G, GV, US, CS, &
1467  int_slope_u, int_slope_v)
1468  type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1469  type(verticalgrid_type), intent(in) :: GV !< Vertical grid structure
1470  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1471  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1472  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: e !< Interface positions [Z ~> m]
1473  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kh_u !< Thickness diffusivity on interfaces
1474  !! at u points [L2 T-1 ~> m2 s-1]
1475  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(inout) :: Kh_v !< Thickness diffusivity on interfaces
1476  !! at v points [L2 T-1 ~> m2 s-1]
1477  real, dimension(SZIB_(G),SZJ_(G)), intent(in) :: Kh_u_CFL !< Maximum stable thickness diffusivity
1478  !! at u points [L2 T-1 ~> m2 s-1]
1479  real, dimension(SZI_(G),SZJB_(G)), intent(in) :: Kh_v_CFL !< Maximum stable thickness diffusivity
1480  !! at v points [L2 T-1 ~> m2 s-1]
1481  type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
1482  real, intent(in) :: dt !< Time increment [T ~> s]
1483  type(thickness_diffuse_cs), pointer :: CS !< Control structure for thickness diffusion
1484  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: int_slope_u !< Ratio that determine how much of
1485  !! the isopycnal slopes are taken directly from
1486  !! the interface slopes without consideration
1487  !! of density gradients.
1488  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(inout) :: int_slope_v !< Ratio that determine how much of
1489  !! the isopycnal slopes are taken directly from
1490  !! the interface slopes without consideration
1491  !! of density gradients.
1492  ! Local variables
1493  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
1494  de_top ! The distances between the top of a layer and the top of the
1495  ! region where the detangling is applied [H ~> m or kg m-2].
1496  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: &
1497  Kh_lay_u ! The tentative interface height diffusivity for each layer at
1498  ! u points [L2 T-1 ~> m2 s-1].
1499  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: &
1500  Kh_lay_v ! The tentative interface height diffusivity for each layer at
1501  ! v points [L2 T-1 ~> m2 s-1].
1502  real, dimension(SZI_(G),SZJ_(G)) :: &
1503  de_bot ! The distances from the bottom of the region where the
1504  ! detangling is applied [H ~> m or kg m-2].
1505  real :: h1, h2 ! The thinner and thicker surrounding thicknesses [H ~> m or kg m-2],
1506  ! with the thinner modified near the boundaries to mask out
1507  ! thickness variations due to topography, etc.
1508  real :: jag_Rat ! The nondimensional jaggedness ratio for a layer, going
1509  ! from 0 (smooth) to 1 (jagged). This is the difference
1510  ! between the arithmetic and harmonic mean thicknesses
1511  ! normalized by the arithmetic mean thickness.
1512  real :: Kh_scale ! A ratio by which Kh_u_CFL is scaled for maximally jagged
1513  ! layers [nondim].
1514  real :: h_neglect ! A thickness that is so small it is usually lost
1515  ! in roundoff and can be neglected [H ~> m or kg m-2].
1516 
1517  real :: I_sl ! The absolute value of the larger in magnitude of the slopes
1518  ! above and below [L Z-1 ~> nondim].
1519  real :: Rsl ! The ratio of the smaller magnitude slope to the larger
1520  ! magnitude one [nondim]. 0 <= Rsl <1.
1521  real :: IRsl ! The (limited) inverse of Rsl [nondim]. 1 < IRsl <= 1e9.
1522  real :: dH ! The thickness gradient divided by the damping timescale
1523  ! and the ratio of the face length to the adjacent cell
1524  ! areas for comparability with the diffusivities [L Z T-1 ~> m2 s-1].
1525  real :: adH ! The absolute value of dH [L Z T-1 ~> m2 s-1].
1526  real :: sign ! 1 or -1, with the same sign as the layer thickness gradient.
1527  real :: sl_K ! The sign-corrected slope of the interface above [Z L-1 ~> nondim].
1528  real :: sl_Kp1 ! The sign-corrected slope of the interface below [Z L-1 ~> nondim].
1529  real :: I_sl_K ! The (limited) inverse of sl_K [L Z-1 ~> nondim].
1530  real :: I_sl_Kp1 ! The (limited) inverse of sl_Kp1 [L Z-1 ~> nondim].
1531  real :: I_4t ! A quarter of a flux scaling factor divided by
1532  ! the damping timescale [T-1 ~> s-1].
1533  real :: Fn_R ! A function of Rsl, such that Rsl < Fn_R < 1.
1534  real :: denom, I_denom ! A denominator and its inverse, various units.
1535  real :: Kh_max ! A local ceiling on the diffusivity [L2 T-1 ~> m2 s-1].
1536  real :: wt1, wt2 ! Nondimensional weights.
1537  ! Variables used only in testing code.
1538  ! real, dimension(SZK_(G)) :: uh_here
1539  ! real, dimension(SZK_(G)+1) :: Sfn
1540  real :: dKh ! An increment in the diffusivity [L2 T-1 ~> m2 s-1].
1541 
1542  real, dimension(SZIB_(G),SZK_(G)+1) :: &
1543  Kh_bg, & ! The background (floor) value of Kh [L2 T-1 ~> m2 s-1].
1544  Kh, & ! The tentative value of Kh [L2 T-1 ~> m2 s-1].
1545  Kh_detangle, & ! The detangling diffusivity that could be used [L2 T-1 ~> m2 s-1].
1546  Kh_min_max_p, & ! The smallest ceiling that can be placed on Kh(I,K)
1547  ! based on the value of Kh(I,K+1) [L2 T-1 ~> m2 s-1].
1548  kh_min_max_m, & ! The smallest ceiling that can be placed on Kh(I,K)
1549  ! based on the value of Kh(I,K-1) [L2 T-1 ~> m2 s-1].
1550  ! The following are variables that define the relationships between
1551  ! successive values of Kh.
1552  ! Search for Kh that satisfy...
1553  ! Kh(I,K) >= Kh_min_m(I,K)*Kh(I,K-1) + Kh0_min_m(I,K)
1554  ! Kh(I,K) >= Kh_min_p(I,K)*Kh(I,K+1) + Kh0_min_p(I,K)
1555  ! Kh(I,K) <= Kh_max_m(I,K)*Kh(I,K-1) + Kh0_max_m(I,K)
1556  ! Kh(I,K) <= Kh_max_p(I,K)*Kh(I,K+1) + Kh0_max_p(I,K)
1557  kh_min_m , & ! See above [nondim].
1558  kh0_min_m , & ! See above [L2 T-1 ~> m2 s-1].
1559  kh_max_m , & ! See above [nondim].
1560  kh0_max_m, & ! See above [L2 T-1 ~> m2 s-1].
1561  kh_min_p , & ! See above [nondim].
1562  kh0_min_p , & ! See above [L2 T-1 ~> m2 s-1].
1563  kh_max_p , & ! See above [nondim].
1564  kh0_max_p ! See above [L2 T-1 ~> m2 s-1].
1565  real, dimension(SZIB_(G)) :: &
1566  Kh_max_max ! The maximum diffusivity permitted in a column [L2 T-1 ~> m2 s-1]..
1567  logical, dimension(SZIB_(G)) :: &
1568  do_i ! If true, work on a column.
1569  integer :: i, j, k, n, ish, jsh, is, ie, js, je, nz, k_top
1570  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1571 
1572  k_top = gv%nk_rho_varies + 1
1573  h_neglect = gv%H_subroundoff
1574  ! The 0.5 is because we are not using uniform weightings, but are
1575  ! distributing the diffusivities more effectively (with wt1 & wt2), but this
1576  ! means that the additions to a single interface can be up to twice as large.
1577  kh_scale = 0.5
1578  if (cs%detangle_time > dt) kh_scale = 0.5 * dt / cs%detangle_time
1579 
1580  do j=js-1,je+1 ; do i=is-1,ie+1
1581  de_top(i,j,k_top) = 0.0 ; de_bot(i,j) = 0.0
1582  enddo ; enddo
1583  do k=k_top+1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
1584  de_top(i,j,k) = de_top(i,j,k-1) + h(i,j,k-1)
1585  enddo ; enddo ; enddo
1586 
1587  do j=js,je ; do i=is-1,ie
1588  kh_lay_u(i,j,nz) = 0.0 ; kh_lay_u(i,j,k_top) = 0.0
1589  enddo ; enddo
1590  do j=js-1,je ; do i=is,ie
1591  kh_lay_v(i,j,nz) = 0.0 ; kh_lay_v(i,j,k_top) = 0.0
1592  enddo ; enddo
1593 
1594  do k=nz-1,k_top+1,-1
1595  ! Find the diffusivities associated with each layer.
1596  do j=js-1,je+1 ; do i=is-1,ie+1
1597  de_bot(i,j) = de_bot(i,j) + h(i,j,k+1)
1598  enddo ; enddo
1599 
1600  do j=js,je ; do i=is-1,ie ; if (g%mask2dCu(i,j) > 0.0) then
1601  if (h(i,j,k) > h(i+1,j,k)) then
1602  h2 = h(i,j,k)
1603  h1 = max( h(i+1,j,k), h2 - min(de_bot(i+1,j), de_top(i+1,j,k)) )
1604  else
1605  h2 = h(i+1,j,k)
1606  h1 = max( h(i,j,k), h2 - min(de_bot(i,j), de_top(i,j,k)) )
1607  endif
1608  jag_rat = (h2 - h1)**2 / (h2 + h1 + h_neglect)**2
1609  kh_lay_u(i,j,k) = (kh_scale * kh_u_cfl(i,j)) * jag_rat**2
1610  endif ; enddo ; enddo
1611 
1612  do j=js-1,je ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.0) then
1613  if (h(i,j,k) > h(i,j+1,k)) then
1614  h2 = h(i,j,k)
1615  h1 = max( h(i,j+1,k), h2 - min(de_bot(i,j+1), de_top(i,j+1,k)) )
1616  else
1617  h2 = h(i,j+1,k)
1618  h1 = max( h(i,j,k), h2 - min(de_bot(i,j), de_top(i,j,k)) )
1619  endif
1620  jag_rat = (h2 - h1)**2 / (h2 + h1 + h_neglect)**2
1621  kh_lay_v(i,j,k) = (kh_scale * kh_v_cfl(i,j)) * jag_rat**2
1622  endif ; enddo ; enddo
1623  enddo
1624 
1625  ! Limit the diffusivities
1626 
1627  i_4t = kh_scale / (4.0 * dt)
1628 
1629  do n=1,2
1630  if (n==1) then ; jsh = js ; ish = is-1
1631  else ; jsh = js-1 ; ish = is ; endif
1632 
1633  do j=jsh,je
1634 
1635  ! First, populate the diffusivities
1636  if (n==1) then ! This is a u-column.
1637  do i=ish,ie
1638  do_i(i) = (g%mask2dCu(i,j) > 0.0)
1639  kh_max_max(i) = kh_u_cfl(i,j)
1640  enddo
1641  do k=1,nz+1 ; do i=ish,ie
1642  kh_bg(i,k) = kh_u(i,j,k) ; kh(i,k) = kh_bg(i,k)
1643  kh_min_max_p(i,k) = kh_bg(i,k) ; kh_min_max_m(i,k) = kh_bg(i,k)
1644  kh_detangle(i,k) = 0.0
1645  enddo ; enddo
1646  else ! This is a v-column.
1647  do i=ish,ie
1648  do_i(i) = (g%mask2dCv(i,j) > 0.0) ; kh_max_max(i) = kh_v_cfl(i,j)
1649  enddo
1650  do k=1,nz+1 ; do i=ish,ie
1651  kh_bg(i,k) = kh_v(i,j,k) ; kh(i,k) = kh_bg(i,k)
1652  kh_min_max_p(i,k) = kh_bg(i,k) ; kh_min_max_m(i,k) = kh_bg(i,k)
1653  kh_detangle(i,k) = 0.0
1654  enddo ; enddo
1655  endif
1656 
1657  ! Determine the limits on the diffusivities.
1658  do k=k_top,nz ; do i=ish,ie ; if (do_i(i)) then
1659  if (n==1) then ! This is a u-column.
1660  dh = 0.0
1661  denom = ((g%IareaT(i+1,j) + g%IareaT(i,j)) * g%dy_Cu(i,j))
1662  ! This expression uses differences in e in place of h for better
1663  ! consistency with the slopes.
1664  if (denom > 0.0) &
1665  dh = i_4t * ((e(i+1,j,k) - e(i+1,j,k+1)) - &
1666  (e(i,j,k) - e(i,j,k+1))) / denom
1667  ! dH = I_4t * (h(i+1,j,k) - h(i,j,k)) / denom
1668 
1669  adh = abs(dh)
1670  sign = 1.0 ; if (dh < 0) sign = -1.0
1671  sl_k = sign * (e(i+1,j,k)-e(i,j,k)) * g%IdxCu(i,j)
1672  sl_kp1 = sign * (e(i+1,j,k+1)-e(i,j,k+1)) * g%IdxCu(i,j)
1673 
1674  ! Add the incremental diffusivites to the surrounding interfaces.
1675  ! Adding more to the more steeply sloping layers (as below) makes
1676  ! the diffusivities more than twice as effective.
1677  denom = (sl_k**2 + sl_kp1**2)
1678  wt1 = 0.5 ; wt2 = 0.5
1679  if (denom > 0.0) then
1680  wt1 = sl_k**2 / denom ; wt2 = sl_kp1**2 / denom
1681  endif
1682  kh_detangle(i,k) = kh_detangle(i,k) + wt1*kh_lay_u(i,j,k)
1683  kh_detangle(i,k+1) = kh_detangle(i,k+1) + wt2*kh_lay_u(i,j,k)
1684  else ! This is a v-column.
1685  dh = 0.0
1686  denom = ((g%IareaT(i,j+1) + g%IareaT(i,j)) * g%dx_Cv(i,j))
1687  if (denom > 0.0) &
1688  dh = i_4t * ((e(i,j+1,k) - e(i,j+1,k+1)) - &
1689  (e(i,j,k) - e(i,j,k+1))) / denom
1690  ! dH = I_4t * (h(i,j+1,k) - h(i,j,k)) / denom
1691 
1692  adh = abs(dh)
1693  sign = 1.0 ; if (dh < 0) sign = -1.0
1694  sl_k = sign * (e(i,j+1,k)-e(i,j,k)) * g%IdyCv(i,j)
1695  sl_kp1 = sign * (e(i,j+1,k+1)-e(i,j,k+1)) * g%IdyCv(i,j)
1696 
1697  ! Add the incremental diffusviites to the surrounding interfaces.
1698  ! Adding more to the more steeply sloping layers (as below) makes
1699  ! the diffusivities more than twice as effective.
1700  denom = (sl_k**2 + sl_kp1**2)
1701  wt1 = 0.5 ; wt2 = 0.5
1702  if (denom > 0.0) then
1703  wt1 = sl_k**2 / denom ; wt2 = sl_kp1**2 / denom
1704  endif
1705  kh_detangle(i,k) = kh_detangle(i,k) + wt1*kh_lay_v(i,j,k)
1706  kh_detangle(i,k+1) = kh_detangle(i,k+1) + wt2*kh_lay_v(i,j,k)
1707  endif
1708 
1709  if (adh == 0.0) then
1710  kh_min_m(i,k+1) = 1.0 ; kh0_min_m(i,k+1) = 0.0
1711  kh_max_m(i,k+1) = 1.0 ; kh0_max_m(i,k+1) = 0.0
1712  kh_min_p(i,k) = 1.0 ; kh0_min_p(i,k) = 0.0
1713  kh_max_p(i,k) = 1.0 ; kh0_max_p(i,k) = 0.0
1714  elseif (adh > 0.0) then
1715  if (sl_k <= sl_kp1) then
1716  ! This case should only arise from nonlinearities in the equation of state.
1717  ! Treat it as though dedx(K) = dedx(K+1) & dH = 0.
1718  kh_min_m(i,k+1) = 1.0 ; kh0_min_m(i,k+1) = 0.0
1719  kh_max_m(i,k+1) = 1.0 ; kh0_max_m(i,k+1) = 0.0
1720  kh_min_p(i,k) = 1.0 ; kh0_min_p(i,k) = 0.0
1721  kh_max_p(i,k) = 1.0 ; kh0_max_p(i,k) = 0.0
1722  elseif (sl_k <= 0.0) then ! Both slopes are opposite to dH
1723  i_sl = -1.0 / sl_kp1
1724  rsl = -sl_k * i_sl ! 0 <= Rsl < 1
1725  irsl = 1e9 ; if (rsl > 1e-9) irsl = 1.0/rsl ! 1 < IRsl <= 1e9
1726 
1727  fn_r = rsl
1728  if (kh_max_max(i) > 0) &
1729  fn_r = min(sqrt(rsl), rsl + (adh * i_sl) / (kh_max_max(i)))
1730 
1731  kh_min_m(i,k+1) = fn_r ; kh0_min_m(i,k+1) = 0.0
1732  kh_max_m(i,k+1) = rsl ; kh0_max_m(i,k+1) = adh * i_sl
1733  kh_min_p(i,k) = irsl ; kh0_min_p(i,k) = -adh * (i_sl*irsl)
1734  kh_max_p(i,k) = 1.0/(fn_r + 1.0e-30) ; kh0_max_p(i,k) = 0.0
1735  elseif (sl_kp1 < 0.0) then ! Opposite (nonzero) signs of slopes.
1736  i_sl_k = 1e18*us%Z_to_L ; if (sl_k > 1e-18*us%L_to_Z) i_sl_k = 1.0 / sl_k
1737  i_sl_kp1 = 1e18*us%Z_to_L ; if (-sl_kp1 > 1e-18*us%L_to_Z) i_sl_kp1 = -1.0 / sl_kp1
1738 
1739  kh_min_m(i,k+1) = 0.0 ; kh0_min_m(i,k+1) = 0.0
1740  kh_max_m(i,k+1) = - sl_k*i_sl_kp1 ; kh0_max_m(i,k+1) = adh*i_sl_kp1
1741  kh_min_p(i,k) = 0.0 ; kh0_min_p(i,k) = 0.0
1742  kh_max_p(i,k) = sl_kp1*i_sl_k ; kh0_max_p(i,k) = adh*i_sl_k
1743 
1744  ! This limit does not use the slope weighting so that potentially
1745  ! sharp gradients in diffusivities are not forced to occur.
1746  kh_max = adh / (sl_k - sl_kp1)
1747  kh_min_max_p(i,k) = max(kh_min_max_p(i,k), kh_max)
1748  kh_min_max_m(i,k+1) = max(kh_min_max_m(i,k+1), kh_max)
1749  else ! Both slopes are of the same sign as dH.
1750  i_sl = 1.0 / sl_k
1751  rsl = sl_kp1 * i_sl ! 0 <= Rsl < 1
1752  irsl = 1e9 ; if (rsl > 1e-9) irsl = 1.0/rsl ! 1 < IRsl <= 1e9
1753 
1754  ! Rsl <= Fn_R <= 1
1755  fn_r = rsl
1756  if (kh_max_max(i) > 0) &
1757  fn_r = min(sqrt(rsl), rsl + (adh * i_sl) / kh_max_max(i))
1758 
1759  kh_min_m(i,k+1) = irsl ; kh0_min_m(i,k+1) = -adh * (i_sl*irsl)
1760  kh_max_m(i,k+1) = 1.0/(fn_r + 1.0e-30) ; kh0_max_m(i,k+1) = 0.0
1761  kh_min_p(i,k) = fn_r ; kh0_min_p(i,k) = 0.0
1762  kh_max_p(i,k) = rsl ; kh0_max_p(i,k) = adh * i_sl
1763  endif
1764  endif
1765  endif ; enddo ; enddo ! I-loop & k-loop
1766 
1767  do k=k_top,nz+1,nz+1-k_top ; do i=ish,ie ; if (do_i(i)) then
1768  ! The diffusivities at k_top and nz+1 are both fixed.
1769  kh_min_m(i,k) = 0.0 ; kh0_min_m(i,k) = 0.0
1770  kh_max_m(i,k) = 0.0 ; kh0_max_m(i,k) = 0.0
1771  kh_min_p(i,k) = 0.0 ; kh0_min_p(i,k) = 0.0
1772  kh_max_p(i,k) = 0.0 ; kh0_max_p(i,k) = 0.0
1773  kh_min_max_p(i,k) = kh_bg(i,k)
1774  kh_min_max_m(i,k) = kh_bg(i,k)
1775  endif ; enddo ; enddo ! I-loop and k_top/nz+1 loop
1776 
1777  ! Search for Kh that satisfy...
1778  ! Kh(I,K) >= Kh_min_m(I,K)*Kh(I,K-1) + Kh0_min_m(I,K)
1779  ! Kh(I,K) >= Kh_min_p(I,K)*Kh(I,K+1) + Kh0_min_p(I,K)
1780  ! Kh(I,K) <= Kh_max_m(I,K)*Kh(I,K-1) + Kh0_max_m(I,K)
1781  ! Kh(I,K) <= Kh_max_p(I,K)*Kh(I,K+1) + Kh0_max_p(I,K)
1782 
1783  ! Increase the diffusivies to satisfy the min constraints.
1784  ! All non-zero min constraints on one diffusivity are max constraints on another.
1785  do k=k_top+1,nz ; do i=ish,ie ; if (do_i(i)) then
1786  kh(i,k) = max(kh_bg(i,k), kh_detangle(i,k), &
1787  min(kh_min_m(i,k)*kh(i,k-1) + kh0_min_m(i,k), kh(i,k-1)))
1788 
1789  if (kh0_max_m(i,k) > kh_bg(i,k)) kh(i,k) = min(kh(i,k), kh0_max_m(i,k))
1790  if (kh0_max_p(i,k) > kh_bg(i,k)) kh(i,k) = min(kh(i,k), kh0_max_p(i,k))
1791  endif ; enddo ; enddo ! I-loop & k-loop
1792  ! This is still true... do i=ish,ie ; Kh(I,nz+1) = Kh_bg(I,nz+1) ; enddo
1793  do k=nz,k_top+1,-1 ; do i=ish,ie ; if (do_i(i)) then
1794  kh(i,k) = max(kh(i,k), min(kh_min_p(i,k)*kh(i,k+1) + kh0_min_p(i,k), kh(i,k+1)))
1795 
1796  kh_max = max(kh_min_max_p(i,k), kh_max_p(i,k)*kh(i,k+1) + kh0_max_p(i,k))
1797  kh(i,k) = min(kh(i,k), kh_max)
1798  endif ; enddo ; enddo ! I-loop & k-loop
1799  ! All non-zero min constraints on one diffusivity are max constraints on
1800  ! another layer, so the min constraints can now be discounted.
1801 
1802  ! Decrease the diffusivities to satisfy the max constraints.
1803  do k=k_top+1,nz ; do i=ish,ie ; if (do_i(i)) then
1804  kh_max = max(kh_min_max_m(i,k), kh_max_m(i,k)*kh(i,k-1) + kh0_max_m(i,k))
1805  if (kh(i,k) > kh_max) kh(i,k) = kh_max
1806  endif ; enddo ; enddo ! i- and K-loops
1807 
1808  ! This code tests the solutions...
1809 ! do i=ish,ie
1810 ! Sfn(:) = 0.0 ; uh_here(:) = 0.0
1811 ! do K=k_top,nz
1812 ! if ((Kh(i,K) > Kh_bg(i,K)) .or. (Kh(i,K+1) > Kh_bg(i,K+1))) then
1813 ! if (n==1) then ! u-point.
1814 ! if ((h(i+1,j,k) - h(i,j,k)) * &
1815 ! ((e(i+1,j,K)-e(i+1,j,K+1)) - (e(i,j,K)-e(i,j,K+1))) > 0.0) then
1816 ! Sfn(K) = -Kh(i,K) * (e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)
1817 ! Sfn(K+1) = -Kh(i,K+1) * (e(i+1,j,K+1)-e(i,j,K+1)) * G%IdxCu(I,j)
1818 ! uh_here(k) = (Sfn(K) - Sfn(K+1))*G%dy_Cu(I,j)
1819 ! if (abs(uh_here(k)) * min(G%IareaT(i,j), G%IareaT(i+1,j)) > &
1820 ! (1e-10*GV%m_to_H)) then
1821 ! if (uh_here(k) * (h(i+1,j,k) - h(i,j,k)) > 0.0) then
1822 ! call MOM_error(WARNING, "Corrective u-transport is up the thickness gradient.", .true.)
1823 ! endif
1824 ! if (((h(i,j,k) - 4.0*dt*G%IareaT(i,j)*uh_here(k)) - &
1825 ! (h(i+1,j,k) + 4.0*dt*G%IareaT(i+1,j)*uh_here(k))) * &
1826 ! (h(i,j,k) - h(i+1,j,k)) < 0.0) then
1827 ! call MOM_error(WARNING, "Corrective u-transport is too large.", .true.)
1828 ! endif
1829 ! endif
1830 ! endif
1831 ! else ! v-point
1832 ! if ((h(i,j+1,k) - h(i,j,k)) * &
1833 ! ((e(i,j+1,K)-e(i,j+1,K+1)) - (e(i,j,K)-e(i,j,K+1))) > 0.0) then
1834 ! Sfn(K) = -Kh(i,K) * (e(i,j+1,K)-e(i,j,K)) * G%IdyCv(i,J)
1835 ! Sfn(K+1) = -Kh(i,K+1) * (e(i,j+1,K+1)-e(i,j,K+1)) * G%IdyCv(i,J)
1836 ! uh_here(k) = (Sfn(K) - Sfn(K+1))*G%dx_Cv(i,J)
1837 ! if (abs(uh_here(K)) * min(G%IareaT(i,j), G%IareaT(i,j+1)) > &
1838 ! (1e-10*GV%m_to_H)) then
1839 ! if (uh_here(K) * (h(i,j+1,k) - h(i,j,k)) > 0.0) then
1840 ! call MOM_error(WARNING, &
1841 ! "Corrective v-transport is up the thickness gradient.", .true.)
1842 ! endif
1843 ! if (((h(i,j,k) - 4.0*dt*G%IareaT(i,j)*uh_here(K)) - &
1844 ! (h(i,j+1,k) + 4.0*dt*G%IareaT(i,j+1)*uh_here(K))) * &
1845 ! (h(i,j,k) - h(i,j+1,k)) < 0.0) then
1846 ! call MOM_error(WARNING, &
1847 ! "Corrective v-transport is too large.", .true.)
1848 ! endif
1849 ! endif
1850 ! endif
1851 ! endif ! u- or v- selection.
1852 ! ! de_dx(I,K) = (e(i+1,j,K)-e(i,j,K)) * G%IdxCu(I,j)
1853 ! endif
1854 ! enddo
1855 ! enddo
1856 
1857  if (n==1) then ! This is a u-column.
1858  do k=k_top+1,nz ; do i=ish,ie
1859  if (kh(i,k) > kh_u(i,j,k)) then
1860  dkh = (kh(i,k) - kh_u(i,j,k))
1861  int_slope_u(i,j,k) = dkh / kh(i,k)
1862  kh_u(i,j,k) = kh(i,k)
1863  endif
1864  enddo ; enddo
1865  else ! This is a v-column.
1866  do k=k_top+1,nz ; do i=ish,ie
1867  if (kh(i,k) > kh_v(i,j,k)) then
1868  dkh = kh(i,k) - kh_v(i,j,k)
1869  int_slope_v(i,j,k) = dkh / kh(i,k)
1870  kh_v(i,j,k) = kh(i,k)
1871  endif
1872  enddo ; enddo
1873  endif
1874 
1875  enddo ! j-loop
1876  enddo ! n-loop over u- and v- directions.
1877 
1878 end subroutine add_detangling_kh
1879 
1880 !> Initialize the thickness diffusion module/structure
1881 subroutine thickness_diffuse_init(Time, G, GV, US, param_file, diag, CDp, CS)
1882  type(time_type), intent(in) :: time !< Current model time
1883  type(ocean_grid_type), intent(in) :: g !< Ocean grid structure
1884  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure
1885  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1886  type(param_file_type), intent(in) :: param_file !< Parameter file handles
1887  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
1888  type(cont_diag_ptrs), intent(inout) :: cdp !< Continuity equation diagnostics
1889  type(thickness_diffuse_cs), pointer :: cs !< Control structure for thickness diffusion
1890 
1891 ! This include declares and sets the variable "version".
1892 #include "version_variable.h"
1893  character(len=40) :: mdl = "MOM_thickness_diffuse" ! This module's name.
1894  real :: omega ! The Earth's rotation rate [T-1 ~> s-1]
1895  real :: strat_floor ! A floor for Brunt-Vasaila frequency in the Ferrari et al. 2010,
1896  ! streamfunction formulation, expressed as a fraction of planetary
1897  ! rotation [nondim].
1898  logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags.
1899 
1900  if (associated(cs)) then
1901  call mom_error(warning, &
1902  "Thickness_diffuse_init called with an associated control structure.")
1903  return
1904  else ; allocate(cs) ; endif
1905 
1906  cs%diag => diag
1907 
1908  ! Read all relevant parameters and write them to the model log.
1909  call log_version(param_file, mdl, version, "")
1910  call get_param(param_file, mdl, "THICKNESSDIFFUSE", cs%thickness_diffuse, &
1911  "If true, interface heights are diffused with a "//&
1912  "coefficient of KHTH.", default=.false.)
1913  call get_param(param_file, mdl, "KHTH", cs%Khth, &
1914  "The background horizontal thickness diffusivity.", &
1915  default=0.0, units="m2 s-1", scale=us%m_to_L**2*us%T_to_s)
1916  call get_param(param_file, mdl, "KHTH_SLOPE_CFF", cs%KHTH_Slope_Cff, &
1917  "The nondimensional coefficient in the Visbeck formula "//&
1918  "for the interface depth diffusivity", units="nondim", &
1919  default=0.0)
1920  call get_param(param_file, mdl, "KHTH_MIN", cs%KHTH_Min, &
1921  "The minimum horizontal thickness diffusivity.", &
1922  default=0.0, units="m2 s-1", scale=us%m_to_L**2*us%T_to_s)
1923  call get_param(param_file, mdl, "KHTH_MAX", cs%KHTH_Max, &
1924  "The maximum horizontal thickness diffusivity.", &
1925  default=0.0, units="m2 s-1", scale=us%m_to_L**2*us%T_to_s)
1926  call get_param(param_file, mdl, "KHTH_MAX_CFL", cs%max_Khth_CFL, &
1927  "The maximum value of the local diffusive CFL ratio that "//&
1928  "is permitted for the thickness diffusivity. 1.0 is the "//&
1929  "marginally unstable value in a pure layered model, but "//&
1930  "much smaller numbers (e.g. 0.1) seem to work better for "//&
1931  "ALE-based models.", units = "nondimensional", default=0.8)
1932 
1933 ! call get_param(param_file, mdl, "USE_QG_LEITH_GM", CS%QG_Leith_GM, &
1934 ! "If true, use the QG Leith viscosity as the GM coefficient.", &
1935 ! default=.false.)
1936 
1937  if (cs%max_Khth_CFL < 0.0) cs%max_Khth_CFL = 0.0
1938  call get_param(param_file, mdl, "DETANGLE_INTERFACES", cs%detangle_interfaces, &
1939  "If defined add 3-d structured enhanced interface height "//&
1940  "diffusivities to horizontally smooth jagged layers.", &
1941  default=.false.)
1942  cs%detangle_time = 0.0
1943  if (cs%detangle_interfaces) &
1944  call get_param(param_file, mdl, "DETANGLE_TIMESCALE", cs%detangle_time, &
1945  "A timescale over which maximally jagged grid-scale "//&
1946  "thickness variations are suppressed. This must be "//&
1947  "longer than DT, or 0 to use DT.", units="s", default=0.0, scale=us%s_to_T)
1948  call get_param(param_file, mdl, "KHTH_SLOPE_MAX", cs%slope_max, &
1949  "A slope beyond which the calculated isopycnal slope is "//&
1950  "not reliable and is scaled away.", units="nondim", default=0.01)
1951  call get_param(param_file, mdl, "KD_SMOOTH", cs%kappa_smooth, &
1952  "A diapycnal diffusivity that is used to interpolate "//&
1953  "more sensible values of T & S into thin layers.", &
1954  units="m2 s-1", default=1.0e-6, scale=us%m_to_Z**2*us%T_to_s)
1955  call get_param(param_file, mdl, "KHTH_USE_FGNV_STREAMFUNCTION", cs%use_FGNV_streamfn, &
1956  "If true, use the streamfunction formulation of "//&
1957  "Ferrari et al., 2010, which effectively emphasizes "//&
1958  "graver vertical modes by smoothing in the vertical.", &
1959  default=.false.)
1960  call get_param(param_file, mdl, "FGNV_FILTER_SCALE", cs%FGNV_scale, &
1961  "A coefficient scaling the vertical smoothing term in the "//&
1962  "Ferrari et al., 2010, streamfunction formulation.", &
1963  units="nondim", default=1., do_not_log=.not.cs%use_FGNV_streamfn)
1964  call get_param(param_file, mdl, "FGNV_C_MIN", cs%FGNV_c_min, &
1965  "A minium wave speed used in the Ferrari et al., 2010, "//&
1966  "streamfunction formulation.", &
1967  default=0., units="m s-1", scale=us%m_s_to_L_T, do_not_log=.not.cs%use_FGNV_streamfn)
1968  call get_param(param_file, mdl, "FGNV_STRAT_FLOOR", strat_floor, &
1969  "A floor for Brunt-Vasaila frequency in the Ferrari et al., 2010, "//&
1970  "streamfunction formulation, expressed as a fraction of planetary "//&
1971  "rotation, OMEGA. This should be tiny but non-zero to avoid degeneracy.", &
1972  default=1.e-15, units="nondim", do_not_log=.not.cs%use_FGNV_streamfn)
1973  call get_param(param_file, mdl, "STANLEY_PRM_DET_COEFF", cs%Stanley_det_coeff, &
1974  "The coefficient correlating SGS temperature variance with the mean "//&
1975  "temperature gradient in the deterministic part of the Stanley parameterization. "//&
1976  "Negative values disable the scheme.", units="nondim", default=-1.0)
1977  call get_param(param_file, mdl, "OMEGA", omega, &
1978  "The rotation rate of the earth.", &
1979  default=7.2921e-5, units="s-1", scale=us%T_to_s, do_not_log=.not.cs%use_FGNV_streamfn)
1980  if (cs%use_FGNV_streamfn) cs%N2_floor = (strat_floor*omega)**2
1981  call get_param(param_file, mdl, "DEBUG", cs%debug, &
1982  "If true, write out verbose debugging data.", &
1983  default=.false., debuggingparam=.true.)
1984 
1985  call get_param(param_file, mdl, "MEKE_GM_SRC_ALT", cs%GM_src_alt, &
1986  "If true, use the GM energy conversion form S^2*N^2*kappa rather "//&
1987  "than the streamfunction for the GM source term.", default=.false.)
1988  call get_param(param_file, mdl, "MEKE_GEOMETRIC", cs%MEKE_GEOMETRIC, &
1989  "If true, uses the GM coefficient formulation from the GEOMETRIC "//&
1990  "framework (Marshall et al., 2012).", default=.false.)
1991  if (cs%MEKE_GEOMETRIC) then
1992  call get_param(param_file, mdl, "MEKE_GEOMETRIC_EPSILON", cs%MEKE_GEOMETRIC_epsilon, &
1993  "Minimum Eady growth rate used in the calculation of GEOMETRIC "//&
1994  "thickness diffusivity.", units="s-1", default=1.0e-7, scale=us%T_to_s)
1995  call get_param(param_file, mdl, "MEKE_GEOMETRIC_ALPHA", cs%MEKE_GEOMETRIC_alpha, &
1996  "The nondimensional coefficient governing the efficiency of the GEOMETRIC "//&
1997  "thickness diffusion.", units="nondim", default=0.05)
1998 
1999  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
2000  "This sets the default value for the various _2018_ANSWERS parameters.", &
2001  default=.false.)
2002  call get_param(param_file, mdl, "MEKE_GEOMETRIC_2018_ANSWERS", cs%MEKE_GEOM_answers_2018, &
2003  "If true, use expressions in the MEKE_GEOMETRIC calculation that recover the "//&
2004  "answers from the original implementation. Otherwise, use expressions that "//&
2005  "satisfy rotational symmetry.", default=default_2018_answers)
2006  endif
2007 
2008  call get_param(param_file, mdl, "USE_KH_IN_MEKE", cs%Use_KH_in_MEKE, &
2009  "If true, uses the thickness diffusivity calculated here to diffuse MEKE.", &
2010  default=.false.)
2011 
2012  call get_param(param_file, mdl, "USE_GME", cs%use_GME_thickness_diffuse, &
2013  "If true, use the GM+E backscatter scheme in association "//&
2014  "with the Gent and McWilliams parameterization.", default=.false.)
2015 
2016  call get_param(param_file, mdl, "USE_GM_WORK_BUG", cs%use_GM_work_bug, &
2017  "If true, compute the top-layer work tendency on the u-grid "//&
2018  "with the incorrect sign, for legacy reproducibility.", &
2019  default=.false.)
2020 
2021  if (cs%use_GME_thickness_diffuse) then
2022  call safe_alloc_ptr(cs%KH_u_GME,g%IsdB,g%IedB,g%jsd,g%jed,g%ke+1)
2023  call safe_alloc_ptr(cs%KH_v_GME,g%isd,g%ied,g%JsdB,g%JedB,g%ke+1)
2024  endif
2025 
2026  cs%id_uhGM = register_diag_field('ocean_model', 'uhGM', diag%axesCuL, time, &
2027  'Time Mean Diffusive Zonal Thickness Flux', &
2028  'kg s-1', conversion=gv%H_to_kg_m2*us%L_to_m**2*us%s_to_T, &
2029  y_cell_method='sum', v_extensive=.true.)
2030  if (cs%id_uhGM > 0) call safe_alloc_ptr(cdp%uhGM,g%IsdB,g%IedB,g%jsd,g%jed,g%ke)
2031  cs%id_vhGM = register_diag_field('ocean_model', 'vhGM', diag%axesCvL, time, &
2032  'Time Mean Diffusive Meridional Thickness Flux', &
2033  'kg s-1', conversion=gv%H_to_kg_m2*us%L_to_m**2*us%s_to_T, &
2034  x_cell_method='sum', v_extensive=.true.)
2035  if (cs%id_vhGM > 0) call safe_alloc_ptr(cdp%vhGM,g%isd,g%ied,g%JsdB,g%JedB,g%ke)
2036 
2037  cs%id_GMwork = register_diag_field('ocean_model', 'GMwork', diag%axesT1, time, &
2038  'Integrated Tendency of Ocean Mesoscale Eddy KE from Parameterized Eddy Advection', &
2039  'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2, cmor_field_name='tnkebto', &
2040  cmor_long_name='Integrated Tendency of Ocean Mesoscale Eddy KE from Parameterized Eddy Advection', &
2041  cmor_standard_name='tendency_of_ocean_eddy_kinetic_energy_content_due_to_parameterized_eddy_advection')
2042  if (cs%id_GMwork > 0) call safe_alloc_ptr(cs%GMwork,g%isd,g%ied,g%jsd,g%jed)
2043 
2044  cs%id_KH_u = register_diag_field('ocean_model', 'KHTH_u', diag%axesCui, time, &
2045  'Parameterized mesoscale eddy advection diffusivity at U-point', &
2046  'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
2047  cs%id_KH_v = register_diag_field('ocean_model', 'KHTH_v', diag%axesCvi, time, &
2048  'Parameterized mesoscale eddy advection diffusivity at V-point', &
2049  'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
2050  cs%id_KH_t = register_diag_field('ocean_model', 'KHTH_t', diag%axesTL, time, &
2051  'Ocean Tracer Diffusivity due to Parameterized Mesoscale Advection', &
2052  'm2 s-1', conversion=us%L_to_m**2*us%s_to_T, &
2053  cmor_field_name='diftrblo', &
2054  cmor_long_name='Ocean Tracer Diffusivity due to Parameterized Mesoscale Advection', &
2055  cmor_units='m2 s-1', &
2056  cmor_standard_name='ocean_tracer_diffusivity_due_to_parameterized_mesoscale_advection')
2057 
2058  cs%id_KH_u1 = register_diag_field('ocean_model', 'KHTH_u1', diag%axesCu1, time, &
2059  'Parameterized mesoscale eddy advection diffusivity at U-points (2-D)', &
2060  'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
2061  cs%id_KH_v1 = register_diag_field('ocean_model', 'KHTH_v1', diag%axesCv1, time, &
2062  'Parameterized mesoscale eddy advection diffusivity at V-points (2-D)', &
2063  'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
2064  cs%id_KH_t1 = register_diag_field('ocean_model', 'KHTH_t1', diag%axesT1, time, &
2065  'Parameterized mesoscale eddy advection diffusivity at T-points (2-D)', &
2066  'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)
2067 
2068  cs%id_slope_x = register_diag_field('ocean_model', 'neutral_slope_x', diag%axesCui, time, &
2069  'Zonal slope of neutral surface', 'nondim')
2070  if (cs%id_slope_x > 0) call safe_alloc_ptr(cs%diagSlopeX,g%IsdB,g%IedB,g%jsd,g%jed,g%ke+1)
2071  cs%id_slope_y = register_diag_field('ocean_model', 'neutral_slope_y', diag%axesCvi, time, &
2072  'Meridional slope of neutral surface', 'nondim')
2073  if (cs%id_slope_y > 0) call safe_alloc_ptr(cs%diagSlopeY,g%isd,g%ied,g%JsdB,g%JedB,g%ke+1)
2074  cs%id_sfn_x = register_diag_field('ocean_model', 'GM_sfn_x', diag%axesCui, time, &
2075  'Parameterized Zonal Overturning Streamfunction', &
2076  'm3 s-1', conversion=gv%H_to_m*us%L_to_m**2*us%s_to_T)
2077  cs%id_sfn_y = register_diag_field('ocean_model', 'GM_sfn_y', diag%axesCvi, time, &
2078  'Parameterized Meridional Overturning Streamfunction', &
2079  'm3 s-1', conversion=gv%H_to_m*us%L_to_m**2*us%s_to_T)
2080  cs%id_sfn_unlim_x = register_diag_field('ocean_model', 'GM_sfn_unlim_x', diag%axesCui, time, &
2081  'Parameterized Zonal Overturning Streamfunction before limiting/smoothing', &
2082  'm3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2083  cs%id_sfn_unlim_y = register_diag_field('ocean_model', 'GM_sfn_unlim_y', diag%axesCvi, time, &
2084  'Parameterized Meridional Overturning Streamfunction before limiting/smoothing', &
2085  'm3 s-1', conversion=us%Z_to_m*us%L_to_m**2*us%s_to_T)
2086 
2087 end subroutine thickness_diffuse_init
2088 
2089 !> Copies ubtav and vbtav from private type into arrays
2090 subroutine thickness_diffuse_get_kh(CS, KH_u_GME, KH_v_GME, G)
2091  type(thickness_diffuse_cs), pointer :: cs !< Control structure for
2092  !! this module
2093  type(ocean_grid_type), intent(in) :: g !< Grid structure
2094  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: kh_u_gme !< interface height
2095  !! diffusivities at u-faces [L2 T-1 ~> m2 s-1]
2096  real, dimension(SZI_(G),SZJB_(G),SZK_(G)+1), intent(inout) :: kh_v_gme !< interface height
2097  !! diffusivities at v-faces [L2 T-1 ~> m2 s-1]
2098  ! Local variables
2099  integer :: i,j,k
2100 
2101  do k=1,g%ke+1 ; do j = g%jsc, g%jec ; do i = g%isc-1, g%iec
2102  kh_u_gme(i,j,k) = cs%KH_u_GME(i,j,k)
2103  enddo ; enddo ; enddo
2104 
2105  do k=1,g%ke+1 ; do j = g%jsc-1, g%jec ; do i = g%isc, g%iec
2106  kh_v_gme(i,j,k) = cs%KH_v_GME(i,j,k)
2107  enddo ; enddo ; enddo
2108 
2109 end subroutine thickness_diffuse_get_kh
2110 
2111 !> Deallocate the thickness diffusion control structure
2112 subroutine thickness_diffuse_end(CS)
2113  type(thickness_diffuse_cs), pointer :: cs !< Control structure for thickness diffusion
2114 
2115  if (associated(cs)) deallocate(cs)
2116 end subroutine thickness_diffuse_end
2117 
2118 !> \namespace mom_thickness_diffuse
2119 !!
2120 !! \section section_gm Thickness diffusion (aka Gent-McWilliams)
2121 !!
2122 !! Thickness diffusion is implemented via along-layer mass fluxes
2123 !! \f[
2124 !! h^\dagger \leftarrow h^n - \Delta t \nabla \cdot ( \vec{uh}^* )
2125 !! \f]
2126 !! where the mass fluxes are cast as the difference in vector streamfunction
2127 !!
2128 !! \f[
2129 !! \vec{uh}^* = \delta_k \vec{\psi} .
2130 !! \f]
2131 !!
2132 !! The GM implementation of thickness diffusion made the streamfunction proportional to the potential density slope
2133 !! \f[
2134 !! \vec{\psi} = - \kappa_h \frac{\nabla_z \rho}{\partial_z \rho}
2135 !! = \frac{g\kappa_h}{\rho_o} \frac{\nabla \rho}{N^2} = \kappa_h \frac{M^2}{N^2}
2136 !! \f]
2137 !! but for robustness the scheme is implemented as
2138 !! \f[
2139 !! \vec{\psi} = \kappa_h \frac{M^2}{\sqrt{N^4 + M^4}}
2140 !! \f]
2141 !! since the quantity \f$\frac{M^2}{\sqrt{N^2 + M^2}}\f$ is bounded between $-1$ and $1$ and does not change sign
2142 !! if \f$N^2<0\f$.
2143 !!
2144 !! Optionally, the method of Ferrari et al, 2010, can be used to obtain the streamfunction which solves the
2145 !! vertically elliptic equation:
2146 !! \f[
2147 !! \gamma_F \partial_z c^2 \partial_z \psi - N_*^2 \psi = ( 1 + \gamma_F ) \kappa_h N_*^2 \frac{M^2}{\sqrt{N^4+M^4}}
2148 !! \f]
2149 !! which recovers the previous streamfunction relation in the limit that \f$ c \rightarrow 0 \f$.
2150 !! Here, \f$c=\max(c_{min},c_g)\f$ is the maximum of either \f$c_{min}\f$ and either the first baroclinic mode
2151 !! wave-speed or the equivalent barotropic mode wave-speed.
2152 !! \f$N_*^2 = \max(N^2,0)\f$ is a non-negative form of the square of the Brunt-Vaisala frequency.
2153 !! The parameter \f$\gamma_F\f$ is used to reduce the vertical smoothing length scale.
2154 !! \f[
2155 !! \kappa_h = \left( \kappa_o + \alpha_{s} L_{s}^2 < S N > + \alpha_{M} \kappa_{M} \right) r(\Delta x,L_d)
2156 !! \f]
2157 !! where \f$ S \f$ is the isoneutral slope magnitude, \f$ N \f$ is the square root of Brunt-Vaisala frequency,
2158 !! \f$\kappa_{M}\f$ is the diffusivity calculated by the MEKE parameterization (mom_meke module) and
2159 !! \f$ r(\Delta x,L_d) \f$ is a function of the local resolution (ratio of grid-spacing, \f$\Delta x\f$,
2160 !! to deformation radius, \f$L_d\f$). The length \f$L_s\f$ is provided by the mom_lateral_mixing_coeffs module
2161 !! (enabled with <code>USE_VARIABLE_MIXING=True</code> and the term \f$<SN>\f$ is the vertical average slope
2162 !! times the Brunt-Vaisala frequency prescribed by Visbeck et al., 1996.
2163 !!
2164 !! The result of the above expression is subsequently bounded by minimum and maximum values, including an upper
2165 !! diffusivity consistent with numerical stability (\f$ \kappa_{cfl} \f$ is calculated internally).
2166 !! \f[
2167 !! \kappa_h \leftarrow \min{\left( \kappa_{max}, \kappa_{cfl}, \max{\left( \kappa_{min}, \kappa_h \right)} \right)}
2168 !! f(c_g,z)
2169 !! \f]
2170 !!
2171 !! where \f$f(c_g,z)\f$ is a vertical structure function.
2172 !! \f$f(c_g,z)\f$ is calculated in module mom_lateral_mixing_coeffs.
2173 !! If <code>KHTH_USE_EBT_STRUCT=True</code> then \f$f(c_g,z)\f$ is set to look like the equivalent barotropic
2174 !! modal velocity structure. Otherwise \f$f(c_g,z)=1\f$ and the diffusivity is independent of depth.
2175 !!
2176 !! In order to calculate meaningful slopes in vanished layers, temporary copies of the thermodynamic variables
2177 !! are passed through a vertical smoother, function vert_fill_ts():
2178 !! \f{eqnarray*}{
2179 !! \left[ 1 + \Delta t \kappa_{smth} \frac{\partial^2}{\partial_z^2} \right] \theta & \leftarrow & \theta \\
2180 !! \left[ 1 + \Delta t \kappa_{smth} \frac{\partial^2}{\partial_z^2} \right] s & \leftarrow & s
2181 !! \f}
2182 !!
2183 !! \subsection section_khth_module_parameters Module mom_thickness_diffuse parameters
2184 !!
2185 !! | Symbol | Module parameter |
2186 !! | ------ | --------------- |
2187 !! | - | <code>THICKNESSDIFFUSE</code> |
2188 !! | \f$ \kappa_o \f$ | <code>KHTH</code> |
2189 !! | \f$ \alpha_{s} \f$ | <code>KHTH_SLOPE_CFF</code> |
2190 !! | \f$ \kappa_{min} \f$ | <code>KHTH_MIN</code> |
2191 !! | \f$ \kappa_{max} \f$ | <code>KHTH_MAX</code> |
2192 !! | - | <code>KHTH_MAX_CFL</code> |
2193 !! | \f$ \kappa_{smth} \f$ | <code>KD_SMOOTH</code> |
2194 !! | \f$ \alpha_{M} \f$ | <code>MEKE_KHTH_FAC</code> (from mom_meke module) |
2195 !! | - | <code>KHTH_USE_EBT_STRUCT</code> (from mom_lateral_mixing_coeffs module) |
2196 !! | - | <code>KHTH_USE_FGNV_STREAMFUNCTION</code> |
2197 !! | \f$ \gamma_F \f$ | <code>FGNV_FILTER_SCALE</code> |
2198 !! | \f$ c_{min} \f$ | <code>FGNV_C_MIN</code> |
2199 !!
2200 !! \subsection section_khth_module_reference References
2201 !!
2202 !! Ferrari, R., S.M. Griffies, A.J.G. Nurser and G.K. Vallis, 2010:
2203 !! A boundary-value problem for the parameterized mesoscale eddy transport.
2204 !! Ocean Modelling, 32, 143-156. http://doi.org/10.1016/j.ocemod.2010.01.004
2205 !!
2206 !! Viscbeck, M., J.C. Marshall, H. Jones, 1996:
2207 !! Dynamics of isolated convective regions in the ocean. J. Phys. Oceangr., 26, 1721-1734.
2208 !! http://dx.doi.org/10.1175/1520-0485(1996)026%3C1721:DOICRI%3E2.0.CO;2
2209 
2210 end module mom_thickness_diffuse
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:80
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:54
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_isopycnal_slopes
Calculations of isoneutral slopes and stratification.
Definition: MOM_isopycnal_slopes.F90:2
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_domains::pass_vector
Do a halo update on a pair of arrays representing the two components of a vector.
Definition: MOM_domains.F90:59
mom_thickness_diffuse::thickness_diffuse_cs
Control structure for thickness diffusion.
Definition: MOM_thickness_diffuse.F90:37
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_meke_types::meke_type
This type is used to exchange information related to the MEKE calculations.
Definition: MOM_MEKE_types.F90:8
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_lateral_mixing_coeffs::varmix_cs
Variable mixing coefficients.
Definition: MOM_lateral_mixing_coeffs.F90:27
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_eos::calculate_density_second_derivs
Calculates the second derivatives of density with various combinations of temperature,...
Definition: MOM_EOS.F90:93
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_variables::cont_diag_ptrs
Pointers to arrays with transports, which can later be used for derived diagnostics,...
Definition: MOM_variables.F90:193
mom_eos::calculate_density_derivs
Calculate the derivatives of density with temperature and salinity from T, S, and P.
Definition: MOM_EOS.F90:81
mom_interface_heights::find_eta
Calculates the heights of sruface or all interfaces from layer thicknesses.
Definition: MOM_interface_heights.F90:21
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_debugging
Provides checksumming functions for debugging.
Definition: MOM_debugging.F90:7
mom_thickness_diffuse
Thickness diffusion (or Gent McWilliams)
Definition: MOM_thickness_diffuse.F90:2
mom_lateral_mixing_coeffs
Variable mixing coefficients.
Definition: MOM_lateral_mixing_coeffs.F90:2
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_interface_heights
Functions for calculating interface heights, including free surface height.
Definition: MOM_interface_heights.F90:2
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
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:240
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:68