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