MOM_thickness_diffuse.F90

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