← back to index

src/parameterizations/lateral/MOM_hor_visc.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

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!> Calculates horizontal viscosity and viscous stresses
6module MOM_hor_visc
7
8use MOM_checksums, only : hchksum, Bchksum, uvchksum
9use MOM_coms, only : min_across_PEs
10use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
11use MOM_diag_mediator, only : post_product_u, post_product_sum_u
12use MOM_diag_mediator, only : post_product_v, post_product_sum_v
13use MOM_diag_mediator, only : diag_ctrl, time_type
14use MOM_domains, only : pass_var, CORNER, pass_vector, AGRID, BGRID_NE
15use MOM_domains, only : To_All, Scalar_Pair
16use MOM_error_handler, only : MOM_error, FATAL, WARNING, is_root_pe
17use MOM_file_parser, only : get_param, log_version, param_file_type
18use MOM_grid, only : ocean_grid_type
19use MOM_interface_heights, only : thickness_to_dz
20use MOM_lateral_mixing_coeffs, only : VarMix_CS, calc_QG_slopes, calc_QG_Leith_viscosity
21use MOM_barotropic, only : barotropic_CS, barotropic_get_tav
22use MOM_thickness_diffuse, only : thickness_diffuse_CS, thickness_diffuse_get_KH
23use MOM_io, only : MOM_read_data, slasher
24use MOM_MEKE_types, only : MEKE_type
25use MOM_open_boundary, only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W
26use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S
27use MOM_open_boundary, only : OBC_STRAIN_NONE, OBC_STRAIN_ZERO, OBC_STRAIN_FREESLIP
28use MOM_open_boundary, only : OBC_STRAIN_COMPUTED, OBC_STRAIN_SPECIFIED
29use MOM_stochastics, only : stochastic_CS
30use MOM_unit_scaling, only : unit_scale_type
31use MOM_verticalGrid, only : verticalGrid_type
32use MOM_variables, only : accel_diag_ptrs, thermo_var_ptrs
33use MOM_Zanna_Bolton, only : ZB2020_lateral_stress, ZB2020_init, ZB2020_end
34use MOM_Zanna_Bolton, only : ZB2020_CS, ZB2020_copy_gradient_and_thickness
35
36implicit none ; private
37
38#include <MOM_memory.h>
39
40public horizontal_viscosity, hor_visc_init, hor_visc_end, hor_visc_vel_stencil
41
42!> Control structure for horizontal viscosity
43type, public :: hor_visc_CS ; private
44 logical :: initialized = .false. !< True if this control structure has been initialized.
45 logical :: Laplacian !< Use a Laplacian horizontal viscosity if true.
46 logical :: biharmonic !< Use a biharmonic horizontal viscosity if true.
47 logical :: debug !< If true, write verbose checksums for debugging purposes.
48 logical :: no_slip !< If true, no slip boundary conditions are used.
49 !! Otherwise free slip boundary conditions are assumed.
50 !! The implementation of the free slip boundary
51 !! conditions on a C-grid is much cleaner than the
52 !! no slip boundary conditions. The use of free slip
53 !! b.c.s is strongly encouraged. The no slip b.c.s
54 !! are not implemented with the biharmonic viscosity.
55 logical :: bound_Kh !< If true, the Laplacian coefficient is locally
56 !! limited to guarantee stability.
57 logical :: EY24_EBT_BS !! If true, use an equivalent barotropic backscatter
58 !! with a stabilizing kill switch in MEKE,
59 !< developed by Yankovsky et al. 2024
60 logical :: bound_Ah !< If true, the biharmonic coefficient is locally
61 !! limited to guarantee stability.
62 real :: Re_Ah !! If nonzero, the biharmonic coefficient is scaled
63 !< so that the biharmonic Reynolds number is equal to this [nondim].
64 real :: bound_coef !< The nondimensional coefficient of the ratio of
65 !! the viscosity bounds to the theoretical maximum
66 !! for stability without considering other terms [nondim].
67 !! The default is 0.8.
68 real :: KS_coef !< A nondimensional coefficient on the biharmonic viscosity that sets the
69 !! kill switch for backscatter. Default is 1.0 [nondim].
70 real :: KS_timescale !< A timescale for computing CFL limit for turning off backscatter [T ~> s].
71 logical :: backscatter_underbound !< If true, the bounds on the biharmonic viscosity are allowed
72 !! to increase where the Laplacian viscosity is negative (due to
73 !! backscatter parameterizations) beyond the largest timestep-dependent
74 !! stable values of biharmonic viscosity when no Laplacian viscosity is
75 !! applied. The default is true for historical reasons, but this option
76 !! probably should not be used as it can lead to numerical instabilities.
77 logical :: Smagorinsky_Kh !< If true, use Smagorinsky nonlinear eddy
78 !! viscosity. KH is the background value.
79 logical :: Smagorinsky_Ah !< If true, use a biharmonic form of Smagorinsky
80 !! nonlinear eddy viscosity. AH is the background.
81 logical :: Leith_Kh !< If true, use 2D Leith nonlinear eddy
82 !! viscosity. KH is the background value.
83 logical :: Modified_Leith !< If true, use extra component of Leith viscosity
84 !! to damp divergent flow. To use, still set Leith_Kh=.TRUE.
85 logical :: use_beta_in_Leith !< If true, includes the beta term in the Leith viscosity
86 logical :: Leith_Ah !< If true, use a biharmonic form of 2D Leith
87 !! nonlinear eddy viscosity. AH is the background.
88 logical :: use_Leithy !< If true, use a biharmonic form of 2D Leith
89 !! nonlinear eddy viscosity with harmonic backscatter.
90 !! Ah is the background. Leithy = Leith+E
91 real :: c_K !< Fraction of energy dissipated by the biharmonic term
92 !! that gets backscattered in the Leith+E scheme. [nondim]
93 logical :: smooth_Ah !< If true (default), then Ah and m_leithy are smoothed.
94 !! This smoothing requires a lot of blocking communication.
95 logical :: use_QG_Leith_visc !< If true, use QG Leith nonlinear eddy viscosity.
96 !! KH is the background value.
97 logical :: bound_Coriolis !< If true & SMAGORINSKY_AH is used, the biharmonic
98 !! viscosity is modified to include a term that
99 !! scales quadratically with the velocity shears.
100 logical :: use_Kh_bg_2d !< Read 2d background viscosity from a file.
101 logical :: Kh_bg_2d_bug !< If true, retain an answer-changing horizontal indexing bug
102 !! in setting the corner-point viscosities when USE_KH_BG_2D=True.
103 real :: Kh_bg_min !< The minimum value allowed for Laplacian horizontal
104 !! viscosity [L2 T-1 ~> m2 s-1]. The default is 0.0.
105 logical :: FrictWork_bug !< If true, retain an answer-changing bug in calculating FrictWork,
106 !! which cancels the h in thickness flux and the h at velocity point.
107 logical :: OBC_strain_bug !< If true, recover a bug that specified shear strain option at open
108 !! boundaries cannot be applied.
109 logical :: use_land_mask !< Use the land mask for the computation of thicknesses
110 !! at velocity locations. This eliminates the dependence on
111 !! arbitrary values over land or outside of the domain.
112 !! Default is False to maintain answers with legacy experiments
113 !! but should be changed to True for new experiments.
114 logical :: anisotropic !< If true, allow anisotropic component to the viscosity.
115 logical :: add_LES_viscosity!< If true, adds the viscosity from Smagorinsky and Leith to
116 !! the background viscosity instead of taking the maximum.
117 real :: Kh_aniso !< The anisotropic viscosity [L2 T-1 ~> m2 s-1].
118 logical :: dynamic_aniso !< If true, the anisotropic viscosity is recomputed as a function
119 !! of state. This is set depending on ANISOTROPIC_MODE.
120 logical :: res_scale_MEKE !< If true, the viscosity contribution from MEKE is scaled by
121 !! the resolution function.
122 logical :: use_GME !< If true, use GME backscatter scheme.
123 integer :: answer_date !< The vintage of the order of arithmetic and expressions in the
124 !! horizontal viscosity calculations. Values below 20190101 recover
125 !! the answers from the end of 2018, while higher values use updated
126 !! and more robust forms of the same expressions.
127 real :: GME_h0 !< The strength of GME tapers quadratically to zero when the bathymetric
128 !! total water column thickness is less than GME_H0 [H ~> m or kg m-2]
129 real :: GME_efficiency !< The nondimensional prefactor multiplying the GME coefficient [nondim]
130 real :: GME_limiter !< The absolute maximum value the GME coefficient is allowed to take [L2 T-1 ~> m2 s-1].
131 real :: min_grid_Kh !< Minimum horizontal Laplacian viscosity used to
132 !! limit the grid Reynolds number [L2 T-1 ~> m2 s-1]
133 real :: min_grid_Ah !< Minimun horizontal biharmonic viscosity used to
134 !! limit grid Reynolds number [L4 T-1 ~> m4 s-1]
135 logical :: use_cont_thick !< If true, thickness at velocity points adopts h[uv] in BT_cont from continuity solver.
136 logical :: use_cont_thick_bug !< If true, retain an answer-changing bug for thickness at velocity points.
137 type(ZB2020_CS) :: ZB2020 !< Zanna-Bolton 2020 control structure.
138 logical :: use_ZB2020 !< If true, use Zanna-Bolton 2020 parameterization.
139 logical :: use_circulation !< If true, use circulation theorem to compute vorticity (for ZB20 or Leith)
140
141 real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: Kh_bg_xx
142 !< The background Laplacian viscosity at h points [L2 T-1 ~> m2 s-1].
143 !! The actual viscosity may be the larger of this
144 !! viscosity and the Smagorinsky and Leith viscosities.
145 real, allocatable :: Kh_bg_2d(:,:)
146 !< The background Laplacian viscosity at h points [L2 T-1 ~> m2 s-1].
147 !! The actual viscosity may be the larger of this
148 !! viscosity and the Smagorinsky and Leith viscosities.
149 real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: Ah_bg_xx
150 !< The background biharmonic viscosity at h points [L4 T-1 ~> m4 s-1].
151 !! The actual viscosity may be the larger of this
152 !! viscosity and the Smagorinsky and Leith viscosities.
153 real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: reduction_xx
154 !< The amount by which stresses through h points are reduced
155 !! due to partial barriers [nondim].
156 real, allocatable :: Kh_Max_xx(:,:) !< The maximum permitted Laplacian viscosity [L2 T-1 ~> m2 s-1].
157 real, allocatable :: Ah_Max_xx(:,:) !< The maximum permitted biharmonic viscosity [L4 T-1 ~> m4 s-1].
158 real, allocatable :: Ah_Max_xx_KS(:,:) !< The maximum permitted biharmonic viscosity for
159 !! the kill switch [L4 T-1 ~> m4 s-1].
160 real, allocatable :: n1n2_h(:,:) !< Factor n1*n2 in the anisotropic direction tensor at h-points [nondim]
161 real, allocatable :: n1n1_m_n2n2_h(:,:) !< Factor n1**2-n2**2 in the anisotropic direction tensor at h-points [nondim]
162 real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: &
163 grid_sp_h2, & !< Harmonic mean of the squares of the grid [L2 ~> m2]
164 grid_sp_h3 !< Harmonic mean of the squares of the grid^(3/2) [L3 ~> m3]
165 real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: Kh_bg_xy
166 !< The background Laplacian viscosity at q points [L2 T-1 ~> m2 s-1].
167 !! The actual viscosity may be the larger of this
168 !! viscosity and the Smagorinsky and Leith viscosities.
169 real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: Ah_bg_xy
170 !< The background biharmonic viscosity at q points [L4 T-1 ~> m4 s-1].
171 !! The actual viscosity may be the larger of this
172 !! viscosity and the Smagorinsky and Leith viscosities.
173 real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: reduction_xy
174 !< The amount by which stresses through q points are reduced
175 !! due to partial barriers [nondim].
176 real, allocatable :: Kh_Max_xy(:,:) !< The maximum permitted Laplacian viscosity [L2 T-1 ~> m2 s-1].
177 real, allocatable :: Ah_Max_xy(:,:) !< The maximum permitted biharmonic viscosity [L4 T-1 ~> m4 s-1].
178 real, allocatable :: Ah_Max_xy_KS(:,:) !< The maximum permitted biharmonic viscosity for
179 !! the kill switch [L4 T-1 ~> m4 s-1].
180 real, allocatable :: n1n2_q(:,:) !< Factor n1*n2 in the anisotropic direction tensor at q-points [nondim]
181 real, allocatable :: n1n1_m_n2n2_q(:,:) !< Factor n1**2-n2**2 in the anisotropic direction tensor at q-points [nondim]
182
183 real ALLOCABLE_, dimension(NIMEM_,NJMEM_) :: &
184 dx2h, & !< Pre-calculated dx^2 at h points [L2 ~> m2]
185 dy2h, & !< Pre-calculated dy^2 at h points [L2 ~> m2]
186 dx_dyT, & !< Pre-calculated dx/dy at h points [nondim]
187 dy_dxT !< Pre-calculated dy/dx at h points [nondim]
188 real, allocatable :: m_const_leithy(:,:) !< Pre-calculated .5*sqrt(c_K)*max{dx,dy} [L ~> m]
189 real, allocatable :: m_leithy_max(:,:) !< Pre-calculated 4./max(dx,dy)^2 at h points [L-2 ~> m-2]
190 real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEMB_PTR_) :: &
191 dx2q, & !< Pre-calculated dx^2 at q points [L2 ~> m2]
192 dy2q, & !< Pre-calculated dy^2 at q points [L2 ~> m2]
193 dx_dyBu, & !< Pre-calculated dx/dy at q points [nondim]
194 dy_dxBu !< Pre-calculated dy/dx at q points [nondim]
195 real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_) :: &
196 Idx2dyCu, & !< 1/(dx^2 dy) at u points [L-3 ~> m-3]
197 Idxdy2u !< 1/(dx dy^2) at u points [L-3 ~> m-3]
198 real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_) :: &
199 Idx2dyCv, & !< 1/(dx^2 dy) at v points [L-3 ~> m-3]
200 Idxdy2v !< 1/(dx dy^2) at v points [L-3 ~> m-3]
201
202 ! The following variables are precalculated time-invariant combinations of
203 ! parameters and metric terms.
204 real, allocatable :: Laplac2_const_xx(:,:) !< Laplacian metric-dependent constants [L2 ~> m2]
205 real, allocatable :: Biharm6_const_xx(:,:) !< Biharmonic metric-dependent constants [L6 ~> m6]
206 real, allocatable :: Laplac3_const_xx(:,:) !< Laplacian metric-dependent constants [L3 ~> m3]
207 real, allocatable :: Biharm_const_xx(:,:) !< Biharmonic metric-dependent constants [L4 ~> m4]
208 real, allocatable :: Biharm_const2_xx(:,:) !< Biharmonic metric-dependent constants [T L4 ~> s m4]
209 real, allocatable :: Re_Ah_const_xx(:,:) !< Biharmonic metric-dependent constants [L3 ~> m3]
210
211 real, allocatable :: Laplac2_const_xy(:,:) !< Laplacian metric-dependent constants [L2 ~> m2]
212 real, allocatable :: Biharm6_const_xy(:,:) !< Biharmonic metric-dependent constants [L6 ~> m6]
213 real, allocatable :: Laplac3_const_xy(:,:) !< Laplacian metric-dependent constants [L3 ~> m3]
214 real, allocatable :: Biharm_const_xy(:,:) !< Biharmonic metric-dependent constants [L4 ~> m4]
215 real, allocatable :: Biharm_const2_xy(:,:) !< Biharmonic metric-dependent constants [T L4 ~> s m4]
216 real, allocatable :: Re_Ah_const_xy(:,:) !< Biharmonic metric-dependent constants [L3 ~> m3]
217
218 type(diag_ctrl), pointer :: diag => NULL() !< structure to regulate diagnostics
219
220 ! real, allocatable :: hf_diffu(:,:,:) ! Zonal horizontal viscous acceleleration times
221 ! ! fractional thickness [L T-2 ~> m s-2].
222 ! real, allocatable :: hf_diffv(:,:,:) ! Meridional horizontal viscous acceleleration times
223 ! ! fractional thickness [L T-2 ~> m s-2].
224 ! 3D diagnostics hf_diffu(diffv) are commented because there is no clarity on proper remapping grid option.
225 ! The code is retained for debugging purposes in the future.
226
227 integer :: num_smooth_gme !< number of smoothing passes for the GME fluxes.
228 !>@{
229 !! Diagnostic id
230 integer :: id_grid_Re_Ah = -1, id_grid_Re_Kh = -1
231 integer :: id_diffu = -1, id_diffv = -1
232 ! integer :: id_hf_diffu = -1, id_hf_diffv = -1
233 integer :: id_h_diffu = -1, id_h_diffv = -1
234 integer :: id_hf_diffu_2d = -1, id_hf_diffv_2d = -1
235 integer :: id_intz_diffu_2d = -1, id_intz_diffv_2d = -1
236 integer :: id_diffu_visc_rem = -1, id_diffv_visc_rem = -1
237 integer :: id_Ah_h = -1, id_Ah_q = -1
238 integer :: id_Kh_h = -1, id_Kh_q = -1
239 integer :: id_GME_coeff_h = -1, id_GME_coeff_q = -1
240 integer :: id_dudx_bt = -1, id_dvdy_bt = -1
241 integer :: id_dudy_bt = -1, id_dvdx_bt = -1
242 integer :: id_vort_xy_q = -1, id_div_xx_h = -1
243 integer :: id_sh_xy_q = -1, id_sh_xx_h = -1
244 integer :: id_FrictWork = -1, id_FrictWorkIntz = -1
245 integer :: id_FrictWork_bh = -1, id_FrictWorkIntz_bh = -1
246 integer :: id_FrictWork_GME = -1
247 integer :: id_normstress = -1, id_shearstress = -1
248 integer :: id_visc_limit_h = -1, id_visc_limit_q = -1
249 integer :: id_visc_limit_h_flag = -1, id_visc_limit_q_flag = -1
250 integer :: id_visc_limit_h_frac = -1, id_visc_limit_q_frac = -1
251 integer :: id_BS_coeff_h = -1, id_BS_coeff_q = -1
252 !>@}
253
254end type hor_visc_CS
255
256contains
257
258!> Calculates the acceleration due to the horizontal viscosity.
259!!
260!! A combination of biharmonic and Laplacian forms can be used. The coefficient
261!! may either be a constant or a shear-dependent form. The biharmonic is
262!! determined by twice taking the divergence of an appropriately defined stress
263!! tensor. The Laplacian is determined by doing so once.
264!!
265!! To work, the following fields must be set outside of the usual
266!! is:ie range before this subroutine is called:
267!! u(is-2:ie+2,js-2:je+2)
268!! v(is-2:ie+2,js-2:je+2)
269!! h(is-1:ie+1,js-1:je+1) or up to h(is-2:ie+2,js-2:je+2) with some Leith options.
27050subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, GV, US, &
27150 CS, tv, dt, OBC, BT, TD, ADp, hu_cont, hv_cont, STOCH)
272 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
273 type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
274 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
275 intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
276 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
277 intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
278 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
279 intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2].
280 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
281 intent(in) :: uh !< The zonal volume transport [H L2 T-1 ~> m3 s-1].
282 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
283 intent(in) :: vh !< The meridional volume transport [H L2 T-1 ~> m3 s-1].
284 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
285 intent(out) :: diffu !< Zonal acceleration due to convergence of
286 !! along-coordinate stress tensor [L T-2 ~> m s-2]
287 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
288 intent(out) :: diffv !< Meridional acceleration due to convergence
289 !! of along-coordinate stress tensor [L T-2 ~> m s-2].
290 type(MEKE_type), intent(inout) :: MEKE !< MEKE fields
291 !! related to Mesoscale Eddy Kinetic Energy.
292 type(VarMix_CS), intent(inout) :: VarMix !< Variable mixing control structure
293 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
294 type(hor_visc_CS), intent(inout) :: CS !< Horizontal viscosity control structure
295 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
296 !! thermodynamic variables
297 real, intent(in) :: dt !< Time increment [T ~> s]
298 type(ocean_OBC_type), optional, pointer :: OBC !< Pointer to an open boundary condition type
299 type(barotropic_CS), optional, intent(in) :: BT !< Barotropic control structure
300 type(thickness_diffuse_CS), optional, intent(in) :: TD !< Thickness diffusion control structure
301 type(accel_diag_ptrs), optional, intent(in) :: ADp !< Acceleration diagnostics
302 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
303 optional, intent(inout) :: hu_cont !< Layer thickness at u-points [H ~> m or kg m-2].
304 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
305 optional, intent(inout) :: hv_cont !< Layer thickness at v-points [H ~> m or kg m-2].
306 type(stochastic_CS), intent(inout), optional :: STOCH !< Stochastic control structure
307
308 ! Local variables
309 real, dimension(SZIB_(G),SZJ_(G)) :: &
31050 Del2u, & ! The u-component of the Laplacian of velocity [L-1 T-1 ~> m-1 s-1]
31150 h_u, & ! Thickness interpolated to u points [H ~> m or kg m-2].
31250 vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
31350 vort_xy_dy_smooth, & ! y-derivative of smoothed vertical vorticity [L-1 T-1 ~> m-1 s-1]
31450 div_xx_dx, & ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
31550 ubtav ! zonal barotropic velocity averaged over a baroclinic time-step [L T-1 ~> m s-1]
316 real, dimension(SZI_(G),SZJB_(G)) :: &
31750 Del2v, & ! The v-component of the Laplacian of velocity [L-1 T-1 ~> m-1 s-1]
31850 h_v, & ! Thickness interpolated to v points [H ~> m or kg m-2].
31950 vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]
32050 vort_xy_dx_smooth, & ! x-derivative of smoothed vertical vorticity [L-1 T-1 ~> m-1 s-1]
32150 div_xx_dy, & ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]
32250 vbtav ! meridional barotropic velocity averaged over a baroclinic time-step [L T-1 ~> m s-1]
323 real, dimension(SZI_(G),SZJ_(G)) :: &
32450 dudx_bt, dvdy_bt, & ! components in the barotropic horizontal tension [T-1 ~> s-1]
32550 div_xx, & ! Estimate of horizontal divergence at h-points [T-1 ~> s-1]
32650 sh_xx, & ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1]
32750 sh_xx_smooth, & ! horizontal tension from smoothed velocity including metric terms [T-1 ~> s-1]
32850 sh_xx_bt, & ! barotropic horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1]
32950 str_xx,& ! str_xx is the diagonal term in the stress tensor [H L2 T-2 ~> m3 s-2 or kg s-2], but
330 ! at some points in the code it is not yet layer integrated, so is in [L2 T-2 ~> m2 s-2].
33150 str_xx_GME,& ! smoothed diagonal term in the stress tensor from GME [L2 T-2 ~> m2 s-2]
33250 bhstr_xx, & ! A copy of str_xx that only contains the biharmonic contribution [H L2 T-2 ~> m3 s-2 or kg s-2]
33350 FrictWorkIntz, & ! depth integrated energy dissipated by lateral friction [R Z L2 T-3 ~> W m-2]
33450 FrictWorkIntz_bh, & ! depth integrated energy dissipated by biharmonic lateral friction [R Z L2 T-3 ~> W m-2]
33550 grad_vort_mag_h, & ! Magnitude of vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1]
33650 grad_vort_mag_h_2d, & ! Magnitude of 2d vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1]
33750 grad_div_mag_h, & ! Magnitude of divergence gradient at h-points [L-1 T-1 ~> m-1 s-1]
33850 dudx, dvdy, & ! components in the horizontal tension [T-1 ~> s-1]
33950 dudx_smooth, dvdy_smooth, & ! components in the horizontal tension from smoothed velocity [T-1 ~> s-1]
34050 GME_effic_h, & ! The filtered efficiency of the GME terms at h points [nondim]
34150 m_leithy, & ! Kh=m_leithy*Ah in Leith+E parameterization [L-2 ~> m-2]
34250 Ah_sq, & ! The square of the biharmonic viscosity [L8 T-2 ~> m8 s-2]
34350 htot, & ! The total thickness of all layers [H ~> m or kg m-2]
34450 str_xx_BS ! The diagonal term in the stress tensor due to backscatter [H L2 T-2 ~> m3 s-2 or kg s-2]
345 real :: Del2vort_h ! Laplacian of vorticity at h-points [L-2 T-1 ~> m-2 s-1]
346 real :: grad_vel_mag_bt_h ! Magnitude of the barotropic velocity gradient tensor squared at h-points [T-2 ~> s-2]
347 real :: boundary_mask_h ! A mask that zeroes out cells with at least one land edge [nondim]
348
349 real, dimension(SZIB_(G),SZJB_(G)) :: &
35050 dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1]
35150 dvdx_smooth, dudy_smooth, & ! components in the shearing strain from smoothed velocity [T-1 ~> s-1]
35250 dDel2vdx, dDel2udy, & ! Components in the biharmonic equivalent of the shearing strain [L-2 T-1 ~> m-2 s-1]
35350 dvdx_bt, dudy_bt, & ! components in the barotropic shearing strain [T-1 ~> s-1]
35450 sh_xy, & ! horizontal shearing strain (du/dy + dv/dx) including metric terms [T-1 ~> s-1]
35550 sh_xy_smooth, & ! horizontal shearing strain from smoothed velocity including metric terms [T-1 ~> s-1]
35650 sh_xy_bt, & ! barotropic horizontal shearing strain (du/dy + dv/dx) inc. metric terms [T-1 ~> s-1]
35750 str_xy, & ! str_xy is the cross term in the stress tensor [H L2 T-2 ~> m3 s-2 or kg s-2], but
358 ! at some points in the code it is not yet layer integrated, so is in [L2 T-2 ~> m2 s-2].
35950 str_xy_GME, & ! smoothed cross term in the stress tensor from GME [L2 T-2 ~> m2 s-2]
36050 bhstr_xy, & ! A copy of str_xy that only contains the biharmonic contribution [H L2 T-2 ~> m3 s-2 or kg s-2]
36150 vort_xy, & ! Vertical vorticity (dv/dx - du/dy) including metric terms [T-1 ~> s-1]
36250 vort_xy_smooth, & ! Vertical vorticity including metric terms, smoothed [T-1 ~> s-1]
36350 grad_vort_mag_q, & ! Magnitude of vorticity gradient at q-points [L-1 T-1 ~> m-1 s-1]
36450 grad_vort_mag_q_2d, & ! Magnitude of 2d vorticity gradient at q-points [L-1 T-1 ~> m-1 s-1]
36550 Del2vort_q, & ! Laplacian of vorticity at q-points [L-2 T-1 ~> m-2 s-1]
36650 grad_div_mag_q, & ! Magnitude of divergence gradient at q-points [L-1 T-1 ~> m-1 s-1]
36750 hq, & ! harmonic mean of the harmonic means of the u- & v point thicknesses [H ~> m or kg m-2]
368 ! This form guarantees that hq/hu < 4.
36950 GME_effic_q, & ! The filtered efficiency of the GME terms at q points [nondim]
37050 str_xy_BS ! The cross term in the stress tensor due to backscatter [H L2 T-2 ~> m3 s-2 or kg s-2]
371 real :: grad_vel_mag_bt_q ! Magnitude of the barotropic velocity gradient tensor squared at q-points [T-2 ~> s-2]
372 real :: boundary_mask_q ! A mask that zeroes out cells with at least one land edge [nondim]
373
374 real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)) :: &
37550 Ah_q, & ! biharmonic viscosity at corner points [L4 T-1 ~> m4 s-1]
37650 Kh_q, & ! Laplacian viscosity at corner points [L2 T-1 ~> m2 s-1]
37750 vort_xy_q, & ! vertical vorticity at corner points [T-1 ~> s-1]
37850 sh_xy_q, & ! horizontal shearing strain at corner points [T-1 ~> s-1]
37950 GME_coeff_q, & !< GME coeff. at q-points [L2 T-1 ~> m2 s-1]
38050 visc_limit_q, & ! used to stabilize the EY24_EBT_BS backscatter [nondim]
38150 visc_limit_q_flag, & ! determines whether backscatter is shut off [nondim]
38250 visc_limit_q_frac, & ! determines how close backscatter is to shutting off [nondim]
38350 BS_coeff_q, & ! A diagnostic array of the backscatter coefficient [L2 T-1 ~> m2 s-1]
38450 ShSt ! A diagnostic array of shear stress [T-1 ~> s-1].
385 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: &
38650 KH_u_GME, & !< Isopycnal height diffusivities in u-columns [L2 T-1 ~> m2 s-1]
38750 slope_x !< Isopycnal slope in i-direction [Z L-1 ~> nondim]
388 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: &
38950 KH_v_GME, & !< Isopycnal height diffusivities in v-columns [L2 T-1 ~> m2 s-1]
39050 slope_y !< Isopycnal slope in j-direction [Z L-1 ~> nondim]
391 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &
39250 Ah_h, & ! biharmonic viscosity at thickness points [L4 T-1 ~> m4 s-1]
39350 Kh_h, & ! Laplacian viscosity at thickness points [L2 T-1 ~> m2 s-1]
39450 dz, & ! Height change across layers [Z ~> m]
39550 FrictWork, & ! work done by MKE dissipation mechanisms [R Z L2 T-3 ~> W m-2]
39650 FrictWork_bh, & ! work done by the biharmonic MKE dissipation mechanisms [R Z L2 T-3 ~> W m-2]
39750 FrictWork_GME, & ! work done by GME [R Z L2 T-3 ~> W m-2]
39850 div_xx_h, & ! horizontal divergence [T-1 ~> s-1]
39950 sh_xx_h, & ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1]
40050 NoSt, & ! A diagnostic array of normal stress [T-1 ~> s-1].
40150 BS_coeff_h ! A diagnostic array of the backscatter coefficient [L2 T-1 ~> m2 s-1]
402 real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
40350 grid_Re_Kh, & ! Grid Reynolds number for Laplacian horizontal viscosity at h points [nondim]
40450 grid_Re_Ah, & ! Grid Reynolds number for Biharmonic horizontal viscosity at h points [nondim]
40550 GME_coeff_h, & ! GME coefficient at h-points [L2 T-1 ~> m2 s-1]
40650 visc_limit_h, & ! Used to stabilize the EY24_EBT_BS backscatter [nondim]
40750 visc_limit_h_flag, & ! determines whether backscatter is shut off [nondim]
40850 visc_limit_h_frac ! determines how close backscatter is to shutting off [nondim]
409 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &
41050 u_smooth ! Zonal velocity, smoothed with a spatial low-pass filter [L T-1 ~> m s-1]
411 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: &
41250 v_smooth ! Meridional velocity, smoothed with a spatial low-pass filter [L T-1 ~> m s-1]
413 real :: AhSm ! Smagorinsky biharmonic viscosity [L4 T-1 ~> m4 s-1]
414 real :: AhLth ! 2D Leith biharmonic viscosity [L4 T-1 ~> m4 s-1]
415 real :: AhLthy ! 2D Leith+E biharmonic viscosity [L4 T-1 ~> m4 s-1]
416 real :: Shear_mag_bc ! Shear_mag value in backscatter [T-1 ~> s-1]
417 real :: sh_xx_sq ! Square of tension (sh_xx) [T-2 ~> s-2]
418 real :: sh_xy_sq ! Square of shearing strain (sh_xy) [T-2 ~> s-2]
419 real :: h2uq, h2vq ! temporary variables [H2 ~> m2 or kg2 m-4].
420 real :: hu, hv ! Thicknesses interpolated by arithmetic means to corner
421 ! points; these are first interpolated to u or v velocity
422 ! points where masks are applied [H ~> m or kg m-2].
423 real :: h_arith_q ! The arithmetic mean total thickness at q points [H ~> m or kg m-2]
424 real :: I_GME_h0 ! The inverse of GME tapering scale [H-1 ~> m-1 or m2 kg-1]
425 real :: h_neglect ! thickness so small it can be lost in roundoff and so neglected [H ~> m or kg m-2]
426 real :: h_neglect3 ! h_neglect^3 [H3 ~> m3 or kg3 m-6]
427 real :: h_min ! Minimum h at the 4 neighboring velocity points [H ~> m]
428 real :: Kh_max_here ! The local maximum Laplacian viscosity for stability [L2 T-1 ~> m2 s-1]
429 real :: RoScl ! The scaling function for MEKE source term [nondim]
430 real :: FatH ! abs(f) at h-point for MEKE source term [T-1 ~> s-1]
431 real :: local_strain ! Local variable for interpolating computed strain rates [T-1 ~> s-1].
432 real :: meke_res_fn ! A copy of the resolution scaling factor if being applied to MEKE [nondim]. Otherwise = 1.
433 real :: GME_coeff ! The GME (negative) viscosity coefficient [L2 T-1 ~> m2 s-1]
434 real :: DY_dxBu ! Ratio of meridional over zonal grid spacing at vertices [nondim]
435 real :: DX_dyBu ! Ratio of zonal over meridional grid spacing at vertices [nondim]
436 real :: Sh_F_pow ! The ratio of shear over the absolute value of f raised to some power and rescaled [nondim]
437 real :: backscat_subround ! The ratio of f over Shear_mag that is so small that the backscatter
438 ! calculation gives the same value as if f were 0 [nondim].
439 real :: KE ! Local kinetic energy [L2 T-2 ~> m2 s-2]
440 real :: d_del2u ! dy-weighted Laplacian(u) diff in x [L-2 T-1 ~> m-2 s-1]
441 real :: d_del2v ! dx-weighted Laplacian(v) diff in y [L-2 T-1 ~> m-2 s-1]
442 real :: d_str ! Stress tensor update [L2 T-2 ~> m2 s-2]
443 real :: grad_vort ! Vorticity gradient magnitude [L-1 T-1 ~> m-1 s-1]
444 real :: grad_vort_qg ! QG-based vorticity gradient magnitude [L-1 T-1 ~> m-1 s-1]
445 real :: grid_Kh ! Laplacian viscosity bound by grid [L2 T-1 ~> m2 s-1]
446 real :: grid_Ah ! Biharmonic viscosity bound by grid [L4 T-1 ~> m4 s-1]
447
448 logical :: rescale_Kh
449 logical :: find_FrictWork
450 logical :: find_FrictWork_bh
451 logical :: apply_OBC = .false.
452 logical :: apply_OBC_strain
453 logical :: use_MEKE_Ku
454 logical :: use_MEKE_Au
455 logical :: skeb_use_frict
456 logical :: use_cont_huv
457 logical :: use_kh_struct
458 integer :: is_vort, ie_vort, js_vort, je_vort ! Loop ranges for vorticity terms
459 integer :: is_Kh, ie_Kh, js_Kh, je_Kh ! Loop ranges for thickness point viscosities
460 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
461 integer :: i, j, k, n
462 real :: inv_PI3, inv_PI2, inv_PI6 ! Powers of the inverse of pi [nondim]
463 real :: tmp
464
465 ! Fields evaluated on active layers, used for constructing 3D stress fields
466 ! NOTE: The position of these declarations can impact performance, due to the
467 ! very large number of stack arrays in this function. Move with caution!
468 ! NOTE: Several of these are declared with the memory extent of q-points, but the
469 ! same arrays are also used at h-points to reduce the memory footprint of this
470 ! module, so they should never be used in halo point or checksum calls.
471 real, dimension(SZIB_(G),SZJB_(G)) :: &
47225 Ah, & ! biharmonic viscosity (h or q) [L4 T-1 ~> m4 s-1]
47350 Kh, & ! Laplacian viscosity (h or q) [L2 T-1 ~> m2 s-1]
47450 Kh_BS, & ! Laplacian antiviscosity [L2 T-1 ~> m2 s-1]
47550 Shear_mag, & ! magnitude of the shear (h or q) [T-1 ~> s-1]
47650 vert_vort_mag, & ! magnitude of the vertical vorticity gradient (h or q) [L-1 T-1 ~> m-1 s-1]
47750 vert_vort_mag_smooth, & ! magnitude of gradient of smoothed vertical vorticity (h or q) [L-1 T-1 ~> m-1 s-1]
47850 hrat_min, & ! h_min divided by the thickness at the stress point (h or q) [nondim]
47925 visc_bound_rem ! fraction of overall viscous bounds that remain to be applied (h or q) [nondim]
480
481 ! New variables: move these up once ready
482 logical :: use_Leith ! True if any Leith parameterizations are enabled
483 logical :: use_vort_xy ! True if vort_xy must be computed
484 logical :: use_Smag ! True if a Smagorinsky viscosity is enabled
485
48625 use_Leith = CS%Leith_Kh .or. CS%Leith_Ah .or. CS%use_Leithy
48725 use_vort_xy = use_Leith .or. CS%id_vort_xy_q > 0 .or. CS%use_ZB2020
48825 use_Smag = CS%Smagorinsky_Kh .or. CS%Smagorinsky_Ah
489
49025 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
49125 Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
492
49325 h_neglect = GV%H_subroundoff
494 !h_neglect3 = h_neglect**3
49525 h_neglect3 = h_neglect*h_neglect*h_neglect
49625 inv_PI3 = 1.0/((4.0*atan(1.0))**3)
49725 inv_PI2 = 1.0/((4.0*atan(1.0))**2)
49825 inv_PI6 = inv_PI3 * inv_PI3
499
50025 if (CS%EY24_EBT_BS) then
5010 visc_limit_h(:,:,:) = 0.
5020 visc_limit_q(:,:,:) = 0.
5030 visc_limit_h_flag(:,:,:) = 0.
5040 visc_limit_q_flag(:,:,:) = 0.
5050 visc_limit_h_frac(:,:,:) = 0.
5060 visc_limit_q_frac(:,:,:) = 0.
507 endif
508
50925 skeb_use_frict = .false.
51025 if (present(STOCH)) skeb_use_frict = STOCH%skeb_use_frict
511
512219325 m_leithy(:,:) = 0.0 ! Initialize
513
51425 if (present(OBC)) then ; if (associated(OBC)) then ; if (OBC%OBC_pe) then
5150 apply_OBC = OBC%Flather_u_BCs_exist_globally .or. OBC%Flather_v_BCs_exist_globally
5160 apply_OBC = .true.
517 endif ; endif ; endif
518
51925 apply_OBC_strain = .false.
52025 if (present(OBC)) then ; if (associated(OBC)) then
521 apply_OBC_strain = (OBC%strain_config /= OBC_STRAIN_NONE) &
5220 .and. ((.not. CS%OBC_strain_bug) .or. (OBC%strain_config /= OBC_STRAIN_SPECIFIED))
523 endif ; endif
524
52525 if (.not.CS%initialized) call MOM_error(FATAL, &
5260 "MOM_hor_visc: Module must be initialized before it is used.")
527
52825 if (.not.(CS%Laplacian .or. CS%biharmonic)) return
529
530 find_FrictWork = CS%id_FrictWork > 0 .or. CS%id_FrictWorkIntz > 0 &
53125 .or. allocated(MEKE%mom_src)
532 find_FrictWork_bh = CS%id_FrictWork_bh > 0 .or. CS%id_FrictWorkIntz_bh > 0 &
53325 .or. allocated(MEKE%mom_src_bh)
534
53525 use_kh_struct = allocated(VarMix%BS_struct)
53625 backscat_subround = 0.0
53725 if (find_FrictWork .and. allocated(MEKE%mom_src) .and. (MEKE%backscatter_Ro_c > 0.0) .and. &
538 (MEKE%backscatter_Ro_Pow /= 0.0)) &
5390 backscat_subround = (1.0e-16/MEKE%backscatter_Ro_c)**(1.0/MEKE%backscatter_Ro_Pow)
540
541 ! Toggle whether to use a Laplacian viscosity derived from MEKE
54225 use_MEKE_Ku = allocated(MEKE%Ku)
54325 use_MEKE_Au = allocated(MEKE%Au)
544
54525 use_cont_huv = CS%use_cont_thick .and. present(hu_cont) .and. present(hv_cont)
54625 if (use_cont_huv .and. .not.CS%use_cont_thick_bug) then
5470 call pass_vector(hu_cont, hv_cont, G%domain, To_All+Scalar_Pair, halo=2)
548 endif
549
55025 rescale_Kh = .false.
55125 if (VarMix%use_variable_mixing) then
55225 rescale_Kh = VarMix%Resoln_scaled_Kh
553 if ((rescale_Kh .or. CS%res_scale_MEKE) &
55425 .and. (.not. allocated(VarMix%Res_fn_h) .or. .not. allocated(VarMix%Res_fn_q))) &
555 call MOM_error(FATAL, "MOM_hor_visc: VarMix%Res_fn_h and VarMix%Res_fn_q "//&
5560 "both need to be associated with Resoln_scaled_Kh or RES_SCALE_MEKE_VISC.")
5570 elseif (CS%res_scale_MEKE) then
558 call MOM_error(FATAL, "MOM_hor_visc: VarMix needs to be associated if "//&
5590 "RES_SCALE_MEKE_VISC is True.")
560 endif
561
562 ! Set the halo sizes used for the thickness-point viscosities.
56325 if (CS%use_Leithy .or. CS%debug) then
5640 js_Kh = js-1 ; je_Kh = je+1 ; is_Kh = is-1 ; ie_Kh = ie+1
565 else
56625 js_Kh = Jsq ; je_Kh = je+1 ; is_Kh = Isq ; ie_Kh = ie+1
567 endif
568
569 ! Set the halo sizes used for the vorticity calculations.
57025 if ((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then
5710 js_vort = js_Kh-2 ; je_vort = Jeq+2 ; is_vort = is_Kh-2 ; ie_vort = Ieq+2
5720 if ((G%isc-G%isd < 3) .or. (G%isc-G%isd < 3)) call MOM_error(FATAL, &
5730 "The minimum halo size is 3 when a Leith viscosity is being used.")
574 else
57525 js_vort = js-2 ; je_vort = Jeq+1 ; is_vort = is-2 ; ie_vort = Ieq+1
576 endif
577
57825 if (CS%use_GME) then
579
580 ! Initialize diagnostic arrays with zeros
5810 GME_coeff_h(:,:,:) = 0.0
5820 GME_coeff_q(:,:,:) = 0.0
5830 str_xx_GME(:,:) = 0.0
5840 str_xy_GME(:,:) = 0.0
585
586 ! Get barotropic velocities and their gradients
5870 call barotropic_get_tav(BT, ubtav, vbtav, G, US)
588
5890 call pass_vector(ubtav, vbtav, G%Domain)
5900 call pass_var(h, G%domain, halo=2)
591
592 ! Calculate the barotropic horizontal tension
5930 do j=js-2,je+2 ; do i=is-2,ie+2
594 dudx_bt(i,j) = CS%DY_dxT(i,j)*((G%IdyCu(I,j) * ubtav(I,j)) - &
5950 (G%IdyCu(I-1,j) * ubtav(I-1,j)))
596 dvdy_bt(i,j) = CS%DX_dyT(i,j)*((G%IdxCv(i,J) * vbtav(i,J)) - &
5970 (G%IdxCv(i,J-1) * vbtav(i,J-1)))
598 enddo ; enddo
5990 do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
6000 sh_xx_bt(i,j) = dudx_bt(i,j) - dvdy_bt(i,j)
601 enddo ; enddo
602
603 ! Components for the barotropic shearing strain
6040 do J=Jsq-2,Jeq+2 ; do I=Isq-2,Ieq+2
605 dvdx_bt(I,J) = CS%DY_dxBu(I,J)*((vbtav(i+1,J)*G%IdyCv(i+1,J)) &
6060 - (vbtav(i,J)*G%IdyCv(i,J)))
607 dudy_bt(I,J) = CS%DX_dyBu(I,J)*((ubtav(I,j+1)*G%IdxCu(I,j+1)) &
6080 - (ubtav(I,j)*G%IdxCu(I,j)))
609 enddo ; enddo
610
6110 if (CS%no_slip) then
6120 do J=js-2,je+1 ; do I=is-2,ie+1
6130 sh_xy_bt(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_bt(I,J) + dudy_bt(I,J) )
614 enddo ; enddo
615 else
6160 do J=js-2,je+1 ; do I=is-2,ie+1
6170 sh_xy_bt(I,J) = G%mask2dBu(I,J) * ( dvdx_bt(I,J) + dudy_bt(I,J) )
618 enddo ; enddo
619 endif
620
6210 do j=js-2,je+2 ; do i=is-2,ie+2
6220 htot(i,j) = 0.0
623 enddo ; enddo
6240 do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2
6250 htot(i,j) = htot(i,j) + h(i,j,k)
626 enddo ; enddo ; enddo
627
6280 I_GME_h0 = 1.0 / CS%GME_h0
6290 do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
6300 boundary_mask_h = (G%mask2dCu(I,j) * G%mask2dCu(I-1,j)) * (G%mask2dCv(i,J) * G%mask2dCv(i,J-1))
631 grad_vel_mag_bt_h = G%mask2dT(I,J) * boundary_mask_h * (dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &
632 (0.25*((dvdx_bt(I,J)+dvdx_bt(I-1,J-1)) + (dvdx_bt(I,J-1)+dvdx_bt(I-1,J))))**2 + &
6330 (0.25*((dudy_bt(I,J)+dudy_bt(I-1,J-1)) + (dudy_bt(I,J-1)+dudy_bt(I-1,J))))**2)
634 ! Probably the following test could be simplified to
635 ! if (boundary_mask_h * G%mask2dT(I,J) > 0.0) then
6360 if (grad_vel_mag_bt_h > 0.0) then
6370 GME_effic_h(i,j) = CS%GME_efficiency * G%mask2dT(I,J) * (MIN(htot(i,j) * I_GME_h0, 1.0)**2)
638 else
6390 GME_effic_h(i,j) = 0.0
640 endif
641 enddo ; enddo
642
6430 do J=js-2,je+1 ; do I=is-2,ie+1
6440 boundary_mask_q = (G%mask2dCv(i,J) * G%mask2dCv(i+1,J)) * (G%mask2dCu(I,j) * G%mask2dCu(I,j+1))
645 grad_vel_mag_bt_q = G%mask2dBu(I,J) * boundary_mask_q * (dvdx_bt(I,J)**2 + dudy_bt(I,J)**2 + &
646 (0.25*((dudx_bt(i,j)+dudx_bt(i+1,j+1)) + (dudx_bt(i,j+1)+dudx_bt(i+1,j))))**2 + &
6470 (0.25*((dvdy_bt(i,j)+dvdy_bt(i+1,j+1)) + (dvdy_bt(i,j+1)+dvdy_bt(i+1,j))))**2)
648 ! Probably the following test could be simplified to
649 ! if (boundary_mask_q * G%mask2dBu(I,J) > 0.0) then
6500 if (grad_vel_mag_bt_q > 0.0) then
6510 h_arith_q = 0.25 * ((htot(i,j) + htot(i+1,j+1)) + (htot(i+1,j) + htot(i,j+1)))
6520 GME_effic_q(I,J) = CS%GME_efficiency * G%mask2dBu(I,J) * (MIN(h_arith_q * I_GME_h0, 1.0)**2)
653 else
6540 GME_effic_q(I,J) = 0.0
655 endif
656 enddo ; enddo
657
6580 call thickness_diffuse_get_KH(TD, KH_u_GME, KH_v_GME, G, GV)
659
6600 call pass_vector(KH_u_GME, KH_v_GME, G%domain, To_All+Scalar_Pair)
661
6620 if (CS%debug) &
6630 call uvchksum("GME KH[u,v]_GME", KH_u_GME, KH_v_GME, G%HI, haloshift=2, unscale=US%L_to_m**2*US%s_to_T)
664
665 endif ! use_GME
666
66725 if (CS%use_Leithy) then
668 ! Smooth the velocity. Right now it happens twice. In the future
669 ! one might make the number of smoothing cycles a user-specified parameter
6700 do k=1,nz
671 ! One call applies the filter twice
6720 u_smooth(:,:,k) = u(:,:,k)
6730 v_smooth(:,:,k) = v(:,:,k)
6740 call smooth_x9_uv(G, u_smooth(:,:,k), v_smooth(:,:,k), zero_land=.false.)
675 enddo
6760 call pass_vector(u_smooth, v_smooth, G%Domain)
677 endif
678
67925 if (CS%use_QG_Leith_visc .and. ((CS%Leith_Kh) .or. (CS%Leith_Ah))) then
6800 call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=2)
681 ! Calculate isopycnal slopes that will be used for some forms of viscosity.
6820 call calc_QG_slopes(h, tv, dt, G, GV, US, slope_x, slope_y, VarMix, OBC)
683 ! If the following halo update is added, the calculations in calc_QG_slopes could work on just
684 ! the computational domains, and some halo updates outside of this routine could be smaller.
685 ! call pass_vector(slope_x, slope_y, G%Domain, halo=2)
686 endif
687
688 !$omp target enter data map(alloc: dudx, dudy, dvdx, dvdy, sh_xx, sh_xy)
689 !$omp target enter data map(alloc: h_u, h_v, hq)
690 !$omp target enter data map(alloc: str_xx, str_xy)
691 !$omp target enter data map(alloc: Del2u, Del2v) if (CS%biharmonic)
692 !$omp target enter data map(alloc: dDel2vdx, dDel2udy) if (CS%biharmonic)
693 !$omp target enter data map(alloc: Shear_mag) if (use_Smag)
694 !$omp target enter data map(alloc: Kh) if (CS%Laplacian)
695 !$omp target enter data map(alloc: Ah) if (CS%biharmonic)
696 ! TODO: Only needed if FrictWork_bh is true, and currently only used on CPU,
697 ! but I do not yet see any benefit to breaking up the calculation.
698 !$omp target enter data map(alloc: bhstr_xx, bhstr_xy) if (CS%biharmonic)
699
700 !$omp target enter data map(alloc: hrat_min) &
701 !$omp if (CS%bound_Kh .or. CS%bound_Ah)
702 !$omp target enter data map(alloc: visc_bound_rem) &
703 !$omp if (CS%bound_Kh .or. CS%bound_Ah)
704 !$omp target enter data map(alloc: sh_xy_q) &
705 !$omp if (CS%id_sh_xy_q > 0)
706
7071900 do k=1,nz
708 ! The following are the forms of the horizontal tension and horizontal
709 ! shearing strain advocated by Smagorinsky (1993) and discussed in
710 ! Griffies and Hallberg (2000).
711
712 ! Calculate horizontal tension
713234375 do concurrent (j=Jsq-1:Jeq+2, i=Isq-1:Ieq+2)
714 dudx(i,j) = CS%DY_dxT(i,j)*((G%IdyCu(I,j) * u(I,j,k)) - &
71515114375 (G%IdyCu(I-1,j) * u(I-1,j,k)))
716 enddo
717
718234375 do concurrent (j=Jsq-1:Jeq+2, i=Isq-1:Ieq+2)
719 dvdy(i,j) = CS%DX_dyT(i,j)*((G%IdxCv(i,J) * v(i,J,k)) - &
72015114375 (G%IdxCv(i,J-1) * v(i,J-1,k)))
721 enddo
722
723234375 do concurrent (j=Jsq-1:Jeq+2, i=Isq-1:Ieq+2)
72415114375 sh_xx(i,j) = dudx(i,j) - dvdy(i,j)
725 enddo
726
727 ! Components for the shearing strain
7281875 do concurrent (J=js_vort:je_vort, I=is_vort:ie_vort)
72914529375 dvdx(I,J) = CS%DY_dxBu(I,J)*((v(i+1,J,k)*G%IdyCv(i+1,J)) - (v(i,J,k)*G%IdyCv(i,J)))
73014761875 dudy(I,J) = CS%DX_dyBu(I,J)*((u(I,j+1,k)*G%IdxCu(I,j+1)) - (u(I,j,k)*G%IdxCu(I,j)))
731 enddo
732
7331875 if (CS%use_Leithy) then
734 ! Calculate horizontal tension from smoothed velocity
7350 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
736 dudx_smooth(i,j) = CS%DY_dxT(i,j)*((G%IdyCu(I,j) * u_smooth(I,j,k)) - &
7370 (G%IdyCu(I-1,j) * u_smooth(I-1,j,k)))
738 dvdy_smooth(i,j) = CS%DX_dyT(i,j)*((G%IdxCv(i,J) * v_smooth(i,J,k)) - &
7390 (G%IdxCv(i,J-1) * v_smooth(i,J-1,k)))
7400 sh_xx_smooth(i,j) = dudx_smooth(i,j) - dvdy_smooth(i,j)
741 enddo ; enddo
742
743 ! Components for the shearing strain from smoothed velocity
7440 do J=js_Kh-1,je_Kh ; do I=is_Kh-1,ie_Kh
745 dvdx_smooth(I,J) = CS%DY_dxBu(I,J) * &
7460 ((v_smooth(i+1,J,k)*G%IdyCv(i+1,J)) - (v_smooth(i,J,k)*G%IdyCv(i,J)))
747 dudy_smooth(I,J) = CS%DX_dyBu(I,J) * &
7480 ((u_smooth(I,j+1,k)*G%IdxCu(I,j+1)) - (u_smooth(I,j,k)*G%IdxCu(I,j)))
749 enddo ; enddo
750 endif ! use Leith+E
751
7521875 if (CS%id_normstress > 0) then
753 !$omp target update from(sh_xx)
7540 do j=js,je ; do i=is,ie
7550 NoSt(i,j,k) = sh_xx(i,j)
756 enddo ; enddo
757 endif
758
759 ! Interpolate the thicknesses to velocity points.
760 ! The extra wide halos are to accommodate the cross-corner-point projections
761 ! in OBCs, which are not ordinarily be necessary, and might not be necessary
762 ! even with OBCs if the accelerations are zeroed at OBC points, in which
763 ! case the j-loop for h_u could collapse to j=js=1,je+1. -RWH
7641875 if (use_cont_huv) then
7650 do concurrent (j=js-2:je+2, I=Isq-1:Ieq+1)
7660 h_u(I,j) = hu_cont(I,j,k)
767 enddo
7680 do concurrent (J=Jsq-1:Jeq+1, i=is-2:ie+2)
7690 h_v(i,J) = hv_cont(i,J,k)
770 enddo
7711875 elseif (CS%use_land_mask) then
7720 do concurrent (j=js-2:je+2, I=is-2:Ieq+1)
7730 h_u(I,j) = 0.5 * (G%mask2dT(i,j)*h(i,j,k) + G%mask2dT(i+1,j)*h(i+1,j,k))
774 enddo
7750 do concurrent (J=js-2:Jeq+1, i=is-2:ie+2)
7760 h_v(i,J) = 0.5 * (G%mask2dT(i,j)*h(i,j,k) + G%mask2dT(i,j+1)*h(i,j+1,k))
777 enddo
778 else
779232500 do concurrent (j=js-2:je+2, I=is-2:Ieq+1)
78014992500 h_u(I,j) = 0.5 * (h(i,j,k) + h(i+1,j,k))
781 enddo
782234375 do concurrent (J=js-2:Jeq+1, i=is-2:ie+2)
78314881875 h_v(i,J) = 0.5 * (h(i,j,k) + h(i,j+1,k))
784 enddo
785 endif
786
787 ! Adjust contributions to shearing strain and interpolated values of
788 ! thicknesses on open boundaries.
7891875 if (apply_OBC) then
790 !$omp target update from(dvdx, dudy, h_u, h_v)
791 ! TODO: Reindent this later
7920 do n=1,OBC%number_of_segments
793
7940 J = OBC%segment(n)%HI%JsdB ; I = OBC%segment(n)%HI%IsdB
7950 if (apply_OBC_strain) then
7960 if (OBC%segment(n)%is_N_or_S .and. (J >= Js_vort) .and. (J <= Je_vort)) then
7970 do I = max(OBC%segment(n)%HI%IsdB,Is_vort), min(OBC%segment(n)%HI%IedB,Ie_vort)
7980 select case (OBC%strain_config)
799 case (OBC_STRAIN_ZERO)
8000 dvdx(I,J) = 0. ; dudy(I,J) = 0.
801 case (OBC_STRAIN_FREESLIP)
8020 dudy(I,J) = 0.
803 case (OBC_STRAIN_COMPUTED)
8040 if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
805 dudy(I,J) = 2.0*CS%DX_dyBu(I,J)* &
8060 (OBC%segment(n)%tangential_vel(I,J,k) - u(I,j,k))*G%IdxCu(I,j)
807 else
808 dudy(I,J) = 2.0*CS%DX_dyBu(I,J)* &
8090 (u(I,j+1,k) - OBC%segment(n)%tangential_vel(I,J,k))*G%IdxCu(I,j+1)
810 endif
811 case (OBC_STRAIN_SPECIFIED)
8120 if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
8130 dudy(I,J) = CS%DX_dyBu(I,J)*OBC%segment(n)%tangential_grad(I,J,k)*G%IdxCu(I,j)*G%dxBu(I,J)
814 else
8150 dudy(I,J) = CS%DX_dyBu(I,J)*OBC%segment(n)%tangential_grad(I,J,k)*G%IdxCu(I,j+1)*G%dxBu(I,J)
816 endif
817 end select
8180 if (CS%use_Leithy) then
8190 dvdx_smooth(I,J) = dvdx(I,J)
8200 dudy_smooth(I,J) = dudy(I,J)
821 endif
822 enddo
8230 elseif (OBC%segment(n)%is_E_or_W .and. (I >= is_vort) .and. (I <= ie_vort)) then
8240 do J = max(OBC%segment(n)%HI%JsdB,js_vort), min(OBC%segment(n)%HI%JedB,je_vort)
8250 select case (OBC%strain_config)
826 case (OBC_STRAIN_ZERO)
8270 dvdx(I,J) = 0. ; dudy(I,J) = 0.
828 case (OBC_STRAIN_FREESLIP)
8290 dvdx(I,J) = 0.
830 case (OBC_STRAIN_COMPUTED)
8310 if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
832 dvdx(I,J) = 2.0*CS%DY_dxBu(I,J)* &
8330 (OBC%segment(n)%tangential_vel(I,J,k) - v(i,J,k))*G%IdyCv(i,J)
834 else
835 dvdx(I,J) = 2.0*CS%DY_dxBu(I,J)* &
8360 (v(i+1,J,k) - OBC%segment(n)%tangential_vel(I,J,k))*G%IdyCv(i+1,J)
837 endif
838 case (OBC_STRAIN_SPECIFIED)
8390 if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
8400 dvdx(I,J) = CS%DY_dxBu(I,J)*OBC%segment(n)%tangential_grad(I,J,k)*G%IdyCv(i,J)*G%dxBu(I,J)
841 else
8420 dvdx(I,J) = CS%DY_dxBu(I,J)*OBC%segment(n)%tangential_grad(I,J,k)*G%IdyCv(i+1,J)*G%dxBu(I,J)
843 endif
844 end select
8450 if (CS%use_Leithy) then
8460 dvdx_smooth(I,J) = dvdx(I,J)
8470 dudy_smooth(I,J) = dudy(I,J)
848 endif
849 enddo
850 endif
851 endif
852
8530 if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
854 ! There are extra wide halos here to accommodate the cross-corner-point
855 ! OBC projections, but they might not be necessary if the accelerations
856 ! are always zeroed out at OBC points, in which case the i-loop below
857 ! becomes do i=is-1,ie+1. -RWH
8580 if ((J >= js-2) .and. (J <= Jeq+1)) then
8590 do i = max(is-2,OBC%segment(n)%HI%isd), min(ie+2,OBC%segment(n)%HI%ied)
8600 h_v(i,J) = h(i,j,k)
861 enddo
862 endif
8630 elseif (OBC%segment(n)%direction == OBC_DIRECTION_S) then
8640 if ((J >= js-2) .and. (J <= Jeq+1)) then
8650 do i = max(is-2,OBC%segment(n)%HI%isd), min(ie+2,OBC%segment(n)%HI%ied)
8660 h_v(i,J) = h(i,j+1,k)
867 enddo
868 endif
8690 elseif (OBC%segment(n)%direction == OBC_DIRECTION_E) then
8700 if ((I >= is-2) .and. (I <= Ieq+1)) then
8710 do j = max(js-2,OBC%segment(n)%HI%jsd), min(je+2,OBC%segment(n)%HI%jed)
8720 h_u(I,j) = h(i,j,k)
873 enddo
874 endif
8750 elseif (OBC%segment(n)%direction == OBC_DIRECTION_W) then
8760 if ((I >= is-2) .and. (I <= Ieq+1)) then
8770 do j = max(js-2,OBC%segment(n)%HI%jsd), min(je+2,OBC%segment(n)%HI%jed)
8780 h_u(I,j) = h(i+1,j,k)
879 enddo
880 endif
881 endif
882 enddo ; endif
883 ! Now project thicknesses across corner points on OBCs.
8841875 if (apply_OBC) then ; do n=1,OBC%number_of_segments
8850 J = OBC%segment(n)%HI%JsdB ; I = OBC%segment(n)%HI%IsdB
8860 if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
8870 if ((J >= js-2) .and. (J <= je)) then
8880 do I = max(is-2,OBC%segment(n)%HI%IsdB), min(Ieq+1,OBC%segment(n)%HI%IedB)
8890 h_u(I,j+1) = h_u(I,j)
890 enddo
891 endif
8920 elseif (OBC%segment(n)%direction == OBC_DIRECTION_S) then
8930 if ((J >= js-1) .and. (J <= je+1)) then
8940 do I = max(is-2,OBC%segment(n)%HI%isd), min(Ieq+1,OBC%segment(n)%HI%ied)
8950 h_u(I,j) = h_u(I,j+1)
896 enddo
897 endif
8980 elseif (OBC%segment(n)%direction == OBC_DIRECTION_E) then
8990 if ((I >= is-2) .and. (I <= ie)) then
9000 do J = max(js-2,OBC%segment(n)%HI%jsd), min(Jeq+1,OBC%segment(n)%HI%jed)
9010 h_v(i+1,J) = h_v(i,J)
902 enddo
903 endif
9040 elseif (OBC%segment(n)%direction == OBC_DIRECTION_W) then
9050 if ((I >= is-1) .and. (I <= ie+1)) then
9060 do J = max(js-2,OBC%segment(n)%HI%jsd), min(Jeq+1,OBC%segment(n)%HI%jed)
9070 h_v(i,J) = h_v(i+1,J)
908 enddo
909 endif
910 endif
911 enddo
912 ! TODO: Fix indentation
913 !$omp target update to(dvdx, dudy, h_u, h_v)
914 endif
915
916 ! Shearing strain (including no-slip boundary conditions at the 2-D land-sea mask).
917 ! dudy and dvdx include modifications at OBCs from above.
9181875 if (CS%no_slip) then
9190 do concurrent (J=js-2:Jeq+1, I=is-2:Ieq+1)
9200 sh_xy(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx(I,J) + dudy(I,J) )
921 enddo
922 else
923232500 do concurrent (J=js-2:Jeq+1, I=is-2:Ieq+1)
92414761875 sh_xy(I,J) = G%mask2dBu(I,J) * ( dvdx(I,J) + dudy(I,J) )
925 enddo
926 endif
927
9281875 if (CS%id_shearstress > 0) then
9290 do concurrent (J=js-2:Jeq+1, I=is-2:Ieq+1)
9300 ShSt(I,J,k) = sh_xy(I,J)
931 enddo
932 endif
933
9341875 if (CS%use_Leithy) then
935 ! Shearing strain (including no-slip boundary conditions at the 2-D land-sea mask).
936 ! dudy_smooth and dvdx_smooth do not (yet) include modifications at OBCs from above.
9370 if (CS%no_slip) then
9380 do J=js-1,Jeq ; do I=is-1,Ieq
9390 sh_xy_smooth(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_smooth(I,J) + dudy_smooth(I,J) )
940 enddo ; enddo
941 else
9420 do J=js-1,Jeq ; do I=is-1,Ieq
9430 sh_xy_smooth(I,J) = G%mask2dBu(I,J) * ( dvdx_smooth(I,J) + dudy_smooth(I,J) )
944 enddo ; enddo
945 endif
946 endif ! use Leith+E
947
948 ! Evaluate Del2u = x.Div(Grad u) and Del2v = y.Div( Grad u)
9491875 if (CS%biharmonic) then
950232500 do concurrent (j=js-1:Jeq+1, I=Isq-1:Ieq+1)
951 Del2u(I,j) = CS%Idx2dyCu(I,j) * ((CS%dx2q(I,J)*sh_xy(I,J)) - (CS%dx2q(I,J-1)*sh_xy(I,J-1))) + &
95214531250 CS%Idxdy2u(I,j) * ((CS%dy2h(i+1,j)*sh_xx(i+1,j)) - (CS%dy2h(i,j)*sh_xx(i,j)))
953 enddo
954
955230625 do concurrent (J=Jsq-1:Jeq+1, i=is-1:Ieq+1)
956 Del2v(i,J) = CS%Idxdy2v(i,J) * ((CS%dy2q(I,J)*sh_xy(I,J)) - (CS%dy2q(I-1,J)*sh_xy(I-1,J))) - &
95714641875 CS%Idx2dyCv(i,J) * ((CS%dx2h(i,j+1)*sh_xx(i,j+1)) - (CS%dx2h(i,j)*sh_xx(i,j)))
958 enddo
959
9601875 if (apply_OBC) then ; if (OBC%zero_biharmonic) then
961 !$omp target update from(Del2u, Del2v)
9620 do n=1,OBC%number_of_segments
9630 I = OBC%segment(n)%HI%IsdB ; J = OBC%segment(n)%HI%JsdB
9640 if (OBC%segment(n)%is_N_or_S .and. (J >= Jsq-1) .and. (J <= Jeq+1)) then
9650 do I=OBC%segment(n)%HI%isd,OBC%segment(n)%HI%ied
9660 Del2v(i,J) = 0.
967 enddo
9680 elseif (OBC%segment(n)%is_E_or_W .and. (I >= Isq-1) .and. (I <= Ieq+1)) then
9690 do j=OBC%segment(n)%HI%jsd,OBC%segment(n)%HI%jed
9700 Del2u(I,j) = 0.
971 enddo
972 endif
973 enddo
974 !$omp target update to(Del2u, Del2v)
975 endif ; endif
976 endif
977
978 ! Vorticity
979
980 ! NOTE: Keep Leith code on CPU for now, but moving it should be
981 ! straightforward.
982
9831875 if (use_vort_xy) then
984 !$omp target update from(dvdx, dudy)
9850 if (CS%no_slip) then
9860 do J=js_vort,je_vort ; do I=is_vort,ie_vort
9870 vort_xy(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx(I,J) - dudy(I,J) )
988 enddo ; enddo
989 else
9900 if (CS%use_circulation) then
9910 do J=js_vort,je_vort ; do I=is_vort,ie_vort
992 vort_xy(I,J) = G%mask2dBu(I,J) * G%IareaBu(I,J) * ( &
993 ((v(i+1,J,k)*G%dyCv(i+1,J)) - (v(i,J,k)*G%dyCv(i,J))) &
994 - ((u(I,j+1,k)*G%dxCu(I,j+1)) - (u(I,j,k)*G%dxCu(I,j))) &
9950 )
996 enddo ; enddo
997 else
9980 do J=js_vort,je_vort ; do I=is_vort,ie_vort
9990 vort_xy(I,J) = G%mask2dBu(I,J) * ( dvdx(I,J) - dudy(I,J) )
1000 enddo ; enddo
1001 endif
1002 endif
1003 endif
1004
10051875 if (CS%use_Leithy) then
10060 if (CS%no_slip) then
10070 do J=js_Kh-1,je_Kh ; do I=is_Kh-1,ie_Kh
10080 vort_xy_smooth(I,J) = (2.0-G%mask2dBu(I,J)) * ( dvdx_smooth(I,J) - dudy_smooth(I,J) )
1009 enddo ; enddo
1010 else
10110 do J=js_Kh-1,je_Kh ; do I=is_Kh-1,ie_Kh
10120 vort_xy_smooth(I,J) = G%mask2dBu(I,J) * ( dvdx_smooth(I,J) - dudy_smooth(I,J) )
1013 enddo ; enddo
1014 endif
1015 endif
1016
10171875 if (use_Leith) then
1018
1019 ! Vorticity gradient
10200 do J=js-2,je_Kh ; do i=is_Kh-1,ie_Kh+1
10210 DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J)
10220 vort_xy_dx(i,J) = DY_dxBu * ((vort_xy(I,J) * G%IdyCu(I,j)) - (vort_xy(I-1,J) * G%IdyCu(I-1,j)))
1023 enddo ; enddo
1024
10250 do j=js_Kh-1,je_Kh+1 ; do I=is-2,ie_Kh
10260 DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J)
10270 vort_xy_dy(I,j) = DX_dyBu * ((vort_xy(I,J) * G%IdxCv(i,J)) - (vort_xy(I,J-1) * G%IdxCv(i,J-1)))
1028 enddo ; enddo
1029
10300 if (CS%use_Leithy) then
1031 ! Gradient of smoothed vorticity
10320 do J=js_Kh-1,je_Kh ; do i=is_Kh,ie_Kh
10330 DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J)
1034 vort_xy_dx_smooth(i,J) = DY_dxBu * &
10350 ((vort_xy_smooth(I,J) * G%IdyCu(I,j)) - (vort_xy_smooth(I-1,J) * G%IdyCu(I-1,j)))
1036 enddo ; enddo
1037
10380 do j=js_Kh,je_Kh ; do I=is_Kh-1,ie_Kh
10390 DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J)
1040 vort_xy_dy_smooth(I,j) = DX_dyBu * &
10410 ((vort_xy_smooth(I,J) * G%IdxCv(i,J)) - (vort_xy_smooth(I,J-1) * G%IdxCv(i,J-1)))
1042 enddo ; enddo
1043 endif ! If Leithy
1044
1045 ! Laplacian of vorticity
1046 ! if (CS%Leith_Ah .or. CS%use_Leithy) then
10470 do J=js_Kh-1,je_Kh ; do I=is_Kh-1,ie_Kh
10480 DY_dxBu = G%dyBu(I,J) * G%IdxBu(I,J)
10490 DX_dyBu = G%dxBu(I,J) * G%IdyBu(I,J)
1050
1051 Del2vort_q(I,J) = DY_dxBu * ((vort_xy_dx(i+1,J) * G%IdyCv(i+1,J)) - (vort_xy_dx(i,J) * G%IdyCv(i,J))) + &
10520 DX_dyBu * ((vort_xy_dy(I,j+1) * G%IdyCu(I,j+1)) - (vort_xy_dy(I,j) * G%IdyCu(I,j)))
1053 enddo ; enddo
1054 ! endif
1055
10560 if (CS%modified_Leith) then
1057 !$omp target update from(dudx, dvdy)
1058
1059 ! Divergence
10600 do j=js_Kh-1,je_Kh+1 ; do i=is_Kh-1,ie_Kh+1
10610 div_xx(i,j) = dudx(i,j) + dvdy(i,j)
1062 enddo ; enddo
1063
1064 ! Divergence gradient
10650 do j=js-1,je+1 ; do I=is_Kh-1,ie_Kh
10660 div_xx_dx(I,j) = G%IdxCu(I,j)*(div_xx(i+1,j) - div_xx(i,j))
1067 enddo ; enddo
10680 do J=js_Kh-1,je_Kh ; do i=is-1,ie+1
10690 div_xx_dy(i,J) = G%IdyCv(i,J)*(div_xx(i,j+1) - div_xx(i,j))
1070 enddo ; enddo
1071
1072 ! Magnitude of divergence gradient
10730 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
1074 grad_div_mag_h(i,j) = sqrt(((0.5*(div_xx_dx(I,j) + div_xx_dx(I-1,j)))**2) + &
10750 ((0.5*(div_xx_dy(i,J) + div_xx_dy(i,J-1)))**2))
1076 enddo ; enddo
10770 do J=js-1,Jeq ; do I=is-1,Ieq
1078 grad_div_mag_q(I,J) = sqrt(((0.5*(div_xx_dx(I,j) + div_xx_dx(I,j+1)))**2) + &
10790 ((0.5*(div_xx_dy(i,J) + div_xx_dy(i+1,J)))**2))
1080 enddo ; enddo
1081
1082 else
1083
10840 do j=js-1,je+1 ; do I=is_Kh-1,ie_Kh
10850 div_xx_dx(I,j) = 0.0
1086 enddo ; enddo
10870 do J=js_Kh-1,je_Kh ; do i=is-1,ie+1
10880 div_xx_dy(i,J) = 0.0
1089 enddo ; enddo
10900 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
10910 grad_div_mag_h(i,j) = 0.0
1092 enddo ; enddo
10930 do J=js-1,Jeq ; do I=is-1,Ieq
10940 grad_div_mag_q(I,J) = 0.0
1095 enddo ; enddo
1096
1097 endif ! CS%modified_Leith
1098
1099 ! Add in beta for the Leith viscosity
1100 ! TODO: Move G%dF_dx, G%dF_dy to GPU
1101
11020 if (CS%use_beta_in_Leith) then
11030 do J=js-2,Jeq+1 ; do i=is-1,ie+1
11040 vort_xy_dx(i,J) = vort_xy_dx(i,J) + 0.5 * ( G%dF_dx(i,j) + G%dF_dx(i,j+1))
1105 enddo ; enddo
11060 do j=js-1,je+1 ; do I=is-2,Ieq+1
11070 vort_xy_dy(I,j) = vort_xy_dy(I,j) + 0.5 * ( G%dF_dy(i,j) + G%dF_dy(i+1,j))
1108 enddo ; enddo
1109 endif ! CS%use_beta_in_Leith
1110
11110 if (CS%use_QG_Leith_visc) then
1112
11130 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
1114 grad_vort_mag_h_2d(i,j) = SQRT(((0.5*(vort_xy_dx(i,J) + vort_xy_dx(i,J-1)))**2) + &
11150 ((0.5*(vort_xy_dy(I,j) + vort_xy_dy(I-1,j)))**2) )
1116 enddo ; enddo
11170 do J=js-1,Jeq ; do I=is-1,Ieq
1118 grad_vort_mag_q_2d(I,J) = SQRT(((0.5*(vort_xy_dx(i,J) + vort_xy_dx(i+1,J)))**2) + &
11190 ((0.5*(vort_xy_dy(I,j) + vort_xy_dy(I,j+1)))**2) )
1120 enddo ; enddo
1121
1122 ! This accumulates terms, some of which are in VarMix.
1123 call calc_QG_Leith_viscosity(VarMix, G, GV, US, h, dz, k, div_xx_dx, div_xx_dy, &
11240 slope_x, slope_y, vort_xy_dx, vort_xy_dy)
1125
1126 endif
1127
11280 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
1129 grad_vort_mag_h(i,j) = SQRT(((0.5*(vort_xy_dx(i,J) + vort_xy_dx(i,J-1)))**2) + &
11300 ((0.5*(vort_xy_dy(I,j) + vort_xy_dy(I-1,j)))**2) )
1131 enddo ; enddo
11320 do J=js-1,Jeq ; do I=is-1,Ieq
1133 grad_vort_mag_q(I,J) = SQRT(((0.5*(vort_xy_dx(i,J) + vort_xy_dx(i+1,J)))**2) + &
11340 ((0.5*(vort_xy_dy(I,j) + vort_xy_dy(I,j+1)))**2) )
1135 enddo ; enddo
1136
11370 if (CS%use_Leithy) then
11380 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
1139 vert_vort_mag_smooth(i,j) = SQRT(((0.5*(vort_xy_dx_smooth(i,J) + &
1140 vort_xy_dx_smooth(i,J-1)))**2) + &
1141 ((0.5*(vort_xy_dy_smooth(I,j) + &
11420 vort_xy_dy_smooth(I-1,j)))**2) )
1143 enddo ; enddo
1144 endif ! Leithy
1145
1146 endif ! CS%Leith_Kh
1147
11481875 if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then
11491875 do concurrent (j=js_Kh:je_Kh, i=is_Kh:ie_Kh)
115014182500 sh_xx_sq = sh_xx(i,j)**2
1151 sh_xy_sq = 0.25 * ( ((sh_xy(I-1,J-1)**2) + (sh_xy(I,J)**2)) &
115214182500 + ((sh_xy(I-1,J)**2) + (sh_xy(I,J-1)**2)) )
115314413125 Shear_mag(i,j) = sqrt(sh_xx_sq + sh_xy_sq)
1154 enddo
1155 endif
1156
11571875 if (CS%bound_Ah .or. CS%bound_Kh) then
11581875 do concurrent (j=js_Kh:je_Kh, i=is_Kh:ie_Kh)
115914182500 h_min = min(h_u(I,j), h_u(I-1,j), h_v(i,J), h_v(i,J-1))
116014413125 hrat_min(i,j) = min(1.0, h_min / (h(i,j,k) + h_neglect))
1161 enddo
1162 endif
1163
11641875 if (CS%Laplacian) then
1165 ! Determine the Laplacian viscosity at h points, using the
1166 ! largest value from several parameterizations. Also get
1167 ! the Laplacian component of str_xx.
1168
11691875 if ((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then
11700 if (CS%use_QG_Leith_visc) then
11710 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
11720 grad_vort = grad_vort_mag_h(i,j) + grad_div_mag_h(i,j)
11730 grad_vort_qg = 3. * grad_vort_mag_h_2d(i,j)
11740 vert_vort_mag(i,j) = min(grad_vort, grad_vort_qg)
1175 enddo ; enddo
1176 else
11770 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
11780 vert_vort_mag(i,j) = grad_vort_mag_h(i,j) + grad_div_mag_h(i,j)
1179 enddo ; enddo
1180 endif
1181 endif
1182
1183 ! Static (pre-computed) background viscosity
11841875 do concurrent (j=js_Kh:je_Kh, i=is_Kh:ie_Kh)
118514413125 Kh(i,j) = CS%Kh_bg_xx(i,j)
1186 enddo
1187
11881875 if (CS%add_LES_viscosity) then
11890 if (CS%Smagorinsky_Kh) then
11900 do concurrent (j=js_Kh:je_Kh, i=is_Kh:ie_Kh)
11910 Kh(i,j) = Kh(i,j) + CS%Laplac2_const_xx(i,j) * Shear_mag(i,j)
1192 enddo
1193 endif
1194
11950 if (CS%Leith_Kh) then
11960 do concurrent (j=js_Kh:je_Kh, i=is_Kh:ie_Kh)
1197 Kh(i,j) = Kh(i,j) &
11980 + CS%Laplac3_const_xx(i,j) * vert_vort_mag(i,j) * inv_PI3
1199 enddo
1200 endif
1201 else
12021875 if (CS%Smagorinsky_Kh) then
12030 do concurrent (j=js_Kh:je_Kh, i=is_Kh:ie_Kh)
12040 Kh(i,j) = max(Kh(i,j), CS%Laplac2_const_xx(i,j) * Shear_mag(i,j))
1205 enddo
1206 endif
1207
12081875 if (CS%Leith_Kh) then
12090 do concurrent (j=js_Kh:je_Kh, i=is_Kh:ie_Kh)
1210 Kh(i,j) = max(Kh(i,j), &
12110 CS%Laplac3_const_xx(i,j) * vert_vort_mag(i,j) * inv_PI3)
1212 enddo
1213 endif
1214 endif
1215
1216 ! All viscosity contributions above are subject to resolution scaling
1217
12181875 if (rescale_Kh) then
1219 !$omp target update from(Kh)
122014300625 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
122114298750 Kh(i,j) = VarMix%Res_fn_h(i,j) * Kh(i,j)
1222 enddo ; enddo
1223 !$omp target update to(Kh)
1224 endif
1225
1226 ! Place a floor on the viscosity, if desired.
12271875 do concurrent (j=js_Kh:je_Kh, i=is_Kh:ie_Kh)
122814413125 Kh(i,j) = max(Kh(i,j), CS%Kh_bg_min)
1229 enddo
1230
12311875 if (use_MEKE_Ku .and. .not. CS%EY24_EBT_BS) then
1232 !$omp target update from(Kh)
1233 ! *Add* the MEKE contribution (which might be negative)
12340 if (use_kh_struct) then
12350 if (CS%res_scale_MEKE) then
12360 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
12370 Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%Res_fn_h(i,j) * VarMix%BS_struct(i,j,k)
1238 enddo ; enddo
1239 else
12400 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
12410 Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k)
1242 enddo ; enddo
1243 endif
1244 else
12450 if (CS%res_scale_MEKE) then
12460 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
12470 Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j) * VarMix%Res_fn_h(i,j)
1248 enddo ; enddo
1249 else
12500 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
12510 Kh(i,j) = Kh(i,j) + MEKE%Ku(i,j)
1252 enddo ; enddo
1253 endif
1254 endif
1255 !$omp target update to(Kh)
1256 endif
1257
12581875 if (CS%anisotropic) then
1259 !$omp target update from(Kh)
12600 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
1261 ! *Add* the tension component of anisotropic viscosity
12620 Kh(i,j) = Kh(i,j) + CS%Kh_aniso * (1. - CS%n1n2_h(i,j)**2)
1263 enddo ; enddo
1264 !$omp target update to(Kh)
1265 endif
1266
1267 !$omp target update to(Kh) &
1268 !$omp if ((use_MEKE_Ku .and. .not. CS%EY24_EBT_BS) .or. CS%anisotropic)
1269
1270 ! Newer method of bounding for stability
12711875 if ((CS%bound_Kh) .and. (CS%bound_Ah)) then
12721875 do concurrent (j=js_Kh:je_Kh, i=is_Kh:ie_Kh)
127314182500 visc_bound_rem(i,j) = 1.0
127414182500 Kh_max_here = hrat_min(i,j) * CS%Kh_Max_xx(i,j)
127514413125 if (Kh(i,j) >= Kh_max_here) then
12764713750 visc_bound_rem(i,j) = 0.0
12774713750 Kh(i,j) = Kh_max_here
12789468750 elseif ((Kh(i,j) > 0.0) .or. (CS%backscatter_underbound .and. (Kh_max_here > 0.0))) then
12790 visc_bound_rem(i,j) = 1.0 - Kh(i,j) / Kh_max_here
1280 endif
1281 enddo
12820 elseif (CS%bound_Kh) then
12830 do concurrent (j=js_Kh:je_Kh, i=is_Kh:ie_Kh)
12840 Kh(i,j) = min(Kh(i,j), hrat_min(i,j) * CS%Kh_Max_xx(i,j))
1285 enddo
1286 endif
1287
1288 ! In Leith+E parameterization Kh is computed after Ah in the biharmonic loop.
1289 ! The harmonic component of str_xx is added in the biharmonic loop.
12901875 if (CS%use_Leithy) then
1291 !$omp target update from(Kh)
12920 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
12930 Kh(i,j) = 0.
1294 enddo ; enddo
1295 !$omp target update to(Kh)
1296 endif
1297
12981875 if (CS%id_Kh_h>0 .or. CS%debug) then
1299 !$omp target update from(Kh)
13000 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
13010 Kh_h(i,j,k) = Kh(i,j)
1302 enddo ; enddo
1303 endif
1304
13051875 if (CS%id_grid_Re_Kh>0) then
1306 !$omp target update from(Kh)
13070 do j=js,je ; do i=is,ie
13080 KE = 0.125*(((u(I,j,k)+u(I-1,j,k))**2) + ((v(i,J,k)+v(i,J-1,k))**2))
13090 grid_Kh = max(Kh(i,j), CS%min_grid_Kh)
13100 grid_Re_Kh(i,j,k) = (sqrt(KE) * sqrt(CS%grid_sp_h2(i,j))) / grid_Kh
1311 enddo ; enddo
1312 endif
1313
13141875 if (CS%id_div_xx_h>0) then
1315 !$omp target update from(dudx, dvdy)
13160 do j=js,je ; do i=is,ie
13170 div_xx_h(i,j,k) = dudx(i,j) + dvdy(i,j)
1318 enddo ; enddo
1319 endif
1320
13211875 if (CS%id_sh_xx_h>0) then
1322 !$omp target update from(sh_xx)
13230 do j=js,je ; do i=is,ie
13240 sh_xx_h(i,j,k) = sh_xx(i,j)
1325 enddo ; enddo
1326 endif
1327
1328230625 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
132914413125 str_xx(i,j) = -Kh(i,j) * sh_xx(i,j)
1330 enddo
1331 else
13320 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
13330 str_xx(i,j) = 0.0
1334 enddo
1335 endif ! Get Kh at h points and get Laplacian component of str_xx
1336
13371875 if (CS%anisotropic) then
1338 !$omp target update from(str_xx, sh_xy)
13390 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
1340 ! Shearing-strain averaged to h-points
13410 local_strain = 0.25 * ( (sh_xy(I,J) + sh_xy(I-1,J-1)) + (sh_xy(I-1,J) + sh_xy(I,J-1)) )
1342 ! *Add* the shear-strain contribution to the xx-component of stress
13430 str_xx(i,j) = str_xx(i,j) - CS%Kh_aniso * CS%n1n2_h(i,j) * CS%n1n1_m_n2n2_h(i,j) * local_strain
1344 enddo ; enddo
1345 !$omp target update to(str_xx)
1346 endif
1347
13481875 if (CS%biharmonic) then
1349 ! Determine the biharmonic viscosity at h points, using the
1350 ! largest value from several parameterizations. Also get the
1351 ! biharmonic component of str_xx.
13521875 do concurrent (j=js_Kh:je_Kh, i=is_Kh:ie_Kh)
135314413125 Ah(i,j) = CS%Ah_bg_xx(i,j)
1354 enddo
1355
13561875 if ((CS%Smagorinsky_Ah) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then
13571875 if (CS%Smagorinsky_Ah) then
13581875 if (CS%bound_Coriolis) then
13591875 do concurrent (j=js_Kh:je_Kh, i=is_Kh:ie_Kh)
1360 AhSm = Shear_mag(i,j) * (CS%Biharm_const_xx(i,j) &
136114182500 + CS%Biharm_const2_xx(i,j) * Shear_mag(i,j))
136214413125 Ah(i,j) = max(Ah(i,j), AhSm)
1363 enddo
1364 else
13650 do concurrent (j=js_Kh:je_Kh, i=is_Kh:ie_Kh)
13660 AhSm = CS%Biharm_const_xx(i,j) * Shear_mag(i,j)
13670 Ah(i,j) = max(Ah(i,j), AhSm)
1368 enddo
1369 endif
1370 endif
1371
13721875 if (CS%Leith_Ah) then
1373 !$omp target update from(Ah)
13740 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
1375 Del2vort_h = 0.25 * ((Del2vort_q(I,J) + Del2vort_q(I-1,J-1)) + &
13760 (Del2vort_q(I-1,J) + Del2vort_q(I,J-1)))
13770 AhLth = CS%Biharm6_const_xx(i,j) * abs(Del2vort_h) * inv_PI6
13780 Ah(i,j) = max(Ah(i,j), AhLth)
1379 enddo ; enddo
1380 !$omp target update to(Ah)
1381 endif
1382
13831875 if (CS%use_Leithy) then
1384 ! TODO: !$omp target update from(...?)
1385
1386 ! Get m_leithy
13870 if (CS%smooth_Ah) m_leithy(:,:) = 0.0 ! This is here to initialize domain edge halo values.
13880 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
1389 Del2vort_h = 0.25 * ((Del2vort_q(I,J) + Del2vort_q(I-1,J-1)) + &
13900 (Del2vort_q(I-1,J) + Del2vort_q(I,J-1)))
13910 AhLth = CS%Biharm6_const_xx(i,j) * inv_PI6 * abs(Del2vort_h)
13920 if (AhLth <= CS%Ah_bg_xx(i,j)) then
13930 m_leithy(i,j) = 0.0
1394 else
13950 if ((CS%m_const_leithy(i,j)*vert_vort_mag(i,j)) < abs(vort_xy_smooth(i,j))) then
13960 m_leithy(i,j) = CS%c_K * (vert_vort_mag(i,j) / vort_xy_smooth(i,j))**2
1397 else
13980 m_leithy(i,j) = CS%m_leithy_max(i,j)
1399 endif
14000 m_leithy(i,j) = G%mask2dBu(i,j) * m_leithy(i,j)
1401 endif
1402 enddo ; enddo
1403
14040 if (CS%smooth_Ah) then
1405 ! Smooth m_leithy. A single call smoothes twice.
14060 call pass_var(m_leithy, G%Domain, halo=2)
14070 call smooth_x9_h(G, m_leithy, zero_land=.true.)
14080 call pass_var(m_leithy, G%Domain)
1409 endif
1410 ! Get Ah
14110 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
1412 Del2vort_h = 0.25 * ((Del2vort_q(I,J) + Del2vort_q(I-1,J-1)) + &
14130 (Del2vort_q(I-1,J) + Del2vort_q(I,J-1)))
1414 AhLthy = CS%Biharm6_const_xx(i,j) * inv_PI6 * &
14150 sqrt(max(0.,Del2vort_h**2 - m_leithy(i,j)*vert_vort_mag_smooth(i,j)**2))
14160 Ah(i,j) = max(CS%Ah_bg_xx(i,j), AhLthy)
1417 enddo ; enddo
14180 if (CS%smooth_Ah) then
1419 ! Smooth Ah before applying upper bound. Square Ah, then smooth, then take its square root.
14200 Ah_sq(:,:) = 0.0 ! This is here to initialize domain edge halo values.
14210 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
14220 Ah_sq(i,j) = Ah(i,j)**2
1423 enddo ; enddo
14240 call pass_var(Ah_sq, G%Domain, halo=2)
1425 ! A single call smoothes twice.
14260 call smooth_x9_h(G, Ah_sq, zero_land=.false.)
14270 call pass_var(Ah_sq, G%Domain)
14280 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
14290 Ah_h(i,j,k) = max(CS%Ah_bg_xx(i,j), sqrt(max(0., Ah_sq(i,j))))
14300 Ah(i,j) = Ah_h(i,j,k)
1431 enddo ; enddo
1432 else
14330 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
14340 Ah_h(i,j,k) = Ah(i,j)
1435 enddo ; enddo
1436 endif
1437 endif
1438 endif ! Smagorinsky_Ah or Leith_Ah or Leith+E
1439
14401875 if (use_MEKE_Au) then
1441 ! *Add* the MEKE contribution
1442 !$omp target update from(Ah)
14430 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
14440 Ah(i,j) = Ah(i,j) + MEKE%Au(i,j)
1445 enddo ; enddo
1446 !$omp target update to(Ah)
1447 endif
1448
14491875 if (CS%Re_Ah > 0.0) then
1450 !$omp target update from(Ah)
14510 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
14520 KE = 0.125*(((u(I,j,k)+u(I-1,j,k))**2) + ((v(i,J,k)+v(i,J-1,k))**2))
14530 Ah(i,j) = sqrt(KE) * CS%Re_Ah_const_xx(i,j)
1454 enddo ; enddo
1455 !$omp target update to(Ah)
1456 endif
1457
14581875 if (CS%bound_Ah) then
14591875 if (CS%bound_Kh) then
14601875 do concurrent (j=js_Kh:je_Kh, i=is_Kh:ie_Kh)
146114413125 Ah(i,j) = min(Ah(i,j), visc_bound_rem(i,j) * hrat_min(i,j) * CS%Ah_Max_xx(i,j))
1462 enddo
1463 else
14640 do concurrent (j=js_Kh:je_Kh, i=is_Kh:ie_Kh)
14650 Ah(i,j) = min(Ah(i,j), hrat_min(i,j) * CS%Ah_Max_xx(i,j))
1466 enddo
1467 endif
1468 endif
1469
14701875 if (CS%EY24_EBT_BS) then
1471 !$omp target update from(Ah)
14720 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
14730 tmp = CS%KS_coef * hrat_min(i,j) * CS%Ah_Max_xx_KS(i,j)
14740 visc_limit_h(i,j,k) = tmp
14750 visc_limit_h_frac(i,j,k) = Ah(i,j) / (CS%KS_coef * hrat_min(i,j) * CS%Ah_Max_xx_KS(i,j))
14760 if (Ah(i,j) >= tmp) then
14770 visc_limit_h_flag(i,j,k) = 1.
1478 endif
1479 enddo ; enddo
1480 endif
1481
14821875 if ((CS%id_Ah_h>0) .or. CS%debug .or. CS%use_Leithy) then
1483 !$omp target update from(Ah)
14840 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
14850 Ah_h(i,j,k) = Ah(i,j)
1486 enddo ; enddo
1487 endif
1488
14891875 if (CS%use_Leithy) then
1490 ! Compute Leith+E Kh after bounds have been applied to Ah
1491 ! and after it has been smoothed. Kh = -m_leithy * Ah
1492 !$omp target update from(Ah, Kh)
14930 do j=js_Kh,je_Kh ; do i=is_Kh,ie_Kh
14940 Kh(i,j) = -m_leithy(i,j) * Ah(i,j)
14950 Kh_h(i,j,k) = Kh(i,j)
1496 enddo ; enddo
1497 !$omp target update to(Kh)
1498 endif
1499
15001875 if (CS%id_grid_Re_Ah > 0) then
1501 !$omp target update from(Ah)
15020 do j=js,je ; do i=is,ie
15030 KE = 0.125 * (((u(I,j,k) + u(I-1,j,k))**2) + ((v(i,J,k) + v(i,J-1,k))**2))
15040 grid_Ah = max(Ah(i,j), CS%min_grid_Ah)
15050 grid_Re_Ah(i,j,k) = (sqrt(KE) * CS%grid_sp_h3(i,j)) / grid_Ah
1506 enddo ; enddo
1507 endif
1508
1509230625 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
151014182500 d_del2u = (G%IdyCu(I,j) * Del2u(I,j)) - (G%IdyCu(I-1,j) * Del2u(I-1,j))
151114182500 d_del2v = (G%IdxCv(i,J) * Del2v(i,J)) - (G%IdxCv(i,J-1) * Del2v(i,J-1))
151214182500 d_str = Ah(i,j) * ((CS%DY_dxT(i,j) * d_del2u) - (CS%DX_dyT(i,j) * d_del2v))
1513
151414182500 str_xx(i,j) = str_xx(i,j) + d_str
1515
1516 ! Keep a copy of the biharmonic contribution for backscatter parameterization
1517 ! XXX: Need to get out of the loop somehow
151814182500 if (find_FrictWork_bh) &
151914413125 bhstr_xx(i,j) = d_str * (h(i,j,k) * CS%reduction_xx(i,j))
1520 enddo
1521
15221875 if (CS%use_Leithy) then
1523 !$omp target update from(Kh)
15240 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
15250 str_xx(i,j) = str_xx(i,j) - Kh(i,j) * sh_xx_smooth(i,j)
1526 enddo ; enddo
1527 endif
1528 endif ! Get biharmonic coefficient at h points and biharmonic part of str_xx
1529
1530 ! Backscatter using MEKE
15311875 if (CS%EY24_EBT_BS) then
1532 !$omp target update from(sh_xx)
15330 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
15340 if (visc_limit_h_flag(i,j,k) > 0) then
15350 Kh_BS(i,j) = 0.
1536 else
15370 if (use_kh_struct) then
15380 Kh_BS(i,j) = MEKE%Ku(i,j) * VarMix%BS_struct(i,j,k)
1539 else
15400 Kh_BS(i,j) = MEKE%Ku(i,j)
1541 endif
1542 endif
1543 enddo ; enddo
1544
15450 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
15460 str_xx_BS(i,j) = -Kh_BS(i,j) * sh_xx(i,j)
1547 enddo ; enddo
1548
15490 if (CS%id_BS_coeff_h>0) then
15500 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
15510 BS_coeff_h(i,j,k) = Kh_BS(i,j)
1552 enddo ; enddo
1553 endif
1554
15550 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
15560 str_xx(i,j) = str_xx(i,j) + str_xx_BS(i,j)
1557 enddo ; enddo
1558 !$omp target update to(str_xx)
1559 endif ! Backscatter
1560
15611875 if (CS%biharmonic) then
1562 ! Gradient of Laplacian, for use in bi-harmonic term
15631875 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
156413839375 dDel2vdx(I,J) = CS%DY_dxBu(I,J)*((Del2v(i+1,J)*G%IdyCv(i+1,J)) - (Del2v(i,J)*G%IdyCv(i,J)))
156514068125 dDel2udy(I,J) = CS%DX_dyBu(I,J)*((Del2u(I,j+1)*G%IdxCu(I,j+1)) - (Del2u(I,j)*G%IdxCu(I,j)))
1566 enddo
1567
1568 ! Adjust contributions to shearing strain on open boundaries.
15691875 if (apply_OBC) then ; if ((OBC%strain_config == OBC_STRAIN_ZERO) .or. &
1570 (OBC%strain_config == OBC_STRAIN_FREESLIP)) then
1571 !$omp target update from(dDel2vdx, dDel2udy)
15720 do n=1,OBC%number_of_segments
15730 J = OBC%segment(n)%HI%JsdB ; I = OBC%segment(n)%HI%IsdB
15740 if (OBC%segment(n)%is_N_or_S .and. (J >= js-1) .and. (J <= Jeq)) then
15750 do I=OBC%segment(n)%HI%IsdB,OBC%segment(n)%HI%IedB
15760 if (OBC%strain_config == OBC_STRAIN_ZERO) then
15770 dDel2vdx(I,J) = 0. ; dDel2udy(I,J) = 0.
15780 elseif (OBC%strain_config == OBC_STRAIN_FREESLIP) then
15790 dDel2udy(I,J) = 0.
1580 endif
1581 enddo
15820 elseif (OBC%segment(n)%is_E_or_W .and. (I >= is-1) .and. (I <= Ieq)) then
15830 do J=OBC%segment(n)%HI%JsdB,OBC%segment(n)%HI%JedB
15840 if (OBC%strain_config == OBC_STRAIN_ZERO) then
15850 dDel2vdx(I,J) = 0. ; dDel2udy(I,J) = 0.
15860 elseif (OBC%strain_config == OBC_STRAIN_FREESLIP) then
15870 dDel2vdx(I,J) = 0.
1588 endif
1589 enddo
1590 endif
1591 enddo
1592 !$omp target update to(dDel2vdx, dDel2udy)
1593 endif ; endif
1594 endif
1595
15961875 if ((CS%Smagorinsky_Kh) .or. (CS%Smagorinsky_Ah)) then
15971875 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
159813839375 sh_xy_sq = sh_xy(I,J)**2
1599 sh_xx_sq = 0.25 * ( ((sh_xx(i,j)**2) + (sh_xx(i+1,j+1)**2)) &
160013839375 + ((sh_xx(i,j+1)**2) + (sh_xx(i+1,j)**2)) )
160114068125 Shear_mag(I,J) = sqrt(sh_xy_sq + sh_xx_sq)
1602 enddo
1603 endif
1604
16051875 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
160613839375 h2uq = 4.0 * (h_u(I,j) * h_u(I,j+1))
160713839375 h2vq = 4.0 * (h_v(i,J) * h_v(i+1,J))
1608 hq(I,J) = (2.0 * (h2uq * h2vq)) &
160914068125 / (h_neglect3 + (h2uq + h2vq) * ((h_u(I,j) + h_u(I,j+1)) + (h_v(i,J) + h_v(i+1,J))))
1610 enddo
1611
16121875 if (CS%bound_Ah .or. CS%bound_Kh) then
16131875 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
161413839375 h_min = min(h_u(I,j), h_u(I,j+1), h_v(i,J), h_v(i+1,J))
161514068125 hrat_min(I,J) = min(1.0, h_min / (hq(I,J) + h_neglect))
1616 enddo
1617 endif
1618
1619 ! TODO: GPU?? Are h_[uv] on CPU? update to hrat_min?
16201875 if (CS%no_slip) then
16210 do J=js-1,Jeq ; do I=is-1,Ieq
16220 if (CS%no_slip .and. (G%mask2dBu(I,J) < 0.5)) then
16230 if ((G%mask2dCu(I,j) + G%mask2dCu(I,j+1)) + &
1624 (G%mask2dCv(i,J) + G%mask2dCv(i+1,J)) > 0.0) then
1625 ! This is a coastal vorticity point, so modify hq and hrat_min.
1626
16270 hu = G%mask2dCu(I,j) * h_u(I,j) + G%mask2dCu(I,j+1) * h_u(I,j+1)
16280 hv = G%mask2dCv(i,J) * h_v(i,J) + G%mask2dCv(i+1,J) * h_v(i+1,J)
16290 if ((G%mask2dCu(I,j) + G%mask2dCu(I,j+1)) * &
1630 (G%mask2dCv(i,J) + G%mask2dCv(i+1,J)) == 0.0) then
1631 ! Only one of hu and hv is nonzero, so just add them.
16320 hq(I,J) = hu + hv
16330 hrat_min(I,J) = 1.0
1634 else
1635 ! Both hu and hv are nonzero, so take the harmonic mean.
16360 hq(I,J) = 2.0 * (hu * hv) / ((hu + hv) + h_neglect)
16370 hrat_min(I,J) = min(1.0, min(hu, hv) / (hq(I,J) + h_neglect) )
1638 endif
1639 endif
1640 endif
1641 enddo ; enddo
1642 endif
1643
1644 ! Pass the velocity gradients and thickness to ZB2020
16451875 if (CS%use_ZB2020) then
1646 !$omp target update to(sh_xx, sh_xy, vort_xy, hq)
1647 call ZB2020_copy_gradient_and_thickness( &
1648 sh_xx, sh_xy, vort_xy, &
1649 hq, &
16500 G, GV, CS%ZB2020, k)
1651 endif
1652
1653 !!$omp target update from(sh_xx, sh_xy)
1654 !!$omp target update from(h_u, h_v, hq)
1655 !!$omp target update from(str_xx)
1656 !!$omp target update from(Shear_mag) if (use_Smag)
1657 !!$omp target update from(Del2u, Del2v) if (CS%biharmonic)
1658 !!$omp target update from(dDel2vdx, dDel2udy) if (CS%biharmonic)
1659
16601875 if (CS%Laplacian) then
1661 ! Determine the Laplacian viscosity at q points, using the
1662 ! largest value from several parameterizations. Also get the
1663 ! Laplacian component of str_xy.
1664
16651875 if ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then
16660 if (CS%use_QG_Leith_visc) then
16670 do J=js-1,Jeq ; do I=is-1,Ieq
16680 grad_vort = grad_vort_mag_q(I,J) + grad_div_mag_q(I,J)
16690 grad_vort_qg = 3. * grad_vort_mag_q_2d(I,J)
16700 vert_vort_mag(I,J) = min(grad_vort, grad_vort_qg)
1671 enddo ; enddo
1672 else
16730 do J=js-1,Jeq ; do I=is-1,Ieq
16740 vert_vort_mag(I,J) = grad_vort_mag_q(I,J) + grad_div_mag_q(I,J)
1675 enddo ; enddo
1676 endif
1677 endif
1678
1679 ! Static (pre-computed) background viscosity
16801875 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
168114068125 Kh(I,J) = CS%Kh_bg_xy(I,J)
1682 enddo
1683
16841875 if (CS%Smagorinsky_Kh) then
16850 if (CS%add_LES_viscosity) then
16860 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
16870 Kh(I,J) = Kh(I,J) + CS%Laplac2_const_xy(I,J) * Shear_mag(I,J)
1688 enddo
1689 else
16900 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
16910 Kh(I,J) = max(Kh(I,J), CS%Laplac2_const_xy(I,J) * Shear_mag(I,J) )
1692 enddo
1693 endif
1694 endif
1695
16961875 if (CS%Leith_Kh) then
1697 !$omp target update from(Kh)
16980 if (CS%add_LES_viscosity) then
16990 do J=js-1,Jeq ; do I=is-1,Ieq
17000 Kh(I,J) = Kh(I,J) + CS%Laplac3_const_xy(I,J) * vert_vort_mag(I,J) * inv_PI3 ! Is this right? -AJA
1701 enddo ; enddo
1702 else
17030 do J=js-1,Jeq ; do I=is-1,Ieq
17040 Kh(I,J) = max(Kh(I,J), CS%Laplac3_const_xy(I,J) * vert_vort_mag(I,J) * inv_PI3)
1705 enddo ; enddo
1706 endif
1707 !$omp target update to(Kh)
1708 endif
1709
1710 ! All viscosity contributions above are subject to resolution scaling
1711
17121875 if (rescale_Kh) then
1713 !$omp target update from(Kh)
171413955625 do J=js-1,Jeq ; do I=is-1,Ieq
171513953750 Kh(I,J) = VarMix%Res_fn_q(I,J) * Kh(I,J)
1716 enddo ; enddo
1717 !$omp target update to(Kh)
1718 endif
1719
17201875 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
172114068125 Kh(I,J) = max(Kh(I,J), CS%Kh_bg_min) ! Place a floor on the viscosity, if desired.
1722 enddo
1723
17241875 if (use_MEKE_Ku .and. .not. CS%EY24_EBT_BS) then
1725 !$omp target update from(Kh)
17260 if (use_kh_struct) then
17270 do J=js-1,Jeq ; do I=is-1,Ieq
17280 meke_res_fn = 1.
17290 if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_q(I,J)
1730
1731 Kh(I,J) = Kh(I,J) + 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + &
1732 (MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + &
1733 ((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + &
17340 (MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) ) * meke_res_fn
1735 enddo ; enddo
1736 else
17370 do J=js-1,Jeq ; do I=is-1,Ieq
17380 meke_res_fn = 1.
17390 if (CS%res_scale_MEKE) meke_res_fn = VarMix%Res_fn_q(I,J)
1740
1741 Kh(I,J) = Kh(I,J) + 0.25 * ( &
1742 (MEKE%Ku(i,j) + MEKE%Ku(i+1,j+1)) + &
1743 (MEKE%Ku(i+1,j) + &
17440 MEKE%Ku(i,j+1)) ) * meke_res_fn
1745 enddo ; enddo
1746 endif
1747 !$omp target update to(Kh)
1748 endif
1749
17501875 if (CS%anisotropic) then
1751 !$omp target update from(Kh)
1752 ! *Add* the shear component of anisotropic viscosity
17530 do J=js-1,Jeq ; do I=is-1,Ieq
17540 Kh(I,J) = Kh(I,J) + CS%Kh_aniso * CS%n1n2_q(I,J)**2
1755 enddo ; enddo
1756 !$omp target update to(Kh)
1757 endif
1758
17591875 if ((CS%bound_Kh) .and. (CS%bound_Ah)) then
1760 ! Newer method of bounding for stability
17611875 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
176213839375 visc_bound_rem(I,J) = 1.0
176313839375 Kh_max_here = hrat_min(I,J) * CS%Kh_Max_xy(I,J)
176414068125 if (Kh(I,J) >= Kh_max_here) then
17654185000 visc_bound_rem(I,J) = 0.0
17664185000 Kh(I,J) = Kh_max_here
17679654375 elseif ((Kh(I,J) > 0.0) .or. (CS%backscatter_underbound .and. (Kh_max_here > 0.0))) then
17680 visc_bound_rem(I,J) = 1.0 - Kh(I,J) / Kh_max_here
1769 endif
1770 enddo
17710 elseif (CS%bound_Kh) then
17720 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
17730 Kh(I,J) = min(Kh(I,J), hrat_min(I,J) * CS%Kh_Max_xy(I,J))
1774 enddo
1775 endif
1776
17771875 if (CS%use_Leithy) then
1778 ! Leith+E doesn't recompute Kh at q points, it just interpolates it from h to q points
17790 do J=js-1,Jeq ; do I=is-1,Ieq
17800 Kh(I,J) = 0.25 * ((Kh_h(i,j,k) + Kh_h(i+1,j+1,k)) + (Kh_h(i,j+1,k) + Kh_h(i+1,j,k)))
1781 enddo ; enddo
1782 !$omp target update to(Kh)
1783 endif
1784
17851875 if (CS%id_Kh_q > 0 .or. CS%debug) then
1786 !$omp target update from (Kh)
17870 do J=js-1,Jeq; do I=is-1,Ieq
17880 Kh_q(I,J,k) = Kh(I,J)
1789 enddo ; enddo
1790 endif
1791
17921875 if (CS%id_vort_xy_q > 0) then
17930 do J=js-1,Jeq ; do I=is-1,Ieq
17940 vort_xy_q(I,J,k) = vort_xy(I,J)
1795 enddo ; enddo
1796 endif
1797
17981875 if (CS%id_sh_xy_q > 0) then
17990 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
18000 sh_xy_q(I,J,k) = sh_xy(I,J)
1801 enddo
1802 endif
1803
18041875 if (.not. CS%use_Leithy) then
18051875 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
180614068125 str_xy(I,J) = -Kh(I,J) * sh_xy(I,J)
1807 enddo
1808 else
1809 !$omp target update from(Kh)
18100 do J=js-1,Jeq ; do I=is-1,Ieq
18110 str_xy(I,J) = -Kh(I,J) * sh_xy_smooth(I,J)
1812 enddo ; enddo
1813 !$omp target update to(str_xy)
1814 endif
1815 else
18160 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
18170 str_xy(I,J) = 0.
1818 enddo
1819 endif ! get harmonic coefficient Kh at q points and harmonic part of str_xy
1820
18211875 if (CS%anisotropic) then
1822 !$omp target update from(sh_xx, str_xy)
18230 do J=js-1,Jeq ; do I=is-1,Ieq
1824 ! Horizontal-tension averaged to q-points
18250 local_strain = 0.25 * ( (sh_xx(i,j) + sh_xx(i+1,j+1)) + (sh_xx(i+1,j) + sh_xx(i,j+1)) )
1826 ! *Add* the tension contribution to the xy-component of stress
18270 str_xy(I,J) = str_xy(I,J) - CS%Kh_aniso * CS%n1n2_q(I,J) * CS%n1n1_m_n2n2_q(I,J) * local_strain
1828 enddo ; enddo
1829 !$omp target update to(str_xy)
1830 endif
1831
18321875 if (CS%biharmonic) then
1833 ! Determine the biharmonic viscosity at q points, using the
1834 ! largest value from several parameterizations. Also get the
1835 ! biharmonic component of str_xy.
18361875 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
183714068125 Ah(I,J) = CS%Ah_bg_xy(I,J)
1838 enddo
1839
18401875 if (CS%Smagorinsky_Ah .or. CS%Leith_Ah) then
18411875 if (CS%Smagorinsky_Ah) then
18421875 if (CS%bound_Coriolis) then
18431875 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
1844 AhSm = Shear_mag(I,J) * (CS%Biharm_const_xy(I,J) &
184513839375 + CS%Biharm_const2_xy(I,J) * Shear_mag(I,J))
184614068125 Ah(I,J) = max(Ah(I,J), AhSm)
1847 enddo
1848 else
18490 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
18500 AhSm = CS%Biharm_const_xy(I,J) * Shear_mag(I,J)
18510 Ah(I,J) = max(Ah(I,J), AhSm)
1852 enddo
1853 endif
1854 endif
1855
18561875 if (CS%Leith_Ah) then
1857 !$omp target update from(Ah)
18580 do J=js-1,Jeq ; do I=is-1,Ieq
18590 AhLth = CS%Biharm6_const_xy(I,J) * abs(Del2vort_q(I,J)) * inv_PI6
18600 Ah(I,J) = max(Ah(I,J), AhLth)
1861 enddo ; enddo
1862 !$omp target update to(Ah)
1863 endif
1864 endif ! Smagorinsky_Ah or Leith_Ah
1865
18661875 if (use_MEKE_Au) then
1867 !$omp target update from(Ah)
1868 ! *Add* the MEKE contribution
18690 do J=js-1,Jeq ; do I=is-1,Ieq
1870 Ah(I,J) = Ah(I,J) + 0.25 * ( &
18710 (MEKE%Au(i,j) + MEKE%Au(i+1,j+1)) + (MEKE%Au(i+1,j) + MEKE%Au(i,j+1)) )
1872 enddo ; enddo
1873 !$omp target update to(Ah)
1874 endif
1875
1876 ! XXX: It is just overwrites the values!
18771875 if (CS%Re_Ah > 0.0) then
18780 do J=js-1,Jeq ; do I=is-1,Ieq
18790 KE = 0.125 * (((u(I,j,k) + u(I,j+1,k))**2) + ((v(i,J,k) + v(i+1,J,k))**2))
18800 Ah(I,J) = sqrt(KE) * CS%Re_Ah_const_xy(I,J)
1881 enddo ; enddo
1882 !$omp target update to(Ah)
1883 endif
1884
18851875 if (CS%bound_Ah) then
18861875 if (CS%bound_Kh) then
18871875 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
188814068125 Ah(I,J) = min(Ah(I,J), visc_bound_rem(I,J) * hrat_min(I,J) * CS%Ah_Max_xy(I,J))
1889 enddo
1890 else
18910 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
18920 Ah(I,J) = min(Ah(I,J), hrat_min(I,J) * CS%Ah_Max_xy(I,J))
1893 enddo
1894 endif
1895 endif
1896
18971875 if (CS%EY24_EBT_BS) then
1898 ! TODO: Fix indent!
1899 !$omp target update from(Ah, hrat_min)
19000 do J=js-1,Jeq ; do I=is-1,Ieq
19010 tmp = CS%KS_coef *hrat_min(I,J) * CS%Ah_Max_xy_KS(I,J)
19020 visc_limit_q(I,J,k) = tmp
19030 visc_limit_q_frac(i,j,k) = Ah(i,j) / (CS%KS_coef * hrat_min(i,j) * CS%Ah_Max_xy_KS(i,j))
19040 if (Ah(I,J) >= tmp) then
19050 visc_limit_q_flag(I,J,k) = 1.
1906 endif
1907 enddo ; enddo
1908 endif
1909
1910 ! Leith+E doesn't recompute Ah at q points, it just interpolates it from h to q points
19111875 if (CS%use_Leithy) then
19120 do J=js-1,Jeq ; do I=is-1,Ieq
19130 Ah(I,J) = 0.25 * ((Ah_h(i,j,k) + Ah_h(i+1,j+1,k)) + (Ah_h(i,j+1,k) + Ah_h(i+1,j,k)))
1914 enddo ; enddo
1915 !$omp target update to(Ah)
1916 endif
1917
19181875 if (CS%id_Ah_q>0 .or. CS%debug) then
1919 !$omp target update from(Ah)
19200 do J=js-1,Jeq ; do I=is-1,Ieq
19210 Ah_q(I,J,k) = Ah(I,J)
1922 enddo ; enddo
1923 endif
1924
1925 ! Again, need to initialize str_xy as if its biharmonic
19261875 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
192713839375 d_str = Ah(I,J) * (dDel2vdx(I,J) + dDel2udy(I,J))
1928
192913839375 str_xy(I,J) = str_xy(I,J) + d_str
1930
1931 ! Keep a copy of the biharmonic contribution for backscatter parameterization
1932 ! NOTE: computing this ought to be conditional! But it uses d_str...
193314068125 bhstr_xy(I,J) = d_str * (hq(I,J) * G%mask2dBu(I,J) * CS%reduction_xy(I,J))
1934 enddo
1935 endif ! Get Ah at q points and biharmonic part of str_xy
1936
1937 ! Backscatter using MEKE
19381875 if (CS%EY24_EBT_BS) then
1939 !$omp target update from(sh_xy, str_xy)
19400 do J=js-1,Jeq ; do I=is-1,Ieq
19410 if (visc_limit_q_flag(I,J,k) > 0) then
19420 Kh_BS(I,J) = 0.
1943 else
19440 if (use_kh_struct) then
1945 Kh_BS(I,J) = 0.25*( ((MEKE%Ku(i,j)*VarMix%BS_struct(i,j,k)) + &
1946 (MEKE%Ku(i+1,j+1)*VarMix%BS_struct(i+1,j+1,k))) + &
1947 ((MEKE%Ku(i+1,j)*VarMix%BS_struct(i+1,j,k)) + &
19480 (MEKE%Ku(i,j+1)*VarMix%BS_struct(i,j+1,k))) )
1949 else
1950 Kh_BS(I,J) = 0.25*( (MEKE%Ku(i,j) + MEKE%Ku(i+1,j+1)) + &
19510 (MEKE%Ku(i+1,j) + MEKE%Ku(i,j+1)) )
1952 endif
1953 endif
1954 enddo ; enddo
1955
19560 do J=js-1,Jeq ; do I=is-1,Ieq
19570 str_xy_BS(I,J) = -Kh_BS(I,J) * (sh_xy(I,J))
1958 enddo ; enddo
1959
19600 if (CS%id_BS_coeff_q>0) then
19610 do J=js-1,Jeq ; do I=is-1,Ieq
19620 BS_coeff_q(I,J,k) = Kh_BS(I,J)
1963 enddo ; enddo
1964 endif
1965
19660 do J=js-1,Jeq ; do I=is-1,Ieq
19670 str_xy(I,J) = str_xy(I,J) + str_xy_BS(I,J)
1968 enddo ; enddo
1969 !$omp target update to(str_xy)
1970 endif ! Backscatter
1971
19721875 if (CS%use_GME) then
1973 !$omp target update from(str_xx, str_xy)
1974 !$omp target update from(hq) if (CS%no_slip)
1975
1976 ! The wider halo here is to permit one pass of smoothing without a halo update.
19770 do j=Jsq-1,Jeq+2 ; do i=Isq-1,Ieq+2
1978 GME_coeff = GME_effic_h(i,j) * 0.25 * &
19790 ((KH_u_GME(I,j,k)+KH_u_GME(I-1,j,k)) + (KH_v_GME(i,J,k)+KH_v_GME(i,J-1,k)))
19800 GME_coeff = MIN(GME_coeff, CS%GME_limiter)
1981
19820 if ((CS%id_GME_coeff_h>0) .or. find_FrictWork) GME_coeff_h(i,j,k) = GME_coeff
19830 str_xx_GME(i,j) = GME_coeff * sh_xx_bt(i,j)
1984 enddo ; enddo
1985
1986 ! The wider halo here is to permit one pass of smoothing without a halo update.
19870 do J=js-2,je+1 ; do I=is-2,ie+1
1988 GME_coeff = GME_effic_q(I,J) * 0.25 * &
19890 ((KH_u_GME(I,j,k)+KH_u_GME(I,j+1,k)) + (KH_v_GME(i,J,k)+KH_v_GME(i+1,J,k)))
19900 GME_coeff = MIN(GME_coeff, CS%GME_limiter)
1991
19920 if (CS%id_GME_coeff_q>0) GME_coeff_q(I,J,k) = GME_coeff
19930 str_xy_GME(I,J) = GME_coeff * sh_xy_bt(I,J)
1994 enddo ; enddo
1995
1996 ! Applying GME diagonal term. This is linear and the arguments can be rescaled.
19970 call smooth_GME(CS, G, GME_flux_h=str_xx_GME)
19980 call smooth_GME(CS, G, GME_flux_q=str_xy_GME)
1999
2000 ! This changes the units of str_xx from [L2 T-2 ~> m2 s-2] to [H L2 T-2 ~> m3 s-2 or kg s-2].
20010 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
20020 str_xx(i,j) = (str_xx(i,j) + str_xx_GME(i,j)) * (h(i,j,k) * CS%reduction_xx(i,j))
2003 enddo ; enddo
2004
2005 ! This adds in GME and changes the units of str_xx from [L2 T-2 ~> m2 s-2] to [H L2 T-2 ~> m3 s-2 or kg s-2].
20060 if (CS%no_slip) then
20070 do J=js-1,Jeq ; do I=is-1,Ieq
20080 str_xy(I,J) = (str_xy(I,J) + str_xy_GME(I,J)) * (hq(I,J) * CS%reduction_xy(I,J))
2009 enddo ; enddo
2010 else
20110 do J=js-1,Jeq ; do I=is-1,Ieq
20120 str_xy(I,J) = (str_xy(I,J) + str_xy_GME(I,J)) * (hq(I,J) * G%mask2dBu(I,J) * CS%reduction_xy(I,J))
2013 enddo ; enddo
2014 endif
2015 !$omp target update to(str_xx, str_xy)
2016 else ! .not. use_GME
2017 ! This changes the units of str_xx from [L2 T-2 ~> m2 s-2] to [H L2 T-2 ~> m3 s-2 or kg s-2].
2018230625 do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
201914413125 str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * CS%reduction_xx(i,j))
2020 enddo
2021
2022 ! This changes the units of str_xy from [L2 T-2 ~> m2 s-2] to [H L2 T-2 ~> m3 s-2 or kg s-2].
20231875 if (CS%no_slip) then
20240 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
20250 str_xy(I,J) = str_xy(I,J) * (hq(I,J) * CS%reduction_xy(I,J))
2026 enddo
2027 else
20281875 do concurrent (J=js-1:Jeq, I=is-1:Ieq)
202914068125 str_xy(I,J) = str_xy(I,J) * (hq(I,J) * G%mask2dBu(I,J) * CS%reduction_xy(I,J))
2030 enddo
2031 endif
2032 endif ! use_GME
2033
2034 ! Evaluate 1/h x.Div(h Grad u) or the biharmonic equivalent.
20351875 do concurrent (j=js:je, I=Isq:Ieq)
2036 diffu(I,j,k) = ((G%IdxCu(I,j)*((CS%dx2q(I,J-1)*str_xy(I,J-1)) - (CS%dx2q(I,J)*str_xy(I,J))) + &
2037 G%IdyCu(I,j)*((CS%dy2h(i,j)*str_xx(i,j)) - (CS%dy2h(i+1,j)*str_xx(i+1,j)))) * &
203813841250 G%IareaCu(I,j)) / (h_u(I,j) + h_neglect)
2039 enddo
2040
20411875 if (apply_OBC) then
2042 !$omp target update from(diffu)
2043 ! This is not the right boundary condition. If all the masking of tendencies are done
2044 ! correctly later then eliminating this block should not change answers.
20450 do n=1,OBC%number_of_segments
20460 if (OBC%segment(n)%is_E_or_W) then
20470 I = OBC%segment(n)%HI%IsdB
20480 do j=OBC%segment(n)%HI%jsd,OBC%segment(n)%HI%jed
20490 diffu(I,j,k) = 0.
2050 enddo
2051 endif
2052 enddo
2053 !$omp target update to(diffu)
2054 endif
2055
2056 ! Evaluate 1/h y.Div(h Grad u) or the biharmonic equivalent.
20571875 do concurrent (J=Jsq:Jeq, i=is:ie)
2058 diffv(i,J,k) = ((G%IdyCv(i,J)*((CS%dy2q(I-1,J)*str_xy(I-1,J)) - (CS%dy2q(I,J)*str_xy(I,J))) - &
2059 G%IdxCv(i,J)*((CS%dx2h(i,j)*str_xx(i,j)) - (CS%dx2h(i,j+1)*str_xx(i,j+1)))) * &
206013951875 G%IareaCv(i,J)) / (h_v(i,J) + h_neglect)
2061 enddo
2062
20631875 if (apply_OBC) then
2064 !$omp target update from(diffv)
2065 ! This is not the right boundary condition. If all the masking of tendencies are done
2066 ! correctly later then eliminating this block should not change answers.
20670 do n=1,OBC%number_of_segments
20680 if (OBC%segment(n)%is_N_or_S) then
20690 J = OBC%segment(n)%HI%JsdB
20700 do i=OBC%segment(n)%HI%isd,OBC%segment(n)%HI%ied
20710 diffv(i,J,k) = 0.
2072 enddo
2073 endif
2074 enddo
2075 !$omp target update to(diffv)
2076 endif
2077
2078 !$omp target update from(h_u, h_v) &
2079 !$omp if ((find_Frictwork .or. find_FrictWork_bh) .and. .not. CS%FrictWork_bug)
2080
20811875 if (find_FrictWork) then
2082 !$omp target update from(str_xx, str_xy)
2083
20840 if (CS%FrictWork_bug) then
2085 ! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v)
2086 ! This is the old formulation that includes energy diffusion
20870 do j=js,je ; do i=is,ie
2088 FrictWork(i,j,k) = GV%H_to_RZ * ( &
2089 ((str_xx(i,j) * (u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j)) &
2090 - (str_xx(i,j) * (v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j))) &
2091 + 0.25*(( (str_xy(I,J) * &
2092 (((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J)) &
2093 + ((v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J)))) &
2094 + (str_xy(I-1,J-1) * &
2095 (((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1)) &
2096 + ((v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)))) ) &
2097 + ( (str_xy(I-1,J) * &
2098 (((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J)) &
2099 + ((v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)))) &
2100 + (str_xy(I,J-1) * &
2101 (((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1)) &
21020 + ((v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)))) ) ) )
2103 enddo ; enddo
2104 else
21050 do j=js,je ; do i=is,ie
2106 FrictWork(i,j,k) = GV%H_to_RZ * G%IareaT(i,j) * ( &
2107 ((str_xx(i,j)*CS%dy2h(i,j) * ( &
2108 (uh(I,j,k)*G%dxCu(I,j)*G%IdyCu(I,j)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) &
2109 - (uh(I-1,j,k)*G%dxCu(I-1,j)*G%IdyCu(I-1,j)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) ) ) &
2110 - (str_xx(i,j)*CS%dx2h(i,j) * ( &
2111 (vh(i,J,k)*G%dyCv(i,J)*G%IdxCv(i,J)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) &
2112 - (vh(i,J-1,k)*G%dyCv(i,J-1)*G%IdxCv(i,J-1)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) ) )) &
2113 + (0.25*(((str_xy(I,J)*( &
2114 (CS%dx2q(I,J)*((uh(I,j+1,k)*G%IareaCu(I,j+1)/(h_u(I,j+1)+h_neglect)) &
2115 - (uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)))) &
2116 + (CS%dy2q(I,J)*((vh(i+1,J,k)*G%IareaCv(i+1,J)/(h_v(i+1,J)+h_neglect)) &
2117 - (vh(i,J,k)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)))) )) &
2118 +(str_xy(I-1,J-1)*( &
2119 (CS%dx2q(I-1,J-1)*((uh(I-1,j,k)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) &
2120 - (uh(I-1,j-1,k)*G%IareaCu(I-1,j-1)/(h_u(I-1,j-1)+h_neglect)))) &
2121 + (CS%dy2q(I-1,J-1)*((vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) &
2122 - (vh(i-1,J-1,k)*G%IareaCv(i-1,J-1)/(h_v(i-1,J-1)+h_neglect)))) )) ) &
2123 +((str_xy(I-1,J)*( &
2124 (CS%dx2q(I-1,J)*((uh(I-1,j+1,k)*G%IareaCu(I-1,j+1)/(h_u(I-1,j+1)+h_neglect)) &
2125 - (uh(I-1,j,k)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)))) &
2126 + (CS%dy2q(I-1,J)*((vh(i,J,k)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) &
2127 - (vh(i-1,J,k)*G%IareaCv(i-1,J)/(h_v(i-1,J)+h_neglect)))) )) &
2128 +(str_xy(I,J-1)*( &
2129 (CS%dx2q(I,J-1)*((uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) &
2130 - (uh(I,j-1,k)*G%IareaCu(I,j-1)/(h_u(I,j-1)+h_neglect)))) &
2131 + (CS%dy2q(I,J-1)*((vh(i+1,J-1,k)*G%IareaCv(i+1,J-1)/(h_v(i+1,J-1)+h_neglect)) &
21320 - (vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)))) )) ) )) )
2133
2134 enddo ; enddo
2135 endif
2136
21370 if (CS%EY24_EBT_BS) then
21380 do j=js,je ; do i=is,ie
21390 FrictWork(i,j,k) = (1. - visc_limit_h_flag(i,j,k)) * FrictWork(i,j,k)
2140 enddo ; enddo
2141 endif
2142 endif
2143
21441875 if (find_FrictWork_bh) then
2145 !$omp target update from(bhstr_xx, bhstr_xy)
21460 if (CS%FrictWork_bug) then
2147 ! Diagnose bhstr_xx*d_x u - bhstr_yy*d_y v + bhstr_xy*(d_y u + d_x v)
2148 ! This is the old formulation that includes energy diffusion !cyc
21490 do j=js,je ; do i=is,ie
2150 FrictWork_bh(i,j,k) = GV%H_to_RZ * ( &
2151 ((bhstr_xx(i,j) * (u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j)) &
2152 - (bhstr_xx(i,j) * (v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j))) &
2153 + 0.25*(( (bhstr_xy(I,J) * &
2154 (((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J)) &
2155 + ((v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J)))) &
2156 + (bhstr_xy(I-1,J-1) * &
2157 (((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1)) &
2158 + ((v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)))) ) &
2159 + ( (bhstr_xy(I-1,J) * &
2160 (((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J)) &
2161 + ((v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)))) &
2162 + (bhstr_xy(I,J-1) * &
2163 (((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1)) &
21640 + ((v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)))) ) ) )
2165 enddo ; enddo
2166 else
21670 do j=js,je ; do i=is,ie
2168 ! Diagnose bhstr_xx*d_x u - bhstr_yy*d_y v + bhstr_xy*(d_y u + d_x v)
2169 FrictWork_bh(i,j,k) = GV%H_to_RZ * G%IareaT(i,j) * ( &
2170 ((bhstr_xx(i,j)*CS%dy2h(i,j) * ( &
2171 (uh(I,j,k)*G%dxCu(I,j)*G%IdyCu(I,j)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) &
2172 - (uh(I-1,j,k)*G%dxCu(I-1,j)*G%IdyCu(I-1,j)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) ) ) &
2173 - (bhstr_xx(i,j)*CS%dx2h(i,j) * ( &
2174 (vh(i,J,k)*G%dyCv(i,J)*G%IdxCv(i,J)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) &
2175 - (vh(i,J-1,k)*G%dyCv(i,J-1)*G%IdxCv(i,J-1)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) ) )) &
2176 + (0.25*(((bhstr_xy(I,J)*( &
2177 (CS%dx2q(I,J)*((uh(I,j+1,k)*G%IareaCu(I,j+1)/(h_u(I,j+1)+h_neglect)) &
2178 - (uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)))) &
2179 + (CS%dy2q(I,J)*((vh(i+1,J,k)*G%IareaCv(i+1,J)/(h_v(i+1,J)+h_neglect)) &
2180 - (vh(i,J,k)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)))) )) &
2181 +(bhstr_xy(I-1,J-1)*( &
2182 (CS%dx2q(I-1,J-1)*((uh(I-1,j,k)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) &
2183 - (uh(I-1,j-1,k)*G%IareaCu(I-1,j-1)/(h_u(I-1,j-1)+h_neglect)))) &
2184 + (CS%dy2q(I-1,J-1)*((vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) &
2185 - (vh(i-1,J-1,k)*G%IareaCv(i-1,J-1)/(h_v(i-1,J-1)+h_neglect)))) )) ) &
2186 +((bhstr_xy(I-1,J)*( &
2187 (CS%dx2q(I-1,J)*((uh(I-1,j+1,k)*G%IareaCu(I-1,j+1)/(h_u(I-1,j+1)+h_neglect)) &
2188 - (uh(I-1,j,k)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)))) &
2189 + (CS%dy2q(I-1,J)*((vh(i,J,k)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) &
2190 - (vh(i-1,J,k)*G%IareaCv(i-1,J)/(h_v(i-1,J)+h_neglect)))) )) &
2191 +(bhstr_xy(I,J-1)*( &
2192 (CS%dx2q(I,J-1)*((uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) &
2193 - (uh(I,j-1,k)*G%IareaCu(I,j-1)/(h_u(I,j-1)+h_neglect)))) &
2194 + (CS%dy2q(I,J-1)*((vh(i+1,J-1,k)*G%IareaCv(i+1,J-1)/(h_v(i+1,J-1)+h_neglect)) &
21950 - (vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)))) )) ) )) )
2196 enddo ; enddo
2197 endif
2198
21990 if (CS%EY24_EBT_BS) then
22000 do j=js,je ; do i=is,ie
22010 FrictWork_bh(i,j,k) = (1. - visc_limit_h_flag(i,j,k)) * FrictWork_bh(i,j,k)
2202 enddo ; enddo
2203 endif
2204 endif
2205
22061875 if (CS%use_GME) then
22070 if (CS%FrictWork_bug) then ; do j=js,je ; do i=is,ie
2208 ! Diagnose str_xx_GME*d_x u - str_yy_GME*d_y v + str_xy_GME*(d_y u + d_x v)
2209 ! This is the old formulation that includes energy diffusion
2210 FrictWork_GME(i,j,k) = GV%H_to_RZ * ( &
2211 ((str_xx_GME(i,j)*(u(I,j,k)-u(I-1,j,k))*G%IdxT(i,j)) &
2212 - (str_xx_GME(i,j)*(v(i,J,k)-v(i,J-1,k))*G%IdyT(i,j))) &
2213 + 0.25*(( (str_xy_GME(I,J) * &
2214 (((u(I,j+1,k)-u(I,j,k))*G%IdyBu(I,J)) &
2215 + ((v(i+1,J,k)-v(i,J,k))*G%IdxBu(I,J)))) &
2216 + (str_xy_GME(I-1,J-1) * &
2217 (((u(I-1,j,k)-u(I-1,j-1,k))*G%IdyBu(I-1,J-1)) &
2218 + ((v(i,J-1,k)-v(i-1,J-1,k))*G%IdxBu(I-1,J-1)))) ) &
2219 + ( (str_xy_GME(I-1,J) * &
2220 (((u(I-1,j+1,k)-u(I-1,j,k))*G%IdyBu(I-1,J)) &
2221 + ((v(i,J,k)-v(i-1,J,k))*G%IdxBu(I-1,J)))) &
2222 + (str_xy_GME(I,J-1) * &
2223 (((u(I,j,k)-u(I,j-1,k))*G%IdyBu(I,J-1)) &
22240 + ((v(i+1,J-1,k)-v(i,J-1,k))*G%IdxBu(I,J-1)))) ) ) )
2225 enddo ; enddo
22260 else ; do j=js,je ; do i=is,ie
2227 FrictWork_GME(i,j,k) = GV%H_to_RZ * G%IareaT(i,j) * ( &
2228 ((str_xx_GME(i,j)*CS%dy2h(i,j) * ( &
2229 (uh(I,j,k)*G%dxCu(I,j)*G%IdyCu(I,j)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) &
2230 - (uh(I-1,j,k)*G%dxCu(I-1,j)*G%IdyCu(I-1,j)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) ) ) &
2231 - (str_xx_GME(i,j)*CS%dx2h(i,j) * ( &
2232 (vh(i,J,k)*G%dyCv(i,J)*G%IdxCv(i,J)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) &
2233 - (vh(i,J-1,k)*G%dyCv(i,J-1)*G%IdxCv(i,J-1)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) ) )) &
2234 + (0.25*(((str_xy_GME(I,J)*( &
2235 (CS%dx2q(I,J)*((uh(I,j+1,k)*G%IareaCu(I,j+1)/(h_u(I,j+1)+h_neglect)) &
2236 - (uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)))) &
2237 + (CS%dy2q(I,J)*((vh(i+1,J,k)*G%IareaCv(i+1,J)/(h_v(i+1,J)+h_neglect)) &
2238 - (vh(i,J,k)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)))) )) &
2239 +(str_xy_GME(I-1,J-1)*( &
2240 (CS%dx2q(I-1,J-1)*((uh(I-1,j,k)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)) &
2241 - (uh(I-1,j-1,k)*G%IareaCu(I-1,j-1)/(h_u(I-1,j-1)+h_neglect)))) &
2242 + (CS%dy2q(I-1,J-1)*((vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)) &
2243 - (vh(i-1,J-1,k)*G%IareaCv(i-1,J-1)/(h_v(i-1,J-1)+h_neglect)))) )) ) &
2244 +((str_xy_GME(I-1,J)*( &
2245 (CS%dx2q(I-1,J)*((uh(I-1,j+1,k)*G%IareaCu(I-1,j+1)/(h_u(I-1,j+1)+h_neglect)) &
2246 - (uh(I-1,j,k)*G%IareaCu(I-1,j)/(h_u(I-1,j)+h_neglect)))) &
2247 + (CS%dy2q(I-1,J)*((vh(i,J,k)*G%IareaCv(i,J)/(h_v(i,J)+h_neglect)) &
2248 - (vh(i-1,J,k)*G%IareaCv(i-1,J)/(h_v(i-1,J)+h_neglect)))) )) &
2249 +(str_xy_GME(I,J-1)*( &
2250 (CS%dx2q(I,J-1)*((uh(I,j,k)*G%IareaCu(I,j)/(h_u(I,j)+h_neglect)) &
2251 - (uh(I,j-1,k)*G%IareaCu(I,j-1)/(h_u(I,j-1)+h_neglect)))) &
2252 + (CS%dy2q(I,J-1)*((vh(i+1,J-1,k)*G%IareaCv(i+1,J-1)/(h_v(i+1,J-1)+h_neglect)) &
22530 - (vh(i,J-1,k)*G%IareaCv(i,J-1)/(h_v(i,J-1)+h_neglect)))) )) ) )) )
2254 enddo ; enddo ; endif
2255 endif
2256
22571875 if (skeb_use_frict) then ; do j=js,je ; do i=is,ie
2258 ! Note that the sign convention is FrictWork < 0 means energy dissipation.
2259 STOCH%skeb_diss(i,j,k) = STOCH%skeb_diss(i,j,k) - STOCH%skeb_frict_coef * &
22600 FrictWork(i,j,k) / (GV%H_to_RZ * (h(i,j,k) + h_neglect))
2261 enddo ; enddo ; endif
2262
2263 ! Make a similar calculation as for FrictWork above but accumulating into
2264 ! the vertically integrated MEKE source term, and adjusting for any
2265 ! energy loss seen as a reduction in the (biharmonic) frictional source term.
22661900 if (find_FrictWork .and. allocated(MEKE%mom_src)) then
22670 if (k==1) then
22680 do j=js,je ; do i=is,ie
22690 MEKE%mom_src(i,j) = 0.
2270 enddo ; enddo
2271
22720 if (allocated(MEKE%mom_src_bh)) then
22730 do j=js,je ; do i=is,ie
22740 MEKE%mom_src_bh(i,j) = 0.
2275 enddo ; enddo
2276 endif
2277
22780 if (allocated(MEKE%GME_snk)) then
22790 do j=js,je ; do i=is,ie
22800 MEKE%GME_snk(i,j) = 0.
2281 enddo ; enddo
2282 endif
2283 endif
22840 if (MEKE%backscatter_Ro_c /= 0.) then
2285 !$omp target update from(sh_xx, sh_xy)
22860 do j=js,je ; do i=is,ie
2287 FatH = 0.25*( (abs(G%CoriolisBu(I-1,J-1)) + abs(G%CoriolisBu(I,J))) + &
22880 (abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J-1))) )
2289 Shear_mag_bc = sqrt(sh_xx(i,j) * sh_xx(i,j) + &
2290 0.25*(((sh_xy(I-1,J-1)*sh_xy(I-1,J-1)) + (sh_xy(I,J)*sh_xy(I,J))) + &
22910 ((sh_xy(I-1,J)*sh_xy(I-1,J)) + (sh_xy(I,J-1)*sh_xy(I,J-1)))))
22920 if ((CS%answer_date > 20190101) .and. (CS%answer_date < 20241201)) then
22930 FatH = (US%s_to_T*FatH)**MEKE%backscatter_Ro_pow ! f^n
2294 ! Note the hard-coded dimensional constant in the following line that can not
2295 ! be rescaled for dimensional consistency.
2296 Shear_mag_bc = (((US%s_to_T * Shear_mag_bc)**MEKE%backscatter_Ro_pow) + 1.e-30) &
22970 * MEKE%backscatter_Ro_c ! c * D^n
2298 ! The Rossby number function is g(Ro) = 1/(1+c.Ro^n)
2299 ! RoScl = 1 - g(Ro)
23000 RoScl = Shear_mag_bc / (FatH + Shear_mag_bc) ! = 1 - f^n/(f^n+c*D^n)
2301 else
23020 if (FatH <= backscat_subround*Shear_mag_bc) then
23030 RoScl = 1.0
2304 else
23050 Sh_F_pow = MEKE%backscatter_Ro_c * (Shear_mag_bc / FatH)**MEKE%backscatter_Ro_pow
23060 RoScl = Sh_F_pow / (1.0 + Sh_F_pow) ! = 1 - f^n/(f^n+c*D^n)
2307 endif
2308 endif
2309
23100 MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + (FrictWork(i,j,k) - RoScl*FrictWork_bh(i,j,k))
2311
23120 if (allocated(MEKE%mom_src_bh)) &
2313 MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) &
23140 + (FrictWork_bh(i,j,k) - RoScl * FrictWork_bh(i,j,k))
2315 enddo ; enddo
2316 else
23170 do j=js,je ; do i=is,ie
23180 MEKE%mom_src(i,j) = MEKE%mom_src(i,j) + FrictWork(i,j,k)
2319 enddo ; enddo
2320
23210 if (allocated(MEKE%mom_src_bh)) then
23220 do j=js,je ; do i=is,ie
23230 MEKE%mom_src_bh(i,j) = MEKE%mom_src_bh(i,j) + FrictWork_bh(i,j,k)
2324 enddo ; enddo
2325 endif
2326 endif ! MEKE%backscatter_Ro_c
2327
23280 if (CS%use_GME .and. allocated(MEKE%GME_snk)) then
23290 do j=js,je ; do i=is,ie
23300 MEKE%GME_snk(i,j) = MEKE%GME_snk(i,j) + FrictWork_GME(i,j,k)
2331 enddo ; enddo
2332 endif
2333 endif ! find_FrictWork and associated(mom_src)
2334 enddo ! end of k loop
2335
2336 !$omp target exit data map(delete: dudx, dudy, dvdx, dvdy, sh_xx, sh_xy)
2337 !$omp target exit data map(delete: h_u, h_v, hq)
2338 !$omp target exit data map(delete: str_xx, str_xy)
2339 !$omp target exit data map(delete: Del2u, Del2v) if (CS%biharmonic)
2340 !$omp target exit data map(delete: dDel2vdx, dDel2udy) if (CS%biharmonic)
2341 !$omp target exit data map(delete: Shear_mag) if (use_Smag)
2342 !$omp target exit data map(delete: Kh) if (CS%Laplacian)
2343 !$omp target exit data map(delete: Ah) if (CS%biharmonic)
2344 !$omp target exit data map(delete: bhstr_xx, bhstr_xy) if (CS%biharmonic)
2345
2346 !$omp target exit data map(delete: hrat_min) &
2347 !$omp if (CS%bound_Kh .or. CS%bound_Ah)
2348 !$omp target exit data map(delete: visc_bound_rem) &
2349 !$omp if (CS%bound_Kh .or. CS%bound_Ah)
2350
2351 ! Offer fields for diagnostic averaging.
235225 if (CS%id_normstress > 0) call post_data(CS%id_normstress, NoSt, CS%diag)
235325 if (CS%id_shearstress > 0) call post_data(CS%id_shearstress, ShSt, CS%diag)
235425 if (CS%id_diffu>0) call post_data(CS%id_diffu, diffu, CS%diag)
235525 if (CS%id_diffv>0) call post_data(CS%id_diffv, diffv, CS%diag)
235625 if (CS%id_FrictWork>0) call post_data(CS%id_FrictWork, FrictWork, CS%diag)
235725 if (CS%id_FrictWork_bh>0) call post_data(CS%id_FrictWork_bh, FrictWork_bh, CS%diag)
235825 if (CS%id_Ah_h>0) call post_data(CS%id_Ah_h, Ah_h, CS%diag)
235925 if (CS%id_grid_Re_Ah>0) call post_data(CS%id_grid_Re_Ah, grid_Re_Ah, CS%diag)
236025 if (CS%id_div_xx_h>0) call post_data(CS%id_div_xx_h, div_xx_h, CS%diag)
236125 if (CS%id_vort_xy_q>0) call post_data(CS%id_vort_xy_q, vort_xy_q, CS%diag)
236225 if (CS%id_sh_xx_h>0) call post_data(CS%id_sh_xx_h, sh_xx_h, CS%diag)
236325 if (CS%id_sh_xy_q>0) call post_data(CS%id_sh_xy_q, sh_xy_q, CS%diag)
236425 if (CS%id_Ah_q>0) call post_data(CS%id_Ah_q, Ah_q, CS%diag)
236525 if (CS%id_Kh_h>0) call post_data(CS%id_Kh_h, Kh_h, CS%diag)
236625 if (CS%id_grid_Re_Kh>0) call post_data(CS%id_grid_Re_Kh, grid_Re_Kh, CS%diag)
236725 if (CS%id_Kh_q>0) call post_data(CS%id_Kh_q, Kh_q, CS%diag)
236825 if (CS%use_GME) then ! post barotropic tension and strain
23690 if (CS%id_GME_coeff_h > 0) call post_data(CS%id_GME_coeff_h, GME_coeff_h, CS%diag)
23700 if (CS%id_GME_coeff_q > 0) call post_data(CS%id_GME_coeff_q, GME_coeff_q, CS%diag)
23710 if (CS%id_FrictWork_GME>0) call post_data(CS%id_FrictWork_GME, FrictWork_GME, CS%diag)
23720 if (CS%id_dudx_bt > 0) call post_data(CS%id_dudx_bt, dudx_bt, CS%diag)
23730 if (CS%id_dvdy_bt > 0) call post_data(CS%id_dvdy_bt, dvdy_bt, CS%diag)
23740 if (CS%id_dudy_bt > 0) call post_data(CS%id_dudy_bt, dudy_bt, CS%diag)
23750 if (CS%id_dvdx_bt > 0) call post_data(CS%id_dvdx_bt, dvdx_bt, CS%diag)
2376 endif
237725 if (CS%EY24_EBT_BS) then
23780 if (CS%id_visc_limit_h>0) call post_data(CS%id_visc_limit_h, visc_limit_h, CS%diag)
23790 if (CS%id_visc_limit_q>0) call post_data(CS%id_visc_limit_q, visc_limit_q, CS%diag)
23800 if (CS%id_visc_limit_h_frac>0) call post_data(CS%id_visc_limit_h_frac, visc_limit_h_frac, CS%diag)
23810 if (CS%id_visc_limit_q_frac>0) call post_data(CS%id_visc_limit_q_frac, visc_limit_q_frac, CS%diag)
23820 if (CS%id_visc_limit_h_flag>0) call post_data(CS%id_visc_limit_h_flag, visc_limit_h_flag, CS%diag)
23830 if (CS%id_visc_limit_q_flag>0) call post_data(CS%id_visc_limit_q_flag, visc_limit_q_flag, CS%diag)
23840 if (CS%id_BS_coeff_h>0) call post_data(CS%id_BS_coeff_h, BS_coeff_h, CS%diag)
23850 if (CS%id_BS_coeff_q>0) call post_data(CS%id_BS_coeff_q, BS_coeff_q, CS%diag)
2386 endif
2387
238825 if (CS%debug) then
23890 if (CS%Laplacian) then
23900 call hchksum(Kh_h, "Kh_h", G%HI, haloshift=1, unscale=US%L_to_m**2*US%s_to_T)
23910 call Bchksum(Kh_q, "Kh_q", G%HI, haloshift=0, symmetric=.true., unscale=US%L_to_m**2*US%s_to_T)
2392 endif
23930 if (CS%biharmonic) then
23940 call hchksum(Ah_h, "Ah_h", G%HI, haloshift=1, unscale=US%L_to_m**4*US%s_to_T)
23950 call Bchksum(Ah_q, "Ah_q", G%HI, haloshift=0, symmetric=.true., unscale=US%L_to_m**4*US%s_to_T)
2396 endif
2397 endif
2398
2399 !$omp target exit data map(delete: sh_xy_q) &
2400 !$omp if (CS%id_sh_xy_q > 0)
2401
240225 if (CS%id_FrictWorkIntz > 0) then
24030 do j=js,je
24040 do i=is,ie ; FrictWorkIntz(i,j) = FrictWork(i,j,1) ; enddo
24050 do k=2,nz ; do i=is,ie
24060 FrictWorkIntz(i,j) = FrictWorkIntz(i,j) + FrictWork(i,j,k)
2407 enddo ; enddo
2408 enddo
24090 call post_data(CS%id_FrictWorkIntz, FrictWorkIntz, CS%diag)
2410 endif
2411
241225 if (CS%id_FrictWorkIntz_bh > 0) then
24130 do j=js,je
24140 do i=is,ie ; FrictWorkIntz_bh(i,j) = FrictWork_bh(i,j,1) ; enddo
24150 do k=2,nz ; do i=is,ie
24160 FrictWorkIntz_bh(i,j) = FrictWorkIntz_bh(i,j) + FrictWork_bh(i,j,k)
2417 enddo ; enddo
2418 enddo
24190 call post_data(CS%id_FrictWorkIntz_bh, FrictWorkIntz_bh, CS%diag)
2420 endif
2421
242225 if (present(ADp)) then
2423 ! Diagnostics of the fractional thicknesses times momentum budget terms
2424 ! 3D diagnostics of hf_diffu(diffv) are commented because there is no clarity on proper remapping grid option.
2425 ! The code is retained for debugging purposes in the future.
2426 !if (CS%id_hf_diffu > 0) call post_product_u(CS%id_hf_diffu, diffu, ADp%diag_hfrac_u, G, nz, CS%diag)
2427 !if (CS%id_hf_diffv > 0) call post_product_v(CS%id_hf_diffv, diffv, ADp%diag_hfrac_v, G, nz, CS%diag)
2428
2429 ! Diagnostics for thickness-weighted vertically averaged momentum budget terms
243024 if (CS%id_hf_diffu_2d > 0) call post_product_sum_u(CS%id_hf_diffu_2d, diffu, ADp%diag_hfrac_u, G, nz, CS%diag)
243124 if (CS%id_hf_diffv_2d > 0) call post_product_sum_v(CS%id_hf_diffv_2d, diffv, ADp%diag_hfrac_v, G, nz, CS%diag)
2432
2433 ! Diagnostics for the vertical sum of layer thickness x momentum budget terms
243424 if (CS%id_intz_diffu_2d > 0) call post_product_sum_u(CS%id_intz_diffu_2d, diffu, ADp%diag_hu, G, nz, CS%diag)
243524 if (CS%id_intz_diffv_2d > 0) call post_product_sum_v(CS%id_intz_diffv_2d, diffv, ADp%diag_hv, G, nz, CS%diag)
2436
2437 ! Diagnostics for thickness x momentum budget terms
243824 if (CS%id_h_diffu > 0) call post_product_u(CS%id_h_diffu, diffu, ADp%diag_hu, G, nz, CS%diag)
243924 if (CS%id_h_diffv > 0) call post_product_v(CS%id_h_diffv, diffv, ADp%diag_hv, G, nz, CS%diag)
2440
2441 ! Diagnostics for momentum budget terms multiplied by visc_rem_[uv],
244224 if (CS%id_diffu_visc_rem > 0) call post_product_u(CS%id_diffu_visc_rem, diffu, ADp%visc_rem_u, G, nz, CS%diag)
244324 if (CS%id_diffv_visc_rem > 0) call post_product_v(CS%id_diffv_visc_rem, diffv, ADp%visc_rem_v, G, nz, CS%diag)
2444 endif
2445
244625 if (CS%use_ZB2020) then
2447 call ZB2020_lateral_stress(u, v, h, diffu, diffv, G, GV, CS%ZB2020, &
24480 CS%dx2h, CS%dy2h, CS%dx2q, CS%dy2q)
2449 endif
2450
245150end subroutine horizontal_viscosity
2452
2453!> Allocates space for and calculates static variables used by horizontal_viscosity.
2454!! hor_visc_init calculates and stores the values of a number of metric functions that
2455!! are used in horizontal_viscosity.
24561subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)
2457 type(time_type), intent(in) :: Time !< Current model time.
2458 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
2459 type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
2460 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2461 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
2462 !! parameters.
2463 type(diag_ctrl), target, intent(inout) :: diag !< Structure to regulate diagnostic output.
2464 type(hor_visc_CS), intent(inout) :: CS !< Horizontal viscosity control structure
2465 type(accel_diag_ptrs), intent(in), optional :: ADp !< Acceleration diagnostics
2466
2467 ! u0v is the Laplacian sensitivities to the v velocities at u points, with u0u, v0u, and v0v defined analogously.
24683 real, dimension(SZIB_(G),SZJ_(G)) :: u0u, u0v ! Laplacian sensitivities at u points [L-2 ~> m-2]
24691 real, dimension(SZI_(G),SZJB_(G)) :: v0u, v0v ! Laplacian sensitivities at v points [L-2 ~> m-2]
2470 real :: grid_sp_h2 ! Harmonic mean of the squares of the grid [L2 ~> m2]
2471 real :: grid_sp_h3 ! Harmonic mean of the squares of the grid^(3/2) [L3 ~> m3]
2472 real :: grid_sp_q2 ! spacings at h and q points [L2 ~> m2]
2473 real :: grid_sp_q3 ! spacings at h and q points^(3/2) [L3 ~> m3]
2474 real :: min_grid_sp_h2 ! Minimum value of grid_sp_h2 [L2 ~> m2]
2475 real :: min_grid_sp_h4 ! Minimum value of grid_sp_h2**2 [L4 ~> m4]
2476 real :: Kh_Limit ! A coefficient [T-1 ~> s-1] used, along with the
2477 ! grid spacing, to limit Laplacian viscosity.
2478 real :: fmax ! maximum absolute value of f at the four
2479 ! vorticity points around a thickness point [T-1 ~> s-1]
2480 real :: BoundCorConst ! A constant used when using viscosity to bound the Coriolis accelerations
2481 ! [T2 L-2 ~> s2 m-2]
2482 real :: Ah_Limit ! coefficient [T-1 ~> s-1] used, along with the
2483 ! grid spacing, to limit biharmonic viscosity
2484 real :: Kh ! Lapacian horizontal viscosity [L2 T-1 ~> m2 s-1]
2485 real :: Ah ! biharmonic horizontal viscosity [L4 T-1 ~> m4 s-1]
2486 real :: Kh_vel_scale ! this speed [L T-1 ~> m s-1] times grid spacing gives Laplacian viscosity
2487 real :: Ah_vel_scale ! this speed [L T-1 ~> m s-1] times grid spacing cubed gives biharmonic viscosity
2488 real :: Ah_time_scale ! damping time-scale for biharmonic visc [T ~> s]
2489 real :: Smag_Lap_const ! nondimensional Laplacian Smagorinsky constant [nondim]
2490 real :: Smag_bi_const ! nondimensional biharmonic Smagorinsky constant [nondim]
2491 real :: Leith_Lap_const ! nondimensional Laplacian Leith constant [nondim]
2492 real :: Leith_bi_const ! nondimensional biharmonic Leith constant [nondim]
2493 real :: dt ! The dynamics time step [T ~> s]
2494 real :: Idt ! The inverse of dt [T-1 ~> s-1]
2495 real :: denom ! work variable; the denominator of a fraction [L-2 ~> m-2] or [L-4 ~> m-4]
2496 real :: maxvel ! largest permitted velocity components [L T-1 ~> m s-1]
2497 real :: bound_Cor_vel ! grid-scale velocity variations at which value
2498 ! the quadratically varying biharmonic viscosity
2499 ! balances Coriolis acceleration [L T-1 ~> m s-1]
2500 real :: Kh_sin_lat ! Amplitude of latitudinally dependent viscosity [L2 T-1 ~> m2 s-1]
2501 real :: Kh_pwr_of_sine ! Power used to raise sin(lat) when using Kh_sin_lat [nondim]
2502 logical :: bound_Cor_def ! parameter setting of BOUND_CORIOLIS
2503 logical :: split ! If true, use the split time stepping scheme.
2504 ! If false and USE_GME = True, issue a FATAL error.
2505 logical :: use_MEKE ! If true, the MEKE parameterization is in use.
2506 logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to
2507 ! recreate the bugs, or if false bugs are only used if actively selected.
2508 real :: backscatter_Ro_c ! Coefficient in Rossby number function for backscatter [nondim]
2509 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags
2510 character(len=200) :: inputdir, filename ! Input file names and paths
2511 character(len=80) :: Kh_var ! Input variable names
2512 real :: deg2rad ! Converts degrees to radians [radians degree-1]
2513 real :: slat_fn ! sin(lat)**Kh_pwr_of_sine [nondim]
2514 real :: aniso_grid_dir(2) ! Vector (n1,n2) for anisotropic direction [nondim]
2515 integer :: aniso_mode ! Selects the mode for setting the anisotropic direction
2516 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
2517 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
2518 integer :: i, j
2519 ! This include declares and sets the variable "version".
2520# include "version_variable.h"
2521 character(len=40) :: mdl = "MOM_hor_visc" ! module name
25221 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
25231 Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
25241 isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
25251 IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
2526
2527 ! init control structure
25281 call ZB2020_init(Time, G, GV, US, param_file, diag, CS%ZB2020, CS%use_ZB2020)
2529
25301 CS%initialized = .true.
2531
25321 CS%diag => diag
2533 ! Read parameters and write them to the model log.
25341 call log_version(param_file, mdl, version, "")
2535
2536 call get_param(param_file, mdl, "USE_CIRCULATION_IN_HORVISC", CS%use_circulation, &
2537 "Use circulation theorem to compute vorticity in horvisc module (for ZB20 or Leith)", &
25381 default=.False.)
2539
2540 ! All parameters are read in all cases to enable parameter spelling checks.
2541 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
2542 "This sets the default value for the various _ANSWER_DATE parameters.", &
25431 default=99991231)
2544
2545 ! Determine whether HOR_VISC_ANSWER_DATE is used, and avoid logging it if it is not used.
2546 call get_param(param_file, mdl, "USE_MEKE", use_MEKE, &
25471 default=.false., do_not_log=.true.)
25481 backscatter_Ro_c = 0.0
25491 if (use_MEKE) call get_param(param_file, mdl, "MEKE_BACKSCAT_RO_C", backscatter_Ro_c, &
2550 "The coefficient in the Rossby number function for scaling the biharmonic "//&
2551 "frictional energy source. Setting to non-zero enables the Rossby number function.", &
25521 units="nondim", default=0.0, do_not_log=.true.)
2553
2554 call get_param(param_file, mdl, "HOR_VISC_ANSWER_DATE", CS%answer_date, &
2555 "The vintage of the order of arithmetic and expressions in the horizontal "//&
2556 "viscosity calculations. Values between 20190102 and 20241201 recover the "//&
2557 "answers from the end of 2018, while higher values use updated and more robust "//&
2558 "forms of the same expressions.", &
25591 default=default_answer_date, do_not_log=(.not.GV%Boussinesq).or.(backscatter_Ro_c==0.0))
25601 if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20241201)
2561
25621 call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.)
2563 call get_param(param_file, mdl, "USE_CONT_THICKNESS", CS%use_cont_thick, &
2564 "If true, use thickness at velocity points from continuity solver. This option "//&
25651 "currently only works with split mode.", default=.false.)
2566 call get_param(param_file, mdl, "USE_CONT_THICKNESS_BUG", CS%use_cont_thick_bug, &
2567 "If true, retain an answer-changing halo update bug when "//&
2568 "USE_CONT_THICKNESS=True. This is not recommended.", &
25691 default=.false., do_not_log=.not.CS%use_cont_thick)
2570
2571 call get_param(param_file, mdl, "LAPLACIAN", CS%Laplacian, &
2572 "If true, use a Laplacian horizontal viscosity.", &
25731 default=.false.)
2574
2575 call get_param(param_file, mdl, "KH", Kh, &
2576 "The background Laplacian horizontal viscosity.", &
2577 units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s, &
25781 do_not_log=.not.CS%Laplacian)
2579 call get_param(param_file, mdl, "KH_BG_MIN", CS%Kh_bg_min, &
2580 "The minimum value allowed for Laplacian horizontal viscosity, KH.", &
2581 units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s, &
25821 do_not_log=.not.CS%Laplacian)
2583 call get_param(param_file, mdl, "KH_VEL_SCALE", Kh_vel_scale, &
2584 "The velocity scale which is multiplied by the grid "//&
2585 "spacing to calculate the Laplacian viscosity. "//&
2586 "The final viscosity is the largest of this scaled "//&
2587 "viscosity, the Smagorinsky and Leith viscosities, and KH.", &
2588 units="m s-1", default=0.0, scale=US%m_s_to_L_T, &
25891 do_not_log=.not.CS%Laplacian)
2590 call get_param(param_file, mdl, "KH_SIN_LAT", Kh_sin_lat, &
2591 "The amplitude of a latitudinally-dependent background "//&
2592 "viscosity of the form KH_SIN_LAT*(SIN(LAT)**KH_PWR_OF_SINE).", &
2593 units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s, &
25941 do_not_log=.not.CS%Laplacian)
2595 call get_param(param_file, mdl, "KH_PWR_OF_SINE", Kh_pwr_of_sine, &
2596 "The power used to raise SIN(LAT) when using a latitudinally "//&
2597 "dependent background viscosity.", &
2598 units="nondim", default=4.0, &
25991 do_not_log=.not.(CS%Laplacian .and. (Kh_sin_lat>0.)) )
2600 call get_param(param_file, mdl, "SMAGORINSKY_KH", CS%Smagorinsky_Kh, &
2601 "If true, use a Smagorinsky nonlinear eddy viscosity.", &
26021 default=.false., do_not_log=.not.CS%Laplacian)
26031 if (.not.CS%Laplacian) CS%Smagorinsky_Kh = .false.
2604 call get_param(param_file, mdl, "SMAG_LAP_CONST", Smag_Lap_const, &
2605 "The nondimensional Laplacian Smagorinsky constant, "//&
2606 "often 0.15.", units="nondim", default=0.0, &
26071 fail_if_missing=CS%Smagorinsky_Kh, do_not_log=.not.CS%Smagorinsky_Kh)
2608 call get_param(param_file, mdl, "LEITH_KH", CS%Leith_Kh, &
2609 "If true, use a Leith nonlinear eddy viscosity.", &
26101 default=.false., do_not_log=.not.CS%Laplacian)
26111 if (.not.CS%Laplacian) CS%Leith_Kh = .false.
2612 call get_param(param_file, mdl, "LEITH_LAP_CONST", Leith_Lap_const, &
2613 "The nondimensional Laplacian Leith constant, "//&
2614 "often set to 1.0", units="nondim", default=0.0, &
26151 fail_if_missing=CS%Leith_Kh, do_not_log=.not.CS%Leith_Kh)
2616 call get_param(param_file, mdl, "RES_SCALE_MEKE_VISC", CS%res_scale_MEKE, &
2617 "If true, the viscosity contribution from MEKE is scaled by "//&
2618 "the resolution function.", default=.false., &
26191 do_not_log=.not.(CS%Laplacian.and.use_MEKE))
26201 if (.not.(CS%Laplacian.and.use_MEKE)) CS%res_scale_MEKE = .false.
2621
2622 call get_param(param_file, mdl, "BOUND_KH", CS%bound_Kh, &
2623 "If true, the Laplacian coefficient is locally limited "//&
26241 "to be stable.", default=.true., do_not_log=.not.CS%Laplacian)
2625 call get_param(param_file, mdl, "EY24_EBT_BS", CS%EY24_EBT_BS, &
2626 "If true, use the backscatter scheme (EBT mode with kill switch) "//&
2627 "developed by Yankovsky et al. (2024). ", &
26281 default=.false., do_not_log=.not.CS%Laplacian)
26291 if (.not.CS%Laplacian) CS%bound_Kh = .false.
26301 if (.not.(CS%Laplacian.and.use_MEKE)) CS%EY24_EBT_BS = .false.
2631 call get_param(param_file, mdl, "ANISOTROPIC_VISCOSITY", CS%anisotropic, &
2632 "If true, allow anistropic viscosity in the Laplacian "//&
2633 "horizontal viscosity.", default=.false., &
26341 do_not_log=.not.CS%Laplacian)
26351 if (.not.CS%Laplacian) CS%anisotropic = .false. ! This replicates the prior code, but is it intended?
2636 call get_param(param_file, mdl, "ADD_LES_VISCOSITY", CS%add_LES_viscosity, &
2637 "If true, adds the viscosity from Smagorinsky and Leith to the "//&
2638 "background viscosity instead of taking the maximum.", default=.false., &
26391 do_not_log=.not.CS%Laplacian)
2640
2641 call get_param(param_file, mdl, "KH_ANISO", CS%Kh_aniso, &
2642 "The background Laplacian anisotropic horizontal viscosity.", &
2643 units="m2 s-1", default=0.0, scale=US%m_to_L**2*US%T_to_s, &
26441 do_not_log=.not.CS%anisotropic)
2645 call get_param(param_file, mdl, "ANISOTROPIC_MODE", aniso_mode, &
2646 "Selects the mode for setting the direction of anisotropy.\n"//&
2647 "\t 0 - Points along the grid i-direction.\n"//&
2648 "\t 1 - Points towards East.\n"//&
2649 "\t 2 - Points along the flow direction, U/|U|.", &
26501 default=0, do_not_log=.not.CS%anisotropic)
26511 if (aniso_mode == 0) then
2652 call get_param(param_file, mdl, "ANISO_GRID_DIR", aniso_grid_dir, &
2653 "The vector pointing in the direction of anisotropy for horizontal viscosity. "//&
2654 "n1,n2 are the i,j components relative to the grid.", &
26551 units="nondim", fail_if_missing=CS%anisotropic, do_not_log=.not.CS%anisotropic)
26560 elseif (aniso_mode == 1) then
2657 call get_param(param_file, mdl, "ANISO_GRID_DIR", aniso_grid_dir, &
2658 "The vector pointing in the direction of anisotropy for horizontal viscosity. "//&
2659 "n1,n2 are the i,j components relative to the spherical coordinates.", &
26600 units="nondim", fail_if_missing=CS%anisotropic, do_not_log=.not.CS%anisotropic)
2661 else
2662 call get_param(param_file, mdl, "ANISO_GRID_DIR", aniso_grid_dir, &
2663 "The vector pointing in the direction of anisotropy for horizontal viscosity.", &
26640 units="nondim", fail_if_missing=.false., do_not_log=.true.)
2665 endif
2666
2667 call get_param(param_file, mdl, "BIHARMONIC", CS%biharmonic, &
2668 "If true, use a biharmonic horizontal viscosity. "//&
2669 "BIHARMONIC may be used with LAPLACIAN.", &
26701 default=.true.)
2671 call get_param(param_file, mdl, "AH", Ah, &
2672 "The background biharmonic horizontal viscosity.", &
2673 units="m4 s-1", default=0.0, scale=US%m_to_L**4*US%T_to_s, &
26741 do_not_log=.not.CS%biharmonic)
2675 call get_param(param_file, mdl, "AH_VEL_SCALE", Ah_vel_scale, &
2676 "The velocity scale which is multiplied by the cube of "//&
2677 "the grid spacing to calculate the biharmonic viscosity. "//&
2678 "The final viscosity is the largest of this scaled "//&
2679 "viscosity, the Smagorinsky and Leith viscosities, and AH.", &
26801 units="m s-1", default=0.0, scale=US%m_s_to_L_T, do_not_log=.not.CS%biharmonic)
2681 call get_param(param_file, mdl, "AH_TIME_SCALE", Ah_time_scale, &
2682 "A time scale whose inverse is multiplied by the fourth "//&
2683 "power of the grid spacing to calculate biharmonic viscosity. "//&
2684 "The final viscosity is the largest of all viscosity "//&
2685 "formulations in use. 0.0 means that it's not used.", &
26861 units="s", default=0.0, scale=US%s_to_T, do_not_log=.not.CS%biharmonic)
2687 call get_param(param_file, mdl, "SMAGORINSKY_AH", CS%Smagorinsky_Ah, &
2688 "If true, use a biharmonic Smagorinsky nonlinear eddy "//&
26891 "viscosity.", default=.false., do_not_log=.not.CS%biharmonic)
26901 if (.not.CS%biharmonic) CS%Smagorinsky_Ah = .false.
2691 call get_param(param_file, mdl, "LEITH_AH", CS%Leith_Ah, &
2692 "If true, use a biharmonic Leith nonlinear eddy "//&
26931 "viscosity.", default=.false., do_not_log=.not.CS%biharmonic)
26941 if (.not.CS%biharmonic) CS%Leith_Ah = .false.
2695 call get_param(param_file, mdl, "USE_LEITHY", CS%use_Leithy, &
2696 "If true, use a biharmonic Leith nonlinear eddy "//&
2697 "viscosity together with a harmonic backscatter.", &
26981 default=.false.)
2699 call get_param(param_file, mdl, "BOUND_AH", CS%bound_Ah, &
2700 "If true, the biharmonic coefficient is locally limited "//&
27011 "to be stable.", default=.true., do_not_log=.not.CS%biharmonic)
27021 if (.not.CS%biharmonic) CS%bound_Ah = .false.
2703 call get_param(param_file, mdl, "RE_AH", CS%Re_Ah, &
2704 "If nonzero, the biharmonic coefficient is scaled "//&
2705 "so that the biharmonic Reynolds number is equal to this.", &
27061 units="nondim", default=0.0, do_not_log=.not.CS%biharmonic)
2707
2708 call get_param(param_file, mdl, "BACKSCATTER_UNDERBOUND", CS%backscatter_underbound, &
2709 "If true, the bounds on the biharmonic viscosity are allowed to "//&
2710 "increase where the Laplacian viscosity is negative (due to backscatter "//&
2711 "parameterizations) beyond the largest timestep-dependent stable values of "//&
2712 "biharmonic viscosity when no Laplacian viscosity is applied. The default "//&
2713 "is true for historical reasons, but this option probably should not be used "//&
2714 "because it can contribute to numerical instabilities.", &
27151 default=.false., do_not_log=.not.((CS%bound_Kh).and.(CS%bound_Ah)))
2716
2717 call get_param(param_file, mdl, "SMAG_BI_CONST",Smag_bi_const, &
2718 "The nondimensional biharmonic Smagorinsky constant, "//&
2719 "typically 0.015 - 0.06.", units="nondim", default=0.0, &
27201 fail_if_missing=CS%Smagorinsky_Ah, do_not_log=.not.CS%Smagorinsky_Ah)
2721
2722 call get_param(param_file, mdl, "USE_BETA_IN_LEITH", CS%use_beta_in_Leith, &
2723 "If true, include the beta term in the Leith nonlinear eddy viscosity.", &
27241 default=CS%Leith_Kh, do_not_log=.not.(CS%Leith_Kh .or. CS%Leith_Ah) )
2725 call get_param(param_file, mdl, "MODIFIED_LEITH", CS%modified_Leith, &
2726 "If true, add a term to Leith viscosity which is "//&
2727 "proportional to the gradient of divergence.", &
27281 default=.false., do_not_log=.not.(CS%Leith_Kh .or. CS%Leith_Ah) )
2729 call get_param(param_file, mdl, "USE_QG_LEITH_VISC", CS%use_QG_Leith_visc, &
2730 "If true, use QG Leith nonlinear eddy viscosity.", &
27311 default=.false., do_not_log=.not.(CS%Leith_Kh .or. CS%Leith_Ah) )
2732! if (CS%use_QG_Leith_visc) then
2733! call MOM_error(FATAL, "USE_QG_LEITH_VISC=True activates code that is a work-in-progress and "//&
2734! "should not be used until a number of bugs are fixed. Specifically it does not "//&
2735! "reproduce across PE count or layout, and may use arrays that have not been properly "//&
2736! "set or allocated. See github.com/mom-ocean/MOM6/issues/1590 for a discussion.")
2737! endif
27381 if (CS%use_QG_Leith_visc .and. .not. (CS%Leith_Kh .or. CS%Leith_Ah) ) then
2739 call MOM_error(FATAL, "MOM_hor_visc.F90, hor_visc_init: "//&
27400 "LEITH_KH or LEITH_AH must be True when USE_QG_LEITH_VISC=True.")
2741 endif
2742
27431 call get_param(param_file, mdl, "BOUND_CORIOLIS", bound_Cor_def, default=.false.)
2744 call get_param(param_file, mdl, "BOUND_CORIOLIS_BIHARM", CS%bound_Coriolis, &
2745 "If true use a viscosity that increases with the square "//&
2746 "of the velocity shears, so that the resulting viscous "//&
2747 "drag is of comparable magnitude to the Coriolis terms "//&
2748 "when the velocity differences between adjacent grid "//&
2749 "points is 0.5*BOUND_CORIOLIS_VEL. The default is the "//&
2750 "value of BOUND_CORIOLIS (or false).", default=bound_Cor_def, &
27511 do_not_log=.not.CS%Smagorinsky_Ah)
27521 if (.not.CS%Smagorinsky_Ah) CS%bound_Coriolis = .false.
2753 call get_param(param_file, mdl, "MAXVEL", maxvel, &
27541 units="m s-1", default=3.0e8, scale=US%m_s_to_L_T)
2755 call get_param(param_file, mdl, "BOUND_CORIOLIS_VEL", bound_Cor_vel, &
2756 "The velocity scale at which BOUND_CORIOLIS_BIHARM causes "//&
2757 "the biharmonic drag to have comparable magnitude to the "//&
2758 "Coriolis acceleration. The default is set by MAXVEL.", &
2759 units="m s-1", default=maxvel*US%L_T_to_m_s, scale=US%m_s_to_L_T, &
27601 do_not_log=.not.(CS%Smagorinsky_Ah .and. CS%bound_Coriolis))
2761 call get_param(param_file, mdl, "LEITH_BI_CONST", Leith_bi_const, &
2762 "The nondimensional biharmonic Leith constant, "//&
2763 "typical values are thus far undetermined.", units="nondim", default=0.0, &
2764 fail_if_missing=(CS%Leith_Ah .or. CS%use_Leithy), &
27651 do_not_log=.not.(CS%Leith_Ah .or. CS%use_Leithy))
2766 call get_param(param_file, mdl, "USE_LAND_MASK_FOR_HVISC", CS%use_land_mask, &
2767 "If true, use the land mask for the computation of thicknesses "//&
2768 "at velocity locations. This eliminates the dependence on arbitrary "//&
27691 "values over land or outside of the domain.", default=.true.)
2770 call get_param(param_file, mdl, "HORVISC_BOUND_COEF", CS%bound_coef, &
2771 "The nondimensional coefficient of the ratio of the "//&
2772 "viscosity bounds to the theoretical maximum for "//&
2773 "stability without considering other terms.", units="nondim", &
27741 default=0.8, do_not_log=.not.(CS%bound_Ah .or. CS%bound_Kh))
2775 call get_param(param_file, mdl, "KILL_SWITCH_COEF", CS%KS_coef, &
2776 "A nondimensional coefficient on the biharmonic viscosity that "// &
2777 "sets the kill switch for backscatter. Default is 1.0.", units="nondim", &
27781 default=1.0, do_not_log=.not.(CS%EY24_EBT_BS))
2779 call get_param(param_file, mdl, "NOSLIP", CS%no_slip, &
2780 "If true, no slip boundary conditions are used; otherwise "//&
2781 "free slip boundary conditions are assumed. The "//&
2782 "implementation of the free slip BCs on a C-grid is much "//&
2783 "cleaner than the no slip BCs. The use of free slip BCs "//&
2784 "is strongly encouraged, and no slip BCs are not used with "//&
27851 "the biharmonic viscosity.", default=.false.)
2786 call get_param(param_file, mdl, "USE_KH_BG_2D", CS%use_Kh_bg_2d, &
2787 "If true, read a file containing 2-d background harmonic "//&
2788 "viscosities. The final viscosity is the maximum of the other "//&
27891 "terms and this background value.", default=.false., do_not_log=.not.CS%Laplacian)
27901 if (.not.CS%Laplacian) CS%use_Kh_bg_2d = .false.
2791 call get_param(param_file, mdl, "KH_BG_2D_BUG", CS%Kh_bg_2d_bug, &
2792 "If true, retain an answer-changing horizontal indexing bug in setting "//&
2793 "the corner-point viscosities when USE_KH_BG_2D=True. This is "//&
27941 "not recommended.", default=.false., do_not_log=.not.CS%use_Kh_bg_2d)
2795 call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
27961 default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
2797 call get_param(param_file, mdl, "FRICTWORK_BUG", CS%FrictWork_bug, &
2798 "If true, retain an answer-changing bug in calculating the FrictWork, "//&
2799 "which cancels the h in thickness flux and the h at velocity point. This is "//&
28001 "not recommended.", default=.false.)
2801 call get_param(param_file, mdl, "OBC_SPECIFIED_STRAIN_BUG", CS%OBC_strain_bug, &
2802 "If true, recover a bug that specified shear strain option at open boundaries "//&
28031 "cannot be applied.", default=.true.)
2804 call get_param(param_file, mdl, "USE_GME", CS%use_GME, &
2805 "If true, use the GM+E backscatter scheme in association \n"//&
28061 "with the Gent and McWilliams parameterization.", default=.false.)
2807 call get_param(param_file, mdl, "SPLIT", split, &
28081 "Use the split time stepping if true.", default=.true., do_not_log=.true.)
28091 if (CS%use_Leithy) then
28100 if (.not.(CS%biharmonic .and. CS%Laplacian)) then
2811 call MOM_error(FATAL, "MOM_hor_visc.F90, hor_visc_init: "//&
28120 "LAPLACIAN and BIHARMONIC must both be True when USE_LEITHY=True.")
2813 endif
2814 endif
2815 call get_param(param_file, mdl, "LEITHY_CK", CS%c_K, &
2816 "Fraction of biharmonic dissipation that gets backscattered, "//&
28171 "in Leith+E.", units="nondim", default=1.0, do_not_log=.not.CS%use_Leithy)
2818 call get_param(param_file, mdl, "SMOOTH_AH", CS%smooth_Ah, &
2819 "If true, Ah and m_leithy are smoothed within Leith+E. This requires "//&
2820 "lots of blocking communications, which can be expensive", &
28211 default=.true., do_not_log=.not.CS%use_Leithy)
2822
28231 if (CS%use_GME .and. .not.split) call MOM_error(FATAL,"ERROR: Currently, USE_GME = True "// &
28240 "cannot be used with SPLIT=False.")
2825
28261 if (CS%use_GME) then
2827 call get_param(param_file, mdl, "GME_NUM_SMOOTHINGS", CS%num_smooth_gme, &
2828 "Number of smoothing passes for the GME fluxes.", &
28290 default=1)
2830 call get_param(param_file, mdl, "GME_H0", CS%GME_h0, &
2831 "The strength of GME tapers quadratically to zero when the bathymetric "//&
2832 "depth is shallower than GME_H0.", &
28330 units="m", scale=GV%m_to_H, default=1000.0)
2834 call get_param(param_file, mdl, "GME_EFFICIENCY", CS%GME_efficiency, &
2835 "The nondimensional prefactor multiplying the GME coefficient.", &
28360 units="nondim", default=1.0)
2837 call get_param(param_file, mdl, "GME_LIMITER", CS%GME_limiter, &
2838 "The absolute maximum value the GME coefficient is allowed to take.", &
28390 units="m2 s-1", scale=US%m_to_L**2*US%T_to_s, default=1.0e7)
2840 endif
2841
28421 if (CS%Laplacian .or. CS%biharmonic) then
2843 call get_param(param_file, mdl, "DT", dt, &
2844 "The (baroclinic) dynamics time step.", units="s", scale=US%s_to_T, &
28451 fail_if_missing=.true.)
28461 Idt = 1.0 / dt
2847 endif
2848 call get_param(param_file, mdl, "KILL_SWITCH_TIMESCALE", CS%KS_timescale, &
2849 "A timescale for computing the CFL limit for viscosity "// &
2850 "that determines when backscatter is shut off. Default is DT.", &
28511 default= dt , units="s", scale=US%s_to_T, do_not_log=.not.(CS%EY24_EBT_BS))
2852
28531 if (CS%no_slip .and. CS%biharmonic) &
2854 call MOM_error(FATAL,"ERROR: NOSLIP and BIHARMONIC cannot be defined "// &
28550 "at the same time in MOM.")
28561 if (.not.(CS%Laplacian .or. CS%biharmonic)) then
2857 ! Only issue inviscid warning if not in single column mode (usually 2x2 domain)
28580 if ( max(G%domain%niglobal, G%domain%njglobal)>2 ) call MOM_error(WARNING, &
2859 "hor_visc_init: It is usually a very bad idea not to use either "//&
28600 "LAPLACIAN or BIHARMONIC viscosity.")
28610 return ! We are not using either Laplacian or Bi-harmonic lateral viscosity
2862 endif
2863
2864 !$omp target update to(CS)
2865
28661 deg2rad = atan(1.0) / 45.
28678773 ALLOC_(CS%dx2h(isd:ied,jsd:jed)) ; CS%dx2h(:,:) = 0.0
28688773 ALLOC_(CS%dy2h(isd:ied,jsd:jed)) ; CS%dy2h(:,:) = 0.0
28698971 ALLOC_(CS%dx2q(IsdB:IedB,JsdB:JedB)) ; CS%dx2q(:,:) = 0.0
28708971 ALLOC_(CS%dy2q(IsdB:IedB,JsdB:JedB)) ; CS%dy2q(:,:) = 0.0
28718773 ALLOC_(CS%dx_dyT(isd:ied,jsd:jed)) ; CS%dx_dyT(:,:) = 0.0
28728773 ALLOC_(CS%dy_dxT(isd:ied,jsd:jed)) ; CS%dy_dxT(:,:) = 0.0
28738971 ALLOC_(CS%dx_dyBu(IsdB:IedB,JsdB:JedB)) ; CS%dx_dyBu(:,:) = 0.0
28748971 ALLOC_(CS%dy_dxBu(IsdB:IedB,JsdB:JedB)) ; CS%dy_dxBu(:,:) = 0.0
2875 !$omp target enter data map(alloc: CS%dx2h, CS%dy2h, CS%dx2q, CS%dy2q)
2876 !$omp target enter data map(alloc: CS%dx_dyT, CS%dy_dxT, CS%dx_dyBu, CS%dy_dxBu)
2877
28781 if (CS%Laplacian) then
28798773 ALLOC_(CS%grid_sp_h2(isd:ied,jsd:jed)) ; CS%grid_sp_h2(:,:) = 0.0
28808773 ALLOC_(CS%Kh_bg_xx(isd:ied,jsd:jed)) ; CS%Kh_bg_xx(:,:) = 0.0
28818971 ALLOC_(CS%Kh_bg_xy(IsdB:IedB,JsdB:JedB)) ; CS%Kh_bg_xy(:,:) = 0.0
28821 if (CS%bound_Kh .or. CS%EY24_EBT_BS) then
28838773 allocate(CS%Kh_Max_xx(Isd:Ied,Jsd:Jed), source=0.0)
28848971 allocate(CS%Kh_Max_xy(IsdB:IedB,JsdB:JedB), source=0.0)
2885 endif
28861 if (CS%Smagorinsky_Kh .or. CS%EY24_EBT_BS) then
28870 allocate(CS%Laplac2_const_xx(isd:ied,jsd:jed), source=0.0)
28880 allocate(CS%Laplac2_const_xy(IsdB:IedB,JsdB:JedB), source=0.0)
2889 endif
28901 if (CS%Leith_Kh) then
28910 allocate(CS%Laplac3_const_xx(isd:ied,jsd:jed), source=0.0)
28920 allocate(CS%Laplac3_const_xy(IsdB:IedB,JsdB:JedB), source=0.0)
2893 endif
2894 endif
28958773 ALLOC_(CS%reduction_xx(isd:ied,jsd:jed)) ; CS%reduction_xx(:,:) = 0.0
28968971 ALLOC_(CS%reduction_xy(IsdB:IedB,JsdB:JedB)) ; CS%reduction_xy(:,:) = 0.0
2897
28981 CS%dynamic_aniso = .false.
28991 if (CS%anisotropic) then
29000 allocate(CS%n1n2_h(isd:ied,jsd:jed), source=0.0)
29010 allocate(CS%n1n1_m_n2n2_h(isd:ied,jsd:jed), source=0.0)
29020 allocate(CS%n1n2_q(IsdB:IedB,JsdB:JedB), source=0.0)
29030 allocate(CS%n1n1_m_n2n2_q(IsdB:IedB,JsdB:JedB), source=0.0)
29040 select case (aniso_mode)
2905 case (0)
29060 call align_aniso_tensor_to_grid(CS, aniso_grid_dir(1), aniso_grid_dir(2))
2907 case (1)
2908 ! call align_aniso_tensor_to_grid(CS, aniso_grid_dir(1), aniso_grid_dir(2))
2909 case (2)
29100 CS%dynamic_aniso = .true.
2911 case default
2912 call MOM_error(FATAL, "MOM_hor_visc: "//&
29130 "Runtime parameter ANISOTROPIC_MODE is out of range.")
2914 end select
2915 endif
2916
2917 call get_param(param_file, mdl, "KH_BG_2D_FILENAME", filename, &
2918 'The filename containing a 2d map of "Kh".', &
29191 default='KH_background_2d.nc', do_not_log=.not.CS%use_Kh_bg_2d)
2920 call get_param(param_file, mdl, "KH_BG_2D_VARNAME", Kh_var, &
2921 'The name in the input file of the horizontal viscosity variable.', &
29221 default='Kh', do_not_log=.not.CS%use_Kh_bg_2d)
2923
29241 if (CS%use_Kh_bg_2d) then
29250 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
29260 inputdir = slasher(inputdir)
29270 allocate(CS%Kh_bg_2d(isd:ied,jsd:jed), source=0.0)
2928 call MOM_read_data(trim(inputdir)//trim(filename), Kh_var, CS%Kh_bg_2d, &
29290 G%domain, timelevel=1, scale=US%m_to_L**2*US%T_to_s)
29300 call pass_var(CS%Kh_bg_2d, G%domain)
2931 endif
29321 if (CS%biharmonic) then
29338841 ALLOC_(CS%Idx2dyCu(IsdB:IedB,jsd:jed)) ; CS%Idx2dyCu(:,:) = 0.0
29348902 ALLOC_(CS%Idx2dyCv(isd:ied,JsdB:JedB)) ; CS%Idx2dyCv(:,:) = 0.0
29358841 ALLOC_(CS%Idxdy2u(IsdB:IedB,jsd:jed)) ; CS%Idxdy2u(:,:) = 0.0
29368902 ALLOC_(CS%Idxdy2v(isd:ied,JsdB:JedB)) ; CS%Idxdy2v(:,:) = 0.0
29378773 ALLOC_(CS%Ah_bg_xx(isd:ied,jsd:jed)) ; CS%Ah_bg_xx(:,:) = 0.0
29388971 ALLOC_(CS%Ah_bg_xy(IsdB:IedB,JsdB:JedB)) ; CS%Ah_bg_xy(:,:) = 0.0
29398773 ALLOC_(CS%grid_sp_h3(isd:ied,jsd:jed)) ; CS%grid_sp_h3(:,:) = 0.0
29401 if (CS%bound_Ah) then
29418773 allocate(CS%Ah_Max_xx(isd:ied,jsd:jed), source=0.0)
29428971 allocate(CS%Ah_Max_xy(IsdB:IedB,JsdB:JedB), source=0.0)
2943 endif
29441 if (CS%EY24_EBT_BS) then
29450 allocate(CS%Ah_Max_xx_KS(isd:ied,jsd:jed), source=0.0)
29460 allocate(CS%Ah_Max_xy_KS(IsdB:IedB,JsdB:JedB), source=0.0)
2947 endif
29481 if (CS%Smagorinsky_Ah) then
29498773 allocate(CS%Biharm_const_xx(isd:ied,jsd:jed), source=0.0)
29508971 allocate(CS%Biharm_const_xy(IsdB:IedB,JsdB:JedB), source=0.0)
29511 if (CS%bound_Coriolis) then
29528773 allocate(CS%Biharm_const2_xx(isd:ied,jsd:jed), source=0.0)
29538971 allocate(CS%Biharm_const2_xy(IsdB:IedB,JsdB:JedB), source=0.0)
2954 endif
2955 endif
29561 if ((CS%Leith_Ah) .or. (CS%use_Leithy)) then
29570 allocate(CS%biharm6_const_xx(isd:ied,jsd:jed), source=0.0)
29580 allocate(CS%biharm6_const_xy(IsdB:IedB,JsdB:JedB), source=0.0)
2959 endif
29601 if (CS%use_Leithy) then
29610 allocate(CS%m_const_leithy(isd:ied,jsd:jed), source=0.0)
29620 allocate(CS%m_leithy_max(isd:ied,jsd:jed), source=0.0)
2963 endif
29641 if (CS%Re_Ah > 0.0) then
29650 allocate(CS%Re_Ah_const_xx(isd:ied,jsd:jed), source=0.0)
29660 allocate(CS%Re_Ah_const_xy(IsdB:IedB,JsdB:JedB), source=0.0)
2967 endif
2968 endif
2969
2970124 do concurrent (J=js-2:Jeq+1, I=is-2:Ieq+1)
29717873 CS%dx2q(I,J) = G%dxBu(I,J)*G%dxBu(I,J) ; CS%dy2q(I,J) = G%dyBu(I,J)*G%dyBu(I,J)
2972 enddo
2973
29741 if (((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) .and. &
2975 ((G%isc-G%isd < 3) .or. (G%isc-G%isd < 3))) call MOM_error(FATAL, &
29760 "The minimum halo size is 3 when a Leith viscosity is being used.")
29771 if (CS%use_Leithy) then
29780 do concurrent (J=js-3:Jeq+2, I=is-3:Ieq+2)
29790 CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J)
2980 enddo
29811 elseif ((CS%Leith_Kh) .or. (CS%Leith_Ah)) then
29820 do concurrent (J=Jsq-2:Jeq+2, I=Isq-2:Ieq+2)
29830 CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J)
2984 enddo
2985 else
2986124 do concurrent (J=js-2:Jeq+1, I=is-2:Ieq+1)
29877873 CS%DX_dyBu(I,J) = G%dxBu(I,J)*G%IdyBu(I,J) ; CS%DY_dxBu(I,J) = G%dyBu(I,J)*G%IdxBu(I,J)
2988 enddo
2989 endif
2990
2991125 do concurrent (j=js-2:Jeq+2, i=is-2:Ieq+2)
29927936 CS%dx2h(i,j) = G%dxT(i,j)*G%dxT(i,j) ; CS%dy2h(i,j) = G%dyT(i,j)*G%dyT(i,j)
29938061 CS%DX_dyT(i,j) = G%dxT(i,j)*G%IdyT(i,j) ; CS%DY_dxT(i,j) = G%dyT(i,j)*G%IdxT(i,j)
2994 enddo
2995
2996 ! TODO: Remove this after every instance has been moved to GPU
2997 !$omp target update from(CS%dx2q, CS%dy2q, CS%dx_dyBu, CS%dy_dxBu)
2998 !$omp target update from(CS%dx2h, CS%dy2h, CS%dx_dyT, CS%dy_dxT)
2999
30007627 do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
30017564 CS%reduction_xx(i,j) = 1.0
30027564 if ((G%dy_Cu(I,j) > 0.0) .and. (G%dy_Cu(I,j) < G%dyCu(I,j)) .and. &
3003 (G%dy_Cu(I,j) < G%dyCu(I,j) * CS%reduction_xx(i,j))) &
30040 CS%reduction_xx(i,j) = G%dy_Cu(I,j) / (G%dyCu(I,j))
30057564 if ((G%dy_Cu(I-1,j) > 0.0) .and. (G%dy_Cu(I-1,j) < G%dyCu(I-1,j)) .and. &
3006 (G%dy_Cu(I-1,j) < G%dyCu(I-1,j) * CS%reduction_xx(i,j))) &
30070 CS%reduction_xx(i,j) = G%dy_Cu(I-1,j) / (G%dyCu(I-1,j))
30087564 if ((G%dx_Cv(i,J) > 0.0) .and. (G%dx_Cv(i,J) < G%dxCv(i,J)) .and. &
3009 (G%dx_Cv(i,J) < G%dxCv(i,J) * CS%reduction_xx(i,j))) &
30100 CS%reduction_xx(i,j) = G%dx_Cv(i,J) / (G%dxCv(i,J))
30117564 if ((G%dx_Cv(i,J-1) > 0.0) .and. (G%dx_Cv(i,J-1) < G%dxCv(i,J-1)) .and. &
3012 (G%dx_Cv(i,J-1) < G%dxCv(i,J-1) * CS%reduction_xx(i,j))) &
301362 CS%reduction_xx(i,j) = G%dx_Cv(i,J-1) / (G%dxCv(i,J-1))
3014 enddo ; enddo
3015
30167443 do J=js-1,Jeq ; do I=is-1,Ieq
30177381 CS%reduction_xy(I,J) = 1.0
30187381 if ((G%dy_Cu(I,j) > 0.0) .and. (G%dy_Cu(I,j) < G%dyCu(I,j)) .and. &
3019 (G%dy_Cu(I,j) < G%dyCu(I,j) * CS%reduction_xy(I,J))) &
30200 CS%reduction_xy(I,J) = G%dy_Cu(I,j) / (G%dyCu(I,j))
30217381 if ((G%dy_Cu(I,j+1) > 0.0) .and. (G%dy_Cu(I,j+1) < G%dyCu(I,j+1)) .and. &
3022 (G%dy_Cu(I,j+1) < G%dyCu(I,j+1) * CS%reduction_xy(I,J))) &
30230 CS%reduction_xy(I,J) = G%dy_Cu(I,j+1) / (G%dyCu(I,j+1))
30247381 if ((G%dx_Cv(i,J) > 0.0) .and. (G%dx_Cv(i,J) < G%dxCv(i,J)) .and. &
3025 (G%dx_Cv(i,J) < G%dxCv(i,J) * CS%reduction_xy(I,J))) &
30260 CS%reduction_xy(I,J) = G%dx_Cv(i,J) / (G%dxCv(i,J))
30277381 if ((G%dx_Cv(i+1,J) > 0.0) .and. (G%dx_Cv(i+1,J) < G%dxCv(i+1,J)) .and. &
3028 (G%dx_Cv(i+1,J) < G%dxCv(i+1,J) * CS%reduction_xy(I,J))) &
302961 CS%reduction_xy(I,J) = G%dx_Cv(i+1,J) / (G%dxCv(i+1,J))
3030 enddo ; enddo
30311 if (CS%Laplacian) then
3032 ! The 0.3 below was 0.4 in MOM1.10. The change in hq requires
3033 ! this to be less than 1/3, rather than 1/2 as before.
30341 if (CS%bound_Kh .or. CS%bound_Ah) Kh_Limit = 0.3 / (dt*4.0)
3035 ! Calculate and store the background viscosity at h-points
3036
30371 min_grid_sp_h2 = huge(1.)
30387627 do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
3039 ! Static factors in the Smagorinsky and Leith schemes
30407564 grid_sp_h2 = (2.0*CS%dx2h(i,j)*CS%dy2h(i,j)) / (CS%dx2h(i,j) + CS%dy2h(i,j))
30417564 CS%grid_sp_h2(i,j) = grid_sp_h2
30427564 grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
30437564 if (CS%Smagorinsky_Kh) CS%Laplac2_const_xx(i,j) = Smag_Lap_const * grid_sp_h2
30447564 if (CS%Leith_Kh) CS%Laplac3_const_xx(i,j) = Leith_Lap_const * grid_sp_h3
3045 ! Maximum of constant background and MICOM viscosity
30467564 CS%Kh_bg_xx(i,j) = MAX(Kh, Kh_vel_scale * sqrt(grid_sp_h2))
3047 ! Use the larger of the above and values read from a file
30487564 if (CS%use_Kh_bg_2d) CS%Kh_bg_xx(i,j) = MAX(CS%Kh_bg_2d(i,j), CS%Kh_bg_xx(i,j))
3049 ! Use the larger of the above and a function of sin(latitude)
30507564 if (Kh_sin_lat>0.) then
30510 slat_fn = abs( sin( deg2rad * G%geoLatT(i,j) ) ) ** Kh_pwr_of_sine
30520 CS%Kh_bg_xx(i,j) = MAX(Kh_sin_lat * slat_fn, CS%Kh_bg_xx(i,j))
3053 endif
30547626 min_grid_sp_h2 = min(grid_sp_h2, min_grid_sp_h2)
3055 enddo ; enddo
30561 call min_across_PEs(min_grid_sp_h2)
3057
3058 ! Calculate and store the background viscosity at q-points
30597443 do J=js-1,Jeq ; do I=is-1,Ieq
3060 ! Static factors in the Smagorinsky and Leith schemes
30617381 grid_sp_q2 = (2.0*CS%dx2q(I,J)*CS%dy2q(I,J)) / (CS%dx2q(I,J) + CS%dy2q(I,J))
30627381 grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)
30637381 if (CS%Smagorinsky_Kh) CS%Laplac2_const_xy(I,J) = Smag_Lap_const * grid_sp_q2
30647381 if (CS%Leith_Kh) CS%Laplac3_const_xy(I,J) = Leith_Lap_const * grid_sp_q3
3065 ! Maximum of constant background and MICOM viscosity
30667381 CS%Kh_bg_xy(I,J) = MAX(Kh, Kh_vel_scale * sqrt(grid_sp_q2))
3067 ! Use the larger of the above and values read from a file
30687381 if (CS%use_Kh_bg_2d) then
30690 if (CS%Kh_bg_2d_bug) then
3070 ! This option is unambiguously wrong but is needed to recover old answers
30710 CS%Kh_bg_xy(I,J) = MAX(CS%Kh_bg_2d(i,j), CS%Kh_bg_xy(I,J))
3072 else
3073 CS%Kh_bg_xy(I,J) = MAX(CS%Kh_bg_xy(I,J), &
3074 0.25*((CS%Kh_bg_2d(i,j) + CS%Kh_bg_2d(i+1,j+1)) + &
30750 (CS%Kh_bg_2d(i+1,j) + CS%Kh_bg_2d(i,j+1))) )
3076 endif
3077 endif
3078
3079 ! Use the larger of the above and a function of sin(latitude)
30807442 if (Kh_sin_lat>0.) then
30810 slat_fn = abs( sin( deg2rad * G%geoLatBu(I,J) ) ) ** Kh_pwr_of_sine
30820 CS%Kh_bg_xy(I,J) = MAX(Kh_sin_lat * slat_fn, CS%Kh_bg_xy(I,J))
3083 endif
3084 enddo ; enddo
3085 endif
30861 if (CS%biharmonic) then
30877689 do j=js-1,Jeq+1 ; do I=is-2,Ieq+1
30887626 CS%Idx2dyCu(I,j) = (G%IdxCu(I,j)*G%IdxCu(I,j)) * G%IdyCu(I,j)
30897688 CS%Idxdy2u(I,j) = G%IdxCu(I,j) * (G%IdyCu(I,j)*G%IdyCu(I,j))
3090 enddo ; enddo
30917750 do J=js-2,Jeq+1 ; do i=is-1,Ieq+1
30927686 CS%Idx2dyCv(i,J) = (G%IdxCv(i,J)*G%IdxCv(i,J)) * G%IdyCv(i,J)
30937749 CS%Idxdy2v(i,J) = G%IdxCv(i,J) * (G%IdyCv(i,J)*G%IdyCv(i,J))
3094 enddo ; enddo
30958971 CS%Ah_bg_xy(:,:) = 0.0
3096 ! The 0.3 below was 0.4 in HIM 1.10. The change in hq requires
3097 ! this to be less than 1/3, rather than 1/2 as before.
30981 if (CS%bound_Ah) Ah_Limit = 0.3 / (dt*64.0)
30991 if (CS%Smagorinsky_Ah .and. CS%bound_Coriolis) &
31001 BoundCorConst = 1.0 / (5.0*(bound_Cor_vel*bound_Cor_vel))
3101
31021 min_grid_sp_h4 = huge(1.)
31037627 do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
31047564 grid_sp_h2 = (2.0*CS%dx2h(i,j)*CS%dy2h(i,j)) / (CS%dx2h(i,j)+CS%dy2h(i,j))
31057564 grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)
31067564 CS%grid_sp_h3(i,j) = grid_sp_h3
31077564 if (CS%Smagorinsky_Ah) then
31087564 CS%Biharm_const_xx(i,j) = Smag_bi_const * (grid_sp_h2 * grid_sp_h2)
31097564 if (CS%bound_Coriolis) then
3110 fmax = MAX(abs(G%CoriolisBu(I-1,J-1)), abs(G%CoriolisBu(I,J-1)), &
31117564 abs(G%CoriolisBu(I-1,J)), abs(G%CoriolisBu(I,J)))
3112 CS%Biharm_const2_xx(i,j) = (grid_sp_h2 * grid_sp_h2 * grid_sp_h2) * &
31137564 (fmax * BoundCorConst)
3114 endif
3115 endif
31167564 if (CS%Leith_Ah) then
31170 CS%biharm6_const_xx(i,j) = Leith_bi_const * (grid_sp_h3 * grid_sp_h3)
3118 endif
31197564 if (CS%use_Leithy) then
31200 CS%biharm6_const_xx(i,j) = Leith_bi_const * max(G%dxT(i,j),G%dyT(i,j))**6
31210 CS%m_const_leithy(i,j) = 0.5 * sqrt(CS%c_K) * max(G%dxT(i,j),G%dyT(i,j))
31220 CS%m_leithy_max(i,j) = 4. / max(G%dxT(i,j),G%dyT(i,j))**2
3123 endif
31247564 CS%Ah_bg_xx(i,j) = MAX(Ah, Ah_vel_scale * grid_sp_h2 * sqrt(grid_sp_h2))
31257564 if (CS%Re_Ah > 0.0) CS%Re_Ah_const_xx(i,j) = grid_sp_h3 / CS%Re_Ah
31267564 if (Ah_time_scale > 0.) CS%Ah_bg_xx(i,j) = &
31270 MAX(CS%Ah_bg_xx(i,j), (grid_sp_h2 * grid_sp_h2) / Ah_time_scale)
31287626 min_grid_sp_h4 = min(grid_sp_h2**2, min_grid_sp_h4)
3129 enddo ; enddo
31301 call min_across_PEs(min_grid_sp_h4)
3131
31327443 do J=js-1,Jeq ; do I=is-1,Ieq
31337381 grid_sp_q2 = (2.0*CS%dx2q(I,J)*CS%dy2q(I,J)) / (CS%dx2q(I,J)+CS%dy2q(I,J))
31347381 grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)
31357381 if (CS%Smagorinsky_Ah) then
31367381 CS%Biharm_const_xy(I,J) = Smag_bi_const * (grid_sp_q2 * grid_sp_q2)
31377381 if (CS%bound_Coriolis) then
3138 CS%Biharm_const2_xy(I,J) = (grid_sp_q2 * grid_sp_q2 * grid_sp_q2) * &
31397381 (abs(G%CoriolisBu(I,J)) * BoundCorConst)
3140 endif
3141 endif
31427381 if ((CS%Leith_Ah) .or. (CS%use_Leithy))then
31430 CS%biharm6_const_xy(I,J) = Leith_bi_const * (grid_sp_q3 * grid_sp_q3)
3144 endif
31457381 CS%Ah_bg_xy(I,J) = MAX(Ah, Ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2))
31467381 if (CS%Re_Ah > 0.0) CS%Re_Ah_const_xy(i,j) = grid_sp_q3 / CS%Re_Ah
31477381 if (Ah_time_scale > 0.) CS%Ah_bg_xy(i,j) = &
314861 MAX(CS%Ah_bg_xy(i,j), (grid_sp_q2 * grid_sp_q2) / Ah_time_scale)
3149 enddo ; enddo
3150 endif
3151 ! The Laplacian bounds should avoid overshoots when CS%bound_coef < 1.
31521 if (CS%Laplacian .and. CS%bound_Kh) then
31537627 do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
3154 denom = max( &
3155 (CS%dy2h(i,j) * CS%DY_dxT(i,j) * (G%IdyCu(I,j) + G%IdyCu(I-1,j)) * &
3156 max(G%IdyCu(I,j)*G%IareaCu(I,j), G%IdyCu(I-1,j)*G%IareaCu(I-1,j)) ), &
3157 (CS%dx2h(i,j) * CS%DX_dyT(i,j) * (G%IdxCv(i,J) + G%IdxCv(i,J-1)) * &
31587564 max(G%IdxCv(i,J)*G%IareaCv(i,J), G%IdxCv(i,J-1)*G%IareaCv(i,J-1)) ) )
31597564 CS%Kh_Max_xx(i,j) = 0.0
31607564 if (denom > 0.0) &
31615112 CS%Kh_Max_xx(i,j) = CS%bound_coef * 0.25 * Idt / denom
3162 enddo ; enddo
31637443 do J=js-1,Jeq ; do I=is-1,Ieq
3164 denom = max( &
3165 (CS%dx2q(I,J) * CS%DX_dyBu(I,J) * (G%IdxCu(I,j+1) + G%IdxCu(I,j)) * &
3166 max(G%IdxCu(I,j)*G%IareaCu(I,j), G%IdxCu(I,j+1)*G%IareaCu(I,j+1)) ), &
3167 (CS%dy2q(I,J) * CS%DY_dxBu(I,J) * (G%IdyCv(i+1,J) + G%IdyCv(i,J)) * &
31687381 max(G%IdyCv(i,J)*G%IareaCv(i,J), G%IdyCv(i+1,J)*G%IareaCv(i+1,J)) ) )
31697381 CS%Kh_Max_xy(I,J) = 0.0
31707381 if (denom > 0.0) &
31715210 CS%Kh_Max_xy(I,J) = CS%bound_coef * 0.25 * Idt / denom
3172 enddo ; enddo
31731 if (CS%debug) then
31740 call hchksum(CS%Kh_Max_xx, "Kh_Max_xx", G%HI, haloshift=0, unscale=US%L_to_m**2*US%s_to_T)
31750 call Bchksum(CS%Kh_Max_xy, "Kh_Max_xy", G%HI, haloshift=0, unscale=US%L_to_m**2*US%s_to_T)
3176 endif
3177 endif
3178 ! The biharmonic bounds should avoid overshoots when CS%bound_coef < 0.5, but
3179 ! empirically work for CS%bound_coef <~ 1.0
31801 if (CS%biharmonic .and. CS%bound_Ah) then
31817689 do j=js-1,Jeq+1 ; do I=is-2,Ieq+1
3182 u0u(I,j) = ((CS%Idxdy2u(I,j)*((CS%dy2h(i+1,j)*CS%DY_dxT(i+1,j)*(G%IdyCu(I+1,j) + G%IdyCu(I,j))) + &
3183 (CS%dy2h(i,j) * CS%DY_dxT(i,j) * (G%IdyCu(I,j) + G%IdyCu(I-1,j))) )) + &
3184 (CS%Idx2dyCu(I,j)*((CS%dx2q(I,J) * CS%DX_dyBu(I,J) * (G%IdxCu(I,j+1) + G%IdxCu(I,j))) + &
31857626 (CS%dx2q(I,J-1)*CS%DX_dyBu(I,J-1)*(G%IdxCu(I,j) + G%IdxCu(I,j-1))) )) )
3186 u0v(I,j) = ((CS%Idxdy2u(I,j)*((CS%dy2h(i+1,j)*CS%DX_dyT(i+1,j)*(G%IdxCv(i+1,J) + G%IdxCv(i+1,J-1))) + &
3187 (CS%dy2h(i,j) * CS%DX_dyT(i,j) * (G%IdxCv(i,J) + G%IdxCv(i,J-1))) )) + &
3188 (CS%Idx2dyCu(I,j)*((CS%dx2q(I,J) * CS%DY_dxBu(I,J) * (G%IdyCv(i+1,J) + G%IdyCv(i,J))) + &
31897688 (CS%dx2q(I,J-1)*CS%DY_dxBu(I,J-1)*(G%IdyCv(i+1,J-1) + G%IdyCv(i,J-1))) )) )
3190 enddo ; enddo
31917750 do J=js-2,Jeq+1 ; do i=is-1,Ieq+1
3192 v0u(i,J) = ((CS%Idxdy2v(i,J)*((CS%dy2q(I,J) * CS%DX_dyBu(I,J) * (G%IdxCu(I,j+1) + G%IdxCu(I,j))) + &
3193 (CS%dy2q(I-1,J)*CS%DX_dyBu(I-1,J)*(G%IdxCu(I-1,j+1) + G%IdxCu(I-1,j))) )) + &
3194 (CS%Idx2dyCv(i,J)*((CS%dx2h(i,j+1)*CS%DY_dxT(i,j+1)*(G%IdyCu(I,j+1) + G%IdyCu(I-1,j+1))) + &
31957686 (CS%dx2h(i,j) * CS%DY_dxT(i,j) * (G%IdyCu(I,j) + G%IdyCu(I-1,j))) ) ))
3196 v0v(i,J) = ((CS%Idxdy2v(i,J)*((CS%dy2q(I,J) * CS%DY_dxBu(I,J) * (G%IdyCv(i+1,J) + G%IdyCv(i,J))) + &
3197 (CS%dy2q(I-1,J)*CS%DY_dxBu(I-1,J)*(G%IdyCv(i,J) + G%IdyCv(i-1,J))) )) + &
3198 (CS%Idx2dyCv(i,J)*((CS%dx2h(i,j+1)*CS%DX_dyT(i,j+1)*(G%IdxCv(i,J+1) + G%IdxCv(i,J))) + &
31997749 (CS%dx2h(i,j) * CS%DX_dyT(i,j) * (G%IdxCv(i,J) + G%IdxCv(i,J-1))) )) )
3200 enddo ; enddo
32017627 do j=js-1,Jeq+1 ; do i=is-1,Ieq+1
3202 denom = max( &
3203 (CS%dy2h(i,j) * &
3204 ((CS%DY_dxT(i,j)*((G%IdyCu(I,j)*u0u(I,j)) + (G%IdyCu(I-1,j)*u0u(I-1,j)))) + &
3205 (CS%DX_dyT(i,j)*((G%IdxCv(i,J)*v0u(i,J)) + (G%IdxCv(i,J-1)*v0u(i,J-1))))) * &
3206 max(G%IdyCu(I,j)*G%IareaCu(I,j), G%IdyCu(I-1,j)*G%IareaCu(I-1,j)) ), &
3207 (CS%dx2h(i,j) * &
3208 ((CS%DY_dxT(i,j)*((G%IdyCu(I,j)*u0v(I,j)) + (G%IdyCu(I-1,j)*u0v(I-1,j)))) + &
3209 (CS%DX_dyT(i,j)*((G%IdxCv(i,J)*v0v(i,J)) + (G%IdxCv(i,J-1)*v0v(i,J-1))))) * &
32107564 max(G%IdxCv(i,J)*G%IareaCv(i,J), G%IdxCv(i,J-1)*G%IareaCv(i,J-1)) ) )
32117564 CS%Ah_Max_xx(I,J) = 0.0
32127626 if (denom > 0.0) then
32135050 CS%Ah_Max_xx(I,J) = CS%bound_coef * 0.5 * Idt / denom
32145050 if (CS%EY24_EBT_BS) then
32150 CS%Ah_Max_xx_KS(i,j) = CS%bound_coef * 0.5 / (CS%KS_timescale * denom)
3216 endif
3217 endif
3218
3219 enddo ; enddo
32207443 do J=js-1,Jeq ; do I=is-1,Ieq
3221 denom = max( &
3222 (CS%dx2q(I,J) * &
3223 ((CS%DX_dyBu(I,J)*((u0u(I,j+1)*G%IdxCu(I,j+1)) + (u0u(I,j)*G%IdxCu(I,j)))) + &
3224 (CS%DY_dxBu(I,J)*((v0u(i+1,J)*G%IdyCv(i+1,J)) + (v0u(i,J)*G%IdyCv(i,J))))) * &
3225 max(G%IdxCu(I,j)*G%IareaCu(I,j), G%IdxCu(I,j+1)*G%IareaCu(I,j+1)) ), &
3226 (CS%dy2q(I,J) * &
3227 ((CS%DX_dyBu(I,J)*((u0v(I,j+1)*G%IdxCu(I,j+1)) + (u0v(I,j)*G%IdxCu(I,j)))) + &
3228 (CS%DY_dxBu(I,J)*((v0v(i+1,J)*G%IdyCv(i+1,J)) + (v0v(i,J)*G%IdyCv(i,J))))) * &
32297381 max(G%IdyCv(i,J)*G%IareaCv(i,J), G%IdyCv(i+1,J)*G%IareaCv(i+1,J)) ) )
32307381 CS%Ah_Max_xy(I,J) = 0.0
32317442 if (denom > 0.0) then
32325149 CS%Ah_Max_xy(I,J) = CS%bound_coef * 0.5 * Idt / denom
32335149 if (CS%EY24_EBT_BS) then
32340 CS%Ah_Max_xy_KS(i,j) = CS%bound_coef * 0.5 / (CS%KS_timescale * denom)
3235 endif
3236 endif
3237
3238 enddo ; enddo
32391 if (CS%debug) then
32400 call hchksum(CS%Ah_Max_xx, "Ah_Max_xx", G%HI, haloshift=0, unscale=US%L_to_m**4*US%s_to_T)
32410 call Bchksum(CS%Ah_Max_xy, "Ah_Max_xy", G%HI, haloshift=0, unscale=US%L_to_m**4*US%s_to_T)
3242 endif
3243 endif
3244 ! Register fields for output from this module.
3245 CS%id_normstress = register_diag_field('ocean_model', 'NoSt', diag%axesTL, Time, &
32461 'Normal Stress', 's-1', conversion=US%s_to_T)
3247 CS%id_shearstress = register_diag_field('ocean_model', 'ShSt', diag%axesBL, Time, &
32481 'Shear Stress', 's-1', conversion=US%s_to_T)
3249 CS%id_diffu = register_diag_field('ocean_model', 'diffu', diag%axesCuL, Time, &
32501 'Zonal Acceleration from Horizontal Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2)
3251 CS%id_diffv = register_diag_field('ocean_model', 'diffv', diag%axesCvL, Time, &
32521 'Meridional Acceleration from Horizontal Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2)
3253
3254 !CS%id_hf_diffu = register_diag_field('ocean_model', 'hf_diffu', diag%axesCuL, Time, &
3255 ! 'Fractional Thickness-weighted Zonal Acceleration from Horizontal Viscosity', &
3256 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
3257 !if ((CS%id_hf_diffu > 0) .and. (present(ADp))) then
3258 ! call safe_alloc_alloc(CS%hf_diffu,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)
3259 ! call safe_alloc_ptr(ADp%diag_hfrac_u,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)
3260 !endif
3261
3262 !CS%id_hf_diffv = register_diag_field('ocean_model', 'hf_diffv', diag%axesCvL, Time, &
3263 ! 'Fractional Thickness-weighted Meridional Acceleration from Horizontal Viscosity', &
3264 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
3265 !if ((CS%id_hf_diffv > 0) .and. (present(ADp))) then
3266 ! call safe_alloc_alloc(CS%hf_diffv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
3267 ! call safe_alloc_ptr(ADp%diag_hfrac_v,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
3268 !endif
3269
3270 CS%id_hf_diffu_2d = register_diag_field('ocean_model', 'hf_diffu_2d', diag%axesCu1, Time, &
3271 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Horizontal Viscosity', &
32721 'm s-2', conversion=US%L_T2_to_m_s2)
32731 if ((CS%id_hf_diffu_2d > 0) .and. (present(ADp))) then
32740 call safe_alloc_ptr(ADp%diag_hfrac_u,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)
3275 endif
3276
3277 CS%id_hf_diffv_2d = register_diag_field('ocean_model', 'hf_diffv_2d', diag%axesCv1, Time, &
3278 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Horizontal Viscosity', &
32791 'm s-2', conversion=US%L_T2_to_m_s2)
32801 if ((CS%id_hf_diffv_2d > 0) .and. (present(ADp))) then
32810 call safe_alloc_ptr(ADp%diag_hfrac_v,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
3282 endif
3283
3284 CS%id_h_diffu = register_diag_field('ocean_model', 'h_diffu', diag%axesCuL, Time, &
3285 'Thickness Multiplied Zonal Acceleration from Horizontal Viscosity', &
32861 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
32871 if ((CS%id_h_diffu > 0) .and. (present(ADp))) then
32880 call safe_alloc_ptr(ADp%diag_hu,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)
3289 endif
3290
3291 CS%id_h_diffv = register_diag_field('ocean_model', 'h_diffv', diag%axesCvL, Time, &
3292 'Thickness Multiplied Meridional Acceleration from Horizontal Viscosity', &
32931 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
32941 if ((CS%id_h_diffv > 0) .and. (present(ADp))) then
32950 call safe_alloc_ptr(ADp%diag_hv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
3296 endif
3297
3298 CS%id_intz_diffu_2d = register_diag_field('ocean_model', 'intz_diffu_2d', diag%axesCu1, Time, &
3299 'Depth-integral of Zonal Acceleration from Horizontal Viscosity', &
33001 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
33011 if ((CS%id_intz_diffu_2d > 0) .and. (present(ADp))) then
33020 call safe_alloc_ptr(ADp%diag_hu,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)
3303 endif
3304
3305 CS%id_intz_diffv_2d = register_diag_field('ocean_model', 'intz_diffv_2d', diag%axesCv1, Time, &
3306 'Depth-integral of Meridional Acceleration from Horizontal Viscosity', &
33071 'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
33081 if ((CS%id_intz_diffv_2d > 0) .and. (present(ADp))) then
33090 call safe_alloc_ptr(ADp%diag_hv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
3310 endif
3311
3312 CS%id_diffu_visc_rem = register_diag_field('ocean_model', 'diffu_visc_rem', diag%axesCuL, Time, &
3313 'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant', &
33141 'm s-2', conversion=US%L_T2_to_m_s2)
33151 if ((CS%id_diffu_visc_rem > 0) .and. (present(ADp))) then
33160 call safe_alloc_ptr(ADp%visc_rem_u,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)
3317 endif
3318
3319 CS%id_diffv_visc_rem = register_diag_field('ocean_model', 'diffv_visc_rem', diag%axesCvL, Time, &
3320 'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant', &
33211 'm s-2', conversion=US%L_T2_to_m_s2)
33221 if ((CS%id_diffv_visc_rem > 0) .and. (present(ADp))) then
33230 call safe_alloc_ptr(ADp%visc_rem_v,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)
3324 endif
3325
33261 if (CS%biharmonic) then
3327 CS%id_Ah_h = register_diag_field('ocean_model', 'Ahh', diag%axesTL, Time, &
3328 'Biharmonic Horizontal Viscosity at h Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T, &
3329 cmor_field_name='difmxybo', &
3330 cmor_long_name='Ocean lateral biharmonic viscosity', &
33311 cmor_standard_name='ocean_momentum_xy_biharmonic_diffusivity')
3332 CS%id_Ah_q = register_diag_field('ocean_model', 'Ahq', diag%axesBL, Time, &
33331 'Biharmonic Horizontal Viscosity at q Points', 'm4 s-1', conversion=US%L_to_m**4*US%s_to_T)
3334 CS%id_grid_Re_Ah = register_diag_field('ocean_model', 'grid_Re_Ah', diag%axesTL, Time, &
33351 'Grid Reynolds number for the Biharmonic horizontal viscosity at h points', 'nondim')
33361 if (CS%EY24_EBT_BS) then
3337 CS%id_visc_limit_h_flag = register_diag_field('ocean_model', 'visc_limit_h_flag', diag%axesTL, Time, &
33380 'Locations where the biharmonic viscosity reached the better_bound limiter at h points', 'nondim')
3339 CS%id_visc_limit_q_flag = register_diag_field('ocean_model', 'visc_limit_q_flag', diag%axesBL, Time, &
33400 'Locations where the biharmonic viscosity reached the better_bound limiter at q points', 'nondim')
3341 CS%id_visc_limit_h = register_diag_field('ocean_model', 'visc_limit_h', diag%axesTL, Time, &
33420 'Value of the biharmonic viscosity limiter at h points', 'nondim')
3343 CS%id_visc_limit_q = register_diag_field('ocean_model', 'visc_limit_q', diag%axesBL, Time, &
33440 'Value of the biharmonic viscosity limiter at q points', 'nondim')
3345 CS%id_visc_limit_h_frac = register_diag_field('ocean_model', 'visc_limit_h_frac', diag%axesTL, Time, &
33460 'Value of the biharmonic viscosity limiter at h points', 'nondim')
3347 CS%id_visc_limit_q_frac = register_diag_field('ocean_model', 'visc_limit_q_frac', diag%axesBL, Time, &
33480 'Value of the biharmonic viscosity limiter at q points', 'nondim')
3349 endif
3350
33511 if (CS%id_grid_Re_Ah > 0) &
3352 ! Compute the smallest biharmonic viscosity capable of modifying the
3353 ! velocity at floating point precision.
33540 CS%min_grid_Ah = spacing(1.) * min_grid_sp_h4 * Idt
3355 endif
33561 if (CS%Laplacian) then
3357 CS%id_Kh_h = register_diag_field('ocean_model', 'Khh', diag%axesTL, Time, &
3358 'Laplacian Horizontal Viscosity at h Points', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T, &
3359 cmor_field_name='difmxylo', &
3360 cmor_long_name='Ocean lateral Laplacian viscosity', &
33611 cmor_standard_name='ocean_momentum_xy_laplacian_diffusivity')
3362 CS%id_Kh_q = register_diag_field('ocean_model', 'Khq', diag%axesBL, Time, &
33631 'Laplacian Horizontal Viscosity at q Points', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T)
3364 CS%id_grid_Re_Kh = register_diag_field('ocean_model', 'grid_Re_Kh', diag%axesTL, Time, &
33651 'Grid Reynolds number for the Laplacian horizontal viscosity at h points', 'nondim')
3366 CS%id_vort_xy_q = register_diag_field('ocean_model', 'vort_xy_q', diag%axesBL, Time, &
33671 'Vertical vorticity at q Points', 's-1', conversion=US%s_to_T)
3368 CS%id_div_xx_h = register_diag_field('ocean_model', 'div_xx_h', diag%axesTL, Time, &
33691 'Horizontal divergence at h Points', 's-1', conversion=US%s_to_T)
3370 CS%id_sh_xy_q = register_diag_field('ocean_model', 'sh_xy_q', diag%axesBL, Time, &
33711 'Shearing strain at q Points', 's-1', conversion=US%s_to_T)
3372 CS%id_sh_xx_h = register_diag_field('ocean_model', 'sh_xx_h', diag%axesTL, Time, &
33731 'Horizontal tension at h Points', 's-1', conversion=US%s_to_T)
3374
33751 if (CS%id_grid_Re_Kh > 0) &
3376 ! Compute a smallest Laplacian viscosity capable of modifying the
3377 ! velocity at floating point precision.
33780 CS%min_grid_Kh = spacing(1.) * min_grid_sp_h2 * Idt
3379 endif
33801 if (CS%use_GME) then
3381 CS%id_dudx_bt = register_diag_field('ocean_model', 'dudx_bt', diag%axesT1, Time, &
3382 'Zonal component of the barotropic shearing strain at h points', 's-1', &
33830 conversion=US%s_to_T)
3384 CS%id_dudy_bt = register_diag_field('ocean_model', 'dudy_bt', diag%axesB1, Time, &
3385 'Zonal component of the barotropic shearing strain at q points', 's-1', &
33860 conversion=US%s_to_T)
3387 CS%id_dvdy_bt = register_diag_field('ocean_model', 'dvdy_bt', diag%axesT1, Time, &
3388 'Meridional component of the barotropic shearing strain at h points', 's-1', &
33890 conversion=US%s_to_T)
3390 CS%id_dvdx_bt = register_diag_field('ocean_model', 'dvdx_bt', diag%axesB1, Time, &
3391 'Meridional component of the barotropic shearing strain at q points', 's-1', &
33920 conversion=US%s_to_T)
3393 CS%id_GME_coeff_h = register_diag_field('ocean_model', 'GME_coeff_h', diag%axesTL, Time, &
33940 'GME coefficient at h Points', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T)
3395 CS%id_GME_coeff_q = register_diag_field('ocean_model', 'GME_coeff_q', diag%axesBL, Time, &
33960 'GME coefficient at q Points', 'm2 s-1', conversion=US%L_to_m**2*US%s_to_T)
3397 CS%id_FrictWork_GME = register_diag_field('ocean_model','FrictWork_GME',diag%axesTL,Time,&
3398 'Integral work done by lateral friction terms in GME (excluding diffusion of energy)', &
33990 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
3400 endif
3401
34021 if (CS%EY24_EBT_BS) then
3403 CS%id_BS_coeff_h = register_diag_field('ocean_model', 'BS_coeff_h', diag%axesTL, Time, &
34040 'Backscatter coefficient at h points', units='m2 s-1', conversion=US%L_to_m**2*US%s_to_T)
3405 CS%id_BS_coeff_q = register_diag_field('ocean_model', 'BS_coeff_q', diag%axesBL, Time, &
34060 'Backscatter coefficient at q points', units='m2 s-1', conversion=US%L_to_m**2*US%s_to_T)
3407 endif
3408
3409 CS%id_FrictWork = register_diag_field('ocean_model','FrictWork',diag%axesTL,Time,&
3410 'Integral work done by lateral friction terms. If GME is turned on, this '//&
3411 'includes the GME contribution.', &
34121 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
3413 CS%id_FrictWorkIntz = register_diag_field('ocean_model','FrictWorkIntz',diag%axesT1,Time, &
3414 'Depth integrated work done by lateral friction', &
3415 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2, &
3416 cmor_field_name='dispkexyfo', &
3417 cmor_long_name='Depth integrated ocean kinetic energy dissipation due to lateral friction',&
34181 cmor_standard_name='ocean_kinetic_energy_dissipation_per_unit_area_due_to_xy_friction')
3419 CS%id_FrictWork_bh = register_diag_field('ocean_model','FrictWork_bh',diag%axesTL,Time,&
3420 'Integral work done by the biharmonic lateral friction terms.', &
34211 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
3422 CS%id_FrictWorkIntz_bh = register_diag_field('ocean_model','FrictWorkIntz_bh',diag%axesT1,Time,&
3423 'Depth integrated work done by the biharmonic lateral friction', &
34241 'W m-2', conversion=US%RZ3_T3_to_W_m2*US%L_to_Z**2)
3425
3426
3427 ! TODO: Position these after their respective loops (and run loops on device)
3428 !$omp target enter data map(to: CS%Idxdy2u, CS%Idxdy2v) if (CS%biharmonic)
3429 !$omp target enter data map(to: CS%Idx2dyCu, CS%Idx2dyCv) if (CS%biharmonic)
3430
3431 !$omp target enter data map(to: CS%Kh_bg_xx, CS%Kh_bg_xy) if (CS%Laplacian)
3432 !$omp target enter data map(to: CS%Kh_max_xx) if (CS%Laplacian)
3433 !$omp target enter data map(to: CS%Kh_max_xy) &
3434 !$omp if (CS%Laplacian .and. CS%bound_Kh)
3435 !$omp target enter data map(to: CS%Laplac2_const_xx) if (CS%Laplacian)
3436 !$omp target enter data map(to: CS%Laplac3_const_xx) if (CS%Laplacian)
3437 !$omp target enter data map(to: CS%Laplac2_const_xy) if (CS%Smagorinsky_Kh)
3438
3439 !$omp target enter data map(to: CS%Ah_bg_xx, CS%Ah_bg_xy) if (CS%biharmonic)
3440 !$omp target enter data map(to: CS%reduction_xx, CS%reduction_xy)
3441 !$omp target enter data map(to: CS%Biharm_const_xx, CS%Biharm_const2_xx) &
3442 !$omp if (CS%Smagorinsky_Ah .or. CS%Leith_Ah .or. CS%use_Leithy)
3443 !$omp target enter data map(to: CS%Biharm_const_xy) &
3444 !$omp if (CS%Smagorinsky_Ah .or. CS%Leith_Ah)
3445 !$omp target enter data map(to: CS%Biharm_const2_xy) &
3446 !$omp if (CS%bound_Coriolis .and. (CS%Smagorinsky_Ah .or. CS%Leith_Ah))
3447 !$omp target enter data map(to: CS%Ah_max_xx) if (CS%bound_Ah)
3448 !$omp target enter data map(to: CS%Ah_max_xy) if (CS%bound_Ah)
3449
3450end subroutine hor_visc_init
3451
3452!> hor_visc_vel_stencil returns the horizontal viscosity input velocity stencil size
345324function hor_visc_vel_stencil(CS) result(stencil)
3454 type(hor_visc_CS), intent(in) :: CS !< Control structure for horizontal viscosity
3455 integer :: stencil !< The horizontal viscosity velocity stencil size with the current settings.
3456
345724 stencil = 2
3458
345924 if ((CS%Leith_Kh) .or. (CS%Leith_Ah) .or. (CS%use_Leithy)) then
34600 stencil = 3
3461 endif
346224end function hor_visc_vel_stencil
3463
3464!> Calculates factors in the anisotropic orientation tensor to be align with the grid.
3465!! With n1=1 and n2=0, this recovers the approach of Large et al, 2001.
34660subroutine align_aniso_tensor_to_grid(CS, n1, n2)
3467 type(hor_visc_CS), intent(inout) :: CS !< Control structure for horizontal viscosity
3468 real, intent(in) :: n1 !< i-component of direction vector [nondim]
3469 real, intent(in) :: n2 !< j-component of direction vector [nondim]
3470 ! Local variables
3471 real :: recip_n2_norm ! The inverse of the squared magnitude of n1 and n2 [nondim]
3472 ! For normalizing n=(n1,n2) in case arguments are not a unit vector
34730 recip_n2_norm = (n1**2) + (n2**2)
34740 if (recip_n2_norm > 0.) recip_n2_norm = 1. / recip_n2_norm
34750 CS%n1n2_h(:,:) = 2. * ( n1 * n2 ) * recip_n2_norm
34760 CS%n1n2_q(:,:) = 2. * ( n1 * n2 ) * recip_n2_norm
34770 CS%n1n1_m_n2n2_h(:,:) = ( (n1 * n1) - (n2 * n2) ) * recip_n2_norm
34780 CS%n1n1_m_n2n2_q(:,:) = ( (n1 * n1) - (n2 * n2) ) * recip_n2_norm
34790end subroutine align_aniso_tensor_to_grid
3480
3481!> Apply a 1-1-4-1-1 Laplacian filter one time on GME diffusive flux to reduce any
3482!! horizontal two-grid-point noise
34830subroutine smooth_GME(CS, G, GME_flux_h, GME_flux_q)
3484 type(hor_visc_CS), intent(in) :: CS !< Control structure
3485 type(ocean_grid_type), intent(in) :: G !< Ocean grid
3486 real, dimension(SZI_(G),SZJ_(G)), optional, intent(inout) :: GME_flux_h!< GME diffusive flux
3487 !! at h points [L2 T-2 ~> m2 s-2]
3488 real, dimension(SZIB_(G),SZJB_(G)), optional, intent(inout) :: GME_flux_q!< GME diffusive flux
3489 !! at q points [L2 T-2 ~> m2 s-2]
3490 ! local variables
34910 real, dimension(SZI_(G),SZJ_(G)) :: GME_flux_h_original ! The previous value of GME_flux_h [L2 T-2 ~> m2 s-2]
34920 real, dimension(SZIB_(G),SZJB_(G)) :: GME_flux_q_original ! The previous value of GME_flux_q [L2 T-2 ~> m2 s-2]
3493 real :: wc, ww, we, wn, ws ! averaging weights for smoothing [nondim]
3494 integer :: i, j, s, halosz
3495 integer :: xh, xq ! The number of valid extra halo points for h and q points.
3496 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq
3497
34980 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
34990 Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
35000 xh = 0 ; xq = 0
3501
35020 do s=1,CS%num_smooth_gme
35030 if (present(GME_flux_h)) then
35040 if (xh < 0) then
3505 ! Update halos if needed, but avoid doing so more often than is needed.
35060 halosz = min(G%isc-G%isd, G%jsc-G%jsd, 2+CS%num_smooth_gme-s)
35070 call pass_var(GME_flux_h, G%Domain, halo=halosz)
35080 xh = halosz - 2
3509 endif
35100 GME_flux_h_original(:,:) = GME_flux_h(:,:)
3511 ! apply smoothing on GME
35120 do j=Jsq-xh,Jeq+1+xh ; do i=Isq-xh,Ieq+1+xh
3513 ! skip land points
35140 if (G%mask2dT(i,j)==0.) cycle
3515 ! compute weights
35160 ww = 0.125 * G%mask2dT(i-1,j)
35170 we = 0.125 * G%mask2dT(i+1,j)
35180 ws = 0.125 * G%mask2dT(i,j-1)
35190 wn = 0.125 * G%mask2dT(i,j+1)
35200 wc = 1.0 - ((ww+we)+(wn+ws))
3521 GME_flux_h(i,j) = wc * GME_flux_h_original(i,j) &
3522 + ((ww * GME_flux_h_original(i-1,j) + we * GME_flux_h_original(i+1,j)) &
35230 + (ws * GME_flux_h_original(i,j-1) + wn * GME_flux_h_original(i,j+1)))
3524 enddo ; enddo
35250 xh = xh - 1
3526 endif
35270 if (present(GME_flux_q)) then
35280 if (xq < 0) then
3529 ! Update halos if needed, but avoid doing so more often than is needed.
35300 halosz = min(G%isc-G%isd, G%jsc-G%jsd, 2+CS%num_smooth_gme-s)
35310 call pass_var(GME_flux_q, G%Domain, position=CORNER, complete=.true., halo=halosz)
35320 xq = halosz - 2
3533 endif
35340 GME_flux_q_original(:,:) = GME_flux_q(:,:)
3535 ! apply smoothing on GME
35360 do J=js-1-xq,je+xq ; do I=is-1-xq,ie+xq
3537 ! skip land points
35380 if (G%mask2dBu(I,J)==0.) cycle
3539 ! compute weights
35400 ww = 0.125 * G%mask2dBu(I-1,J)
35410 we = 0.125 * G%mask2dBu(I+1,J)
35420 ws = 0.125 * G%mask2dBu(I,J-1)
35430 wn = 0.125 * G%mask2dBu(I,J+1)
35440 wc = 1.0 - ((ww+we)+(wn+ws))
3545 GME_flux_q(I,J) = wc * GME_flux_q_original(I,J) &
3546 + ((ww * GME_flux_q_original(I-1,J) + we * GME_flux_q_original(I+1,J)) &
35470 + (ws * GME_flux_q_original(I,J-1) + wn * GME_flux_q_original(I,J+1)))
3548 enddo ; enddo
35490 xq = xq - 1
3550 endif
3551 enddo ! s-loop
35520end subroutine smooth_GME
3553
3554!> Apply a 9-point smoothing filter twice to a field staggered at a thickness point to reduce
3555!! horizontal two-grid-point noise.
3556!! Note that this subroutine does not conserve mass, so don't use it in situations where you
3557!! need conservation. Also note that it assumes that the input field has valid values in the
3558!! first two halo points upon entry.
35590subroutine smooth_x9_h(G, field_h, zero_land)
3560 type(ocean_grid_type), intent(in) :: G !< Ocean grid
3561 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: field_h !< h-point field to be smoothed [arbitrary]
3562 logical, optional, intent(in) :: zero_land !< If present and false, return the average
3563 !! of the surrounding ocean points when
3564 !! smoothing, otherwise use a value of 0 for
3565 !! land points and include them in the averages.
3566
3567 ! Local variables
35680 real :: fh_prev(SZI_(G),SZJ_(G)) ! The value of the h-point field at the previous iteration [arbitrary]
3569 real :: Iwts ! The inverse of the sum of the weights [nondim]
3570 logical :: zero_land_val ! The value of the zero_land optional argument or .true. if it is absent.
3571 integer :: i, j, s, is, ie, js, je
3572
35730 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
3574
35750 zero_land_val = .true. ; if (present(zero_land)) zero_land_val = zero_land
3576
35770 do s=1,0,-1
35780 fh_prev(:,:) = field_h(:,:)
3579 ! apply smoothing on field_h using rotationally symmetric expressions.
35800 do j=js-s,je+s ; do i=is-s,ie+s ; if (G%mask2dT(i,j) > 0.0) then
35810 Iwts = 0.0625
35820 if (.not. zero_land_val) &
3583 Iwts = 1.0 / ( (4.0*G%mask2dT(i,j) + &
3584 ( 2.0*((G%mask2dT(i-1,j) + G%mask2dT(i+1,j)) + &
3585 (G%mask2dT(i,j-1) + G%mask2dT(i,j+1))) + &
3586 ((G%mask2dT(i-1,j-1) + G%mask2dT(i+1,j+1)) + &
35870 (G%mask2dT(i-1,j+1) + G%mask2dT(i+1,j-1))) ) ) + 1.0e-16 )
3588 field_h(i,j) = Iwts * ( 4.0*G%mask2dT(i,j) * fh_prev(i,j) &
3589 + (2.0*((G%mask2dT(i-1,j) * fh_prev(i-1,j) + G%mask2dT(i+1,j) * fh_prev(i+1,j)) + &
3590 (G%mask2dT(i,j-1) * fh_prev(i,j-1) + G%mask2dT(i,j+1) * fh_prev(i,j+1))) &
3591 + ((G%mask2dT(i-1,j-1) * fh_prev(i-1,j-1) + G%mask2dT(i+1,j+1) * fh_prev(i+1,j+1)) + &
35920 (G%mask2dT(i-1,j+1) * fh_prev(i-1,j+1) + G%mask2dT(i+1,j-1) * fh_prev(i-1,j-1))) ))
3593 endif ; enddo ; enddo
3594 enddo
3595
35960end subroutine smooth_x9_h
3597
3598!> Apply a 9-point smoothing filter twice to a pair of velocity components to reduce
3599!! horizontal two-grid-point noise.
3600!! Note that this subroutine does not conserve angular momentum, so don't use it
3601!! in situations where you need conservation. Also note that it assumes that the
3602!! input fields have valid values in the first two halo points upon entry.
36030subroutine smooth_x9_uv(G, field_u, field_v, zero_land)
3604 type(ocean_grid_type), intent(in) :: G !< Ocean grid
3605 real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: field_u !< u-point field to be smoothed [arbitrary]
3606 real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: field_v !< v-point field to be smoothed [arbitrary]
3607 logical, optional, intent(in) :: zero_land !< If present and false, return the average
3608 !! of the surrounding ocean points when
3609 !! smoothing, otherwise use a value of 0 for
3610 !! land points and include them in the averages.
3611
3612 ! Local variables.
36130 real :: fu_prev(SZIB_(G),SZJ_(G)) ! The value of the u-point field at the previous iteration [arbitrary]
36140 real :: fv_prev(SZI_(G),SZJB_(G)) ! The value of the v-point field at the previous iteration [arbitrary]
3615 real :: Iwts ! The inverse of the sum of the weights [nondim]
3616 logical :: zero_land_val ! The value of the zero_land optional argument or .true. if it is absent.
3617 integer :: i, j, s, is, ie, js, je, Isq, Ieq, Jsq, Jeq
3618
36190 is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
36200 Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
3621
36220 zero_land_val = .true. ; if (present(zero_land)) zero_land_val = zero_land
3623
36240 do s=1,0,-1
36250 fu_prev(:,:) = field_u(:,:)
3626 ! apply smoothing on field_u using the original non-rotationally symmetric expressions.
36270 do j=js-s,je+s ; do I=Isq-s,Ieq+s ; if (G%mask2dCu(I,j) > 0.0) then
36280 Iwts = 0.0625
36290 if (.not. zero_land_val) &
3630 Iwts = 1.0 / ( (4.0*G%mask2dCu(I,j) + &
3631 ( 2.0*((G%mask2dCu(I-1,j) + G%mask2dCu(I+1,j)) + &
3632 (G%mask2dCu(I,j-1) + G%mask2dCu(I,j+1))) + &
3633 ((G%mask2dCu(I-1,j-1) + G%mask2dCu(I+1,j+1)) + &
36340 (G%mask2dCu(I-1,j+1) + G%mask2dCu(I+1,j-1))) ) ) + 1.0e-16 )
3635 field_u(I,j) = Iwts * ( 4.0*G%mask2dCu(I,j) * fu_prev(I,j) &
3636 + (2.0*((G%mask2dCu(I-1,j) * fu_prev(I-1,j) + G%mask2dCu(I+1,j) * fu_prev(I+1,j)) + &
3637 (G%mask2dCu(I,j-1) * fu_prev(I,j-1) + G%mask2dCu(I,j+1) * fu_prev(I,j+1))) &
3638 + ((G%mask2dCu(I-1,j-1) * fu_prev(I-1,j-1) + G%mask2dCu(I+1,j+1) * fu_prev(I+1,j+1)) + &
36390 (G%mask2dCu(I-1,j+1) * fu_prev(I-1,j+1) + G%mask2dCu(I+1,j-1) * fu_prev(I-1,j-1))) ))
3640 endif ; enddo ; enddo
3641
36420 fv_prev(:,:) = field_v(:,:)
3643 ! apply smoothing on field_v using the original non-rotationally symmetric expressions.
36440 do J=Jsq-s,Jeq+s ; do i=is-s,ie+s ; if (G%mask2dCv(i,J) > 0.0) then
36450 Iwts = 0.0625
36460 if (.not. zero_land_val) &
3647 Iwts = 1.0 / ( (4.0*G%mask2dCv(i,J) + &
3648 ( 2.0*((G%mask2dCv(i-1,J) + G%mask2dCv(i+1,J)) + &
3649 (G%mask2dCv(i,J-1) + G%mask2dCv(i,J+1))) + &
3650 ((G%mask2dCv(i-1,J-1) + G%mask2dCv(i+1,J+1)) + &
36510 (G%mask2dCv(i-1,J+1) + G%mask2dCv(i+1,J-1))) ) ) + 1.0e-16 )
3652 field_v(i,J) = Iwts * ( 4.0*G%mask2dCv(i,J) * fv_prev(i,J) &
3653 + (2.0*((G%mask2dCv(i-1,J) * fv_prev(i-1,J) + G%mask2dCv(i+1,J) * fv_prev(i+1,J)) + &
3654 (G%mask2dCv(i,J-1) * fv_prev(i,J-1) + G%mask2dCv(i,J+1) * fv_prev(i,J+1))) &
3655 + ((G%mask2dCv(i-1,J-1) * fv_prev(i-1,J-1) + G%mask2dCv(i+1,J+1) * fv_prev(i+1,J+1)) + &
36560 (G%mask2dCv(i-1,J+1) * fv_prev(i-1,J+1) + G%mask2dCv(i+1,J-1) * fv_prev(i-1,J-1))) ))
3657 endif ; enddo ; enddo
3658 enddo
3659
36600end subroutine smooth_x9_uv
3661
3662!> Deallocates any variables allocated in hor_visc_init.
36631subroutine hor_visc_end(CS)
3664 type(hor_visc_CS), intent(inout) :: CS !< Horizontal viscosity control structure
3665
3666 !$omp target exit data map(delete: CS%DX_dyT, CS%DY_dxT)
3667 !$omp target exit data map(delete: CS%Dx_dyBu, CS%DY_dxBu)
3668
3669 !$omp target exit data map(delete: CS%Idxdy2u, CS%Idxdy2v) if (CS%biharmonic)
3670 !$omp target exit data map(delete: CS%Idx2dyCu, CS%Idx2dyCv) if (CS%biharmonic)
3671 !$omp target exit data map(delete: CS%dx2q, CS%dy2q)
3672 !$omp target exit data map(delete: CS%dx2h, CS%dy2h)
3673
3674 !$omp target exit data map(delete: CS%Kh_bg_xx, CS%Kh_bg_xy) if (CS%Laplacian)
3675 !$omp target exit data map(delete: CS%Kh_Max_xx) if (CS%Laplacian)
3676 !$omp target exit data map(delete: CS%Kh_max_xy) &
3677 !$omp if (CS%Laplacian .and. CS%bound_Kh)
3678 !$omp target exit data map(delete: CS%Laplac2_const_xx) if (CS%Laplacian)
3679 !$omp target exit data map(delete: CS%Laplac3_const_xx) if (CS%Laplacian)
3680 !$omp target exit data map(delete: CS%Laplac2_const_xy) if (CS%Smagorinsky_Kh)
3681
3682 !$omp target exit data map(delete: CS%Ah_bg_xx, CS%Ah_bg_xy) if (CS%biharmonic)
3683 !$omp target exit data map(delete: CS%reduction_xx, CS%reduction_xy)
3684 !$omp target exit data map(delete: CS%Biharm_const_xx, CS%Biharm_const2_xx) &
3685 !$omp if (CS%Smagorinsky_Ah .or. CS%Leith_Ah .or. CS%use_Leithy)
3686 !$omp target exit data map(delete: CS%Biharm_const_xy) &
3687 !$omp if (CS%Smagorinsky_Ah .or. CS%Leith_Ah)
3688 !$omp target exit data map(delete: CS%Biharm_const2_xy) &
3689 !$omp if (CS%bound_Coriolis .and. (CS%Smagorinsky_Ah .or. CS%Leith_Ah))
3690 !$omp target exit data map(delete: CS%Ah_max_xx) if (CS%bound_Ah)
3691 !$omp target exit data map(delete: CS%Ah_max_xy) if (CS%bound_Ah)
3692
36931 if (CS%Laplacian .or. CS%biharmonic) then
36941 DEALLOC_(CS%dx2h) ; DEALLOC_(CS%dx2q) ; DEALLOC_(CS%dy2h) ; DEALLOC_(CS%dy2q)
36951 DEALLOC_(CS%dx_dyT) ; DEALLOC_(CS%dy_dxT) ; DEALLOC_(CS%dx_dyBu) ; DEALLOC_(CS%dy_dxBu)
36961 DEALLOC_(CS%reduction_xx) ; DEALLOC_(CS%reduction_xy)
3697 endif
36981 if (CS%Laplacian) then
36991 DEALLOC_(CS%Kh_bg_xx) ; DEALLOC_(CS%Kh_bg_xy)
37001 DEALLOC_(CS%grid_sp_h2)
37011 if (allocated(CS%Kh_bg_2d)) deallocate(CS%Kh_bg_2d)
3702
37031 if (allocated(CS%Kh_Max_xx)) deallocate(CS%Kh_Max_xx)
37041 if (allocated(CS%Kh_Max_xy)) deallocate(CS%Kh_Max_xy)
37051 if (allocated(CS%Laplac2_const_xx)) deallocate(CS%Laplac2_const_xx)
37061 if (allocated(CS%Laplac2_const_xy)) deallocate(CS%Laplac2_const_xy)
37071 if (allocated(CS%Laplac3_const_xx)) deallocate(CS%Laplac3_const_xx)
37081 if (allocated(CS%Laplac3_const_xy)) deallocate(CS%Laplac3_const_xy)
3709 endif
37101 if (CS%biharmonic) then
37111 DEALLOC_(CS%grid_sp_h3)
37121 DEALLOC_(CS%Idx2dyCu) ; DEALLOC_(CS%Idx2dyCv)
37131 DEALLOC_(CS%Idxdy2u) ; DEALLOC_(CS%Idxdy2v)
37141 DEALLOC_(CS%Ah_bg_xx) ; DEALLOC_(CS%Ah_bg_xy)
3715
37161 if (allocated(CS%Ah_Max_xx)) deallocate(CS%Ah_Max_xx)
37171 if (allocated(CS%Ah_Max_xy)) deallocate(CS%Ah_Max_xy)
37181 if (allocated(CS%Ah_Max_xx_KS)) deallocate(CS%Ah_Max_xx_KS)
37191 if (allocated(CS%Ah_Max_xy_KS)) deallocate(CS%Ah_Max_xy_KS)
37201 if (allocated(CS%Biharm_const_xx)) deallocate(CS%Biharm_const_xx)
37211 if (allocated(CS%Biharm_const_xy)) deallocate(CS%Biharm_const_xy)
37221 if (allocated(CS%Biharm_const2_xx)) deallocate(CS%Biharm_const2_xx)
37231 if (allocated(CS%Biharm_const2_xy)) deallocate(CS%Biharm_const2_xy)
37241 if (allocated(CS%Biharm6_const_xx)) deallocate(CS%Biharm6_const_xx)
37251 if (allocated(CS%Biharm6_const_xy)) deallocate(CS%Biharm6_const_xy)
37261 if (allocated(CS%m_const_leithy)) deallocate(CS%m_const_leithy)
37271 if (allocated(CS%m_leithy_max)) deallocate(CS%m_leithy_max)
37281 if (allocated(CS%Re_Ah_const_xx)) deallocate(CS%Re_Ah_const_xx)
37291 if (allocated(CS%Re_Ah_const_xy)) deallocate(CS%Re_Ah_const_xy)
3730 endif
3731
37321 if (allocated(CS%n1n2_h)) deallocate(CS%n1n2_h)
37331 if (allocated(CS%n1n2_q)) deallocate(CS%n1n2_q)
37341 if (allocated(CS%n1n1_m_n2n2_h)) deallocate(CS%n1n1_m_n2n2_h)
37351 if (allocated(CS%n1n1_m_n2n2_q)) deallocate(CS%n1n1_m_n2n2_q)
3736
37371 if (CS%use_ZB2020) then
37380 call ZB2020_end(CS%ZB2020)
3739 endif
3740
37411end subroutine hor_visc_end
3742!> \namespace mom_hor_visc
3743!!
3744!! \section section_horizontal_viscosity Horizontal viscosity in MOM
3745!!
3746!! This module contains the subroutine horizontal_viscosity that calculates the
3747!! effects of horizontal viscosity, including parameterizations of the value of
3748!! the viscosity itself. Subroutine horizontal_viscosity calculates the acceleration due to
3749!! some combination of a biharmonic viscosity and a Laplacian viscosity. Either or
3750!! both may use a coefficient that depends on the shear and strain of the flow.
3751!! All metric terms are retained. The Laplacian is calculated as the divergence of
3752!! a stress tensor, using the form suggested by \cite Smagorinsky1993. The biharmonic
3753!! is calculated by twice applying the divergence of the stress tensor that is
3754!! used to calculate the Laplacian, but without the dependence on thickness in the
3755!! first pass. This form permits a variable viscosity, and indicates no
3756!! acceleration for either resting fluid or solid body rotation.
3757!!
3758!! The form of the viscous accelerations is discussed extensively in \cite griffies2000,
3759!! and the implementation here follows that discussion closely.
3760!! We use the notation of \cite Smith2003 with the exception that the
3761!! isotropic viscosity is \f$\kappa_h\f$.
3762!!
3763!! In general, the horizontal stress tensor can be written as
3764!! \f[
3765!! {\bf \sigma} =
3766!! \begin{pmatrix}
3767!! \frac{1}{2} \left( \sigma_D + \sigma_T \right) & \frac{1}{2} \sigma_S \\\\
3768!! \frac{1}{2} \sigma_S & \frac{1}{2} \left( \sigma_D - \sigma_T \right)
3769!! \end{pmatrix}
3770!! \f]
3771!! where \f$\sigma_D\f$, \f$\sigma_T\f$ and \f$\sigma_S\f$ are stresses associated with
3772!! invariant factors in the strain-rate tensor. For a Newtonian fluid, the stress
3773!! tensor is usually linearly related to the strain-rate tensor. The horizontal
3774!! strain-rate tensor is
3775!! \f[
3776!! \dot{\bf e} =
3777!! \begin{pmatrix}
3778!! \frac{1}{2} \left( \dot{e}_D + \dot{e}_T \right) & \frac{1}{2} \dot{e}_S \\\\
3779!! \frac{1}{2} \dot{e}_S & \frac{1}{2} \left( \dot{e}_D - \dot{e}_T \right)
3780!! \end{pmatrix}
3781!! \f]
3782!! where \f$\dot{e}_D = \partial_x u + \partial_y v\f$ is the horizontal divergence,
3783!! \f$\dot{e}_T = \partial_x u - \partial_y v\f$ is the horizontal tension, and
3784!! \f$\dot{e}_S = \partial_y u + \partial_x v\f$ is the horizontal shear strain.
3785!!
3786!! The trace of the stress tensor, \f$tr(\bf \sigma) = \sigma_D\f$, is usually
3787!! absorbed into the pressure and only the deviatoric stress tensor considered.
3788!! From here on, we drop \f$\sigma_D\f$. The trace of the strain tensor, \f$tr(\bf e) =
3789!! \dot{e}_D\f$ is non-zero for horizontally divergent flow but only enters the
3790!! stress tensor through \f$\sigma_D\f$ and so we will drop \f$\sigma_D\f$ from
3791!! calculations of the strain tensor in the code. Therefore the horizontal stress
3792!! tensor can be considered to be
3793!! \f[
3794!! {\bf \sigma} =
3795!! \begin{pmatrix}
3796!! \frac{1}{2} \sigma_T & \frac{1}{2} \sigma_S \\\\
3797!! \frac{1}{2} \sigma_S & - \frac{1}{2} \sigma_T
3798!! \end{pmatrix}
3799!! .\f]
3800!!
3801!! The stresses above are linearly related to the strain through a viscosity
3802!! coefficient, \f$\kappa_h\f$:
3803!! \f{eqnarray*}{
3804!! \sigma_T & = & 2 \kappa_h \dot{e}_T \\\\
3805!! \sigma_S & = & 2 \kappa_h \dot{e}_S
3806!! .
3807!! \f}
3808!!
3809!! The viscosity \f$\kappa_h\f$ may either be a constant or variable. For example,
3810!! \f$\kappa_h\f$ may vary with the shear, as proposed by \cite Smagorinsky1993.
3811!!
3812!! The accelerations resulting form the divergence of the stress tensor are
3813!! \f{eqnarray*}{
3814!! \hat{\bf x} \cdot \left( \nabla \cdot {\bf \sigma} \right)
3815!! & = &
3816!! \partial_x \left( \frac{1}{2} \sigma_T \right)
3817!! + \partial_y \left( \frac{1}{2} \sigma_S \right)
3818!! \\\\
3819!! & = &
3820!! \partial_x \left( \kappa_h \dot{e}_T \right)
3821!! + \partial_y \left( \kappa_h \dot{e}_S \right)
3822!! \\\\
3823!! \hat{\bf y} \cdot \left( \nabla \cdot {\bf \sigma} \right)
3824!! & = &
3825!! \partial_x \left( \frac{1}{2} \sigma_S \right)
3826!! + \partial_y \left( - \frac{1}{2} \sigma_T \right)
3827!! \\\\
3828!! & = &
3829!! \partial_x \left( \kappa_h \dot{e}_S \right)
3830!! + \partial_y \left( - \kappa_h \dot{e}_T \right)
3831!! .
3832!! \f}
3833!!
3834!! The form of the Laplacian viscosity in general coordinates is:
3835!! \f{eqnarray*}{
3836!! \hat{\bf x} \cdot \left( \nabla \cdot \sigma \right)
3837!! & = &
3838!! \frac{1}{h} \left[ \partial_x \left( \kappa_h h \dot{e}_T \right)
3839!! + \partial_y \left( \kappa_h h \dot{e}_S \right) \right]
3840!! \\\\
3841!! \hat{\bf y} \cdot \left( \nabla \cdot \sigma \right)
3842!! & = &
3843!! \frac{1}{h} \left[ \partial_x \left( \kappa_h h \dot{e}_S \right)
3844!! - \partial_y \left( \kappa_h h \dot{e}_T \right) \right]
3845!! .
3846!! \f}
3847!!
3848!! \subsection section_laplacian_viscosity_coefficient Laplacian viscosity coefficient
3849!!
3850!! The horizontal viscosity coefficient, \f$\kappa_h\f$, can have multiple components.
3851!! The isotropic components are:
3852!! - A uniform background component, \f$\kappa_{bg}\f$.
3853!! - A constant but spatially variable 2D map, \f$\kappa_{2d}(x,y)\f$.
3854!! - A ''MICOM'' viscosity, \f$U_\nu \Delta(x,y)\f$, which uses a constant
3855!! velocity scale, \f$U_\nu\f$ and a measure of the grid-spacing \f$\Delta(x,y)^2 =
3856!! \frac{2 \Delta x^2 \Delta y^2}{\Delta x^2 + \Delta y^2}\f$.
3857!! - A function of
3858!! latitude, \f$\kappa_{\phi}(x,y) = \kappa_{\pi/2} |\sin(\phi)|^n\f$.
3859!! - A dynamic Smagorinsky viscosity, \f$\kappa_{Sm}(x,y,t) = C_{Sm} \Delta^2 \sqrt{\dot{e}_T^2 + \dot{e}_S^2}\f$.
3860!! - A dynamic Leith viscosity, \f$\kappa_{Lth}(x,y,t) =
3861!! C_{Lth} \Delta^3 \sqrt{|\nabla \zeta|^2 + |\nabla \dot{e}_D|^2}\f$.
3862!!
3863!! A maximum stable viscosity, \f$\kappa_{max}(x,y)\f$ is calculated based on the
3864!! grid-spacing and time-step and used to clip calculated viscosities.
3865!!
3866!! The static components of \f$\kappa_h\f$ are first combined as follows:
3867!! \f[
3868!! \kappa_{static} = \min \left[ \max\left(
3869!! \kappa_{bg},
3870!! U_\nu \Delta(x,y),
3871!! \kappa_{2d}(x,y),
3872!! \kappa_\phi(x,y)
3873!! \right)
3874!! , \kappa_{max}(x,y) \right]
3875!! \f]
3876!! and stored in the module control structure as variables <code>Kh_bg_xx</code> and
3877!! <code>Kh_bg_xy</code> for the tension (h-points) and shear (q-points) components
3878!! respectively.
3879!!
3880!! The full viscosity includes the dynamic components as follows:
3881!! \f[
3882!! \kappa_h(x,y,t) = r(\Delta,L_d)
3883!! \max \left( \kappa_{static}, \kappa_{Sm}, \kappa_{Lth} \right)
3884!! \f]
3885!! where \f$r(\Delta,L_d)\f$ is a resolution function.
3886!!
3887!! The dynamic Smagorinsky and Leith viscosity schemes are exclusive with each
3888!! other.
3889!!
3890!! \subsection section_viscous_boundary_conditions Viscous boundary conditions
3891!!
3892!! Free slip boundary conditions have been coded, although no slip boundary
3893!! conditions can be used with the Laplacian viscosity based on the 2D land-sea
3894!! mask. For a western boundary, for example, the boundary conditions with the
3895!! biharmonic operator would be written as:
3896!! \f[
3897!! \partial_x v = 0 ; \partial_x^3 v = 0 ; u = 0 ; \partial_x^2 u = 0 ,
3898!! \f]
3899!! while for a Laplacian operator, they are simply
3900!! \f[
3901!! \partial_x v = 0 ; u = 0 .
3902!! \f]
3903!! These boundary conditions are largely dictated by the use of an Arakawa
3904!! C-grid and by the varying layer thickness.
3905!!
3906!! \subsection section_anisotropic_viscosity Anisotropic viscosity
3907!!
3908!! \cite Large2001 proposed enhancing viscosity in a particular direction and the
3909!! approach was generalized in \cite Smith2003. We use the second form of their
3910!! two coefficient anisotropic viscosity (section 4.3). We also replace their
3911!! \f$A^\prime\f$ and $D$ such that \f$2A^\prime = 2 \kappa_h + D\f$ and
3912!! \f$\kappa_a = D\f$ so that \f$\kappa_h\f$ can be considered the isotropic
3913!! viscosity and \f$\kappa_a=D\f$ can be consider the anisotropic viscosity. The
3914!! direction of anisotropy is defined by a unit vector \f$\hat{\bf
3915!! n}=(n_1,n_2)\f$.
3916!!
3917!! The contributions to the stress tensor are
3918!! \f[
3919!! \begin{pmatrix}
3920!! \sigma_T \\\\ \sigma_S
3921!! \end{pmatrix}
3922!! =
3923!! \left[
3924!! \begin{pmatrix}
3925!! 2 \kappa_h + \kappa_a & 0 \\\\
3926!! 0 & 2 \kappa_h
3927!! \end{pmatrix}
3928!! + 2 \kappa_a n_1 n_2
3929!! \begin{pmatrix}
3930!! - 2 n_1 n_2 & n_1^2 - n_2^2 \\\\
3931!! n_1^2 - n_2^2 & 2 n_1 n_2
3932!! \end{pmatrix}
3933!! \right]
3934!! \begin{pmatrix}
3935!! \dot{e}_T \\\\ \dot{e}_S
3936!! \end{pmatrix}
3937!! \f]
3938!! Dissipation of kinetic energy requires \f$\kappa_h \geq 0\f$ and \f$2 \kappa_h + \kappa_a \geq 0\f$.
3939!! Note that when anisotropy is aligned with the x-direction, \f$n_1 = \pm 1\f$, then
3940!! \f$n_2 = 0\f$ and the cross terms vanish. The accelerations in this aligned limit
3941!! with constant coefficients become
3942!! \f{eqnarray*}{
3943!! \hat{\bf x} \cdot \nabla \cdot {\bf \sigma}
3944!! & = &
3945!! \partial_x \left( \left( \kappa_h + \frac{1}{2} \kappa_a \right) \dot{e}_T \right)
3946!! + \partial_y \left( \kappa_h \dot{e}_S \right)
3947!! \\\\
3948!! & = &
3949!! \left( \kappa_h + \kappa_a \right) \partial_{xx} u
3950!! + \kappa_h \partial_{yy} u
3951!! - \frac{1}{2} \kappa_a \partial_x \left( \partial_x u + \partial_y v \right)
3952!! \\\\
3953!! \hat{\bf y} \cdot \nabla \cdot {\bf \sigma}
3954!! & = &
3955!! \partial_x \left( \kappa_h \dot{e}_S \right)
3956!! - \partial_y \left( \left( \kappa_h + \frac{1}{2} \kappa_a \right) \dot{e}_T \right)
3957!! \\\\
3958!! & = &
3959!! \kappa_h \partial_{xx} v
3960!! + \left( \kappa_h + \kappa_a \right) \partial_{yy} v
3961!! - \frac{1}{2} \kappa_a \partial_y \left( \partial_x u + \partial_y v \right)
3962!! \f}
3963!! which has contributions akin to a negative divergence damping (a divergence
3964!! enhancement?) but which is weaker than the enhanced tension terms by half.
3965!!
3966!! \subsection section_viscous_discretization Discretization
3967!!
3968!! The horizontal tension, \f$\dot{e}_T\f$, is stored in variable <code>sh_xx</code> and
3969!! discretized as
3970!! \f[
3971!! \dot{e}_T
3972!! = \frac{\Delta y}{\Delta x} \delta_i \left( \frac{1}{\Delta y} u \right)
3973!! - \frac{\Delta x}{\Delta y} \delta_j \left( \frac{1}{\Delta x} v \right)
3974!! .
3975!! \f]
3976!! The horizontal divergent strain, \f$\dot{e}_D\f$, is stored in variable
3977!! <code>div_xx</code> and discretized as
3978!! \f[
3979!! \dot{e}_D
3980!! = \frac{1}{h A} \left( \delta_i \left( \overline{h}^i \Delta y \, u \right)
3981!! + \delta_j \left( \overline{h}^j\Delta x \, v \right) \right)
3982!! .
3983!! \f]
3984!! Note that for expediency this is the exact discretization used in the
3985!! continuity equation.
3986!!
3987!! The horizontal shear strain, \f$\dot{e}_S\f$, is stored in variable <code>sh_xy</code>
3988!! and discretized as
3989!! \f[
3990!! \dot{e}_S = v_x + u_y
3991!! \f]
3992!! where
3993!! \f{align*}{
3994!! v_x &= \frac{\Delta y}{\Delta x} \delta_i \left( \frac{1}{\Delta y} v \right) \\\\
3995!! u_y &= \frac{\Delta x}{\Delta y} \delta_j \left( \frac{1}{\Delta x} u \right)
3996!! \f}
3997!! which are calculated separately so that no-slip or free-slip boundary
3998!! conditions can be applied to \f$v_x\f$ and \f$u_y\f$ where appropriate.
3999!!
4000!! The tendency for the x-component of the divergence of stress is stored in
4001!! variable <code>diffu</code> and discretized as
4002!! \f[
4003!! \hat{\bf x} \cdot \left( \nabla \cdot {\bf \sigma} \right) =
4004!! \frac{1}{A \overline{h}^i} \left(
4005!! \frac{1}{\Delta y} \delta_i \left( h \Delta y^2 \kappa_h \dot{e}_T \right) +
4006!! \frac{1}{\Delta x} \delta_j \left( \tilde{h}^{ij} \Delta x^2 \kappa_h \dot{e}_S \right)
4007!! \right)
4008!! .
4009!! \f]
4010!!
4011!! The tendency for the y-component of the divergence of stress is stored in
4012!! variable <code>diffv</code> and discretized as
4013!! \f[
4014!! \hat{\bf y} \cdot \left( \nabla \cdot {\bf \sigma} \right) =
4015!! \frac{1}{A \overline{h}^j} \left(
4016!! \frac{1}{\Delta y} \delta_i \left( \tilde{h}^{ij} \Delta y^2 A_M \dot{e}_S \right)
4017!! - \frac{1}{\Delta x} \delta_j \left( h \Delta x^2 A_M \dot{e}_T \right)
4018!! \right)
4019!! .
4020!! \f]
4021!!
4022!! \subsection section_viscous_refs References
4023!!
4024!! Griffies, S.M., and Hallberg, R.W., 2000: Biharmonic friction with a
4025!! Smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models.
4026!! Monthly Weather Review, 128(8), 2935-2946.
4027!! https://doi.org/10.1175/1520-0493(2000)128%3C2935:BFWASL%3E2.0.CO;2
4028!!
4029!! Large, W.G., Danabasoglu, G., McWilliams, J.C., Gent, P.R. and Bryan, F.O.,
4030!! 2001: Equatorial circulation of a global ocean climate model with
4031!! anisotropic horizontal viscosity.
4032!! Journal of Physical Oceanography, 31(2), pp.518-536.
4033!! https://doi.org/10.1175/1520-0485(2001)031%3C0518:ECOAGO%3E2.0.CO;2
4034!!
4035!! Smagorinsky, J., 1993: Some historical remarks on the use of nonlinear
4036!! viscosities. Large eddy simulation of complex engineering and geophysical
4037!! flows, 1, 69-106.
4038!!
4039!! Smith, R.D., and McWilliams, J.C., 2003: Anisotropic horizontal viscosity for
4040!! ocean models. Ocean Modelling, 5(2), 129-156.
4041!! https://doi.org/10.1016/S1463-5003(02)00016-1
40420end module MOM_hor_visc