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