MOM_hor_visc.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.

2! See the LICENSE file for licensing information.

3! SPDX-License-Identifier: Apache-2.0

4

5!> 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.

270subroutine horizontal_viscosity(u, v, h, uh, vh, diffu, diffv, MEKE, VarMix, G, GV, US, &

271 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)) :: &

310 del2u, & ! The u-component of the Laplacian of velocity [L-1 T-1 ~> m-1 s-1]

311 h_u, & ! Thickness interpolated to u points [H ~> m or kg m-2].

312 vort_xy_dy, & ! y-derivative of vertical vorticity (d/dy(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]

313 vort_xy_dy_smooth, & ! y-derivative of smoothed vertical vorticity [L-1 T-1 ~> m-1 s-1]

314 div_xx_dx, & ! x-derivative of horizontal divergence (d/dx(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]

315 ubtav ! zonal barotropic velocity averaged over a baroclinic time-step [L T-1 ~> m s-1]

316 real, dimension(SZI_(G),SZJB_(G)) :: &

317 del2v, & ! The v-component of the Laplacian of velocity [L-1 T-1 ~> m-1 s-1]

318 h_v, & ! Thickness interpolated to v points [H ~> m or kg m-2].

319 vort_xy_dx, & ! x-derivative of vertical vorticity (d/dx(dv/dx - du/dy)) [L-1 T-1 ~> m-1 s-1]

320 vort_xy_dx_smooth, & ! x-derivative of smoothed vertical vorticity [L-1 T-1 ~> m-1 s-1]

321 div_xx_dy, & ! y-derivative of horizontal divergence (d/dy(du/dx + dv/dy)) [L-1 T-1 ~> m-1 s-1]

322 vbtav ! meridional barotropic velocity averaged over a baroclinic time-step [L T-1 ~> m s-1]

323 real, dimension(SZI_(G),SZJ_(G)) :: &

324 dudx_bt, dvdy_bt, & ! components in the barotropic horizontal tension [T-1 ~> s-1]

325 div_xx, & ! Estimate of horizontal divergence at h-points [T-1 ~> s-1]

326 sh_xx, & ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1]

327 sh_xx_smooth, & ! horizontal tension from smoothed velocity including metric terms [T-1 ~> s-1]

328 sh_xx_bt, & ! barotropic horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1]

329 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].

331 str_xx_gme,& ! smoothed diagonal term in the stress tensor from GME [L2 T-2 ~> m2 s-2]

332 bhstr_xx, & ! A copy of str_xx that only contains the biharmonic contribution [H L2 T-2 ~> m3 s-2 or kg s-2]

333 frictworkintz, & ! depth integrated energy dissipated by lateral friction [R Z L2 T-3 ~> W m-2]

334 frictworkintz_bh, & ! depth integrated energy dissipated by biharmonic lateral friction [R Z L2 T-3 ~> W m-2]

335 grad_vort_mag_h, & ! Magnitude of vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1]

336 grad_vort_mag_h_2d, & ! Magnitude of 2d vorticity gradient at h-points [L-1 T-1 ~> m-1 s-1]

337 grad_div_mag_h, & ! Magnitude of divergence gradient at h-points [L-1 T-1 ~> m-1 s-1]

338 dudx, dvdy, & ! components in the horizontal tension [T-1 ~> s-1]

339 dudx_smooth, dvdy_smooth, & ! components in the horizontal tension from smoothed velocity [T-1 ~> s-1]

340 gme_effic_h, & ! The filtered efficiency of the GME terms at h points [nondim]

341 m_leithy, & ! Kh=m_leithy*Ah in Leith+E parameterization [L-2 ~> m-2]

342 ah_sq, & ! The square of the biharmonic viscosity [L8 T-2 ~> m8 s-2]

343 htot, & ! The total thickness of all layers [H ~> m or kg m-2]

344 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)) :: &

350 dvdx, dudy, & ! components in the shearing strain [T-1 ~> s-1]

351 dvdx_smooth, dudy_smooth, & ! components in the shearing strain from smoothed velocity [T-1 ~> s-1]

352 ddel2vdx, ddel2udy, & ! Components in the biharmonic equivalent of the shearing strain [L-2 T-1 ~> m-2 s-1]

353 dvdx_bt, dudy_bt, & ! components in the barotropic shearing strain [T-1 ~> s-1]

354 sh_xy, & ! horizontal shearing strain (du/dy + dv/dx) including metric terms [T-1 ~> s-1]

355 sh_xy_smooth, & ! horizontal shearing strain from smoothed velocity including metric terms [T-1 ~> s-1]

356 sh_xy_bt, & ! barotropic horizontal shearing strain (du/dy + dv/dx) inc. metric terms [T-1 ~> s-1]

357 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].

359 str_xy_gme, & ! smoothed cross term in the stress tensor from GME [L2 T-2 ~> m2 s-2]

360 bhstr_xy, & ! A copy of str_xy that only contains the biharmonic contribution [H L2 T-2 ~> m3 s-2 or kg s-2]

361 vort_xy, & ! Vertical vorticity (dv/dx - du/dy) including metric terms [T-1 ~> s-1]

362 vort_xy_smooth, & ! Vertical vorticity including metric terms, smoothed [T-1 ~> s-1]

363 grad_vort_mag_q, & ! Magnitude of vorticity gradient at q-points [L-1 T-1 ~> m-1 s-1]

364 grad_vort_mag_q_2d, & ! Magnitude of 2d vorticity gradient at q-points [L-1 T-1 ~> m-1 s-1]

365 del2vort_q, & ! Laplacian of vorticity at q-points [L-2 T-1 ~> m-2 s-1]

366 grad_div_mag_q, & ! Magnitude of divergence gradient at q-points [L-1 T-1 ~> m-1 s-1]

367 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.

369 gme_effic_q, & ! The filtered efficiency of the GME terms at q points [nondim]

370 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)) :: &

375 ah_q, & ! biharmonic viscosity at corner points [L4 T-1 ~> m4 s-1]

376 kh_q, & ! Laplacian viscosity at corner points [L2 T-1 ~> m2 s-1]

377 vort_xy_q, & ! vertical vorticity at corner points [T-1 ~> s-1]

378 sh_xy_q, & ! horizontal shearing strain at corner points [T-1 ~> s-1]

379 gme_coeff_q, & !< GME coeff. at q-points [L2 T-1 ~> m2 s-1]

380 visc_limit_q, & ! used to stabilize the EY24_EBT_BS backscatter [nondim]

381 visc_limit_q_flag, & ! determines whether backscatter is shut off [nondim]

382 visc_limit_q_frac, & ! determines how close backscatter is to shutting off [nondim]

383 bs_coeff_q, & ! A diagnostic array of the backscatter coefficient [L2 T-1 ~> m2 s-1]

384 shst ! A diagnostic array of shear stress [T-1 ~> s-1].

385 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1) :: &

386 kh_u_gme, & !< Isopycnal height diffusivities in u-columns [L2 T-1 ~> m2 s-1]

387 slope_x !< Isopycnal slope in i-direction [Z L-1 ~> nondim]

388 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1) :: &

389 kh_v_gme, & !< Isopycnal height diffusivities in v-columns [L2 T-1 ~> m2 s-1]

390 slope_y !< Isopycnal slope in j-direction [Z L-1 ~> nondim]

391 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: &

392 ah_h, & ! biharmonic viscosity at thickness points [L4 T-1 ~> m4 s-1]

393 kh_h, & ! Laplacian viscosity at thickness points [L2 T-1 ~> m2 s-1]

394 dz, & ! Height change across layers [Z ~> m]

395 frictwork, & ! work done by MKE dissipation mechanisms [R Z L2 T-3 ~> W m-2]

396 frictwork_bh, & ! work done by the biharmonic MKE dissipation mechanisms [R Z L2 T-3 ~> W m-2]

397 frictwork_gme, & ! work done by GME [R Z L2 T-3 ~> W m-2]

398 div_xx_h, & ! horizontal divergence [T-1 ~> s-1]

399 sh_xx_h, & ! horizontal tension (du/dx - dv/dy) including metric terms [T-1 ~> s-1]

400 nost, & ! A diagnostic array of normal stress [T-1 ~> s-1].

401 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)) :: &

403 grid_re_kh, & ! Grid Reynolds number for Laplacian horizontal viscosity at h points [nondim]

404 grid_re_ah, & ! Grid Reynolds number for Biharmonic horizontal viscosity at h points [nondim]

405 gme_coeff_h, & ! GME coefficient at h-points [L2 T-1 ~> m2 s-1]

406 visc_limit_h, & ! Used to stabilize the EY24_EBT_BS backscatter [nondim]

407 visc_limit_h_flag, & ! determines whether backscatter is shut off [nondim]

408 visc_limit_h_frac ! determines how close backscatter is to shutting off [nondim]

409 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: &

410 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)) :: &

412 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 :: apply_obc = .false.

451 logical :: apply_obc_strain

452 logical :: use_meke_ku

453 logical :: use_meke_au

454 logical :: skeb_use_frict

455 logical :: use_cont_huv

456 logical :: use_kh_struct

457 integer :: is_vort, ie_vort, js_vort, je_vort ! Loop ranges for vorticity terms

458 integer :: is_kh, ie_kh, js_kh, je_kh ! Loop ranges for thickness point viscosities

459 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz

460 integer :: i, j, k, n

461 real :: inv_pi3, inv_pi2, inv_pi6 ! Powers of the inverse of pi [nondim]

462 real :: tmp

463

464 ! Fields evaluated on active layers, used for constructing 3D stress fields

465 ! NOTE: The position of these declarations can impact performance, due to the

466 ! very large number of stack arrays in this function. Move with caution!

467 ! NOTE: Several of these are declared with the memory extent of q-points, but the

468 ! same arrays are also used at h-points to reduce the memory footprint of this

469 ! module, so they should never be used in halo point or checksum calls.

470 real, dimension(SZIB_(G),SZJB_(G)) :: &

471 ah, & ! biharmonic viscosity (h or q) [L4 T-1 ~> m4 s-1]

472 kh, & ! Laplacian viscosity (h or q) [L2 T-1 ~> m2 s-1]

473 kh_bs, & ! Laplacian antiviscosity [L2 T-1 ~> m2 s-1]

474 shear_mag, & ! magnitude of the shear (h or q) [T-1 ~> s-1]

475 vert_vort_mag, & ! magnitude of the vertical vorticity gradient (h or q) [L-1 T-1 ~> m-1 s-1]

476 vert_vort_mag_smooth, & ! magnitude of gradient of smoothed vertical vorticity (h or q) [L-1 T-1 ~> m-1 s-1]

477 hrat_min, & ! h_min divided by the thickness at the stress point (h or q) [nondim]

478 visc_bound_rem ! fraction of overall viscous bounds that remain to be applied (h or q) [nondim]

479

480 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

481 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

482

483 h_neglect = gv%H_subroundoff

484 !h_neglect3 = h_neglect**3

485 h_neglect3 = h_neglect*h_neglect*h_neglect

486 inv_pi3 = 1.0/((4.0*atan(1.0))**3)

487 inv_pi2 = 1.0/((4.0*atan(1.0))**2)

488 inv_pi6 = inv_pi3 * inv_pi3

489

490 if (cs%EY24_EBT_BS) then

491 visc_limit_h(:,:,:) = 0.

492 visc_limit_q(:,:,:) = 0.

493 visc_limit_h_flag(:,:,:) = 0.

494 visc_limit_q_flag(:,:,:) = 0.

495 visc_limit_h_frac(:,:,:) = 0.

496 visc_limit_q_frac(:,:,:) = 0.

497 endif

498

499 skeb_use_frict = .false.

500 if (present(stoch)) skeb_use_frict = stoch%skeb_use_frict

501

502 m_leithy(:,:) = 0.0 ! Initialize

503

504 if (present(obc)) then ; if (associated(obc)) then ; if (obc%OBC_pe) then

505 apply_obc = obc%Flather_u_BCs_exist_globally .or. obc%Flather_v_BCs_exist_globally

506 apply_obc = .true.

507 endif ; endif ; endif

508

509 apply_obc_strain = .false.

510 if (present(obc)) then ; if (associated(obc)) then

511 apply_obc_strain = (obc%strain_config /= obc_strain_none) &

512 .and. ((.not. cs%OBC_strain_bug) .or. (obc%strain_config /= obc_strain_specified))

513 endif ; endif

514

515 if (.not.cs%initialized) call mom_error(fatal, &

516 "MOM_hor_visc: Module must be initialized before it is used.")

517

518 if (.not.(cs%Laplacian .or. cs%biharmonic)) return

519

520 find_frictwork = (cs%id_FrictWork > 0)

521 if (cs%id_FrictWorkIntz > 0) find_frictwork = .true.

522

523 if (allocated(meke%mom_src)) find_frictwork = .true.

524 use_kh_struct = allocated(varmix%BS_struct)

525 backscat_subround = 0.0

526 if (find_frictwork .and. allocated(meke%mom_src) .and. (meke%backscatter_Ro_c > 0.0) .and. &

527 (meke%backscatter_Ro_Pow /= 0.0)) &

528 backscat_subround = (1.0e-16/meke%backscatter_Ro_c)**(1.0/meke%backscatter_Ro_Pow)

529

530 ! Toggle whether to use a Laplacian viscosity derived from MEKE

531 use_meke_ku = allocated(meke%Ku)

532 use_meke_au = allocated(meke%Au)

533

534 use_cont_huv = cs%use_cont_thick .and. present(hu_cont) .and. present(hv_cont)

535 if (use_cont_huv .and. .not.cs%use_cont_thick_bug) then

536 call pass_vector(hu_cont, hv_cont, g%domain, to_all+scalar_pair, halo=2)

537 endif

538

539 rescale_kh = .false.

540 if (varmix%use_variable_mixing) then

541 rescale_kh = varmix%Resoln_scaled_Kh

542 if ((rescale_kh .or. cs%res_scale_MEKE) &

543 .and. (.not. allocated(varmix%Res_fn_h) .or. .not. allocated(varmix%Res_fn_q))) &

544 call mom_error(fatal, "MOM_hor_visc: VarMix%Res_fn_h and VarMix%Res_fn_q "//&

545 "both need to be associated with Resoln_scaled_Kh or RES_SCALE_MEKE_VISC.")

546 elseif (cs%res_scale_MEKE) then

547 call mom_error(fatal, "MOM_hor_visc: VarMix needs to be associated if "//&

548 "RES_SCALE_MEKE_VISC is True.")

549 endif

550

551 ! Set the halo sizes used for the thickness-point viscosities.

552 if (cs%use_Leithy .or. cs%debug) then

553 js_kh = js-1 ; je_kh = je+1 ; is_kh = is-1 ; ie_kh = ie+1

554 else

555 js_kh = jsq ; je_kh = je+1 ; is_kh = isq ; ie_kh = ie+1

556 endif

557

558 ! Set the halo sizes used for the vorticity calculations.

559 if ((cs%Leith_Kh) .or. (cs%Leith_Ah) .or. (cs%use_Leithy)) then

560 js_vort = js_kh-2 ; je_vort = jeq+2 ; is_vort = is_kh-2 ; ie_vort = ieq+2

561 if ((g%isc-g%isd < 3) .or. (g%isc-g%isd < 3)) call mom_error(fatal, &

562 "The minimum halo size is 3 when a Leith viscosity is being used.")

563 else

564 js_vort = js-2 ; je_vort = jeq+1 ; is_vort = is-2 ; ie_vort = ieq+1

565 endif

566

567 if (cs%use_GME) then

568

569 ! Initialize diagnostic arrays with zeros

570 gme_coeff_h(:,:,:) = 0.0

571 gme_coeff_q(:,:,:) = 0.0

572 str_xx_gme(:,:) = 0.0

573 str_xy_gme(:,:) = 0.0

574

575 ! Get barotropic velocities and their gradients

576 call barotropic_get_tav(bt, ubtav, vbtav, g, us)

577

578 call pass_vector(ubtav, vbtav, g%Domain)

579 call pass_var(h, g%domain, halo=2)

580

581 ! Calculate the barotropic horizontal tension

582 do j=js-2,je+2 ; do i=is-2,ie+2

583 dudx_bt(i,j) = cs%DY_dxT(i,j)*((g%IdyCu(i,j) * ubtav(i,j)) - &

584 (g%IdyCu(i-1,j) * ubtav(i-1,j)))

585 dvdy_bt(i,j) = cs%DX_dyT(i,j)*((g%IdxCv(i,j) * vbtav(i,j)) - &

586 (g%IdxCv(i,j-1) * vbtav(i,j-1)))

587 enddo ; enddo

588 do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2

589 sh_xx_bt(i,j) = dudx_bt(i,j) - dvdy_bt(i,j)

590 enddo ; enddo

591

592 ! Components for the barotropic shearing strain

593 do j=jsq-2,jeq+2 ; do i=isq-2,ieq+2

594 dvdx_bt(i,j) = cs%DY_dxBu(i,j)*((vbtav(i+1,j)*g%IdyCv(i+1,j)) &

595 - (vbtav(i,j)*g%IdyCv(i,j)))

596 dudy_bt(i,j) = cs%DX_dyBu(i,j)*((ubtav(i,j+1)*g%IdxCu(i,j+1)) &

597 - (ubtav(i,j)*g%IdxCu(i,j)))

598 enddo ; enddo

599

600 if (cs%no_slip) then

601 do j=js-2,je+1 ; do i=is-2,ie+1

602 sh_xy_bt(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx_bt(i,j) + dudy_bt(i,j) )

603 enddo ; enddo

604 else

605 do j=js-2,je+1 ; do i=is-2,ie+1

606 sh_xy_bt(i,j) = g%mask2dBu(i,j) * ( dvdx_bt(i,j) + dudy_bt(i,j) )

607 enddo ; enddo

608 endif

609

610 do j=js-2,je+2 ; do i=is-2,ie+2

611 htot(i,j) = 0.0

612 enddo ; enddo

613 do k=1,nz ; do j=js-2,je+2 ; do i=is-2,ie+2

614 htot(i,j) = htot(i,j) + h(i,j,k)

615 enddo ; enddo ; enddo

616

617 i_gme_h0 = 1.0 / cs%GME_h0

618 do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2

619 boundary_mask_h = (g%mask2dCu(i,j) * g%mask2dCu(i-1,j)) * (g%mask2dCv(i,j) * g%mask2dCv(i,j-1))

620 grad_vel_mag_bt_h = g%mask2dT(i,j) * boundary_mask_h * (dudx_bt(i,j)**2 + dvdy_bt(i,j)**2 + &

621 (0.25*((dvdx_bt(i,j)+dvdx_bt(i-1,j-1)) + (dvdx_bt(i,j-1)+dvdx_bt(i-1,j))))**2 + &

622 (0.25*((dudy_bt(i,j)+dudy_bt(i-1,j-1)) + (dudy_bt(i,j-1)+dudy_bt(i-1,j))))**2)

623 ! Probably the following test could be simplified to

624 ! if (boundary_mask_h * G%mask2dT(I,J) > 0.0) then

625 if (grad_vel_mag_bt_h > 0.0) then

626 gme_effic_h(i,j) = cs%GME_efficiency * g%mask2dT(i,j) * (min(htot(i,j) * i_gme_h0, 1.0)**2)

627 else

628 gme_effic_h(i,j) = 0.0

629 endif

630 enddo ; enddo

631

632 do j=js-2,je+1 ; do i=is-2,ie+1

633 boundary_mask_q = (g%mask2dCv(i,j) * g%mask2dCv(i+1,j)) * (g%mask2dCu(i,j) * g%mask2dCu(i,j+1))

634 grad_vel_mag_bt_q = g%mask2dBu(i,j) * boundary_mask_q * (dvdx_bt(i,j)**2 + dudy_bt(i,j)**2 + &

635 (0.25*((dudx_bt(i,j)+dudx_bt(i+1,j+1)) + (dudx_bt(i,j+1)+dudx_bt(i+1,j))))**2 + &

636 (0.25*((dvdy_bt(i,j)+dvdy_bt(i+1,j+1)) + (dvdy_bt(i,j+1)+dvdy_bt(i+1,j))))**2)

637 ! Probably the following test could be simplified to

638 ! if (boundary_mask_q * G%mask2dBu(I,J) > 0.0) then

639 if (grad_vel_mag_bt_q > 0.0) then

640 h_arith_q = 0.25 * ((htot(i,j) + htot(i+1,j+1)) + (htot(i+1,j) + htot(i,j+1)))

641 gme_effic_q(i,j) = cs%GME_efficiency * g%mask2dBu(i,j) * (min(h_arith_q * i_gme_h0, 1.0)**2)

642 else

643 gme_effic_q(i,j) = 0.0

644 endif

645 enddo ; enddo

646

647 call thickness_diffuse_get_kh(td, kh_u_gme, kh_v_gme, g, gv)

648

649 call pass_vector(kh_u_gme, kh_v_gme, g%domain, to_all+scalar_pair)

650

651 if (cs%debug) &

652 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)

653

654 endif ! use_GME

655

656 if (cs%use_Leithy) then

657 ! Smooth the velocity. Right now it happens twice. In the future

658 ! one might make the number of smoothing cycles a user-specified parameter

659 do k=1,nz

660 ! One call applies the filter twice

661 u_smooth(:,:,k) = u(:,:,k)

662 v_smooth(:,:,k) = v(:,:,k)

663 call smooth_x9_uv(g, u_smooth(:,:,k), v_smooth(:,:,k), zero_land=.false.)

664 enddo

665 call pass_vector(u_smooth, v_smooth, g%Domain)

666 endif

667

668 if (cs%use_QG_Leith_visc .and. ((cs%Leith_Kh) .or. (cs%Leith_Ah))) then

669 call thickness_to_dz(h, tv, dz, g, gv, us, halo_size=2)

670 ! Calculate isopycnal slopes that will be used for some forms of viscosity.

671 call calc_qg_slopes(h, tv, dt, g, gv, us, slope_x, slope_y, varmix, obc)

672 ! If the following halo update is added, the calculations in calc_QG_slopes could work on just

673 ! the computational domains, and some halo updates outside of this routine could be smaller.

674 ! call pass_vector(slope_x, slope_y, G%Domain, halo=2)

675 endif

676

677 !$OMP parallel do default(none) if (.not. CS%smooth_AH) &

678 !$OMP shared( &

679 !$OMP CS, G, GV, US, OBC, VarMix, MEKE, u, v, h, uh, vh, &

680 !$OMP is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, &

681 !$OMP is_vort, ie_vort, js_vort, je_vort, &

682 !$OMP is_Kh, ie_Kh, js_Kh, je_Kh, &

683 !$OMP apply_OBC, apply_OBC_strain, rescale_Kh, find_FrictWork, use_kh_struct, skeb_use_frict, &

684 !$OMP use_MEKE_Ku, use_MEKE_Au, u_smooth, v_smooth, use_cont_huv, slope_x, slope_y, dz, &

685 !$OMP backscat_subround, GME_effic_h, GME_effic_q, &

686 !$OMP h_neglect, h_neglect3, inv_PI3, inv_PI6, &

687 !$OMP diffu, diffv, Kh_h, Kh_q, Ah_h, Ah_q, FrictWork, FrictWork_bh, FrictWork_GME, &

688 !$OMP div_xx_h, sh_xx_h, vort_xy_q, sh_xy_q, GME_coeff_h, GME_coeff_q, &

689 !$OMP KH_u_GME, KH_v_GME, grid_Re_Kh, grid_Re_Ah, NoSt, ShSt, hu_cont, hv_cont, STOCH &

690 !$OMP ) &

691 !$OMP private( &

692 !$OMP i, j, k, n, tmp, &

693 !$OMP dudx, dudy, dvdx, dvdy, sh_xx, sh_xy, h_u, h_v, &

694 !$OMP Del2u, Del2v, DY_dxBu, DX_dyBu, sh_xx_bt, sh_xy_bt, &

695 !$OMP str_xx, str_xy, bhstr_xx, bhstr_xy, str_xx_GME, str_xy_GME, &

696 !$OMP vort_xy, vort_xy_dx, vort_xy_dy, div_xx, div_xx_dx, div_xx_dy, &

697 !$OMP grad_div_mag_h, grad_div_mag_q, grad_vort_mag_h, grad_vort_mag_q, &

698 !$OMP grad_vort, grad_vort_qg, grad_vort_mag_h_2d, grad_vort_mag_q_2d, &

699 !$OMP sh_xx_sq, sh_xy_sq, meke_res_fn, Shear_mag, Shear_mag_bc, vert_vort_mag, &

700 !$OMP h_min, hrat_min, visc_bound_rem, Kh_max_here, &

701 !$OMP grid_Ah, grid_Kh, d_Del2u, d_Del2v, d_str, &

702 !$OMP Kh, Ah, AhSm, AhLth, local_strain, Sh_F_pow, &

703 !$OMP dDel2vdx, dDel2udy, Del2vort_q, Del2vort_h, KE, &

704 !$OMP h2uq, h2vq, hu, hv, hq, FatH, RoScl, GME_coeff, &

705 !$OMP dudx_smooth, dudy_smooth, dvdx_smooth, dvdy_smooth, &

706 !$OMP vort_xy_smooth, vort_xy_dx_smooth, vort_xy_dy_smooth, &

707 !$OMP sh_xx_smooth, sh_xy_smooth, &

708 !$OMP vert_vort_mag_smooth, m_leithy, Ah_sq, AhLthy, &

709 !$OMP Kh_BS, str_xx_bs, str_xy_bs, bs_coeff_h, bs_coeff_q &

710 !$OMP ) &

711 !$OMP firstprivate( &

712 !$OMP visc_limit_h, visc_limit_h_frac, visc_limit_h_flag, &

713 !$OMP visc_limit_q, visc_limit_q_frac, visc_limit_q_flag &

714 !$OMP )

715 do k=1,nz

716

717 ! The following are the forms of the horizontal tension and horizontal

718 ! shearing strain advocated by Smagorinsky (1993) and discussed in

719 ! Griffies and Hallberg (2000).

720

721 ! NOTE: There is a ~1% speedup when the tension and shearing loops below

722 ! are fused (presumably due to shared access of Id[xy]C[uv]). However,

723 ! this breaks the center/vertex index case convention, and also evaluates

724 ! the dudx and dvdy terms beyond their valid bounds.

725 ! TODO: Explore methods for retaining both the syntax and speedup.

726

727 ! Calculate horizontal tension

728 do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2

729 dudx(i,j) = cs%DY_dxT(i,j)*((g%IdyCu(i,j) * u(i,j,k)) - &

730 (g%IdyCu(i-1,j) * u(i-1,j,k)))

731 dvdy(i,j) = cs%DX_dyT(i,j)*((g%IdxCv(i,j) * v(i,j,k)) - &

732 (g%IdxCv(i,j-1) * v(i,j-1,k)))

733 sh_xx(i,j) = dudx(i,j) - dvdy(i,j)

734 enddo ; enddo

735

736 ! Components for the shearing strain

737 do j=js_vort,je_vort ; do i=is_vort,ie_vort

738 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)))

739 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)))

740 enddo ; enddo

741

742 if (cs%use_Leithy) then

743 ! Calculate horizontal tension from smoothed velocity

744 do j=jsq,jeq+1 ; do i=isq,ieq+1

745 dudx_smooth(i,j) = cs%DY_dxT(i,j)*((g%IdyCu(i,j) * u_smooth(i,j,k)) - &

746 (g%IdyCu(i-1,j) * u_smooth(i-1,j,k)))

747 dvdy_smooth(i,j) = cs%DX_dyT(i,j)*((g%IdxCv(i,j) * v_smooth(i,j,k)) - &

748 (g%IdxCv(i,j-1) * v_smooth(i,j-1,k)))

749 sh_xx_smooth(i,j) = dudx_smooth(i,j) - dvdy_smooth(i,j)

750 enddo ; enddo

751

752 ! Components for the shearing strain from smoothed velocity

753 do j=js_kh-1,je_kh ; do i=is_kh-1,ie_kh

754 dvdx_smooth(i,j) = cs%DY_dxBu(i,j) * &

755 ((v_smooth(i+1,j,k)*g%IdyCv(i+1,j)) - (v_smooth(i,j,k)*g%IdyCv(i,j)))

756 dudy_smooth(i,j) = cs%DX_dyBu(i,j) * &

757 ((u_smooth(i,j+1,k)*g%IdxCu(i,j+1)) - (u_smooth(i,j,k)*g%IdxCu(i,j)))

758 enddo ; enddo

759 endif ! use Leith+E

760

761 if (cs%id_normstress > 0) then

762 do j=js,je ; do i=is,ie

763 nost(i,j,k) = sh_xx(i,j)

764 enddo ; enddo

765 endif

766

767 ! Interpolate the thicknesses to velocity points.

768 ! The extra wide halos are to accommodate the cross-corner-point projections

769 ! in OBCs, which are not ordinarily be necessary, and might not be necessary

770 ! even with OBCs if the accelerations are zeroed at OBC points, in which

771 ! case the j-loop for h_u could collapse to j=js=1,je+1. -RWH

772 if (use_cont_huv) then

773 do j=js-2,je+2 ; do i=isq-1,ieq+1

774 h_u(i,j) = hu_cont(i,j,k)

775 enddo ; enddo

776 do j=jsq-1,jeq+1 ; do i=is-2,ie+2

777 h_v(i,j) = hv_cont(i,j,k)

778 enddo ; enddo

779 elseif (cs%use_land_mask) then

780 do j=js-2,je+2 ; do i=is-2,ieq+1

781 h_u(i,j) = 0.5 * (g%mask2dT(i,j)*h(i,j,k) + g%mask2dT(i+1,j)*h(i+1,j,k))

782 enddo ; enddo

783 do j=js-2,jeq+1 ; do i=is-2,ie+2

784 h_v(i,j) = 0.5 * (g%mask2dT(i,j)*h(i,j,k) + g%mask2dT(i,j+1)*h(i,j+1,k))

785 enddo ; enddo

786 else

787 do j=js-2,je+2 ; do i=is-2,ieq+1

788 h_u(i,j) = 0.5 * (h(i,j,k) + h(i+1,j,k))

789 enddo ; enddo

790 do j=js-2,jeq+1 ; do i=is-2,ie+2

791 h_v(i,j) = 0.5 * (h(i,j,k) + h(i,j+1,k))

792 enddo ; enddo

793 endif

794

795 ! Adjust contributions to shearing strain and interpolated values of

796 ! thicknesses on open boundaries.

797 if (apply_obc) then ; do n=1,obc%number_of_segments

798 j = obc%segment(n)%HI%JsdB ; i = obc%segment(n)%HI%IsdB

799 if (apply_obc_strain) then

800 if (obc%segment(n)%is_N_or_S .and. (j >= js_vort) .and. (j <= je_vort)) then

801 do i = max(obc%segment(n)%HI%IsdB,is_vort), min(obc%segment(n)%HI%IedB,ie_vort)

802 select case (obc%strain_config)

803 case (obc_strain_zero)

804 dvdx(i,j) = 0. ; dudy(i,j) = 0.

805 case (obc_strain_freeslip)

806 dudy(i,j) = 0.

807 case (obc_strain_computed)

808 if (obc%segment(n)%direction == obc_direction_n) then

809 dudy(i,j) = 2.0*cs%DX_dyBu(i,j)* &

810 (obc%segment(n)%tangential_vel(i,j,k) - u(i,j,k))*g%IdxCu(i,j)

811 else

812 dudy(i,j) = 2.0*cs%DX_dyBu(i,j)* &

813 (u(i,j+1,k) - obc%segment(n)%tangential_vel(i,j,k))*g%IdxCu(i,j+1)

814 endif

815 case (obc_strain_specified)

816 if (obc%segment(n)%direction == obc_direction_n) then

817 dudy(i,j) = cs%DX_dyBu(i,j)*obc%segment(n)%tangential_grad(i,j,k)*g%IdxCu(i,j)*g%dxBu(i,j)

818 else

819 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)

820 endif

821 end select

822 if (cs%use_Leithy) then

823 dvdx_smooth(i,j) = dvdx(i,j)

824 dudy_smooth(i,j) = dudy(i,j)

825 endif

826 enddo

827 elseif (obc%segment(n)%is_E_or_W .and. (i >= is_vort) .and. (i <= ie_vort)) then

828 do j = max(obc%segment(n)%HI%JsdB,js_vort), min(obc%segment(n)%HI%JedB,je_vort)

829 select case (obc%strain_config)

830 case (obc_strain_zero)

831 dvdx(i,j) = 0. ; dudy(i,j) = 0.

832 case (obc_strain_freeslip)

833 dvdx(i,j) = 0.

834 case (obc_strain_computed)

835 if (obc%segment(n)%direction == obc_direction_e) then

836 dvdx(i,j) = 2.0*cs%DY_dxBu(i,j)* &

837 (obc%segment(n)%tangential_vel(i,j,k) - v(i,j,k))*g%IdyCv(i,j)

838 else

839 dvdx(i,j) = 2.0*cs%DY_dxBu(i,j)* &

840 (v(i+1,j,k) - obc%segment(n)%tangential_vel(i,j,k))*g%IdyCv(i+1,j)

841 endif

842 case (obc_strain_specified)

843 if (obc%segment(n)%direction == obc_direction_e) then

844 dvdx(i,j) = cs%DY_dxBu(i,j)*obc%segment(n)%tangential_grad(i,j,k)*g%IdyCv(i,j)*g%dxBu(i,j)

845 else

846 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)

847 endif

848 end select

849 if (cs%use_Leithy) then

850 dvdx_smooth(i,j) = dvdx(i,j)

851 dudy_smooth(i,j) = dudy(i,j)

852 endif

853 enddo

854 endif

855 endif

856

857 if (obc%segment(n)%direction == obc_direction_n) then

858 ! There are extra wide halos here to accommodate the cross-corner-point

859 ! OBC projections, but they might not be necessary if the accelerations

860 ! are always zeroed out at OBC points, in which case the i-loop below

861 ! becomes do i=is-1,ie+1. -RWH

862 if ((j >= js-2) .and. (j <= jeq+1)) then

863 do i = max(is-2,obc%segment(n)%HI%isd), min(ie+2,obc%segment(n)%HI%ied)

864 h_v(i,j) = h(i,j,k)

865 enddo

866 endif

867 elseif (obc%segment(n)%direction == obc_direction_s) then

868 if ((j >= js-2) .and. (j <= jeq+1)) then

869 do i = max(is-2,obc%segment(n)%HI%isd), min(ie+2,obc%segment(n)%HI%ied)

870 h_v(i,j) = h(i,j+1,k)

871 enddo

872 endif

873 elseif (obc%segment(n)%direction == obc_direction_e) then

874 if ((i >= is-2) .and. (i <= ieq+1)) then

875 do j = max(js-2,obc%segment(n)%HI%jsd), min(je+2,obc%segment(n)%HI%jed)

876 h_u(i,j) = h(i,j,k)

877 enddo

878 endif

879 elseif (obc%segment(n)%direction == obc_direction_w) then

880 if ((i >= is-2) .and. (i <= ieq+1)) then

881 do j = max(js-2,obc%segment(n)%HI%jsd), min(je+2,obc%segment(n)%HI%jed)

882 h_u(i,j) = h(i+1,j,k)

883 enddo

884 endif

885 endif

886 enddo ; endif

887 ! Now project thicknesses across corner points on OBCs.

888 if (apply_obc) then ; do n=1,obc%number_of_segments

889 j = obc%segment(n)%HI%JsdB ; i = obc%segment(n)%HI%IsdB

890 if (obc%segment(n)%direction == obc_direction_n) then

891 if ((j >= js-2) .and. (j <= je)) then

892 do i = max(is-2,obc%segment(n)%HI%IsdB), min(ieq+1,obc%segment(n)%HI%IedB)

893 h_u(i,j+1) = h_u(i,j)

894 enddo

895 endif

896 elseif (obc%segment(n)%direction == obc_direction_s) then

897 if ((j >= js-1) .and. (j <= je+1)) then

898 do i = max(is-2,obc%segment(n)%HI%isd), min(ieq+1,obc%segment(n)%HI%ied)

899 h_u(i,j) = h_u(i,j+1)

900 enddo

901 endif

902 elseif (obc%segment(n)%direction == obc_direction_e) then

903 if ((i >= is-2) .and. (i <= ie)) then

904 do j = max(js-2,obc%segment(n)%HI%jsd), min(jeq+1,obc%segment(n)%HI%jed)

905 h_v(i+1,j) = h_v(i,j)

906 enddo

907 endif

908 elseif (obc%segment(n)%direction == obc_direction_w) then

909 if ((i >= is-1) .and. (i <= ie+1)) then

910 do j = max(js-2,obc%segment(n)%HI%jsd), min(jeq+1,obc%segment(n)%HI%jed)

911 h_v(i,j) = h_v(i+1,j)

912 enddo

913 endif

914 endif

915 enddo ; endif

916

917 ! Shearing strain (including no-slip boundary conditions at the 2-D land-sea mask).

918 ! dudy and dvdx include modifications at OBCs from above.

919 if (cs%no_slip) then

920 do j=js-2,jeq+1 ; do i=is-2,ieq+1

921 sh_xy(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx(i,j) + dudy(i,j) )

922 if (cs%id_shearstress > 0) shst(i,j,k) = sh_xy(i,j)

923 enddo ; enddo

924 else

925 do j=js-2,jeq+1 ; do i=is-2,ieq+1

926 sh_xy(i,j) = g%mask2dBu(i,j) * ( dvdx(i,j) + dudy(i,j) )

927 if (cs%id_shearstress > 0) shst(i,j,k) = sh_xy(i,j)

928 enddo ; enddo

929 endif

930

931 if (cs%use_Leithy) then

932 ! Shearing strain (including no-slip boundary conditions at the 2-D land-sea mask).

933 ! dudy_smooth and dvdx_smooth do not (yet) include modifications at OBCs from above.

934 if (cs%no_slip) then

935 do j=js-1,jeq ; do i=is-1,ieq

936 sh_xy_smooth(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx_smooth(i,j) + dudy_smooth(i,j) )

937 enddo ; enddo

938 else

939 do j=js-1,jeq ; do i=is-1,ieq

940 sh_xy_smooth(i,j) = g%mask2dBu(i,j) * ( dvdx_smooth(i,j) + dudy_smooth(i,j) )

941 enddo ; enddo

942 endif

943 endif ! use Leith+E

944

945 ! Evaluate Del2u = x.Div(Grad u) and Del2v = y.Div( Grad u)

946 if (cs%biharmonic) then

947 do j=js-1,jeq+1 ; do i=isq-1,ieq+1

948 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))) + &

949 cs%Idxdy2u(i,j) * ((cs%dy2h(i+1,j)*sh_xx(i+1,j)) - (cs%dy2h(i,j)*sh_xx(i,j)))

950 enddo ; enddo

951 do j=jsq-1,jeq+1 ; do i=is-1,ieq+1

952 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))) - &

953 cs%Idx2dyCv(i,j) * ((cs%dx2h(i,j+1)*sh_xx(i,j+1)) - (cs%dx2h(i,j)*sh_xx(i,j)))

954 enddo ; enddo

955 if (apply_obc) then ; if (obc%zero_biharmonic) then

956 do n=1,obc%number_of_segments

957 i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB

958 if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1)) then

959 do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied

960 del2v(i,j) = 0.

961 enddo

962 elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1)) then

963 do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed

964 del2u(i,j) = 0.

965 enddo

966 endif

967 enddo

968 endif ; endif

969 endif

970

971 ! Vorticity

972 if ((cs%Leith_Kh) .or. (cs%Leith_Ah) .or. (cs%use_Leithy) .or. (cs%id_vort_xy_q>0) .or. cs%use_ZB2020) then

973 if (cs%no_slip) then

974 do j=js_vort,je_vort ; do i=is_vort,ie_vort

975 vort_xy(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx(i,j) - dudy(i,j) )

976 enddo ; enddo

977 else

978 if (cs%use_circulation) then

979 do j=js_vort,je_vort ; do i=is_vort,ie_vort

980 vort_xy(i,j) = g%mask2dBu(i,j) * g%IareaBu(i,j) * ( &

981 ((v(i+1,j,k)*g%dyCv(i+1,j)) - (v(i,j,k)*g%dyCv(i,j))) &

982 - ((u(i,j+1,k)*g%dxCu(i,j+1)) - (u(i,j,k)*g%dxCu(i,j))) &

983 )

984 enddo ; enddo

985 else

986 do j=js_vort,je_vort ; do i=is_vort,ie_vort

987 vort_xy(i,j) = g%mask2dBu(i,j) * ( dvdx(i,j) - dudy(i,j) )

988 enddo ; enddo

989 endif

990 endif

991 endif

992

993 if (cs%use_Leithy) then

994 if (cs%no_slip) then

995 do j=js_kh-1,je_kh ; do i=is_kh-1,ie_kh

996 vort_xy_smooth(i,j) = (2.0-g%mask2dBu(i,j)) * ( dvdx_smooth(i,j) - dudy_smooth(i,j) )

997 enddo ; enddo

998 else

999 do j=js_kh-1,je_kh ; do i=is_kh-1,ie_kh

1000 vort_xy_smooth(i,j) = g%mask2dBu(i,j) * ( dvdx_smooth(i,j) - dudy_smooth(i,j) )

1001 enddo ; enddo

1002 endif

1003 endif

1004

1005

1006 if ((cs%Leith_Kh) .or. (cs%Leith_Ah) .or. (cs%use_Leithy)) then

1007

1008 ! Vorticity gradient

1009 do j=js-2,je_kh ; do i=is_kh-1,ie_kh+1

1010 dy_dxbu = g%dyBu(i,j) * g%IdxBu(i,j)

1011 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)))

1012 enddo ; enddo

1013

1014 do j=js_kh-1,je_kh+1 ; do i=is-2,ie_kh

1015 dx_dybu = g%dxBu(i,j) * g%IdyBu(i,j)

1016 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)))

1017 enddo ; enddo

1018

1019 if (cs%use_Leithy) then

1020 ! Gradient of smoothed vorticity

1021 do j=js_kh-1,je_kh ; do i=is_kh,ie_kh

1022 dy_dxbu = g%dyBu(i,j) * g%IdxBu(i,j)

1023 vort_xy_dx_smooth(i,j) = dy_dxbu * &

1024 ((vort_xy_smooth(i,j) * g%IdyCu(i,j)) - (vort_xy_smooth(i-1,j) * g%IdyCu(i-1,j)))

1025 enddo ; enddo

1026

1027 do j=js_kh,je_kh ; do i=is_kh-1,ie_kh

1028 dx_dybu = g%dxBu(i,j) * g%IdyBu(i,j)

1029 vort_xy_dy_smooth(i,j) = dx_dybu * &

1030 ((vort_xy_smooth(i,j) * g%IdxCv(i,j)) - (vort_xy_smooth(i,j-1) * g%IdxCv(i,j-1)))

1031 enddo ; enddo

1032 endif ! If Leithy

1033

1034 ! Laplacian of vorticity

1035 ! if (CS%Leith_Ah .or. CS%use_Leithy) then

1036 do j=js_kh-1,je_kh ; do i=is_kh-1,ie_kh

1037 dy_dxbu = g%dyBu(i,j) * g%IdxBu(i,j)

1038 dx_dybu = g%dxBu(i,j) * g%IdyBu(i,j)

1039

1040 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))) + &

1041 dx_dybu * ((vort_xy_dy(i,j+1) * g%IdyCu(i,j+1)) - (vort_xy_dy(i,j) * g%IdyCu(i,j)))

1042 enddo ; enddo

1043 ! endif

1044

1045 if (cs%modified_Leith) then

1046

1047 ! Divergence

1048 do j=js_kh-1,je_kh+1 ; do i=is_kh-1,ie_kh+1

1049 div_xx(i,j) = dudx(i,j) + dvdy(i,j)

1050 enddo ; enddo

1051

1052 ! Divergence gradient

1053 do j=js-1,je+1 ; do i=is_kh-1,ie_kh

1054 div_xx_dx(i,j) = g%IdxCu(i,j)*(div_xx(i+1,j) - div_xx(i,j))

1055 enddo ; enddo

1056 do j=js_kh-1,je_kh ; do i=is-1,ie+1

1057 div_xx_dy(i,j) = g%IdyCv(i,j)*(div_xx(i,j+1) - div_xx(i,j))

1058 enddo ; enddo

1059

1060 ! Magnitude of divergence gradient

1061 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1062 grad_div_mag_h(i,j) = sqrt(((0.5*(div_xx_dx(i,j) + div_xx_dx(i-1,j)))**2) + &

1063 ((0.5*(div_xx_dy(i,j) + div_xx_dy(i,j-1)))**2))

1064 enddo ; enddo

1065 do j=js-1,jeq ; do i=is-1,ieq

1066 grad_div_mag_q(i,j) = sqrt(((0.5*(div_xx_dx(i,j) + div_xx_dx(i,j+1)))**2) + &

1067 ((0.5*(div_xx_dy(i,j) + div_xx_dy(i+1,j)))**2))

1068 enddo ; enddo

1069

1070 else

1071

1072 do j=js-1,je+1 ; do i=is_kh-1,ie_kh

1073 div_xx_dx(i,j) = 0.0

1074 enddo ; enddo

1075 do j=js_kh-1,je_kh ; do i=is-1,ie+1

1076 div_xx_dy(i,j) = 0.0

1077 enddo ; enddo

1078 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1079 grad_div_mag_h(i,j) = 0.0

1080 enddo ; enddo

1081 do j=js-1,jeq ; do i=is-1,ieq

1082 grad_div_mag_q(i,j) = 0.0

1083 enddo ; enddo

1084

1085 endif ! CS%modified_Leith

1086

1087 ! Add in beta for the Leith viscosity

1088 if (cs%use_beta_in_Leith) then

1089 do j=js-2,jeq+1 ; do i=is-1,ie+1

1090 vort_xy_dx(i,j) = vort_xy_dx(i,j) + 0.5 * ( g%dF_dx(i,j) + g%dF_dx(i,j+1))

1091 enddo ; enddo

1092 do j=js-1,je+1 ; do i=is-2,ieq+1

1093 vort_xy_dy(i,j) = vort_xy_dy(i,j) + 0.5 * ( g%dF_dy(i,j) + g%dF_dy(i+1,j))

1094 enddo ; enddo

1095 endif ! CS%use_beta_in_Leith

1096

1097 if (cs%use_QG_Leith_visc) then

1098

1099 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1100 grad_vort_mag_h_2d(i,j) = sqrt(((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i,j-1)))**2) + &

1101 ((0.5*(vort_xy_dy(i,j) + vort_xy_dy(i-1,j)))**2) )

1102 enddo ; enddo

1103 do j=js-1,jeq ; do i=is-1,ieq

1104 grad_vort_mag_q_2d(i,j) = sqrt(((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i+1,j)))**2) + &

1105 ((0.5*(vort_xy_dy(i,j) + vort_xy_dy(i,j+1)))**2) )

1106 enddo ; enddo

1107

1108 ! This accumulates terms, some of which are in VarMix.

1109 call calc_qg_leith_viscosity(varmix, g, gv, us, h, dz, k, div_xx_dx, div_xx_dy, &

1110 slope_x, slope_y, vort_xy_dx, vort_xy_dy)

1111

1112 endif

1113

1114 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1115 grad_vort_mag_h(i,j) = sqrt(((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i,j-1)))**2) + &

1116 ((0.5*(vort_xy_dy(i,j) + vort_xy_dy(i-1,j)))**2) )

1117 enddo ; enddo

1118 do j=js-1,jeq ; do i=is-1,ieq

1119 grad_vort_mag_q(i,j) = sqrt(((0.5*(vort_xy_dx(i,j) + vort_xy_dx(i+1,j)))**2) + &

1120 ((0.5*(vort_xy_dy(i,j) + vort_xy_dy(i,j+1)))**2) )

1121 enddo ; enddo

1122

1123 if (cs%use_Leithy) then

1124 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1125 vert_vort_mag_smooth(i,j) = sqrt(((0.5*(vort_xy_dx_smooth(i,j) + &

1126 vort_xy_dx_smooth(i,j-1)))**2) + &

1127 ((0.5*(vort_xy_dy_smooth(i,j) + &

1128 vort_xy_dy_smooth(i-1,j)))**2) )

1129 enddo ; enddo

1130 endif ! Leithy

1131

1132 endif ! CS%Leith_Kh

1133

1134 if ((cs%Smagorinsky_Kh) .or. (cs%Smagorinsky_Ah)) then

1135 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1136 sh_xx_sq = sh_xx(i,j)**2

1137 sh_xy_sq = 0.25 * ( ((sh_xy(i-1,j-1)**2) + (sh_xy(i,j)**2)) &

1138 + ((sh_xy(i-1,j)**2) + (sh_xy(i,j-1)**2)) )

1139 shear_mag(i,j) = sqrt(sh_xx_sq + sh_xy_sq)

1140 enddo ; enddo

1141 endif

1142

1143 if (cs%bound_Ah .or. cs%bound_Kh) then

1144 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1145 h_min = min(h_u(i,j), h_u(i-1,j), h_v(i,j), h_v(i,j-1))

1146 hrat_min(i,j) = min(1.0, h_min / (h(i,j,k) + h_neglect))

1147 enddo ; enddo

1148 endif

1149

1150 if (cs%Laplacian) then

1151 ! Determine the Laplacian viscosity at h points, using the

1152 ! largest value from several parameterizations. Also get

1153 ! the Laplacian component of str_xx.

1154

1155 if ((cs%Leith_Kh) .or. (cs%Leith_Ah) .or. (cs%use_Leithy)) then

1156 if (cs%use_QG_Leith_visc) then

1157 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1158 grad_vort = grad_vort_mag_h(i,j) + grad_div_mag_h(i,j)

1159 grad_vort_qg = 3. * grad_vort_mag_h_2d(i,j)

1160 vert_vort_mag(i,j) = min(grad_vort, grad_vort_qg)

1161 enddo ; enddo

1162 else

1163 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1164 vert_vort_mag(i,j) = grad_vort_mag_h(i,j) + grad_div_mag_h(i,j)

1165 enddo ; enddo

1166 endif

1167 endif

1168

1169 ! Static (pre-computed) background viscosity

1170 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1171 kh(i,j) = cs%Kh_bg_xx(i,j)

1172 enddo ; enddo

1173

1174 ! NOTE: The following do-block can be decomposed and vectorized after the

1175 ! stack size has been reduced.

1176 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1177 if (cs%add_LES_viscosity) then

1178 if (cs%Smagorinsky_Kh) &

1179 kh(i,j) = kh(i,j) + cs%Laplac2_const_xx(i,j) * shear_mag(i,j)

1180 if (cs%Leith_Kh) &

1181 kh(i,j) = kh(i,j) + cs%Laplac3_const_xx(i,j) * vert_vort_mag(i,j) * inv_pi3

1182 else

1183 if (cs%Smagorinsky_Kh) &

1184 kh(i,j) = max(kh(i,j), cs%Laplac2_const_xx(i,j) * shear_mag(i,j))

1185 if (cs%Leith_Kh) &

1186 kh(i,j) = max(kh(i,j), cs%Laplac3_const_xx(i,j) * vert_vort_mag(i,j) * inv_pi3)

1187 endif

1188 enddo ; enddo

1189

1190 ! All viscosity contributions above are subject to resolution scaling

1191

1192 if (rescale_kh) then

1193 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1194 kh(i,j) = varmix%Res_fn_h(i,j) * kh(i,j)

1195 enddo ; enddo

1196 endif

1197

1198 ! Place a floor on the viscosity, if desired.

1199 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1200 kh(i,j) = max(kh(i,j), cs%Kh_bg_min)

1201 enddo ; enddo

1202

1203 if (use_meke_ku .and. .not. cs%EY24_EBT_BS) then

1204 ! *Add* the MEKE contribution (which might be negative)

1205 if (use_kh_struct) then

1206 if (cs%res_scale_MEKE) then

1207 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1208 kh(i,j) = kh(i,j) + meke%Ku(i,j) * varmix%Res_fn_h(i,j) * varmix%BS_struct(i,j,k)

1209 enddo ; enddo

1210 else

1211 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1212 kh(i,j) = kh(i,j) + meke%Ku(i,j) * varmix%BS_struct(i,j,k)

1213 enddo ; enddo

1214 endif

1215 else

1216 if (cs%res_scale_MEKE) then

1217 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1218 kh(i,j) = kh(i,j) + meke%Ku(i,j) * varmix%Res_fn_h(i,j)

1219 enddo ; enddo

1220 else

1221 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1222 kh(i,j) = kh(i,j) + meke%Ku(i,j)

1223 enddo ; enddo

1224 endif

1225 endif

1226 endif

1227

1228 if (cs%anisotropic) then

1229 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1230 ! *Add* the tension component of anisotropic viscosity

1231 kh(i,j) = kh(i,j) + cs%Kh_aniso * (1. - cs%n1n2_h(i,j)**2)

1232 enddo ; enddo

1233 endif

1234

1235 ! Newer method of bounding for stability

1236 if ((cs%bound_Kh) .and. (cs%bound_Ah)) then

1237 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1238 visc_bound_rem(i,j) = 1.0

1239 kh_max_here = hrat_min(i,j) * cs%Kh_Max_xx(i,j)

1240 if (kh(i,j) >= kh_max_here) then

1241 visc_bound_rem(i,j) = 0.0

1242 kh(i,j) = kh_max_here

1243 elseif ((kh(i,j) > 0.0) .or. (cs%backscatter_underbound .and. (kh_max_here > 0.0))) then

1244 visc_bound_rem(i,j) = 1.0 - kh(i,j) / kh_max_here

1245 endif

1246 enddo ; enddo

1247 elseif (cs%bound_Kh) then

1248 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1249 kh(i,j) = min(kh(i,j), hrat_min(i,j) * cs%Kh_Max_xx(i,j))

1250 enddo ; enddo

1251 endif

1252

1253 ! In Leith+E parameterization Kh is computed after Ah in the biharmonic loop.

1254 ! The harmonic component of str_xx is added in the biharmonic loop.

1255 if (cs%use_Leithy) then

1256 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1257 kh(i,j) = 0.

1258 enddo ; enddo

1259 endif

1260

1261 if (cs%id_Kh_h>0 .or. cs%debug) then

1262 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1263 kh_h(i,j,k) = kh(i,j)

1264 enddo ; enddo

1265 endif

1266

1267 if (cs%id_grid_Re_Kh>0) then

1268 do j=js,je ; do i=is,ie

1269 ke = 0.125*(((u(i,j,k)+u(i-1,j,k))**2) + ((v(i,j,k)+v(i,j-1,k))**2))

1270 grid_kh = max(kh(i,j), cs%min_grid_Kh)

1271 grid_re_kh(i,j,k) = (sqrt(ke) * sqrt(cs%grid_sp_h2(i,j))) / grid_kh

1272 enddo ; enddo

1273 endif

1274

1275 if (cs%id_div_xx_h>0) then

1276 do j=js,je ; do i=is,ie

1277 div_xx_h(i,j,k) = dudx(i,j) + dvdy(i,j)

1278 enddo ; enddo

1279 endif

1280

1281 if (cs%id_sh_xx_h>0) then

1282 do j=js,je ; do i=is,ie

1283 sh_xx_h(i,j,k) = sh_xx(i,j)

1284 enddo ; enddo

1285 endif

1286

1287 do j=jsq,jeq+1 ; do i=isq,ieq+1

1288 str_xx(i,j) = -kh(i,j) * sh_xx(i,j)

1289 enddo ; enddo

1290 else

1291 do j=jsq,jeq+1 ; do i=isq,ieq+1

1292 str_xx(i,j) = 0.0

1293 enddo ; enddo

1294 endif ! Get Kh at h points and get Laplacian component of str_xx

1295

1296 if (cs%anisotropic) then

1297 do j=jsq,jeq+1 ; do i=isq,ieq+1

1298 ! Shearing-strain averaged to h-points

1299 local_strain = 0.25 * ( (sh_xy(i,j) + sh_xy(i-1,j-1)) + (sh_xy(i-1,j) + sh_xy(i,j-1)) )

1300 ! *Add* the shear-strain contribution to the xx-component of stress

1301 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

1302 enddo ; enddo

1303 endif

1304

1305 if (cs%biharmonic) then

1306 ! Determine the biharmonic viscosity at h points, using the

1307 ! largest value from several parameterizations. Also get the

1308 ! biharmonic component of str_xx.

1309 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1310 ah(i,j) = cs%Ah_bg_xx(i,j)

1311 enddo ; enddo

1312

1313 if ((cs%Smagorinsky_Ah) .or. (cs%Leith_Ah) .or. (cs%use_Leithy)) then

1314 if (cs%Smagorinsky_Ah) then

1315 if (cs%bound_Coriolis) then

1316 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1317 ahsm = shear_mag(i,j) * (cs%Biharm_const_xx(i,j) &

1318 + cs%Biharm_const2_xx(i,j) * shear_mag(i,j))

1319 ah(i,j) = max(ah(i,j), ahsm)

1320 enddo ; enddo

1321 else

1322 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1323 ahsm = cs%Biharm_const_xx(i,j) * shear_mag(i,j)

1324 ah(i,j) = max(ah(i,j), ahsm)

1325 enddo ; enddo

1326 endif

1327 endif

1328

1329 if (cs%Leith_Ah) then

1330 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1331 del2vort_h = 0.25 * ((del2vort_q(i,j) + del2vort_q(i-1,j-1)) + &

1332 (del2vort_q(i-1,j) + del2vort_q(i,j-1)))

1333 ahlth = cs%Biharm6_const_xx(i,j) * abs(del2vort_h) * inv_pi6

1334 ah(i,j) = max(ah(i,j), ahlth)

1335 enddo ; enddo

1336 endif

1337

1338 if (cs%use_Leithy) then

1339 ! Get m_leithy

1340 if (cs%smooth_Ah) m_leithy(:,:) = 0.0 ! This is here to initialize domain edge halo values.

1341 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1342 del2vort_h = 0.25 * ((del2vort_q(i,j) + del2vort_q(i-1,j-1)) + &

1343 (del2vort_q(i-1,j) + del2vort_q(i,j-1)))

1344 ahlth = cs%Biharm6_const_xx(i,j) * inv_pi6 * abs(del2vort_h)

1345 if (ahlth <= cs%Ah_bg_xx(i,j)) then

1346 m_leithy(i,j) = 0.0

1347 else

1348 if ((cs%m_const_leithy(i,j)*vert_vort_mag(i,j)) < abs(vort_xy_smooth(i,j))) then

1349 m_leithy(i,j) = cs%c_K * (vert_vort_mag(i,j) / vort_xy_smooth(i,j))**2

1350 else

1351 m_leithy(i,j) = cs%m_leithy_max(i,j)

1352 endif

1353 m_leithy(i,j) = g%mask2dBu(i,j) * m_leithy(i,j)

1354 endif

1355 enddo ; enddo

1356

1357 if (cs%smooth_Ah) then

1358 ! Smooth m_leithy. A single call smoothes twice.

1359 call pass_var(m_leithy, g%Domain, halo=2)

1360 call smooth_x9_h(g, m_leithy, zero_land=.true.)

1361 call pass_var(m_leithy, g%Domain)

1362 endif

1363 ! Get Ah

1364 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1365 del2vort_h = 0.25 * ((del2vort_q(i,j) + del2vort_q(i-1,j-1)) + &

1366 (del2vort_q(i-1,j) + del2vort_q(i,j-1)))

1367 ahlthy = cs%Biharm6_const_xx(i,j) * inv_pi6 * &

1368 sqrt(max(0.,del2vort_h**2 - m_leithy(i,j)*vert_vort_mag_smooth(i,j)**2))

1369 ah(i,j) = max(cs%Ah_bg_xx(i,j), ahlthy)

1370 enddo ; enddo

1371 if (cs%smooth_Ah) then

1372 ! Smooth Ah before applying upper bound. Square Ah, then smooth, then take its square root.

1373 ah_sq(:,:) = 0.0 ! This is here to initialize domain edge halo values.

1374 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1375 ah_sq(i,j) = ah(i,j)**2

1376 enddo ; enddo

1377 call pass_var(ah_sq, g%Domain, halo=2)

1378 ! A single call smoothes twice.

1379 call smooth_x9_h(g, ah_sq, zero_land=.false.)

1380 call pass_var(ah_sq, g%Domain)

1381 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1382 ah_h(i,j,k) = max(cs%Ah_bg_xx(i,j), sqrt(max(0., ah_sq(i,j))))

1383 ah(i,j) = ah_h(i,j,k)

1384 enddo ; enddo

1385 else

1386 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1387 ah_h(i,j,k) = ah(i,j)

1388 enddo ; enddo

1389 endif

1390 endif

1391

1392 endif ! Smagorinsky_Ah or Leith_Ah or Leith+E

1393

1394 if (use_meke_au) then

1395 ! *Add* the MEKE contribution

1396 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1397 ah(i,j) = ah(i,j) + meke%Au(i,j)

1398 enddo ; enddo

1399 endif

1400

1401 if (cs%Re_Ah > 0.0) then

1402 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1403 ke = 0.125*(((u(i,j,k)+u(i-1,j,k))**2) + ((v(i,j,k)+v(i,j-1,k))**2))

1404 ah(i,j) = sqrt(ke) * cs%Re_Ah_const_xx(i,j)

1405 enddo ; enddo

1406 endif

1407

1408 if (cs%bound_Ah) then

1409 if (cs%bound_Kh) then

1410 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1411 ah(i,j) = min(ah(i,j), visc_bound_rem(i,j) * hrat_min(i,j) * cs%Ah_Max_xx(i,j))

1412 enddo ; enddo

1413 else

1414 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1415 ah(i,j) = min(ah(i,j), hrat_min(i,j) * cs%Ah_Max_xx(i,j))

1416 enddo ; enddo

1417 endif

1418 endif

1419

1420 if (cs%EY24_EBT_BS) then

1421 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1422 tmp = cs%KS_coef * hrat_min(i,j) * cs%Ah_Max_xx_KS(i,j)

1423 visc_limit_h(i,j,k) = tmp

1424 visc_limit_h_frac(i,j,k) = ah(i,j) / (cs%KS_coef * hrat_min(i,j) * cs%Ah_Max_xx_KS(i,j))

1425 if (ah(i,j) >= tmp) then

1426 visc_limit_h_flag(i,j,k) = 1.

1427 endif

1428 enddo ; enddo

1429 endif

1430

1431 if ((cs%id_Ah_h>0) .or. cs%debug .or. cs%use_Leithy) then

1432 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1433 ah_h(i,j,k) = ah(i,j)

1434 enddo ; enddo

1435 endif

1436

1437 if (cs%use_Leithy) then

1438 ! Compute Leith+E Kh after bounds have been applied to Ah

1439 ! and after it has been smoothed. Kh = -m_leithy * Ah

1440 do j=js_kh,je_kh ; do i=is_kh,ie_kh

1441 kh(i,j) = -m_leithy(i,j) * ah(i,j)

1442 kh_h(i,j,k) = kh(i,j)

1443 enddo ; enddo

1444 endif

1445

1446 if (cs%id_grid_Re_Ah > 0) then

1447 do j=js,je ; do i=is,ie

1448 ke = 0.125 * (((u(i,j,k) + u(i-1,j,k))**2) + ((v(i,j,k) + v(i,j-1,k))**2))

1449 grid_ah = max(ah(i,j), cs%min_grid_Ah)

1450 grid_re_ah(i,j,k) = (sqrt(ke) * cs%grid_sp_h3(i,j)) / grid_ah

1451 enddo ; enddo

1452 endif

1453

1454 do j=jsq,jeq+1 ; do i=isq,ieq+1

1455 d_del2u = (g%IdyCu(i,j) * del2u(i,j)) - (g%IdyCu(i-1,j) * del2u(i-1,j))

1456 d_del2v = (g%IdxCv(i,j) * del2v(i,j)) - (g%IdxCv(i,j-1) * del2v(i,j-1))

1457 d_str = ah(i,j) * ((cs%DY_dxT(i,j) * d_del2u) - (cs%DX_dyT(i,j) * d_del2v))

1458

1459 str_xx(i,j) = str_xx(i,j) + d_str

1460

1461 if (cs%use_Leithy) str_xx(i,j) = str_xx(i,j) - kh(i,j) * sh_xx_smooth(i,j)

1462

1463 ! Keep a copy of the biharmonic contribution for backscatter parameterization

1464 bhstr_xx(i,j) = d_str * (h(i,j,k) * cs%reduction_xx(i,j))

1465 enddo ; enddo

1466 endif ! Get biharmonic coefficient at h points and biharmonic part of str_xx

1467

1468 ! Backscatter using MEKE

1469 if (cs%EY24_EBT_BS) then

1470 do j=jsq,jeq+1 ; do i=isq,ieq+1

1471 if (visc_limit_h_flag(i,j,k) > 0) then

1472 kh_bs(i,j) = 0.

1473 else

1474 if (use_kh_struct) then

1475 kh_bs(i,j) = meke%Ku(i,j) * varmix%BS_struct(i,j,k)

1476 else

1477 kh_bs(i,j) = meke%Ku(i,j)

1478 endif

1479 endif

1480 enddo ; enddo

1481

1482 do j=jsq,jeq+1 ; do i=isq,ieq+1

1483 str_xx_bs(i,j) = -kh_bs(i,j) * sh_xx(i,j)

1484 enddo ; enddo

1485

1486 if (cs%id_BS_coeff_h>0) then

1487 do j=jsq,jeq+1 ; do i=isq,ieq+1

1488 bs_coeff_h(i,j,k) = kh_bs(i,j)

1489 enddo ; enddo

1490 endif

1491

1492 do j=jsq,jeq+1 ; do i=isq,ieq+1

1493 str_xx(i,j) = str_xx(i,j) + str_xx_bs(i,j)

1494 enddo ; enddo

1495 endif ! Backscatter

1496

1497 if (cs%biharmonic) then

1498 ! Gradient of Laplacian, for use in bi-harmonic term

1499 do j=js-1,jeq ; do i=is-1,ieq

1500 ddel2vdx(i,j) = cs%DY_dxBu(i,j)*((del2v(i+1,j)*g%IdyCv(i+1,j)) - (del2v(i,j)*g%IdyCv(i,j)))

1501 ddel2udy(i,j) = cs%DX_dyBu(i,j)*((del2u(i,j+1)*g%IdxCu(i,j+1)) - (del2u(i,j)*g%IdxCu(i,j)))

1502 enddo ; enddo

1503 ! Adjust contributions to shearing strain on open boundaries.

1504 if (apply_obc) then ; if ((obc%strain_config == obc_strain_zero) .or. &

1505 (obc%strain_config == obc_strain_freeslip)) then

1506 do n=1,obc%number_of_segments

1507 j = obc%segment(n)%HI%JsdB ; i = obc%segment(n)%HI%IsdB

1508 if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= jeq)) then

1509 do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB

1510 if (obc%strain_config == obc_strain_zero) then

1511 ddel2vdx(i,j) = 0. ; ddel2udy(i,j) = 0.

1512 elseif (obc%strain_config == obc_strain_freeslip) then

1513 ddel2udy(i,j) = 0.

1514 endif

1515 enddo

1516 elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ieq)) then

1517 do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB

1518 if (obc%strain_config == obc_strain_zero) then

1519 ddel2vdx(i,j) = 0. ; ddel2udy(i,j) = 0.

1520 elseif (obc%strain_config == obc_strain_freeslip) then

1521 ddel2vdx(i,j) = 0.

1522 endif

1523 enddo

1524 endif

1525 enddo

1526 endif ; endif

1527 endif

1528

1529 if ((cs%Smagorinsky_Kh) .or. (cs%Smagorinsky_Ah)) then

1530 do j=js-1,jeq ; do i=is-1,ieq

1531 sh_xy_sq = sh_xy(i,j)**2

1532 sh_xx_sq = 0.25 * ( ((sh_xx(i,j)**2) + (sh_xx(i+1,j+1)**2)) &

1533 + ((sh_xx(i,j+1)**2) + (sh_xx(i+1,j)**2)) )

1534 shear_mag(i,j) = sqrt(sh_xy_sq + sh_xx_sq)

1535 enddo ; enddo

1536 endif

1537

1538 do j=js-1,jeq ; do i=is-1,ieq

1539 h2uq = 4.0 * (h_u(i,j) * h_u(i,j+1))

1540 h2vq = 4.0 * (h_v(i,j) * h_v(i+1,j))

1541 hq(i,j) = (2.0 * (h2uq * h2vq)) &

1542 / (h_neglect3 + (h2uq + h2vq) * ((h_u(i,j) + h_u(i,j+1)) + (h_v(i,j) + h_v(i+1,j))))

1543 enddo ; enddo

1544

1545 if (cs%bound_Ah .or. cs%bound_Kh) then

1546 do j=js-1,jeq ; do i=is-1,ieq

1547 h_min = min(h_u(i,j), h_u(i,j+1), h_v(i,j), h_v(i+1,j))

1548 hrat_min(i,j) = min(1.0, h_min / (hq(i,j) + h_neglect))

1549 enddo ; enddo

1550

1551 endif

1552

1553 if (cs%no_slip) then

1554 do j=js-1,jeq ; do i=is-1,ieq

1555 if (cs%no_slip .and. (g%mask2dBu(i,j) < 0.5)) then

1556 if ((g%mask2dCu(i,j) + g%mask2dCu(i,j+1)) + &

1557 (g%mask2dCv(i,j) + g%mask2dCv(i+1,j)) > 0.0) then

1558 ! This is a coastal vorticity point, so modify hq and hrat_min.

1559

1560 hu = g%mask2dCu(i,j) * h_u(i,j) + g%mask2dCu(i,j+1) * h_u(i,j+1)

1561 hv = g%mask2dCv(i,j) * h_v(i,j) + g%mask2dCv(i+1,j) * h_v(i+1,j)

1562 if ((g%mask2dCu(i,j) + g%mask2dCu(i,j+1)) * &

1563 (g%mask2dCv(i,j) + g%mask2dCv(i+1,j)) == 0.0) then

1564 ! Only one of hu and hv is nonzero, so just add them.

1565 hq(i,j) = hu + hv

1566 hrat_min(i,j) = 1.0

1567 else

1568 ! Both hu and hv are nonzero, so take the harmonic mean.

1569 hq(i,j) = 2.0 * (hu * hv) / ((hu + hv) + h_neglect)

1570 hrat_min(i,j) = min(1.0, min(hu, hv) / (hq(i,j) + h_neglect) )

1571 endif

1572 endif

1573 endif

1574 enddo ; enddo

1575 endif

1576

1577 ! Pass the velocity gradients and thickness to ZB2020

1578 if (cs%use_ZB2020) then

1579 call zb2020_copy_gradient_and_thickness(sh_xx, sh_xy, vort_xy, hq, g, gv, cs%ZB2020, k)

1580 endif

1581

1582 if (cs%Laplacian) then

1583 ! Determine the Laplacian viscosity at q points, using the

1584 ! largest value from several parameterizations. Also get the

1585 ! Laplacian component of str_xy.

1586

1587 if ((cs%Leith_Kh) .or. (cs%Leith_Ah)) then

1588 if (cs%use_QG_Leith_visc) then

1589 do j=js-1,jeq ; do i=is-1,ieq

1590 grad_vort = grad_vort_mag_q(i,j) + grad_div_mag_q(i,j)

1591 grad_vort_qg = 3. * grad_vort_mag_q_2d(i,j)

1592 vert_vort_mag(i,j) = min(grad_vort, grad_vort_qg)

1593 enddo ; enddo

1594 else

1595 do j=js-1,jeq ; do i=is-1,ieq

1596 vert_vort_mag(i,j) = grad_vort_mag_q(i,j) + grad_div_mag_q(i,j)

1597 enddo ; enddo

1598 endif

1599 endif

1600

1601 ! Static (pre-computed) background viscosity

1602 do j=js-1,jeq ; do i=is-1,ieq

1603 kh(i,j) = cs%Kh_bg_xy(i,j)

1604 enddo ; enddo

1605

1606 if (cs%Smagorinsky_Kh) then

1607 if (cs%add_LES_viscosity) then

1608 do j=js-1,jeq ; do i=is-1,ieq

1609 kh(i,j) = kh(i,j) + cs%Laplac2_const_xy(i,j) * shear_mag(i,j)

1610 enddo ; enddo

1611 else

1612 do j=js-1,jeq ; do i=is-1,ieq

1613 kh(i,j) = max(kh(i,j), cs%Laplac2_const_xy(i,j) * shear_mag(i,j) )

1614 enddo ; enddo

1615 endif

1616 endif

1617

1618 if (cs%Leith_Kh) then

1619 if (cs%add_LES_viscosity) then

1620 do j=js-1,jeq ; do i=is-1,ieq

1621 kh(i,j) = kh(i,j) + cs%Laplac3_const_xy(i,j) * vert_vort_mag(i,j) * inv_pi3 ! Is this right? -AJA

1622 enddo ; enddo

1623 else

1624 do j=js-1,jeq ; do i=is-1,ieq

1625 kh(i,j) = max(kh(i,j), cs%Laplac3_const_xy(i,j) * vert_vort_mag(i,j) * inv_pi3)

1626 enddo ; enddo

1627 endif

1628 endif

1629

1630 ! All viscosity contributions above are subject to resolution scaling

1631

1632 if (rescale_kh) then

1633 do j=js-1,jeq ; do i=is-1,ieq

1634 kh(i,j) = varmix%Res_fn_q(i,j) * kh(i,j)

1635 enddo ; enddo

1636 endif

1637

1638 do j=js-1,jeq ; do i=is-1,ieq

1639 kh(i,j) = max(kh(i,j), cs%Kh_bg_min) ! Place a floor on the viscosity, if desired.

1640 enddo ; enddo

1641

1642 if (use_meke_ku .and. .not. cs%EY24_EBT_BS) then

1643 if (use_kh_struct) then

1644 do j=js-1,jeq ; do i=is-1,ieq

1645 meke_res_fn = 1.

1646 if (cs%res_scale_MEKE) meke_res_fn = varmix%Res_fn_q(i,j)

1647

1648 kh(i,j) = kh(i,j) + 0.25*( ((meke%Ku(i,j)*varmix%BS_struct(i,j,k)) + &

1649 (meke%Ku(i+1,j+1)*varmix%BS_struct(i+1,j+1,k))) + &

1650 ((meke%Ku(i+1,j)*varmix%BS_struct(i+1,j,k)) + &

1651 (meke%Ku(i,j+1)*varmix%BS_struct(i,j+1,k))) ) * meke_res_fn

1652 enddo ; enddo

1653 else

1654 do j=js-1,jeq ; do i=is-1,ieq

1655 meke_res_fn = 1.

1656 if (cs%res_scale_MEKE) meke_res_fn = varmix%Res_fn_q(i,j)

1657

1658 kh(i,j) = kh(i,j) + 0.25 * ( &

1659 (meke%Ku(i,j) + meke%Ku(i+1,j+1)) + &

1660 (meke%Ku(i+1,j) + &

1661 meke%Ku(i,j+1)) ) * meke_res_fn

1662 enddo ; enddo

1663 endif

1664 endif

1665

1666 if (cs%anisotropic) then

1667 ! *Add* the shear component of anisotropic viscosity

1668 do j=js-1,jeq ; do i=is-1,ieq

1669 kh(i,j) = kh(i,j) + cs%Kh_aniso * cs%n1n2_q(i,j)**2

1670 enddo ; enddo

1671 endif

1672

1673 do j=js-1,jeq ; do i=is-1,ieq

1674 ! Newer method of bounding for stability

1675 if ((cs%bound_Kh) .and. (cs%bound_Ah)) then

1676 visc_bound_rem(i,j) = 1.0

1677 kh_max_here = hrat_min(i,j) * cs%Kh_Max_xy(i,j)

1678 if (kh(i,j) >= kh_max_here) then

1679 visc_bound_rem(i,j) = 0.0

1680 kh(i,j) = kh_max_here

1681 elseif ((kh(i,j) > 0.0) .or. (cs%backscatter_underbound .and. (kh_max_here > 0.0))) then

1682 visc_bound_rem(i,j) = 1.0 - kh(i,j) / kh_max_here

1683 endif

1684 elseif (cs%bound_Kh) then

1685 kh(i,j) = min(kh(i,j), hrat_min(i,j) * cs%Kh_Max_xy(i,j))

1686 endif

1687 enddo ; enddo

1688

1689 if (cs%use_Leithy) then

1690 ! Leith+E doesn't recompute Kh at q points, it just interpolates it from h to q points

1691 do j=js-1,jeq ; do i=is-1,ieq

1692 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)))

1693 enddo ; enddo

1694 endif

1695

1696 if (cs%id_Kh_q > 0 .or. cs%debug) then

1697 do j=js-1,jeq ; do i=is-1,ieq

1698 kh_q(i,j,k) = kh(i,j)

1699 enddo ; enddo

1700 endif

1701

1702 if (cs%id_vort_xy_q > 0) then

1703 do j=js-1,jeq ; do i=is-1,ieq

1704 vort_xy_q(i,j,k) = vort_xy(i,j)

1705 enddo ; enddo

1706 endif

1707

1708 if (cs%id_sh_xy_q > 0) then

1709 do j=js-1,jeq ; do i=is-1,ieq

1710 sh_xy_q(i,j,k) = sh_xy(i,j)

1711 enddo ; enddo

1712 endif

1713

1714 if (.not. cs%use_Leithy) then

1715 do j=js-1,jeq ; do i=is-1,ieq

1716 str_xy(i,j) = -kh(i,j) * sh_xy(i,j)

1717 enddo ; enddo

1718 else

1719 do j=js-1,jeq ; do i=is-1,ieq

1720 str_xy(i,j) = -kh(i,j) * sh_xy_smooth(i,j)

1721 enddo ; enddo

1722 endif

1723 else

1724 do j=js-1,jeq ; do i=is-1,ieq

1725 str_xy(i,j) = 0.

1726 enddo ; enddo

1727 endif ! get harmonic coefficient Kh at q points and harmonic part of str_xy

1728

1729 if (cs%anisotropic) then

1730 do j=js-1,jeq ; do i=is-1,ieq

1731 ! Horizontal-tension averaged to q-points

1732 local_strain = 0.25 * ( (sh_xx(i,j) + sh_xx(i+1,j+1)) + (sh_xx(i+1,j) + sh_xx(i,j+1)) )

1733 ! *Add* the tension contribution to the xy-component of stress

1734 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

1735 enddo ; enddo

1736 endif

1737

1738 if (cs%biharmonic) then

1739 ! Determine the biharmonic viscosity at q points, using the

1740 ! largest value from several parameterizations. Also get the

1741 ! biharmonic component of str_xy.

1742 do j=js-1,jeq ; do i=is-1,ieq

1743 ah(i,j) = cs%Ah_bg_xy(i,j)

1744 enddo ; enddo

1745

1746 if (cs%Smagorinsky_Ah .or. cs%Leith_Ah) then

1747 if (cs%Smagorinsky_Ah) then

1748 if (cs%bound_Coriolis) then

1749 do j=js-1,jeq ; do i=is-1,ieq

1750 ahsm = shear_mag(i,j) * (cs%Biharm_const_xy(i,j) &

1751 + cs%Biharm_const2_xy(i,j) * shear_mag(i,j))

1752 ah(i,j) = max(ah(i,j), ahsm)

1753 enddo ; enddo

1754 else

1755 do j=js-1,jeq ; do i=is-1,ieq

1756 ahsm = cs%Biharm_const_xy(i,j) * shear_mag(i,j)

1757 ah(i,j) = max(ah(i,j), ahsm)

1758 enddo ; enddo

1759 endif

1760 endif

1761

1762 if (cs%Leith_Ah) then

1763 do j=js-1,jeq ; do i=is-1,ieq

1764 ahlth = cs%Biharm6_const_xy(i,j) * abs(del2vort_q(i,j)) * inv_pi6

1765 ah(i,j) = max(ah(i,j), ahlth)

1766 enddo ; enddo

1767 endif

1768 endif ! Smagorinsky_Ah or Leith_Ah

1769

1770 if (use_meke_au) then

1771 ! *Add* the MEKE contribution

1772 do j=js-1,jeq ; do i=is-1,ieq

1773 ah(i,j) = ah(i,j) + 0.25 * ( &

1774 (meke%Au(i,j) + meke%Au(i+1,j+1)) + (meke%Au(i+1,j) + meke%Au(i,j+1)) )

1775 enddo ; enddo

1776 endif

1777

1778 if (cs%Re_Ah > 0.0) then

1779 do j=js-1,jeq ; do i=is-1,ieq

1780 ke = 0.125 * (((u(i,j,k) + u(i,j+1,k))**2) + ((v(i,j,k) + v(i+1,j,k))**2))

1781 ah(i,j) = sqrt(ke) * cs%Re_Ah_const_xy(i,j)

1782 enddo ; enddo

1783 endif

1784

1785 if (cs%bound_Ah) then

1786 if (cs%bound_Kh) then

1787 do j=js-1,jeq ; do i=is-1,ieq

1788 ah(i,j) = min(ah(i,j), visc_bound_rem(i,j) * hrat_min(i,j) * cs%Ah_Max_xy(i,j))

1789 enddo ; enddo

1790 else

1791 do j=js-1,jeq ; do i=is-1,ieq

1792 ah(i,j) = min(ah(i,j), hrat_min(i,j) * cs%Ah_Max_xy(i,j))

1793 enddo ; enddo

1794 endif

1795 endif

1796

1797 if (cs%EY24_EBT_BS) then

1798 do j=js-1,jeq ; do i=is-1,ieq

1799 tmp = cs%KS_coef *hrat_min(i,j) * cs%Ah_Max_xy_KS(i,j)

1800 visc_limit_q(i,j,k) = tmp

1801 visc_limit_q_frac(i,j,k) = ah(i,j) / (cs%KS_coef * hrat_min(i,j) * cs%Ah_Max_xy_KS(i,j))

1802 if (ah(i,j) >= tmp) then

1803 visc_limit_q_flag(i,j,k) = 1.

1804 endif

1805 enddo ; enddo

1806 endif

1807

1808 ! Leith+E doesn't recompute Ah at q points, it just interpolates it from h to q points

1809 if (cs%use_Leithy) then

1810 do j=js-1,jeq ; do i=is-1,ieq

1811 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)))

1812 enddo ; enddo

1813 endif

1814

1815 if (cs%id_Ah_q>0 .or. cs%debug) then

1816 do j=js-1,jeq ; do i=is-1,ieq

1817 ah_q(i,j,k) = ah(i,j)

1818 enddo ; enddo

1819 endif

1820

1821 ! Again, need to initialize str_xy as if its biharmonic

1822 do j=js-1,jeq ; do i=is-1,ieq

1823 d_str = ah(i,j) * (ddel2vdx(i,j) + ddel2udy(i,j))

1824

1825 str_xy(i,j) = str_xy(i,j) + d_str

1826

1827 ! Keep a copy of the biharmonic contribution for backscatter parameterization

1828 bhstr_xy(i,j) = d_str * (hq(i,j) * g%mask2dBu(i,j) * cs%reduction_xy(i,j))

1829 enddo ; enddo

1830 endif ! Get Ah at q points and biharmonic part of str_xy

1831

1832 ! Backscatter using MEKE

1833 if (cs%EY24_EBT_BS) then

1834 do j=js-1,jeq ; do i=is-1,ieq

1835 if (visc_limit_q_flag(i,j,k) > 0) then

1836 kh_bs(i,j) = 0.

1837 else

1838 if (use_kh_struct) then

1839 kh_bs(i,j) = 0.25*( ((meke%Ku(i,j)*varmix%BS_struct(i,j,k)) + &

1840 (meke%Ku(i+1,j+1)*varmix%BS_struct(i+1,j+1,k))) + &

1841 ((meke%Ku(i+1,j)*varmix%BS_struct(i+1,j,k)) + &

1842 (meke%Ku(i,j+1)*varmix%BS_struct(i,j+1,k))) )

1843 else

1844 kh_bs(i,j) = 0.25*( (meke%Ku(i,j) + meke%Ku(i+1,j+1)) + &

1845 (meke%Ku(i+1,j) + meke%Ku(i,j+1)) )

1846 endif

1847 endif

1848 enddo ; enddo

1849

1850 do j=js-1,jeq ; do i=is-1,ieq

1851 str_xy_bs(i,j) = -kh_bs(i,j) * (sh_xy(i,j))

1852 enddo ; enddo

1853

1854 if (cs%id_BS_coeff_q>0) then

1855 do j=js-1,jeq ; do i=is-1,ieq

1856 bs_coeff_q(i,j,k) = kh_bs(i,j)

1857 enddo ; enddo

1858 endif

1859

1860 do j=js-1,jeq ; do i=is-1,ieq

1861 str_xy(i,j) = str_xy(i,j) + str_xy_bs(i,j)

1862 enddo ; enddo

1863 endif ! Backscatter

1864

1865 if (cs%use_GME) then

1866 ! The wider halo here is to permit one pass of smoothing without a halo update.

1867 do j=jsq-1,jeq+2 ; do i=isq-1,ieq+2

1868 gme_coeff = gme_effic_h(i,j) * 0.25 * &

1869 ((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)))

1870 gme_coeff = min(gme_coeff, cs%GME_limiter)

1871

1872 if ((cs%id_GME_coeff_h>0) .or. find_frictwork) gme_coeff_h(i,j,k) = gme_coeff

1873 str_xx_gme(i,j) = gme_coeff * sh_xx_bt(i,j)

1874 enddo ; enddo

1875

1876 ! The wider halo here is to permit one pass of smoothing without a halo update.

1877 do j=js-2,je+1 ; do i=is-2,ie+1

1878 gme_coeff = gme_effic_q(i,j) * 0.25 * &

1879 ((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)))

1880 gme_coeff = min(gme_coeff, cs%GME_limiter)

1881

1882 if (cs%id_GME_coeff_q>0) gme_coeff_q(i,j,k) = gme_coeff

1883 str_xy_gme(i,j) = gme_coeff * sh_xy_bt(i,j)

1884 enddo ; enddo

1885

1886 ! Applying GME diagonal term. This is linear and the arguments can be rescaled.

1887 call smooth_gme(cs, g, gme_flux_h=str_xx_gme)

1888 call smooth_gme(cs, g, gme_flux_q=str_xy_gme)

1889

1890 ! 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].

1891 do j=jsq,jeq+1 ; do i=isq,ieq+1

1892 str_xx(i,j) = (str_xx(i,j) + str_xx_gme(i,j)) * (h(i,j,k) * cs%reduction_xx(i,j))

1893 enddo ; enddo

1894

1895 ! 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].

1896 if (cs%no_slip) then

1897 do j=js-1,jeq ; do i=is-1,ieq

1898 str_xy(i,j) = (str_xy(i,j) + str_xy_gme(i,j)) * (hq(i,j) * cs%reduction_xy(i,j))

1899 enddo ; enddo

1900 else

1901 do j=js-1,jeq ; do i=is-1,ieq

1902 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))

1903 enddo ; enddo

1904 endif

1905

1906 else ! .not. use_GME

1907 ! 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].

1908 do j=jsq,jeq+1 ; do i=isq,ieq+1

1909 str_xx(i,j) = str_xx(i,j) * (h(i,j,k) * cs%reduction_xx(i,j))

1910 enddo ; enddo

1911

1912 ! 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].

1913 if (cs%no_slip) then

1914 do j=js-1,jeq ; do i=is-1,ieq

1915 str_xy(i,j) = str_xy(i,j) * (hq(i,j) * cs%reduction_xy(i,j))

1916 enddo ; enddo

1917 else

1918 do j=js-1,jeq ; do i=is-1,ieq

1919 str_xy(i,j) = str_xy(i,j) * (hq(i,j) * g%mask2dBu(i,j) * cs%reduction_xy(i,j))

1920 enddo ; enddo

1921 endif

1922 endif ! use_GME

1923

1924 ! Evaluate 1/h x.Div(h Grad u) or the biharmonic equivalent.

1925 do j=js,je ; do i=isq,ieq

1926 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))) + &

1927 g%IdyCu(i,j)*((cs%dy2h(i,j)*str_xx(i,j)) - (cs%dy2h(i+1,j)*str_xx(i+1,j)))) * &

1928 g%IareaCu(i,j)) / (h_u(i,j) + h_neglect)

1929 enddo ; enddo

1930

1931 if (apply_obc) then

1932 ! This is not the right boundary condition. If all the masking of tendencies are done

1933 ! correctly later then eliminating this block should not change answers.

1934 do n=1,obc%number_of_segments

1935 if (obc%segment(n)%is_E_or_W) then

1936 i = obc%segment(n)%HI%IsdB

1937 do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed

1938 diffu(i,j,k) = 0.

1939 enddo

1940 endif

1941 enddo

1942 endif

1943

1944 ! Evaluate 1/h y.Div(h Grad u) or the biharmonic equivalent.

1945 do j=jsq,jeq ; do i=is,ie

1946 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))) - &

1947 g%IdxCv(i,j)*((cs%dx2h(i,j)*str_xx(i,j)) - (cs%dx2h(i,j+1)*str_xx(i,j+1)))) * &

1948 g%IareaCv(i,j)) / (h_v(i,j) + h_neglect)

1949 enddo ; enddo

1950

1951 if (apply_obc) then

1952 ! This is not the right boundary condition. If all the masking of tendencies are done

1953 ! correctly later then eliminating this block should not change answers.

1954 do n=1,obc%number_of_segments

1955 if (obc%segment(n)%is_N_or_S) then

1956 j = obc%segment(n)%HI%JsdB

1957 do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied

1958 diffv(i,j,k) = 0.

1959 enddo

1960 endif

1961 enddo

1962 endif

1963

1964 if (find_frictwork) then

1965 if (cs%FrictWork_bug) then

1966 ! Diagnose str_xx*d_x u - str_yy*d_y v + str_xy*(d_y u + d_x v)

1967 ! This is the old formulation that includes energy diffusion

1968 do j=js,je ; do i=is,ie

1969 frictwork(i,j,k) = gv%H_to_RZ * ( &

1970 ((str_xx(i,j) * (u(i,j,k)-u(i-1,j,k))*g%IdxT(i,j)) &

1971 - (str_xx(i,j) * (v(i,j,k)-v(i,j-1,k))*g%IdyT(i,j))) &

1972 + 0.25*(( (str_xy(i,j) * &

1973 (((u(i,j+1,k)-u(i,j,k))*g%IdyBu(i,j)) &

1974 + ((v(i+1,j,k)-v(i,j,k))*g%IdxBu(i,j)))) &

1975 + (str_xy(i-1,j-1) * &

1976 (((u(i-1,j,k)-u(i-1,j-1,k))*g%IdyBu(i-1,j-1)) &

1977 + ((v(i,j-1,k)-v(i-1,j-1,k))*g%IdxBu(i-1,j-1)))) ) &

1978 + ( (str_xy(i-1,j) * &

1979 (((u(i-1,j+1,k)-u(i-1,j,k))*g%IdyBu(i-1,j)) &

1980 + ((v(i,j,k)-v(i-1,j,k))*g%IdxBu(i-1,j)))) &

1981 + (str_xy(i,j-1) * &

1982 (((u(i,j,k)-u(i,j-1,k))*g%IdyBu(i,j-1)) &

1983 + ((v(i+1,j-1,k)-v(i,j-1,k))*g%IdxBu(i,j-1)))) ) ) )

1984 enddo ; enddo

1985 else

1986 do j=js,je ; do i=is,ie

1987 frictwork(i,j,k) = gv%H_to_RZ * g%IareaT(i,j) * ( &

1988 ((str_xx(i,j)*cs%dy2h(i,j) * ( &

1989 (uh(i,j,k)*g%dxCu(i,j)*g%IdyCu(i,j)*g%IareaCu(i,j)/(h_u(i,j)+h_neglect)) &

1990 - (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)) ) ) &

1991 - (str_xx(i,j)*cs%dx2h(i,j) * ( &

1992 (vh(i,j,k)*g%dyCv(i,j)*g%IdxCv(i,j)*g%IareaCv(i,j)/(h_v(i,j)+h_neglect)) &

1993 - (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)) ) )) &

1994 + (0.25*(((str_xy(i,j)*( &

1995 (cs%dx2q(i,j)*((uh(i,j+1,k)*g%IareaCu(i,j+1)/(h_u(i,j+1)+h_neglect)) &

1996 - (uh(i,j,k)*g%IareaCu(i,j)/(h_u(i,j)+h_neglect)))) &

1997 + (cs%dy2q(i,j)*((vh(i+1,j,k)*g%IareaCv(i+1,j)/(h_v(i+1,j)+h_neglect)) &

1998 - (vh(i,j,k)*g%IareaCv(i,j)/(h_v(i,j)+h_neglect)))) )) &

1999 +(str_xy(i-1,j-1)*( &

2000 (cs%dx2q(i-1,j-1)*((uh(i-1,j,k)*g%IareaCu(i-1,j)/(h_u(i-1,j)+h_neglect)) &

2001 - (uh(i-1,j-1,k)*g%IareaCu(i-1,j-1)/(h_u(i-1,j-1)+h_neglect)))) &

2002 + (cs%dy2q(i-1,j-1)*((vh(i,j-1,k)*g%IareaCv(i,j-1)/(h_v(i,j-1)+h_neglect)) &

2003 - (vh(i-1,j-1,k)*g%IareaCv(i-1,j-1)/(h_v(i-1,j-1)+h_neglect)))) )) ) &

2004 +((str_xy(i-1,j)*( &

2005 (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)) &

2006 - (uh(i-1,j,k)*g%IareaCu(i-1,j)/(h_u(i-1,j)+h_neglect)))) &

2007 + (cs%dy2q(i-1,j)*((vh(i,j,k)*g%IareaCv(i,j)/(h_v(i,j)+h_neglect)) &

2008 - (vh(i-1,j,k)*g%IareaCv(i-1,j)/(h_v(i-1,j)+h_neglect)))) )) &

2009 +(str_xy(i,j-1)*( &

2010 (cs%dx2q(i,j-1)*((uh(i,j,k)*g%IareaCu(i,j)/(h_u(i,j)+h_neglect)) &

2011 - (uh(i,j-1,k)*g%IareaCu(i,j-1)/(h_u(i,j-1)+h_neglect)))) &

2012 + (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)) &

2013 - (vh(i,j-1,k)*g%IareaCv(i,j-1)/(h_v(i,j-1)+h_neglect)))) )) ) )) )

2014

2015 enddo ; enddo

2016 endif

2017

2018 if (cs%EY24_EBT_BS) then

2019 do j=js,je ; do i=is,ie

2020 frictwork(i,j,k) = (1. - visc_limit_h_flag(i,j,k)) * frictwork(i,j,k)

2021 enddo ; enddo

2022 endif

2023 endif

2024

2025 if (cs%id_FrictWork_bh>0 .or. cs%id_FrictWorkIntz_bh > 0 .or. allocated(meke%mom_src_bh)) then

2026 if (cs%FrictWork_bug) then

2027 ! Diagnose bhstr_xx*d_x u - bhstr_yy*d_y v + bhstr_xy*(d_y u + d_x v)

2028 ! This is the old formulation that includes energy diffusion !cyc

2029 do j=js,je ; do i=is,ie

2030 frictwork_bh(i,j,k) = gv%H_to_RZ * ( &

2031 ((bhstr_xx(i,j) * (u(i,j,k)-u(i-1,j,k))*g%IdxT(i,j)) &

2032 - (bhstr_xx(i,j) * (v(i,j,k)-v(i,j-1,k))*g%IdyT(i,j))) &

2033 + 0.25*(( (bhstr_xy(i,j) * &

2034 (((u(i,j+1,k)-u(i,j,k))*g%IdyBu(i,j)) &

2035 + ((v(i+1,j,k)-v(i,j,k))*g%IdxBu(i,j)))) &

2036 + (bhstr_xy(i-1,j-1) * &

2037 (((u(i-1,j,k)-u(i-1,j-1,k))*g%IdyBu(i-1,j-1)) &

2038 + ((v(i,j-1,k)-v(i-1,j-1,k))*g%IdxBu(i-1,j-1)))) ) &

2039 + ( (bhstr_xy(i-1,j) * &

2040 (((u(i-1,j+1,k)-u(i-1,j,k))*g%IdyBu(i-1,j)) &

2041 + ((v(i,j,k)-v(i-1,j,k))*g%IdxBu(i-1,j)))) &

2042 + (bhstr_xy(i,j-1) * &

2043 (((u(i,j,k)-u(i,j-1,k))*g%IdyBu(i,j-1)) &

2044 + ((v(i+1,j-1,k)-v(i,j-1,k))*g%IdxBu(i,j-1)))) ) ) )

2045 enddo ; enddo

2046 else

2047 do j=js,je ; do i=is,ie

2048 ! Diagnose bhstr_xx*d_x u - bhstr_yy*d_y v + bhstr_xy*(d_y u + d_x v)

2049 frictwork_bh(i,j,k) = gv%H_to_RZ * g%IareaT(i,j) * ( &

2050 ((bhstr_xx(i,j)*cs%dy2h(i,j) * ( &

2051 (uh(i,j,k)*g%dxCu(i,j)*g%IdyCu(i,j)*g%IareaCu(i,j)/(h_u(i,j)+h_neglect)) &

2052 - (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)) ) ) &

2053 - (bhstr_xx(i,j)*cs%dx2h(i,j) * ( &

2054 (vh(i,j,k)*g%dyCv(i,j)*g%IdxCv(i,j)*g%IareaCv(i,j)/(h_v(i,j)+h_neglect)) &

2055 - (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)) ) )) &

2056 + (0.25*(((bhstr_xy(i,j)*( &

2057 (cs%dx2q(i,j)*((uh(i,j+1,k)*g%IareaCu(i,j+1)/(h_u(i,j+1)+h_neglect)) &

2058 - (uh(i,j,k)*g%IareaCu(i,j)/(h_u(i,j)+h_neglect)))) &

2059 + (cs%dy2q(i,j)*((vh(i+1,j,k)*g%IareaCv(i+1,j)/(h_v(i+1,j)+h_neglect)) &

2060 - (vh(i,j,k)*g%IareaCv(i,j)/(h_v(i,j)+h_neglect)))) )) &

2061 +(bhstr_xy(i-1,j-1)*( &

2062 (cs%dx2q(i-1,j-1)*((uh(i-1,j,k)*g%IareaCu(i-1,j)/(h_u(i-1,j)+h_neglect)) &

2063 - (uh(i-1,j-1,k)*g%IareaCu(i-1,j-1)/(h_u(i-1,j-1)+h_neglect)))) &

2064 + (cs%dy2q(i-1,j-1)*((vh(i,j-1,k)*g%IareaCv(i,j-1)/(h_v(i,j-1)+h_neglect)) &

2065 - (vh(i-1,j-1,k)*g%IareaCv(i-1,j-1)/(h_v(i-1,j-1)+h_neglect)))) )) ) &

2066 +((bhstr_xy(i-1,j)*( &

2067 (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)) &

2068 - (uh(i-1,j,k)*g%IareaCu(i-1,j)/(h_u(i-1,j)+h_neglect)))) &

2069 + (cs%dy2q(i-1,j)*((vh(i,j,k)*g%IareaCv(i,j)/(h_v(i,j)+h_neglect)) &

2070 - (vh(i-1,j,k)*g%IareaCv(i-1,j)/(h_v(i-1,j)+h_neglect)))) )) &

2071 +(bhstr_xy(i,j-1)*( &

2072 (cs%dx2q(i,j-1)*((uh(i,j,k)*g%IareaCu(i,j)/(h_u(i,j)+h_neglect)) &

2073 - (uh(i,j-1,k)*g%IareaCu(i,j-1)/(h_u(i,j-1)+h_neglect)))) &

2074 + (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)) &

2075 - (vh(i,j-1,k)*g%IareaCv(i,j-1)/(h_v(i,j-1)+h_neglect)))) )) ) )) )

2076 enddo ; enddo

2077 endif

2078

2079 if (cs%EY24_EBT_BS) then

2080 do j=js,je ; do i=is,ie

2081 frictwork_bh(i,j,k) = (1. - visc_limit_h_flag(i,j,k)) * frictwork_bh(i,j,k)

2082 enddo ; enddo

2083 endif

2084 endif

2085

2086 if (cs%use_GME) then

2087 if (cs%FrictWork_bug) then ; do j=js,je ; do i=is,ie

2088 ! Diagnose str_xx_GME*d_x u - str_yy_GME*d_y v + str_xy_GME*(d_y u + d_x v)

2089 ! This is the old formulation that includes energy diffusion

2090 frictwork_gme(i,j,k) = gv%H_to_RZ * ( &

2091 ((str_xx_gme(i,j)*(u(i,j,k)-u(i-1,j,k))*g%IdxT(i,j)) &

2092 - (str_xx_gme(i,j)*(v(i,j,k)-v(i,j-1,k))*g%IdyT(i,j))) &

2093 + 0.25*(( (str_xy_gme(i,j) * &

2094 (((u(i,j+1,k)-u(i,j,k))*g%IdyBu(i,j)) &

2095 + ((v(i+1,j,k)-v(i,j,k))*g%IdxBu(i,j)))) &

2096 + (str_xy_gme(i-1,j-1) * &

2097 (((u(i-1,j,k)-u(i-1,j-1,k))*g%IdyBu(i-1,j-1)) &

2098 + ((v(i,j-1,k)-v(i-1,j-1,k))*g%IdxBu(i-1,j-1)))) ) &

2099 + ( (str_xy_gme(i-1,j) * &

2100 (((u(i-1,j+1,k)-u(i-1,j,k))*g%IdyBu(i-1,j)) &

2101 + ((v(i,j,k)-v(i-1,j,k))*g%IdxBu(i-1,j)))) &

2102 + (str_xy_gme(i,j-1) * &

2103 (((u(i,j,k)-u(i,j-1,k))*g%IdyBu(i,j-1)) &

2104 + ((v(i+1,j-1,k)-v(i,j-1,k))*g%IdxBu(i,j-1)))) ) ) )

2105 enddo ; enddo

2106 else ; do j=js,je ; do i=is,ie

2107 frictwork_gme(i,j,k) = gv%H_to_RZ * g%IareaT(i,j) * ( &

2108 ((str_xx_gme(i,j)*cs%dy2h(i,j) * ( &

2109 (uh(i,j,k)*g%dxCu(i,j)*g%IdyCu(i,j)*g%IareaCu(i,j)/(h_u(i,j)+h_neglect)) &

2110 - (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)) ) ) &

2111 - (str_xx_gme(i,j)*cs%dx2h(i,j) * ( &

2112 (vh(i,j,k)*g%dyCv(i,j)*g%IdxCv(i,j)*g%IareaCv(i,j)/(h_v(i,j)+h_neglect)) &

2113 - (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)) ) )) &

2114 + (0.25*(((str_xy_gme(i,j)*( &

2115 (cs%dx2q(i,j)*((uh(i,j+1,k)*g%IareaCu(i,j+1)/(h_u(i,j+1)+h_neglect)) &

2116 - (uh(i,j,k)*g%IareaCu(i,j)/(h_u(i,j)+h_neglect)))) &

2117 + (cs%dy2q(i,j)*((vh(i+1,j,k)*g%IareaCv(i+1,j)/(h_v(i+1,j)+h_neglect)) &

2118 - (vh(i,j,k)*g%IareaCv(i,j)/(h_v(i,j)+h_neglect)))) )) &

2119 +(str_xy_gme(i-1,j-1)*( &

2120 (cs%dx2q(i-1,j-1)*((uh(i-1,j,k)*g%IareaCu(i-1,j)/(h_u(i-1,j)+h_neglect)) &

2121 - (uh(i-1,j-1,k)*g%IareaCu(i-1,j-1)/(h_u(i-1,j-1)+h_neglect)))) &

2122 + (cs%dy2q(i-1,j-1)*((vh(i,j-1,k)*g%IareaCv(i,j-1)/(h_v(i,j-1)+h_neglect)) &

2123 - (vh(i-1,j-1,k)*g%IareaCv(i-1,j-1)/(h_v(i-1,j-1)+h_neglect)))) )) ) &

2124 +((str_xy_gme(i-1,j)*( &

2125 (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)) &

2126 - (uh(i-1,j,k)*g%IareaCu(i-1,j)/(h_u(i-1,j)+h_neglect)))) &

2127 + (cs%dy2q(i-1,j)*((vh(i,j,k)*g%IareaCv(i,j)/(h_v(i,j)+h_neglect)) &

2128 - (vh(i-1,j,k)*g%IareaCv(i-1,j)/(h_v(i-1,j)+h_neglect)))) )) &

2129 +(str_xy_gme(i,j-1)*( &

2130 (cs%dx2q(i,j-1)*((uh(i,j,k)*g%IareaCu(i,j)/(h_u(i,j)+h_neglect)) &

2131 - (uh(i,j-1,k)*g%IareaCu(i,j-1)/(h_u(i,j-1)+h_neglect)))) &

2132 + (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)) &

2133 - (vh(i,j-1,k)*g%IareaCv(i,j-1)/(h_v(i,j-1)+h_neglect)))) )) ) )) )

2134 enddo ; enddo ; endif

2135 endif

2136

2137 if (skeb_use_frict) then ; do j=js,je ; do i=is,ie

2138 ! Note that the sign convention is FrictWork < 0 means energy dissipation.

2139 stoch%skeb_diss(i,j,k) = stoch%skeb_diss(i,j,k) - stoch%skeb_frict_coef * &

2140 frictwork(i,j,k) / (gv%H_to_RZ * (h(i,j,k) + h_neglect))

2141 enddo ; enddo ; endif

2142

2143 ! Make a similar calculation as for FrictWork above but accumulating into

2144 ! the vertically integrated MEKE source term, and adjusting for any

2145 ! energy loss seen as a reduction in the (biharmonic) frictional source term.

2146 if (find_frictwork .and. allocated(meke%mom_src)) then

2147 if (k==1) then

2148 do j=js,je ; do i=is,ie

2149 meke%mom_src(i,j) = 0.

2150 enddo ; enddo

2151

2152 if (allocated(meke%mom_src_bh)) then

2153 do j=js,je ; do i=is,ie

2154 meke%mom_src_bh(i,j) = 0.

2155 enddo ; enddo

2156 endif

2157

2158 if (allocated(meke%GME_snk)) then

2159 do j=js,je ; do i=is,ie

2160 meke%GME_snk(i,j) = 0.

2161 enddo ; enddo

2162 endif

2163 endif

2164 if (meke%backscatter_Ro_c /= 0.) then

2165 do j=js,je ; do i=is,ie

2166 fath = 0.25*( (abs(g%CoriolisBu(i-1,j-1)) + abs(g%CoriolisBu(i,j))) + &

2167 (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j-1))) )

2168 shear_mag_bc = sqrt(sh_xx(i,j) * sh_xx(i,j) + &

2169 0.25*(((sh_xy(i-1,j-1)*sh_xy(i-1,j-1)) + (sh_xy(i,j)*sh_xy(i,j))) + &

2170 ((sh_xy(i-1,j)*sh_xy(i-1,j)) + (sh_xy(i,j-1)*sh_xy(i,j-1)))))

2171 if ((cs%answer_date > 20190101) .and. (cs%answer_date < 20241201)) then

2172 fath = (us%s_to_T*fath)**meke%backscatter_Ro_pow ! f^n

2173 ! Note the hard-coded dimensional constant in the following line that can not

2174 ! be rescaled for dimensional consistency.

2175 shear_mag_bc = (((us%s_to_T * shear_mag_bc)**meke%backscatter_Ro_pow) + 1.e-30) &

2176 * meke%backscatter_Ro_c ! c * D^n

2177 ! The Rossby number function is g(Ro) = 1/(1+c.Ro^n)

2178 ! RoScl = 1 - g(Ro)

2179 roscl = shear_mag_bc / (fath + shear_mag_bc) ! = 1 - f^n/(f^n+c*D^n)

2180 else

2181 if (fath <= backscat_subround*shear_mag_bc) then

2182 roscl = 1.0

2183 else

2184 sh_f_pow = meke%backscatter_Ro_c * (shear_mag_bc / fath)**meke%backscatter_Ro_pow

2185 roscl = sh_f_pow / (1.0 + sh_f_pow) ! = 1 - f^n/(f^n+c*D^n)

2186 endif

2187 endif

2188

2189 meke%mom_src(i,j) = meke%mom_src(i,j) + (frictwork(i,j,k) - roscl*frictwork_bh(i,j,k))

2190

2191 if (allocated(meke%mom_src_bh)) &

2192 meke%mom_src_bh(i,j) = meke%mom_src_bh(i,j) &

2193 + (frictwork_bh(i,j,k) - roscl * frictwork_bh(i,j,k))

2194 enddo ; enddo

2195 else

2196 do j=js,je ; do i=is,ie

2197 meke%mom_src(i,j) = meke%mom_src(i,j) + frictwork(i,j,k)

2198 enddo ; enddo

2199

2200 if (allocated(meke%mom_src_bh)) then

2201 do j=js,je ; do i=is,ie

2202 meke%mom_src_bh(i,j) = meke%mom_src_bh(i,j) + frictwork_bh(i,j,k)

2203 enddo ; enddo

2204 endif

2205 endif ! MEKE%backscatter_Ro_c

2206

2207 if (cs%use_GME .and. allocated(meke%GME_snk)) then

2208 do j=js,je ; do i=is,ie

2209 meke%GME_snk(i,j) = meke%GME_snk(i,j) + frictwork_gme(i,j,k)

2210 enddo ; enddo

2211 endif

2212 endif ! find_FrictWork and associated(mom_src)

2213 enddo ! end of k loop

2214

2215 ! Offer fields for diagnostic averaging.

2216 if (cs%id_normstress > 0) call post_data(cs%id_normstress, nost, cs%diag)

2217 if (cs%id_shearstress > 0) call post_data(cs%id_shearstress, shst, cs%diag)

2218 if (cs%id_diffu>0) call post_data(cs%id_diffu, diffu, cs%diag)

2219 if (cs%id_diffv>0) call post_data(cs%id_diffv, diffv, cs%diag)

2220 if (cs%id_FrictWork>0) call post_data(cs%id_FrictWork, frictwork, cs%diag)

2221 if (cs%id_FrictWork_bh>0) call post_data(cs%id_FrictWork_bh, frictwork_bh, cs%diag)

2222 if (cs%id_Ah_h>0) call post_data(cs%id_Ah_h, ah_h, cs%diag)

2223 if (cs%id_grid_Re_Ah>0) call post_data(cs%id_grid_Re_Ah, grid_re_ah, cs%diag)

2224 if (cs%id_div_xx_h>0) call post_data(cs%id_div_xx_h, div_xx_h, cs%diag)

2225 if (cs%id_vort_xy_q>0) call post_data(cs%id_vort_xy_q, vort_xy_q, cs%diag)

2226 if (cs%id_sh_xx_h>0) call post_data(cs%id_sh_xx_h, sh_xx_h, cs%diag)

2227 if (cs%id_sh_xy_q>0) call post_data(cs%id_sh_xy_q, sh_xy_q, cs%diag)

2228 if (cs%id_Ah_q>0) call post_data(cs%id_Ah_q, ah_q, cs%diag)

2229 if (cs%id_Kh_h>0) call post_data(cs%id_Kh_h, kh_h, cs%diag)

2230 if (cs%id_grid_Re_Kh>0) call post_data(cs%id_grid_Re_Kh, grid_re_kh, cs%diag)

2231 if (cs%id_Kh_q>0) call post_data(cs%id_Kh_q, kh_q, cs%diag)

2232 if (cs%use_GME) then ! post barotropic tension and strain

2233 if (cs%id_GME_coeff_h > 0) call post_data(cs%id_GME_coeff_h, gme_coeff_h, cs%diag)

2234 if (cs%id_GME_coeff_q > 0) call post_data(cs%id_GME_coeff_q, gme_coeff_q, cs%diag)

2235 if (cs%id_FrictWork_GME>0) call post_data(cs%id_FrictWork_GME, frictwork_gme, cs%diag)

2236 if (cs%id_dudx_bt > 0) call post_data(cs%id_dudx_bt, dudx_bt, cs%diag)

2237 if (cs%id_dvdy_bt > 0) call post_data(cs%id_dvdy_bt, dvdy_bt, cs%diag)

2238 if (cs%id_dudy_bt > 0) call post_data(cs%id_dudy_bt, dudy_bt, cs%diag)

2239 if (cs%id_dvdx_bt > 0) call post_data(cs%id_dvdx_bt, dvdx_bt, cs%diag)

2240 endif

2241 if (cs%EY24_EBT_BS) then

2242 if (cs%id_visc_limit_h>0) call post_data(cs%id_visc_limit_h, visc_limit_h, cs%diag)

2243 if (cs%id_visc_limit_q>0) call post_data(cs%id_visc_limit_q, visc_limit_q, cs%diag)

2244 if (cs%id_visc_limit_h_frac>0) call post_data(cs%id_visc_limit_h_frac, visc_limit_h_frac, cs%diag)

2245 if (cs%id_visc_limit_q_frac>0) call post_data(cs%id_visc_limit_q_frac, visc_limit_q_frac, cs%diag)

2246 if (cs%id_visc_limit_h_flag>0) call post_data(cs%id_visc_limit_h_flag, visc_limit_h_flag, cs%diag)

2247 if (cs%id_visc_limit_q_flag>0) call post_data(cs%id_visc_limit_q_flag, visc_limit_q_flag, cs%diag)

2248 if (cs%id_BS_coeff_h>0) call post_data(cs%id_BS_coeff_h, bs_coeff_h, cs%diag)

2249 if (cs%id_BS_coeff_q>0) call post_data(cs%id_BS_coeff_q, bs_coeff_q, cs%diag)

2250 endif

2251

2252 if (cs%debug) then

2253 if (cs%Laplacian) then

2254 call hchksum(kh_h, "Kh_h", g%HI, haloshift=1, unscale=us%L_to_m**2*us%s_to_T)

2255 call bchksum(kh_q, "Kh_q", g%HI, haloshift=0, symmetric=.true., unscale=us%L_to_m**2*us%s_to_T)

2256 endif

2257 if (cs%biharmonic) then

2258 call hchksum(ah_h, "Ah_h", g%HI, haloshift=1, unscale=us%L_to_m**4*us%s_to_T)

2259 call bchksum(ah_q, "Ah_q", g%HI, haloshift=0, symmetric=.true., unscale=us%L_to_m**4*us%s_to_T)

2260 endif

2261 endif

2262

2263 if (cs%id_FrictWorkIntz > 0) then

2264 do j=js,je

2265 do i=is,ie ; frictworkintz(i,j) = frictwork(i,j,1) ; enddo

2266 do k=2,nz ; do i=is,ie

2267 frictworkintz(i,j) = frictworkintz(i,j) + frictwork(i,j,k)

2268 enddo ; enddo

2269 enddo

2270 call post_data(cs%id_FrictWorkIntz, frictworkintz, cs%diag)

2271 endif

2272

2273 if (cs%id_FrictWorkIntz_bh > 0) then

2274 do j=js,je

2275 do i=is,ie ; frictworkintz_bh(i,j) = frictwork_bh(i,j,1) ; enddo

2276 do k=2,nz ; do i=is,ie

2277 frictworkintz_bh(i,j) = frictworkintz_bh(i,j) + frictwork_bh(i,j,k)

2278 enddo ; enddo

2279 enddo

2280 call post_data(cs%id_FrictWorkIntz_bh, frictworkintz_bh, cs%diag)

2281 endif

2282

2283 if (present(adp)) then

2284 ! Diagnostics of the fractional thicknesses times momentum budget terms

2285 ! 3D diagnostics of hf_diffu(diffv) are commented because there is no clarity on proper remapping grid option.

2286 ! The code is retained for debugging purposes in the future.

2287 !if (CS%id_hf_diffu > 0) call post_product_u(CS%id_hf_diffu, diffu, ADp%diag_hfrac_u, G, nz, CS%diag)

2288 !if (CS%id_hf_diffv > 0) call post_product_v(CS%id_hf_diffv, diffv, ADp%diag_hfrac_v, G, nz, CS%diag)

2289

2290 ! Diagnostics for thickness-weighted vertically averaged momentum budget terms

2291 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)

2292 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)

2293

2294 ! Diagnostics for the vertical sum of layer thickness x momentum budget terms

2295 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)

2296 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)

2297

2298 ! Diagnostics for thickness x momentum budget terms

2299 if (cs%id_h_diffu > 0) call post_product_u(cs%id_h_diffu, diffu, adp%diag_hu, g, nz, cs%diag)

2300 if (cs%id_h_diffv > 0) call post_product_v(cs%id_h_diffv, diffv, adp%diag_hv, g, nz, cs%diag)

2301

2302 ! Diagnostics for momentum budget terms multiplied by visc_rem_[uv],

2303 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)

2304 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)

2305 endif

2306

2307 if (cs%use_ZB2020) then

2308 call zb2020_lateral_stress(u, v, h, diffu, diffv, g, gv, cs%ZB2020, &

2309 cs%dx2h, cs%dy2h, cs%dx2q, cs%dy2q)

2310 endif

2311

2312end subroutine horizontal_viscosity

2313

2314!> Allocates space for and calculates static variables used by horizontal_viscosity.

2315!! hor_visc_init calculates and stores the values of a number of metric functions that

2316!! are used in horizontal_viscosity.

2317subroutine hor_visc_init(Time, G, GV, US, param_file, diag, CS, ADp)

2318 type(time_type), intent(in) :: time !< Current model time.

2319 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.

2320 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure

2321 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

2322 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time

2323 !! parameters.

2324 type(diag_ctrl), target, intent(inout) :: diag !< Structure to regulate diagnostic output.

2325 type(hor_visc_cs), intent(inout) :: cs !< Horizontal viscosity control structure

2326 type(accel_diag_ptrs), intent(in), optional :: adp !< Acceleration diagnostics

2327

2328 ! u0v is the Laplacian sensitivities to the v velocities at u points, with u0u, v0u, and v0v defined analogously.

2329 real, dimension(SZIB_(G),SZJ_(G)) :: u0u, u0v ! Laplacian sensitivities at u points [L-2 ~> m-2]

2330 real, dimension(SZI_(G),SZJB_(G)) :: v0u, v0v ! Laplacian sensitivities at v points [L-2 ~> m-2]

2331 real :: grid_sp_h2 ! Harmonic mean of the squares of the grid [L2 ~> m2]

2332 real :: grid_sp_h3 ! Harmonic mean of the squares of the grid^(3/2) [L3 ~> m3]

2333 real :: grid_sp_q2 ! spacings at h and q points [L2 ~> m2]

2334 real :: grid_sp_q3 ! spacings at h and q points^(3/2) [L3 ~> m3]

2335 real :: min_grid_sp_h2 ! Minimum value of grid_sp_h2 [L2 ~> m2]

2336 real :: min_grid_sp_h4 ! Minimum value of grid_sp_h2**2 [L4 ~> m4]

2337 real :: kh_limit ! A coefficient [T-1 ~> s-1] used, along with the

2338 ! grid spacing, to limit Laplacian viscosity.

2339 real :: fmax ! maximum absolute value of f at the four

2340 ! vorticity points around a thickness point [T-1 ~> s-1]

2341 real :: boundcorconst ! A constant used when using viscosity to bound the Coriolis accelerations

2342 ! [T2 L-2 ~> s2 m-2]

2343 real :: ah_limit ! coefficient [T-1 ~> s-1] used, along with the

2344 ! grid spacing, to limit biharmonic viscosity

2345 real :: kh ! Lapacian horizontal viscosity [L2 T-1 ~> m2 s-1]

2346 real :: ah ! biharmonic horizontal viscosity [L4 T-1 ~> m4 s-1]

2347 real :: kh_vel_scale ! this speed [L T-1 ~> m s-1] times grid spacing gives Laplacian viscosity

2348 real :: ah_vel_scale ! this speed [L T-1 ~> m s-1] times grid spacing cubed gives biharmonic viscosity

2349 real :: ah_time_scale ! damping time-scale for biharmonic visc [T ~> s]

2350 real :: smag_lap_const ! nondimensional Laplacian Smagorinsky constant [nondim]

2351 real :: smag_bi_const ! nondimensional biharmonic Smagorinsky constant [nondim]

2352 real :: leith_lap_const ! nondimensional Laplacian Leith constant [nondim]

2353 real :: leith_bi_const ! nondimensional biharmonic Leith constant [nondim]

2354 real :: dt ! The dynamics time step [T ~> s]

2355 real :: idt ! The inverse of dt [T-1 ~> s-1]

2356 real :: denom ! work variable; the denominator of a fraction [L-2 ~> m-2] or [L-4 ~> m-4]

2357 real :: maxvel ! largest permitted velocity components [L T-1 ~> m s-1]

2358 real :: bound_cor_vel ! grid-scale velocity variations at which value

2359 ! the quadratically varying biharmonic viscosity

2360 ! balances Coriolis acceleration [L T-1 ~> m s-1]

2361 real :: kh_sin_lat ! Amplitude of latitudinally dependent viscosity [L2 T-1 ~> m2 s-1]

2362 real :: kh_pwr_of_sine ! Power used to raise sin(lat) when using Kh_sin_lat [nondim]

2363 logical :: bound_cor_def ! parameter setting of BOUND_CORIOLIS

2364 logical :: split ! If true, use the split time stepping scheme.

2365 ! If false and USE_GME = True, issue a FATAL error.

2366 logical :: use_meke ! If true, the MEKE parameterization is in use.

2367 logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to

2368 ! recreate the bugs, or if false bugs are only used if actively selected.

2369 real :: backscatter_ro_c ! Coefficient in Rossby number function for backscatter [nondim]

2370 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags

2371 character(len=200) :: inputdir, filename ! Input file names and paths

2372 character(len=80) :: kh_var ! Input variable names

2373 real :: deg2rad ! Converts degrees to radians [radians degree-1]

2374 real :: slat_fn ! sin(lat)**Kh_pwr_of_sine [nondim]

2375 real :: aniso_grid_dir(2) ! Vector (n1,n2) for anisotropic direction [nondim]

2376 integer :: aniso_mode ! Selects the mode for setting the anisotropic direction

2377 integer :: is, ie, js, je, isq, ieq, jsq, jeq, nz

2378 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb

2379 integer :: i, j

2380 ! This include declares and sets the variable "version".

2381# include "version_variable.h"

2382 character(len=40) :: mdl = "MOM_hor_visc" ! module name

2383 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

2384 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

2385 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

2386 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

2387

2388 ! init control structure

2389 call zb2020_init(time, g, gv, us, param_file, diag, cs%ZB2020, cs%use_ZB2020)

2390

2391 cs%initialized = .true.

2392

2393 cs%diag => diag

2394 ! Read parameters and write them to the model log.

2395 call log_version(param_file, mdl, version, "")

2396

2397 call get_param(param_file, mdl, "USE_CIRCULATION_IN_HORVISC", cs%use_circulation, &

2398 "Use circulation theorem to compute vorticity in horvisc module (for ZB20 or Leith)", &

2399 default=.false.)

2400

2401 ! All parameters are read in all cases to enable parameter spelling checks.

2402 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &

2403 "This sets the default value for the various _ANSWER_DATE parameters.", &

2404 default=99991231)

2405

2406 ! Determine whether HOR_VISC_ANSWER_DATE is used, and avoid logging it if it is not used.

2407 call get_param(param_file, mdl, "USE_MEKE", use_meke, &

2408 default=.false., do_not_log=.true.)

2409 backscatter_ro_c = 0.0

2410 if (use_meke) call get_param(param_file, mdl, "MEKE_BACKSCAT_RO_C", backscatter_ro_c, &

2411 "The coefficient in the Rossby number function for scaling the biharmonic "//&

2412 "frictional energy source. Setting to non-zero enables the Rossby number function.", &

2413 units="nondim", default=0.0, do_not_log=.true.)

2414

2415 call get_param(param_file, mdl, "HOR_VISC_ANSWER_DATE", cs%answer_date, &

2416 "The vintage of the order of arithmetic and expressions in the horizontal "//&

2417 "viscosity calculations. Values between 20190102 and 20241201 recover the "//&

2418 "answers from the end of 2018, while higher values use updated and more robust "//&

2419 "forms of the same expressions.", &

2420 default=default_answer_date, do_not_log=(.not.gv%Boussinesq).or.(backscatter_ro_c==0.0))

2421 if (.not.gv%Boussinesq) cs%answer_date = max(cs%answer_date, 20241201)

2422

2423 call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)

2424 call get_param(param_file, mdl, "USE_CONT_THICKNESS", cs%use_cont_thick, &

2425 "If true, use thickness at velocity points from continuity solver. This option "//&

2426 "currently only works with split mode.", default=.false.)

2427 call get_param(param_file, mdl, "USE_CONT_THICKNESS_BUG", cs%use_cont_thick_bug, &

2428 "If true, retain an answer-changing halo update bug when "//&

2429 "USE_CONT_THICKNESS=True. This is not recommended.", &

2430 default=.false., do_not_log=.not.cs%use_cont_thick)

2431

2432 call get_param(param_file, mdl, "LAPLACIAN", cs%Laplacian, &

2433 "If true, use a Laplacian horizontal viscosity.", &

2434 default=.false.)

2435

2436 call get_param(param_file, mdl, "KH", kh, &

2437 "The background Laplacian horizontal viscosity.", &

2438 units="m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s, &

2439 do_not_log=.not.cs%Laplacian)

2440 call get_param(param_file, mdl, "KH_BG_MIN", cs%Kh_bg_min, &

2441 "The minimum value allowed for Laplacian horizontal viscosity, KH.", &

2442 units="m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s, &

2443 do_not_log=.not.cs%Laplacian)

2444 call get_param(param_file, mdl, "KH_VEL_SCALE", kh_vel_scale, &

2445 "The velocity scale which is multiplied by the grid "//&

2446 "spacing to calculate the Laplacian viscosity. "//&

2447 "The final viscosity is the largest of this scaled "//&

2448 "viscosity, the Smagorinsky and Leith viscosities, and KH.", &

2449 units="m s-1", default=0.0, scale=us%m_s_to_L_T, &

2450 do_not_log=.not.cs%Laplacian)

2451 call get_param(param_file, mdl, "KH_SIN_LAT", kh_sin_lat, &

2452 "The amplitude of a latitudinally-dependent background "//&

2453 "viscosity of the form KH_SIN_LAT*(SIN(LAT)**KH_PWR_OF_SINE).", &

2454 units="m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s, &

2455 do_not_log=.not.cs%Laplacian)

2456 call get_param(param_file, mdl, "KH_PWR_OF_SINE", kh_pwr_of_sine, &

2457 "The power used to raise SIN(LAT) when using a latitudinally "//&

2458 "dependent background viscosity.", &

2459 units="nondim", default=4.0, &

2460 do_not_log=.not.(cs%Laplacian .and. (kh_sin_lat>0.)) )

2461 call get_param(param_file, mdl, "SMAGORINSKY_KH", cs%Smagorinsky_Kh, &

2462 "If true, use a Smagorinsky nonlinear eddy viscosity.", &

2463 default=.false., do_not_log=.not.cs%Laplacian)

2464 if (.not.cs%Laplacian) cs%Smagorinsky_Kh = .false.

2465 call get_param(param_file, mdl, "SMAG_LAP_CONST", smag_lap_const, &

2466 "The nondimensional Laplacian Smagorinsky constant, "//&

2467 "often 0.15.", units="nondim", default=0.0, &

2468 fail_if_missing=cs%Smagorinsky_Kh, do_not_log=.not.cs%Smagorinsky_Kh)

2469 call get_param(param_file, mdl, "LEITH_KH", cs%Leith_Kh, &

2470 "If true, use a Leith nonlinear eddy viscosity.", &

2471 default=.false., do_not_log=.not.cs%Laplacian)

2472 if (.not.cs%Laplacian) cs%Leith_Kh = .false.

2473 call get_param(param_file, mdl, "LEITH_LAP_CONST", leith_lap_const, &

2474 "The nondimensional Laplacian Leith constant, "//&

2475 "often set to 1.0", units="nondim", default=0.0, &

2476 fail_if_missing=cs%Leith_Kh, do_not_log=.not.cs%Leith_Kh)

2477 call get_param(param_file, mdl, "RES_SCALE_MEKE_VISC", cs%res_scale_MEKE, &

2478 "If true, the viscosity contribution from MEKE is scaled by "//&

2479 "the resolution function.", default=.false., &

2480 do_not_log=.not.(cs%Laplacian.and.use_meke))

2481 if (.not.(cs%Laplacian.and.use_meke)) cs%res_scale_MEKE = .false.

2482

2483 call get_param(param_file, mdl, "BOUND_KH", cs%bound_Kh, &

2484 "If true, the Laplacian coefficient is locally limited "//&

2485 "to be stable.", default=.true., do_not_log=.not.cs%Laplacian)

2486 call get_param(param_file, mdl, "EY24_EBT_BS", cs%EY24_EBT_BS, &

2487 "If true, use the backscatter scheme (EBT mode with kill switch) "//&

2488 "developed by Yankovsky et al. (2024). ", &

2489 default=.false., do_not_log=.not.cs%Laplacian)

2490 if (.not.cs%Laplacian) cs%bound_Kh = .false.

2491 if (.not.(cs%Laplacian.and.use_meke)) cs%EY24_EBT_BS = .false.

2492 call get_param(param_file, mdl, "ANISOTROPIC_VISCOSITY", cs%anisotropic, &

2493 "If true, allow anistropic viscosity in the Laplacian "//&

2494 "horizontal viscosity.", default=.false., &

2495 do_not_log=.not.cs%Laplacian)

2496 if (.not.cs%Laplacian) cs%anisotropic = .false. ! This replicates the prior code, but is it intended?

2497 call get_param(param_file, mdl, "ADD_LES_VISCOSITY", cs%add_LES_viscosity, &

2498 "If true, adds the viscosity from Smagorinsky and Leith to the "//&

2499 "background viscosity instead of taking the maximum.", default=.false., &

2500 do_not_log=.not.cs%Laplacian)

2501

2502 call get_param(param_file, mdl, "KH_ANISO", cs%Kh_aniso, &

2503 "The background Laplacian anisotropic horizontal viscosity.", &

2504 units="m2 s-1", default=0.0, scale=us%m_to_L**2*us%T_to_s, &

2505 do_not_log=.not.cs%anisotropic)

2506 call get_param(param_file, mdl, "ANISOTROPIC_MODE", aniso_mode, &

2507 "Selects the mode for setting the direction of anisotropy.\n"//&

2508 "\t 0 - Points along the grid i-direction.\n"//&

2509 "\t 1 - Points towards East.\n"//&

2510 "\t 2 - Points along the flow direction, U/|U|.", &

2511 default=0, do_not_log=.not.cs%anisotropic)

2512 if (aniso_mode == 0) then

2513 call get_param(param_file, mdl, "ANISO_GRID_DIR", aniso_grid_dir, &

2514 "The vector pointing in the direction of anisotropy for horizontal viscosity. "//&

2515 "n1,n2 are the i,j components relative to the grid.", &

2516 units="nondim", fail_if_missing=cs%anisotropic, do_not_log=.not.cs%anisotropic)

2517 elseif (aniso_mode == 1) then

2518 call get_param(param_file, mdl, "ANISO_GRID_DIR", aniso_grid_dir, &

2519 "The vector pointing in the direction of anisotropy for horizontal viscosity. "//&

2520 "n1,n2 are the i,j components relative to the spherical coordinates.", &

2521 units="nondim", fail_if_missing=cs%anisotropic, do_not_log=.not.cs%anisotropic)

2522 else

2523 call get_param(param_file, mdl, "ANISO_GRID_DIR", aniso_grid_dir, &

2524 "The vector pointing in the direction of anisotropy for horizontal viscosity.", &

2525 units="nondim", fail_if_missing=.false., do_not_log=.true.)

2526 endif

2527

2528 call get_param(param_file, mdl, "BIHARMONIC", cs%biharmonic, &

2529 "If true, use a biharmonic horizontal viscosity. "//&

2530 "BIHARMONIC may be used with LAPLACIAN.", &

2531 default=.true.)

2532 call get_param(param_file, mdl, "AH", ah, &

2533 "The background biharmonic horizontal viscosity.", &

2534 units="m4 s-1", default=0.0, scale=us%m_to_L**4*us%T_to_s, &

2535 do_not_log=.not.cs%biharmonic)

2536 call get_param(param_file, mdl, "AH_VEL_SCALE", ah_vel_scale, &

2537 "The velocity scale which is multiplied by the cube of "//&

2538 "the grid spacing to calculate the biharmonic viscosity. "//&

2539 "The final viscosity is the largest of this scaled "//&

2540 "viscosity, the Smagorinsky and Leith viscosities, and AH.", &

2541 units="m s-1", default=0.0, scale=us%m_s_to_L_T, do_not_log=.not.cs%biharmonic)

2542 call get_param(param_file, mdl, "AH_TIME_SCALE", ah_time_scale, &

2543 "A time scale whose inverse is multiplied by the fourth "//&

2544 "power of the grid spacing to calculate biharmonic viscosity. "//&

2545 "The final viscosity is the largest of all viscosity "//&

2546 "formulations in use. 0.0 means that it's not used.", &

2547 units="s", default=0.0, scale=us%s_to_T, do_not_log=.not.cs%biharmonic)

2548 call get_param(param_file, mdl, "SMAGORINSKY_AH", cs%Smagorinsky_Ah, &

2549 "If true, use a biharmonic Smagorinsky nonlinear eddy "//&

2550 "viscosity.", default=.false., do_not_log=.not.cs%biharmonic)

2551 if (.not.cs%biharmonic) cs%Smagorinsky_Ah = .false.

2552 call get_param(param_file, mdl, "LEITH_AH", cs%Leith_Ah, &

2553 "If true, use a biharmonic Leith nonlinear eddy "//&

2554 "viscosity.", default=.false., do_not_log=.not.cs%biharmonic)

2555 if (.not.cs%biharmonic) cs%Leith_Ah = .false.

2556 call get_param(param_file, mdl, "USE_LEITHY", cs%use_Leithy, &

2557 "If true, use a biharmonic Leith nonlinear eddy "//&

2558 "viscosity together with a harmonic backscatter.", &

2559 default=.false.)

2560 call get_param(param_file, mdl, "BOUND_AH", cs%bound_Ah, &

2561 "If true, the biharmonic coefficient is locally limited "//&

2562 "to be stable.", default=.true., do_not_log=.not.cs%biharmonic)

2563 if (.not.cs%biharmonic) cs%bound_Ah = .false.

2564 call get_param(param_file, mdl, "RE_AH", cs%Re_Ah, &

2565 "If nonzero, the biharmonic coefficient is scaled "//&

2566 "so that the biharmonic Reynolds number is equal to this.", &

2567 units="nondim", default=0.0, do_not_log=.not.cs%biharmonic)

2568

2569 call get_param(param_file, mdl, "BACKSCATTER_UNDERBOUND", cs%backscatter_underbound, &

2570 "If true, the bounds on the biharmonic viscosity are allowed to "//&

2571 "increase where the Laplacian viscosity is negative (due to backscatter "//&

2572 "parameterizations) beyond the largest timestep-dependent stable values of "//&

2573 "biharmonic viscosity when no Laplacian viscosity is applied. The default "//&

2574 "is true for historical reasons, but this option probably should not be used "//&

2575 "because it can contribute to numerical instabilities.", &

2576 default=.false., do_not_log=.not.((cs%bound_Kh).and.(cs%bound_Ah)))

2577

2578 call get_param(param_file, mdl, "SMAG_BI_CONST",smag_bi_const, &

2579 "The nondimensional biharmonic Smagorinsky constant, "//&

2580 "typically 0.015 - 0.06.", units="nondim", default=0.0, &

2581 fail_if_missing=cs%Smagorinsky_Ah, do_not_log=.not.cs%Smagorinsky_Ah)

2582

2583 call get_param(param_file, mdl, "USE_BETA_IN_LEITH", cs%use_beta_in_Leith, &

2584 "If true, include the beta term in the Leith nonlinear eddy viscosity.", &

2585 default=cs%Leith_Kh, do_not_log=.not.(cs%Leith_Kh .or. cs%Leith_Ah) )

2586 call get_param(param_file, mdl, "MODIFIED_LEITH", cs%modified_Leith, &

2587 "If true, add a term to Leith viscosity which is "//&

2588 "proportional to the gradient of divergence.", &

2589 default=.false., do_not_log=.not.(cs%Leith_Kh .or. cs%Leith_Ah) )

2590 call get_param(param_file, mdl, "USE_QG_LEITH_VISC", cs%use_QG_Leith_visc, &

2591 "If true, use QG Leith nonlinear eddy viscosity.", &

2592 default=.false., do_not_log=.not.(cs%Leith_Kh .or. cs%Leith_Ah) )

2593! if (CS%use_QG_Leith_visc) then

2594! call MOM_error(FATAL, "USE_QG_LEITH_VISC=True activates code that is a work-in-progress and "//&

2595! "should not be used until a number of bugs are fixed. Specifically it does not "//&

2596! "reproduce across PE count or layout, and may use arrays that have not been properly "//&

2597! "set or allocated. See github.com/mom-ocean/MOM6/issues/1590 for a discussion.")

2598! endif

2599 if (cs%use_QG_Leith_visc .and. .not. (cs%Leith_Kh .or. cs%Leith_Ah) ) then

2600 call mom_error(fatal, "MOM_hor_visc.F90, hor_visc_init: "//&

2601 "LEITH_KH or LEITH_AH must be True when USE_QG_LEITH_VISC=True.")

2602 endif

2603

2604 call get_param(param_file, mdl, "BOUND_CORIOLIS", bound_cor_def, default=.false.)

2605 call get_param(param_file, mdl, "BOUND_CORIOLIS_BIHARM", cs%bound_Coriolis, &

2606 "If true use a viscosity that increases with the square "//&

2607 "of the velocity shears, so that the resulting viscous "//&

2608 "drag is of comparable magnitude to the Coriolis terms "//&

2609 "when the velocity differences between adjacent grid "//&

2610 "points is 0.5*BOUND_CORIOLIS_VEL. The default is the "//&

2611 "value of BOUND_CORIOLIS (or false).", default=bound_cor_def, &

2612 do_not_log=.not.cs%Smagorinsky_Ah)

2613 if (.not.cs%Smagorinsky_Ah) cs%bound_Coriolis = .false.

2614 call get_param(param_file, mdl, "MAXVEL", maxvel, &

2615 units="m s-1", default=3.0e8, scale=us%m_s_to_L_T)

2616 call get_param(param_file, mdl, "BOUND_CORIOLIS_VEL", bound_cor_vel, &

2617 "The velocity scale at which BOUND_CORIOLIS_BIHARM causes "//&

2618 "the biharmonic drag to have comparable magnitude to the "//&

2619 "Coriolis acceleration. The default is set by MAXVEL.", &

2620 units="m s-1", default=maxvel*us%L_T_to_m_s, scale=us%m_s_to_L_T, &

2621 do_not_log=.not.(cs%Smagorinsky_Ah .and. cs%bound_Coriolis))

2622 call get_param(param_file, mdl, "LEITH_BI_CONST", leith_bi_const, &

2623 "The nondimensional biharmonic Leith constant, "//&

2624 "typical values are thus far undetermined.", units="nondim", default=0.0, &

2625 fail_if_missing=(cs%Leith_Ah .or. cs%use_Leithy), &

2626 do_not_log=.not.(cs%Leith_Ah .or. cs%use_Leithy))

2627 call get_param(param_file, mdl, "USE_LAND_MASK_FOR_HVISC", cs%use_land_mask, &

2628 "If true, use the land mask for the computation of thicknesses "//&

2629 "at velocity locations. This eliminates the dependence on arbitrary "//&

2630 "values over land or outside of the domain.", default=.true.)

2631 call get_param(param_file, mdl, "HORVISC_BOUND_COEF", cs%bound_coef, &

2632 "The nondimensional coefficient of the ratio of the "//&

2633 "viscosity bounds to the theoretical maximum for "//&

2634 "stability without considering other terms.", units="nondim", &

2635 default=0.8, do_not_log=.not.(cs%bound_Ah .or. cs%bound_Kh))

2636 call get_param(param_file, mdl, "KILL_SWITCH_COEF", cs%KS_coef, &

2637 "A nondimensional coefficient on the biharmonic viscosity that "// &

2638 "sets the kill switch for backscatter. Default is 1.0.", units="nondim", &

2639 default=1.0, do_not_log=.not.(cs%EY24_EBT_BS))

2640 call get_param(param_file, mdl, "NOSLIP", cs%no_slip, &

2641 "If true, no slip boundary conditions are used; otherwise "//&

2642 "free slip boundary conditions are assumed. The "//&

2643 "implementation of the free slip BCs on a C-grid is much "//&

2644 "cleaner than the no slip BCs. The use of free slip BCs "//&

2645 "is strongly encouraged, and no slip BCs are not used with "//&

2646 "the biharmonic viscosity.", default=.false.)

2647 call get_param(param_file, mdl, "USE_KH_BG_2D", cs%use_Kh_bg_2d, &

2648 "If true, read a file containing 2-d background harmonic "//&

2649 "viscosities. The final viscosity is the maximum of the other "//&

2650 "terms and this background value.", default=.false., do_not_log=.not.cs%Laplacian)

2651 if (.not.cs%Laplacian) cs%use_Kh_bg_2d = .false.

2652 call get_param(param_file, mdl, "KH_BG_2D_BUG", cs%Kh_bg_2d_bug, &

2653 "If true, retain an answer-changing horizontal indexing bug in setting "//&

2654 "the corner-point viscosities when USE_KH_BG_2D=True. This is "//&

2655 "not recommended.", default=.false., do_not_log=.not.cs%use_Kh_bg_2d)

2656 call get_param(param_file, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &

2657 default=.true., do_not_log=.true.) ! This is logged from MOM.F90.

2658 call get_param(param_file, mdl, "FRICTWORK_BUG", cs%FrictWork_bug, &

2659 "If true, retain an answer-changing bug in calculating the FrictWork, "//&

2660 "which cancels the h in thickness flux and the h at velocity point. This is "//&

2661 "not recommended.", default=.false.)

2662 call get_param(param_file, mdl, "OBC_SPECIFIED_STRAIN_BUG", cs%OBC_strain_bug, &

2663 "If true, recover a bug that specified shear strain option at open boundaries "//&

2664 "cannot be applied.", default=.true.)

2665 call get_param(param_file, mdl, "USE_GME", cs%use_GME, &

2666 "If true, use the GM+E backscatter scheme in association \n"//&

2667 "with the Gent and McWilliams parameterization.", default=.false.)

2668 call get_param(param_file, mdl, "SPLIT", split, &

2669 "Use the split time stepping if true.", default=.true., do_not_log=.true.)

2670 if (cs%use_Leithy) then

2671 if (.not.(cs%biharmonic .and. cs%Laplacian)) then

2672 call mom_error(fatal, "MOM_hor_visc.F90, hor_visc_init: "//&

2673 "LAPLACIAN and BIHARMONIC must both be True when USE_LEITHY=True.")

2674 endif

2675 endif

2676 call get_param(param_file, mdl, "LEITHY_CK", cs%c_K, &

2677 "Fraction of biharmonic dissipation that gets backscattered, "//&

2678 "in Leith+E.", units="nondim", default=1.0, do_not_log=.not.cs%use_Leithy)

2679 call get_param(param_file, mdl, "SMOOTH_AH", cs%smooth_Ah, &

2680 "If true, Ah and m_leithy are smoothed within Leith+E. This requires "//&

2681 "lots of blocking communications, which can be expensive", &

2682 default=.true., do_not_log=.not.cs%use_Leithy)

2683

2684 if (cs%use_GME .and. .not.split) call mom_error(fatal,"ERROR: Currently, USE_GME = True "// &

2685 "cannot be used with SPLIT=False.")

2686

2687 if (cs%use_GME) then

2688 call get_param(param_file, mdl, "GME_NUM_SMOOTHINGS", cs%num_smooth_gme, &

2689 "Number of smoothing passes for the GME fluxes.", &

2690 default=1)

2691 call get_param(param_file, mdl, "GME_H0", cs%GME_h0, &

2692 "The strength of GME tapers quadratically to zero when the bathymetric "//&

2693 "depth is shallower than GME_H0.", &

2694 units="m", scale=gv%m_to_H, default=1000.0)

2695 call get_param(param_file, mdl, "GME_EFFICIENCY", cs%GME_efficiency, &

2696 "The nondimensional prefactor multiplying the GME coefficient.", &

2697 units="nondim", default=1.0)

2698 call get_param(param_file, mdl, "GME_LIMITER", cs%GME_limiter, &

2699 "The absolute maximum value the GME coefficient is allowed to take.", &

2700 units="m2 s-1", scale=us%m_to_L**2*us%T_to_s, default=1.0e7)

2701 endif

2702

2703 if (cs%Laplacian .or. cs%biharmonic) then

2704 call get_param(param_file, mdl, "DT", dt, &

2705 "The (baroclinic) dynamics time step.", units="s", scale=us%s_to_T, &

2706 fail_if_missing=.true.)

2707 idt = 1.0 / dt

2708 endif

2709 call get_param(param_file, mdl, "KILL_SWITCH_TIMESCALE", cs%KS_timescale, &

2710 "A timescale for computing the CFL limit for viscosity "// &

2711 "that determines when backscatter is shut off. Default is DT.", &

2712 default= dt , units="s", scale=us%s_to_T, do_not_log=.not.(cs%EY24_EBT_BS))

2713

2714 if (cs%no_slip .and. cs%biharmonic) &

2715 call mom_error(fatal,"ERROR: NOSLIP and BIHARMONIC cannot be defined "// &

2716 "at the same time in MOM.")

2717 if (.not.(cs%Laplacian .or. cs%biharmonic)) then

2718 ! Only issue inviscid warning if not in single column mode (usually 2x2 domain)

2719 if ( max(g%domain%niglobal, g%domain%njglobal)>2 ) call mom_error(warning, &

2720 "hor_visc_init: It is usually a very bad idea not to use either "//&

2721 "LAPLACIAN or BIHARMONIC viscosity.")

2722 return ! We are not using either Laplacian or Bi-harmonic lateral viscosity

2723 endif

2724 deg2rad = atan(1.0) / 45.

2725 alloc_(cs%dx2h(isd:ied,jsd:jed)) ; cs%dx2h(:,:) = 0.0

2726 alloc_(cs%dy2h(isd:ied,jsd:jed)) ; cs%dy2h(:,:) = 0.0

2727 alloc_(cs%dx2q(isdb:iedb,jsdb:jedb)) ; cs%dx2q(:,:) = 0.0

2728 alloc_(cs%dy2q(isdb:iedb,jsdb:jedb)) ; cs%dy2q(:,:) = 0.0

2729 alloc_(cs%dx_dyT(isd:ied,jsd:jed)) ; cs%dx_dyT(:,:) = 0.0

2730 alloc_(cs%dy_dxT(isd:ied,jsd:jed)) ; cs%dy_dxT(:,:) = 0.0

2731 alloc_(cs%dx_dyBu(isdb:iedb,jsdb:jedb)) ; cs%dx_dyBu(:,:) = 0.0

2732 alloc_(cs%dy_dxBu(isdb:iedb,jsdb:jedb)) ; cs%dy_dxBu(:,:) = 0.0

2733 if (cs%Laplacian) then

2734 alloc_(cs%grid_sp_h2(isd:ied,jsd:jed)) ; cs%grid_sp_h2(:,:) = 0.0

2735 alloc_(cs%Kh_bg_xx(isd:ied,jsd:jed)) ; cs%Kh_bg_xx(:,:) = 0.0

2736 alloc_(cs%Kh_bg_xy(isdb:iedb,jsdb:jedb)) ; cs%Kh_bg_xy(:,:) = 0.0

2737 if (cs%bound_Kh .or. cs%EY24_EBT_BS) then

2738 allocate(cs%Kh_Max_xx(isd:ied,jsd:jed), source=0.0)

2739 allocate(cs%Kh_Max_xy(isdb:iedb,jsdb:jedb), source=0.0)

2740 endif

2741 if (cs%Smagorinsky_Kh .or. cs%EY24_EBT_BS) then

2742 allocate(cs%Laplac2_const_xx(isd:ied,jsd:jed), source=0.0)

2743 allocate(cs%Laplac2_const_xy(isdb:iedb,jsdb:jedb), source=0.0)

2744 endif

2745 if (cs%Leith_Kh) then

2746 allocate(cs%Laplac3_const_xx(isd:ied,jsd:jed), source=0.0)

2747 allocate(cs%Laplac3_const_xy(isdb:iedb,jsdb:jedb), source=0.0)

2748 endif

2749 endif

2750 alloc_(cs%reduction_xx(isd:ied,jsd:jed)) ; cs%reduction_xx(:,:) = 0.0

2751 alloc_(cs%reduction_xy(isdb:iedb,jsdb:jedb)) ; cs%reduction_xy(:,:) = 0.0

2752

2753 cs%dynamic_aniso = .false.

2754 if (cs%anisotropic) then

2755 allocate(cs%n1n2_h(isd:ied,jsd:jed), source=0.0)

2756 allocate(cs%n1n1_m_n2n2_h(isd:ied,jsd:jed), source=0.0)

2757 allocate(cs%n1n2_q(isdb:iedb,jsdb:jedb), source=0.0)

2758 allocate(cs%n1n1_m_n2n2_q(isdb:iedb,jsdb:jedb), source=0.0)

2759 select case (aniso_mode)

2760 case (0)

2761 call align_aniso_tensor_to_grid(cs, aniso_grid_dir(1), aniso_grid_dir(2))

2762 case (1)

2763 ! call align_aniso_tensor_to_grid(CS, aniso_grid_dir(1), aniso_grid_dir(2))

2764 case (2)

2765 cs%dynamic_aniso = .true.

2766 case default

2767 call mom_error(fatal, "MOM_hor_visc: "//&

2768 "Runtime parameter ANISOTROPIC_MODE is out of range.")

2769 end select

2770 endif

2771

2772 call get_param(param_file, mdl, "KH_BG_2D_FILENAME", filename, &

2773 'The filename containing a 2d map of "Kh".', &

2774 default='KH_background_2d.nc', do_not_log=.not.cs%use_Kh_bg_2d)

2775 call get_param(param_file, mdl, "KH_BG_2D_VARNAME", kh_var, &

2776 'The name in the input file of the horizontal viscosity variable.', &

2777 default='Kh', do_not_log=.not.cs%use_Kh_bg_2d)

2778

2779 if (cs%use_Kh_bg_2d) then

2780 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")

2781 inputdir = slasher(inputdir)

2782 allocate(cs%Kh_bg_2d(isd:ied,jsd:jed), source=0.0)

2783 call mom_read_data(trim(inputdir)//trim(filename), kh_var, cs%Kh_bg_2d, &

2784 g%domain, timelevel=1, scale=us%m_to_L**2*us%T_to_s)

2785 call pass_var(cs%Kh_bg_2d, g%domain)

2786 endif

2787 if (cs%biharmonic) then

2788 alloc_(cs%Idx2dyCu(isdb:iedb,jsd:jed)) ; cs%Idx2dyCu(:,:) = 0.0

2789 alloc_(cs%Idx2dyCv(isd:ied,jsdb:jedb)) ; cs%Idx2dyCv(:,:) = 0.0

2790 alloc_(cs%Idxdy2u(isdb:iedb,jsd:jed)) ; cs%Idxdy2u(:,:) = 0.0

2791 alloc_(cs%Idxdy2v(isd:ied,jsdb:jedb)) ; cs%Idxdy2v(:,:) = 0.0

2792 alloc_(cs%Ah_bg_xx(isd:ied,jsd:jed)) ; cs%Ah_bg_xx(:,:) = 0.0

2793 alloc_(cs%Ah_bg_xy(isdb:iedb,jsdb:jedb)) ; cs%Ah_bg_xy(:,:) = 0.0

2794 alloc_(cs%grid_sp_h3(isd:ied,jsd:jed)) ; cs%grid_sp_h3(:,:) = 0.0

2795 if (cs%bound_Ah) then

2796 allocate(cs%Ah_Max_xx(isd:ied,jsd:jed), source=0.0)

2797 allocate(cs%Ah_Max_xy(isdb:iedb,jsdb:jedb), source=0.0)

2798 endif

2799 if (cs%EY24_EBT_BS) then

2800 allocate(cs%Ah_Max_xx_KS(isd:ied,jsd:jed), source=0.0)

2801 allocate(cs%Ah_Max_xy_KS(isdb:iedb,jsdb:jedb), source=0.0)

2802 endif

2803 if (cs%Smagorinsky_Ah) then

2804 allocate(cs%Biharm_const_xx(isd:ied,jsd:jed), source=0.0)

2805 allocate(cs%Biharm_const_xy(isdb:iedb,jsdb:jedb), source=0.0)

2806 if (cs%bound_Coriolis) then

2807 allocate(cs%Biharm_const2_xx(isd:ied,jsd:jed), source=0.0)

2808 allocate(cs%Biharm_const2_xy(isdb:iedb,jsdb:jedb), source=0.0)

2809 endif

2810 endif

2811 if ((cs%Leith_Ah) .or. (cs%use_Leithy)) then

2812 allocate(cs%biharm6_const_xx(isd:ied,jsd:jed), source=0.0)

2813 allocate(cs%biharm6_const_xy(isdb:iedb,jsdb:jedb), source=0.0)

2814 endif

2815 if (cs%use_Leithy) then

2816 allocate(cs%m_const_leithy(isd:ied,jsd:jed), source=0.0)

2817 allocate(cs%m_leithy_max(isd:ied,jsd:jed), source=0.0)

2818 endif

2819 if (cs%Re_Ah > 0.0) then

2820 allocate(cs%Re_Ah_const_xx(isd:ied,jsd:jed), source=0.0)

2821 allocate(cs%Re_Ah_const_xy(isdb:iedb,jsdb:jedb), source=0.0)

2822 endif

2823 endif

2824 do j=js-2,jeq+1 ; do i=is-2,ieq+1

2825 cs%dx2q(i,j) = g%dxBu(i,j)*g%dxBu(i,j) ; cs%dy2q(i,j) = g%dyBu(i,j)*g%dyBu(i,j)

2826 enddo ; enddo

2827

2828 if (((cs%Leith_Kh) .or. (cs%Leith_Ah) .or. (cs%use_Leithy)) .and. &

2829 ((g%isc-g%isd < 3) .or. (g%isc-g%isd < 3))) call mom_error(fatal, &

2830 "The minimum halo size is 3 when a Leith viscosity is being used.")

2831 if (cs%use_Leithy) then

2832 do j=js-3,jeq+2 ; do i=is-3,ieq+2

2833 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)

2834 enddo ; enddo

2835 elseif ((cs%Leith_Kh) .or. (cs%Leith_Ah)) then

2836 do j=jsq-2,jeq+2 ; do i=isq-2,ieq+2

2837 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)

2838 enddo ; enddo

2839 else

2840 do j=js-2,jeq+1 ; do i=is-2,ieq+1

2841 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)

2842 enddo ; enddo

2843 endif

2844

2845 do j=js-2,jeq+2 ; do i=is-2,ieq+2

2846 cs%dx2h(i,j) = g%dxT(i,j)*g%dxT(i,j) ; cs%dy2h(i,j) = g%dyT(i,j)*g%dyT(i,j)

2847 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)

2848 enddo ; enddo

2849 do j=jsq,jeq+1 ; do i=isq,ieq+1

2850 cs%reduction_xx(i,j) = 1.0

2851 if ((g%dy_Cu(i,j) > 0.0) .and. (g%dy_Cu(i,j) < g%dyCu(i,j)) .and. &

2852 (g%dy_Cu(i,j) < g%dyCu(i,j) * cs%reduction_xx(i,j))) &

2853 cs%reduction_xx(i,j) = g%dy_Cu(i,j) / (g%dyCu(i,j))

2854 if ((g%dy_Cu(i-1,j) > 0.0) .and. (g%dy_Cu(i-1,j) < g%dyCu(i-1,j)) .and. &

2855 (g%dy_Cu(i-1,j) < g%dyCu(i-1,j) * cs%reduction_xx(i,j))) &

2856 cs%reduction_xx(i,j) = g%dy_Cu(i-1,j) / (g%dyCu(i-1,j))

2857 if ((g%dx_Cv(i,j) > 0.0) .and. (g%dx_Cv(i,j) < g%dxCv(i,j)) .and. &

2858 (g%dx_Cv(i,j) < g%dxCv(i,j) * cs%reduction_xx(i,j))) &

2859 cs%reduction_xx(i,j) = g%dx_Cv(i,j) / (g%dxCv(i,j))

2860 if ((g%dx_Cv(i,j-1) > 0.0) .and. (g%dx_Cv(i,j-1) < g%dxCv(i,j-1)) .and. &

2861 (g%dx_Cv(i,j-1) < g%dxCv(i,j-1) * cs%reduction_xx(i,j))) &

2862 cs%reduction_xx(i,j) = g%dx_Cv(i,j-1) / (g%dxCv(i,j-1))

2863 enddo ; enddo

2864 do j=js-1,jeq ; do i=is-1,ieq

2865 cs%reduction_xy(i,j) = 1.0

2866 if ((g%dy_Cu(i,j) > 0.0) .and. (g%dy_Cu(i,j) < g%dyCu(i,j)) .and. &

2867 (g%dy_Cu(i,j) < g%dyCu(i,j) * cs%reduction_xy(i,j))) &

2868 cs%reduction_xy(i,j) = g%dy_Cu(i,j) / (g%dyCu(i,j))

2869 if ((g%dy_Cu(i,j+1) > 0.0) .and. (g%dy_Cu(i,j+1) < g%dyCu(i,j+1)) .and. &

2870 (g%dy_Cu(i,j+1) < g%dyCu(i,j+1) * cs%reduction_xy(i,j))) &

2871 cs%reduction_xy(i,j) = g%dy_Cu(i,j+1) / (g%dyCu(i,j+1))

2872 if ((g%dx_Cv(i,j) > 0.0) .and. (g%dx_Cv(i,j) < g%dxCv(i,j)) .and. &

2873 (g%dx_Cv(i,j) < g%dxCv(i,j) * cs%reduction_xy(i,j))) &

2874 cs%reduction_xy(i,j) = g%dx_Cv(i,j) / (g%dxCv(i,j))

2875 if ((g%dx_Cv(i+1,j) > 0.0) .and. (g%dx_Cv(i+1,j) < g%dxCv(i+1,j)) .and. &

2876 (g%dx_Cv(i+1,j) < g%dxCv(i+1,j) * cs%reduction_xy(i,j))) &

2877 cs%reduction_xy(i,j) = g%dx_Cv(i+1,j) / (g%dxCv(i+1,j))

2878 enddo ; enddo

2879 if (cs%Laplacian) then

2880 ! The 0.3 below was 0.4 in MOM1.10. The change in hq requires

2881 ! this to be less than 1/3, rather than 1/2 as before.

2882 if (cs%bound_Kh .or. cs%bound_Ah) kh_limit = 0.3 / (dt*4.0)

2883 ! Calculate and store the background viscosity at h-points

2884

2885 min_grid_sp_h2 = huge(1.)

2886 do j=js-1,jeq+1 ; do i=is-1,ieq+1

2887 ! Static factors in the Smagorinsky and Leith schemes

2888 grid_sp_h2 = (2.0*cs%dx2h(i,j)*cs%dy2h(i,j)) / (cs%dx2h(i,j) + cs%dy2h(i,j))

2889 cs%grid_sp_h2(i,j) = grid_sp_h2

2890 grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)

2891 if (cs%Smagorinsky_Kh) cs%Laplac2_const_xx(i,j) = smag_lap_const * grid_sp_h2

2892 if (cs%Leith_Kh) cs%Laplac3_const_xx(i,j) = leith_lap_const * grid_sp_h3

2893 ! Maximum of constant background and MICOM viscosity

2894 cs%Kh_bg_xx(i,j) = max(kh, kh_vel_scale * sqrt(grid_sp_h2))

2895 ! Use the larger of the above and values read from a file

2896 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))

2897 ! Use the larger of the above and a function of sin(latitude)

2898 if (kh_sin_lat>0.) then

2899 slat_fn = abs( sin( deg2rad * g%geoLatT(i,j) ) ) ** kh_pwr_of_sine

2900 cs%Kh_bg_xx(i,j) = max(kh_sin_lat * slat_fn, cs%Kh_bg_xx(i,j))

2901 endif

2902 min_grid_sp_h2 = min(grid_sp_h2, min_grid_sp_h2)

2903 enddo ; enddo

2904 call min_across_pes(min_grid_sp_h2)

2905

2906 ! Calculate and store the background viscosity at q-points

2907 do j=js-1,jeq ; do i=is-1,ieq

2908 ! Static factors in the Smagorinsky and Leith schemes

2909 grid_sp_q2 = (2.0*cs%dx2q(i,j)*cs%dy2q(i,j)) / (cs%dx2q(i,j) + cs%dy2q(i,j))

2910 grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)

2911 if (cs%Smagorinsky_Kh) cs%Laplac2_const_xy(i,j) = smag_lap_const * grid_sp_q2

2912 if (cs%Leith_Kh) cs%Laplac3_const_xy(i,j) = leith_lap_const * grid_sp_q3

2913 ! Maximum of constant background and MICOM viscosity

2914 cs%Kh_bg_xy(i,j) = max(kh, kh_vel_scale * sqrt(grid_sp_q2))

2915 ! Use the larger of the above and values read from a file

2916 if (cs%use_Kh_bg_2d) then

2917 if (cs%Kh_bg_2d_bug) then

2918 ! This option is unambiguously wrong but is needed to recover old answers

2919 cs%Kh_bg_xy(i,j) = max(cs%Kh_bg_2d(i,j), cs%Kh_bg_xy(i,j))

2920 else

2921 cs%Kh_bg_xy(i,j) = max(cs%Kh_bg_xy(i,j), &

2922 0.25*((cs%Kh_bg_2d(i,j) + cs%Kh_bg_2d(i+1,j+1)) + &

2923 (cs%Kh_bg_2d(i+1,j) + cs%Kh_bg_2d(i,j+1))) )

2924 endif

2925 endif

2926

2927 ! Use the larger of the above and a function of sin(latitude)

2928 if (kh_sin_lat>0.) then

2929 slat_fn = abs( sin( deg2rad * g%geoLatBu(i,j) ) ) ** kh_pwr_of_sine

2930 cs%Kh_bg_xy(i,j) = max(kh_sin_lat * slat_fn, cs%Kh_bg_xy(i,j))

2931 endif

2932 enddo ; enddo

2933 endif

2934 if (cs%biharmonic) then

2935 do j=js-1,jeq+1 ; do i=is-2,ieq+1

2936 cs%Idx2dyCu(i,j) = (g%IdxCu(i,j)*g%IdxCu(i,j)) * g%IdyCu(i,j)

2937 cs%Idxdy2u(i,j) = g%IdxCu(i,j) * (g%IdyCu(i,j)*g%IdyCu(i,j))

2938 enddo ; enddo

2939 do j=js-2,jeq+1 ; do i=is-1,ieq+1

2940 cs%Idx2dyCv(i,j) = (g%IdxCv(i,j)*g%IdxCv(i,j)) * g%IdyCv(i,j)

2941 cs%Idxdy2v(i,j) = g%IdxCv(i,j) * (g%IdyCv(i,j)*g%IdyCv(i,j))

2942 enddo ; enddo

2943 cs%Ah_bg_xy(:,:) = 0.0

2944 ! The 0.3 below was 0.4 in HIM 1.10. The change in hq requires

2945 ! this to be less than 1/3, rather than 1/2 as before.

2946 if (cs%bound_Ah) ah_limit = 0.3 / (dt*64.0)

2947 if (cs%Smagorinsky_Ah .and. cs%bound_Coriolis) &

2948 boundcorconst = 1.0 / (5.0*(bound_cor_vel*bound_cor_vel))

2949

2950 min_grid_sp_h4 = huge(1.)

2951 do j=js-1,jeq+1 ; do i=is-1,ieq+1

2952 grid_sp_h2 = (2.0*cs%dx2h(i,j)*cs%dy2h(i,j)) / (cs%dx2h(i,j)+cs%dy2h(i,j))

2953 grid_sp_h3 = grid_sp_h2*sqrt(grid_sp_h2)

2954 cs%grid_sp_h3(i,j) = grid_sp_h3

2955 if (cs%Smagorinsky_Ah) then

2956 cs%Biharm_const_xx(i,j) = smag_bi_const * (grid_sp_h2 * grid_sp_h2)

2957 if (cs%bound_Coriolis) then

2958 fmax = max(abs(g%CoriolisBu(i-1,j-1)), abs(g%CoriolisBu(i,j-1)), &

2959 abs(g%CoriolisBu(i-1,j)), abs(g%CoriolisBu(i,j)))

2960 cs%Biharm_const2_xx(i,j) = (grid_sp_h2 * grid_sp_h2 * grid_sp_h2) * &

2961 (fmax * boundcorconst)

2962 endif

2963 endif

2964 if (cs%Leith_Ah) then

2965 cs%biharm6_const_xx(i,j) = leith_bi_const * (grid_sp_h3 * grid_sp_h3)

2966 endif

2967 if (cs%use_Leithy) then

2968 cs%biharm6_const_xx(i,j) = leith_bi_const * max(g%dxT(i,j),g%dyT(i,j))**6

2969 cs%m_const_leithy(i,j) = 0.5 * sqrt(cs%c_K) * max(g%dxT(i,j),g%dyT(i,j))

2970 cs%m_leithy_max(i,j) = 4. / max(g%dxT(i,j),g%dyT(i,j))**2

2971 endif

2972 cs%Ah_bg_xx(i,j) = max(ah, ah_vel_scale * grid_sp_h2 * sqrt(grid_sp_h2))

2973 if (cs%Re_Ah > 0.0) cs%Re_Ah_const_xx(i,j) = grid_sp_h3 / cs%Re_Ah

2974 if (ah_time_scale > 0.) cs%Ah_bg_xx(i,j) = &

2975 max(cs%Ah_bg_xx(i,j), (grid_sp_h2 * grid_sp_h2) / ah_time_scale)

2976 min_grid_sp_h4 = min(grid_sp_h2**2, min_grid_sp_h4)

2977 enddo ; enddo

2978 call min_across_pes(min_grid_sp_h4)

2979

2980 do j=js-1,jeq ; do i=is-1,ieq

2981 grid_sp_q2 = (2.0*cs%dx2q(i,j)*cs%dy2q(i,j)) / (cs%dx2q(i,j)+cs%dy2q(i,j))

2982 grid_sp_q3 = grid_sp_q2*sqrt(grid_sp_q2)

2983 if (cs%Smagorinsky_Ah) then

2984 cs%Biharm_const_xy(i,j) = smag_bi_const * (grid_sp_q2 * grid_sp_q2)

2985 if (cs%bound_Coriolis) then

2986 cs%Biharm_const2_xy(i,j) = (grid_sp_q2 * grid_sp_q2 * grid_sp_q2) * &

2987 (abs(g%CoriolisBu(i,j)) * boundcorconst)

2988 endif

2989 endif

2990 if ((cs%Leith_Ah) .or. (cs%use_Leithy))then

2991 cs%biharm6_const_xy(i,j) = leith_bi_const * (grid_sp_q3 * grid_sp_q3)

2992 endif

2993 cs%Ah_bg_xy(i,j) = max(ah, ah_vel_scale * grid_sp_q2 * sqrt(grid_sp_q2))

2994 if (cs%Re_Ah > 0.0) cs%Re_Ah_const_xy(i,j) = grid_sp_q3 / cs%Re_Ah

2995 if (ah_time_scale > 0.) cs%Ah_bg_xy(i,j) = &

2996 max(cs%Ah_bg_xy(i,j), (grid_sp_q2 * grid_sp_q2) / ah_time_scale)

2997 enddo ; enddo

2998 endif

2999 ! The Laplacian bounds should avoid overshoots when CS%bound_coef < 1.

3000 if (cs%Laplacian .and. cs%bound_Kh) then

3001 do j=js-1,jeq+1 ; do i=is-1,ieq+1

3002 denom = max( &

3003 (cs%dy2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j)) * &

3004 max(g%IdyCu(i,j)*g%IareaCu(i,j), g%IdyCu(i-1,j)*g%IareaCu(i-1,j)) ), &

3005 (cs%dx2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1)) * &

3006 max(g%IdxCv(i,j)*g%IareaCv(i,j), g%IdxCv(i,j-1)*g%IareaCv(i,j-1)) ) )

3007 cs%Kh_Max_xx(i,j) = 0.0

3008 if (denom > 0.0) &

3009 cs%Kh_Max_xx(i,j) = cs%bound_coef * 0.25 * idt / denom

3010 enddo ; enddo

3011 do j=js-1,jeq ; do i=is-1,ieq

3012 denom = max( &

3013 (cs%dx2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j)) * &

3014 max(g%IdxCu(i,j)*g%IareaCu(i,j), g%IdxCu(i,j+1)*g%IareaCu(i,j+1)) ), &

3015 (cs%dy2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j)) * &

3016 max(g%IdyCv(i,j)*g%IareaCv(i,j), g%IdyCv(i+1,j)*g%IareaCv(i+1,j)) ) )

3017 cs%Kh_Max_xy(i,j) = 0.0

3018 if (denom > 0.0) &

3019 cs%Kh_Max_xy(i,j) = cs%bound_coef * 0.25 * idt / denom

3020 enddo ; enddo

3021 if (cs%debug) then

3022 call hchksum(cs%Kh_Max_xx, "Kh_Max_xx", g%HI, haloshift=0, unscale=us%L_to_m**2*us%s_to_T)

3023 call bchksum(cs%Kh_Max_xy, "Kh_Max_xy", g%HI, haloshift=0, unscale=us%L_to_m**2*us%s_to_T)

3024 endif

3025 endif

3026 ! The biharmonic bounds should avoid overshoots when CS%bound_coef < 0.5, but

3027 ! empirically work for CS%bound_coef <~ 1.0

3028 if (cs%biharmonic .and. cs%bound_Ah) then

3029 do j=js-1,jeq+1 ; do i=is-2,ieq+1

3030 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))) + &

3031 (cs%dy2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j))) )) + &

3032 (cs%Idx2dyCu(i,j)*((cs%dx2q(i,j) * cs%DX_dyBu(i,j) * (g%IdxCu(i,j+1) + g%IdxCu(i,j))) + &

3033 (cs%dx2q(i,j-1)*cs%DX_dyBu(i,j-1)*(g%IdxCu(i,j) + g%IdxCu(i,j-1))) )) )

3034 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))) + &

3035 (cs%dy2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1))) )) + &

3036 (cs%Idx2dyCu(i,j)*((cs%dx2q(i,j) * cs%DY_dxBu(i,j) * (g%IdyCv(i+1,j) + g%IdyCv(i,j))) + &

3037 (cs%dx2q(i,j-1)*cs%DY_dxBu(i,j-1)*(g%IdyCv(i+1,j-1) + g%IdyCv(i,j-1))) )) )

3038 enddo ; enddo

3039 do j=js-2,jeq+1 ; do i=is-1,ieq+1

3040 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))) + &

3041 (cs%dy2q(i-1,j)*cs%DX_dyBu(i-1,j)*(g%IdxCu(i-1,j+1) + g%IdxCu(i-1,j))) )) + &

3042 (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))) + &

3043 (cs%dx2h(i,j) * cs%DY_dxT(i,j) * (g%IdyCu(i,j) + g%IdyCu(i-1,j))) ) ))

3044 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))) + &

3045 (cs%dy2q(i-1,j)*cs%DY_dxBu(i-1,j)*(g%IdyCv(i,j) + g%IdyCv(i-1,j))) )) + &

3046 (cs%Idx2dyCv(i,j)*((cs%dx2h(i,j+1)*cs%DX_dyT(i,j+1)*(g%IdxCv(i,j+1) + g%IdxCv(i,j))) + &

3047 (cs%dx2h(i,j) * cs%DX_dyT(i,j) * (g%IdxCv(i,j) + g%IdxCv(i,j-1))) )) )

3048 enddo ; enddo

3049 do j=js-1,jeq+1 ; do i=is-1,ieq+1

3050 denom = max( &

3051 (cs%dy2h(i,j) * &

3052 ((cs%DY_dxT(i,j)*((g%IdyCu(i,j)*u0u(i,j)) + (g%IdyCu(i-1,j)*u0u(i-1,j)))) + &

3053 (cs%DX_dyT(i,j)*((g%IdxCv(i,j)*v0u(i,j)) + (g%IdxCv(i,j-1)*v0u(i,j-1))))) * &

3054 max(g%IdyCu(i,j)*g%IareaCu(i,j), g%IdyCu(i-1,j)*g%IareaCu(i-1,j)) ), &

3055 (cs%dx2h(i,j) * &

3056 ((cs%DY_dxT(i,j)*((g%IdyCu(i,j)*u0v(i,j)) + (g%IdyCu(i-1,j)*u0v(i-1,j)))) + &

3057 (cs%DX_dyT(i,j)*((g%IdxCv(i,j)*v0v(i,j)) + (g%IdxCv(i,j-1)*v0v(i,j-1))))) * &

3058 max(g%IdxCv(i,j)*g%IareaCv(i,j), g%IdxCv(i,j-1)*g%IareaCv(i,j-1)) ) )

3059 cs%Ah_Max_xx(i,j) = 0.0

3060 if (denom > 0.0) then

3061 cs%Ah_Max_xx(i,j) = cs%bound_coef * 0.5 * idt / denom

3062 if (cs%EY24_EBT_BS) then

3063 cs%Ah_Max_xx_KS(i,j) = cs%bound_coef * 0.5 / (cs%KS_timescale * denom)

3064 endif

3065 endif

3066

3067 enddo ; enddo

3068 do j=js-1,jeq ; do i=is-1,ieq

3069 denom = max( &

3070 (cs%dx2q(i,j) * &

3071 ((cs%DX_dyBu(i,j)*((u0u(i,j+1)*g%IdxCu(i,j+1)) + (u0u(i,j)*g%IdxCu(i,j)))) + &

3072 (cs%DY_dxBu(i,j)*((v0u(i+1,j)*g%IdyCv(i+1,j)) + (v0u(i,j)*g%IdyCv(i,j))))) * &

3073 max(g%IdxCu(i,j)*g%IareaCu(i,j), g%IdxCu(i,j+1)*g%IareaCu(i,j+1)) ), &

3074 (cs%dy2q(i,j) * &

3075 ((cs%DX_dyBu(i,j)*((u0v(i,j+1)*g%IdxCu(i,j+1)) + (u0v(i,j)*g%IdxCu(i,j)))) + &

3076 (cs%DY_dxBu(i,j)*((v0v(i+1,j)*g%IdyCv(i+1,j)) + (v0v(i,j)*g%IdyCv(i,j))))) * &

3077 max(g%IdyCv(i,j)*g%IareaCv(i,j), g%IdyCv(i+1,j)*g%IareaCv(i+1,j)) ) )

3078 cs%Ah_Max_xy(i,j) = 0.0

3079 if (denom > 0.0) then

3080 cs%Ah_Max_xy(i,j) = cs%bound_coef * 0.5 * idt / denom

3081 if (cs%EY24_EBT_BS) then

3082 cs%Ah_Max_xy_KS(i,j) = cs%bound_coef * 0.5 / (cs%KS_timescale * denom)

3083 endif

3084 endif

3085

3086 enddo ; enddo

3087 if (cs%debug) then

3088 call hchksum(cs%Ah_Max_xx, "Ah_Max_xx", g%HI, haloshift=0, unscale=us%L_to_m**4*us%s_to_T)

3089 call bchksum(cs%Ah_Max_xy, "Ah_Max_xy", g%HI, haloshift=0, unscale=us%L_to_m**4*us%s_to_T)

3090 endif

3091 endif

3092 ! Register fields for output from this module.

3093 cs%id_normstress = register_diag_field('ocean_model', 'NoSt', diag%axesTL, time, &

3094 'Normal Stress', 's-1', conversion=us%s_to_T)

3095 cs%id_shearstress = register_diag_field('ocean_model', 'ShSt', diag%axesBL, time, &

3096 'Shear Stress', 's-1', conversion=us%s_to_T)

3097 cs%id_diffu = register_diag_field('ocean_model', 'diffu', diag%axesCuL, time, &

3098 'Zonal Acceleration from Horizontal Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)

3099 cs%id_diffv = register_diag_field('ocean_model', 'diffv', diag%axesCvL, time, &

3100 'Meridional Acceleration from Horizontal Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)

3101

3102 !CS%id_hf_diffu = register_diag_field('ocean_model', 'hf_diffu', diag%axesCuL, Time, &

3103 ! 'Fractional Thickness-weighted Zonal Acceleration from Horizontal Viscosity', &

3104 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)

3105 !if ((CS%id_hf_diffu > 0) .and. (present(ADp))) then

3106 ! call safe_alloc_alloc(CS%hf_diffu,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)

3107 ! call safe_alloc_ptr(ADp%diag_hfrac_u,G%IsdB,G%IedB,G%jsd,G%jed,GV%ke)

3108 !endif

3109

3110 !CS%id_hf_diffv = register_diag_field('ocean_model', 'hf_diffv', diag%axesCvL, Time, &

3111 ! 'Fractional Thickness-weighted Meridional Acceleration from Horizontal Viscosity', &

3112 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)

3113 !if ((CS%id_hf_diffv > 0) .and. (present(ADp))) then

3114 ! call safe_alloc_alloc(CS%hf_diffv,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)

3115 ! call safe_alloc_ptr(ADp%diag_hfrac_v,G%isd,G%ied,G%JsdB,G%JedB,GV%ke)

3116 !endif

3117

3118 cs%id_hf_diffu_2d = register_diag_field('ocean_model', 'hf_diffu_2d', diag%axesCu1, time, &

3119 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Horizontal Viscosity', &

3120 'm s-2', conversion=us%L_T2_to_m_s2)

3121 if ((cs%id_hf_diffu_2d > 0) .and. (present(adp))) then

3122 call safe_alloc_ptr(adp%diag_hfrac_u,g%IsdB,g%IedB,g%jsd,g%jed,gv%ke)

3123 endif

3124

3125 cs%id_hf_diffv_2d = register_diag_field('ocean_model', 'hf_diffv_2d', diag%axesCv1, time, &

3126 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Horizontal Viscosity', &

3127 'm s-2', conversion=us%L_T2_to_m_s2)

3128 if ((cs%id_hf_diffv_2d > 0) .and. (present(adp))) then

3129 call safe_alloc_ptr(adp%diag_hfrac_v,g%isd,g%ied,g%JsdB,g%JedB,gv%ke)

3130 endif

3131

3132 cs%id_h_diffu = register_diag_field('ocean_model', 'h_diffu', diag%axesCuL, time, &

3133 'Thickness Multiplied Zonal Acceleration from Horizontal Viscosity', &

3134 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

3135 if ((cs%id_h_diffu > 0) .and. (present(adp))) then

3136 call safe_alloc_ptr(adp%diag_hu,g%IsdB,g%IedB,g%jsd,g%jed,gv%ke)

3137 endif

3138

3139 cs%id_h_diffv = register_diag_field('ocean_model', 'h_diffv', diag%axesCvL, time, &

3140 'Thickness Multiplied Meridional Acceleration from Horizontal Viscosity', &

3141 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

3142 if ((cs%id_h_diffv > 0) .and. (present(adp))) then

3143 call safe_alloc_ptr(adp%diag_hv,g%isd,g%ied,g%JsdB,g%JedB,gv%ke)

3144 endif

3145

3146 cs%id_intz_diffu_2d = register_diag_field('ocean_model', 'intz_diffu_2d', diag%axesCu1, time, &

3147 'Depth-integral of Zonal Acceleration from Horizontal Viscosity', &

3148 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

3149 if ((cs%id_intz_diffu_2d > 0) .and. (present(adp))) then

3150 call safe_alloc_ptr(adp%diag_hu,g%IsdB,g%IedB,g%jsd,g%jed,gv%ke)

3151 endif

3152

3153 cs%id_intz_diffv_2d = register_diag_field('ocean_model', 'intz_diffv_2d', diag%axesCv1, time, &

3154 'Depth-integral of Meridional Acceleration from Horizontal Viscosity', &

3155 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

3156 if ((cs%id_intz_diffv_2d > 0) .and. (present(adp))) then

3157 call safe_alloc_ptr(adp%diag_hv,g%isd,g%ied,g%JsdB,g%JedB,gv%ke)

3158 endif

3159

3160 cs%id_diffu_visc_rem = register_diag_field('ocean_model', 'diffu_visc_rem', diag%axesCuL, time, &

3161 'Zonal Acceleration from Horizontal Viscosity multiplied by viscous remnant', &

3162 'm s-2', conversion=us%L_T2_to_m_s2)

3163 if ((cs%id_diffu_visc_rem > 0) .and. (present(adp))) then

3164 call safe_alloc_ptr(adp%visc_rem_u,g%IsdB,g%IedB,g%jsd,g%jed,gv%ke)

3165 endif

3166

3167 cs%id_diffv_visc_rem = register_diag_field('ocean_model', 'diffv_visc_rem', diag%axesCvL, time, &

3168 'Meridional Acceleration from Horizontal Viscosity multiplied by viscous remnant', &

3169 'm s-2', conversion=us%L_T2_to_m_s2)

3170 if ((cs%id_diffv_visc_rem > 0) .and. (present(adp))) then

3171 call safe_alloc_ptr(adp%visc_rem_v,g%isd,g%ied,g%JsdB,g%JedB,gv%ke)

3172 endif

3173

3174 if (cs%biharmonic) then

3175 cs%id_Ah_h = register_diag_field('ocean_model', 'Ahh', diag%axesTL, time, &

3176 'Biharmonic Horizontal Viscosity at h Points', 'm4 s-1', conversion=us%L_to_m**4*us%s_to_T, &

3177 cmor_field_name='difmxybo', &

3178 cmor_long_name='Ocean lateral biharmonic viscosity', &

3179 cmor_standard_name='ocean_momentum_xy_biharmonic_diffusivity')

3180 cs%id_Ah_q = register_diag_field('ocean_model', 'Ahq', diag%axesBL, time, &

3181 'Biharmonic Horizontal Viscosity at q Points', 'm4 s-1', conversion=us%L_to_m**4*us%s_to_T)

3182 cs%id_grid_Re_Ah = register_diag_field('ocean_model', 'grid_Re_Ah', diag%axesTL, time, &

3183 'Grid Reynolds number for the Biharmonic horizontal viscosity at h points', 'nondim')

3184 if (cs%EY24_EBT_BS) then

3185 cs%id_visc_limit_h_flag = register_diag_field('ocean_model', 'visc_limit_h_flag', diag%axesTL, time, &

3186 'Locations where the biharmonic viscosity reached the better_bound limiter at h points', 'nondim')

3187 cs%id_visc_limit_q_flag = register_diag_field('ocean_model', 'visc_limit_q_flag', diag%axesBL, time, &

3188 'Locations where the biharmonic viscosity reached the better_bound limiter at q points', 'nondim')

3189 cs%id_visc_limit_h = register_diag_field('ocean_model', 'visc_limit_h', diag%axesTL, time, &

3190 'Value of the biharmonic viscosity limiter at h points', 'nondim')

3191 cs%id_visc_limit_q = register_diag_field('ocean_model', 'visc_limit_q', diag%axesBL, time, &

3192 'Value of the biharmonic viscosity limiter at q points', 'nondim')

3193 cs%id_visc_limit_h_frac = register_diag_field('ocean_model', 'visc_limit_h_frac', diag%axesTL, time, &

3194 'Value of the biharmonic viscosity limiter at h points', 'nondim')

3195 cs%id_visc_limit_q_frac = register_diag_field('ocean_model', 'visc_limit_q_frac', diag%axesBL, time, &

3196 'Value of the biharmonic viscosity limiter at q points', 'nondim')

3197 endif

3198

3199 if (cs%id_grid_Re_Ah > 0) &

3200 ! Compute the smallest biharmonic viscosity capable of modifying the

3201 ! velocity at floating point precision.

3202 cs%min_grid_Ah = spacing(1.) * min_grid_sp_h4 * idt

3203 endif

3204 if (cs%Laplacian) then

3205 cs%id_Kh_h = register_diag_field('ocean_model', 'Khh', diag%axesTL, time, &

3206 'Laplacian Horizontal Viscosity at h Points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T, &

3207 cmor_field_name='difmxylo', &

3208 cmor_long_name='Ocean lateral Laplacian viscosity', &

3209 cmor_standard_name='ocean_momentum_xy_laplacian_diffusivity')

3210 cs%id_Kh_q = register_diag_field('ocean_model', 'Khq', diag%axesBL, time, &

3211 'Laplacian Horizontal Viscosity at q Points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)

3212 cs%id_grid_Re_Kh = register_diag_field('ocean_model', 'grid_Re_Kh', diag%axesTL, time, &

3213 'Grid Reynolds number for the Laplacian horizontal viscosity at h points', 'nondim')

3214 cs%id_vort_xy_q = register_diag_field('ocean_model', 'vort_xy_q', diag%axesBL, time, &

3215 'Vertical vorticity at q Points', 's-1', conversion=us%s_to_T)

3216 cs%id_div_xx_h = register_diag_field('ocean_model', 'div_xx_h', diag%axesTL, time, &

3217 'Horizontal divergence at h Points', 's-1', conversion=us%s_to_T)

3218 cs%id_sh_xy_q = register_diag_field('ocean_model', 'sh_xy_q', diag%axesBL, time, &

3219 'Shearing strain at q Points', 's-1', conversion=us%s_to_T)

3220 cs%id_sh_xx_h = register_diag_field('ocean_model', 'sh_xx_h', diag%axesTL, time, &

3221 'Horizontal tension at h Points', 's-1', conversion=us%s_to_T)

3222

3223 if (cs%id_grid_Re_Kh > 0) &

3224 ! Compute a smallest Laplacian viscosity capable of modifying the

3225 ! velocity at floating point precision.

3226 cs%min_grid_Kh = spacing(1.) * min_grid_sp_h2 * idt

3227 endif

3228 if (cs%use_GME) then

3229 cs%id_dudx_bt = register_diag_field('ocean_model', 'dudx_bt', diag%axesT1, time, &

3230 'Zonal component of the barotropic shearing strain at h points', 's-1', &

3231 conversion=us%s_to_T)

3232 cs%id_dudy_bt = register_diag_field('ocean_model', 'dudy_bt', diag%axesB1, time, &

3233 'Zonal component of the barotropic shearing strain at q points', 's-1', &

3234 conversion=us%s_to_T)

3235 cs%id_dvdy_bt = register_diag_field('ocean_model', 'dvdy_bt', diag%axesT1, time, &

3236 'Meridional component of the barotropic shearing strain at h points', 's-1', &

3237 conversion=us%s_to_T)

3238 cs%id_dvdx_bt = register_diag_field('ocean_model', 'dvdx_bt', diag%axesB1, time, &

3239 'Meridional component of the barotropic shearing strain at q points', 's-1', &

3240 conversion=us%s_to_T)

3241 cs%id_GME_coeff_h = register_diag_field('ocean_model', 'GME_coeff_h', diag%axesTL, time, &

3242 'GME coefficient at h Points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)

3243 cs%id_GME_coeff_q = register_diag_field('ocean_model', 'GME_coeff_q', diag%axesBL, time, &

3244 'GME coefficient at q Points', 'm2 s-1', conversion=us%L_to_m**2*us%s_to_T)

3245 cs%id_FrictWork_GME = register_diag_field('ocean_model','FrictWork_GME',diag%axesTL,time,&

3246 'Integral work done by lateral friction terms in GME (excluding diffusion of energy)', &

3247 'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2)

3248 endif

3249

3250 if (cs%EY24_EBT_BS) then

3251 cs%id_BS_coeff_h = register_diag_field('ocean_model', 'BS_coeff_h', diag%axesTL, time, &

3252 'Backscatter coefficient at h points', units='m2 s-1', conversion=us%L_to_m**2*us%s_to_T)

3253 cs%id_BS_coeff_q = register_diag_field('ocean_model', 'BS_coeff_q', diag%axesBL, time, &

3254 'Backscatter coefficient at q points', units='m2 s-1', conversion=us%L_to_m**2*us%s_to_T)

3255 endif

3256

3257 cs%id_FrictWork = register_diag_field('ocean_model','FrictWork',diag%axesTL,time,&

3258 'Integral work done by lateral friction terms. If GME is turned on, this '//&

3259 'includes the GME contribution.', &

3260 'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2)

3261 cs%id_FrictWorkIntz = register_diag_field('ocean_model','FrictWorkIntz',diag%axesT1,time, &

3262 'Depth integrated work done by lateral friction', &

3263 'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2, &

3264 cmor_field_name='dispkexyfo', &

3265 cmor_long_name='Depth integrated ocean kinetic energy dissipation due to lateral friction',&

3266 cmor_standard_name='ocean_kinetic_energy_dissipation_per_unit_area_due_to_xy_friction')

3267 cs%id_FrictWork_bh = register_diag_field('ocean_model','FrictWork_bh',diag%axesTL,time,&

3268 'Integral work done by the biharmonic lateral friction terms.', &

3269 'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2)

3270 cs%id_FrictWorkIntz_bh = register_diag_field('ocean_model','FrictWorkIntz_bh',diag%axesT1,time,&

3271 'Depth integrated work done by the biharmonic lateral friction', &

3272 'W m-2', conversion=us%RZ3_T3_to_W_m2*us%L_to_Z**2)

3273

3274end subroutine hor_visc_init

3275

3276!> hor_visc_vel_stencil returns the horizontal viscosity input velocity stencil size

3277function hor_visc_vel_stencil(CS) result(stencil)

3278 type(hor_visc_cs), intent(in) :: cs !< Control structure for horizontal viscosity

3279 integer :: stencil !< The horizontal viscosity velocity stencil size with the current settings.

3280

3281 stencil = 2

3282

3283 if ((cs%Leith_Kh) .or. (cs%Leith_Ah) .or. (cs%use_Leithy)) then

3284 stencil = 3

3285 endif

3286end function hor_visc_vel_stencil

3287

3288!> Calculates factors in the anisotropic orientation tensor to be align with the grid.

3289!! With n1=1 and n2=0, this recovers the approach of Large et al, 2001.

3290subroutine align_aniso_tensor_to_grid(CS, n1, n2)

3291 type(hor_visc_cs), intent(inout) :: CS !< Control structure for horizontal viscosity

3292 real, intent(in) :: n1 !< i-component of direction vector [nondim]

3293 real, intent(in) :: n2 !< j-component of direction vector [nondim]

3294 ! Local variables

3295 real :: recip_n2_norm ! The inverse of the squared magnitude of n1 and n2 [nondim]

3296 ! For normalizing n=(n1,n2) in case arguments are not a unit vector

3297 recip_n2_norm = (n1**2) + (n2**2)

3298 if (recip_n2_norm > 0.) recip_n2_norm = 1. / recip_n2_norm

3299 cs%n1n2_h(:,:) = 2. * ( n1 * n2 ) * recip_n2_norm

3300 cs%n1n2_q(:,:) = 2. * ( n1 * n2 ) * recip_n2_norm

3301 cs%n1n1_m_n2n2_h(:,:) = ( (n1 * n1) - (n2 * n2) ) * recip_n2_norm

3302 cs%n1n1_m_n2n2_q(:,:) = ( (n1 * n1) - (n2 * n2) ) * recip_n2_norm

3303end subroutine align_aniso_tensor_to_grid

3304

3305!> Apply a 1-1-4-1-1 Laplacian filter one time on GME diffusive flux to reduce any

3306!! horizontal two-grid-point noise

3307subroutine smooth_gme(CS, G, GME_flux_h, GME_flux_q)

3308 type(hor_visc_cs), intent(in) :: CS !< Control structure

3309 type(ocean_grid_type), intent(in) :: G !< Ocean grid

3310 real, dimension(SZI_(G),SZJ_(G)), optional, intent(inout) :: GME_flux_h!< GME diffusive flux

3311 !! at h points [L2 T-2 ~> m2 s-2]

3312 real, dimension(SZIB_(G),SZJB_(G)), optional, intent(inout) :: GME_flux_q!< GME diffusive flux

3313 !! at q points [L2 T-2 ~> m2 s-2]

3314 ! local variables

3315 real, dimension(SZI_(G),SZJ_(G)) :: GME_flux_h_original ! The previous value of GME_flux_h [L2 T-2 ~> m2 s-2]

3316 real, dimension(SZIB_(G),SZJB_(G)) :: GME_flux_q_original ! The previous value of GME_flux_q [L2 T-2 ~> m2 s-2]

3317 real :: wc, ww, we, wn, ws ! averaging weights for smoothing [nondim]

3318 integer :: i, j, s, halosz

3319 integer :: xh, xq ! The number of valid extra halo points for h and q points.

3320 integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq

3321

3322 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

3323 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

3324 xh = 0 ; xq = 0

3325

3326 do s=1,cs%num_smooth_gme

3327 if (present(gme_flux_h)) then

3328 if (xh < 0) then

3329 ! Update halos if needed, but avoid doing so more often than is needed.

3330 halosz = min(g%isc-g%isd, g%jsc-g%jsd, 2+cs%num_smooth_gme-s)

3331 call pass_var(gme_flux_h, g%Domain, halo=halosz)

3332 xh = halosz - 2

3333 endif

3334 gme_flux_h_original(:,:) = gme_flux_h(:,:)

3335 ! apply smoothing on GME

3336 do j=jsq-xh,jeq+1+xh ; do i=isq-xh,ieq+1+xh

3337 ! skip land points

3338 if (g%mask2dT(i,j)==0.) cycle

3339 ! compute weights

3340 ww = 0.125 * g%mask2dT(i-1,j)

3341 we = 0.125 * g%mask2dT(i+1,j)

3342 ws = 0.125 * g%mask2dT(i,j-1)

3343 wn = 0.125 * g%mask2dT(i,j+1)

3344 wc = 1.0 - ((ww+we)+(wn+ws))

3345 gme_flux_h(i,j) = wc * gme_flux_h_original(i,j) &

3346 + ((ww * gme_flux_h_original(i-1,j) + we * gme_flux_h_original(i+1,j)) &

3347 + (ws * gme_flux_h_original(i,j-1) + wn * gme_flux_h_original(i,j+1)))

3348 enddo ; enddo

3349 xh = xh - 1

3350 endif

3351 if (present(gme_flux_q)) then

3352 if (xq < 0) then

3353 ! Update halos if needed, but avoid doing so more often than is needed.

3354 halosz = min(g%isc-g%isd, g%jsc-g%jsd, 2+cs%num_smooth_gme-s)

3355 call pass_var(gme_flux_q, g%Domain, position=corner, complete=.true., halo=halosz)

3356 xq = halosz - 2

3357 endif

3358 gme_flux_q_original(:,:) = gme_flux_q(:,:)

3359 ! apply smoothing on GME

3360 do j=js-1-xq,je+xq ; do i=is-1-xq,ie+xq

3361 ! skip land points

3362 if (g%mask2dBu(i,j)==0.) cycle

3363 ! compute weights

3364 ww = 0.125 * g%mask2dBu(i-1,j)

3365 we = 0.125 * g%mask2dBu(i+1,j)

3366 ws = 0.125 * g%mask2dBu(i,j-1)

3367 wn = 0.125 * g%mask2dBu(i,j+1)

3368 wc = 1.0 - ((ww+we)+(wn+ws))

3369 gme_flux_q(i,j) = wc * gme_flux_q_original(i,j) &

3370 + ((ww * gme_flux_q_original(i-1,j) + we * gme_flux_q_original(i+1,j)) &

3371 + (ws * gme_flux_q_original(i,j-1) + wn * gme_flux_q_original(i,j+1)))

3372 enddo ; enddo

3373 xq = xq - 1

3374 endif

3375 enddo ! s-loop

3376end subroutine smooth_gme

3377

3378!> Apply a 9-point smoothing filter twice to a field staggered at a thickness point to reduce

3379!! horizontal two-grid-point noise.

3380!! Note that this subroutine does not conserve mass, so don't use it in situations where you

3381!! need conservation. Also note that it assumes that the input field has valid values in the

3382!! first two halo points upon entry.

3383subroutine smooth_x9_h(G, field_h, zero_land)

3384 type(ocean_grid_type), intent(in) :: G !< Ocean grid

3385 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: field_h !< h-point field to be smoothed [arbitrary]

3386 logical, optional, intent(in) :: zero_land !< If present and false, return the average

3387 !! of the surrounding ocean points when

3388 !! smoothing, otherwise use a value of 0 for

3389 !! land points and include them in the averages.

3390

3391 ! Local variables

3392 real :: fh_prev(SZI_(G),SZJ_(G)) ! The value of the h-point field at the previous iteration [arbitrary]

3393 real :: Iwts ! The inverse of the sum of the weights [nondim]

3394 logical :: zero_land_val ! The value of the zero_land optional argument or .true. if it is absent.

3395 integer :: i, j, s, is, ie, js, je

3396

3397 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

3398

3399 zero_land_val = .true. ; if (present(zero_land)) zero_land_val = zero_land

3400

3401 do s=1,0,-1

3402 fh_prev(:,:) = field_h(:,:)

3403 ! apply smoothing on field_h using rotationally symmetric expressions.

3404 do j=js-s,je+s ; do i=is-s,ie+s ; if (g%mask2dT(i,j) > 0.0) then

3405 iwts = 0.0625

3406 if (.not. zero_land_val) &

3407 iwts = 1.0 / ( (4.0*g%mask2dT(i,j) + &

3408 ( 2.0*((g%mask2dT(i-1,j) + g%mask2dT(i+1,j)) + &

3409 (g%mask2dT(i,j-1) + g%mask2dT(i,j+1))) + &

3410 ((g%mask2dT(i-1,j-1) + g%mask2dT(i+1,j+1)) + &

3411 (g%mask2dT(i-1,j+1) + g%mask2dT(i+1,j-1))) ) ) + 1.0e-16 )

3412 field_h(i,j) = iwts * ( 4.0*g%mask2dT(i,j) * fh_prev(i,j) &

3413 + (2.0*((g%mask2dT(i-1,j) * fh_prev(i-1,j) + g%mask2dT(i+1,j) * fh_prev(i+1,j)) + &

3414 (g%mask2dT(i,j-1) * fh_prev(i,j-1) + g%mask2dT(i,j+1) * fh_prev(i,j+1))) &

3415 + ((g%mask2dT(i-1,j-1) * fh_prev(i-1,j-1) + g%mask2dT(i+1,j+1) * fh_prev(i+1,j+1)) + &

3416 (g%mask2dT(i-1,j+1) * fh_prev(i-1,j+1) + g%mask2dT(i+1,j-1) * fh_prev(i-1,j-1))) ))

3417 endif ; enddo ; enddo

3418 enddo

3419

3420end subroutine smooth_x9_h

3421

3422!> Apply a 9-point smoothing filter twice to a pair of velocity components to reduce

3423!! horizontal two-grid-point noise.

3424!! Note that this subroutine does not conserve angular momentum, so don't use it

3425!! in situations where you need conservation. Also note that it assumes that the

3426!! input fields have valid values in the first two halo points upon entry.

3427subroutine smooth_x9_uv(G, field_u, field_v, zero_land)

3428 type(ocean_grid_type), intent(in) :: G !< Ocean grid

3429 real, dimension(SZIB_(G),SZJ_(G)), intent(inout) :: field_u !< u-point field to be smoothed [arbitrary]

3430 real, dimension(SZI_(G),SZJB_(G)), intent(inout) :: field_v !< v-point field to be smoothed [arbitrary]

3431 logical, optional, intent(in) :: zero_land !< If present and false, return the average

3432 !! of the surrounding ocean points when

3433 !! smoothing, otherwise use a value of 0 for

3434 !! land points and include them in the averages.

3435

3436 ! Local variables.

3437 real :: fu_prev(SZIB_(G),SZJ_(G)) ! The value of the u-point field at the previous iteration [arbitrary]

3438 real :: fv_prev(SZI_(G),SZJB_(G)) ! The value of the v-point field at the previous iteration [arbitrary]

3439 real :: Iwts ! The inverse of the sum of the weights [nondim]

3440 logical :: zero_land_val ! The value of the zero_land optional argument or .true. if it is absent.

3441 integer :: i, j, s, is, ie, js, je, Isq, Ieq, Jsq, Jeq

3442

3443 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

3444 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

3445

3446 zero_land_val = .true. ; if (present(zero_land)) zero_land_val = zero_land

3447

3448 do s=1,0,-1

3449 fu_prev(:,:) = field_u(:,:)

3450 ! apply smoothing on field_u using the original non-rotationally symmetric expressions.

3451 do j=js-s,je+s ; do i=isq-s,ieq+s ; if (g%mask2dCu(i,j) > 0.0) then

3452 iwts = 0.0625

3453 if (.not. zero_land_val) &

3454 iwts = 1.0 / ( (4.0*g%mask2dCu(i,j) + &

3455 ( 2.0*((g%mask2dCu(i-1,j) + g%mask2dCu(i+1,j)) + &

3456 (g%mask2dCu(i,j-1) + g%mask2dCu(i,j+1))) + &

3457 ((g%mask2dCu(i-1,j-1) + g%mask2dCu(i+1,j+1)) + &

3458 (g%mask2dCu(i-1,j+1) + g%mask2dCu(i+1,j-1))) ) ) + 1.0e-16 )

3459 field_u(i,j) = iwts * ( 4.0*g%mask2dCu(i,j) * fu_prev(i,j) &

3460 + (2.0*((g%mask2dCu(i-1,j) * fu_prev(i-1,j) + g%mask2dCu(i+1,j) * fu_prev(i+1,j)) + &

3461 (g%mask2dCu(i,j-1) * fu_prev(i,j-1) + g%mask2dCu(i,j+1) * fu_prev(i,j+1))) &

3462 + ((g%mask2dCu(i-1,j-1) * fu_prev(i-1,j-1) + g%mask2dCu(i+1,j+1) * fu_prev(i+1,j+1)) + &

3463 (g%mask2dCu(i-1,j+1) * fu_prev(i-1,j+1) + g%mask2dCu(i+1,j-1) * fu_prev(i-1,j-1))) ))

3464 endif ; enddo ; enddo

3465

3466 fv_prev(:,:) = field_v(:,:)

3467 ! apply smoothing on field_v using the original non-rotationally symmetric expressions.

3468 do j=jsq-s,jeq+s ; do i=is-s,ie+s ; if (g%mask2dCv(i,j) > 0.0) then

3469 iwts = 0.0625

3470 if (.not. zero_land_val) &

3471 iwts = 1.0 / ( (4.0*g%mask2dCv(i,j) + &

3472 ( 2.0*((g%mask2dCv(i-1,j) + g%mask2dCv(i+1,j)) + &

3473 (g%mask2dCv(i,j-1) + g%mask2dCv(i,j+1))) + &

3474 ((g%mask2dCv(i-1,j-1) + g%mask2dCv(i+1,j+1)) + &

3475 (g%mask2dCv(i-1,j+1) + g%mask2dCv(i+1,j-1))) ) ) + 1.0e-16 )

3476 field_v(i,j) = iwts * ( 4.0*g%mask2dCv(i,j) * fv_prev(i,j) &

3477 + (2.0*((g%mask2dCv(i-1,j) * fv_prev(i-1,j) + g%mask2dCv(i+1,j) * fv_prev(i+1,j)) + &

3478 (g%mask2dCv(i,j-1) * fv_prev(i,j-1) + g%mask2dCv(i,j+1) * fv_prev(i,j+1))) &

3479 + ((g%mask2dCv(i-1,j-1) * fv_prev(i-1,j-1) + g%mask2dCv(i+1,j+1) * fv_prev(i+1,j+1)) + &

3480 (g%mask2dCv(i-1,j+1) * fv_prev(i-1,j+1) + g%mask2dCv(i+1,j-1) * fv_prev(i-1,j-1))) ))

3481 endif ; enddo ; enddo

3482 enddo

3483

3484end subroutine smooth_x9_uv

3485

3486!> Deallocates any variables allocated in hor_visc_init.

3487subroutine hor_visc_end(CS)

3488 type(hor_visc_cs), intent(inout) :: cs !< Horizontal viscosity control structure

3489 if (cs%Laplacian .or. cs%biharmonic) then

3490 dealloc_(cs%dx2h) ; dealloc_(cs%dx2q) ; dealloc_(cs%dy2h) ; dealloc_(cs%dy2q)

3491 dealloc_(cs%dx_dyT) ; dealloc_(cs%dy_dxT) ; dealloc_(cs%dx_dyBu) ; dealloc_(cs%dy_dxBu)

3492 dealloc_(cs%reduction_xx) ; dealloc_(cs%reduction_xy)

3493 endif

3494 if (cs%Laplacian) then

3495 dealloc_(cs%Kh_bg_xx) ; dealloc_(cs%Kh_bg_xy)

3496 dealloc_(cs%grid_sp_h2)

3497 if (allocated(cs%Kh_bg_2d)) deallocate(cs%Kh_bg_2d)

3498

3499 if (allocated(cs%Kh_Max_xx)) deallocate(cs%Kh_Max_xx)

3500 if (allocated(cs%Kh_Max_xy)) deallocate(cs%Kh_Max_xy)

3501 if (allocated(cs%Laplac2_const_xx)) deallocate(cs%Laplac2_const_xx)

3502 if (allocated(cs%Laplac2_const_xy)) deallocate(cs%Laplac2_const_xy)

3503 if (allocated(cs%Laplac3_const_xx)) deallocate(cs%Laplac3_const_xx)

3504 if (allocated(cs%Laplac3_const_xy)) deallocate(cs%Laplac3_const_xy)

3505 endif

3506 if (cs%biharmonic) then

3507 dealloc_(cs%grid_sp_h3)

3508 dealloc_(cs%Idx2dyCu) ; dealloc_(cs%Idx2dyCv)

3509 dealloc_(cs%Idxdy2u) ; dealloc_(cs%Idxdy2v)

3510 dealloc_(cs%Ah_bg_xx) ; dealloc_(cs%Ah_bg_xy)

3511

3512 if (allocated(cs%Ah_Max_xx)) deallocate(cs%Ah_Max_xx)

3513 if (allocated(cs%Ah_Max_xy)) deallocate(cs%Ah_Max_xy)

3514 if (allocated(cs%Ah_Max_xx_KS)) deallocate(cs%Ah_Max_xx_KS)

3515 if (allocated(cs%Ah_Max_xy_KS)) deallocate(cs%Ah_Max_xy_KS)

3516 if (allocated(cs%Biharm_const_xx)) deallocate(cs%Biharm_const_xx)

3517 if (allocated(cs%Biharm_const_xy)) deallocate(cs%Biharm_const_xy)

3518 if (allocated(cs%Biharm_const2_xx)) deallocate(cs%Biharm_const2_xx)

3519 if (allocated(cs%Biharm_const2_xy)) deallocate(cs%Biharm_const2_xy)

3520 if (allocated(cs%Biharm6_const_xx)) deallocate(cs%Biharm6_const_xx)

3521 if (allocated(cs%Biharm6_const_xy)) deallocate(cs%Biharm6_const_xy)

3522 if (allocated(cs%m_const_leithy)) deallocate(cs%m_const_leithy)

3523 if (allocated(cs%m_leithy_max)) deallocate(cs%m_leithy_max)

3524 if (allocated(cs%Re_Ah_const_xx)) deallocate(cs%Re_Ah_const_xx)

3525 if (allocated(cs%Re_Ah_const_xy)) deallocate(cs%Re_Ah_const_xy)

3526 endif

3527

3528 if (allocated(cs%n1n2_h)) deallocate(cs%n1n2_h)

3529 if (allocated(cs%n1n2_q)) deallocate(cs%n1n2_q)

3530 if (allocated(cs%n1n1_m_n2n2_h)) deallocate(cs%n1n1_m_n2n2_h)

3531 if (allocated(cs%n1n1_m_n2n2_q)) deallocate(cs%n1n1_m_n2n2_q)

3532

3533 if (cs%use_ZB2020) then

3534 call zb2020_end(cs%ZB2020)

3535 endif

3536

3537end subroutine hor_visc_end

3538!> \namespace mom_hor_visc

3539!!

3540!! \section section_horizontal_viscosity Horizontal viscosity in MOM

3541!!

3542!! This module contains the subroutine horizontal_viscosity that calculates the

3543!! effects of horizontal viscosity, including parameterizations of the value of

3544!! the viscosity itself. Subroutine horizontal_viscosity calculates the acceleration due to

3545!! some combination of a biharmonic viscosity and a Laplacian viscosity. Either or

3546!! both may use a coefficient that depends on the shear and strain of the flow.

3547!! All metric terms are retained. The Laplacian is calculated as the divergence of

3548!! a stress tensor, using the form suggested by \cite Smagorinsky1993. The biharmonic

3549!! is calculated by twice applying the divergence of the stress tensor that is

3550!! used to calculate the Laplacian, but without the dependence on thickness in the

3551!! first pass. This form permits a variable viscosity, and indicates no

3552!! acceleration for either resting fluid or solid body rotation.

3553!!

3554!! The form of the viscous accelerations is discussed extensively in \cite griffies2000,

3555!! and the implementation here follows that discussion closely.

3556!! We use the notation of \cite Smith2003 with the exception that the

3557!! isotropic viscosity is \f$\kappa_h\f$.

3558!!

3559!! In general, the horizontal stress tensor can be written as

3560!! \f[

3561!! {\bf \sigma} =

3562!! \begin{pmatrix}

3563!! \frac{1}{2} \left( \sigma_D + \sigma_T \right) & \frac{1}{2} \sigma_S \\\\

3564!! \frac{1}{2} \sigma_S & \frac{1}{2} \left( \sigma_D - \sigma_T \right)

3565!! \end{pmatrix}

3566!! \f]

3567!! where \f$\sigma_D\f$, \f$\sigma_T\f$ and \f$\sigma_S\f$ are stresses associated with

3568!! invariant factors in the strain-rate tensor. For a Newtonian fluid, the stress

3569!! tensor is usually linearly related to the strain-rate tensor. The horizontal

3570!! strain-rate tensor is

3571!! \f[

3572!! \dot{\bf e} =

3573!! \begin{pmatrix}

3574!! \frac{1}{2} \left( \dot{e}_D + \dot{e}_T \right) & \frac{1}{2} \dot{e}_S \\\\

3575!! \frac{1}{2} \dot{e}_S & \frac{1}{2} \left( \dot{e}_D - \dot{e}_T \right)

3576!! \end{pmatrix}

3577!! \f]

3578!! where \f$\dot{e}_D = \partial_x u + \partial_y v\f$ is the horizontal divergence,

3579!! \f$\dot{e}_T = \partial_x u - \partial_y v\f$ is the horizontal tension, and

3580!! \f$\dot{e}_S = \partial_y u + \partial_x v\f$ is the horizontal shear strain.

3581!!

3582!! The trace of the stress tensor, \f$tr(\bf \sigma) = \sigma_D\f$, is usually

3583!! absorbed into the pressure and only the deviatoric stress tensor considered.

3584!! From here on, we drop \f$\sigma_D\f$. The trace of the strain tensor, \f$tr(\bf e) =

3585!! \dot{e}_D\f$ is non-zero for horizontally divergent flow but only enters the

3586!! stress tensor through \f$\sigma_D\f$ and so we will drop \f$\sigma_D\f$ from

3587!! calculations of the strain tensor in the code. Therefore the horizontal stress

3588!! tensor can be considered to be

3589!! \f[

3590!! {\bf \sigma} =

3591!! \begin{pmatrix}

3592!! \frac{1}{2} \sigma_T & \frac{1}{2} \sigma_S \\\\

3593!! \frac{1}{2} \sigma_S & - \frac{1}{2} \sigma_T

3594!! \end{pmatrix}

3595!! .\f]

3596!!

3597!! The stresses above are linearly related to the strain through a viscosity

3598!! coefficient, \f$\kappa_h\f$:

3599!! \f{eqnarray*}{

3600!! \sigma_T & = & 2 \kappa_h \dot{e}_T \\\\

3601!! \sigma_S & = & 2 \kappa_h \dot{e}_S

3602!! .

3603!! \f}

3604!!

3605!! The viscosity \f$\kappa_h\f$ may either be a constant or variable. For example,

3606!! \f$\kappa_h\f$ may vary with the shear, as proposed by \cite Smagorinsky1993.

3607!!

3608!! The accelerations resulting form the divergence of the stress tensor are

3609!! \f{eqnarray*}{

3610!! \hat{\bf x} \cdot \left( \nabla \cdot {\bf \sigma} \right)

3611!! & = &

3612!! \partial_x \left( \frac{1}{2} \sigma_T \right)

3613!! + \partial_y \left( \frac{1}{2} \sigma_S \right)

3614!! \\\\

3615!! & = &

3616!! \partial_x \left( \kappa_h \dot{e}_T \right)

3617!! + \partial_y \left( \kappa_h \dot{e}_S \right)

3618!! \\\\

3619!! \hat{\bf y} \cdot \left( \nabla \cdot {\bf \sigma} \right)

3620!! & = &

3621!! \partial_x \left( \frac{1}{2} \sigma_S \right)

3622!! + \partial_y \left( - \frac{1}{2} \sigma_T \right)

3623!! \\\\

3624!! & = &

3625!! \partial_x \left( \kappa_h \dot{e}_S \right)

3626!! + \partial_y \left( - \kappa_h \dot{e}_T \right)

3627!! .

3628!! \f}

3629!!

3630!! The form of the Laplacian viscosity in general coordinates is:

3631!! \f{eqnarray*}{

3632!! \hat{\bf x} \cdot \left( \nabla \cdot \sigma \right)

3633!! & = &

3634!! \frac{1}{h} \left[ \partial_x \left( \kappa_h h \dot{e}_T \right)

3635!! + \partial_y \left( \kappa_h h \dot{e}_S \right) \right]

3636!! \\\\

3637!! \hat{\bf y} \cdot \left( \nabla \cdot \sigma \right)

3638!! & = &

3639!! \frac{1}{h} \left[ \partial_x \left( \kappa_h h \dot{e}_S \right)

3640!! - \partial_y \left( \kappa_h h \dot{e}_T \right) \right]

3641!! .

3642!! \f}

3643!!

3644!! \subsection section_laplacian_viscosity_coefficient Laplacian viscosity coefficient

3645!!

3646!! The horizontal viscosity coefficient, \f$\kappa_h\f$, can have multiple components.

3647!! The isotropic components are:

3648!! - A uniform background component, \f$\kappa_{bg}\f$.

3649!! - A constant but spatially variable 2D map, \f$\kappa_{2d}(x,y)\f$.

3650!! - A ''MICOM'' viscosity, \f$U_\nu \Delta(x,y)\f$, which uses a constant

3651!! velocity scale, \f$U_\nu\f$ and a measure of the grid-spacing \f$\Delta(x,y)^2 =

3652!! \frac{2 \Delta x^2 \Delta y^2}{\Delta x^2 + \Delta y^2}\f$.

3653!! - A function of

3654!! latitude, \f$\kappa_{\phi}(x,y) = \kappa_{\pi/2} |\sin(\phi)|^n\f$.

3655!! - A dynamic Smagorinsky viscosity, \f$\kappa_{Sm}(x,y,t) = C_{Sm} \Delta^2 \sqrt{\dot{e}_T^2 + \dot{e}_S^2}\f$.

3656!! - A dynamic Leith viscosity, \f$\kappa_{Lth}(x,y,t) =

3657!! C_{Lth} \Delta^3 \sqrt{|\nabla \zeta|^2 + |\nabla \dot{e}_D|^2}\f$.

3658!!

3659!! A maximum stable viscosity, \f$\kappa_{max}(x,y)\f$ is calculated based on the

3660!! grid-spacing and time-step and used to clip calculated viscosities.

3661!!

3662!! The static components of \f$\kappa_h\f$ are first combined as follows:

3663!! \f[

3664!! \kappa_{static} = \min \left[ \max\left(

3665!! \kappa_{bg},

3666!! U_\nu \Delta(x,y),

3667!! \kappa_{2d}(x,y),

3668!! \kappa_\phi(x,y)

3669!! \right)

3670!! , \kappa_{max}(x,y) \right]

3671!! \f]

3672!! and stored in the module control structure as variables <code>Kh_bg_xx</code> and

3673!! <code>Kh_bg_xy</code> for the tension (h-points) and shear (q-points) components

3674!! respectively.

3675!!

3676!! The full viscosity includes the dynamic components as follows:

3677!! \f[

3678!! \kappa_h(x,y,t) = r(\Delta,L_d)

3679!! \max \left( \kappa_{static}, \kappa_{Sm}, \kappa_{Lth} \right)

3680!! \f]

3681!! where \f$r(\Delta,L_d)\f$ is a resolution function.

3682!!

3683!! The dynamic Smagorinsky and Leith viscosity schemes are exclusive with each

3684!! other.

3685!!

3686!! \subsection section_viscous_boundary_conditions Viscous boundary conditions

3687!!

3688!! Free slip boundary conditions have been coded, although no slip boundary

3689!! conditions can be used with the Laplacian viscosity based on the 2D land-sea

3690!! mask. For a western boundary, for example, the boundary conditions with the

3691!! biharmonic operator would be written as:

3692!! \f[

3693!! \partial_x v = 0 ; \partial_x^3 v = 0 ; u = 0 ; \partial_x^2 u = 0 ,

3694!! \f]

3695!! while for a Laplacian operator, they are simply

3696!! \f[

3697!! \partial_x v = 0 ; u = 0 .

3698!! \f]

3699!! These boundary conditions are largely dictated by the use of an Arakawa

3700!! C-grid and by the varying layer thickness.

3701!!

3702!! \subsection section_anisotropic_viscosity Anisotropic viscosity

3703!!

3704!! \cite Large2001 proposed enhancing viscosity in a particular direction and the

3705!! approach was generalized in \cite Smith2003. We use the second form of their

3706!! two coefficient anisotropic viscosity (section 4.3). We also replace their

3707!! \f$A^\prime\f$ and $D$ such that \f$2A^\prime = 2 \kappa_h + D\f$ and

3708!! \f$\kappa_a = D\f$ so that \f$\kappa_h\f$ can be considered the isotropic

3709!! viscosity and \f$\kappa_a=D\f$ can be consider the anisotropic viscosity. The

3710!! direction of anisotropy is defined by a unit vector \f$\hat{\bf

3711!! n}=(n_1,n_2)\f$.

3712!!

3713!! The contributions to the stress tensor are

3714!! \f[

3715!! \begin{pmatrix}

3716!! \sigma_T \\\\ \sigma_S

3717!! \end{pmatrix}

3718!! =

3719!! \left[

3720!! \begin{pmatrix}

3721!! 2 \kappa_h + \kappa_a & 0 \\\\

3722!! 0 & 2 \kappa_h

3723!! \end{pmatrix}

3724!! + 2 \kappa_a n_1 n_2

3725!! \begin{pmatrix}

3726!! - 2 n_1 n_2 & n_1^2 - n_2^2 \\\\

3727!! n_1^2 - n_2^2 & 2 n_1 n_2

3728!! \end{pmatrix}

3729!! \right]

3730!! \begin{pmatrix}

3731!! \dot{e}_T \\\\ \dot{e}_S

3732!! \end{pmatrix}

3733!! \f]

3734!! Dissipation of kinetic energy requires \f$\kappa_h \geq 0\f$ and \f$2 \kappa_h + \kappa_a \geq 0\f$.

3735!! Note that when anisotropy is aligned with the x-direction, \f$n_1 = \pm 1\f$, then

3736!! \f$n_2 = 0\f$ and the cross terms vanish. The accelerations in this aligned limit

3737!! with constant coefficients become

3738!! \f{eqnarray*}{

3739!! \hat{\bf x} \cdot \nabla \cdot {\bf \sigma}

3740!! & = &

3741!! \partial_x \left( \left( \kappa_h + \frac{1}{2} \kappa_a \right) \dot{e}_T \right)

3742!! + \partial_y \left( \kappa_h \dot{e}_S \right)

3743!! \\\\

3744!! & = &

3745!! \left( \kappa_h + \kappa_a \right) \partial_{xx} u

3746!! + \kappa_h \partial_{yy} u

3747!! - \frac{1}{2} \kappa_a \partial_x \left( \partial_x u + \partial_y v \right)

3748!! \\\\

3749!! \hat{\bf y} \cdot \nabla \cdot {\bf \sigma}

3750!! & = &

3751!! \partial_x \left( \kappa_h \dot{e}_S \right)

3752!! - \partial_y \left( \left( \kappa_h + \frac{1}{2} \kappa_a \right) \dot{e}_T \right)

3753!! \\\\

3754!! & = &

3755!! \kappa_h \partial_{xx} v

3756!! + \left( \kappa_h + \kappa_a \right) \partial_{yy} v

3757!! - \frac{1}{2} \kappa_a \partial_y \left( \partial_x u + \partial_y v \right)

3758!! \f}

3759!! which has contributions akin to a negative divergence damping (a divergence

3760!! enhancement?) but which is weaker than the enhanced tension terms by half.

3761!!

3762!! \subsection section_viscous_discretization Discretization

3763!!

3764!! The horizontal tension, \f$\dot{e}_T\f$, is stored in variable <code>sh_xx</code> and

3765!! discretized as

3766!! \f[

3767!! \dot{e}_T

3768!! = \frac{\Delta y}{\Delta x} \delta_i \left( \frac{1}{\Delta y} u \right)

3769!! - \frac{\Delta x}{\Delta y} \delta_j \left( \frac{1}{\Delta x} v \right)

3770!! .

3771!! \f]

3772!! The horizontal divergent strain, \f$\dot{e}_D\f$, is stored in variable

3773!! <code>div_xx</code> and discretized as

3774!! \f[

3775!! \dot{e}_D

3776!! = \frac{1}{h A} \left( \delta_i \left( \overline{h}^i \Delta y \, u \right)

3777!! + \delta_j \left( \overline{h}^j\Delta x \, v \right) \right)

3778!! .

3779!! \f]

3780!! Note that for expediency this is the exact discretization used in the

3781!! continuity equation.

3782!!

3783!! The horizontal shear strain, \f$\dot{e}_S\f$, is stored in variable <code>sh_xy</code>

3784!! and discretized as

3785!! \f[

3786!! \dot{e}_S = v_x + u_y

3787!! \f]

3788!! where

3789!! \f{align*}{

3790!! v_x &= \frac{\Delta y}{\Delta x} \delta_i \left( \frac{1}{\Delta y} v \right) \\\\

3791!! u_y &= \frac{\Delta x}{\Delta y} \delta_j \left( \frac{1}{\Delta x} u \right)

3792!! \f}

3793!! which are calculated separately so that no-slip or free-slip boundary

3794!! conditions can be applied to \f$v_x\f$ and \f$u_y\f$ where appropriate.

3795!!

3796!! The tendency for the x-component of the divergence of stress is stored in

3797!! variable <code>diffu</code> and discretized as

3798!! \f[

3799!! \hat{\bf x} \cdot \left( \nabla \cdot {\bf \sigma} \right) =

3800!! \frac{1}{A \overline{h}^i} \left(

3801!! \frac{1}{\Delta y} \delta_i \left( h \Delta y^2 \kappa_h \dot{e}_T \right) +

3802!! \frac{1}{\Delta x} \delta_j \left( \tilde{h}^{ij} \Delta x^2 \kappa_h \dot{e}_S \right)

3803!! \right)

3804!! .

3805!! \f]

3806!!

3807!! The tendency for the y-component of the divergence of stress is stored in

3808!! variable <code>diffv</code> and discretized as

3809!! \f[

3810!! \hat{\bf y} \cdot \left( \nabla \cdot {\bf \sigma} \right) =

3811!! \frac{1}{A \overline{h}^j} \left(

3812!! \frac{1}{\Delta y} \delta_i \left( \tilde{h}^{ij} \Delta y^2 A_M \dot{e}_S \right)

3813!! - \frac{1}{\Delta x} \delta_j \left( h \Delta x^2 A_M \dot{e}_T \right)

3814!! \right)

3815!! .

3816!! \f]

3817!!

3818!! \subsection section_viscous_refs References

3819!!

3820!! Griffies, S.M., and Hallberg, R.W., 2000: Biharmonic friction with a

3821!! Smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models.

3822!! Monthly Weather Review, 128(8), 2935-2946.

3823!! https://doi.org/10.1175/1520-0493(2000)128%3C2935:BFWASL%3E2.0.CO;2

3824!!

3825!! Large, W.G., Danabasoglu, G., McWilliams, J.C., Gent, P.R. and Bryan, F.O.,

3826!! 2001: Equatorial circulation of a global ocean climate model with

3827!! anisotropic horizontal viscosity.

3828!! Journal of Physical Oceanography, 31(2), pp.518-536.

3829!! https://doi.org/10.1175/1520-0485(2001)031%3C0518:ECOAGO%3E2.0.CO;2

3830!!

3831!! Smagorinsky, J., 1993: Some historical remarks on the use of nonlinear

3832!! viscosities. Large eddy simulation of complex engineering and geophysical

3833!! flows, 1, 69-106.

3834!!

3835!! Smith, R.D., and McWilliams, J.C., 2003: Anisotropic horizontal viscosity for

3836!! ocean models. Ocean Modelling, 5(2), 129-156.

3837!! https://doi.org/10.1016/S1463-5003(02)00016-1

3838end module mom_hor_visc