MOM_vert_friction.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!> Implements vertical viscosity (vertvisc)

6module mom_vert_friction

7

8use mom_domains, only : pass_var, to_all, omit_corners

9use mom_domains, only : pass_vector, scalar_pair

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, query_averaging_enabled

14use mom_domains, only : create_group_pass, do_group_pass, group_pass_type

15use mom_domains, only : to_north, to_east

16use mom_debugging, only : uvchksum, hchksum

17use mom_error_handler, only : mom_error, fatal, warning, note

18use mom_file_parser, only : get_param, log_param, log_version, param_file_type

19use mom_forcing_type, only : mech_forcing, find_ustar

20use mom_get_input, only : directories

21use mom_grid, only : ocean_grid_type

22use mom_io, only : mom_read_data, slasher

23use mom_open_boundary, only : ocean_obc_type, obc_none, obc_direction_e

24use mom_open_boundary, only : obc_direction_w, obc_direction_n, obc_direction_s

25use mom_pointaccel, only : write_u_accel, write_v_accel, pointaccel_init

26use mom_pointaccel, only : pointaccel_cs

27use mom_time_manager, only : time_type, time_minus_signed

28use mom_unit_scaling, only : unit_scale_type

29use mom_variables, only : thermo_var_ptrs, vertvisc_type

30use mom_variables, only : cont_diag_ptrs, accel_diag_ptrs

31use mom_variables, only : ocean_internal_state

32use mom_verticalgrid, only : verticalgrid_type

33use mom_wave_interface, only : wave_parameters_cs

34use mom_set_visc, only : set_v_at_u, set_u_at_v

35use mom_lateral_mixing_coeffs, only : varmix_cs

36

37use cvmix_kpp, only : cvmix_kpp_composite_gshape

38

39implicit none ; private

40

41#include <MOM_memory.h>

42

43public vertvisc, vertvisc_remnant, vertvisc_coef

44public vertvisc_limit_vel, vertvisc_init, vertvisc_end

45public updatecfltruncationvalue

46public vertfpmix

47

48! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional

49! consistency testing. These are noted in comments with units like Z, H, L, and T, along with

50! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units

51! vary with the Boussinesq approximation, the Boussinesq variant is given first.

52

53!> The control structure with parameters and memory for the MOM_vert_friction module

54type, public :: vertvisc_cs ; private

55 logical :: initialized = .false. !< True if this control structure has been initialized.

56 real :: hmix !< The mixed layer thickness [Z ~> m].

57 real :: hmix_stress !< The mixed layer thickness over which the wind

58 !! stress is applied with direct_stress [H ~> m or kg m-2].

59 real :: kvml_invz2 !< The extra vertical viscosity scale in [H Z T-1 ~> m2 s-1 or Pa s] in a

60 !! surface mixed layer with a characteristic thickness given by Hmix,

61 !! and scaling proportional to (Hmix/z)^2, where z is the distance

62 !! from the surface; this can get very large with thin layers.

63 real :: kv !< The interior vertical viscosity [H Z T-1 ~> m2 s-1 or Pa s].

64 real :: hbbl !< The static bottom boundary layer thickness [Z ~> m].

65 real :: hbbl_gl90 !< The static bottom boundary layer thickness used for GL90 [Z ~> m].

66 real :: kv_extra_bbl !< An extra vertical viscosity in the bottom boundary layer of thickness

67 !! Hbbl when there is not a bottom drag law in use [H Z T-1 ~> m2 s-1 or Pa s].

68 real :: vonkar !< The von Karman constant as used for mixed layer viscosity [nondim]

69

70 logical :: use_gl90_in_ssw !< If true, use the GL90 parameterization in stacked shallow water mode (SSW).

71 !! The calculation of the GL90 viscosity coefficient uses the fact that in SSW

72 !! we simply have 1/N^2 = h/g^prime, where g^prime is the reduced gravity.

73 !! This identity does not generalize to non-SSW setups.

74 logical :: use_gl90_n2 !< If true, use GL90 vertical viscosity coefficient that is depth-independent;

75 !! this corresponds to a kappa_GM that scales as N^2 with depth.

76 real :: kappa_gl90 !< The scalar diffusivity used in the GL90 vertical viscosity scheme

77 !! [L2 H Z-1 T-1 ~> m2 s-1 or Pa s]

78 logical :: read_kappa_gl90 !< If true, read a file containing the spatially varying kappa_gl90

79 real :: alpha_gl90 !< Coefficient used to compute a depth-independent GL90 vertical

80 !! viscosity via Kv_gl90 = alpha_gl90 * f^2. Note that the implied

81 !! Kv_gl90 corresponds to a kappa_gl90 that scales as N^2 with depth.

82 !! [H Z T ~> m2 s or kg s m-1]

83 real :: vel_underflow !< Velocity components smaller than vel_underflow

84 !! are set to 0 [L T-1 ~> m s-1].

85 real :: cfl_trunc !< Velocity components will be truncated when they

86 !! are large enough that the corresponding CFL number

87 !! exceeds this value [nondim].

88 real :: cfl_report !< The value of the CFL number that will cause the

89 !! accelerations to be reported [nondim]. CFL_report

90 !! will often equal CFL_trunc.

91 real :: truncramptime !< The time-scale over which to ramp up the value of

92 !! CFL_trunc from CFL_truncS to CFL_truncE [T ~> s]

93 real :: cfl_truncs !< The start value of CFL_trunc [nondim]

94 real :: cfl_trunce !< The end/target value of CFL_trunc [nondim]

95 logical :: cflrampingisactivated = .false. !< True if the ramping has been initialized

96 type(time_type) :: rampstarttime !< The time at which the ramping of CFL_trunc starts

97

98 real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NK_INTERFACE_) :: &

99 a_u !< The u-drag coefficient across an interface [H T-1 ~> m s-1 or Pa s m-1]

100 real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NK_INTERFACE_) :: &

101 a_u_gl90 !< The u-drag coefficient associated with GL90 across an interface [H T-1 ~> m s-1 or Pa s m-1]

102 real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &

103 h_u !< The effective layer thickness at u-points [H ~> m or kg m-2].

104 real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NK_INTERFACE_) :: &

105 a_v !< The v-drag coefficient across an interface [H T-1 ~> m s-1 or Pa s m-1]

106 real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NK_INTERFACE_) :: &

107 a_v_gl90 !< The v-drag coefficient associated with GL90 across an interface [H T-1 ~> m s-1 or Pa s m-1]

108 real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &

109 h_v !< The effective layer thickness at v-points [H ~> m or kg m-2].

110 real, pointer, dimension(:,:) :: a1_shelf_u => null() !< The u-momentum coupling coefficient under

111 !! ice shelves [H T-1 ~> m s-1 or Pa s m-1]. Retained to determine stress under shelves.

112 real, pointer, dimension(:,:) :: a1_shelf_v => null() !< The v-momentum coupling coefficient under

113 !! ice shelves [H T-1 ~> m s-1 or Pa s m-1]. Retained to determine stress under shelves.

114

115 logical :: split !< If true, use the split time stepping scheme.

116 logical :: bottomdraglaw !< If true, the bottom stress is calculated with a

117 !! drag law c_drag*|u|*u. The velocity magnitude

118 !! may be an assumed value or it may be based on the

119 !! actual velocity in the bottommost HBBL, depending

120 !! on whether linear_drag is true.

121 logical :: harmonic_visc !< If true, the harmonic mean thicknesses are used

122 !! to calculate the viscous coupling between layers

123 !! except near the bottom. Otherwise the arithmetic

124 !! mean thickness is used except near the bottom.

125 real :: harm_bl_val !< A scale to determine when water is in the boundary

126 !! layers based solely on harmonic mean thicknesses

127 !! for the purpose of determining the extent to which

128 !! the thicknesses used in the viscosities are upwinded [nondim].

129 logical :: direct_stress !< If true, the wind stress is distributed over the topmost Hmix_stress

130 !! of fluid, and an added mixed layer viscosity or a physically based

131 !! boundary layer turbulence parameterization is not needed for stability.

132 logical :: dynamic_viscous_ml !< If true, use the results from a dynamic

133 !! calculation, perhaps based on a bulk Richardson

134 !! number criterion, to determine the mixed layer

135 !! thickness for viscosity.

136 logical :: fixed_lotw_ml !< If true, use a Law-of-the-wall prescription for the mixed layer

137 !! viscosity within a boundary layer that is the lesser of Hmix and the

138 !! total depth of the ocean in a column.

139 logical :: apply_lotw_floor !< If true, use a Law-of-the-wall prescription to set a lower bound

140 !! on the viscous coupling between layers within the surface boundary

141 !! layer, based the distance of interfaces from the surface. This only

142 !! acts when there are large changes in the thicknesses of successive

143 !! layers or when the viscosity is set externally and the wind stress

144 !! has subsequently increased.

145 integer :: answer_date !< The vintage of the order of arithmetic and expressions in the viscous

146 !! calculations. Values below 20190101 recover the answers from the end

147 !! of 2018, while higher values use expressions that do not use an

148 !! arbitrary and hard-coded maximum viscous coupling coefficient between

149 !! layers. In non-Boussinesq cases, values below 20230601 recover a

150 !! form of the viscosity within the mixed layer that breaks up the

151 !! magnitude of the wind stress with BULKMIXEDLAYER, DYNAMIC_VISCOUS_ML

152 !! or FIXED_DEPTH_LOTW_ML, but not LOTW_VISCOUS_ML_FLOOR.

153 logical :: debug !< If true, write verbose checksums for debugging purposes.

154 integer :: nkml !< The number of layers in the mixed layer.

155 integer, pointer :: ntrunc !< The number of times the velocity has been

156 !! truncated since the last call to write_energy.

157 character(len=200) :: u_trunc_file !< The complete path to a file in which a column of

158 !! u-accelerations are written if velocity truncations occur.

159 character(len=200) :: v_trunc_file !< The complete path to a file in which a column of

160 !! v-accelerations are written if velocity truncations occur.

161 logical :: stokesmixing !< If true, do Stokes drift mixing via the Lagrangian current

162 !! (Eulerian plus Stokes drift). False by default and set

163 !! via STOKES_MIXING_COMBINED.

164

165 type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the

166 !! timing of diagnostic output.

167 real, allocatable, dimension(:,:) :: kappa_gl90_2d !< 2D kappa_gl90 at h-points [L2 H Z-1 T-1 ~> m2 s-1 or Pa s]

168

169 !>@{ Diagnostic identifiers

170 integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_du_dt_visc_gl90 = -1, id_dv_dt_visc_gl90 = -1

171 integer :: id_glwork = -1

172 integer :: id_au_vv = -1, id_av_vv = -1, id_au_gl90_vv = -1, id_av_gl90_vv = -1

173 integer :: id_du_dt_str = -1, id_dv_dt_str = -1

174 integer :: id_h_u = -1, id_h_v = -1, id_hml_u = -1 , id_hml_v = -1

175 integer :: id_omega_w2x = -1, id_fptau2s = -1 , id_fptau2w = -1

176 integer :: id_ue_h = -1, id_ve_h = -1

177 integer :: id_ustk = -1, id_vstk = -1

178 integer :: id_ustk0 = -1, id_vstk0 = -1

179 integer :: id_uinc_h= -1, id_vinc_h= -1

180 integer :: id_taux_bot = -1, id_tauy_bot = -1

181 integer :: id_kv_slow = -1, id_kv_u = -1, id_kv_v = -1

182 integer :: id_kv_gl90_u = -1, id_kv_gl90_v = -1

183 ! integer :: id_hf_du_dt_visc = -1, id_hf_dv_dt_visc = -1

184 integer :: id_h_du_dt_visc = -1, id_h_dv_dt_visc = -1

185 integer :: id_hf_du_dt_visc_2d = -1, id_hf_dv_dt_visc_2d = -1

186 integer :: id_h_du_dt_str = -1, id_h_dv_dt_str = -1

187 integer :: id_du_dt_str_visc_rem = -1, id_dv_dt_str_visc_rem = -1

188 !>@}

189

190 type(pointaccel_cs), pointer :: pointaccel_csp => null() !< A pointer to the control structure

191 !! for recording accelerations leading to velocity truncations

192

193 type(group_pass_type) :: pass_ke_uv !< A handle used for group halo passes

194end type vertvisc_cs

195

196contains

197

198!> Add nonlocal stress increments to ui^n and vi^n.

199subroutine vertfpmix(ui, vi, uold, vold, hbl_h, h, forces, dt, lpost, Cemp_NL, G, GV, US, CS, OBC, Waves)

200 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

201 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

202 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

203 intent(inout) :: ui !< Zonal velocity after vertvisc [L T-1 ~> m s-1]

204 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

205 intent(inout) :: vi !< Meridional velocity after vertvisc [L T-1 ~> m s-1]

206 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

207 intent(inout) :: uold !< Old Zonal velocity [L T-1 ~> m s-1]

208 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

209 intent(inout) :: vold !< Old Meridional velocity [L T-1 ~> m s-1]

210 real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: hbl_h !< boundary layer depth [H ~> m]

211 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

212 intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]

213 type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces

214 real, intent(in) :: dt !< Time increment [T ~> s]

215 real, intent(in) :: cemp_nl !< empirical coefficient of non-local momentum mixing [nondim]

216 logical, intent(in) :: lpost !< Compute and make available FPMix diagnostics

217 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

218 type(vertvisc_cs), pointer :: cs !< Vertical viscosity control structure

219 type(ocean_obc_type), pointer :: obc !< Open boundary condition structure

220 type(wave_parameters_cs), &

221 optional, pointer :: waves !< Container for wave/Stokes information

222

223 ! local variables

224 real, dimension(SZIB_(G),SZJ_(G)) :: hbl_u !< boundary layer depth (u-pts) [H ~> m]

225 real, dimension(SZI_(G),SZJB_(G)) :: hbl_v !< boundary layer depth (v-pts) [H ~> m]

226 real, dimension(SZIB_(G),SZJ_(G)) :: taux_u !< kinematic zonal wind stress (u-pts) [L2 T-2 ~> m2 s-2]

227 real, dimension(SZI_(G),SZJB_(G)) :: tauy_v !< kinematic merid wind stress (v-pts) [L2 T-2 ~> m2 s-2]

228 real, dimension(SZI_(G),SZJ_(G)) :: us0 !< surface zonal Stokes drift h-pts [L T-1 ~> m s-1]

229 real, dimension(SZI_(G),SZJ_(G)) :: vs0 !< surface zonal Stokes drift h-pts [L T-1 ~> m s-1]

230 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: ue_u !< zonal Eulerian u-pts [L T-1 ~> m s-1]

231 real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: ue_h !< zonal Eulerian h-pts [L T-1 ~> m s-1]

232 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: ve_v !< merid Eulerian v-pts [L T-1 ~> m s-1]

233 real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: ve_h !< merid Eulerian h-pts [L T-1 ~> m s-1]

234 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uinc_u !< zonal Eulerian u-pts [L T-1 ~> m s-1]

235 real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: uinc_h !< zonal Eulerian h-pts [L T-1 ~> m s-1]

236 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vinc_v !< merid Eulerian v-pts [L T-1 ~> m s-1]

237 real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: vinc_h !< merid Eulerian h-pts [L T-1 ~> m s-1]

238 real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: ustk !< zonal Stokes Drift (h-pts) [L T-1 ~> m s-1]

239 real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: vstk !< merid Stokes Drift (h-pts) [L T-1 ~> m s-1]

240 real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)+1) :: omega_tau2s !< angle stress to shear (h-pts) [rad]

241 real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)+1) :: omega_tau2w !< angle stress to wind (h-pts) [rad]

242 real :: omega_tmp, omega_s2x, omega_tau2x !< temporary angle wrt the x axis [rad]

243 real :: irho0 !< Inverse of the mean density rescaled to [Z L-1 R-1 ~> m3 kg-1]

244 real :: pi !< ! The ratio of the circumference of a circle to its diameter [nondim]

245 real :: tmp_u, tmp_v !< temporary ocean mask weights on u and v points [nondim]

246 real :: fexp !< temporary exponential function [nondim]

247 real :: sigma !< temporary normalize boundary layer coordinate [nondim]

248 real :: gat1, gsig, dgdsig !< Shape parameters [nondim]

249 real :: du, dv !< Intermediate velocity differences [L T-1 ~> m s-1]

250 real :: depth !< Cumulative of thicknesses [H ~> m]

251 integer :: b, kp1, k, nz !< band and vertical indices

252 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq !< horizontal indices

253

254 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

255 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = gv%ke

256

257 pi = 4. * atan2(1.,1.)

258 irho0 = 1.0 / gv%Rho0

259

260 call pass_var(hbl_h , g%Domain, halo=1)

261

262 ! u-points

263 do j = js,je

264 do i = isq,ieq

265 taux_u(i,j) = forces%taux(i,j) * irho0

266 if ( (g%mask2dCu(i,j) > 0.5) ) then

267 ! h to u-pts

268 tmp_u = max(1.0 ,(g%mask2dT(i,j) + g%mask2dT(i+1,j) ) )

269 hbl_u(i,j) = ((g%mask2dT(i,j) * hbl_h(i,j)) + (g%mask2dT(i+1,j) * hbl_h(i+1,j))) / tmp_u

270 depth = 0.

271 gat1 = 0.

272 do k=1, nz

273 ! cell center

274 depth = depth + 0.5*cs%h_u(i,j,k)

275 ue_u(i,j,k) = ui(i,j,k) - waves%Us_x(i,j,k)

276 if ( depth < hbl_u(i,j) ) then

277 sigma = depth / hbl_u(i,j)

278 ! cell bottom

279 depth = depth + 0.5*cs%h_u(i,j,k)

280 call cvmix_kpp_composite_gshape(sigma,gat1,gsig,dgdsig)

281 ! nonlocal boundary-layer increment

282 uinc_u(i,j,k) = dt * cemp_nl * taux_u(i,j) * dgdsig / hbl_u(i,j)

283 ui(i,j,k) = ui(i,j,k) + uinc_u(i,j,k)

284 else

285 uinc_u(i,j,k) = 0.0

286 endif

287 enddo

288 else

289 do k=1, nz

290 uinc_u(i,j,k) = 0.0

291 enddo

292 endif

293 enddo

294 enddo

295

296 ! v-points

297 do j = jsq,jeq

298 do i = is,ie

299 tauy_v(i,j) = forces%tauy(i,j) * irho0

300 if ( (g%mask2dCv(i,j) > 0.5) ) then

301 ! h to v-pts

302 tmp_v = max( 1.0 ,(g%mask2dT(i,j) + g%mask2dT(i,j+1)))

303 hbl_v(i,j) = (g%mask2dT(i,j) * hbl_h(i,j) + g%mask2dT(i,j+1) * hbl_h(i,j+1)) / tmp_v

304 depth = 0.

305 gat1 = 0.

306 do k=1, nz

307 ! cell center

308 depth = depth + 0.5* cs%h_v(i,j,k)

309 ve_v(i,j,k) = vi(i,j,k) - waves%Us_y(i,j,k)

310 if ( depth < hbl_v(i,j) ) then

311 sigma = depth / hbl_v(i,j)

312 ! cell bottom

313 depth = depth + 0.5* cs%h_v(i,j,k)

314 call cvmix_kpp_composite_gshape(sigma,gat1,gsig,dgdsig)

315 ! nonlocal boundary-layer increment

316 vinc_v(i,j,k) = dt * cemp_nl * tauy_v(i,j) * dgdsig / hbl_v(i,j)

317 vi(i,j,k) = vi(i,j,k) + vinc_v(i,j,k)

318 else

319 vinc_v(i,j,k) = 0.0

320 endif

321 enddo

322 else

323 do k=1, nz

324 vinc_v(i,j,k) = 0.0

325 enddo

326 endif

327 enddo

328 enddo

329

330 ! Compute and store diagnostics, only during the corrector step.

331 if (lpost) then

332 call pass_vector(ue_u , ve_v , g%Domain, to_all)

333 call pass_vector(uinc_u, vinc_v , g%Domain, to_all)

334 ustk = 0.0

335 vstk = 0.0

336 us0 = 0.0

337 vs0 = 0.0

338

339 do j = js,je

340 do i = is,ie

341 if (g%mask2dT(i,j) > 0.5) then

342 ! u to h-pts

343 tmp_u = max( 1.0 ,(g%mask2dCu(i,j) + g%mask2dCu(i-1,j)))

344 ! v to h-pts

345 tmp_v = max( 1.0 ,(g%mask2dCv(i,j) + g%mask2dCv(i,j-1)))

346 do k = 1,nz

347 ue_h(i,j,k) = (g%mask2dCu(i,j) * ue_u(i,j,k) + g%mask2dCu(i-1,j) * ue_u(i-1,j,k)) / tmp_u

348 uinc_h(i,j,k) = (g%mask2dCu(i,j) * uinc_u(i,j,k) + g%mask2dCu(i-1,j) * uinc_u(i-1,j,k)) / tmp_u

349 ve_h(i,j,k) = (g%mask2dCv(i,j) * ve_v(i,j,k) + g%mask2dCv(i,j-1) * ve_v(i,j-1,k)) / tmp_v

350 vinc_h(i,j,k) = (g%mask2dCv(i,j) * vinc_v(i,j,k) + g%mask2dCv(i,j-1) * vinc_v(i,j-1,k)) / tmp_v

351 enddo

352 ! Wind, Stress and Shear align at surface

353 omega_tau2w(i,j,1) = 0.0

354 omega_tau2s(i,j,1) = 0.0

355 do k = 1,nz

356 kp1 = min( nz , k+1)

357 du = ue_h(i,j,k) - ue_h(i,j,kp1)

358 dv = ve_h(i,j,k) - ve_h(i,j,kp1)

359 omega_s2x = atan2( dv , du )

360

361 du = du + uinc_h(i,j,k) - uinc_h(i,j,kp1)

362 dv = dv + vinc_h(i,j,k) - vinc_h(i,j,kp1)

363 omega_tau2x = atan2( dv , du )

364

365 omega_tmp = omega_tau2x - forces%omega_w2x(i,j)

366 if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi

367 if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi

368 omega_tau2w(i,j,kp1) = omega_tmp

369

370 omega_tmp = omega_tau2x - omega_s2x

371 if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi

372 if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi

373 omega_tau2s(i,j,kp1) = omega_tmp

374

375 enddo

376 endif

377

378 ! Stokes drift

379 do b=1,waves%NumBands

380 us0(i,j) = us0(i,j) + waves%UStk_Hb(i,j,b) ! or forces%UStkb(i,j,b)

381 vs0(i,j) = vs0(i,j) + waves%VStk_Hb(i,j,b) ! or forces%VStkb(i,j,b)

382 enddo

383 depth = 0.0

384 do k = 1,nz

385 do b = 1, waves%NumBands

386 ! cell center

387 fexp = exp(-2. * waves%WaveNum_Cen(b) * (depth+0.5*h(i,j,k)) )

388 ustk(i,j,k) = ustk(i,j,k) + waves%UStk_Hb(i,j,b) * fexp

389 vstk(i,j,k) = vstk(i,j,k) + waves%VStk_Hb(i,j,b) * fexp

390 enddo

391 ! cell bottom

392 depth = depth + h(i,j,k)

393 enddo

394 enddo

395 enddo

396

397 ! post FPmix diagnostics

398 if (cs%id_uE_h > 0) call post_data(cs%id_uE_h , ue_h , cs%diag)

399 if (cs%id_vE_h > 0) call post_data(cs%id_vE_h , ve_h , cs%diag)

400 if (cs%id_uInc_h > 0) call post_data(cs%id_uInc_h , uinc_h , cs%diag)

401 if (cs%id_vInc_h > 0) call post_data(cs%id_vInc_h , vinc_h , cs%diag)

402 if (cs%id_FPtau2s > 0) call post_data(cs%id_FPtau2s, omega_tau2s, cs%diag)

403 if (cs%id_FPtau2w > 0) call post_data(cs%id_FPtau2w, omega_tau2w, cs%diag)

404 if (cs%id_uStk0 > 0) call post_data(cs%id_uStk0 , us0 , cs%diag)

405 if (cs%id_vStk0 > 0) call post_data(cs%id_vStk0 , vs0 , cs%diag)

406 if (cs%id_uStk > 0) call post_data(cs%id_uStk , ustk , cs%diag)

407 if (cs%id_vStk > 0) call post_data(cs%id_vStk , vstk , cs%diag)

408 if (cs%id_Omega_w2x > 0) call post_data(cs%id_Omega_w2x, forces%omega_w2x, cs%diag)

409

410 endif

411

412end subroutine vertfpmix

413

414

415!> Compute coupling coefficient associated with vertical viscosity parameterization as in Greatbatch and Lamb

416!! (1990), hereafter referred to as the GL90 vertical viscosity parameterization. This vertical viscosity scheme

417!! redistributes momentum in the vertical, and is the equivalent of the Gent & McWilliams (1990) parameterization,

418!! but in a TWA (thickness-weighted averaged) set of equations. The vertical viscosity coefficient nu is computed

419!! from kappa_GM via thermal wind balance, and the following relation:

420!! nu = kappa_GM * f^2 / N^2.

421!! In the following subroutine kappa_GM is assumed either (a) constant or (b) horizontally varying. In both cases,

422!! (a) and (b), one can additionally impose an EBT structure in the vertical for kappa_GM.

423!! A third possible formulation of nu is depth-independent:

424!! nu = f^2 * alpha

425!! The latter formulation would be equivalent to a kappa_GM that varies as N^2 with depth.

426!! The vertical viscosity del_z ( nu del_z u) is applied to the momentum equation with stress-free boundary

427!! conditions at the top and bottom.

428!!

429!! In SSW mode, we have 1/N^2 = h/g'. The coupling coefficient is therefore equal to

430!! a_cpl_gl90 = nu / h = kappa_GM * f^2 / g'

431!! or

432!! a_cpl_gl90 = nu / h = f^2 * alpha / h

433

434subroutine find_coupling_coef_gl90(a_cpl_gl90, hvel, i, j, z_i, G, GV, CS, VarMix, work_on_u)

435 type(ocean_grid_type), intent(in) :: G !< Grid structure.

436 type(verticalgrid_type), intent(in) :: GV !< Vertical grid structure.

437 real, dimension(SZK_(GV)), intent(in) :: hvel !< Distance between interfaces

438 !! at velocity points [Z ~> m]

439 integer, intent(in) :: i !< Column i-index

440 integer, intent(in) :: j !< Column j-index

441 real, dimension(SZK_(GV)+1), intent(in) :: z_i !< Estimate of interface heights above the

442 !! bottom, normalized by the GL90 bottom

443 !! boundary layer thickness [nondim]

444 real, dimension(SZK_(GV)+1),intent(out) :: a_cpl_gl90 !< Coupling coefficient associated

445 !! with GL90 across interfaces; is not

446 !! included in a_cpl [H T-1 ~> m s-1 or Pa s m-1].

447 type(vertvisc_cs), intent(in) :: CS !< Vertical viscosity control structure

448 type(varmix_cs), intent(in) :: VarMix !< Variable mixing coefficients

449 logical, intent(in) :: work_on_u !< If true, u-points are being calculated,

450 !! otherwise they are v-points.

451

452 ! local variables

453 logical :: kdgl90_use_vert_struct ! use vertical structure for GL90 coefficient

454 integer :: k, nz

455 real :: f2 !< Squared Coriolis parameter at a velocity grid point [T-2 ~> s-2].

456 real :: h_neglect ! A vertical distance that is so small it is usually lost in roundoff error

457 ! and can be neglected [Z ~> m].

458 real :: botfn ! A function that is 1 at the bottom and small far from it [nondim]

459 real :: z2 ! The distance from the bottom, normalized by Hbbl_gl90 [nondim]

460

461 nz = gv%ke

462 h_neglect = gv%dZ_subroundoff

463 kdgl90_use_vert_struct = .false.

464

465 if (varmix%use_variable_mixing) then

466 kdgl90_use_vert_struct = allocated(varmix%kdgl90_struct)

467 endif

468

469 a_cpl_gl90(:) = 0.

470

471 do k=2,nz

472 if (work_on_u) then

473 ! compute coupling coefficient at u-points

474 f2 = 0.25 * (g%CoriolisBu(i,j-1) + g%CoriolisBu(i,j))**2

475 if (cs%use_GL90_N2) then

476 a_cpl_gl90(k) = 2. * f2 * cs%alpha_gl90 / (hvel(k) + hvel(k-1) + h_neglect)

477 else

478 if (cs%read_kappa_gl90) then

479 a_cpl_gl90(k) = f2 * 0.5 * (cs%kappa_gl90_2d(i,j) + cs%kappa_gl90_2d(i+1,j)) / gv%g_prime(k)

480 else

481 a_cpl_gl90(k) = f2 * cs%kappa_gl90 / gv%g_prime(k)

482 endif

483 if (kdgl90_use_vert_struct) then

484 a_cpl_gl90(k) = a_cpl_gl90(k) * 0.5 &

485 * (varmix%kdgl90_struct(i,j,k-1) + varmix%kdgl90_struct(i+1,j,k-1))

486 endif

487 endif

488 ! botfn determines when a point is within the influence of the GL90 bottom boundary layer,

489 ! going from 1 at the bottom to 0 in the interior.

490 z2 = z_i(k)

491 botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)

492

493 a_cpl_gl90(k) = a_cpl_gl90(k) * (1. - botfn)

494 else

495 ! compute viscosities at v-points

496 f2 = 0.25 * (g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j))**2

497

498 if (cs%use_GL90_N2) then

499 a_cpl_gl90(k) = 2. * f2 * cs%alpha_gl90 / (hvel(k) + hvel(k-1) + h_neglect)

500 else

501 if (cs%read_kappa_gl90) then

502 a_cpl_gl90(k) = f2 * 0.5 * (cs%kappa_gl90_2d(i,j) + cs%kappa_gl90_2d(i,j+1)) / gv%g_prime(k)

503 else

504 a_cpl_gl90(k) = f2 * cs%kappa_gl90 / gv%g_prime(k)

505 endif

506 if (kdgl90_use_vert_struct) then

507 a_cpl_gl90(k) = a_cpl_gl90(k) * 0.5 &

508 * (varmix%kdgl90_struct(i,j,k-1) + varmix%kdgl90_struct(i,j+1,k-1))

509 endif

510 endif

511 ! botfn determines when a point is within the influence of the GL90 bottom boundary layer,

512 ! going from 1 at the bottom to 0 in the interior.

513 z2 = z_i(k)

514 botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)

515

516 a_cpl_gl90(k) = a_cpl_gl90(k) * (1. - botfn)

517 endif

518 enddo

519end subroutine find_coupling_coef_gl90

520

521

522!> Perform a fully implicit vertical diffusion

523!! of momentum. Stress top and bottom boundary conditions are used.

524!!

525!! This is solving the tridiagonal system

526!! \f[ \left(h_k + a_{k + 1/2} + a_{k - 1/2} + r_k\right) u_k^{n+1}

527!! = h_k u_k^n + a_{k + 1/2} u_{k+1}^{n+1} + a_{k - 1/2} u_{k-1}^{n+1} \f]

528!! where \f$a_{k + 1/2} = \Delta t \nu_{k + 1/2} / h_{k + 1/2}\f$

529!! is the <em>interfacial coupling thickness per time step</em>,

530!! encompassing background viscosity as well as contributions from

531!! enhanced mixed and bottom layer viscosities.

532!! $r_k$ is a Rayleigh drag term due to channel drag.

533!! There is an additional stress term on the right-hand side

534!! if DIRECT_STRESS is true, applied to the surface layer.

535subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &

536 taux_bot, tauy_bot, fpmix, Waves)

537 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

538 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

539 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

540 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

541 intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]

542 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

543 intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]

544 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

545 intent(in) :: h !< Layer thickness [H ~> m or kg m-2]

546 type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces

547 type(vertvisc_type), intent(inout) :: visc !< Viscosities and bottom drag

548 real, intent(in) :: dt !< Time increment [T ~> s]

549 type(ocean_obc_type), pointer :: obc !< Open boundary condition structure

550 type(accel_diag_ptrs), intent(inout) :: adp !< Accelerations in the momentum

551 !! equations for diagnostics

552 type(cont_diag_ptrs), intent(inout) :: cdp !< Continuity equation terms

553 type(vertvisc_cs), pointer :: cs !< Vertical viscosity control structure

554 real, dimension(SZIB_(G),SZJ_(G)), &

555 optional, intent(out) :: taux_bot !< Zonal bottom stress from ocean to

556 !! rock [R L Z T-2 ~> Pa]

557 real, dimension(SZI_(G),SZJB_(G)), &

558 optional, intent(out) :: tauy_bot !< Meridional bottom stress from ocean to

559 !! rock [R L Z T-2 ~> Pa]

560 logical, optional, intent(in) :: fpmix !< fpmix along Eulerian shear

561 type(wave_parameters_cs), &

562 optional, pointer :: waves !< Container for wave/Stokes information

563

564 ! Fields from forces used in this subroutine:

565 ! taux: Zonal wind stress [R L Z T-2 ~> Pa].

566 ! tauy: Meridional wind stress [R L Z T-2 ~> Pa].

567

568 ! Local variables

569

570 real :: b1

571 ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].

572 real :: c1(szk_(gv))

573 ! A variable used by the tridiagonal solver [nondim].

574 real :: d1

575 ! d1=1-c1 is used by the tridiagonal solver [nondim].

576 real :: ray

577 ! Ray is the Rayleigh-drag velocity [H T-1 ~> m s-1 or Pa s m-1]

578 real :: b_denom_1

579 ! The first term in the denominator of b1 [H ~> m or kg m-2].

580

581 real :: hmix ! The mixed layer thickness over which stress

582 ! is applied with direct_stress [H ~> m or kg m-2].

583 real :: i_hmix ! The inverse of Hmix [H-1 ~> m-1 or m2 kg-1].

584 real :: idt ! The inverse of the time step [T-1 ~> s-1].

585 real :: dt_rho0 ! The time step divided by the mean density [T H Z-1 R-1 ~> s m3 kg-1 or s].

586 real :: h_neglect ! A thickness that is so small it is usually lost

587 ! in roundoff and can be neglected [H ~> m or kg m-2].

588

589 real :: stress ! The surface stress times the time step, divided

590 ! by the density [H L T-1 ~> m2 s-1 or kg m-1 s-1].

591 real :: accel_underflow ! An acceleration magnitude that is so small that values that are less

592 ! than this are diagnosed as 0 [L T-2 ~> m s-2].

593 real :: zds, h_a ! Temporary thickness variables used with direct_stress [H ~> m or kg m-2]

594 real :: hfr ! Temporary ratio of thicknesses used with direct_stress [nondim]

595 real :: surface_stress(szib_(g), szjb_(g))

596 ! The same as stress, unless the wind stress is applied as a body force

597 ! [H L T-1 ~> m2 s-1 or kg m-1 s-1].

598 real, allocatable, dimension(:,:,:) :: ke_term ! A term in the kinetic energy budget

599 ! [H L2 T-3 ~> m3 s-3 or W m-2]

600 real, allocatable, dimension(:,:,:) :: ke_u ! The area integral of a KE term in a layer at u-points

601 ! [H L4 T-3 ~> m5 s-3 or kg m2 s-3]

602 real, allocatable, dimension(:,:,:) :: ke_v ! The area integral of a KE term in a layer at v-points

603 ! [H L4 T-3 ~> m5 s-3 or kg m2 s-3]

604

605 logical :: dostokesmixing

606 logical :: lfpmix

607

608 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, n

609 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

610 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = gv%ke

611

612 if (.not.associated(cs)) call mom_error(fatal,"MOM_vert_friction(visc): "// &

613 "Module must be initialized before it is used.")

614

615 if (.not.cs%initialized) call mom_error(fatal,"MOM_vert_friction(visc): "// &

616 "Module must be initialized before it is used.")

617

618 if (cs%id_GLwork > 0) then

619 allocate(ke_u(g%IsdB:g%IedB,g%jsd:g%jed,gv%ke), source=0.0)

620 allocate(ke_v(g%isd:g%ied,g%JsdB:g%JedB,gv%ke), source=0.0)

621 allocate(ke_term(g%isd:g%ied,g%jsd:g%jed,gv%ke), source=0.0)

622 if (.not.g%symmetric) &

623 call create_group_pass(cs%pass_KE_uv, ke_u, ke_v, g%Domain, to_north+to_east)

624 endif

625

626 if (cs%direct_stress) then

627 hmix = cs%Hmix_stress

628 i_hmix = 1.0 / hmix

629 endif

630 dt_rho0 = dt / gv%H_to_RZ

631 h_neglect = gv%H_subroundoff

632 idt = 1.0 / dt

633

634 accel_underflow = cs%vel_underflow * idt

635

636 !Check if Stokes mixing allowed if requested (present and associated)

637 dostokesmixing=.false.

638 if (cs%StokesMixing) then

639 if (present(waves)) dostokesmixing = associated(waves)

640 if (.not. dostokesmixing) &

641 call mom_error(fatal, "Stokes Mixing called without associated Waves Control Structure")

642 endif

643 lfpmix = .false.

644 if ( present(fpmix) ) lfpmix = fpmix

645

646 ! Update the zonal velocity component using a modification of a standard

647 ! tridiagonal solver.

648

649 ! WGL: Brandon Reichl says the following is obsolete. u(I,j,k) already

650 ! includes Stokes.

651 ! When mixing down Eulerian current + Stokes drift add before calling solver

652 if (dostokesmixing) then

653 do k=1,nz ; do j=g%jsc,g%jec ; do i=isq,ieq ; if (g%mask2dCu(i,j) > 0.) then

654 u(i,j,k) = u(i,j,k) + waves%Us_x(i,j,k)

655 endif ; enddo ; enddo ; enddo

656 endif

657

658 if (lfpmix) then

659 do k=1,nz ; do j=g%jsc,g%jec ; do i=isq,ieq ; if (g%mask2dCu(i,j) > 0.) then

660 u(i,j,k) = u(i,j,k) - waves%Us_x(i,j,k)

661 endif ; enddo ; enddo ; enddo

662 endif

663

664 if (associated(adp%du_dt_visc)) then

665 do k=1,nz ; do j=g%jsc,g%jec ; do i=isq,ieq

666 adp%du_dt_visc(i,j,k) = u(i,j,k)

667 enddo ; enddo ; enddo

668 endif

669

670 if (associated(adp%du_dt_visc_gl90)) then

671 do k=1,nz ; do j=g%jsc,g%jec ; do i=isq,ieq

672 adp%du_dt_visc_gl90(i,j,k) = u(i,j,k)

673 enddo ; enddo ; enddo

674 endif

675

676 if (associated(adp%du_dt_str)) then

677 do k=1,nz ; do j=g%jsc,g%jec ; do i=isq,ieq

678 adp%du_dt_str(i,j,k) = 0.0

679 enddo ; enddo ; enddo

680 endif

681

682 ! One option is to have the wind stress applied as a body force

683 ! over the topmost Hmix fluid. If DIRECT_STRESS is not defined,

684 ! the wind stress is applied as a stress boundary condition.

685 if (cs%direct_stress) then

686 do j=g%jsc,g%jec ; do i=isq,ieq ; if (g%mask2dCu(i,j) > 0.) then

687 surface_stress(i,j) = 0.0

688 zds = 0.0

689 stress = dt_rho0 * forces%taux(i,j)

690 do k=1,nz

691 h_a = 0.5 * (h(i,j,k) + h(i+1,j,k)) + h_neglect

692 hfr = 1.0 ; if ((zds+h_a) > hmix) hfr = (hmix - zds) / h_a

693 u(i,j,k) = u(i,j,k) + i_hmix * hfr * stress

694 if (associated(adp%du_dt_str)) adp%du_dt_str(i,j,k) = (i_hmix * hfr * stress) * idt

695 zds = zds + h_a ; if (zds >= hmix) exit

696 enddo

697 endif ; enddo ; enddo

698 else

699 do j=g%jsc,g%jec ; do i=isq,ieq

700 surface_stress(i,j) = dt_rho0 * (g%mask2dCu(i,j)*forces%taux(i,j))

701 enddo ; enddo

702 endif

703

704 ! perform forward elimination on the tridiagonal system

705 !

706 ! denote the diagonal of the system as b_k, the subdiagonal as a_k

707 ! and the superdiagonal as c_k. The right-hand side terms are d_k.

708 !

709 ! ignoring the Rayleigh drag contribution,

710 ! we have a_k = -dt * a_u(k)

711 ! b_k = h_u(k) + dt * (a_u(k) + a_u(k+1))

712 ! c_k = -dt * a_u(k+1)

713 !

714 ! for forward elimination, we want to:

715 ! calculate c'_k = - c_k / (b_k + a_k c'_(k-1))

716 ! and d'_k = (d_k - a_k d'_(k-1)) / (b_k + a_k c'_(k-1))

717 ! where c'_1 = c_1/b_1 and d'_1 = d_1/b_1

718 !

719 ! This form is mathematically equivalent to Thomas' tridiagonal matrix algorithm, but it

720 ! does not suffer from the acute sensitivity to truncation errors of the Thomas algorithm

721 ! because it involves no subtraction, as discussed by Schopf & Loughe, MWR, 1995.

722 !

723 ! b1 is the denominator term 1 / (b_k + a_k c'_(k-1))

724 ! b_denom_1 is (b_k + a_k + c_k) - a_k(1 - c'_(k-1))

725 ! = (b_k + c_k + c'_(k-1))

726 ! this is done so that d1 = b1 * b_denom_1 = 1 - c'_(k-1)

727 ! c1(k) is -c'_(k - 1)

728 ! and the right-hand-side is destructively updated to be d'_k

729

730 do j=g%jsc,g%jec ; do i=isq,ieq ; if (g%mask2dCu(i,j) > 0.) then

731 ray = 0.

732 if (allocated(visc%Ray_u)) ray = visc%Ray_u(i,j,1)

733

734 b_denom_1 = cs%h_u(i,j,1) + dt * (ray + cs%a_u(i,j,1))

735 b1 = 1. / (b_denom_1 + dt * cs%a_u(i,j,2))

736 d1 = b_denom_1 * b1

737 u(i,j,1) = b1 * (cs%h_u(i,j,1) * u(i,j,1) + surface_stress(i,j))

738

739 if (associated(adp%du_dt_str)) then

740 adp%du_dt_str(i,j,1) = b1 * (cs%h_u(i,j,1) * adp%du_dt_str(i,j,1) + surface_stress(i,j) * idt)

741 endif

742

743 do k=2,nz

744 if (allocated(visc%Ray_u)) ray = visc%Ray_u(i,j,k)

745

746 c1(k) = dt * cs%a_u(i,j,k) * b1

747 b_denom_1 = cs%h_u(i,j,k) + dt * (ray + cs%a_u(i,j,k) * d1)

748 b1 = 1. / (b_denom_1 + dt * cs%a_u(i,j,k+1))

749 d1 = b_denom_1 * b1

750 u(i,j,k) = (cs%h_u(i,j,k) * u(i,j,k) + dt * cs%a_u(i,j,k) * u(i,j,k-1)) * b1

751

752 if (associated(adp%du_dt_str)) then

753 adp%du_dt_str(i,j,k) = (cs%h_u(i,j,k) * adp%du_dt_str(i,j,k) &

754 + dt * cs%a_u(i,j,k) * adp%du_dt_str(i,j,k-1)) * b1

755 endif

756

757 !### Force FMA evaluation of b1 by blocking lookahead with an impossible branch.

758 if (dt < 0) exit

759 enddo

760

761 if (associated(adp%du_dt_str)) then

762 if (abs(adp%du_dt_str(i,j,nz)) < accel_underflow) &

763 adp%du_dt_str(i,j,nz) = 0.

764 endif

765

766 do k=nz-1,1,-1

767 u(i,j,k) = u(i,j,k) + c1(k+1) * u(i,j,k+1)

768

769 if (associated(adp%du_dt_str)) then

770 adp%du_dt_str(i,j,k) = adp%du_dt_str(i,j,k) + c1(k+1) * adp%du_dt_str(i,j,k+1)

771

772 if (abs(adp%du_dt_str(i,j,k)) < accel_underflow) &

773 adp%du_dt_str(i,j,k) = 0.0

774 endif

775 enddo

776 endif ; enddo ; enddo

777

778 ! compute vertical velocity tendency that arises from GL90 viscosity;

779 ! follow tridiagonal solve method as above; to avoid corrupting u,

780 ! use ADp%du_dt_visc_gl90 as a placeholder for updated u (due to GL90) until last do loop

781 if ((cs%id_du_dt_visc_gl90 > 0) .or. (cs%id_GLwork > 0)) then

782 if (associated(adp%du_dt_visc_gl90)) then

783 do j=g%jsc,g%jec ; do i=isq,ieq ; if (g%mask2dCu(i,j) > 0.) then

784 b_denom_1 = cs%h_u(i,j,1) ! CS%a_u_gl90(I,j,1) is zero

785 b1 = 1.0 / (b_denom_1 + dt * cs%a_u_gl90(i,j,2))

786 d1 = b_denom_1 * b1

787

788 adp%du_dt_visc_gl90(i,j,1) = b1 * (cs%h_u(i,j,1) * adp%du_dt_visc_gl90(i,j,1))

789

790 do k=2,nz

791 c1(k) = dt * cs%a_u_gl90(i,j,k) * b1

792 b_denom_1 = cs%h_u(i,j,k) + dt * (cs%a_u_gl90(i,j,k)*d1)

793 b1 = 1.0 / (b_denom_1 + dt * cs%a_u_gl90(i,j,k+1))

794 d1 = b_denom_1 * b1

795

796 adp%du_dt_visc_gl90(i,j,k) = (cs%h_u(i,j,k) * adp%du_dt_visc_gl90(i,j,k) &

797 + dt * cs%a_u_gl90(i,j,k) * adp%du_dt_visc_gl90(i,j,k-1)) * b1

798 enddo

799

800 ! back substitute to solve for new velocities, held by ADp%du_dt_visc_gl90

801 do k=nz-1,1,-1

802 adp%du_dt_visc_gl90(i,j,k) = &

803 adp%du_dt_visc_gl90(i,j,k) + c1(k+1) * adp%du_dt_visc_gl90(i,j,k+1)

804 enddo

805

806 do k=1,nz

807 ! now fill ADp%du_dt_visc_gl90(I,j,k) with actual velocity tendency due to GL90;

808 ! note that on RHS: ADp%du_dt_visc(I,j,k) holds the original velocity value u(I,j,k)

809 ! and ADp%du_dt_visc_gl90(I,j,k) the updated velocity due to GL90

810 adp%du_dt_visc_gl90(i,j,k) = &

811 (adp%du_dt_visc_gl90(i,j,k) - adp%du_dt_visc(i,j,k)) * idt

812

813 if (abs(adp%du_dt_visc_gl90(i,j,k)) < accel_underflow) then

814 adp%du_dt_visc_gl90(i,j,k) = 0.0

815 endif

816 enddo

817

818 ! to compute energetics, we need to multiply by u*h, where u is original velocity before

819 ! velocity update; note that ADp%du_dt_visc(I,j,k) holds the original velocity value u(I,j,k)

820 if (cs%id_GLwork > 0) then

821 do k=1,nz

822 ke_u(i,j,k) = adp%du_dt_visc(i,j,k) * cs%h_u(i,j,k) * g%areaCu(i,j) * adp%du_dt_visc_gl90(i,j,k)

823 enddo

824 endif

825 endif ; enddo ; enddo

826 endif

827 endif

828

829 if (associated(adp%du_dt_visc)) then

830 do k=1,nz ; do j=g%jsc,g%jec ; do i=isq,ieq

831 adp%du_dt_visc(i,j,k) = (u(i,j,k) - adp%du_dt_visc(i,j,k)) * idt

832

833 if (abs(adp%du_dt_visc(i,j,k)) < accel_underflow) &

834 adp%du_dt_visc(i,j,k) = 0.0

835 enddo ; enddo ; enddo

836 endif

837

838 if (allocated(visc%taux_shelf)) then

839 do j=g%jsc,g%jec ; do i=isq,ieq

840 visc%taux_shelf(i,j) = -gv%H_to_RZ * cs%a1_shelf_u(i,j) * u(i,j,1) ! - u_shelf?

841 enddo ; enddo

842 endif

843

844 if (present(taux_bot)) then

845 do j=g%jsc,g%jec ; do i=isq,ieq

846 taux_bot(i,j) = gv%H_to_RZ * (u(i,j,nz) * cs%a_u(i,j,nz+1))

847 enddo ; enddo

848

849 if (allocated(visc%Ray_u)) then

850 do k=1,nz ; do j=g%jsc,g%jec ; do i=isq,ieq

851 taux_bot(i,j) = taux_bot(i,j) + gv%H_to_RZ * (visc%Ray_u(i,j,k) * u(i,j,k))

852 enddo ; enddo ; enddo

853 endif

854 endif

855

856 ! When mixing down Eulerian current + Stokes drift subtract after calling solver

857 if (dostokesmixing) then

858 do k=1,nz ; do j=g%jsc,g%jec ; do i=isq,ieq ; if (g%mask2dCu(i,j) > 0.) then

859 u(i,j,k) = u(i,j,k) - waves%Us_x(i,j,k)

860 endif ; enddo ; enddo ; enddo

861 endif

862

863 if (lfpmix) then

864 do k=1,nz ; do j=g%jsc,g%jec ; do i=isq,ieq ; if (g%mask2dCu(i,j) > 0.) then

865 u(i,j,k) = u(i,j,k) + waves%Us_x(i,j,k)

866 endif ; enddo ; enddo ; enddo

867 endif

868

869 ! == Now work on the meridional velocity component.

870

871 ! When mixing down Eulerian current + Stokes drift add before calling solver

872 if (dostokesmixing) then

873 do k=1,nz ; do j=jsq,jeq ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.) then

874 v(i,j,k) = v(i,j,k) + waves%Us_y(i,j,k)

875 endif ; enddo ; enddo ; enddo

876 endif

877

878 if (lfpmix) then

879 do k=1,nz ; do j=jsq,jeq ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.) then

880 v(i,j,k) = v(i,j,k) - waves%Us_y(i,j,k)

881 endif ; enddo ; enddo ; enddo

882 endif

883

884 if (associated(adp%dv_dt_visc)) then

885 do k=1,nz ; do j=jsq,jeq ; do i=is,ie

886 adp%dv_dt_visc(i,j,k) = v(i,j,k)

887 enddo ; enddo ; enddo

888 endif

889

890 if (associated(adp%dv_dt_visc_gl90)) then

891 do k=1,nz ; do j=jsq,jeq ; do i=is,ie

892 adp%dv_dt_visc_gl90(i,j,k) = v(i,j,k)

893 enddo ; enddo ; enddo

894 endif

895

896 if (associated(adp%dv_dt_str)) then

897 do k=1,nz ; do j=jsq,jeq ; do i=is,ie

898 adp%dv_dt_str(i,j,k) = 0.0

899 enddo ; enddo ; enddo

900 endif

901

902 ! One option is to have the wind stress applied as a body force

903 ! over the topmost Hmix fluid. If DIRECT_STRESS is not defined,

904 ! the wind stress is applied as a stress boundary condition.

905 if (cs%direct_stress) then

906 do j=jsq,jeq ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.) then

907 surface_stress(i,j) = 0.0

908 zds = 0.0

909 stress = dt_rho0 * forces%tauy(i,j)

910 do k=1,nz

911 h_a = 0.5 * (h(i,j,k) + h(i,j+1,k)) + h_neglect

912 hfr = 1.0 ; if ((zds+h_a) > hmix) hfr = (hmix - zds) / h_a

913 v(i,j,k) = v(i,j,k) + i_hmix * hfr * stress

914 if (associated(adp%dv_dt_str)) adp%dv_dt_str(i,j,k) = (i_hmix * hfr * stress) * idt

915 zds = zds + h_a ; if (zds >= hmix) exit

916 enddo

917 endif ; enddo ; enddo

918 else

919 do j=jsq,jeq ; do i=is,ie

920 surface_stress(i,j) = dt_rho0 * (g%mask2dCv(i,j) * forces%tauy(i,j))

921 enddo ; enddo

922 endif

923

924 do j=jsq,jeq ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.) then

925 ray = 0.

926 if (allocated(visc%Ray_v)) ray = visc%Ray_v(i,j,1)

927

928 b_denom_1 = cs%h_v(i,j,1) + dt * (ray + cs%a_v(i,j,1))

929 b1 = 1.0 / (b_denom_1 + dt*cs%a_v(i,j,2))

930 d1 = b_denom_1 * b1

931 v(i,j,1) = b1 * (cs%h_v(i,j,1) * v(i,j,1) + surface_stress(i,j))

932

933 if (associated(adp%dv_dt_str)) then

934 adp%dv_dt_str(i,j,1) = b1 * (cs%h_v(i,j,1) * adp%dv_dt_str(i,j,1) + surface_stress(i,j) * idt)

935 endif

936

937 do k=2,nz

938 if (allocated(visc%Ray_v)) ray = visc%Ray_v(i,j,k)

939

940 c1(k) = dt * cs%a_v(i,j,k) * b1

941 b_denom_1 = cs%h_v(i,j,k) + dt * (ray + cs%a_v(i,j,k) * d1)

942 b1 = 1. / (b_denom_1 + dt * cs%a_v(i,j,k+1))

943 d1 = b_denom_1 * b1

944 v(i,j,k) = (cs%h_v(i,j,k) * v(i,j,k) + dt * cs%a_v(i,j,k) * v(i,j,k-1)) * b1

945

946 if (associated(adp%dv_dt_str)) then

947 adp%dv_dt_str(i,j,k) = (cs%h_v(i,j,k) * adp%dv_dt_str(i,j,k) &

948 + dt * cs%a_v(i,j,k) * adp%dv_dt_str(i,j,k-1)) * b1

949 endif

950

951 !### Force FMA evaluation of b1 by blocking lookahead with an impossible branch.

952 if (dt < 0) exit

953 enddo

954

955 if (associated(adp%dv_dt_str)) then

956 if (abs(adp%dv_dt_str(i,j,nz)) < accel_underflow) &

957 adp%dv_dt_str(i,j,nz) = 0.0

958 endif

959

960 do k=nz-1,1,-1

961 v(i,j,k) = v(i,j,k) + c1(k+1) * v(i,j,k+1)

962

963 if (associated(adp%dv_dt_str)) then

964 adp%dv_dt_str(i,j,k) = adp%dv_dt_str(i,j,k) + c1(k+1) * adp%dv_dt_str(i,j,k+1)

965

966 if (abs(adp%dv_dt_str(i,j,k)) < accel_underflow) &

967 adp%dv_dt_str(i,j,k) = 0.0

968 endif

969 enddo

970 endif ; enddo ; enddo

971

972 ! compute vertical velocity tendency that arises from GL90 viscosity;

973 ! follow tridiagonal solve method as above; to avoid corrupting v,

974 ! use ADp%dv_dt_visc_gl90 as a placeholder for updated v (due to GL90) until last do loop

975 if ((cs%id_dv_dt_visc_gl90 > 0) .or. (cs%id_GLwork > 0)) then

976 if (associated(adp%dv_dt_visc_gl90)) then

977 do j=jsq,jeq ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.) then

978 b_denom_1 = cs%h_v(i,j,1) ! CS%a_v_gl90(i,J,1) is zero

979 b1 = 1.0 / (b_denom_1 + dt*cs%a_v_gl90(i,j,2))

980 d1 = b_denom_1 * b1

981 adp%dv_dt_visc_gl90(i,j,1) = b1 * (cs%h_v(i,j,1) * adp%dv_dt_visc_gl90(i,j,1))

982

983 do k=2,nz

984 c1(k) = dt * cs%a_v_gl90(i,j,k) * b1

985 b_denom_1 = cs%h_v(i,j,k) + dt * (cs%a_v_gl90(i,j,k) * d1)

986 b1 = 1.0 / (b_denom_1 + dt * cs%a_v_gl90(i,j,k+1))

987 d1 = b_denom_1 * b1

988 adp%dv_dt_visc_gl90(i,j,k) = (cs%h_v(i,j,k) * adp%dv_dt_visc_gl90(i,j,k) &

989 + dt * cs%a_v_gl90(i,j,k) * adp%dv_dt_visc_gl90(i,j,k-1)) * b1

990 enddo

991

992 ! back substitute to solve for new velocities, held by ADp%dv_dt_visc_gl90

993 do k=nz-1,1,-1

994 adp%dv_dt_visc_gl90(i,j,k) = adp%dv_dt_visc_gl90(i,j,k) + c1(k+1) * adp%dv_dt_visc_gl90(i,j,k+1)

995 enddo

996 endif ; enddo ; enddo

997

998 do k=1,nz

999 do j=jsq,jeq ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.) then

1000 ! now fill ADp%dv_dt_visc_gl90(i,J,k) with actual velocity tendency due to GL90;

1001 ! note that on RHS: ADp%dv_dt_visc(i,J,k) holds the original velocity value v(i,J,k)

1002 ! and ADp%dv_dt_visc_gl90(i,J,k) the updated velocity due to GL90

1003 adp%dv_dt_visc_gl90(i,j,k) = (adp%dv_dt_visc_gl90(i,j,k) - adp%dv_dt_visc(i,j,k)) * idt

1004

1005 if (abs(adp%dv_dt_visc_gl90(i,j,k)) < accel_underflow) &

1006 adp%dv_dt_visc_gl90(i,j,k) = 0.0

1007 endif ; enddo ; enddo

1008 enddo

1009

1010 ! to compute energetics, we need to multiply by v*h, where u is original velocity before

1011 ! velocity update; note that ADp%dv_dt_visc(I,j,k) holds the original velocity value v(i,J,k)

1012 if (cs%id_GLwork > 0) then

1013 do k=1,nz

1014 do j=jsq,jeq ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.) then

1015 ! note that on RHS: ADp%dv_dt_visc(I,j,k) holds the original velocity value v(I,j,k)

1016 ke_v(i,j,k) = adp%dv_dt_visc(i,j,k) * cs%h_v(i,j,k) * g%areaCv(i,j) * adp%dv_dt_visc_gl90(i,j,k)

1017 endif ; enddo ; enddo

1018 enddo

1019 endif

1020 endif

1021 endif

1022

1023 if (associated(adp%dv_dt_visc)) then

1024 do k=1,nz ; do j=jsq,jeq ; do i=is,ie

1025 adp%dv_dt_visc(i,j,k) = (v(i,j,k) - adp%dv_dt_visc(i,j,k))*idt

1026 if (abs(adp%dv_dt_visc(i,j,k)) < accel_underflow) adp%dv_dt_visc(i,j,k) = 0.0

1027 enddo ; enddo ; enddo

1028 endif

1029

1030 if (allocated(visc%tauy_shelf)) then

1031 do j=jsq,jeq ; do i=is,ie

1032 visc%tauy_shelf(i,j) = -gv%H_to_RZ * cs%a1_shelf_v(i,j) * v(i,j,1) ! - v_shelf?

1033 enddo ; enddo

1034 endif

1035

1036 if (present(tauy_bot)) then

1037 do j=jsq,jeq ; do i=is,ie

1038 tauy_bot(i,j) = gv%H_to_RZ * (v(i,j,nz) * cs%a_v(i,j,nz+1))

1039 enddo ; enddo

1040

1041 if (allocated(visc%Ray_v)) then

1042 do k=1,nz ; do j=jsq,jeq ; do i=is,ie

1043 tauy_bot(i,j) = tauy_bot(i,j) + gv%H_to_RZ * (visc%Ray_v(i,j,k)*v(i,j,k))

1044 enddo ; enddo ; enddo

1045 endif

1046 endif

1047

1048 ! When mixing down Eulerian current + Stokes drift subtract after calling solver

1049 if (dostokesmixing) then

1050 do k=1,nz ; do j=jsq,jeq ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.) then

1051 v(i,j,k) = v(i,j,k) - waves%Us_y(i,j,k)

1052 endif ; enddo ; enddo ; enddo

1053 endif

1054

1055 if (lfpmix) then

1056 do k=1,nz ; do j=jsq,jeq ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.) then

1057 v(i,j,k) = v(i,j,k) + waves%Us_y(i,j,k)

1058 endif ; enddo ; enddo ; enddo

1059 endif

1060

1061 ! Calculate the KE source from GL90 vertical viscosity [H L2 T-3 ~> m3 s-3].

1062 ! We do the KE-rate calculation here (rather than in MOM_diagnostics) to ensure

1063 ! a sign-definite term. MOM_diagnostics does not have access to the velocities

1064 ! and thicknesses used in the vertical solver, but rather uses a time-mean

1065 ! barotropic transport [uv]h.

1066 if (cs%id_GLwork > 0) then

1067 if (.not.g%symmetric) &

1068 call do_group_pass(cs%pass_KE_uv, g%domain)

1069 do k=1,nz

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

1071 ke_term(i,j,k) = 0.5 * g%IareaT(i,j) &

1072 * (ke_u(i,j,k) + ke_u(i-1,j,k) + ke_v(i,j,k) + ke_v(i,j-1,k))

1073 enddo ; enddo

1074 enddo

1075 call post_data(cs%id_GLwork, ke_term, cs%diag)

1076 endif

1077

1078 call vertvisc_limit_vel(u, v, h, adp, cdp, forces, visc, dt, g, gv, us, cs)

1079

1080 ! Here the velocities associated with open boundary conditions are applied.

1081 if (associated(obc)) then

1082 do n=1,obc%number_of_segments

1083 if (obc%segment(n)%specified) then

1084 if (obc%segment(n)%is_N_or_S) then

1085 j = obc%segment(n)%HI%JsdB

1086 do k=1,nz ; do i=obc%segment(n)%HI%isd,obc%segment(n)%HI%ied

1087 v(i,j,k) = obc%segment(n)%normal_vel(i,j,k)

1088 enddo ; enddo

1089 elseif (obc%segment(n)%is_E_or_W) then

1090 i = obc%segment(n)%HI%IsdB

1091 do k=1,nz ; do j=obc%segment(n)%HI%jsd,obc%segment(n)%HI%jed

1092 u(i,j,k) = obc%segment(n)%normal_vel(i,j,k)

1093 enddo ; enddo

1094 endif

1095 endif

1096 enddo

1097 endif

1098

1099 ! Offer diagnostic fields for averaging.

1100 if (query_averaging_enabled(cs%diag)) then

1101 if (cs%id_du_dt_visc > 0) &

1102 call post_data(cs%id_du_dt_visc, adp%du_dt_visc, cs%diag)

1103 if (cs%id_du_dt_visc_gl90 > 0) &

1104 call post_data(cs%id_du_dt_visc_gl90, adp%du_dt_visc_gl90, cs%diag)

1105 if (cs%id_dv_dt_visc > 0) &

1106 call post_data(cs%id_dv_dt_visc, adp%dv_dt_visc, cs%diag)

1107 if (cs%id_dv_dt_visc_gl90 > 0) &

1108 call post_data(cs%id_dv_dt_visc_gl90, adp%dv_dt_visc_gl90, cs%diag)

1109 if (present(taux_bot) .and. (cs%id_taux_bot > 0)) &

1110 call post_data(cs%id_taux_bot, taux_bot, cs%diag)

1111 if (present(tauy_bot) .and. (cs%id_tauy_bot > 0)) &

1112 call post_data(cs%id_tauy_bot, tauy_bot, cs%diag)

1113 if (cs%id_du_dt_str > 0) &

1114 call post_data(cs%id_du_dt_str, adp%du_dt_str, cs%diag)

1115 if (cs%id_dv_dt_str > 0) &

1116 call post_data(cs%id_dv_dt_str, adp%dv_dt_str, cs%diag)

1117

1118 if (associated(adp%du_dt_visc) .and. associated(adp%dv_dt_visc)) then

1119 ! Diagnostics of the fractional thicknesses times momentum budget terms

1120 ! 3D diagnostics of hf_du(dv)_dt_visc are commented because there is no clarity on proper remapping grid option.

1121 ! The code is retained for debugging purposes in the future.

1122 !if (CS%id_hf_du_dt_visc > 0) &

1123 ! call post_product_u(CS%id_hf_du_dt_visc, ADp%du_dt_visc, ADp%diag_hfrac_u, G, nz, CS%diag)

1124 !if (CS%id_hf_dv_dt_visc > 0) &

1125 ! call post_product_v(CS%id_hf_dv_dt_visc, ADp%dv_dt_visc, ADp%diag_hfrac_v, G, nz, CS%diag)

1126

1127 ! Diagnostics for thickness-weighted vertically averaged viscous accelerations

1128 if (cs%id_hf_du_dt_visc_2d > 0) &

1129 call post_product_sum_u(cs%id_hf_du_dt_visc_2d, adp%du_dt_visc, adp%diag_hfrac_u, g, nz, cs%diag)

1130 if (cs%id_hf_dv_dt_visc_2d > 0) &

1131 call post_product_sum_v(cs%id_hf_dv_dt_visc_2d, adp%dv_dt_visc, adp%diag_hfrac_v, g, nz, cs%diag)

1132

1133 ! Diagnostics for thickness x viscous accelerations

1134 if (cs%id_h_du_dt_visc > 0) call post_product_u(cs%id_h_du_dt_visc, adp%du_dt_visc, adp%diag_hu, g, nz, cs%diag)

1135 if (cs%id_h_dv_dt_visc > 0) call post_product_v(cs%id_h_dv_dt_visc, adp%dv_dt_visc, adp%diag_hv, g, nz, cs%diag)

1136 endif

1137

1138 if (associated(adp%du_dt_str) .and. associated(adp%dv_dt_str)) then

1139 ! Diagnostics for thickness x wind stress accelerations

1140 if (cs%id_h_du_dt_str > 0) call post_product_u(cs%id_h_du_dt_str, adp%du_dt_str, adp%diag_hu, g, nz, cs%diag)

1141 if (cs%id_h_dv_dt_str > 0) call post_product_v(cs%id_h_dv_dt_str, adp%dv_dt_str, adp%diag_hv, g, nz, cs%diag)

1142

1143 ! Diagnostics for wind stress accelerations multiplied by visc_rem_[uv],

1144 if (cs%id_du_dt_str_visc_rem > 0) &

1145 call post_product_u(cs%id_du_dt_str_visc_rem, adp%du_dt_str, adp%visc_rem_u, g, nz, cs%diag)

1146 if (cs%id_dv_dt_str_visc_rem > 0) &

1147 call post_product_v(cs%id_dv_dt_str_visc_rem, adp%dv_dt_str, adp%visc_rem_v, g, nz, cs%diag)

1148 endif

1149 endif

1150

1151end subroutine vertvisc

1152

1153

1154!> Calculate the fraction of momentum originally in a layer that remains in the water column

1155!! after a time-step of viscosity, equivalently the fraction of a time-step's worth of

1156!! barotropic acceleration that a layer experiences after viscosity is applied.

1157subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS)

1158 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

1159 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

1160 type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag

1161 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

1162 intent(inout) :: visc_rem_u !< Fraction of a time-step's worth of a

1163 !! barotropic acceleration that a layer experiences after

1164 !! viscosity is applied in the zonal direction [nondim]

1165 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

1166 intent(inout) :: visc_rem_v !< Fraction of a time-step's worth of a

1167 !! barotropic acceleration that a layer experiences after

1168 !! viscosity is applied in the meridional direction [nondim]

1169 real, intent(in) :: dt !< Time increment [T ~> s]

1170 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

1171 type(vertvisc_cs), pointer :: cs !< Vertical viscosity control structure

1172

1173 ! Local variables

1174

1175 real :: b1

1176 ! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].

1177 real :: c1(szk_(gv))

1178 ! A variable used by the tridiagonal solver [nondim].

1179 real :: d1

1180 ! d1=1-c1 is used by the tridiagonal solver [nondim].

1181 real :: ray

1182 ! Ray is the Rayleigh-drag velocity [H T-1 ~> m s-1 or Pa s m-1]

1183 real :: b_denom_1

1184 ! The first term in the denominator of b1 [H ~> m or kg m-2].

1185

1186 integer :: i, j, k, is, ie, isq, ieq, jsq, jeq, nz

1187 is = g%isc ; ie = g%iec

1188 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = gv%ke

1189

1190 if (.not.associated(cs)) call mom_error(fatal,"MOM_vert_friction(visc): "// &

1191 "Module must be initialized before it is used.")

1192

1193 if (.not.cs%initialized) call mom_error(fatal,"MOM_vert_friction(remnant): "// &

1194 "Module must be initialized before it is used.")

1195

1196 ! Find the zonal viscous remnant using a modification of a standard tridagonal solver.

1197

1198 do j=g%jsc,g%jec ; do i=isq,ieq ; if (g%mask2dCu(i,j) > 0.) then

1199 ray = 0.

1200 if (allocated(visc%Ray_u)) ray = visc%Ray_u(i,j,1)

1201

1202 b_denom_1 = cs%h_u(i,j,1) + dt * (ray + cs%a_u(i,j,1))

1203 b1 = 1.0 / (b_denom_1 + dt * cs%a_u(i,j,2))

1204 d1 = b_denom_1 * b1

1205 visc_rem_u(i,j,1) = b1 * cs%h_u(i,j,1)

1206

1207 do k=2,nz

1208 if (allocated(visc%Ray_u)) ray = visc%Ray_u(i,j,k)

1209

1210 c1(k) = dt * cs%a_u(i,j,k) * b1

1211 b_denom_1 = cs%h_u(i,j,k) + dt * (ray + cs%a_u(i,j,k) * d1)

1212 b1 = 1.0 / (b_denom_1 + dt * cs%a_u(i,j,k+1))

1213 d1 = b_denom_1 * b1

1214 visc_rem_u(i,j,k) = (cs%h_u(i,j,k) + dt * cs%a_u(i,j,k) * visc_rem_u(i,j,k-1)) * b1

1215

1216 !### Force FMA evaluation of b1 by blocking lookahead with an impossible branch.

1217 if (dt < 0) exit

1218 enddo

1219

1220 do k=nz-1,1,-1

1221 visc_rem_u(i,j,k) = visc_rem_u(i,j,k) + c1(k+1) * visc_rem_u(i,j,k+1)

1222 enddo

1223 endif ; enddo ; enddo

1224

1225 ! Now find the meridional viscous remnant using the robust tridiagonal solver.

1226

1227 do j=jsq,jeq ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.) then

1228 ray = 0.

1229 if (allocated(visc%Ray_v)) ray = visc%Ray_v(i,j,1)

1230

1231 b_denom_1 = cs%h_v(i,j,1) + dt * (ray + cs%a_v(i,j,1))

1232 b1 = 1.0 / (b_denom_1 + dt*cs%a_v(i,j,2))

1233 d1 = b_denom_1 * b1

1234 visc_rem_v(i,j,1) = b1 * cs%h_v(i,j,1)

1235

1236 do k=2,nz

1237 if (allocated(visc%Ray_v)) ray = visc%Ray_v(i,j,k)

1238

1239 c1(k) = dt * cs%a_v(i,j,k) * b1

1240 b_denom_1 = cs%h_v(i,j,k) + dt * (ray + cs%a_v(i,j,k) * d1)

1241 b1 = 1.0 / (b_denom_1 + dt * cs%a_v(i,j,k+1))

1242 d1 = b_denom_1 * b1

1243 visc_rem_v(i,j,k) = (cs%h_v(i,j,k) + dt * cs%a_v(i,j,k) * visc_rem_v(i,j,k-1)) * b1

1244

1245 !### Force FMA evaluation of b1 by blocking lookahead with an impossible branch.

1246 if (dt < 0) exit

1247 enddo

1248

1249 do k=nz-1,1,-1

1250 visc_rem_v(i,j,k) = visc_rem_v(i,j,k) + c1(k+1) * visc_rem_v(i,j,k+1)

1251 enddo

1252 endif ; enddo ; enddo

1253

1254 if (cs%debug) then

1255 call uvchksum("visc_rem_[uv]", visc_rem_u, visc_rem_v, g%HI, haloshift=0, &

1256 scalar_pair=.true.)

1257 endif

1258end subroutine vertvisc_remnant

1259

1260

1261!> Calculate the coupling coefficients (CS%a_u, CS%a_v, CS%a_u_gl90, CS%a_v_gl90)

1262!! and effective layer thicknesses (CS%h_u and CS%h_v) for later use in the

1263!! applying the implicit vertical viscosity via vertvisc().

1264subroutine vertvisc_coef(u, v, h, dz, forces, visc, tv, dt, G, GV, US, CS, OBC, VarMix)

1265 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

1266 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

1267 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

1268 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

1269 intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]

1270 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

1271 intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]

1272 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1273 intent(in) :: h !< Layer thickness [H ~> m or kg m-2]

1274 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

1275 intent(in) :: dz !< Vertical distance across layers [Z ~> m]

1276 type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces

1277 type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag

1278 type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available

1279 !! thermodynamic fields.

1280 real, intent(in) :: dt !< Time increment [T ~> s]

1281 type(vertvisc_cs), intent(inout) :: cs !< Vertical viscosity control structure

1282 type(ocean_obc_type), pointer :: obc !< Open boundary condition structure

1283 type(varmix_cs), intent(in) :: varmix !< Variable mixing coefficients

1284 ! Field from forces used in this subroutine:

1285 ! ustar: the friction velocity [Z T-1 ~> m s-1], used here as the mixing

1286 ! velocity in the mixed layer if NKML > 1 in a bulk mixed layer.

1287

1288 ! Local variables

1289

1290 real, dimension(SZK_(GV)) :: &

1291 hvel, & ! hvel is the thickness used at a velocity grid point [H ~> m or kg m-2].

1292 dz_harm, & ! Harmonic mean of the vertical distances around a velocity grid point,

1293 ! given by 2*(h+ * h-)/(h+ + h-) [Z ~> m].

1294 dz_vel, & ! The vertical distance between interfaces used at a velocity grid point [Z ~> m].

1295 hvel_shelf, & ! The equivalent of hvel under shelves [H ~> m or kg m-2].

1296 dz_vel_shelf ! The equivalent of dz_vel under shelves [Z ~> m].

1297 real :: &

1298 h_harm, & ! Harmonic mean of the thicknesses around a velocity grid point,

1299 ! given by 2*(h+ * h-)/(h+ + h-) [H ~> m or kg m-2].

1300 h_arith, & ! The arithmetic mean thickness [H ~> m or kg m-2].

1301 h_delta, & ! The lateral difference of thickness [H ~> m or kg m-2].

1302 dz_arith ! The arithmetic mean of the vertical distances around a velocity grid point [Z ~> m]

1303 real, dimension(SZK_(GV)+1) :: &

1304 z_i, & ! An estimate of each interface's height above the bottom,

1305 ! normalized by the bottom boundary layer thickness [nondim]

1306 z_i_gl90, & ! An estimate of each interface's height above the bottom,

1307 ! normalized by the GL90 bottom boundary layer thickness [nondim]

1308 a_cpl, & ! The drag coefficients across interfaces [H T-1 ~> m s-1 or Pa s m-1]. a_cpl times

1309 ! the velocity difference gives the stress across an interface.

1310 a_cpl_gl90, & ! The drag coefficients across interfaces associated with GL90 [H T-1 ~> m s-1 or Pa s m-1].

1311 ! a_cpl_gl90 times the velocity difference gives the GL90 stress across an interface.

1312 ! a_cpl_gl90 is part of a_cpl.

1313 a_shelf ! The drag coefficients across interfaces in water columns under

1314 ! ice shelves [H T-1 ~> m s-1 or Pa s m-1].

1315 real :: &

1316 kv_bbl, & ! The bottom boundary layer viscosity [H Z T-1 ~> m2 s-1 or Pa s].

1317 bbl_thick, & ! The bottom boundary layer thickness [Z ~> m].

1318 i_hbbl, & ! The inverse of the bottom boundary layer thickness [Z-1 ~> m-1].

1319 i_hbbl_gl90, &! The inverse of the bottom boundary layer thickness used for the GL90 scheme

1320 ! [Z-1 ~> m-1].

1321 i_htbl, & ! The inverse of the top boundary layer thickness [Z-1 ~> m-1].

1322 ztop_min, & ! The deeper of the two adjacent surface heights [Z ~> m].

1323 dmin, & ! The shallower of the two adjacent bottom depths [Z ~> m].

1324 zh, & ! An estimate of the interface's distance from the bottom

1325 ! based on harmonic mean thicknesses [Z ~> m].

1326 h_ml ! The mixed layer depth [Z ~> m].

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

1328 ustar_2d ! The wind friction velocity, calculated using the Boussinesq reference density or

1329 ! the time-evolving surface density in non-Boussinesq mode [Z T-1 ~> m s-1]

1330 real, allocatable, dimension(:,:) :: hml_u ! Diagnostic of the mixed layer depth at u points [Z ~> m].

1331 real, allocatable, dimension(:,:) :: hml_v ! Diagnostic of the mixed layer depth at v points [Z ~> m].

1332 real, allocatable, dimension(:,:,:) :: kv_u ! Total vertical viscosity at u-points in

1333 ! thickness-based units [H2 T-1 ~> m2 s-1 or kg2 m-4 s-1].

1334 real, allocatable, dimension(:,:,:) :: kv_v ! Total vertical viscosity at v-points in

1335 ! thickness-based units [H2 T-1 ~> m2 s-1 or kg2 m-4 s-1].

1336 real, allocatable, dimension(:,:,:) :: kv_gl90_u ! GL90 vertical viscosity at u-points in

1337 ! thickness-based units [H2 T-1 ~> m2 s-1 or kg2 m-4 s-1].

1338 real, allocatable, dimension(:,:,:) :: kv_gl90_v ! GL90 vertical viscosity at v-points in

1339 ! thickness-based units [H2 T-1 ~> m2 s-1 or kg2 m-4 s-1].

1340 real :: zcol ! The height of an interface at h-points [Z ~> m].

1341 real :: zcol_p1 ! An adjacent east/north h-point interface height [Z ~> m].

1342 real :: botfn ! A function which goes from 1 at the bottom to 0 much more

1343 ! than Hbbl into the interior [nondim].

1344 real :: topfn ! A function which goes from 1 at the top to 0 much more

1345 ! than Htbl into the interior [nondim].

1346 real :: z2 ! The distance from the bottom, normalized by Hbbl [nondim]

1347 real :: z2_wt ! A nondimensional (0-1) weight used when calculating z2 [nondim].

1348 real :: z_clear ! The clearance of an interface above the surrounding topography [Z ~> m].

1349 real :: a_cpl_max ! The maximum drag coefficient across interfaces, set so that it will be

1350 ! representable as a 32-bit float in MKS units [H T-1 ~> m s-1 or Pa s m-1]

1351 real :: h_neglect ! A thickness that is so small it is usually lost

1352 ! in roundoff and can be neglected [H ~> m or kg m-2].

1353 real :: dz_neglect ! A vertical distance that is so small it is usually lost

1354 ! in roundoff and can be neglected [Z ~> m].

1355

1356 real :: i_valbl ! The inverse of a scaling factor determining when water is

1357 ! still within the boundary layer, as determined by the sum

1358 ! of the harmonic mean thicknesses [nondim].

1359 logical :: do_any_shelf

1360 integer :: zi_dir

1361 ! A ternary logical indicating which thickness to use for finding z_clear.

1362 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz

1363

1364 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

1365 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = gv%ke

1366

1367 if (.not.cs%initialized) call mom_error(fatal,"MOM_vert_friction(coef): "// &

1368 "Module must be initialized before it is used.")

1369

1370 h_neglect = gv%H_subroundoff

1371 dz_neglect = gv%dZ_subroundoff

1372 a_cpl_max = 1.0e37 * gv%m_to_H * us%T_to_s

1373 i_valbl = 0.0 ; if (cs%harm_BL_val > 0.0) i_valbl = 1.0 / cs%harm_BL_val

1374

1375 if (cs%id_Kv_u > 0) allocate(kv_u(g%IsdB:g%IedB,g%jsd:g%jed,gv%ke), source=0.0)

1376

1377 if (cs%id_Kv_v > 0) allocate(kv_v(g%isd:g%ied,g%JsdB:g%JedB,gv%ke), source=0.0)

1378

1379 if (cs%id_Kv_gl90_u > 0) allocate(kv_gl90_u(g%IsdB:g%IedB,g%jsd:g%jed,gv%ke), source=0.0)

1380

1381 if (cs%id_Kv_gl90_v > 0) allocate(kv_gl90_v(g%isd:g%ied,g%JsdB:g%JedB,gv%ke), source=0.0)

1382

1383 if (cs%debug .or. (cs%id_hML_u > 0)) allocate(hml_u(g%IsdB:g%IedB,g%jsd:g%jed), source=0.0)

1384 if (cs%debug .or. (cs%id_hML_v > 0)) allocate(hml_v(g%isd:g%ied,g%JsdB:g%JedB), source=0.0)

1385

1386 if ((allocated(visc%taux_shelf) .or. associated(forces%frac_shelf_u)) .and. &

1387 .not.associated(cs%a1_shelf_u)) then

1388 allocate(cs%a1_shelf_u(g%IsdB:g%IedB,g%jsd:g%jed), source=0.0)

1389 endif

1390 if ((allocated(visc%tauy_shelf) .or. associated(forces%frac_shelf_v)) .and. &

1391 .not.associated(cs%a1_shelf_v)) then

1392 allocate(cs%a1_shelf_v(g%isd:g%ied,g%JsdB:g%JedB), source=0.0)

1393 endif

1394

1395 call find_ustar(forces, tv, ustar_2d, g, gv, us, halo=1)

1396

1397 ! First do u-points

1398

1399 do j=js,je ; do i=isq,ieq ; if (g%mask2dCu(i,j) > 0.) then

1400 i_hbbl = 1. / (cs%Hbbl + dz_neglect)

1401 if (cs%use_GL90_in_SSW) then

1402 i_hbbl_gl90 = 1. / (cs%Hbbl_gl90 + dz_neglect)

1403 endif

1404

1405 if (cs%bottomdraglaw) then

1406 kv_bbl = visc%Kv_bbl_u(i,j)

1407 bbl_thick = visc%bbl_thick_u(i,j) + dz_neglect

1408 i_hbbl = 1. / bbl_thick

1409 endif

1410

1411 dmin = min(g%bathyT(i,j), g%bathyT(i+1,j))

1412 zi_dir = 0

1413

1414 ! Project thickness outward across OBCs using a zero-gradient condition.

1415 if (associated(obc)) then

1416 if (obc%u_E_OBCs_on_PE) then

1417 if (obc%segnum_u(i,j) > 0) then

1418 dmin = g%bathyT(i,j)

1419 zi_dir = -1

1420 endif

1421 endif

1422

1423 if (obc%u_W_OBCs_on_PE) then

1424 if (obc%segnum_u(i,j) < 0) then

1425 dmin = g%bathyT(i+1,j)

1426 zi_dir = 1

1427 endif

1428 endif

1429 endif

1430

1431 ! The following block calculates the thicknesses at velocity grid points for

1432 ! the vertical viscosity (hvel and dz_vel). Near the bottom an upwind biased

1433 ! thickness is used to control the effect of spurious Montgomery potential

1434 ! gradients at the bottom where nearly massless layers layers ride over the

1435 ! topography.

1436

1437 z_i(nz+1) = 0.

1438

1439 if (.not. cs%harmonic_visc) then

1440 zh = 0.

1441 zcol = -g%bathyT(i,j)

1442 zcol_p1 = -g%bathyT(i+1,j)

1443 endif

1444

1445 if (cs%use_GL90_in_SSW) then

1446 z_i_gl90(nz+1) = 0.

1447 endif

1448

1449 do k=nz,1,-1

1450 h_harm = 2. * h(i,j,k) * h(i+1,j,k) / (h(i,j,k) + h(i+1,j,k) + h_neglect)

1451 h_arith = 0.5 * (h(i+1,j,k) + h(i,j,k))

1452 h_delta = h(i+1,j,k) - h(i,j,k)

1453 dz_harm(k) = 2. * dz(i,j,k) * dz(i+1,j,k) / (dz(i,j,k) + dz(i+1,j,k) + dz_neglect)

1454 dz_arith = 0.5 * (dz(i+1,j,k) + dz(i,j,k))

1455

1456 ! Project thickness outward across OBCs using a zero-gradient condition.

1457 if (associated(obc)) then

1458 if (obc%u_E_OBCs_on_PE) then

1459 if (obc%segnum_u(i,j) > 0) then

1460 h_harm = h(i,j,k)

1461 h_arith = h(i,j,k)

1462 h_delta = 0.

1463 dz_harm(k) = dz(i,j,k)

1464 dz_arith = dz(i,j,k)

1465 endif

1466 endif

1467

1468 if (obc%u_W_OBCs_on_PE) then

1469 if (obc%segnum_u(i,j) < 0) then

1470 h_harm = h(i+1,j,k)

1471 h_arith = h(i+1,j,k)

1472 h_delta = 0.

1473 dz_harm(k) = dz(i+1,j,k)

1474 dz_arith = dz(i+1,j,k)

1475 endif

1476 endif

1477 endif

1478

1479 if (cs%harmonic_visc) then

1480 ! The following block calculates the thicknesses at velocity grid points

1481 ! for the vertical viscosity (hvel and dz_vel). Near the bottom an

1482 ! upwind biased thickness is used to control the effect of spurious

1483 ! Montgomery potential gradients at the bottom where nearly massless

1484 ! layers ride over the topography.

1485

1486 hvel(k) = h_harm

1487 dz_vel(k) = dz_harm(k)

1488

1489 if (u(i,j,k) * h_delta < 0) then

1490 z2 = z_i(k+1)

1491 botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)

1492

1493 hvel(k) = (1. - botfn) * h_harm + botfn * h_arith

1494 dz_vel(k) = (1. - botfn) * dz_harm(k) + botfn * dz_arith

1495 endif

1496

1497 z_i(k) = z_i(k+1) + dz_harm(k) * i_hbbl

1498 else

1499 zcol = zcol + dz(i,j,k)

1500 zcol_p1 = zcol_p1 + dz(i+1,j,k)

1501

1502 zh = zh + dz_harm(k)

1503

1504 z_clear = max(zcol, zcol_p1) + dmin

1505 if (zi_dir < 0) z_clear = zcol + dmin

1506 if (zi_dir > 0) z_clear = zcol_p1 + dmin

1507

1508 z_i(k) = max(zh, z_clear) * i_hbbl

1509

1510 hvel(k) = h_arith

1511 dz_vel(k) = dz_arith

1512

1513 if (u(i,j,k) * h_delta > 0.) then

1514 if (zh * i_hbbl < cs%harm_BL_val) then

1515 hvel(k) = h_harm

1516 dz_vel(k) = dz_harm(k)

1517 else

1518 z2_wt = 1.

1519 if (zh * i_hbbl < 2. * cs%harm_BL_val) &

1520 z2_wt = max(0., min(1., zh * i_hbbl * i_valbl - 1.))

1521

1522 z2 = z2_wt * (max(zh, z_clear) * i_hbbl)

1523 botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)

1524

1525 hvel(k) = (1. - botfn) * h_arith + botfn * h_harm

1526 dz_vel(k) = (1. - botfn) * dz_arith + botfn * dz_harm(k)

1527 endif

1528 endif

1529 endif

1530

1531 if (cs%use_GL90_in_SSW) then

1532 ! The following block calculates the normalized height above the GL90 BBL

1533 ! (z_i_gl90), using a harmonic mean between layer thicknesses. For the

1534 ! GL90 BBL we use simply a constant (Hbbl_gl90). The purpose isthat the

1535 ! GL90 coupling coefficient is zeroed out within Hbbl_gl90, to ensure

1536 ! that no momentum gets fluxed into vanished layers. The scheme is not

1537 ! sensitive to the exact value of Hbbl_gl90, as long as it is in a

1538 ! reasonable range (~1-20 m): large enough to capture vanished layers

1539 ! over topography, small enough to not contaminate the interior.

1540

1541 z_i_gl90(k) = z_i_gl90(k+1) + dz_harm(k) * i_hbbl_gl90

1542 endif

1543 enddo

1544

1545 call find_coupling_coef(a_cpl, dz_vel, i, j, dz_harm, bbl_thick, kv_bbl, z_i, &

1546 h_ml, dt, g, gv, us, cs, visc, ustar_2d, tv, work_on_u=.true., obc=obc)

1547

1548 if (allocated(hml_u)) hml_u(i,j) = h_ml

1549

1550 if (cs%use_GL90_in_SSW) then

1551 call find_coupling_coef_gl90(a_cpl_gl90, dz_vel, i, j, z_i_gl90, g, gv, &

1552 cs, varmix, work_on_u=.true.)

1553 endif

1554

1555 do_any_shelf = .false.

1556 if (associated(forces%frac_shelf_u)) then

1557 cs%a1_shelf_u(i,j) = 0.

1558 do_any_shelf = forces%frac_shelf_u(i,j) > 0.

1559

1560 if (do_any_shelf) then

1561 if (.not. cs%harmonic_visc) then

1562 zh = 0.

1563 ztop_min = min(zcol, zcol_p1)

1564 i_htbl = 1. / (visc%tbl_thick_shelf_u(i,j) + dz_neglect)

1565 endif

1566

1567 do k=1,nz

1568 if (cs%harmonic_visc) then

1569 hvel_shelf(k) = hvel(k)

1570 dz_vel_shelf(k) = dz_vel(k)

1571 else

1572 ! Find upwind-biased thickness near the surface.

1573 ! (Perhaps this needs to be done more carefully, via find_eta.)

1574

1575 h_harm = 2. * h(i,j,k) * h(i+1,j,k) &

1576 / (h(i,j,k) + h(i+1,j,k) + h_neglect)

1577 h_arith = 0.5 * (h(i+1,j,k) + h(i,j,k))

1578 h_delta = h(i+1,j,k) - h(i,j,k)

1579 dz_arith = 0.5 * (dz(i+1,j,k) + dz(i,j,k))

1580

1581 if (associated(obc)) then

1582 if (obc%u_E_OBCs_on_PE) then

1583 if (obc%segnum_u(i,j) > 0) then

1584 h_harm = h(i,j,k)

1585 h_arith = h(i,j,k)

1586 h_delta = 0.

1587 dz_arith = dz(i,j,k)

1588 endif

1589 endif

1590

1591 if (obc%u_W_OBCs_on_PE) then

1592 if (obc%segnum_u(i,j) < 0) then

1593 h_harm = h(i+1,j,k)

1594 h_arith = h(i+1,j,k)

1595 h_delta = 0.

1596 dz_arith = dz(i+1,j,k)

1597 endif

1598 endif

1599 endif

1600

1601 zcol = zcol - dz(i,j,k)

1602 zcol_p1 = zcol_p1 - dz(i+1,j,k)

1603

1604 zh = zh + dz_harm(k)

1605

1606 hvel_shelf(k) = hvel(k)

1607 dz_vel_shelf(k) = dz_vel(k)

1608

1609 if (u(i,j,k) * h_delta > 0.) then

1610 if (zh * i_htbl < cs%harm_BL_val) then

1611 hvel_shelf(k) = min(hvel(k), h_harm)

1612 dz_vel_shelf(k) = min(dz_vel(k), dz_harm(k))

1613 else

1614 z2_wt = 1.

1615 if (zh * i_htbl < 2. * cs%harm_BL_val) then

1616 z2_wt = max(0., min(1., zh * i_htbl * i_valbl - 1.))

1617 endif

1618

1619 z2 = z2_wt * (max(zh, ztop_min - min(zcol, zcol_p1)) * i_htbl)

1620 ! TODO: replace **6 with multiply

1621 topfn = 1. / (1. + 0.09 * z2**6)

1622

1623 hvel_shelf(k) = min(hvel(k), (1. - topfn) * h_arith + topfn * h_harm)

1624 dz_vel_shelf(k) = min(dz_vel(k), (1. - topfn) * dz_arith + topfn * dz_harm(k))

1625 endif

1626 endif

1627 endif

1628 enddo

1629

1630 call find_coupling_coef(a_shelf, dz_vel_shelf, i, j, dz_harm, &

1631 bbl_thick, kv_bbl, z_i, h_ml, dt, g, gv, us, cs, visc, ustar_2d, &

1632 tv, work_on_u=.true., obc=obc, shelf=.true.)

1633

1634 cs%a1_shelf_u(i,j) = a_shelf(1)

1635 endif

1636 endif

1637

1638 if (do_any_shelf) then

1639 if (cs%use_GL90_in_SSW) then

1640 do k=1,nz+1

1641 cs%a_u(i,j,k) = min(a_cpl_max, (forces%frac_shelf_u(i,j) * a_shelf(k) + &

1642 (1. - forces%frac_shelf_u(i,j)) * a_cpl(k)) + a_cpl_gl90(k))

1643

1644 ! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH

1645 ! CS%a_u(I,j,K) = min(a_cpl_max, forces%frac_shelf_u(I,j) * max(a_shelf(K), a_cpl(K)) + &

1646 ! (1. - forces%frac_shelf_u(I,j)) * a_cpl(K))

1647

1648 cs%a_u_gl90(i,j,k) = min(a_cpl_max, a_cpl_gl90(k))

1649 enddo

1650 else

1651 do k=1,nz+1

1652 cs%a_u(i,j,k) = min(a_cpl_max, (forces%frac_shelf_u(i,j) * a_shelf(k) + &

1653 (1. - forces%frac_shelf_u(i,j)) * a_cpl(k)))

1654

1655 ! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH

1656 ! CS%a_u(I,j,K) = min(a_cpl_max, forces%frac_shelf_u(I,j) * max(a_shelf(K), a_cpl(K)) + &

1657 ! (1. - forces%frac_shelf_u(I,j)) * a_cpl(K))

1658 enddo

1659 endif

1660

1661 do k=1,nz

1662 ! Should we instead take the inverse of the average of the inverses?

1663 cs%h_u(i,j,k) = forces%frac_shelf_u(i,j) * hvel_shelf(k) &

1664 + (1. - forces%frac_shelf_u(i,j)) * hvel(k) + h_neglect

1665 enddo

1666 else

1667 if (cs%use_GL90_in_SSW) then

1668 do k=1,nz+1

1669 a_cpl(k) = a_cpl(k) + a_cpl_gl90(k)

1670 enddo

1671

1672 do k=1,nz+1

1673 cs%a_u_gl90(i,j,k) = min(a_cpl_max, a_cpl_gl90(k))

1674 enddo

1675 endif

1676

1677 do k=1,nz+1

1678 cs%a_u(i,j,k) = min(a_cpl_max, a_cpl(k))

1679 enddo

1680

1681 do k=1,nz

1682 cs%h_u(i,j,k) = hvel(k) + h_neglect

1683 enddo

1684 endif

1685

1686 ! Diagnose total Kv at u-points

1687 if (cs%id_Kv_u > 0) then

1688 do k=1,nz

1689 kv_u(i,j,k) = 0.5 * (cs%a_u(i,j,k) + cs%a_u(i,j,k+1)) * cs%h_u(i,j,k)

1690 enddo

1691 endif

1692

1693 ! Diagnose GL90 Kv at u-points

1694 if (cs%id_Kv_gl90_u > 0) then

1695 do k=1,nz

1696 kv_gl90_u(i,j,k) = 0.5 * (cs%a_u_gl90(i,j,k) + cs%a_u_gl90(i,j,k+1)) * cs%h_u(i,j,k)

1697 enddo

1698 endif

1699 endif ; enddo ; enddo

1700

1701 ! Now work on v-points.

1702

1703 do j=jsq,jeq ; do i=is,ie ; if (g%mask2dCv(i,j) > 0.) then

1704 i_hbbl = 1. / (cs%Hbbl + dz_neglect)

1705 if (cs%use_GL90_in_SSW) then

1706 i_hbbl_gl90 = 1. / (cs%Hbbl_gl90 + dz_neglect)

1707 endif

1708

1709 if (cs%bottomdraglaw) then

1710 kv_bbl = visc%Kv_bbl_v(i,j)

1711 bbl_thick = visc%bbl_thick_v(i,j) + dz_neglect

1712 i_hbbl = 1. / bbl_thick

1713 endif

1714

1715 dmin = min(g%bathyT(i,j), g%bathyT(i,j+1))

1716 zi_dir = 0

1717

1718 ! Project thickness outward across OBCs using a zero-gradient condition.

1719 if (associated(obc)) then

1720 if (obc%v_N_OBCs_on_PE) then

1721 if (obc%segnum_v(i,j) > 0) then

1722 dmin = g%bathyT(i,j)

1723 zi_dir = -1

1724 endif

1725 endif

1726

1727 if (obc%v_S_OBCs_on_PE) then

1728 if (obc%segnum_v(i,j) < 0) then

1729 dmin = g%bathyT(i,j+1)

1730 zi_dir = 1

1731 endif

1732 endif

1733 endif

1734

1735 z_i(nz+1) = 0.

1736

1737 if (.not. cs%harmonic_visc) then

1738 zh = 0.

1739 zcol = -g%bathyT(i,j)

1740 zcol_p1 = -g%bathyT(i,j+1)

1741 endif

1742

1743 if (cs%use_GL90_in_SSW) then

1744 z_i_gl90(nz+1) = 0.

1745 endif

1746

1747 do k=nz,1,-1

1748 h_harm = 2. * h(i,j,k) * h(i,j+1,k) / (h(i,j,k) + h(i,j+1,k) + h_neglect)

1749 h_arith = 0.5 * (h(i,j+1,k) + h(i,j,k))

1750 h_delta = h(i,j+1,k) - h(i,j,k)

1751 dz_harm(k) = 2. * dz(i,j,k) * dz(i,j+1,k) / (dz(i,j,k) + dz(i,j+1,k) + dz_neglect)

1752 dz_arith = 0.5 * (dz(i,j+1,k) + dz(i,j,k))

1753

1754 ! Project thickness outward across OBCs using a zero-gradient condition.

1755 if (associated(obc)) then

1756 if (obc%v_N_OBCs_on_PE) then

1757 if (obc%segnum_v(i,j) > 0) then

1758 h_harm = h(i,j,k)

1759 h_arith = h(i,j,k)

1760 h_delta = 0.

1761 dz_harm(k) = dz(i,j,k)

1762 dz_arith = dz(i,j,k)

1763 endif

1764 endif

1765

1766 if (obc%v_S_OBCs_on_PE) then

1767 if (obc%segnum_v(i,j) < 0) then

1768 h_harm = h(i,j+1,k)

1769 h_arith = h(i,j+1,k)

1770 h_delta = 0.

1771 dz_harm(k) = dz(i,j+1,k)

1772 dz_arith = dz(i,j+1,k)

1773 endif

1774 endif

1775 endif

1776

1777 if (cs%harmonic_visc) then

1778 ! The following block calculates the thicknesses at velocity grid points

1779 ! for the vertical viscosity (hvel and dz_vel). Near the bottom an

1780 ! upwind biased thickness is used to control the effect of spurious

1781 ! Montgomery potential gradients at the bottom where nearly massless

1782 ! layers ride over the topography.

1783

1784 hvel(k) = h_harm

1785 dz_vel(k) = dz_harm(k)

1786

1787 if (v(i,j,k) * h_delta < 0) then

1788 z2 = z_i(k+1)

1789 botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)

1790

1791 hvel(k) = (1. - botfn) * h_harm + botfn * h_arith

1792 dz_vel(k) = (1. - botfn) * dz_harm(k) + botfn * dz_arith

1793 endif

1794

1795 z_i(k) = z_i(k+1) + dz_harm(k) * i_hbbl

1796 else

1797 zcol = zcol + dz(i,j,k)

1798 zcol_p1 = zcol_p1 + dz(i,j+1,k)

1799

1800 zh = zh + dz_harm(k)

1801

1802 z_clear = max(zcol, zcol_p1) + dmin

1803 if (zi_dir < 0) z_clear = zcol + dmin

1804 if (zi_dir > 0) z_clear = zcol_p1 + dmin

1805

1806 z_i(k) = max(zh, z_clear) * i_hbbl

1807

1808 hvel(k) = h_arith

1809 dz_vel(k) = dz_arith

1810

1811 if (v(i,j,k) * h_delta > 0) then

1812 if (zh * i_hbbl < cs%harm_BL_val) then

1813 hvel(k) = h_harm

1814 dz_vel(k) = dz_harm(k)

1815 else

1816 z2_wt = 1.

1817 if (zh * i_hbbl < 2. * cs%harm_BL_val) &

1818 z2_wt = max(0., min(1., zh * i_hbbl * i_valbl - 1.))

1819

1820 ! TODO: should z_clear be used here?

1821 z2 = z2_wt * (max(zh, max(zcol, zcol_p1) + dmin) * i_hbbl)

1822 botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)

1823

1824 hvel(k) = (1. - botfn) * h_arith + botfn * h_harm

1825 dz_vel(k) = (1. - botfn) * dz_arith + botfn * dz_harm(k)

1826 endif

1827 endif

1828 endif

1829

1830 if (cs%use_GL90_in_SSW) then

1831 ! The following block calculates the normalized height above the GL90 BBL

1832 ! (z_i_gl90), using a harmonic mean between layer thicknesses. For the

1833 ! GL90 BBL we use simply a constant (Hbbl_gl90). The purpose is that the

1834 ! GL90 coupling coefficient is zeroed out within Hbbl_gl90, to ensure

1835 ! that no momentum gets fluxed into vanished layers. The scheme is not

1836 ! sensitive to the exact value of Hbbl_gl90, as long as it is in a

1837 ! reasonable range (~1-20 m): large enough to capture vanished layers

1838 ! over topography, small enough to not contaminate the interior.

1839

1840 z_i_gl90(k) = z_i_gl90(k+1) + dz_harm(k) * i_hbbl_gl90

1841 endif

1842 enddo

1843

1844 call find_coupling_coef(a_cpl, dz_vel, i, j, dz_harm, bbl_thick, kv_bbl, z_i, &

1845 h_ml, dt, g, gv, us, cs, visc, ustar_2d, tv, work_on_u=.false., obc=obc)

1846

1847 if (allocated(hml_v)) hml_v(i,j) = h_ml

1848

1849 if (cs%use_GL90_in_SSW) then

1850 call find_coupling_coef_gl90(a_cpl_gl90, dz_vel, i, j, z_i_gl90, g, gv, &

1851 cs, varmix, work_on_u=.false.)

1852 endif

1853

1854 do_any_shelf = .false.

1855 if (associated(forces%frac_shelf_v)) then

1856 cs%a1_shelf_v(i,j) = 0.

1857 do_any_shelf = forces%frac_shelf_v(i,j) > 0.

1858

1859 if (do_any_shelf) then

1860 ! Initialize non-harmonic depths

1861 if (.not. cs%harmonic_visc) then

1862 zh = 0.

1863 ztop_min = min(zcol, zcol_p1)

1864 i_htbl = 1. / (visc%tbl_thick_shelf_v(i,j) + dz_neglect)

1865 endif

1866

1867 do k=1,nz

1868 if (cs%harmonic_visc) then

1869 hvel_shelf(k) = hvel(k)

1870 dz_vel_shelf(k) = dz_vel(k)

1871 else

1872 ! Find upwind-biased thickness near the surface.

1873 ! Perhaps this needs to be done more carefully, via find_eta.

1874 h_harm = 2. * h(i,j,k) * h(i,j+1,k) &

1875 / (h(i,j,k) + h(i,j+1,k) + h_neglect)

1876 h_arith = 0.5 * (h(i,j+1,k) + h(i,j,k))

1877 h_delta = h(i,j+1,k) - h(i,j,k)

1878 dz_arith = 0.5 * (dz(i,j+1,k) + dz(i,j,k))

1879

1880 ! Project thickness outward across OBCs using a zero-gradient condition.

1881 if (associated(obc)) then

1882 if (obc%v_N_OBCs_on_PE) then

1883 if (obc%segnum_v(i,j) > 0) then

1884 h_harm = h(i,j,k)

1885 h_arith = h(i,j,k)

1886 h_delta = 0.

1887 dz_arith = dz(i,j,k)

1888 endif

1889 endif

1890

1891 if (obc%v_S_OBCs_on_PE) then

1892 if (obc%segnum_v(i,j) < 0) then

1893 h_harm = h(i,j+1,k)

1894 h_arith = h(i,j+1,k)

1895 h_delta = 0.

1896 dz_arith = dz(i,j+1,k)

1897 endif

1898 endif

1899 endif

1900

1901 zcol = zcol - dz(i,j,k)

1902 zcol_p1 = zcol_p1 - dz(i,j+1,k)

1903

1904 zh = zh + dz_harm(k)

1905

1906 hvel_shelf(k) = hvel(k)

1907 dz_vel_shelf(k) = dz_vel(k)

1908

1909 if (v(i,j,k) * h_delta > 0.) then

1910 if (zh * i_htbl < cs%harm_BL_val) then

1911 hvel_shelf(k) = min(hvel(k), h_harm)

1912 dz_vel_shelf(k) = min(dz_vel(k), dz_harm(k))

1913 else

1914 z2_wt = 1.

1915 if (zh * i_htbl < 2. * cs%harm_BL_val) &

1916 z2_wt = max(0., min(1., zh * i_htbl * i_valbl - 1.))

1917

1918 z2 = z2_wt * (max(zh, ztop_min - min(zcol, zcol_p1)) * i_htbl)

1919 ! TODO: Replace **6

1920 topfn = 1. / (1. + 0.09 * z2**6)

1921

1922 hvel_shelf(k) = min(hvel(k), (1. - topfn) * h_arith + topfn * h_harm)

1923 dz_vel_shelf(k) = min(dz_vel(k), (1. - topfn) * dz_arith + topfn * dz_harm(k))

1924 endif

1925 endif

1926 endif

1927 enddo

1928

1929 call find_coupling_coef(a_shelf, dz_vel_shelf, i, j, dz_harm, &

1930 bbl_thick, kv_bbl, z_i, h_ml, dt, g, gv, us, cs, visc, ustar_2d, &

1931 tv, work_on_u=.false., obc=obc, shelf=.true.)

1932

1933 cs%a1_shelf_v(i,j) = a_shelf(1)

1934 endif

1935 endif

1936

1937 if (do_any_shelf) then

1938 if (cs%use_GL90_in_SSW) then

1939 do k=1,nz+1

1940 cs%a_v(i,j,k) = min(a_cpl_max, (forces%frac_shelf_v(i,j) * a_shelf(k) + &

1941 (1. - forces%frac_shelf_v(i,j)) * a_cpl(k)) + a_cpl_gl90(k))

1942

1943 ! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH

1944 ! CS%a_v(I,j,K) = min(a_cpl_max, forces%frac_shelf_v(I,j) * max(a_shelf(K), a_cpl(K)) + &

1945 ! (1. - forces%frac_shelf_v(I,j)) * a_cpl(K))

1946

1947 cs%a_v_gl90(i,j,k) = min(a_cpl_max, a_cpl_gl90(k))

1948 enddo

1949 else

1950 do k=1,nz+1

1951 cs%a_v(i,j,k) = min(a_cpl_max, (forces%frac_shelf_v(i,j) * a_shelf(k) + &

1952 (1. - forces%frac_shelf_v(i,j)) * a_cpl(k)))

1953 ! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH

1954 ! CS%a_v(I,j,K) = min(a_cpl_max, forces%frac_shelf_v(I,j) * max(a_shelf(K), a_cpl(K)) + &

1955 ! (1. - forces%frac_shelf_v(I,j)) * a_cpl(K))

1956 enddo

1957 endif

1958

1959 do k=1,nz

1960 ! Should we instead take the inverse of the average of the inverses?

1961 cs%h_v(i,j,k) = forces%frac_shelf_v(i,j) * hvel_shelf(k) &

1962 + (1. - forces%frac_shelf_v(i,j)) * hvel(k) + h_neglect

1963 enddo

1964 else

1965 if (cs%use_GL90_in_SSW) then

1966 do k=1,nz+1

1967 a_cpl(k) = a_cpl(k) + a_cpl_gl90(k)

1968 enddo

1969

1970 do k=1,nz+1

1971 cs%a_v_gl90(i,j,k) = min(a_cpl_max, a_cpl_gl90(k))

1972 enddo

1973 endif

1974

1975 do k=1,nz+1

1976 cs%a_v(i,j,k) = min(a_cpl_max, a_cpl(k))

1977 enddo

1978

1979 do k=1,nz

1980 cs%h_v(i,j,k) = hvel(k) + h_neglect

1981 enddo

1982 endif

1983

1984 ! Diagnose total Kv at v-points

1985 if (cs%id_Kv_v > 0) then

1986 do k=1,nz

1987 kv_v(i,j,k) = 0.5 * (cs%a_v(i,j,k) + cs%a_v(i,j,k+1)) * cs%h_v(i,j,k)

1988 enddo

1989 endif

1990

1991 ! Diagnose GL90 Kv at v-points

1992 if (cs%id_Kv_gl90_v > 0) then

1993 do k=1,nz

1994 kv_gl90_v(i,j,k) = 0.5 * (cs%a_v_gl90(i,j,k) + cs%a_v_gl90(i,j,k+1)) * cs%h_v(i,j,k)

1995 enddo

1996 endif

1997 endif ; enddo ; enddo

1998

1999 if (cs%debug) then

2000 call uvchksum("vertvisc_coef h_[uv]", cs%h_u, cs%h_v, g%HI, haloshift=0, &

2001 unscale=gv%H_to_m, scalar_pair=.true.)

2002 call uvchksum("vertvisc_coef a_[uv]", cs%a_u, cs%a_v, g%HI, haloshift=0, &

2003 unscale=gv%H_to_m*us%s_to_T, scalar_pair=.true.)

2004 if (allocated(hml_u) .and. allocated(hml_v)) &

2005 call uvchksum("vertvisc_coef hML_[uv]", hml_u, hml_v, g%HI, &

2006 haloshift=0, unscale=us%Z_to_m, scalar_pair=.true.)

2007 endif

2008

2009! Offer diagnostic fields for averaging.

2010 if (query_averaging_enabled(cs%diag)) then

2011 if (associated(visc%Kv_slow) .and. (cs%id_Kv_slow > 0)) &

2012 call post_data(cs%id_Kv_slow, visc%Kv_slow, cs%diag)

2013 if (cs%id_Kv_u > 0) call post_data(cs%id_Kv_u, kv_u, cs%diag)

2014 if (cs%id_Kv_v > 0) call post_data(cs%id_Kv_v, kv_v, cs%diag)

2015 if (cs%id_Kv_gl90_u > 0) call post_data(cs%id_Kv_gl90_u, kv_gl90_u, cs%diag)

2016 if (cs%id_Kv_gl90_v > 0) call post_data(cs%id_Kv_gl90_v, kv_gl90_v, cs%diag)

2017 if (cs%id_au_vv > 0) call post_data(cs%id_au_vv, cs%a_u, cs%diag)

2018 if (cs%id_av_vv > 0) call post_data(cs%id_av_vv, cs%a_v, cs%diag)

2019 if (cs%id_au_gl90_vv > 0) call post_data(cs%id_au_gl90_vv, cs%a_u_gl90, cs%diag)

2020 if (cs%id_av_gl90_vv > 0) call post_data(cs%id_av_gl90_vv, cs%a_v_gl90, cs%diag)

2021 if (cs%id_h_u > 0) call post_data(cs%id_h_u, cs%h_u, cs%diag)

2022 if (cs%id_h_v > 0) call post_data(cs%id_h_v, cs%h_v, cs%diag)

2023 if (cs%id_hML_u > 0) call post_data(cs%id_hML_u, hml_u, cs%diag)

2024 if (cs%id_hML_v > 0) call post_data(cs%id_hML_v, hml_v, cs%diag)

2025 endif

2026

2027 if (allocated(hml_u)) deallocate(hml_u)

2028 if (allocated(hml_v)) deallocate(hml_v)

2029

2030end subroutine vertvisc_coef

2031

2032

2033!> Calculate the 'coupling coefficient' (a_cpl) at the interfaces.

2034!! If BOTTOMDRAGLAW is defined, the minimum of Hbbl and half the adjacent

2035!! layer thicknesses are used to calculate a_cpl near the bottom.

2036subroutine find_coupling_coef(a_cpl, hvel, i, j, h_harm, bbl_thick, kv_bbl, z_i, h_ml, &

2037 dt, G, GV, US, CS, visc, Ustar_2d, tv, work_on_u, OBC, shelf)

2038 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure

2039 type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure

2040 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

2041 real, dimension(SZK_(GV)+1), &

2042 intent(out) :: a_cpl !< Coupling coefficient across interfaces [H T-1 ~> m s-1 or Pa s m-1]

2043 real, dimension(SZK_(GV)), &

2044 intent(in) :: hvel !< Distance between interfaces at velocity points [Z ~> m]

2045 integer, intent(in) :: i !< Column i-index

2046 integer, intent(in) :: j !< Column j-index

2047 real, dimension(SZK_(GV)), &

2048 intent(in) :: h_harm !< Harmonic mean of thicknesses around a velocity

2049 !! grid point [Z ~> m]

2050 real, intent(in) :: bbl_thick !< Bottom boundary layer thickness [Z ~> m]

2051 real, intent(in) :: kv_bbl !< Bottom boundary layer viscosity, exclusive of

2052 !! any depth-dependent contributions from

2053 !! visc%Kv_shear [H Z T-1 ~> m2 s-1 or Pa s]

2054 real, dimension(SZK_(GV)+1), &

2055 intent(in) :: z_i !< Estimate of interface heights above the bottom,

2056 !! normalized by the bottom boundary layer thickness [nondim]

2057 real, intent(out) :: h_ml !< Mixed layer depth [Z ~> m]

2058 real, intent(in) :: dt !< Time increment [T ~> s]

2059 type(vertvisc_cs), intent(in) :: CS !< Vertical viscosity control structure

2060 type(vertvisc_type), intent(in) :: visc !< Structure containing viscosities and bottom drag

2061 real, dimension(SZI_(G),SZJ_(G)), &

2062 intent(in) :: Ustar_2d !< The wind friction velocity, calculated using

2063 !! the Boussinesq reference density or the

2064 !! time-evolving surface density in non-Boussinesq

2065 !! mode [Z T-1 ~> m s-1]

2066 type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available

2067 !! thermodynamic fields.

2068 logical, intent(in) :: work_on_u !< If true, u-points are being calculated,

2069 !! otherwise they are v-points

2070 type(ocean_obc_type), pointer :: OBC !< Open boundary condition structure

2071 logical, optional, intent(in) :: shelf !< If present and true, use a surface boundary

2072 !! condition appropriate for an ice shelf.

2073

2074 ! Local variables

2075

2076 real :: &

2077 u_star, & ! ustar at a velocity point [Z T-1 ~> m s-1]

2078 tau_mag, & ! The magnitude of the wind stress at a velocity point including gustiness [H Z T-2 ~> m2 s-2 or Pa]

2079 absf, & ! The average of the neighboring absolute values of f [T-1 ~> s-1].

2080 rho_av1, & ! The harmonic mean surface layer density at velocity points [R ~> kg m-3]

2081 z_t, & ! The distance from the top, sometimes normalized

2082 ! by Hmix, [Z ~> m] or [nondim].

2083 kv_tbl, & ! The viscosity in a top boundary layer under ice [H Z T-1 ~> m2 s-1 or Pa s]

2084 tbl_thick, &! The thickness of the top boundary layer [Z ~> m]

2085 kv_add, & ! A viscosity to add [H Z T-1 ~> m2 s-1 or Pa s]

2086 kv_tot ! The total viscosity at an interface [H Z T-1 ~> m2 s-1 or Pa s]

2087 integer :: &

2088 nk_in_ml ! The index of the deepest interface in the mixed layer.

2089 real :: h_shear ! The distance over which shears occur [Z ~> m].

2090 real :: dhc ! The distance between the center of adjacent layers [Z ~> m].

2091 real :: visc_ml ! The mixed layer viscosity [H Z T-1 ~> m2 s-1 or Pa s].

2092 real :: I_Hmix ! The inverse of the mixed layer thickness [Z-1 ~> m-1].

2093 real :: a_ml ! The layer coupling coefficient across an interface in

2094 ! the mixed layer [H T-1 ~> m s-1 or Pa s m-1].

2095 real :: a_floor ! A lower bound on the layer coupling coefficient across an interface in

2096 ! the mixed layer [H T-1 ~> m s-1 or Pa s m-1].

2097 real :: I_amax ! The inverse of the maximum coupling coefficient [T H-1 ~> s m-1 or s m2 kg-1].

2098 real :: temp1 ! A temporary variable [Z2 ~> m2]

2099 real :: ustar2_denom ! A temporary variable in the surface boundary layer turbulence

2100 ! calculations [H Z-1 T-1 ~> s-1 or kg m-3 s-1]

2101 real :: h_neglect ! A vertical distance that is so small it is usually lost

2102 ! in roundoff and can be neglected [Z ~> m].

2103 real :: z2 ! A copy of z_i [nondim]

2104 real :: botfn ! A function that is 1 at the bottom and small far from it [nondim]

2105 real :: topfn ! A function that is 1 at the top and small far from it [nondim]

2106 real :: kv_top ! A viscosity associated with the top boundary layer [H Z T-1 ~> m2 s-1 or Pa s]

2107 logical :: do_shelf, do_OBCs, can_exit

2108 integer :: k

2109 integer :: nz, max_nk

2110

2111 nz = gv%ke

2112

2113 h_neglect = gv%dZ_subroundoff

2114

2115 if (cs%answer_date < 20190101) then

2116 ! The maximum coupling coefficient was originally introduced to avoid

2117 ! truncation error problems in the tridiagonal solver. Effectively, the 1e-10

2118 ! sets the maximum coupling coefficient increment to 1e10 m per timestep.

2119 i_amax = (1.0e-10*gv%H_to_m) * dt

2120 else

2121 i_amax = 0.0

2122 endif

2123

2124 do_shelf = .false. ; if (present(shelf)) do_shelf = shelf

2125

2126 do_obcs = .false.

2127 if (associated(obc)) then

2128 if (work_on_u) then

2129 do_obcs = obc%u_E_OBCs_on_PE .or. obc%u_W_OBCs_on_PE

2130 else

2131 do_obcs = obc%v_N_OBCs_on_PE .or. obc%v_S_OBCs_on_PE

2132 endif

2133 endif

2134

2135 a_cpl(:) = 0.

2136 h_ml = 0.

2137

2138 if (cs%Kvml_invZ2 > 0. .and. .not. do_shelf) then

2139 i_hmix = 1. / (cs%Hmix + h_neglect)

2140 z_t = h_neglect * i_hmix

2141 endif

2142

2143 do k=2,nz

2144 kv_tot = cs%Kv

2145

2146 if (cs%Kvml_invZ2 > 0. .and. .not. do_shelf) then

2147 ! This is an older (vintage ~1997) way to prevent wind stresses from driving very

2148 ! large flows in nearly massless near-surface layers when there is not a physically-

2149 ! based surface boundary layer parameterization. It does not have a plausible

2150 ! physical basis, and probably should not be used.

2151 z_t = z_t + h_harm(k-1) * i_hmix

2152 kv_tot = cs%Kv + cs%Kvml_invZ2 / ((z_t * z_t) * &

2153 (1. + 0.09 * z_t * z_t * z_t * z_t * z_t * z_t))

2154 endif

2155

2156 if (associated(visc%Kv_shear)) then

2157 ! Add in viscosities that are determined by physical processes that are handled in

2158 ! other modules, and which do not respond immediately to the changing layer thicknesses.

2159 ! These processes may include shear-driven mixing or contributions from some boundary

2160 ! layer turbulence schemes. Other viscosity contributions that respond to the evolving

2161 ! layer thicknesses or the surface wind stresses are added later.

2162 if (work_on_u) then

2163 kv_add = 0.5 * (visc%Kv_shear(i,j,k) + visc%Kv_shear(i+1,j,k))

2164

2165 if (do_obcs) then

2166 if (obc%u_E_OBCs_on_PE) then

2167 if (obc%segnum_u(i,j) > 0) then

2168 kv_add = visc%Kv_shear(i,j,k)

2169 endif

2170 endif

2171

2172 if (obc%u_W_OBCs_on_PE) then

2173 if (obc%segnum_u(i,j) < 0) then

2174 kv_add = visc%Kv_shear(i+1,j,k)

2175 endif

2176 endif

2177 endif

2178

2179 kv_tot = kv_tot + kv_add

2180 else

2181 kv_add = 0.5 * (visc%Kv_shear(i,j,k) + visc%Kv_shear(i,j+1,k))

2182

2183 if (do_obcs) then

2184 if (obc%v_N_OBCs_on_PE) then

2185 if (obc%segnum_v(i,j) > 0) then

2186 kv_add = visc%Kv_shear(i,j,k)

2187 endif

2188 endif

2189

2190 if (obc%v_S_OBCs_on_PE) then

2191 if (obc%segnum_v(i,j) < 0) then

2192 kv_add = visc%Kv_shear(i,j+1,k)

2193 endif

2194 endif

2195 endif

2196

2197 kv_tot = kv_tot + kv_add

2198 endif

2199 endif

2200

2201 if (associated(visc%Kv_shear_Bu)) then

2202 ! This is similar to what was done above, but for contributions coming from the corner

2203 ! (vorticity) points. Because OBCs run through the faces and corners there is no need

2204 ! to further modify these viscosities here to take OBCs into account.

2205 if (work_on_u) then

2206 kv_tot = kv_tot + 0.5 * (visc%Kv_shear_Bu(i,j-1,k) + visc%Kv_shear_Bu(i,j,k))

2207 else

2208 kv_tot = kv_tot + 0.5 * (visc%Kv_shear_Bu(i-1,j,k) + visc%Kv_shear_Bu(i,j,k))

2209 endif

2210 endif

2211

2212 ! Set the viscous coupling coefficients, excluding surface mixed layer contributions

2213 ! for now, but including viscous bottom drag, working up from the bottom.

2214 if (cs%bottomdraglaw) then

2215 ! botfn determines when a point is within the influence of the bottom

2216 ! boundary layer, going from 1 at the bottom to 0 in the interior.

2217 z2 = z_i(k)

2218 botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)

2219

2220 kv_tot = kv_tot + (kv_bbl - cs%Kv) * botfn

2221 dhc = 0.5 * (hvel(k) + hvel(k-1))

2222 if (dhc > bbl_thick) then

2223 h_shear = ((1. - botfn) * dhc + botfn * bbl_thick) + h_neglect

2224 else

2225 h_shear = dhc + h_neglect

2226 endif

2227

2228 ! Calculate the coupling coefficients from the viscosities.

2229 a_cpl(k) = kv_tot / (h_shear + (i_amax * kv_tot))

2230 elseif (abs(cs%Kv_extra_bbl) > 0.0) then

2231 ! There is a simple enhancement of the near-bottom viscosities, but no

2232 ! adjustment of the viscous coupling length scales to give a particular

2233 ! bottom stress.

2234

2235 ! botfn determines when a point is within the influence of the bottom

2236 ! boundary layer, going from 1 at the bottom to 0 in the interior.

2237 z2 = z_i(k)

2238 botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)

2239

2240 kv_tot = kv_tot + cs%Kv_extra_bbl * botfn

2241 h_shear = 0.5 * (hvel(k) + hvel(k-1) + h_neglect)

2242

2243 ! Calculate the coupling coefficients from the viscosities.

2244 a_cpl(k) = kv_tot / (h_shear + i_amax * kv_tot)

2245 else

2246 ! Any near-bottom viscous enhancements were already incorporated into

2247 ! Kv_tot, and there is no adjustment of the viscous coupling length

2248 ! scales to give a particular bottom stress.

2249

2250 h_shear = 0.5 * (hvel(k) + hvel(k-1) + h_neglect)

2251 ! Calculate the coupling coefficients from the viscosities.

2252 a_cpl(k) = kv_tot / (h_shear + i_amax * kv_tot)

2253 endif

2254 enddo

2255

2256 ! Assign the bottom coupling coefficients

2257 if (cs%bottomdraglaw) then

2258 dhc = hvel(nz) * 0.5

2259 a_cpl(nz+1) = kv_bbl / ((min(dhc, bbl_thick) + h_neglect) + i_amax * kv_bbl)

2260 elseif (abs(cs%Kv_extra_bbl) > 0.0) then

2261 a_cpl(nz+1) = (cs%Kv + cs%Kv_extra_bbl) &

2262 / ((0.5 * hvel(nz) + h_neglect) + i_amax * (cs%Kv + cs%Kv_extra_bbl))

2263 else

2264 a_cpl(nz+1) = cs%Kv / ((0.5 * hvel(nz) + h_neglect) + i_amax * cs%Kv)

2265 endif

2266

2267 ! Add surface intensified viscous coupling, either as a no-slip boundary condition under a

2268 ! rigid ice-shelf, or due to wind-stress driven surface boundary layer mixing that has not

2269 ! already been added via visc%Kv_shear.

2270 if (do_shelf) then

2271 ! Set the coefficients to include the no-slip surface stress.

2272 if (work_on_u) then

2273 kv_tbl = visc%Kv_tbl_shelf_u(i,j)

2274 tbl_thick = visc%tbl_thick_shelf_u(i,j) + h_neglect

2275 else

2276 kv_tbl = visc%Kv_tbl_shelf_v(i,j)

2277 tbl_thick = visc%tbl_thick_shelf_v(i,j) + h_neglect

2278 endif

2279

2280 z_t = 0.0

2281

2282 ! If a_cpl(1) were not already 0, it would be added here.

2283 if (0.5 * hvel(1) > tbl_thick) then

2284 a_cpl(1) = kv_tbl / (tbl_thick + i_amax * kv_tbl)

2285 else

2286 a_cpl(1) = kv_tbl / ((0.5 * hvel(1) + h_neglect) + i_amax * kv_tbl)

2287 endif

2288

2289 do k=2,nz

2290 z_t = z_t + hvel(k-1) / tbl_thick

2291 topfn = 1. / (1. + 0.09 * z_t**6)

2292

2293 dhc = 0.5 * (hvel(k) + hvel(k-1))

2294 if (dhc > tbl_thick) then

2295 h_shear = ((1. - topfn) * dhc + topfn * tbl_thick) + h_neglect

2296 else

2297 h_shear = dhc + h_neglect

2298 endif

2299

2300 kv_top = topfn * kv_tbl

2301 a_cpl(k) = a_cpl(k) + kv_top / (h_shear + i_amax * kv_top)

2302 enddo

2303 elseif (cs%dynamic_viscous_ML .or. (gv%nkml>0) .or. cs%fixed_LOTW_ML .or. cs%apply_LOTW_floor) then

2304

2305 ! Find the friction velocity and the absolute value of the Coriolis parameter at this point.

2306 u_star = 0. ! Zero out the friction velocity on land points.

2307 tau_mag = 0. ! Zero out the friction velocity on land points.

2308

2309 if (allocated(tv%SpV_avg)) then

2310 rho_av1 = 0.

2311

2312 if (work_on_u) then

2313 u_star = 0.5 * (ustar_2d(i,j) + ustar_2d(i+1,j))

2314 rho_av1 = 2. / (tv%SpV_avg(i,j,1) + tv%SpV_avg(i+1,j,1))

2315 absf = 0.5 * (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i,j)))

2316

2317 if (do_obcs) then

2318 if (obc%u_E_OBCs_on_PE) then

2319 if (obc%segnum_u(i,j) > 0) then

2320 u_star = ustar_2d(i,j)

2321 rho_av1 = 1. / tv%SpV_avg(i,j,1)

2322 endif

2323 endif

2324

2325 if (obc%u_W_OBCs_on_PE) then

2326 if (obc%segnum_u(i,j) < 0) then

2327 u_star = ustar_2d(i+1,j)

2328 rho_av1 = 1. / tv%SpV_avg(i+1,j,1)

2329 endif

2330 endif

2331 endif

2332 else

2333 u_star = 0.5 * (ustar_2d(i,j) + ustar_2d(i,j+1))

2334 rho_av1 = 2. / (tv%SpV_avg(i,j,1) + tv%SpV_avg(i,j+1,1))

2335 absf = 0.5 * (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))

2336

2337 if (do_obcs) then

2338 if (obc%v_N_OBCs_on_PE) then

2339 if (obc%segnum_v(i,j) > 0) then

2340 u_star = ustar_2d(i,j)

2341 rho_av1 = 1. / tv%SpV_avg(i,j,1)

2342 endif

2343 endif

2344

2345 if (obc%v_S_OBCs_on_PE) then

2346 if (obc%segnum_v(i,j) < 0) then

2347 u_star = ustar_2d(i,j+1)

2348 rho_av1 = 1. / tv%SpV_avg(i,j+1,1)

2349 endif

2350 endif

2351 endif

2352 endif

2353

2354 tau_mag = gv%RZ_to_H * rho_av1 * u_star**2

2355 else ! (.not.allocated(tv%SpV_avg))

2356 if (work_on_u) then

2357 u_star = 0.5 * (ustar_2d(i,j) + ustar_2d(i+1,j))

2358 absf = 0.5 * (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i,j)))

2359

2360 if (do_obcs) then

2361 if (obc%u_E_OBCs_on_PE) then

2362 if (obc%segnum_u(i,j) > 0) then

2363 u_star = ustar_2d(i,j)

2364 endif

2365 endif

2366

2367 if (obc%u_W_OBCs_on_PE) then

2368 if (obc%segnum_u(i,j) < 0) then

2369 u_star = ustar_2d(i+1,j)

2370 endif

2371 endif

2372 endif

2373 else

2374 u_star = 0.5 * (ustar_2d(i,j) + ustar_2d(i,j+1))

2375 absf = 0.5 * (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))

2376

2377 if (do_obcs) then

2378 if (obc%v_N_OBCs_on_PE) then

2379 if (obc%segnum_v(i,j) > 0) then

2380 u_star = ustar_2d(i,j)

2381 endif

2382 endif

2383

2384 if (obc%v_S_OBCs_on_PE) then

2385 if (obc%segnum_v(i,j) < 0) then

2386 u_star = ustar_2d(i,j+1)

2387 endif

2388 endif

2389 endif

2390 endif

2391

2392 tau_mag = gv%Z_to_H * u_star**2

2393 endif

2394

2395 ! Determine the thickness of the surface ocean boundary layer and its extent in index space.

2396 nk_in_ml = 0

2397 if (cs%dynamic_viscous_ML) then

2398 ! The fractional number of layers that are within the viscous boundary layer were

2399 ! previously stored in visc%nkml_visc_[uv].

2400 h_ml = h_neglect

2401 max_nk = 0

2402

2403 if (work_on_u) then

2404 nk_in_ml = ceiling(visc%nkml_visc_u(i,j))

2405 max_nk = max(max_nk, nk_in_ml)

2406

2407 do k=1,max_nk

2408 if (k <= visc%nkml_visc_u(i,j)) then ! This layer is all in the ML.

2409 h_ml = h_ml + hvel(k)

2410 elseif (k < visc%nkml_visc_u(i,j) + 1.) then ! Part of this layer is in the ML.

2411 h_ml = h_ml + ((visc%nkml_visc_u(i,j) + 1.) - k) * hvel(k)

2412 endif

2413 enddo

2414 else

2415 nk_in_ml = ceiling(visc%nkml_visc_v(i,j))

2416 max_nk = max(max_nk, nk_in_ml)

2417

2418 do k=1,max_nk

2419 if (k <= visc%nkml_visc_v(i,j)) then ! This layer is all in the ML.

2420 h_ml = h_ml + hvel(k)

2421 elseif (k < visc%nkml_visc_v(i,j) + 1.) then ! Part of this layer is in the ML.

2422 h_ml = h_ml + ((visc%nkml_visc_v(i,j) + 1.) - k) * hvel(k)

2423 endif

2424 enddo

2425 endif

2426 elseif (gv%nkml>0) then

2427 ! This is a simple application of a refined-bulk mixed layer with GV%nkml sublayers.

2428 max_nk = gv%nkml

2429 nk_in_ml = gv%nkml

2430

2431 h_ml = h_neglect

2432

2433 do k=1,gv%nkml

2434 h_ml = h_ml + hvel(k)

2435 enddo

2436 elseif (cs%fixed_LOTW_ML .or. cs%apply_LOTW_floor) then

2437 ! Determine which interfaces are within CS%Hmix of the surface, and set the viscous

2438 ! boundary layer thickness to the smaller of CS%Hmix and the depth of the ocean.

2439 h_ml = 0.0

2440 do k=1,nz

2441 can_exit = .true.

2442 if (h_ml < cs%Hmix) then

2443 nk_in_ml = k

2444

2445 if (h_ml + hvel(k) < cs%Hmix) then

2446 h_ml = h_ml + hvel(k)

2447 can_exit = .false. ! Part of the next deeper layer is also in the mixed layer.

2448 else

2449 h_ml = cs%Hmix

2450 endif

2451 endif

2452

2453 if (can_exit) exit ! All remaining layers in this row are below the mixed layer depth.

2454 enddo

2455

2456 max_nk = max(0, nk_in_ml)

2457 endif

2458

2459 ! Avoid working on columns where the viscous coupling could not be increased.

2460 if (u_star <= 0.) nk_in_ml = 0

2461

2462 ! Set the viscous coupling at the interfaces as the larger of what was previously

2463 ! set and the contributions from the surface boundary layer.

2464 z_t = 0.

2465 if (cs%apply_LOTW_floor .and. &

2466 (cs%dynamic_viscous_ML .or. gv%nkml > 0 .or. cs%fixed_LOTW_ML)) then

2467 do k=2,max_nk

2468 if (k <= nk_in_ml) then

2469 z_t = z_t + hvel(k-1)

2470

2471 ! The viscosity in visc_ml is set to go to 0 at the mixed layer top and bottom

2472 ! (in a log-layer) and be further limited by rotation to give the natural Ekman length.

2473 temp1 = (z_t * h_ml - z_t * z_t)

2474 if (gv%Boussinesq) then

2475 ustar2_denom = (cs%vonKar * gv%Z_to_H * u_star**2) &

2476 / (absf * temp1 + (h_ml + h_neglect) * u_star)

2477 else

2478 ustar2_denom = (cs%vonKar * tau_mag) &

2479 / (absf * temp1 + (h_ml + h_neglect) * u_star)

2480 endif

2481

2482 visc_ml = temp1 * ustar2_denom

2483 ! Set the viscous coupling based on the model's vertical resolution. The omission of

2484 ! the I_amax factor here is consistent with answer dates above 20190101.

2485 a_ml = visc_ml / (0.25 * (hvel(k) + hvel(k-1) + h_neglect))

2486

2487 ! As a floor on the viscous coupling, assume that the length scale in the denominator can

2488 ! not be larger than the distance from the surface, consistent with a logarithmic velocity

2489 ! profile. This is consistent with visc_ml, but cancels out common factors of z_t.

2490 a_floor = (h_ml - z_t) * ustar2_denom

2491

2492 ! Choose the largest estimate of a_cpl.

2493 a_cpl(k) = max(a_cpl(k), a_ml, a_floor)

2494 ! An option could be added to change this to: a_cpl(i,K) = max(a_cpl(i,K) + a_ml, a_floor)

2495 endif

2496 enddo

2497 elseif (cs%apply_LOTW_floor) then

2498 do k=2,max_nk

2499 if (k <= nk_in_ml) then

2500 z_t = z_t + hvel(k-1)

2501

2502 temp1 = (z_t * h_ml - z_t * z_t)

2503 if (gv%Boussinesq) then

2504 ustar2_denom = (cs%vonKar * gv%Z_to_H * u_star**2) &

2505 / (absf * temp1 + (h_ml + h_neglect) * u_star)

2506 else

2507 ustar2_denom = (cs%vonKar * tau_mag) &

2508 / (absf * temp1 + (h_ml + h_neglect) * u_star)

2509 endif

2510

2511 ! As a floor on the viscous coupling, assume that the length scale in the denominator can not

2512 ! be larger than the distance from the surface, consistent with a logarithmic velocity profile.

2513 a_cpl(k) = max(a_cpl(k), (h_ml - z_t) * ustar2_denom)

2514 endif

2515 enddo

2516 else

2517 do k=2,max_nk

2518 if (k <= nk_in_ml) then

2519 z_t = z_t + hvel(k-1)

2520

2521 temp1 = (z_t * h_ml - z_t * z_t)

2522 ! This viscosity is set to go to 0 at the mixed layer top and bottom (in a log-layer)

2523 ! and be further limited by rotation to give the natural Ekman length.

2524 ! The following expressions are mathematically equivalent.

2525 if (gv%Boussinesq .or. (cs%answer_date < 20230601)) then

2526 visc_ml = u_star * cs%vonKar * (gv%Z_to_H * temp1 * u_star) &

2527 / (absf * temp1 + (h_ml + h_neglect) * u_star)

2528 else

2529 visc_ml = cs%vonKar * (temp1 * tau_mag) &

2530 / (absf * temp1 + (h_ml + h_neglect) * u_star)

2531 endif

2532 a_ml = visc_ml / (0.25 * (hvel(k) + hvel(k-1) + h_neglect) + 0.5 * i_amax * visc_ml)

2533

2534 ! Choose the largest estimate of a_cpl, but these could be changed to be additive.

2535 a_cpl(k) = max(a_cpl(k), a_ml)

2536 ! An option could be added to change this to: a_cpl(i,K) = a_cpl(i,K) + a_ml

2537 endif

2538 enddo

2539 endif

2540 endif

2541end subroutine find_coupling_coef

2542

2543

2544!> Velocity components which exceed a threshold for physically reasonable values are truncated,

2545!! and the running sum of the number of trunctionas within the non-symmetric memory computational

2546!! domain is incremented. Optionally, any column with excessive velocities may be sent

2547!! to a diagnostic reporting subroutine.

2548subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS)

2549 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

2550 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

2551 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

2552 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &

2553 intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]

2554 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &

2555 intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]

2556 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &

2557 intent(in) :: h !< Layer thickness [H ~> m or kg m-2]

2558 type(accel_diag_ptrs), intent(in) :: adp !< Acceleration diagnostic pointers

2559 type(cont_diag_ptrs), intent(in) :: cdp !< Continuity diagnostic pointers

2560 type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces

2561 type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag

2562 real, intent(in) :: dt !< Time increment [T ~> s]

2563 type(vertvisc_cs), pointer :: cs !< Vertical viscosity control structure

2564

2565 ! Local variables

2566 real :: cfl ! The local CFL number [nondim]

2567 real :: h_report ! A thickness below which not to report truncations [H ~> m or kg m-2]

2568 real :: vel_report(szib_(g),szjb_(g)) ! The velocity to report [L T-1 ~> m s-1]

2569 real :: u_old(szib_(g),szj_(g),szk_(gv)) ! The previous u-velocity [L T-1 ~> m s-1]

2570 real :: v_old(szi_(g),szjb_(g),szk_(gv)) ! The previous v-velocity [L T-1 ~> m s-1]

2571 logical :: trunc_any, dowrite(szib_(g),szjb_(g))

2572 logical :: do_any_write

2573 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz

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

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

2576

2577 h_report = 3.0 * gv%Angstrom_H

2578

2579 if (len_trim(cs%u_trunc_file) > 0) then

2580 do_any_write = .false.

2581 trunc_any = .false.

2582

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

2584 dowrite(i,j) = .false.

2585 vel_report(i,j) = 3.0e8 * us%m_s_to_L_T

2586 enddo ; enddo

2587

2588 do k=1,nz ; do j=js,je ; do i=isq,ieq

2589 if (abs(u(i,j,k)) < cs%vel_underflow) u(i,j,k) = 0.0

2590 if (u(i,j,k) < 0.0) then

2591 cfl = (-u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))

2592 else

2593 cfl = (u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))

2594 endif

2595 if (cfl > cs%CFL_trunc) trunc_any = .true.

2596 if (cfl > cs%CFL_report) then

2597 dowrite(i,j) = .true.

2598 do_any_write = .true.

2599 vel_report(i,j) = min(vel_report(i,j), abs(u(i,j,k)))

2600 endif

2601 enddo ; enddo ; enddo

2602

2603 do j=js,je ; do i=isq,ieq ; if (dowrite(i,j)) then

2604 u_old(i,j,:) = u(i,j,:)

2605 endif ; enddo ; enddo

2606

2607 if (trunc_any) then

2608 do k=1,nz ; do j=js,je ; do i=isq,ieq

2609 if ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL_trunc) then

2610 u(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i+1,j) / (dt * g%dy_Cu(i,j)))

2611 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2612 (cs%h_u(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2613 elseif ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then

2614 u(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dy_Cu(i,j)))

2615 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2616 (cs%h_u(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2617 endif

2618 enddo ; enddo ; enddo

2619 endif

2620

2621 if (do_any_write) then

2622 do j=js,je ; do i=isq,ieq ; if (dowrite(i,j)) then

2623 ! Call a diagnostic reporting subroutines are called if unphysically large values are found.

2624 call write_u_accel(i, j, u_old, h, adp, cdp, dt, g, gv, us, cs%PointAccel_CSp, &

2625 vel_report(i,j), forces%taux(i,j), a=cs%a_u, hv=cs%h_u)

2626 endif ; enddo ; enddo

2627 endif

2628 else ! Do not report accelerations leading to large velocities.

2629 do k=1,nz ; do j=js,je ; do i=isq,ieq

2630 if (abs(u(i,j,k)) < cs%vel_underflow) then ; u(i,j,k) = 0.0

2631 elseif ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL_trunc) then

2632 u(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i+1,j) / (dt * g%dy_Cu(i,j)))

2633 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2634 (cs%h_u(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2635 elseif ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then

2636 u(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dy_Cu(i,j)))

2637 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2638 (cs%h_u(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2639 endif

2640 enddo ; enddo ; enddo

2641 endif

2642

2643 if (len_trim(cs%v_trunc_file) > 0) then

2644 do_any_write =.false.

2645 trunc_any = .false.

2646

2647

2648 do j=jsq,jeq ; do i=is,ie

2649 dowrite(i,j) = .false.

2650 vel_report(i,j) = 3.0e8 * us%m_s_to_L_T

2651 enddo ; enddo

2652

2653 do k=1,nz ; do j=jsq,jeq ; do i=is,ie

2654 if (abs(v(i,j,k)) < cs%vel_underflow) v(i,j,k) = 0.0

2655 if (v(i,j,k) < 0.0) then

2656 cfl = (-v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))

2657 else

2658 cfl = (v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))

2659 endif

2660 if (cfl > cs%CFL_trunc) trunc_any = .true.

2661 if (cfl > cs%CFL_report) then

2662 dowrite(i,j) = .true.

2663 do_any_write = .true.

2664 vel_report(i,j) = min(vel_report(i,j), abs(v(i,j,k)))

2665 endif

2666 enddo ; enddo ; enddo

2667

2668 do j=jsq,jeq ; do i=is,ie ; if (dowrite(i,j)) then

2669 v_old(i,j,:) = v(i,j,:)

2670 endif ; enddo ; enddo

2671

2672 if (trunc_any) then

2673 do k=1,nz ; do j=jsq,jeq ; do i=is,ie

2674 if ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL_trunc) then

2675 v(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i,j+1) / (dt * g%dx_Cv(i,j)))

2676 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2677 (cs%h_v(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2678 elseif ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then

2679 v(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dx_Cv(i,j)))

2680 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2681 (cs%h_v(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2682 endif

2683 enddo ; enddo ; enddo

2684 endif

2685

2686 if (do_any_write) then

2687 do j=jsq,jeq ; do i=is,ie ; if (dowrite(i,j)) then

2688 ! Call a diagnostic reporting subroutines are called if unphysically large values are found.

2689 call write_v_accel(i, j, v_old, h, adp, cdp, dt, g, gv, us, cs%PointAccel_CSp, &

2690 vel_report(i,j), forces%tauy(i,j), a=cs%a_v, hv=cs%h_v)

2691 endif ; enddo ; enddo

2692 endif

2693 else ! Do not report accelerations leading to large velocities.

2694 do k=1,nz ; do j=jsq,jeq ; do i=is,ie

2695 if (abs(v(i,j,k)) < cs%vel_underflow) then ; v(i,j,k) = 0.0

2696 elseif ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL_trunc) then

2697 v(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i,j+1) / (dt * g%dx_Cv(i,j)))

2698 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2699 (cs%h_v(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2700 elseif ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then

2701 v(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dx_Cv(i,j)))

2702 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2703 (cs%h_v(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2704 endif

2705 enddo ; enddo ; enddo

2706 endif

2707

2708end subroutine vertvisc_limit_vel

2709

2710

2711!> Initialize the vertical friction module

2712subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &

2713 ntrunc, CS, fpmix)

2714 type(ocean_internal_state), &

2715 target, intent(in) :: mis !< The "MOM Internal State", a set of pointers

2716 !! to the fields and accelerations that make

2717 !! up the ocean's physical state

2718 type(time_type), target, intent(in) :: time !< Current model time

2719 type(ocean_grid_type), intent(in) :: g !< Ocean grid structure

2720 type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid structure

2721 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

2722 type(param_file_type), intent(in) :: param_file !< File to parse for parameters

2723 type(diag_ctrl), target, intent(inout) :: diag !< Diagnostic control structure

2724 type(accel_diag_ptrs), intent(inout) :: adp !< Acceleration diagnostic pointers

2725 type(directories), intent(in) :: dirs !< Relevant directory paths

2726 integer, target, intent(inout) :: ntrunc !< Number of velocity truncations

2727 type(vertvisc_cs), pointer :: cs !< Vertical viscosity control structure

2728 logical, optional, intent(in) :: fpmix !< Nonlocal momentum mixing

2729

2730 ! Local variables

2731

2732 real :: kv_bbl ! A viscosity in the bottom boundary layer with a simple scheme [H Z T-1 ~> m2 s-1 or Pa s]

2733 real :: kv_back_z ! A background kinematic viscosity [Z2 T-1 ~> m2 s-1]

2734 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.

2735 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz

2736 logical :: lfpmix

2737 character(len=200) :: kappa_gl90_file, inputdir, kdgl90_varname

2738 ! This include declares and sets the variable "version".

2739# include "version_variable.h"

2740 character(len=40) :: mdl = "MOM_vert_friction" ! This module's name.

2741 character(len=40) :: thickness_units

2742 real :: kv_mks ! KVML in MKS [m2 s-1]

2743

2744 if (associated(cs)) then

2745 call mom_error(warning, "vertvisc_init called with an associated "// &

2746 "control structure.")

2747 return

2748 endif

2749 allocate(cs)

2750

2751 cs%initialized = .true.

2752

2753 if (gv%Boussinesq) then ; thickness_units = "m"

2754 else ; thickness_units = "kg m-2" ; endif

2755

2756 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke

2757 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

2758

2759 cs%diag => diag ; cs%ntrunc => ntrunc ; ntrunc = 0

2760

2761 lfpmix = .false.

2762 if (present(fpmix)) lfpmix = fpmix

2763

2764! Default, read and log parameters

2765 call log_version(param_file, mdl, version, "", log_to_all=.true., debugging=.true.)

2766 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &

2767 "This sets the default value for the various _ANSWER_DATE parameters.", &

2768 default=99991231)

2769 call get_param(param_file, mdl, "VERT_FRICTION_ANSWER_DATE", cs%answer_date, &

2770 "The vintage of the order of arithmetic and expressions in the viscous "//&

2771 "calculations. Values below 20190101 recover the answers from the end of 2018, "//&

2772 "while higher values use expressions that do not use an arbitrary hard-coded "//&

2773 "maximum viscous coupling coefficient between layers. Values below 20230601 "//&

2774 "recover a form of the viscosity within the mixed layer that breaks up the "//&

2775 "magnitude of the wind stress in some non-Boussinesq cases.", &

2776 default=default_answer_date, do_not_log=.not.gv%Boussinesq)

2777 if (.not.gv%Boussinesq) cs%answer_date = max(cs%answer_date, 20230701)

2778

2779 call get_param(param_file, mdl, "BOTTOMDRAGLAW", cs%bottomdraglaw, &

2780 "If true, the bottom stress is calculated with a drag "//&

2781 "law of the form c_drag*|u|*u. The velocity magnitude "//&

2782 "may be an assumed value or it may be based on the "//&

2783 "actual velocity in the bottommost HBBL, depending on "//&

2784 "LINEAR_DRAG.", default=.true.)

2785 call get_param(param_file, mdl, "DIRECT_STRESS", cs%direct_stress, &

2786 "If true, the wind stress is distributed over the topmost HMIX_STRESS of fluid "//&

2787 "(like in HYCOM), and an added mixed layer viscosity or a physically based "//&

2788 "boundary layer turbulence parameterization is not needed for stability.", &

2789 default=.false.)

2790 call get_param(param_file, mdl, "DYNAMIC_VISCOUS_ML", cs%dynamic_viscous_ML, &

2791 "If true, use a bulk Richardson number criterion to "//&

2792 "determine the mixed layer thickness for viscosity.", &

2793 default=.false.)

2794 call get_param(param_file, mdl, "FIXED_DEPTH_LOTW_ML", cs%fixed_LOTW_ML, &

2795 "If true, use a Law-of-the-wall prescription for the mixed layer viscosity "//&

2796 "within a boundary layer that is the lesser of HMIX_FIXED and the total "//&

2797 "depth of the ocean in a column.", default=.false.)

2798 call get_param(param_file, mdl, "LOTW_VISCOUS_ML_FLOOR", cs%apply_LOTW_floor, &

2799 "If true, use a Law-of-the-wall prescription to set a lower bound on the "//&

2800 "viscous coupling between layers within the surface boundary layer, based "//&

2801 "the distance of interfaces from the surface. This only acts when there "//&

2802 "are large changes in the thicknesses of successive layers or when the "//&

2803 "viscosity is set externally and the wind stress has subsequently increased.", &

2804 default=.false.)

2805 call get_param(param_file, mdl, 'VON_KARMAN_CONST', cs%vonKar, &

2806 'The value the von Karman constant as used for mixed layer viscosity.', &

2807 units='nondim', default=0.41)

2808 call get_param(param_file, mdl, "U_TRUNC_FILE", cs%u_trunc_file, &

2809 "The absolute path to a file into which the accelerations "//&

2810 "leading to zonal velocity truncations are written. "//&

2811 "Undefine this for efficiency if this diagnostic is not needed.", &

2812 default=" ", debuggingparam=.true.)

2813 call get_param(param_file, mdl, "V_TRUNC_FILE", cs%v_trunc_file, &

2814 "The absolute path to a file into which the accelerations "//&

2815 "leading to meridional velocity truncations are written. "//&

2816 "Undefine this for efficiency if this diagnostic is not needed.", &

2817 default=" ", debuggingparam=.true.)

2818 call get_param(param_file, mdl, "HARMONIC_VISC", cs%harmonic_visc, &

2819 "If true, use the harmonic mean thicknesses for "//&

2820 "calculating the vertical viscosity.", default=.false.)

2821 call get_param(param_file, mdl, "HARMONIC_BL_SCALE", cs%harm_BL_val, &

2822 "A scale to determine when water is in the boundary "//&

2823 "layers based solely on harmonic mean thicknesses for "//&

2824 "the purpose of determining the extent to which the "//&

2825 "thicknesses used in the viscosities are upwinded.", &

2826 default=0.0, units="nondim")

2827 call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)

2828

2829 if (gv%nkml < 1) then

2830 call get_param(param_file, mdl, "HMIX_FIXED", cs%Hmix, &

2831 "The prescribed depth over which the near-surface viscosity and "//&

2832 "diffusivity are elevated when the bulk mixed layer is not used.", &

2833 units="m", scale=us%m_to_Z, fail_if_missing=.true.)

2834 endif

2835 if (cs%direct_stress) then

2836 if (gv%nkml < 1) then

2837 call get_param(param_file, mdl, "HMIX_STRESS", cs%Hmix_stress, &

2838 "The depth over which the wind stress is applied if DIRECT_STRESS is true.", &

2839 units="m", default=us%Z_to_m*cs%Hmix, scale=gv%m_to_H)

2840 else

2841 call get_param(param_file, mdl, "HMIX_STRESS", cs%Hmix_stress, &

2842 "The depth over which the wind stress is applied if DIRECT_STRESS is true.", &

2843 units="m", fail_if_missing=.true., scale=gv%m_to_H)

2844 endif

2845 if (cs%Hmix_stress <= 0.0) call mom_error(fatal, "vertvisc_init: " // &

2846 "HMIX_STRESS must be set to a positive value if DIRECT_STRESS is true.")

2847 endif

2848 call get_param(param_file, mdl, "KV", kv_back_z, &

2849 "The background kinematic viscosity in the interior. "//&

2850 "The molecular value, ~1e-6 m2 s-1, may be used.", &

2851 units="m2 s-1", fail_if_missing=.true., scale=us%m2_s_to_Z2_T)

2852 ! Convert input kinematic viscosity to dynamic viscosity when non-Boussinesq.

2853 cs%Kv = (us%Z2_T_to_m2_s*gv%m2_s_to_HZ_T) * kv_back_z

2854

2855 call get_param(param_file, mdl, "USE_GL90_IN_SSW", cs%use_GL90_in_SSW, &

2856 "If true, use simpler method to calculate 1/N^2 in GL90 vertical "// &

2857 "viscosity coefficient. This method is valid in stacked shallow water mode.", &

2858 default=.false.)

2859 call get_param(param_file, mdl, "KD_GL90", cs%kappa_gl90, &

2860 "The scalar diffusivity used in GL90 vertical viscosity scheme.", &

2861 units="m2 s-1", default=0.0, scale=us%m_to_L*us%Z_to_L*gv%m_to_H*us%T_to_s, &

2862 do_not_log=.not.cs%use_GL90_in_SSW)

2863 call get_param(param_file, mdl, "READ_KD_GL90", cs%read_kappa_gl90, &

2864 "If true, read a file (given by KD_GL90_FILE) containing the "//&

2865 "spatially varying diffusivity KD_GL90 used in the GL90 scheme.", default=.false., &

2866 do_not_log=.not.cs%use_GL90_in_SSW)

2867 if (cs%read_kappa_gl90) then

2868 if (cs%kappa_gl90 > 0) then

2869 call mom_error(fatal, "MOM_vert_friction.F90, vertvisc_init: KD_GL90 > 0 "// &

2870 "is not compatible with READ_KD_GL90 = .TRUE. ")

2871 endif

2872 call get_param(param_file, mdl, "INPUTDIR", inputdir, &

2873 "The directory in which all input files are found.", &

2874 default=".", do_not_log=.true.)

2875 inputdir = slasher(inputdir)

2876 call get_param(param_file, mdl, "KD_GL90_FILE", kappa_gl90_file, &

2877 "The file containing the spatially varying diffusivity used in the "// &

2878 "GL90 scheme.", default="kd_gl90.nc", do_not_log=.not.cs%use_GL90_in_SSW)

2879 call get_param(param_file, mdl, "KD_GL90_VARIABLE", kdgl90_varname, &

2880 "The name of the GL90 diffusivity variable to read "//&

2881 "from KD_GL90_FILE.", default="kd_gl90", do_not_log=.not.cs%use_GL90_in_SSW)

2882 kappa_gl90_file = trim(inputdir) // trim(kappa_gl90_file)

2883

2884 allocate(cs%kappa_gl90_2d(g%isd:g%ied, g%jsd:g%jed), source=0.0)

2885 call mom_read_data(kappa_gl90_file, kdgl90_varname, cs%kappa_gl90_2d(:,:), g%domain, &

2886 scale=us%m_to_L*us%Z_to_L*gv%m_to_H*us%T_to_s)

2887 call pass_var(cs%kappa_gl90_2d, g%domain)

2888 endif

2889 call get_param(param_file, mdl, "USE_GL90_N2", cs%use_GL90_N2, &

2890 "If true, use GL90 vertical viscosity coefficient that is depth-independent; "// &

2891 "this corresponds to a kappa_GM that scales as N^2 with depth.", &

2892 default=.false., do_not_log=.not.cs%use_GL90_in_SSW)

2893 if (cs%use_GL90_N2) then

2894 if (.not. cs%use_GL90_in_SSW) call mom_error(fatal, &

2895 "MOM_vert_friction.F90, vertvisc_init: "//&

2896 "When USE_GL90_N2=True, USE_GL90_in_SSW must also be True.")

2897 if (cs%kappa_gl90 > 0) then

2898 call mom_error(fatal, "MOM_vert_friction.F90, vertvisc_init: KD_GL90 > 0 "// &

2899 "is not compatible with USE_GL90_N2 = .TRUE. ")

2900 endif

2901 if (cs%read_kappa_gl90) call mom_error(fatal, &

2902 "MOM_vert_friction.F90, vertvisc_init: "//&

2903 "READ_KD_GL90 = .TRUE. is not compatible with USE_GL90_N2 = .TRUE.")

2904 call get_param(param_file, mdl, "alpha_GL90", cs%alpha_gl90, &

2905 "Coefficient used to compute a depth-independent GL90 vertical "//&

2906 "viscosity via Kv_GL90 = alpha_GL90 * f2. Is only used "// &

2907 "if USE_GL90_N2 is true. Note that the implied Kv_GL90 "// &

2908 "corresponds to a KD_GL90 that scales as N^2 with depth.", &

2909 units="m2 s", default=0.0, scale=gv%m_to_H*us%m_to_Z*us%s_to_T, &

2910 do_not_log=.not.cs%use_GL90_in_SSW)

2911 endif

2912 call get_param(param_file, mdl, "HBBL_GL90", cs%Hbbl_gl90, &

2913 "The thickness of the GL90 bottom boundary layer, "//&

2914 "which defines the range over which the GL90 coupling "//&

2915 "coefficient is zeroed out, in order to avoid fluxing "//&

2916 "momentum into vanished layers over steep topography.", &

2917 units="m", default=5.0, scale=us%m_to_Z, do_not_log=.not.cs%use_GL90_in_SSW)

2918

2919 cs%Kvml_invZ2 = 0.0

2920 if (gv%nkml < 1) then

2921 call get_param(param_file, mdl, "KVML", kv_mks, &

2922 "The scale for an extra kinematic viscosity in the mixed layer", &

2923 units="m2 s-1", default=-1.0, do_not_log=.true.)

2924 if (kv_mks >= 0.0) then

2925 call mom_error(warning, "KVML is a deprecated parameter. Use KV_ML_INVZ2 instead.")

2926 else

2927 kv_mks = 0.0

2928 endif

2929 call get_param(param_file, mdl, "KV_ML_INVZ2", cs%Kvml_invZ2, &

2930 "An extra kinematic viscosity in a mixed layer of thickness HMIX_FIXED, "//&

2931 "with the actual viscosity scaling as 1/(z*HMIX_FIXED)^2, where z is the "//&

2932 "distance from the surface, to allow for finite wind stresses to be "//&

2933 "transmitted through infinitesimally thin surface layers. This is an "//&

2934 "older option for numerical convenience without a strong physical basis, "//&

2935 "and its use is now discouraged.", &

2936 units="m2 s-1", default=kv_mks, scale=gv%m2_s_to_HZ_T)

2937 endif

2938

2939 if (.not.cs%bottomdraglaw) then

2940 call get_param(param_file, mdl, "KV_EXTRA_BBL", cs%Kv_extra_bbl, &

2941 "An extra kinematic viscosity in the benthic boundary layer. "//&

2942 "KV_EXTRA_BBL is not used if BOTTOMDRAGLAW is true.", &

2943 units="m2 s-1", default=0.0, scale=gv%m2_s_to_HZ_T, do_not_log=.true.)

2944 if (cs%Kv_extra_bbl == 0.0) then

2945 call get_param(param_file, mdl, "KVBBL", kv_bbl, &

2946 "An extra kinematic viscosity in the benthic boundary layer. "//&

2947 "KV_EXTRA_BBL is not used if BOTTOMDRAGLAW is true.", &

2948 units="m2 s-1", default=us%Z2_T_to_m2_s*kv_back_z, scale=gv%m2_s_to_HZ_T, &

2949 do_not_log=.true.)

2950 if (abs(kv_bbl - cs%Kv) > 1.0e-15*abs(cs%Kv)) then

2951 call mom_error(warning, "KVBBL is a deprecated parameter. Use KV_EXTRA_BBL instead.")

2952 cs%Kv_extra_bbl = kv_bbl - cs%Kv

2953 endif

2954 endif

2955 call log_param(param_file, mdl, "KV_EXTRA_BBL", cs%Kv_extra_bbl, &

2956 "An extra kinematic viscosity in the benthic boundary layer. "//&

2957 "KV_EXTRA_BBL is not used if BOTTOMDRAGLAW is true.", &

2958 units="m2 s-1", default=0.0, unscale=gv%HZ_T_to_m2_s)

2959 endif

2960 call get_param(param_file, mdl, "HBBL", cs%Hbbl, &

2961 "The thickness of a bottom boundary layer with a viscosity increased by "//&

2962 "KV_EXTRA_BBL if BOTTOMDRAGLAW is not defined, or the thickness over which "//&

2963 "near-bottom velocities are averaged for the drag law if BOTTOMDRAGLAW is "//&

2964 "defined but LINEAR_DRAG is not.", &

2965 units="m", fail_if_missing=.true., scale=us%m_to_Z)

2966 call get_param(param_file, mdl, "CFL_TRUNCATE", cs%CFL_trunc, &

2967 "The value of the CFL number that will cause velocity "//&

2968 "components to be truncated; instability can occur past 0.5.", &

2969 units="nondim", default=0.5)

2970 call get_param(param_file, mdl, "CFL_REPORT", cs%CFL_report, &

2971 "The value of the CFL number that causes accelerations "//&

2972 "to be reported; the default is CFL_TRUNCATE.", &

2973 units="nondim", default=cs%CFL_trunc)

2974 call get_param(param_file, mdl, "CFL_TRUNCATE_RAMP_TIME", cs%truncRampTime, &

2975 "The time over which the CFL truncation value is ramped "//&

2976 "up at the beginning of the run.", &

2977 units="s", default=0., scale=us%s_to_T)

2978 cs%CFL_truncE = cs%CFL_trunc

2979 call get_param(param_file, mdl, "CFL_TRUNCATE_START", cs%CFL_truncS, &

2980 "The start value of the truncation CFL number used when "//&

2981 "ramping up CFL_TRUNC.", &

2982 units="nondim", default=0.)

2983 call get_param(param_file, mdl, "STOKES_MIXING_COMBINED", cs%StokesMixing, &

2984 "Flag to use Stokes drift Mixing via the Lagrangian "//&

2985 " current (Eulerian plus Stokes drift). "//&

2986 " Still needs work and testing, so not recommended for use.",&

2987 default=.false.)

2988 !BGR 04/04/2018{

2989 ! StokesMixing is required for MOM6 for some Langmuir mixing parameterization.

2990 ! The code used here has not been developed for vanishing layers or in

2991 ! conjunction with any bottom friction. Therefore, the following line is

2992 ! added so this functionality cannot be used without user intervention in

2993 ! the code. This will prevent general use of this functionality until proper

2994 ! care is given to the previously mentioned issues. Comment out the following

2995 ! MOM_error to use, but do so at your own risk and with these points in mind.

2996 !}

2997 if (cs%StokesMixing) then

2998 call mom_error(fatal, "Stokes mixing requires user intervention in the code.\n"//&

2999 " Model now exiting. See MOM_vert_friction.F90 for \n"//&

3000 " details (search 'BGR 04/04/2018' to locate comment).")

3001 endif

3002 call get_param(param_file, mdl, "VEL_UNDERFLOW", cs%vel_underflow, &

3003 "A negligibly small velocity magnitude below which velocity "//&

3004 "components are set to 0. A reasonable value might be "//&

3005 "1e-30 m/s, which is less than an Angstrom divided by "//&

3006 "the age of the universe.", units="m s-1", default=0.0, scale=us%m_s_to_L_T)

3007

3008 alloc_(cs%a_u(isdb:iedb,jsd:jed,nz+1)) ; cs%a_u(:,:,:) = 0.0

3009 alloc_(cs%a_u_gl90(isdb:iedb,jsd:jed,nz+1)) ; cs%a_u_gl90(:,:,:) = 0.0

3010 alloc_(cs%h_u(isdb:iedb,jsd:jed,nz)) ; cs%h_u(:,:,:) = 0.0

3011 alloc_(cs%a_v(isd:ied,jsdb:jedb,nz+1)) ; cs%a_v(:,:,:) = 0.0

3012 alloc_(cs%a_v_gl90(isd:ied,jsdb:jedb,nz+1)) ; cs%a_v_gl90(:,:,:) = 0.0

3013 alloc_(cs%h_v(isd:ied,jsdb:jedb,nz)) ; cs%h_v(:,:,:) = 0.0

3014

3015 cs%id_Kv_slow = register_diag_field('ocean_model', 'Kv_slow', diag%axesTi, time, &

3016 'Slow varying vertical viscosity', 'm2 s-1', conversion=gv%HZ_T_to_m2_s)

3017

3018 cs%id_Kv_u = register_diag_field('ocean_model', 'Kv_u', diag%axesCuL, time, &

3019 'Total vertical viscosity at u-points', 'm2 s-1', conversion=gv%H_to_m**2*us%s_to_T)

3020

3021 cs%id_Kv_v = register_diag_field('ocean_model', 'Kv_v', diag%axesCvL, time, &

3022 'Total vertical viscosity at v-points', 'm2 s-1', conversion=gv%H_to_m**2*us%s_to_T)

3023

3024 cs%id_Kv_gl90_u = register_diag_field('ocean_model', 'Kv_gl90_u', diag%axesCuL, time, &

3025 'GL90 vertical viscosity at u-points', 'm2 s-1', conversion=gv%H_to_m**2*us%s_to_T)

3026

3027 cs%id_Kv_gl90_v = register_diag_field('ocean_model', 'Kv_gl90_v', diag%axesCvL, time, &

3028 'GL90 vertical viscosity at v-points', 'm2 s-1', conversion=gv%H_to_m**2*us%s_to_T)

3029

3030 cs%id_au_vv = register_diag_field('ocean_model', 'au_visc', diag%axesCui, time, &

3031 'Zonal Viscous Vertical Coupling Coefficient', 'm s-1', conversion=gv%H_to_m*us%s_to_T)

3032

3033 cs%id_av_vv = register_diag_field('ocean_model', 'av_visc', diag%axesCvi, time, &

3034 'Meridional Viscous Vertical Coupling Coefficient', 'm s-1', conversion=gv%H_to_m*us%s_to_T)

3035

3036 cs%id_au_gl90_vv = register_diag_field('ocean_model', 'au_gl90_visc', diag%axesCui, time, &

3037 'Zonal Viscous Vertical GL90 Coupling Coefficient', 'm s-1', conversion=gv%H_to_m*us%s_to_T)

3038

3039 cs%id_av_gl90_vv = register_diag_field('ocean_model', 'av_gl90_visc', diag%axesCvi, time, &

3040 'Meridional Viscous Vertical GL90 Coupling Coefficient', 'm s-1', conversion=gv%H_to_m*us%s_to_T)

3041

3042 cs%id_h_u = register_diag_field('ocean_model', 'Hu_visc', diag%axesCuL, time, &

3043 'Thickness at Zonal Velocity Points for Viscosity', &

3044 thickness_units, conversion=gv%H_to_MKS)

3045 ! Alternately, to always give this variable in 'm' use the following line instead:

3046 ! 'm', conversion=GV%H_to_m)

3047

3048 cs%id_h_v = register_diag_field('ocean_model', 'Hv_visc', diag%axesCvL, time, &

3049 'Thickness at Meridional Velocity Points for Viscosity', &

3050 thickness_units, conversion=gv%H_to_MKS)

3051

3052 cs%id_hML_u = register_diag_field('ocean_model', 'HMLu_visc', diag%axesCu1, time, &

3053 'Mixed Layer Thickness at Zonal Velocity Points for Viscosity', &

3054 thickness_units, conversion=us%Z_to_m)

3055

3056 cs%id_hML_v = register_diag_field('ocean_model', 'HMLv_visc', diag%axesCv1, time, &

3057 'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', &

3058 thickness_units, conversion=us%Z_to_m)

3059

3060 if (lfpmix) then

3061 cs%id_uE_h = register_diag_field('ocean_model', 'uE_h' , cs%diag%axesTL, &

3062 time, 'x-zonal Eulerian' , 'm s-1', conversion=us%L_T_to_m_s)

3063 cs%id_vE_h = register_diag_field('ocean_model', 'vE_h' , cs%diag%axesTL, &

3064 time, 'y-merid Eulerian' , 'm s-1', conversion=us%L_T_to_m_s)

3065 cs%id_uInc_h = register_diag_field('ocean_model','uInc_h',cs%diag%axesTL, &

3066 time, 'x-zonal Eulerian' , 'm s-1', conversion=us%L_T_to_m_s)

3067 cs%id_vInc_h = register_diag_field('ocean_model','vInc_h',cs%diag%axesTL, &

3068 time, 'x-zonal Eulerian' , 'm s-1', conversion=us%L_T_to_m_s)

3069 cs%id_uStk = register_diag_field('ocean_model', 'uStk' , cs%diag%axesTL, &

3070 time, 'x-FP du increment' , 'm s-1', conversion=us%L_T_to_m_s)

3071 cs%id_vStk = register_diag_field('ocean_model', 'vStk' , cs%diag%axesTL, &

3072 time, 'y-FP dv increment' , 'm s-1', conversion=us%L_T_to_m_s)

3073

3074 cs%id_FPtau2s = register_diag_field('ocean_model','Omega_tau2s',cs%diag%axesTi, &

3075 time, 'Stress direction from shear','radians')

3076 cs%id_FPtau2w = register_diag_field('ocean_model','Omega_tau2w',cs%diag%axesTi, &

3077 time, 'Stress direction from wind','radians')

3078 cs%id_uStk0 = register_diag_field('ocean_model', 'uStk0' , diag%axesT1, &

3079 time, 'Zonal Surface Stokes', 'm s-1', conversion=us%L_T_to_m_s)

3080 cs%id_vStk0 = register_diag_field('ocean_model', 'vStk0' , diag%axesT1, &

3081 time, 'Merid Surface Stokes', 'm s-1', conversion=us%L_T_to_m_s)

3082 endif

3083

3084 cs%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', diag%axesCuL, time, &

3085 'Zonal Acceleration from Vertical Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)

3086 if (cs%id_du_dt_visc > 0) call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)

3087 cs%id_dv_dt_visc = register_diag_field('ocean_model', 'dv_dt_visc', diag%axesCvL, time, &

3088 'Meridional Acceleration from Vertical Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)

3089 if (cs%id_dv_dt_visc > 0) call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)

3090 cs%id_GLwork = register_diag_field('ocean_model', 'GLwork', diag%axesTL, time, &

3091 'Sign-definite Kinetic Energy Source from GL90 Vertical Viscosity', &

3092 'm3 s-3', conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)

3093 cs%id_du_dt_visc_gl90 = register_diag_field('ocean_model', 'du_dt_visc_gl90', diag%axesCuL, time, &

3094 'Zonal Acceleration from GL90 Vertical Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)

3095 if ((cs%id_du_dt_visc_gl90 > 0) .or. (cs%id_GLwork > 0)) then

3096 call safe_alloc_ptr(adp%du_dt_visc_gl90,isdb,iedb,jsd,jed,nz)

3097 call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)

3098 endif

3099 cs%id_dv_dt_visc_gl90 = register_diag_field('ocean_model', 'dv_dt_visc_gl90', diag%axesCvL, time, &

3100 'Meridional Acceleration from GL90 Vertical Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)

3101 if ((cs%id_dv_dt_visc_gl90 > 0) .or. (cs%id_GLwork > 0)) then

3102 call safe_alloc_ptr(adp%dv_dt_visc_gl90,isd,ied,jsdb,jedb,nz)

3103 call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)

3104 endif

3105 cs%id_du_dt_str = register_diag_field('ocean_model', 'du_dt_str', diag%axesCuL, time, &

3106 'Zonal Acceleration from Surface Wind Stresses', 'm s-2', conversion=us%L_T2_to_m_s2)

3107 if (cs%id_du_dt_str > 0) call safe_alloc_ptr(adp%du_dt_str,isdb,iedb,jsd,jed,nz)

3108 cs%id_dv_dt_str = register_diag_field('ocean_model', 'dv_dt_str', diag%axesCvL, time, &

3109 'Meridional Acceleration from Surface Wind Stresses', 'm s-2', conversion=us%L_T2_to_m_s2)

3110 if (cs%id_dv_dt_str > 0) call safe_alloc_ptr(adp%dv_dt_str,isd,ied,jsdb,jedb,nz)

3111

3112 cs%id_taux_bot = register_diag_field('ocean_model', 'taux_bot', diag%axesCu1, &

3113 time, 'Zonal Bottom Stress from Ocean to Earth', &

3114 'Pa', conversion=us%RZ_to_kg_m2*us%L_T2_to_m_s2)

3115 cs%id_tauy_bot = register_diag_field('ocean_model', 'tauy_bot', diag%axesCv1, &

3116 time, 'Meridional Bottom Stress from Ocean to Earth', &

3117 'Pa', conversion=us%RZ_to_kg_m2*us%L_T2_to_m_s2)

3118

3119 !CS%id_hf_du_dt_visc = register_diag_field('ocean_model', 'hf_du_dt_visc', diag%axesCuL, Time, &

3120 ! 'Fractional Thickness-weighted Zonal Acceleration from Vertical Viscosity', &

3121 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)

3122 !if (CS%id_hf_du_dt_visc > 0) then

3123 ! call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)

3124 ! call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)

3125 !endif

3126

3127 !CS%id_hf_dv_dt_visc = register_diag_field('ocean_model', 'hf_dv_dt_visc', diag%axesCvL, Time, &

3128 ! 'Fractional Thickness-weighted Meridional Acceleration from Vertical Viscosity', &

3129 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)

3130 !if (CS%id_hf_dv_dt_visc > 0) then

3131 ! call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)

3132 ! call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)

3133 !endif

3134

3135 cs%id_hf_du_dt_visc_2d = register_diag_field('ocean_model', 'hf_du_dt_visc_2d', diag%axesCu1, time, &

3136 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Vertical Viscosity', &

3137 'm s-2', conversion=us%L_T2_to_m_s2)

3138 if (cs%id_hf_du_dt_visc_2d > 0) then

3139 call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)

3140 call safe_alloc_ptr(adp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)

3141 endif

3142

3143 cs%id_hf_dv_dt_visc_2d = register_diag_field('ocean_model', 'hf_dv_dt_visc_2d', diag%axesCv1, time, &

3144 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Vertical Viscosity', &

3145 'm s-2', conversion=us%L_T2_to_m_s2)

3146 if (cs%id_hf_dv_dt_visc_2d > 0) then

3147 call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)

3148 call safe_alloc_ptr(adp%diag_hfrac_v,isd,ied,jsdb,jedb,nz)

3149 endif

3150

3151 cs%id_h_du_dt_visc = register_diag_field('ocean_model', 'h_du_dt_visc', diag%axesCuL, time, &

3152 'Thickness Multiplied Zonal Acceleration from Horizontal Viscosity', &

3153 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

3154 if (cs%id_h_du_dt_visc > 0) then

3155 call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)

3156 call safe_alloc_ptr(adp%diag_hu,isdb,iedb,jsd,jed,nz)

3157 endif

3158

3159 cs%id_h_dv_dt_visc = register_diag_field('ocean_model', 'h_dv_dt_visc', diag%axesCvL, time, &

3160 'Thickness Multiplied Meridional Acceleration from Horizontal Viscosity', &

3161 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

3162 if (cs%id_h_dv_dt_visc > 0) then

3163 call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)

3164 call safe_alloc_ptr(adp%diag_hv,isd,ied,jsdb,jedb,nz)

3165 endif

3166

3167 cs%id_h_du_dt_str = register_diag_field('ocean_model', 'h_du_dt_str', diag%axesCuL, time, &

3168 'Thickness Multiplied Zonal Acceleration from Surface Wind Stresses', &

3169 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

3170 if (cs%id_h_du_dt_str > 0) then

3171 call safe_alloc_ptr(adp%du_dt_str,isdb,iedb,jsd,jed,nz)

3172 call safe_alloc_ptr(adp%diag_hu,isdb,iedb,jsd,jed,nz)

3173 endif

3174

3175 cs%id_h_dv_dt_str = register_diag_field('ocean_model', 'h_dv_dt_str', diag%axesCvL, time, &

3176 'Thickness Multiplied Meridional Acceleration from Surface Wind Stresses', &

3177 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

3178 if (cs%id_h_dv_dt_str > 0) then

3179 call safe_alloc_ptr(adp%dv_dt_str,isd,ied,jsdb,jedb,nz)

3180 call safe_alloc_ptr(adp%diag_hv,isd,ied,jsdb,jedb,nz)

3181 endif

3182

3183 cs%id_du_dt_str_visc_rem = register_diag_field('ocean_model', 'du_dt_str_visc_rem', diag%axesCuL, time, &

3184 'Zonal Acceleration from Surface Wind Stresses multiplied by viscous remnant', &

3185 'm s-2', conversion=us%L_T2_to_m_s2)

3186 if (cs%id_du_dt_str_visc_rem > 0) then

3187 call safe_alloc_ptr(adp%du_dt_str,isdb,iedb,jsd,jed,nz)

3188 call safe_alloc_ptr(adp%visc_rem_u,isdb,iedb,jsd,jed,nz)

3189 endif

3190

3191 cs%id_dv_dt_str_visc_rem = register_diag_field('ocean_model', 'dv_dt_str_visc_rem', diag%axesCvL, time, &

3192 'Meridional Acceleration from Surface Wind Stresses multiplied by viscous remnant', &

3193 'm s-2', conversion=us%L_T2_to_m_s2)

3194 if (cs%id_dv_dt_str_visc_rem > 0) then

3195 call safe_alloc_ptr(adp%dv_dt_str,isd,ied,jsdb,jedb,nz)

3196 call safe_alloc_ptr(adp%visc_rem_v,isd,ied,jsdb,jedb,nz)

3197 endif

3198

3199 if ((len_trim(cs%u_trunc_file) > 0) .or. (len_trim(cs%v_trunc_file) > 0)) &

3200 call pointaccel_init(mis, time, g, param_file, diag, dirs, cs%PointAccel_CSp)

3201

3202end subroutine vertvisc_init

3203

3204!> Update the CFL truncation value as a function of time.

3205!! If called with the optional argument activate=.true., record the

3206!! value of Time as the beginning of the ramp period.

3207subroutine updatecfltruncationvalue(Time, CS, US, activate)

3208 type(time_type), target, intent(in) :: time !< Current model time

3209 type(vertvisc_cs), pointer :: cs !< Vertical viscosity control structure

3210 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type

3211 logical, optional, intent(in) :: activate !< Specify whether to record the value of

3212 !! Time as the beginning of the ramp period

3213

3214 ! Local variables

3215 real :: deltatime ! The time since CS%rampStartTime [T ~> s], which may be negative.

3216 real :: wghta ! The relative weight of the final value [nondim]

3217 character(len=12) :: msg

3218

3219 if (cs%truncRampTime==0.) return ! This indicates to ramping is turned off

3220

3221 ! We use the optional argument to indicate this Time should be recorded as the

3222 ! beginning of the ramp-up period.

3223 if (present(activate)) then

3224 if (activate) then

3225 cs%rampStartTime = time ! Record the current time

3226 cs%CFLrampingIsActivated = .true.

3227 endif

3228 endif

3229 if (.not.cs%CFLrampingIsActivated) return

3230 deltatime = max(0., time_minus_signed(time, cs%rampStartTime, scale=us%s_to_T))

3231 if (deltatime >= cs%truncRampTime) then

3232 cs%CFL_trunc = cs%CFL_truncE

3233 cs%truncRampTime = 0. ! This turns off ramping after this call

3234 else

3235 wghta = min( 1., deltatime / cs%truncRampTime ) ! Linear profile in time

3236 !wghtA = wghtA*wghtA ! Convert linear profile to parabolic profile in time

3237 !wghtA = wghtA*wghtA*(3. - 2.*wghtA) ! Convert linear profile to cosine profile

3238 wghta = 1. - ( (1. - wghta)**2 ) ! Convert linear profile to inverted parabolic profile

3239 cs%CFL_trunc = cs%CFL_truncS + wghta * ( cs%CFL_truncE - cs%CFL_truncS )

3240 endif

3241 write(msg(1:12),'(es12.3)') cs%CFL_trunc

3242 call mom_error(note, "MOM_vert_friction: updateCFLtruncationValue set CFL limit to "//trim(msg))

3243end subroutine updatecfltruncationvalue

3244

3245!> Clean up and deallocate the vertical friction module

3246subroutine vertvisc_end(CS)

3247 type(vertvisc_cs), intent(inout) :: cs !< Vertical viscosity control structure that

3248 !! will be deallocated in this subroutine.

3249

3250 if ((len_trim(cs%u_trunc_file) > 0) .or. (len_trim(cs%v_trunc_file) > 0)) &

3251 deallocate(cs%PointAccel_CSp)

3252

3253 dealloc_(cs%a_u) ; dealloc_(cs%h_u)

3254 dealloc_(cs%a_v) ; dealloc_(cs%h_v)

3255 if (associated(cs%a1_shelf_u)) deallocate(cs%a1_shelf_u)

3256 if (associated(cs%a1_shelf_v)) deallocate(cs%a1_shelf_v)

3257 if (allocated(cs%kappa_gl90_2d)) deallocate(cs%kappa_gl90_2d)

3258end subroutine vertvisc_end

3259

3260!> \namespace mom_vert_friction

3261!! \author Robert Hallberg

3262!! \date April 1994 - October 2006

3263!!

3264!! The vertical diffusion of momentum is fully implicit. This is

3265!! necessary to allow for vanishingly small layers. The coupling

3266!! is based on the distance between the centers of adjacent layers,

3267!! except where a layer is close to the bottom compared with a

3268!! bottom boundary layer thickness when a bottom drag law is used.

3269!! A stress top b.c. and a no slip bottom b.c. are used. There

3270!! is no limit on the time step for vertvisc.

3271!!

3272!! Near the bottom, the horizontal thickness interpolation scheme

3273!! changes to an upwind biased estimate to control the effect of

3274!! spurious Montgomery potential gradients at the bottom where

3275!! nearly massless layers layers ride over the topography. Within a

3276!! few boundary layer depths of the bottom, the harmonic mean

3277!! thickness (i.e. (2 h+ h-) / (h+ + h-) ) is used if the velocity

3278!! is from the thinner side and the arithmetic mean thickness

3279!! (i.e. (h+ + h-)/2) is used if the velocity is from the thicker

3280!! side. Both of these thickness estimates are second order

3281!! accurate. Above this the arithmetic mean thickness is used.

3282!!

3283!! In addition, vertvisc truncates any velocity component that exceeds a

3284!! maximum CFL number to a fraction of this value. This basically keeps

3285!! instabilities spatially localized. The number of times the velocity is

3286!! truncated is reported each time the energies are saved, and if

3287!! exceeds CS%Maxtrunc the model will stop itself and change the time

3288!! to a large value. This has proven very useful in (1) diagnosing

3289!! model failures and (2) letting the model settle down to a

3290!! meaningful integration from a poorly specified initial condition.

3291!!

3292!! The same code is used for the two velocity components, by

3293!! indirectly referencing the velocities and defining a handful of

3294!! direction-specific defined variables.

3295!!

3296!! Macros written all in capital letters are defined in MOM_memory.h.

3297!!

3298!! A small fragment of the grid is shown below:

3299!! \verbatim

3300!! j+1 x ^ x ^ x At x: q

3301!! j+1 > o > o > At ^: v, frhatv, tauy

3302!! j x ^ x ^ x At >: u, frhatu, taux

3303!! j > o > o > At o: h

3304!! j-1 x ^ x ^ x

3305!! i-1 i i+1 At x & ^:

3306!! i i+1 At > & o:

3307!! \endverbatim

3308!!

3309!! The boundaries always run through q grid points (x).

3310end module mom_vert_friction