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