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) [L Z T-2 ~> m2 s-2]

227 real, dimension(SZI_(G),SZJB_(G)) :: tauy_v !< kinematic merid wind stress (v-pts) [L Z 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 [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 ! initialize arrays

261 ue_h(:,:,:) = 0.0

262 ve_h(:,:,:) = 0.0

263 ue_u(:,:,:) = 0.0

264 ve_v(:,:,:) = 0.0

265 vinc_v(:,:,:) = 0.0

266 uinc_u(:,:,:) = 0.0

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

268

269 ! u-points

270 do j = js,je

271 do i = isq,ieq

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

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

274 ! h to u-pts

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

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

277 depth = 0.

278 gat1 = 0.

279 do k=1, nz

280 ! cell center

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

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

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

284 sigma = depth / hbl_u(i,j)

285 ! cell bottom

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

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

288 ! nonlocal boundary-layer increment

289 uinc_u(i,j,k) = dt * cemp_nl * taux_u(i,j) * dgdsig / (hbl_u(i,j) * gv%H_to_Z)

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

291 else

292 uinc_u(i,j,k) = 0.0

293 endif

294 enddo

295 else

296 do k=1, nz

297 uinc_u(i,j,k) = 0.0

298 enddo

299 endif

300 enddo

301 enddo

302

303 ! v-points

304 do j = jsq,jeq

305 do i = is,ie

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

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

308 ! h to v-pts

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

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

311 depth = 0.

312 gat1 = 0.

313 do k=1, nz

314 ! cell center

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

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

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

318 sigma = depth / hbl_v(i,j)

319 ! cell bottom

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

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

322 ! nonlocal boundary-layer increment

323 vinc_v(i,j,k) = dt * cemp_nl * tauy_v(i,j) * dgdsig / (hbl_v(i,j) * gv%H_to_Z)

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

325 else

326 vinc_v(i,j,k) = 0.0

327 endif

328 enddo

329 else

330 do k=1, nz

331 vinc_v(i,j,k) = 0.0

332 enddo

333 endif

334 enddo

335 enddo

336

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

338 if (lpost) then

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

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

341 ustk = 0.0

342 vstk = 0.0

343 us0 = 0.0

344 vs0 = 0.0

345

346 do j = js,je

347 do i = is,ie

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

349 ! u to h-pts

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

351 ! v to h-pts

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

353 do k = 1,nz

354 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

355 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

356 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

357 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

358 enddo

359 ! Wind, Stress and Shear align at surface

360 omega_tau2w(i,j,:) = 0.0

361 omega_tau2s(i,j,:) = 0.0

362 do k = 1,nz

363 kp1 = min( nz , k+1)

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

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

366 omega_s2x = atan2(dv, du)

367

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

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

370 omega_tau2x = atan2(dv, du)

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

372

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

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

375 omega_tau2w(i,j,kp1) = omega_tmp

376

377 omega_tmp = omega_tau2x - omega_s2x

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

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

380 omega_tau2s(i,j,kp1) = omega_tmp

381

382 enddo

383 endif

384

385 ! Stokes drift

386 do b=1,waves%NumBands

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

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

389 enddo

390 depth = 0.0

391 do k = 1,nz

392 do b = 1, waves%NumBands

393 ! cell center

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

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

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

397 enddo

398 ! cell bottom

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

400 enddo

401 enddo

402 enddo

403

404 ! post FPmix diagnostics

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

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

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

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

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

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

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

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

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

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

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

416

417 endif

418

419end subroutine vertfpmix

420

421

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

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

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

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

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

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

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

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

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

431!! nu = f^2 * alpha

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

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

434!! conditions at the top and bottom.

435!!

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

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

438!! or

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

440

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

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

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

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

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

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

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

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

449 !! bottom, normalized by the GL90 bottom

450 !! boundary layer thickness [nondim]

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

452 !! with GL90 across interfaces; is not

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

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

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

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

457 !! otherwise they are v-points.

458

459 ! local variables

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

461 integer :: k, nz

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

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

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

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

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

467

468 nz = gv%ke

469 h_neglect = gv%dZ_subroundoff

470 kdgl90_use_vert_struct = .false.

471

472 if (varmix%use_variable_mixing) then

473 kdgl90_use_vert_struct = allocated(varmix%kdgl90_struct)

474 endif

475

476 a_cpl_gl90(:) = 0.

477

478 do k=2,nz

479 if (work_on_u) then

480 ! compute coupling coefficient at u-points

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

482 if (cs%use_GL90_N2) then

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

484 else

485 if (cs%read_kappa_gl90) then

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

487 else

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

489 endif

490 if (kdgl90_use_vert_struct) then

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

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

493 endif

494 endif

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

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

497 z2 = z_i(k)

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

499

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

501 else

502 ! compute viscosities at v-points

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

504

505 if (cs%use_GL90_N2) then

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

507 else

508 if (cs%read_kappa_gl90) then

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

510 else

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

512 endif

513 if (kdgl90_use_vert_struct) then

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

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

516 endif

517 endif

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

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

520 z2 = z_i(k)

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

522

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

524 endif

525 enddo

526end subroutine find_coupling_coef_gl90

527

528

529!> Perform a fully implicit vertical diffusion

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

531!!

532!! This is solving the tridiagonal system

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

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

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

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

537!! encompassing background viscosity as well as contributions from

538!! enhanced mixed and bottom layer viscosities.

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

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

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

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

543 taux_bot, tauy_bot, fpmix, Waves)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

558 !! equations for diagnostics

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

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

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

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

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

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

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

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

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

568 type(wave_parameters_cs), &

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

570

571 ! Fields from forces used in this subroutine:

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

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

574

575 ! Local variables

576

577 real :: b1

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

579 real :: c1(szk_(gv))

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

581 real :: d1

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

583 real :: ray

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

585 real :: b_denom_1

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

587

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

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

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

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

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

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

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

595

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

611

612 logical :: dostokesmixing

613 logical :: lfpmix

614

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

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

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

618

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

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

621

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

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

624

625 if (cs%id_GLwork > 0) then

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

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

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

629 if (.not.g%symmetric) &

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

631 endif

632

633 if (cs%direct_stress) then

634 hmix = cs%Hmix_stress

635 i_hmix = 1.0 / hmix

636 endif

637 dt_rho0 = dt / gv%H_to_RZ

638 h_neglect = gv%H_subroundoff

639 idt = 1.0 / dt

640

641 accel_underflow = cs%vel_underflow * idt

642

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

644 dostokesmixing=.false.

645 if (cs%StokesMixing) then

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

647 if (.not. dostokesmixing) &

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

649 endif

650 lfpmix = .false.

651 if ( present(fpmix) ) lfpmix = fpmix

652

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

654 ! tridiagonal solver.

655

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

657 ! includes Stokes.

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

659 if (dostokesmixing) then

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

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

662 endif ; enddo ; enddo ; enddo

663 endif

664

665 if (lfpmix) then

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

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

668 endif ; enddo ; enddo ; enddo

669 endif

670

671 if (associated(adp%du_dt_visc)) then

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

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

674 enddo ; enddo ; enddo

675 endif

676

677 if (associated(adp%du_dt_visc_gl90)) then

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

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

680 enddo ; enddo ; enddo

681 endif

682

683 if (associated(adp%du_dt_str)) then

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

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

686 enddo ; enddo ; enddo

687 endif

688

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

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

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

692 if (cs%direct_stress) then

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

694 surface_stress(i,j) = 0.0

695 zds = 0.0

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

697 do k=1,nz

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

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

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

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

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

703 enddo

704 endif ; enddo ; enddo

705 else

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

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

708 enddo ; enddo

709 endif

710

711 ! perform forward elimination on the tridiagonal system

712 !

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

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

715 !

716 ! ignoring the Rayleigh drag contribution,

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

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

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

720 !

721 ! for forward elimination, we want to:

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

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

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

725 !

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

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

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

729 !

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

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

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

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

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

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

736

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

738 ray = 0.

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

740

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

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

743 d1 = b_denom_1 * b1

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

745

746 if (associated(adp%du_dt_str)) then

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

748 endif

749

750 do k=2,nz

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

752

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

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

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

756 d1 = b_denom_1 * b1

757 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

758

759 if (associated(adp%du_dt_str)) then

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

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

762 endif

763

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

765 if (dt < 0) exit

766 enddo

767

768 if (associated(adp%du_dt_str)) then

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

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

771 endif

772

773 do k=nz-1,1,-1

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

775

776 if (associated(adp%du_dt_str)) then

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

778

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

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

781 endif

782 enddo

783 endif ; enddo ; enddo

784

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

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

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

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

789 if (associated(adp%du_dt_visc_gl90)) then

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

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

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

793 d1 = b_denom_1 * b1

794

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

796

797 do k=2,nz

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

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

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

801 d1 = b_denom_1 * b1

802

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

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

805 enddo

806

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

808 do k=nz-1,1,-1

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

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

811 enddo

812

813 do k=1,nz

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

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

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

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

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

819

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

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

822 endif

823 enddo

824

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

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

827 if (cs%id_GLwork > 0) then

828 do k=1,nz

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

830 enddo

831 endif

832 endif ; enddo ; enddo

833 endif

834 endif

835

836 if (associated(adp%du_dt_visc)) then

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

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

839

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

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

842 enddo ; enddo ; enddo

843 endif

844

845 if (allocated(visc%taux_shelf)) then

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

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

848 enddo ; enddo

849 endif

850

851 if (present(taux_bot)) then

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

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

854 enddo ; enddo

855

856 if (allocated(visc%Ray_u)) then

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

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

859 enddo ; enddo ; enddo

860 endif

861 endif

862

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

864 if (dostokesmixing) then

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

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

867 endif ; enddo ; enddo ; enddo

868 endif

869

870 if (lfpmix) then

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

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

873 endif ; enddo ; enddo ; enddo

874 endif

875

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

877

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

879 if (dostokesmixing) then

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

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

882 endif ; enddo ; enddo ; enddo

883 endif

884

885 if (lfpmix) then

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

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

888 endif ; enddo ; enddo ; enddo

889 endif

890

891 if (associated(adp%dv_dt_visc)) then

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

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

894 enddo ; enddo ; enddo

895 endif

896

897 if (associated(adp%dv_dt_visc_gl90)) then

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

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

900 enddo ; enddo ; enddo

901 endif

902

903 if (associated(adp%dv_dt_str)) then

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

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

906 enddo ; enddo ; enddo

907 endif

908

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

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

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

912 if (cs%direct_stress) then

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

914 surface_stress(i,j) = 0.0

915 zds = 0.0

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

917 do k=1,nz

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

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

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

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

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

923 enddo

924 endif ; enddo ; enddo

925 else

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

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

928 enddo ; enddo

929 endif

930

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

932 ray = 0.

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

934

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

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

937 d1 = b_denom_1 * b1

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

939

940 if (associated(adp%dv_dt_str)) then

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

942 endif

943

944 do k=2,nz

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

946

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

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

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

950 d1 = b_denom_1 * b1

951 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

952

953 if (associated(adp%dv_dt_str)) then

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

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

956 endif

957

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

959 if (dt < 0) exit

960 enddo

961

962 if (associated(adp%dv_dt_str)) then

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

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

965 endif

966

967 do k=nz-1,1,-1

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

969

970 if (associated(adp%dv_dt_str)) then

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

972

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

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

975 endif

976 enddo

977 endif ; enddo ; enddo

978

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

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

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

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

983 if (associated(adp%dv_dt_visc_gl90)) then

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

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

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

987 d1 = b_denom_1 * b1

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

989

990 do k=2,nz

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

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

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

994 d1 = b_denom_1 * b1

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

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

997 enddo

998

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

1000 do k=nz-1,1,-1

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

1002 enddo

1003 endif ; enddo ; enddo

1004

1005 do k=1,nz

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

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

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

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

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

1011

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

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

1014 endif ; enddo ; enddo

1015 enddo

1016

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

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

1019 if (cs%id_GLwork > 0) then

1020 do k=1,nz

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

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

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

1024 endif ; enddo ; enddo

1025 enddo

1026 endif

1027 endif

1028 endif

1029

1030 if (associated(adp%dv_dt_visc)) then

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

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

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

1034 enddo ; enddo ; enddo

1035 endif

1036

1037 if (allocated(visc%tauy_shelf)) then

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

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

1040 enddo ; enddo

1041 endif

1042

1043 if (present(tauy_bot)) then

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

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

1046 enddo ; enddo

1047

1048 if (allocated(visc%Ray_v)) then

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

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

1051 enddo ; enddo ; enddo

1052 endif

1053 endif

1054

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

1056 if (dostokesmixing) then

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

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

1059 endif ; enddo ; enddo ; enddo

1060 endif

1061

1062 if (lfpmix) then

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

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

1065 endif ; enddo ; enddo ; enddo

1066 endif

1067

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

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

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

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

1072 ! barotropic transport [uv]h.

1073 if (cs%id_GLwork > 0) then

1074 if (.not.g%symmetric) &

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

1076 do k=1,nz

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

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

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

1080 enddo ; enddo

1081 enddo

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

1083 endif

1084

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

1086

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

1088 if (associated(obc)) then

1089 do n=1,obc%number_of_segments

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

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

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

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

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

1095 enddo ; enddo

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

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

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

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

1100 enddo ; enddo

1101 endif

1102 endif

1103 enddo

1104 endif

1105

1106 ! Offer diagnostic fields for averaging.

1107 if (query_averaging_enabled(cs%diag)) then

1108 if (cs%id_du_dt_visc > 0) &

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

1110 if (cs%id_du_dt_visc_gl90 > 0) &

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

1112 if (cs%id_dv_dt_visc > 0) &

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

1114 if (cs%id_dv_dt_visc_gl90 > 0) &

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

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

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

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

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

1120 if (cs%id_du_dt_str > 0) &

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

1122 if (cs%id_dv_dt_str > 0) &

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

1124

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

1126 ! Diagnostics of the fractional thicknesses times momentum budget terms

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

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

1129 !if (CS%id_hf_du_dt_visc > 0) &

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

1131 !if (CS%id_hf_dv_dt_visc > 0) &

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

1133

1134 ! Diagnostics for thickness-weighted vertically averaged viscous accelerations

1135 if (cs%id_hf_du_dt_visc_2d > 0) &

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

1137 if (cs%id_hf_dv_dt_visc_2d > 0) &

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

1139

1140 ! Diagnostics for thickness x viscous accelerations

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

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

1143 endif

1144

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

1146 ! Diagnostics for thickness x wind stress accelerations

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

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

1149

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

1151 if (cs%id_du_dt_str_visc_rem > 0) &

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

1153 if (cs%id_dv_dt_str_visc_rem > 0) &

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

1155 endif

1156 endif

1157

1158end subroutine vertvisc

1159

1160

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

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

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

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

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

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

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

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

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

1170 !! barotropic acceleration that a layer experiences after

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

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

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

1174 !! barotropic acceleration that a layer experiences after

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

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

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

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

1179

1180 ! Local variables

1181

1182 real :: b1

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

1184 real :: c1(szk_(gv))

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

1186 real :: d1

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

1188 real :: ray

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

1190 real :: b_denom_1

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

1192

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

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

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

1196

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

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

1199

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

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

1202

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

1204

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

1206 ray = 0.

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

1208

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

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

1211 d1 = b_denom_1 * b1

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

1213

1214 do k=2,nz

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

1216

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

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

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

1220 d1 = b_denom_1 * b1

1221 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

1222

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

1224 if (dt < 0) exit

1225 enddo

1226

1227 do k=nz-1,1,-1

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

1229 enddo

1230 endif ; enddo ; enddo

1231

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

1233

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

1235 ray = 0.

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

1237

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

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

1240 d1 = b_denom_1 * b1

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

1242

1243 do k=2,nz

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

1245

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

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

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

1249 d1 = b_denom_1 * b1

1250 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

1251

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

1253 if (dt < 0) exit

1254 enddo

1255

1256 do k=nz-1,1,-1

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

1258 enddo

1259 endif ; enddo ; enddo

1260

1261 if (cs%debug) then

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

1263 scalar_pair=.true.)

1264 endif

1265end subroutine vertvisc_remnant

1266

1267

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1286 !! thermodynamic fields.

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

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

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

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

1291 ! Field from forces used in this subroutine:

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

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

1294

1295 ! Local variables

1296

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

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

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

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

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

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

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

1304 real :: &

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1319 ! a_cpl_gl90 is part of a_cpl.

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

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

1322 real :: &

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1350 ! than Hbbl into the interior [nondim].

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

1352 ! than Htbl into the interior [nondim].

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

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

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

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

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

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

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

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

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

1362

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

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

1365 ! of the harmonic mean thicknesses [nondim].

1366 logical :: do_any_shelf

1367 integer :: zi_dir

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

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

1370

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

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

1373

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

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

1376

1377 h_neglect = gv%H_subroundoff

1378 dz_neglect = gv%dZ_subroundoff

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

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

1381

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

1383

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

1385

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

1387

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

1389

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

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

1392

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

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

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

1396 endif

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

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

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

1400 endif

1401

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

1403

1404 ! First do u-points

1405

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

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

1408 if (cs%use_GL90_in_SSW) then

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

1410 endif

1411

1412 if (cs%bottomdraglaw) then

1413 kv_bbl = visc%Kv_bbl_u(i,j)

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

1415 i_hbbl = 1. / bbl_thick

1416 endif

1417

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

1419 zi_dir = 0

1420

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

1422 if (associated(obc)) then

1423 if (obc%u_E_OBCs_on_PE) then

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

1425 dmin = g%bathyT(i,j)

1426 zi_dir = -1

1427 endif

1428 endif

1429

1430 if (obc%u_W_OBCs_on_PE) then

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

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

1433 zi_dir = 1

1434 endif

1435 endif

1436 endif

1437

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

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

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

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

1442 ! topography.

1443

1444 z_i(nz+1) = 0.

1445

1446 if (.not. cs%harmonic_visc) then

1447 zh = 0.

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

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

1450 endif

1451

1452 if (cs%use_GL90_in_SSW) then

1453 z_i_gl90(nz+1) = 0.

1454 endif

1455

1456 do k=nz,1,-1

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

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

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

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

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

1462

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

1464 if (associated(obc)) then

1465 if (obc%u_E_OBCs_on_PE) then

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

1467 h_harm = h(i,j,k)

1468 h_arith = h(i,j,k)

1469 h_delta = 0.

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

1471 dz_arith = dz(i,j,k)

1472 endif

1473 endif

1474

1475 if (obc%u_W_OBCs_on_PE) then

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

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

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

1479 h_delta = 0.

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

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

1482 endif

1483 endif

1484 endif

1485

1486 if (cs%harmonic_visc) then

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

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

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

1490 ! Montgomery potential gradients at the bottom where nearly massless

1491 ! layers ride over the topography.

1492

1493 hvel(k) = h_harm

1494 dz_vel(k) = dz_harm(k)

1495

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

1497 z2 = z_i(k+1)

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

1499

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

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

1502 endif

1503

1504 z_i(k) = z_i(k+1) + dz_harm(k) * i_hbbl

1505 else

1506 zcol = zcol + dz(i,j,k)

1507 zcol_p1 = zcol_p1 + dz(i+1,j,k)

1508

1509 zh = zh + dz_harm(k)

1510

1511 z_clear = max(zcol, zcol_p1) + dmin

1512 if (zi_dir < 0) z_clear = zcol + dmin

1513 if (zi_dir > 0) z_clear = zcol_p1 + dmin

1514

1515 z_i(k) = max(zh, z_clear) * i_hbbl

1516

1517 hvel(k) = h_arith

1518 dz_vel(k) = dz_arith

1519

1520 if (u(i,j,k) * h_delta > 0.) then

1521 if (zh * i_hbbl < cs%harm_BL_val) then

1522 hvel(k) = h_harm

1523 dz_vel(k) = dz_harm(k)

1524 else

1525 z2_wt = 1.

1526 if (zh * i_hbbl < 2. * cs%harm_BL_val) &

1527 z2_wt = max(0., min(1., zh * i_hbbl * i_valbl - 1.))

1528

1529 z2 = z2_wt * (max(zh, z_clear) * i_hbbl)

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

1531

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

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

1534 endif

1535 endif

1536 endif

1537

1538 if (cs%use_GL90_in_SSW) then

1539 ! The following block calculates the normalized height above the GL90 BBL

1540 ! (z_i_gl90), using a harmonic mean between layer thicknesses. For the

1541 ! GL90 BBL we use simply a constant (Hbbl_gl90). The purpose isthat the

1542 ! GL90 coupling coefficient is zeroed out within Hbbl_gl90, to ensure

1543 ! that no momentum gets fluxed into vanished layers. The scheme is not

1544 ! sensitive to the exact value of Hbbl_gl90, as long as it is in a

1545 ! reasonable range (~1-20 m): large enough to capture vanished layers

1546 ! over topography, small enough to not contaminate the interior.

1547

1548 z_i_gl90(k) = z_i_gl90(k+1) + dz_harm(k) * i_hbbl_gl90

1549 endif

1550 enddo

1551

1552 call find_coupling_coef(a_cpl, dz_vel, i, j, dz_harm, bbl_thick, kv_bbl, z_i, &

1553 h_ml, dt, g, gv, us, cs, visc, ustar_2d, tv, work_on_u=.true., obc=obc)

1554

1555 if (allocated(hml_u)) hml_u(i,j) = h_ml

1556

1557 if (cs%use_GL90_in_SSW) then

1558 call find_coupling_coef_gl90(a_cpl_gl90, dz_vel, i, j, z_i_gl90, g, gv, &

1559 cs, varmix, work_on_u=.true.)

1560 endif

1561

1562 do_any_shelf = .false.

1563 if (associated(forces%frac_shelf_u)) then

1564 cs%a1_shelf_u(i,j) = 0.

1565 do_any_shelf = forces%frac_shelf_u(i,j) > 0.

1566

1567 if (do_any_shelf) then

1568 if (.not. cs%harmonic_visc) then

1569 zh = 0.

1570 ztop_min = min(zcol, zcol_p1)

1571 i_htbl = 1. / (visc%tbl_thick_shelf_u(i,j) + dz_neglect)

1572 endif

1573

1574 do k=1,nz

1575 if (cs%harmonic_visc) then

1576 hvel_shelf(k) = hvel(k)

1577 dz_vel_shelf(k) = dz_vel(k)

1578 else

1579 ! Find upwind-biased thickness near the surface.

1580 ! (Perhaps this needs to be done more carefully, via find_eta.)

1581

1582 h_harm = 2. * h(i,j,k) * h(i+1,j,k) &

1583 / (h(i,j,k) + h(i+1,j,k) + h_neglect)

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

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

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

1587

1588 if (associated(obc)) then

1589 if (obc%u_E_OBCs_on_PE) then

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

1591 h_harm = h(i,j,k)

1592 h_arith = h(i,j,k)

1593 h_delta = 0.

1594 dz_arith = dz(i,j,k)

1595 endif

1596 endif

1597

1598 if (obc%u_W_OBCs_on_PE) then

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

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

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

1602 h_delta = 0.

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

1604 endif

1605 endif

1606 endif

1607

1608 zcol = zcol - dz(i,j,k)

1609 zcol_p1 = zcol_p1 - dz(i+1,j,k)

1610

1611 zh = zh + dz_harm(k)

1612

1613 hvel_shelf(k) = hvel(k)

1614 dz_vel_shelf(k) = dz_vel(k)

1615

1616 if (u(i,j,k) * h_delta > 0.) then

1617 if (zh * i_htbl < cs%harm_BL_val) then

1618 hvel_shelf(k) = min(hvel(k), h_harm)

1619 dz_vel_shelf(k) = min(dz_vel(k), dz_harm(k))

1620 else

1621 z2_wt = 1.

1622 if (zh * i_htbl < 2. * cs%harm_BL_val) then

1623 z2_wt = max(0., min(1., zh * i_htbl * i_valbl - 1.))

1624 endif

1625

1626 z2 = z2_wt * (max(zh, ztop_min - min(zcol, zcol_p1)) * i_htbl)

1627 ! TODO: replace **6 with multiply

1628 topfn = 1. / (1. + 0.09 * z2**6)

1629

1630 hvel_shelf(k) = min(hvel(k), (1. - topfn) * h_arith + topfn * h_harm)

1631 dz_vel_shelf(k) = min(dz_vel(k), (1. - topfn) * dz_arith + topfn * dz_harm(k))

1632 endif

1633 endif

1634 endif

1635 enddo

1636

1637 call find_coupling_coef(a_shelf, dz_vel_shelf, i, j, dz_harm, &

1638 bbl_thick, kv_bbl, z_i, h_ml, dt, g, gv, us, cs, visc, ustar_2d, &

1639 tv, work_on_u=.true., obc=obc, shelf=.true.)

1640

1641 cs%a1_shelf_u(i,j) = a_shelf(1)

1642 endif

1643 endif

1644

1645 if (do_any_shelf) then

1646 if (cs%use_GL90_in_SSW) then

1647 do k=1,nz+1

1648 cs%a_u(i,j,k) = min(a_cpl_max, (forces%frac_shelf_u(i,j) * a_shelf(k) + &

1649 (1. - forces%frac_shelf_u(i,j)) * a_cpl(k)) + a_cpl_gl90(k))

1650

1651 ! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH

1652 ! CS%a_u(I,j,K) = min(a_cpl_max, forces%frac_shelf_u(I,j) * max(a_shelf(K), a_cpl(K)) + &

1653 ! (1. - forces%frac_shelf_u(I,j)) * a_cpl(K))

1654

1655 cs%a_u_gl90(i,j,k) = min(a_cpl_max, a_cpl_gl90(k))

1656 enddo

1657 else

1658 do k=1,nz+1

1659 cs%a_u(i,j,k) = min(a_cpl_max, (forces%frac_shelf_u(i,j) * a_shelf(k) + &

1660 (1. - forces%frac_shelf_u(i,j)) * a_cpl(k)))

1661

1662 ! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH

1663 ! CS%a_u(I,j,K) = min(a_cpl_max, forces%frac_shelf_u(I,j) * max(a_shelf(K), a_cpl(K)) + &

1664 ! (1. - forces%frac_shelf_u(I,j)) * a_cpl(K))

1665 enddo

1666 endif

1667

1668 do k=1,nz

1669 ! Should we instead take the inverse of the average of the inverses?

1670 cs%h_u(i,j,k) = forces%frac_shelf_u(i,j) * hvel_shelf(k) &

1671 + (1. - forces%frac_shelf_u(i,j)) * hvel(k) + h_neglect

1672 enddo

1673 else

1674 if (cs%use_GL90_in_SSW) then

1675 do k=1,nz+1

1676 a_cpl(k) = a_cpl(k) + a_cpl_gl90(k)

1677 enddo

1678

1679 do k=1,nz+1

1680 cs%a_u_gl90(i,j,k) = min(a_cpl_max, a_cpl_gl90(k))

1681 enddo

1682 endif

1683

1684 do k=1,nz+1

1685 cs%a_u(i,j,k) = min(a_cpl_max, a_cpl(k))

1686 enddo

1687

1688 do k=1,nz

1689 cs%h_u(i,j,k) = hvel(k) + h_neglect

1690 enddo

1691 endif

1692

1693 ! Diagnose total Kv at u-points

1694 if (cs%id_Kv_u > 0) then

1695 do k=1,nz

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

1697 enddo

1698 endif

1699

1700 ! Diagnose GL90 Kv at u-points

1701 if (cs%id_Kv_gl90_u > 0) then

1702 do k=1,nz

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

1704 enddo

1705 endif

1706 endif ; enddo ; enddo

1707

1708 ! Now work on v-points.

1709

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

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

1712 if (cs%use_GL90_in_SSW) then

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

1714 endif

1715

1716 if (cs%bottomdraglaw) then

1717 kv_bbl = visc%Kv_bbl_v(i,j)

1718 bbl_thick = visc%bbl_thick_v(i,j) + dz_neglect

1719 i_hbbl = 1. / bbl_thick

1720 endif

1721

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

1723 zi_dir = 0

1724

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

1726 if (associated(obc)) then

1727 if (obc%v_N_OBCs_on_PE) then

1728 if (obc%segnum_v(i,j) > 0) then

1729 dmin = g%bathyT(i,j)

1730 zi_dir = -1

1731 endif

1732 endif

1733

1734 if (obc%v_S_OBCs_on_PE) then

1735 if (obc%segnum_v(i,j) < 0) then

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

1737 zi_dir = 1

1738 endif

1739 endif

1740 endif

1741

1742 z_i(nz+1) = 0.

1743

1744 if (.not. cs%harmonic_visc) then

1745 zh = 0.

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

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

1748 endif

1749

1750 if (cs%use_GL90_in_SSW) then

1751 z_i_gl90(nz+1) = 0.

1752 endif

1753

1754 do k=nz,1,-1

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

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

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

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

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

1760

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

1762 if (associated(obc)) then

1763 if (obc%v_N_OBCs_on_PE) then

1764 if (obc%segnum_v(i,j) > 0) then

1765 h_harm = h(i,j,k)

1766 h_arith = h(i,j,k)

1767 h_delta = 0.

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

1769 dz_arith = dz(i,j,k)

1770 endif

1771 endif

1772

1773 if (obc%v_S_OBCs_on_PE) then

1774 if (obc%segnum_v(i,j) < 0) then

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

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

1777 h_delta = 0.

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

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

1780 endif

1781 endif

1782 endif

1783

1784 if (cs%harmonic_visc) then

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

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

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

1788 ! Montgomery potential gradients at the bottom where nearly massless

1789 ! layers ride over the topography.

1790

1791 hvel(k) = h_harm

1792 dz_vel(k) = dz_harm(k)

1793

1794 if (v(i,j,k) * h_delta < 0) then

1795 z2 = z_i(k+1)

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

1797

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

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

1800 endif

1801

1802 z_i(k) = z_i(k+1) + dz_harm(k) * i_hbbl

1803 else

1804 zcol = zcol + dz(i,j,k)

1805 zcol_p1 = zcol_p1 + dz(i,j+1,k)

1806

1807 zh = zh + dz_harm(k)

1808

1809 z_clear = max(zcol, zcol_p1) + dmin

1810 if (zi_dir < 0) z_clear = zcol + dmin

1811 if (zi_dir > 0) z_clear = zcol_p1 + dmin

1812

1813 z_i(k) = max(zh, z_clear) * i_hbbl

1814

1815 hvel(k) = h_arith

1816 dz_vel(k) = dz_arith

1817

1818 if (v(i,j,k) * h_delta > 0) then

1819 if (zh * i_hbbl < cs%harm_BL_val) then

1820 hvel(k) = h_harm

1821 dz_vel(k) = dz_harm(k)

1822 else

1823 z2_wt = 1.

1824 if (zh * i_hbbl < 2. * cs%harm_BL_val) &

1825 z2_wt = max(0., min(1., zh * i_hbbl * i_valbl - 1.))

1826

1827 ! TODO: should z_clear be used here?

1828 z2 = z2_wt * (max(zh, max(zcol, zcol_p1) + dmin) * i_hbbl)

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

1830

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

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

1833 endif

1834 endif

1835 endif

1836

1837 if (cs%use_GL90_in_SSW) then

1838 ! The following block calculates the normalized height above the GL90 BBL

1839 ! (z_i_gl90), using a harmonic mean between layer thicknesses. For the

1840 ! GL90 BBL we use simply a constant (Hbbl_gl90). The purpose is that the

1841 ! GL90 coupling coefficient is zeroed out within Hbbl_gl90, to ensure

1842 ! that no momentum gets fluxed into vanished layers. The scheme is not

1843 ! sensitive to the exact value of Hbbl_gl90, as long as it is in a

1844 ! reasonable range (~1-20 m): large enough to capture vanished layers

1845 ! over topography, small enough to not contaminate the interior.

1846

1847 z_i_gl90(k) = z_i_gl90(k+1) + dz_harm(k) * i_hbbl_gl90

1848 endif

1849 enddo

1850

1851 call find_coupling_coef(a_cpl, dz_vel, i, j, dz_harm, bbl_thick, kv_bbl, z_i, &

1852 h_ml, dt, g, gv, us, cs, visc, ustar_2d, tv, work_on_u=.false., obc=obc)

1853

1854 if (allocated(hml_v)) hml_v(i,j) = h_ml

1855

1856 if (cs%use_GL90_in_SSW) then

1857 call find_coupling_coef_gl90(a_cpl_gl90, dz_vel, i, j, z_i_gl90, g, gv, &

1858 cs, varmix, work_on_u=.false.)

1859 endif

1860

1861 do_any_shelf = .false.

1862 if (associated(forces%frac_shelf_v)) then

1863 cs%a1_shelf_v(i,j) = 0.

1864 do_any_shelf = forces%frac_shelf_v(i,j) > 0.

1865

1866 if (do_any_shelf) then

1867 ! Initialize non-harmonic depths

1868 if (.not. cs%harmonic_visc) then

1869 zh = 0.

1870 ztop_min = min(zcol, zcol_p1)

1871 i_htbl = 1. / (visc%tbl_thick_shelf_v(i,j) + dz_neglect)

1872 endif

1873

1874 do k=1,nz

1875 if (cs%harmonic_visc) then

1876 hvel_shelf(k) = hvel(k)

1877 dz_vel_shelf(k) = dz_vel(k)

1878 else

1879 ! Find upwind-biased thickness near the surface.

1880 ! Perhaps this needs to be done more carefully, via find_eta.

1881 h_harm = 2. * h(i,j,k) * h(i,j+1,k) &

1882 / (h(i,j,k) + h(i,j+1,k) + h_neglect)

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

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

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

1886

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

1888 if (associated(obc)) then

1889 if (obc%v_N_OBCs_on_PE) then

1890 if (obc%segnum_v(i,j) > 0) then

1891 h_harm = h(i,j,k)

1892 h_arith = h(i,j,k)

1893 h_delta = 0.

1894 dz_arith = dz(i,j,k)

1895 endif

1896 endif

1897

1898 if (obc%v_S_OBCs_on_PE) then

1899 if (obc%segnum_v(i,j) < 0) then

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

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

1902 h_delta = 0.

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

1904 endif

1905 endif

1906 endif

1907

1908 zcol = zcol - dz(i,j,k)

1909 zcol_p1 = zcol_p1 - dz(i,j+1,k)

1910

1911 zh = zh + dz_harm(k)

1912

1913 hvel_shelf(k) = hvel(k)

1914 dz_vel_shelf(k) = dz_vel(k)

1915

1916 if (v(i,j,k) * h_delta > 0.) then

1917 if (zh * i_htbl < cs%harm_BL_val) then

1918 hvel_shelf(k) = min(hvel(k), h_harm)

1919 dz_vel_shelf(k) = min(dz_vel(k), dz_harm(k))

1920 else

1921 z2_wt = 1.

1922 if (zh * i_htbl < 2. * cs%harm_BL_val) &

1923 z2_wt = max(0., min(1., zh * i_htbl * i_valbl - 1.))

1924

1925 z2 = z2_wt * (max(zh, ztop_min - min(zcol, zcol_p1)) * i_htbl)

1926 ! TODO: Replace **6

1927 topfn = 1. / (1. + 0.09 * z2**6)

1928

1929 hvel_shelf(k) = min(hvel(k), (1. - topfn) * h_arith + topfn * h_harm)

1930 dz_vel_shelf(k) = min(dz_vel(k), (1. - topfn) * dz_arith + topfn * dz_harm(k))

1931 endif

1932 endif

1933 endif

1934 enddo

1935

1936 call find_coupling_coef(a_shelf, dz_vel_shelf, i, j, dz_harm, &

1937 bbl_thick, kv_bbl, z_i, h_ml, dt, g, gv, us, cs, visc, ustar_2d, &

1938 tv, work_on_u=.false., obc=obc, shelf=.true.)

1939

1940 cs%a1_shelf_v(i,j) = a_shelf(1)

1941 endif

1942 endif

1943

1944 if (do_any_shelf) then

1945 if (cs%use_GL90_in_SSW) then

1946 do k=1,nz+1

1947 cs%a_v(i,j,k) = min(a_cpl_max, (forces%frac_shelf_v(i,j) * a_shelf(k) + &

1948 (1. - forces%frac_shelf_v(i,j)) * a_cpl(k)) + a_cpl_gl90(k))

1949

1950 ! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH

1951 ! CS%a_v(I,j,K) = min(a_cpl_max, forces%frac_shelf_v(I,j) * max(a_shelf(K), a_cpl(K)) + &

1952 ! (1. - forces%frac_shelf_v(I,j)) * a_cpl(K))

1953

1954 cs%a_v_gl90(i,j,k) = min(a_cpl_max, a_cpl_gl90(k))

1955 enddo

1956 else

1957 do k=1,nz+1

1958 cs%a_v(i,j,k) = min(a_cpl_max, (forces%frac_shelf_v(i,j) * a_shelf(k) + &

1959 (1. - forces%frac_shelf_v(i,j)) * a_cpl(k)))

1960 ! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH

1961 ! CS%a_v(I,j,K) = min(a_cpl_max, forces%frac_shelf_v(I,j) * max(a_shelf(K), a_cpl(K)) + &

1962 ! (1. - forces%frac_shelf_v(I,j)) * a_cpl(K))

1963 enddo

1964 endif

1965

1966 do k=1,nz

1967 ! Should we instead take the inverse of the average of the inverses?

1968 cs%h_v(i,j,k) = forces%frac_shelf_v(i,j) * hvel_shelf(k) &

1969 + (1. - forces%frac_shelf_v(i,j)) * hvel(k) + h_neglect

1970 enddo

1971 else

1972 if (cs%use_GL90_in_SSW) then

1973 do k=1,nz+1

1974 a_cpl(k) = a_cpl(k) + a_cpl_gl90(k)

1975 enddo

1976

1977 do k=1,nz+1

1978 cs%a_v_gl90(i,j,k) = min(a_cpl_max, a_cpl_gl90(k))

1979 enddo

1980 endif

1981

1982 do k=1,nz+1

1983 cs%a_v(i,j,k) = min(a_cpl_max, a_cpl(k))

1984 enddo

1985

1986 do k=1,nz

1987 cs%h_v(i,j,k) = hvel(k) + h_neglect

1988 enddo

1989 endif

1990

1991 ! Diagnose total Kv at v-points

1992 if (cs%id_Kv_v > 0) then

1993 do k=1,nz

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

1995 enddo

1996 endif

1997

1998 ! Diagnose GL90 Kv at v-points

1999 if (cs%id_Kv_gl90_v > 0) then

2000 do k=1,nz

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

2002 enddo

2003 endif

2004 endif ; enddo ; enddo

2005

2006 if (cs%debug) then

2007 call uvchksum("vertvisc_coef h_[uv]", cs%h_u, cs%h_v, g%HI, haloshift=0, &

2008 unscale=gv%H_to_m, scalar_pair=.true.)

2009 call uvchksum("vertvisc_coef a_[uv]", cs%a_u, cs%a_v, g%HI, haloshift=0, &

2010 unscale=gv%H_to_m*us%s_to_T, scalar_pair=.true.)

2011 if (allocated(hml_u) .and. allocated(hml_v)) &

2012 call uvchksum("vertvisc_coef hML_[uv]", hml_u, hml_v, g%HI, &

2013 haloshift=0, unscale=us%Z_to_m, scalar_pair=.true.)

2014 endif

2015

2016! Offer diagnostic fields for averaging.

2017 if (query_averaging_enabled(cs%diag)) then

2018 if (associated(visc%Kv_slow) .and. (cs%id_Kv_slow > 0)) &

2019 call post_data(cs%id_Kv_slow, visc%Kv_slow, cs%diag)

2020 if (cs%id_Kv_u > 0) call post_data(cs%id_Kv_u, kv_u, cs%diag)

2021 if (cs%id_Kv_v > 0) call post_data(cs%id_Kv_v, kv_v, cs%diag)

2022 if (cs%id_Kv_gl90_u > 0) call post_data(cs%id_Kv_gl90_u, kv_gl90_u, cs%diag)

2023 if (cs%id_Kv_gl90_v > 0) call post_data(cs%id_Kv_gl90_v, kv_gl90_v, cs%diag)

2024 if (cs%id_au_vv > 0) call post_data(cs%id_au_vv, cs%a_u, cs%diag)

2025 if (cs%id_av_vv > 0) call post_data(cs%id_av_vv, cs%a_v, cs%diag)

2026 if (cs%id_au_gl90_vv > 0) call post_data(cs%id_au_gl90_vv, cs%a_u_gl90, cs%diag)

2027 if (cs%id_av_gl90_vv > 0) call post_data(cs%id_av_gl90_vv, cs%a_v_gl90, cs%diag)

2028 if (cs%id_h_u > 0) call post_data(cs%id_h_u, cs%h_u, cs%diag)

2029 if (cs%id_h_v > 0) call post_data(cs%id_h_v, cs%h_v, cs%diag)

2030 if (cs%id_hML_u > 0) call post_data(cs%id_hML_u, hml_u, cs%diag)

2031 if (cs%id_hML_v > 0) call post_data(cs%id_hML_v, hml_v, cs%diag)

2032 endif

2033

2034 if (allocated(hml_u)) deallocate(hml_u)

2035 if (allocated(hml_v)) deallocate(hml_v)

2036

2037end subroutine vertvisc_coef

2038

2039

2040!> Calculate the 'coupling coefficient' (a_cpl) at the interfaces.

2041!! If BOTTOMDRAGLAW is defined, the minimum of Hbbl and half the adjacent

2042!! layer thicknesses are used to calculate a_cpl near the bottom.

2043subroutine find_coupling_coef(a_cpl, hvel, i, j, h_harm, bbl_thick, kv_bbl, z_i, h_ml, &

2044 dt, G, GV, US, CS, visc, Ustar_2d, tv, work_on_u, OBC, shelf)

2045 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure

2046 type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid structure

2047 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type

2048 real, dimension(SZK_(GV)+1), &

2049 intent(out) :: a_cpl !< Coupling coefficient across interfaces [H T-1 ~> m s-1 or Pa s m-1]

2050 real, dimension(SZK_(GV)), &

2051 intent(in) :: hvel !< Distance between interfaces at velocity points [Z ~> m]

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

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

2054 real, dimension(SZK_(GV)), &

2055 intent(in) :: h_harm !< Harmonic mean of thicknesses around a velocity

2056 !! grid point [Z ~> m]

2057 real, intent(in) :: bbl_thick !< Bottom boundary layer thickness [Z ~> m]

2058 real, intent(in) :: kv_bbl !< Bottom boundary layer viscosity, exclusive of

2059 !! any depth-dependent contributions from

2060 !! visc%Kv_shear [H Z T-1 ~> m2 s-1 or Pa s]

2061 real, dimension(SZK_(GV)+1), &

2062 intent(in) :: z_i !< Estimate of interface heights above the bottom,

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

2064 real, intent(out) :: h_ml !< Mixed layer depth [Z ~> m]

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

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

2067 type(vertvisc_type), intent(in) :: visc !< Structure containing viscosities and bottom drag

2068 real, dimension(SZI_(G),SZJ_(G)), &

2069 intent(in) :: Ustar_2d !< The wind friction velocity, calculated using

2070 !! the Boussinesq reference density or the

2071 !! time-evolving surface density in non-Boussinesq

2072 !! mode [Z T-1 ~> m s-1]

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

2074 !! thermodynamic fields.

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

2076 !! otherwise they are v-points

2077 type(ocean_obc_type), pointer :: OBC !< Open boundary condition structure

2078 logical, optional, intent(in) :: shelf !< If present and true, use a surface boundary

2079 !! condition appropriate for an ice shelf.

2080

2081 ! Local variables

2082

2083 real :: &

2084 u_star, & ! ustar at a velocity point [Z T-1 ~> m s-1]

2085 tau_mag, & ! The magnitude of the wind stress at a velocity point including gustiness [H Z T-2 ~> m2 s-2 or Pa]

2086 absf, & ! The average of the neighboring absolute values of f [T-1 ~> s-1].

2087 rho_av1, & ! The harmonic mean surface layer density at velocity points [R ~> kg m-3]

2088 z_t, & ! The distance from the top, sometimes normalized

2089 ! by Hmix, [Z ~> m] or [nondim].

2090 kv_tbl, & ! The viscosity in a top boundary layer under ice [H Z T-1 ~> m2 s-1 or Pa s]

2091 tbl_thick, &! The thickness of the top boundary layer [Z ~> m]

2092 kv_add, & ! A viscosity to add [H Z T-1 ~> m2 s-1 or Pa s]

2093 kv_tot ! The total viscosity at an interface [H Z T-1 ~> m2 s-1 or Pa s]

2094 integer :: &

2095 nk_in_ml ! The index of the deepest interface in the mixed layer.

2096 real :: h_shear ! The distance over which shears occur [Z ~> m].

2097 real :: dhc ! The distance between the center of adjacent layers [Z ~> m].

2098 real :: visc_ml ! The mixed layer viscosity [H Z T-1 ~> m2 s-1 or Pa s].

2099 real :: I_Hmix ! The inverse of the mixed layer thickness [Z-1 ~> m-1].

2100 real :: a_ml ! The layer coupling coefficient across an interface in

2101 ! the mixed layer [H T-1 ~> m s-1 or Pa s m-1].

2102 real :: a_floor ! A lower bound on the layer coupling coefficient across an interface in

2103 ! the mixed layer [H T-1 ~> m s-1 or Pa s m-1].

2104 real :: I_amax ! The inverse of the maximum coupling coefficient [T H-1 ~> s m-1 or s m2 kg-1].

2105 real :: temp1 ! A temporary variable [Z2 ~> m2]

2106 real :: ustar2_denom ! A temporary variable in the surface boundary layer turbulence

2107 ! calculations [H Z-1 T-1 ~> s-1 or kg m-3 s-1]

2108 real :: h_neglect ! A vertical distance that is so small it is usually lost

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

2110 real :: z2 ! A copy of z_i [nondim]

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

2112 real :: topfn ! A function that is 1 at the top and small far from it [nondim]

2113 real :: kv_top ! A viscosity associated with the top boundary layer [H Z T-1 ~> m2 s-1 or Pa s]

2114 logical :: do_shelf, do_OBCs, can_exit

2115 integer :: k

2116 integer :: nz, max_nk

2117

2118 nz = gv%ke

2119

2120 h_neglect = gv%dZ_subroundoff

2121

2122 if (cs%answer_date < 20190101) then

2123 ! The maximum coupling coefficient was originally introduced to avoid

2124 ! truncation error problems in the tridiagonal solver. Effectively, the 1e-10

2125 ! sets the maximum coupling coefficient increment to 1e10 m per timestep.

2126 i_amax = (1.0e-10*gv%H_to_m) * dt

2127 else

2128 i_amax = 0.0

2129 endif

2130

2131 do_shelf = .false. ; if (present(shelf)) do_shelf = shelf

2132

2133 do_obcs = .false.

2134 if (associated(obc)) then

2135 if (work_on_u) then

2136 do_obcs = obc%u_E_OBCs_on_PE .or. obc%u_W_OBCs_on_PE

2137 else

2138 do_obcs = obc%v_N_OBCs_on_PE .or. obc%v_S_OBCs_on_PE

2139 endif

2140 endif

2141

2142 a_cpl(:) = 0.

2143 h_ml = 0.

2144

2145 if (cs%Kvml_invZ2 > 0. .and. .not. do_shelf) then

2146 i_hmix = 1. / (cs%Hmix + h_neglect)

2147 z_t = h_neglect * i_hmix

2148 endif

2149

2150 do k=2,nz

2151 kv_tot = cs%Kv

2152

2153 if (cs%Kvml_invZ2 > 0. .and. .not. do_shelf) then

2154 ! This is an older (vintage ~1997) way to prevent wind stresses from driving very

2155 ! large flows in nearly massless near-surface layers when there is not a physically-

2156 ! based surface boundary layer parameterization. It does not have a plausible

2157 ! physical basis, and probably should not be used.

2158 z_t = z_t + h_harm(k-1) * i_hmix

2159 kv_tot = cs%Kv + cs%Kvml_invZ2 / ((z_t * z_t) * &

2160 (1. + 0.09 * z_t * z_t * z_t * z_t * z_t * z_t))

2161 endif

2162

2163 if (associated(visc%Kv_shear)) then

2164 ! Add in viscosities that are determined by physical processes that are handled in

2165 ! other modules, and which do not respond immediately to the changing layer thicknesses.

2166 ! These processes may include shear-driven mixing or contributions from some boundary

2167 ! layer turbulence schemes. Other viscosity contributions that respond to the evolving

2168 ! layer thicknesses or the surface wind stresses are added later.

2169 if (work_on_u) then

2170 kv_add = 0.5 * (visc%Kv_shear(i,j,k) + visc%Kv_shear(i+1,j,k))

2171

2172 if (do_obcs) then

2173 if (obc%u_E_OBCs_on_PE) then

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

2175 kv_add = visc%Kv_shear(i,j,k)

2176 endif

2177 endif

2178

2179 if (obc%u_W_OBCs_on_PE) then

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

2181 kv_add = visc%Kv_shear(i+1,j,k)

2182 endif

2183 endif

2184 endif

2185

2186 kv_tot = kv_tot + kv_add

2187 else

2188 kv_add = 0.5 * (visc%Kv_shear(i,j,k) + visc%Kv_shear(i,j+1,k))

2189

2190 if (do_obcs) then

2191 if (obc%v_N_OBCs_on_PE) then

2192 if (obc%segnum_v(i,j) > 0) then

2193 kv_add = visc%Kv_shear(i,j,k)

2194 endif

2195 endif

2196

2197 if (obc%v_S_OBCs_on_PE) then

2198 if (obc%segnum_v(i,j) < 0) then

2199 kv_add = visc%Kv_shear(i,j+1,k)

2200 endif

2201 endif

2202 endif

2203

2204 kv_tot = kv_tot + kv_add

2205 endif

2206 endif

2207

2208 if (associated(visc%Kv_shear_Bu)) then

2209 ! This is similar to what was done above, but for contributions coming from the corner

2210 ! (vorticity) points. Because OBCs run through the faces and corners there is no need

2211 ! to further modify these viscosities here to take OBCs into account.

2212 if (work_on_u) then

2213 kv_tot = kv_tot + 0.5 * (visc%Kv_shear_Bu(i,j-1,k) + visc%Kv_shear_Bu(i,j,k))

2214 else

2215 kv_tot = kv_tot + 0.5 * (visc%Kv_shear_Bu(i-1,j,k) + visc%Kv_shear_Bu(i,j,k))

2216 endif

2217 endif

2218

2219 ! Set the viscous coupling coefficients, excluding surface mixed layer contributions

2220 ! for now, but including viscous bottom drag, working up from the bottom.

2221 if (cs%bottomdraglaw) then

2222 ! botfn determines when a point is within the influence of the bottom

2223 ! boundary layer, going from 1 at the bottom to 0 in the interior.

2224 z2 = z_i(k)

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

2226

2227 kv_tot = kv_tot + (kv_bbl - cs%Kv) * botfn

2228 dhc = 0.5 * (hvel(k) + hvel(k-1))

2229 if (dhc > bbl_thick) then

2230 h_shear = ((1. - botfn) * dhc + botfn * bbl_thick) + h_neglect

2231 else

2232 h_shear = dhc + h_neglect

2233 endif

2234

2235 ! Calculate the coupling coefficients from the viscosities.

2236 a_cpl(k) = kv_tot / (h_shear + (i_amax * kv_tot))

2237 elseif (abs(cs%Kv_extra_bbl) > 0.0) then

2238 ! There is a simple enhancement of the near-bottom viscosities, but no

2239 ! adjustment of the viscous coupling length scales to give a particular

2240 ! bottom stress.

2241

2242 ! botfn determines when a point is within the influence of the bottom

2243 ! boundary layer, going from 1 at the bottom to 0 in the interior.

2244 z2 = z_i(k)

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

2246

2247 kv_tot = kv_tot + cs%Kv_extra_bbl * botfn

2248 h_shear = 0.5 * (hvel(k) + hvel(k-1) + h_neglect)

2249

2250 ! Calculate the coupling coefficients from the viscosities.

2251 a_cpl(k) = kv_tot / (h_shear + i_amax * kv_tot)

2252 else

2253 ! Any near-bottom viscous enhancements were already incorporated into

2254 ! Kv_tot, and there is no adjustment of the viscous coupling length

2255 ! scales to give a particular bottom stress.

2256

2257 h_shear = 0.5 * (hvel(k) + hvel(k-1) + h_neglect)

2258 ! Calculate the coupling coefficients from the viscosities.

2259 a_cpl(k) = kv_tot / (h_shear + i_amax * kv_tot)

2260 endif

2261 enddo

2262

2263 ! Assign the bottom coupling coefficients

2264 if (cs%bottomdraglaw) then

2265 dhc = hvel(nz) * 0.5

2266 a_cpl(nz+1) = kv_bbl / ((min(dhc, bbl_thick) + h_neglect) + i_amax * kv_bbl)

2267 elseif (abs(cs%Kv_extra_bbl) > 0.0) then

2268 a_cpl(nz+1) = (cs%Kv + cs%Kv_extra_bbl) &

2269 / ((0.5 * hvel(nz) + h_neglect) + i_amax * (cs%Kv + cs%Kv_extra_bbl))

2270 else

2271 a_cpl(nz+1) = cs%Kv / ((0.5 * hvel(nz) + h_neglect) + i_amax * cs%Kv)

2272 endif

2273

2274 ! Add surface intensified viscous coupling, either as a no-slip boundary condition under a

2275 ! rigid ice-shelf, or due to wind-stress driven surface boundary layer mixing that has not

2276 ! already been added via visc%Kv_shear.

2277 if (do_shelf) then

2278 ! Set the coefficients to include the no-slip surface stress.

2279 if (work_on_u) then

2280 kv_tbl = visc%Kv_tbl_shelf_u(i,j)

2281 tbl_thick = visc%tbl_thick_shelf_u(i,j) + h_neglect

2282 else

2283 kv_tbl = visc%Kv_tbl_shelf_v(i,j)

2284 tbl_thick = visc%tbl_thick_shelf_v(i,j) + h_neglect

2285 endif

2286

2287 z_t = 0.0

2288

2289 ! If a_cpl(1) were not already 0, it would be added here.

2290 if (0.5 * hvel(1) > tbl_thick) then

2291 a_cpl(1) = kv_tbl / (tbl_thick + i_amax * kv_tbl)

2292 else

2293 a_cpl(1) = kv_tbl / ((0.5 * hvel(1) + h_neglect) + i_amax * kv_tbl)

2294 endif

2295

2296 do k=2,nz

2297 z_t = z_t + hvel(k-1) / tbl_thick

2298 topfn = 1. / (1. + 0.09 * z_t**6)

2299

2300 dhc = 0.5 * (hvel(k) + hvel(k-1))

2301 if (dhc > tbl_thick) then

2302 h_shear = ((1. - topfn) * dhc + topfn * tbl_thick) + h_neglect

2303 else

2304 h_shear = dhc + h_neglect

2305 endif

2306

2307 kv_top = topfn * kv_tbl

2308 a_cpl(k) = a_cpl(k) + kv_top / (h_shear + i_amax * kv_top)

2309 enddo

2310 elseif (cs%dynamic_viscous_ML .or. (gv%nkml>0) .or. cs%fixed_LOTW_ML .or. cs%apply_LOTW_floor) then

2311

2312 ! Find the friction velocity and the absolute value of the Coriolis parameter at this point.

2313 u_star = 0. ! Zero out the friction velocity on land points.

2314 tau_mag = 0. ! Zero out the friction velocity on land points.

2315

2316 if (allocated(tv%SpV_avg)) then

2317 rho_av1 = 0.

2318

2319 if (work_on_u) then

2320 u_star = 0.5 * (ustar_2d(i,j) + ustar_2d(i+1,j))

2321 rho_av1 = 2. / (tv%SpV_avg(i,j,1) + tv%SpV_avg(i+1,j,1))

2322 absf = 0.5 * (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i,j)))

2323

2324 if (do_obcs) then

2325 if (obc%u_E_OBCs_on_PE) then

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

2327 u_star = ustar_2d(i,j)

2328 rho_av1 = 1. / tv%SpV_avg(i,j,1)

2329 endif

2330 endif

2331

2332 if (obc%u_W_OBCs_on_PE) then

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

2334 u_star = ustar_2d(i+1,j)

2335 rho_av1 = 1. / tv%SpV_avg(i+1,j,1)

2336 endif

2337 endif

2338 endif

2339 else

2340 u_star = 0.5 * (ustar_2d(i,j) + ustar_2d(i,j+1))

2341 rho_av1 = 2. / (tv%SpV_avg(i,j,1) + tv%SpV_avg(i,j+1,1))

2342 absf = 0.5 * (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))

2343

2344 if (do_obcs) then

2345 if (obc%v_N_OBCs_on_PE) then

2346 if (obc%segnum_v(i,j) > 0) then

2347 u_star = ustar_2d(i,j)

2348 rho_av1 = 1. / tv%SpV_avg(i,j,1)

2349 endif

2350 endif

2351

2352 if (obc%v_S_OBCs_on_PE) then

2353 if (obc%segnum_v(i,j) < 0) then

2354 u_star = ustar_2d(i,j+1)

2355 rho_av1 = 1. / tv%SpV_avg(i,j+1,1)

2356 endif

2357 endif

2358 endif

2359 endif

2360

2361 tau_mag = gv%RZ_to_H * rho_av1 * u_star**2

2362 else ! (.not.allocated(tv%SpV_avg))

2363 if (work_on_u) then

2364 u_star = 0.5 * (ustar_2d(i,j) + ustar_2d(i+1,j))

2365 absf = 0.5 * (abs(g%CoriolisBu(i,j-1)) + abs(g%CoriolisBu(i,j)))

2366

2367 if (do_obcs) then

2368 if (obc%u_E_OBCs_on_PE) then

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

2370 u_star = ustar_2d(i,j)

2371 endif

2372 endif

2373

2374 if (obc%u_W_OBCs_on_PE) then

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

2376 u_star = ustar_2d(i+1,j)

2377 endif

2378 endif

2379 endif

2380 else

2381 u_star = 0.5 * (ustar_2d(i,j) + ustar_2d(i,j+1))

2382 absf = 0.5 * (abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))

2383

2384 if (do_obcs) then

2385 if (obc%v_N_OBCs_on_PE) then

2386 if (obc%segnum_v(i,j) > 0) then

2387 u_star = ustar_2d(i,j)

2388 endif

2389 endif

2390

2391 if (obc%v_S_OBCs_on_PE) then

2392 if (obc%segnum_v(i,j) < 0) then

2393 u_star = ustar_2d(i,j+1)

2394 endif

2395 endif

2396 endif

2397 endif

2398

2399 tau_mag = gv%Z_to_H * u_star**2

2400 endif

2401

2402 ! Determine the thickness of the surface ocean boundary layer and its extent in index space.

2403 nk_in_ml = 0

2404 if (cs%dynamic_viscous_ML) then

2405 ! The fractional number of layers that are within the viscous boundary layer were

2406 ! previously stored in visc%nkml_visc_[uv].

2407 h_ml = h_neglect

2408 max_nk = 0

2409

2410 if (work_on_u) then

2411 nk_in_ml = ceiling(visc%nkml_visc_u(i,j))

2412 max_nk = max(max_nk, nk_in_ml)

2413

2414 do k=1,max_nk

2415 if (k <= visc%nkml_visc_u(i,j)) then ! This layer is all in the ML.

2416 h_ml = h_ml + hvel(k)

2417 elseif (k < visc%nkml_visc_u(i,j) + 1.) then ! Part of this layer is in the ML.

2418 h_ml = h_ml + ((visc%nkml_visc_u(i,j) + 1.) - k) * hvel(k)

2419 endif

2420 enddo

2421 else

2422 nk_in_ml = ceiling(visc%nkml_visc_v(i,j))

2423 max_nk = max(max_nk, nk_in_ml)

2424

2425 do k=1,max_nk

2426 if (k <= visc%nkml_visc_v(i,j)) then ! This layer is all in the ML.

2427 h_ml = h_ml + hvel(k)

2428 elseif (k < visc%nkml_visc_v(i,j) + 1.) then ! Part of this layer is in the ML.

2429 h_ml = h_ml + ((visc%nkml_visc_v(i,j) + 1.) - k) * hvel(k)

2430 endif

2431 enddo

2432 endif

2433 elseif (gv%nkml>0) then

2434 ! This is a simple application of a refined-bulk mixed layer with GV%nkml sublayers.

2435 max_nk = gv%nkml

2436 nk_in_ml = gv%nkml

2437

2438 h_ml = h_neglect

2439

2440 do k=1,gv%nkml

2441 h_ml = h_ml + hvel(k)

2442 enddo

2443 elseif (cs%fixed_LOTW_ML .or. cs%apply_LOTW_floor) then

2444 ! Determine which interfaces are within CS%Hmix of the surface, and set the viscous

2445 ! boundary layer thickness to the smaller of CS%Hmix and the depth of the ocean.

2446 h_ml = 0.0

2447 do k=1,nz

2448 can_exit = .true.

2449 if (h_ml < cs%Hmix) then

2450 nk_in_ml = k

2451

2452 if (h_ml + hvel(k) < cs%Hmix) then

2453 h_ml = h_ml + hvel(k)

2454 can_exit = .false. ! Part of the next deeper layer is also in the mixed layer.

2455 else

2456 h_ml = cs%Hmix

2457 endif

2458 endif

2459

2460 if (can_exit) exit ! All remaining layers in this row are below the mixed layer depth.

2461 enddo

2462

2463 max_nk = max(0, nk_in_ml)

2464 endif

2465

2466 ! Avoid working on columns where the viscous coupling could not be increased.

2467 if (u_star <= 0.) nk_in_ml = 0

2468

2469 ! Set the viscous coupling at the interfaces as the larger of what was previously

2470 ! set and the contributions from the surface boundary layer.

2471 z_t = 0.

2472 if (cs%apply_LOTW_floor .and. &

2473 (cs%dynamic_viscous_ML .or. gv%nkml > 0 .or. cs%fixed_LOTW_ML)) then

2474 do k=2,max_nk

2475 if (k <= nk_in_ml) then

2476 z_t = z_t + hvel(k-1)

2477

2478 ! The viscosity in visc_ml is set to go to 0 at the mixed layer top and bottom

2479 ! (in a log-layer) and be further limited by rotation to give the natural Ekman length.

2480 temp1 = (z_t * h_ml - z_t * z_t)

2481 if (gv%Boussinesq) then

2482 ustar2_denom = (cs%vonKar * gv%Z_to_H * u_star**2) &

2483 / (absf * temp1 + (h_ml + h_neglect) * u_star)

2484 else

2485 ustar2_denom = (cs%vonKar * tau_mag) &

2486 / (absf * temp1 + (h_ml + h_neglect) * u_star)

2487 endif

2488

2489 visc_ml = temp1 * ustar2_denom

2490 ! Set the viscous coupling based on the model's vertical resolution. The omission of

2491 ! the I_amax factor here is consistent with answer dates above 20190101.

2492 a_ml = visc_ml / (0.25 * (hvel(k) + hvel(k-1) + h_neglect))

2493

2494 ! As a floor on the viscous coupling, assume that the length scale in the denominator can

2495 ! not be larger than the distance from the surface, consistent with a logarithmic velocity

2496 ! profile. This is consistent with visc_ml, but cancels out common factors of z_t.

2497 a_floor = (h_ml - z_t) * ustar2_denom

2498

2499 ! Choose the largest estimate of a_cpl.

2500 a_cpl(k) = max(a_cpl(k), a_ml, a_floor)

2501 ! An option could be added to change this to: a_cpl(i,K) = max(a_cpl(i,K) + a_ml, a_floor)

2502 endif

2503 enddo

2504 elseif (cs%apply_LOTW_floor) then

2505 do k=2,max_nk

2506 if (k <= nk_in_ml) then

2507 z_t = z_t + hvel(k-1)

2508

2509 temp1 = (z_t * h_ml - z_t * z_t)

2510 if (gv%Boussinesq) then

2511 ustar2_denom = (cs%vonKar * gv%Z_to_H * u_star**2) &

2512 / (absf * temp1 + (h_ml + h_neglect) * u_star)

2513 else

2514 ustar2_denom = (cs%vonKar * tau_mag) &

2515 / (absf * temp1 + (h_ml + h_neglect) * u_star)

2516 endif

2517

2518 ! As a floor on the viscous coupling, assume that the length scale in the denominator can not

2519 ! be larger than the distance from the surface, consistent with a logarithmic velocity profile.

2520 a_cpl(k) = max(a_cpl(k), (h_ml - z_t) * ustar2_denom)

2521 endif

2522 enddo

2523 else

2524 do k=2,max_nk

2525 if (k <= nk_in_ml) then

2526 z_t = z_t + hvel(k-1)

2527

2528 temp1 = (z_t * h_ml - z_t * z_t)

2529 ! This viscosity is set to go to 0 at the mixed layer top and bottom (in a log-layer)

2530 ! and be further limited by rotation to give the natural Ekman length.

2531 ! The following expressions are mathematically equivalent.

2532 if (gv%Boussinesq .or. (cs%answer_date < 20230601)) then

2533 visc_ml = u_star * cs%vonKar * (gv%Z_to_H * temp1 * u_star) &

2534 / (absf * temp1 + (h_ml + h_neglect) * u_star)

2535 else

2536 visc_ml = cs%vonKar * (temp1 * tau_mag) &

2537 / (absf * temp1 + (h_ml + h_neglect) * u_star)

2538 endif

2539 a_ml = visc_ml / (0.25 * (hvel(k) + hvel(k-1) + h_neglect) + 0.5 * i_amax * visc_ml)

2540

2541 ! Choose the largest estimate of a_cpl, but these could be changed to be additive.

2542 a_cpl(k) = max(a_cpl(k), a_ml)

2543 ! An option could be added to change this to: a_cpl(i,K) = a_cpl(i,K) + a_ml

2544 endif

2545 enddo

2546 endif

2547 endif

2548end subroutine find_coupling_coef

2549

2550

2551!> Velocity components which exceed a threshold for physically reasonable values are truncated,

2552!! and the running sum of the number of trunctionas within the non-symmetric memory computational

2553!! domain is incremented. Optionally, any column with excessive velocities may be sent

2554!! to a diagnostic reporting subroutine.

2555subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS)

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

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

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

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

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

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

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

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

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

2565 type(accel_diag_ptrs), intent(in) :: adp !< Acceleration diagnostic pointers

2566 type(cont_diag_ptrs), intent(in) :: cdp !< Continuity diagnostic pointers

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

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

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

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

2571

2572 ! Local variables

2573 real :: cfl ! The local CFL number [nondim]

2574 real :: h_report ! A thickness below which not to report truncations [H ~> m or kg m-2]

2575 real :: vel_report(szib_(g),szjb_(g)) ! The velocity to report [L T-1 ~> m s-1]

2576 real :: u_old(szib_(g),szj_(g),szk_(gv)) ! The previous u-velocity [L T-1 ~> m s-1]

2577 real :: v_old(szi_(g),szjb_(g),szk_(gv)) ! The previous v-velocity [L T-1 ~> m s-1]

2578 logical :: trunc_any, dowrite(szib_(g),szjb_(g))

2579 logical :: do_any_write

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

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

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

2583

2584 h_report = 3.0 * gv%Angstrom_H

2585

2586 if (len_trim(cs%u_trunc_file) > 0) then

2587 do_any_write = .false.

2588 trunc_any = .false.

2589

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

2591 dowrite(i,j) = .false.

2592 vel_report(i,j) = 3.0e8 * us%m_s_to_L_T

2593 enddo ; enddo

2594

2595 do k=1,nz ; do j=js,je ; do i=isq,ieq

2596 if (abs(u(i,j,k)) < cs%vel_underflow) u(i,j,k) = 0.0

2597 if (u(i,j,k) < 0.0) then

2598 cfl = (-u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))

2599 else

2600 cfl = (u(i,j,k) * dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))

2601 endif

2602 if (cfl > cs%CFL_trunc) trunc_any = .true.

2603 if (cfl > cs%CFL_report) then

2604 dowrite(i,j) = .true.

2605 do_any_write = .true.

2606 vel_report(i,j) = min(vel_report(i,j), abs(u(i,j,k)))

2607 endif

2608 enddo ; enddo ; enddo

2609

2610 do j=js,je ; do i=isq,ieq ; if (dowrite(i,j)) then

2611 u_old(i,j,:) = u(i,j,:)

2612 endif ; enddo ; enddo

2613

2614 if (trunc_any) then

2615 do k=1,nz ; do j=js,je ; do i=isq,ieq

2616 if ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL_trunc) then

2617 u(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i+1,j) / (dt * g%dy_Cu(i,j)))

2618 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2619 (cs%h_u(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2620 elseif ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then

2621 u(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dy_Cu(i,j)))

2622 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2623 (cs%h_u(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2624 endif

2625 enddo ; enddo ; enddo

2626 endif

2627

2628 if (do_any_write) then

2629 do j=js,je ; do i=isq,ieq ; if (dowrite(i,j)) then

2630 ! Call a diagnostic reporting subroutines are called if unphysically large values are found.

2631 call write_u_accel(i, j, u_old, h, adp, cdp, dt, g, gv, us, cs%PointAccel_CSp, &

2632 vel_report(i,j), forces%taux(i,j), a=cs%a_u, hv=cs%h_u)

2633 endif ; enddo ; enddo

2634 endif

2635 else ! Do not report accelerations leading to large velocities.

2636 do k=1,nz ; do j=js,je ; do i=isq,ieq

2637 if (abs(u(i,j,k)) < cs%vel_underflow) then ; u(i,j,k) = 0.0

2638 elseif ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i+1,j) < -cs%CFL_trunc) then

2639 u(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i+1,j) / (dt * g%dy_Cu(i,j)))

2640 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2641 (cs%h_u(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2642 elseif ((u(i,j,k) * (dt * g%dy_Cu(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then

2643 u(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dy_Cu(i,j)))

2644 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2645 (cs%h_u(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2646 endif

2647 enddo ; enddo ; enddo

2648 endif

2649

2650 if (len_trim(cs%v_trunc_file) > 0) then

2651 do_any_write =.false.

2652 trunc_any = .false.

2653

2654

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

2656 dowrite(i,j) = .false.

2657 vel_report(i,j) = 3.0e8 * us%m_s_to_L_T

2658 enddo ; enddo

2659

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

2661 if (abs(v(i,j,k)) < cs%vel_underflow) v(i,j,k) = 0.0

2662 if (v(i,j,k) < 0.0) then

2663 cfl = (-v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))

2664 else

2665 cfl = (v(i,j,k) * dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))

2666 endif

2667 if (cfl > cs%CFL_trunc) trunc_any = .true.

2668 if (cfl > cs%CFL_report) then

2669 dowrite(i,j) = .true.

2670 do_any_write = .true.

2671 vel_report(i,j) = min(vel_report(i,j), abs(v(i,j,k)))

2672 endif

2673 enddo ; enddo ; enddo

2674

2675 do j=jsq,jeq ; do i=is,ie ; if (dowrite(i,j)) then

2676 v_old(i,j,:) = v(i,j,:)

2677 endif ; enddo ; enddo

2678

2679 if (trunc_any) then

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

2681 if ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL_trunc) then

2682 v(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i,j+1) / (dt * g%dx_Cv(i,j)))

2683 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2684 (cs%h_v(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2685 elseif ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then

2686 v(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dx_Cv(i,j)))

2687 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2688 (cs%h_v(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2689 endif

2690 enddo ; enddo ; enddo

2691 endif

2692

2693 if (do_any_write) then

2694 do j=jsq,jeq ; do i=is,ie ; if (dowrite(i,j)) then

2695 ! Call a diagnostic reporting subroutines are called if unphysically large values are found.

2696 call write_v_accel(i, j, v_old, h, adp, cdp, dt, g, gv, us, cs%PointAccel_CSp, &

2697 vel_report(i,j), forces%tauy(i,j), a=cs%a_v, hv=cs%h_v)

2698 endif ; enddo ; enddo

2699 endif

2700 else ! Do not report accelerations leading to large velocities.

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

2702 if (abs(v(i,j,k)) < cs%vel_underflow) then ; v(i,j,k) = 0.0

2703 elseif ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j+1) < -cs%CFL_trunc) then

2704 v(i,j,k) = (-0.9*cs%CFL_trunc) * (g%areaT(i,j+1) / (dt * g%dx_Cv(i,j)))

2705 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2706 (cs%h_v(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2707 elseif ((v(i,j,k) * (dt * g%dx_Cv(i,j))) * g%IareaT(i,j) > cs%CFL_trunc) then

2708 v(i,j,k) = (0.9*cs%CFL_trunc) * (g%areaT(i,j) / (dt * g%dx_Cv(i,j)))

2709 if (((i >= g%isc) .and. (i <= g%iec) .and. (j >= g%jsc) .and. (j <= g%jec)) .and. &

2710 (cs%h_v(i,j,k) > h_report)) cs%ntrunc = cs%ntrunc + 1

2711 endif

2712 enddo ; enddo ; enddo

2713 endif

2714

2715end subroutine vertvisc_limit_vel

2716

2717

2718!> Initialize the vertical friction module

2719subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &

2720 ntrunc, CS, fpmix)

2721 type(ocean_internal_state), &

2722 target, intent(in) :: mis !< The "MOM Internal State", a set of pointers

2723 !! to the fields and accelerations that make

2724 !! up the ocean's physical state

2725 type(time_type), target, intent(in) :: time !< Current model time

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

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

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

2729 type(param_file_type), intent(in) :: param_file !< File to parse for parameters

2730 type(diag_ctrl), target, intent(inout) :: diag !< Diagnostic control structure

2731 type(accel_diag_ptrs), intent(inout) :: adp !< Acceleration diagnostic pointers

2732 type(directories), intent(in) :: dirs !< Relevant directory paths

2733 integer, target, intent(inout) :: ntrunc !< Number of velocity truncations

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

2735 logical, optional, intent(in) :: fpmix !< Nonlocal momentum mixing

2736

2737 ! Local variables

2738

2739 real :: kv_bbl ! A viscosity in the bottom boundary layer with a simple scheme [H Z T-1 ~> m2 s-1 or Pa s]

2740 real :: kv_back_z ! A background kinematic viscosity [Z2 T-1 ~> m2 s-1]

2741 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.

2742 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz

2743 logical :: lfpmix

2744 character(len=200) :: kappa_gl90_file, inputdir, kdgl90_varname

2745 ! This include declares and sets the variable "version".

2746# include "version_variable.h"

2747 character(len=40) :: mdl = "MOM_vert_friction" ! This module's name.

2748 character(len=40) :: thickness_units

2749 real :: kv_mks ! KVML in MKS [m2 s-1]

2750

2751 if (associated(cs)) then

2752 call mom_error(warning, "vertvisc_init called with an associated "// &

2753 "control structure.")

2754 return

2755 endif

2756 allocate(cs)

2757

2758 cs%initialized = .true.

2759

2760 if (gv%Boussinesq) then ; thickness_units = "m"

2761 else ; thickness_units = "kg m-2" ; endif

2762

2763 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke

2764 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

2765

2766 cs%diag => diag ; cs%ntrunc => ntrunc ; ntrunc = 0

2767

2768 lfpmix = .false.

2769 if (present(fpmix)) lfpmix = fpmix

2770

2771! Default, read and log parameters

2772 call log_version(param_file, mdl, version, "", log_to_all=.true., debugging=.true.)

2773 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &

2774 "This sets the default value for the various _ANSWER_DATE parameters.", &

2775 default=99991231)

2776 call get_param(param_file, mdl, "VERT_FRICTION_ANSWER_DATE", cs%answer_date, &

2777 "The vintage of the order of arithmetic and expressions in the viscous "//&

2778 "calculations. Values below 20190101 recover the answers from the end of 2018, "//&

2779 "while higher values use expressions that do not use an arbitrary hard-coded "//&

2780 "maximum viscous coupling coefficient between layers. Values below 20230601 "//&

2781 "recover a form of the viscosity within the mixed layer that breaks up the "//&

2782 "magnitude of the wind stress in some non-Boussinesq cases.", &

2783 default=default_answer_date, do_not_log=.not.gv%Boussinesq)

2784 if (.not.gv%Boussinesq) cs%answer_date = max(cs%answer_date, 20230701)

2785

2786 call get_param(param_file, mdl, "BOTTOMDRAGLAW", cs%bottomdraglaw, &

2787 "If true, the bottom stress is calculated with a drag "//&

2788 "law of the form c_drag*|u|*u. The velocity magnitude "//&

2789 "may be an assumed value or it may be based on the "//&

2790 "actual velocity in the bottommost HBBL, depending on "//&

2791 "LINEAR_DRAG.", default=.true.)

2792 call get_param(param_file, mdl, "DIRECT_STRESS", cs%direct_stress, &

2793 "If true, the wind stress is distributed over the topmost HMIX_STRESS of fluid "//&

2794 "(like in HYCOM), and an added mixed layer viscosity or a physically based "//&

2795 "boundary layer turbulence parameterization is not needed for stability.", &

2796 default=.false.)

2797 call get_param(param_file, mdl, "DYNAMIC_VISCOUS_ML", cs%dynamic_viscous_ML, &

2798 "If true, use a bulk Richardson number criterion to "//&

2799 "determine the mixed layer thickness for viscosity.", &

2800 default=.false.)

2801 call get_param(param_file, mdl, "FIXED_DEPTH_LOTW_ML", cs%fixed_LOTW_ML, &

2802 "If true, use a Law-of-the-wall prescription for the mixed layer viscosity "//&

2803 "within a boundary layer that is the lesser of HMIX_FIXED and the total "//&

2804 "depth of the ocean in a column.", default=.false.)

2805 call get_param(param_file, mdl, "LOTW_VISCOUS_ML_FLOOR", cs%apply_LOTW_floor, &

2806 "If true, use a Law-of-the-wall prescription to set a lower bound on the "//&

2807 "viscous coupling between layers within the surface boundary layer, based "//&

2808 "the distance of interfaces from the surface. This only acts when there "//&

2809 "are large changes in the thicknesses of successive layers or when the "//&

2810 "viscosity is set externally and the wind stress has subsequently increased.", &

2811 default=.false.)

2812 call get_param(param_file, mdl, 'VON_KARMAN_CONST', cs%vonKar, &

2813 'The value the von Karman constant as used for mixed layer viscosity.', &

2814 units='nondim', default=0.41)

2815 call get_param(param_file, mdl, "U_TRUNC_FILE", cs%u_trunc_file, &

2816 "The absolute path to a file into which the accelerations "//&

2817 "leading to zonal velocity truncations are written. "//&

2818 "Undefine this for efficiency if this diagnostic is not needed.", &

2819 default=" ", debuggingparam=.true.)

2820 call get_param(param_file, mdl, "V_TRUNC_FILE", cs%v_trunc_file, &

2821 "The absolute path to a file into which the accelerations "//&

2822 "leading to meridional velocity truncations are written. "//&

2823 "Undefine this for efficiency if this diagnostic is not needed.", &

2824 default=" ", debuggingparam=.true.)

2825 call get_param(param_file, mdl, "HARMONIC_VISC", cs%harmonic_visc, &

2826 "If true, use the harmonic mean thicknesses for "//&

2827 "calculating the vertical viscosity.", default=.false.)

2828 call get_param(param_file, mdl, "HARMONIC_BL_SCALE", cs%harm_BL_val, &

2829 "A scale to determine when water is in the boundary "//&

2830 "layers based solely on harmonic mean thicknesses for "//&

2831 "the purpose of determining the extent to which the "//&

2832 "thicknesses used in the viscosities are upwinded.", &

2833 default=0.0, units="nondim")

2834 call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)

2835

2836 if (gv%nkml < 1) then

2837 call get_param(param_file, mdl, "HMIX_FIXED", cs%Hmix, &

2838 "The prescribed depth over which the near-surface viscosity and "//&

2839 "diffusivity are elevated when the bulk mixed layer is not used.", &

2840 units="m", scale=us%m_to_Z, fail_if_missing=.true.)

2841 endif

2842 if (cs%direct_stress) then

2843 if (gv%nkml < 1) then

2844 call get_param(param_file, mdl, "HMIX_STRESS", cs%Hmix_stress, &

2845 "The depth over which the wind stress is applied if DIRECT_STRESS is true.", &

2846 units="m", default=us%Z_to_m*cs%Hmix, scale=gv%m_to_H)

2847 else

2848 call get_param(param_file, mdl, "HMIX_STRESS", cs%Hmix_stress, &

2849 "The depth over which the wind stress is applied if DIRECT_STRESS is true.", &

2850 units="m", fail_if_missing=.true., scale=gv%m_to_H)

2851 endif

2852 if (cs%Hmix_stress <= 0.0) call mom_error(fatal, "vertvisc_init: " // &

2853 "HMIX_STRESS must be set to a positive value if DIRECT_STRESS is true.")

2854 endif

2855 call get_param(param_file, mdl, "KV", kv_back_z, &

2856 "The background kinematic viscosity in the interior. "//&

2857 "The molecular value, ~1e-6 m2 s-1, may be used.", &

2858 units="m2 s-1", fail_if_missing=.true., scale=us%m2_s_to_Z2_T)

2859 ! Convert input kinematic viscosity to dynamic viscosity when non-Boussinesq.

2860 cs%Kv = (us%Z2_T_to_m2_s*gv%m2_s_to_HZ_T) * kv_back_z

2861

2862 call get_param(param_file, mdl, "USE_GL90_IN_SSW", cs%use_GL90_in_SSW, &

2863 "If true, use simpler method to calculate 1/N^2 in GL90 vertical "// &

2864 "viscosity coefficient. This method is valid in stacked shallow water mode.", &

2865 default=.false.)

2866 call get_param(param_file, mdl, "KD_GL90", cs%kappa_gl90, &

2867 "The scalar diffusivity used in GL90 vertical viscosity scheme.", &

2868 units="m2 s-1", default=0.0, scale=us%m_to_L*us%Z_to_L*gv%m_to_H*us%T_to_s, &

2869 do_not_log=.not.cs%use_GL90_in_SSW)

2870 call get_param(param_file, mdl, "READ_KD_GL90", cs%read_kappa_gl90, &

2871 "If true, read a file (given by KD_GL90_FILE) containing the "//&

2872 "spatially varying diffusivity KD_GL90 used in the GL90 scheme.", default=.false., &

2873 do_not_log=.not.cs%use_GL90_in_SSW)

2874 if (cs%read_kappa_gl90) then

2875 if (cs%kappa_gl90 > 0) then

2876 call mom_error(fatal, "MOM_vert_friction.F90, vertvisc_init: KD_GL90 > 0 "// &

2877 "is not compatible with READ_KD_GL90 = .TRUE. ")

2878 endif

2879 call get_param(param_file, mdl, "INPUTDIR", inputdir, &

2880 "The directory in which all input files are found.", &

2881 default=".", do_not_log=.true.)

2882 inputdir = slasher(inputdir)

2883 call get_param(param_file, mdl, "KD_GL90_FILE", kappa_gl90_file, &

2884 "The file containing the spatially varying diffusivity used in the "// &

2885 "GL90 scheme.", default="kd_gl90.nc", do_not_log=.not.cs%use_GL90_in_SSW)

2886 call get_param(param_file, mdl, "KD_GL90_VARIABLE", kdgl90_varname, &

2887 "The name of the GL90 diffusivity variable to read "//&

2888 "from KD_GL90_FILE.", default="kd_gl90", do_not_log=.not.cs%use_GL90_in_SSW)

2889 kappa_gl90_file = trim(inputdir) // trim(kappa_gl90_file)

2890

2891 allocate(cs%kappa_gl90_2d(g%isd:g%ied, g%jsd:g%jed), source=0.0)

2892 call mom_read_data(kappa_gl90_file, kdgl90_varname, cs%kappa_gl90_2d(:,:), g%domain, &

2893 scale=us%m_to_L*us%Z_to_L*gv%m_to_H*us%T_to_s)

2894 call pass_var(cs%kappa_gl90_2d, g%domain)

2895 endif

2896 call get_param(param_file, mdl, "USE_GL90_N2", cs%use_GL90_N2, &

2897 "If true, use GL90 vertical viscosity coefficient that is depth-independent; "// &

2898 "this corresponds to a kappa_GM that scales as N^2 with depth.", &

2899 default=.false., do_not_log=.not.cs%use_GL90_in_SSW)

2900 if (cs%use_GL90_N2) then

2901 if (.not. cs%use_GL90_in_SSW) call mom_error(fatal, &

2902 "MOM_vert_friction.F90, vertvisc_init: "//&

2903 "When USE_GL90_N2=True, USE_GL90_in_SSW must also be True.")

2904 if (cs%kappa_gl90 > 0) then

2905 call mom_error(fatal, "MOM_vert_friction.F90, vertvisc_init: KD_GL90 > 0 "// &

2906 "is not compatible with USE_GL90_N2 = .TRUE. ")

2907 endif

2908 if (cs%read_kappa_gl90) call mom_error(fatal, &

2909 "MOM_vert_friction.F90, vertvisc_init: "//&

2910 "READ_KD_GL90 = .TRUE. is not compatible with USE_GL90_N2 = .TRUE.")

2911 call get_param(param_file, mdl, "alpha_GL90", cs%alpha_gl90, &

2912 "Coefficient used to compute a depth-independent GL90 vertical "//&

2913 "viscosity via Kv_GL90 = alpha_GL90 * f2. Is only used "// &

2914 "if USE_GL90_N2 is true. Note that the implied Kv_GL90 "// &

2915 "corresponds to a KD_GL90 that scales as N^2 with depth.", &

2916 units="m2 s", default=0.0, scale=gv%m_to_H*us%m_to_Z*us%s_to_T, &

2917 do_not_log=.not.cs%use_GL90_in_SSW)

2918 endif

2919 call get_param(param_file, mdl, "HBBL_GL90", cs%Hbbl_gl90, &

2920 "The thickness of the GL90 bottom boundary layer, "//&

2921 "which defines the range over which the GL90 coupling "//&

2922 "coefficient is zeroed out, in order to avoid fluxing "//&

2923 "momentum into vanished layers over steep topography.", &

2924 units="m", default=5.0, scale=us%m_to_Z, do_not_log=.not.cs%use_GL90_in_SSW)

2925

2926 cs%Kvml_invZ2 = 0.0

2927 if (gv%nkml < 1) then

2928 call get_param(param_file, mdl, "KVML", kv_mks, &

2929 "The scale for an extra kinematic viscosity in the mixed layer", &

2930 units="m2 s-1", default=-1.0, do_not_log=.true.)

2931 if (kv_mks >= 0.0) then

2932 call mom_error(warning, "KVML is a deprecated parameter. Use KV_ML_INVZ2 instead.")

2933 else

2934 kv_mks = 0.0

2935 endif

2936 call get_param(param_file, mdl, "KV_ML_INVZ2", cs%Kvml_invZ2, &

2937 "An extra kinematic viscosity in a mixed layer of thickness HMIX_FIXED, "//&

2938 "with the actual viscosity scaling as 1/(z*HMIX_FIXED)^2, where z is the "//&

2939 "distance from the surface, to allow for finite wind stresses to be "//&

2940 "transmitted through infinitesimally thin surface layers. This is an "//&

2941 "older option for numerical convenience without a strong physical basis, "//&

2942 "and its use is now discouraged.", &

2943 units="m2 s-1", default=kv_mks, scale=gv%m2_s_to_HZ_T)

2944 endif

2945

2946 if (.not.cs%bottomdraglaw) then

2947 call get_param(param_file, mdl, "KV_EXTRA_BBL", cs%Kv_extra_bbl, &

2948 "An extra kinematic viscosity in the benthic boundary layer. "//&

2949 "KV_EXTRA_BBL is not used if BOTTOMDRAGLAW is true.", &

2950 units="m2 s-1", default=0.0, scale=gv%m2_s_to_HZ_T, do_not_log=.true.)

2951 if (cs%Kv_extra_bbl == 0.0) then

2952 call get_param(param_file, mdl, "KVBBL", kv_bbl, &

2953 "An extra kinematic viscosity in the benthic boundary layer. "//&

2954 "KV_EXTRA_BBL is not used if BOTTOMDRAGLAW is true.", &

2955 units="m2 s-1", default=us%Z2_T_to_m2_s*kv_back_z, scale=gv%m2_s_to_HZ_T, &

2956 do_not_log=.true.)

2957 if (abs(kv_bbl - cs%Kv) > 1.0e-15*abs(cs%Kv)) then

2958 call mom_error(warning, "KVBBL is a deprecated parameter. Use KV_EXTRA_BBL instead.")

2959 cs%Kv_extra_bbl = kv_bbl - cs%Kv

2960 endif

2961 endif

2962 call log_param(param_file, mdl, "KV_EXTRA_BBL", cs%Kv_extra_bbl, &

2963 "An extra kinematic viscosity in the benthic boundary layer. "//&

2964 "KV_EXTRA_BBL is not used if BOTTOMDRAGLAW is true.", &

2965 units="m2 s-1", default=0.0, unscale=gv%HZ_T_to_m2_s)

2966 endif

2967 call get_param(param_file, mdl, "HBBL", cs%Hbbl, &

2968 "The thickness of a bottom boundary layer with a viscosity increased by "//&

2969 "KV_EXTRA_BBL if BOTTOMDRAGLAW is not defined, or the thickness over which "//&

2970 "near-bottom velocities are averaged for the drag law if BOTTOMDRAGLAW is "//&

2971 "defined but LINEAR_DRAG is not.", &

2972 units="m", fail_if_missing=.true., scale=us%m_to_Z)

2973 call get_param(param_file, mdl, "CFL_TRUNCATE", cs%CFL_trunc, &

2974 "The value of the CFL number that will cause velocity "//&

2975 "components to be truncated; instability can occur past 0.5.", &

2976 units="nondim", default=0.5)

2977 call get_param(param_file, mdl, "CFL_REPORT", cs%CFL_report, &

2978 "The value of the CFL number that causes accelerations "//&

2979 "to be reported; the default is CFL_TRUNCATE.", &

2980 units="nondim", default=cs%CFL_trunc)

2981 call get_param(param_file, mdl, "CFL_TRUNCATE_RAMP_TIME", cs%truncRampTime, &

2982 "The time over which the CFL truncation value is ramped "//&

2983 "up at the beginning of the run.", &

2984 units="s", default=0., scale=us%s_to_T)

2985 cs%CFL_truncE = cs%CFL_trunc

2986 call get_param(param_file, mdl, "CFL_TRUNCATE_START", cs%CFL_truncS, &

2987 "The start value of the truncation CFL number used when "//&

2988 "ramping up CFL_TRUNC.", &

2989 units="nondim", default=0.)

2990 call get_param(param_file, mdl, "STOKES_MIXING_COMBINED", cs%StokesMixing, &

2991 "Flag to use Stokes drift Mixing via the Lagrangian "//&

2992 " current (Eulerian plus Stokes drift). "//&

2993 " Still needs work and testing, so not recommended for use.",&

2994 default=.false.)

2995 !BGR 04/04/2018{

2996 ! StokesMixing is required for MOM6 for some Langmuir mixing parameterization.

2997 ! The code used here has not been developed for vanishing layers or in

2998 ! conjunction with any bottom friction. Therefore, the following line is

2999 ! added so this functionality cannot be used without user intervention in

3000 ! the code. This will prevent general use of this functionality until proper

3001 ! care is given to the previously mentioned issues. Comment out the following

3002 ! MOM_error to use, but do so at your own risk and with these points in mind.

3003 !}

3004 if (cs%StokesMixing) then

3005 call mom_error(fatal, "Stokes mixing requires user intervention in the code.\n"//&

3006 " Model now exiting. See MOM_vert_friction.F90 for \n"//&

3007 " details (search 'BGR 04/04/2018' to locate comment).")

3008 endif

3009 call get_param(param_file, mdl, "VEL_UNDERFLOW", cs%vel_underflow, &

3010 "A negligibly small velocity magnitude below which velocity "//&

3011 "components are set to 0. A reasonable value might be "//&

3012 "1e-30 m/s, which is less than an Angstrom divided by "//&

3013 "the age of the universe.", units="m s-1", default=0.0, scale=us%m_s_to_L_T)

3014

3015 alloc_(cs%a_u(isdb:iedb,jsd:jed,nz+1)) ; cs%a_u(:,:,:) = 0.0

3016 alloc_(cs%a_u_gl90(isdb:iedb,jsd:jed,nz+1)) ; cs%a_u_gl90(:,:,:) = 0.0

3017 alloc_(cs%h_u(isdb:iedb,jsd:jed,nz)) ; cs%h_u(:,:,:) = 0.0

3018 alloc_(cs%a_v(isd:ied,jsdb:jedb,nz+1)) ; cs%a_v(:,:,:) = 0.0

3019 alloc_(cs%a_v_gl90(isd:ied,jsdb:jedb,nz+1)) ; cs%a_v_gl90(:,:,:) = 0.0

3020 alloc_(cs%h_v(isd:ied,jsdb:jedb,nz)) ; cs%h_v(:,:,:) = 0.0

3021

3022 cs%id_Kv_slow = register_diag_field('ocean_model', 'Kv_slow', diag%axesTi, time, &

3023 'Slow varying vertical viscosity', 'm2 s-1', conversion=gv%HZ_T_to_m2_s)

3024

3025 cs%id_Kv_u = register_diag_field('ocean_model', 'Kv_u', diag%axesCuL, time, &

3026 'Total vertical viscosity at u-points', 'm2 s-1', conversion=gv%H_to_m**2*us%s_to_T)

3027

3028 cs%id_Kv_v = register_diag_field('ocean_model', 'Kv_v', diag%axesCvL, time, &

3029 'Total vertical viscosity at v-points', 'm2 s-1', conversion=gv%H_to_m**2*us%s_to_T)

3030

3031 cs%id_Kv_gl90_u = register_diag_field('ocean_model', 'Kv_gl90_u', diag%axesCuL, time, &

3032 'GL90 vertical viscosity at u-points', 'm2 s-1', conversion=gv%H_to_m**2*us%s_to_T)

3033

3034 cs%id_Kv_gl90_v = register_diag_field('ocean_model', 'Kv_gl90_v', diag%axesCvL, time, &

3035 'GL90 vertical viscosity at v-points', 'm2 s-1', conversion=gv%H_to_m**2*us%s_to_T)

3036

3037 cs%id_au_vv = register_diag_field('ocean_model', 'au_visc', diag%axesCui, time, &

3038 'Zonal Viscous Vertical Coupling Coefficient', 'm s-1', conversion=gv%H_to_m*us%s_to_T)

3039

3040 cs%id_av_vv = register_diag_field('ocean_model', 'av_visc', diag%axesCvi, time, &

3041 'Meridional Viscous Vertical Coupling Coefficient', 'm s-1', conversion=gv%H_to_m*us%s_to_T)

3042

3043 cs%id_au_gl90_vv = register_diag_field('ocean_model', 'au_gl90_visc', diag%axesCui, time, &

3044 'Zonal Viscous Vertical GL90 Coupling Coefficient', 'm s-1', conversion=gv%H_to_m*us%s_to_T)

3045

3046 cs%id_av_gl90_vv = register_diag_field('ocean_model', 'av_gl90_visc', diag%axesCvi, time, &

3047 'Meridional Viscous Vertical GL90 Coupling Coefficient', 'm s-1', conversion=gv%H_to_m*us%s_to_T)

3048

3049 cs%id_h_u = register_diag_field('ocean_model', 'Hu_visc', diag%axesCuL, time, &

3050 'Thickness at Zonal Velocity Points for Viscosity', &

3051 thickness_units, conversion=gv%H_to_MKS)

3052 ! Alternately, to always give this variable in 'm' use the following line instead:

3053 ! 'm', conversion=GV%H_to_m)

3054

3055 cs%id_h_v = register_diag_field('ocean_model', 'Hv_visc', diag%axesCvL, time, &

3056 'Thickness at Meridional Velocity Points for Viscosity', &

3057 thickness_units, conversion=gv%H_to_MKS)

3058

3059 cs%id_hML_u = register_diag_field('ocean_model', 'HMLu_visc', diag%axesCu1, time, &

3060 'Mixed Layer Thickness at Zonal Velocity Points for Viscosity', &

3061 thickness_units, conversion=us%Z_to_m)

3062

3063 cs%id_hML_v = register_diag_field('ocean_model', 'HMLv_visc', diag%axesCv1, time, &

3064 'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', &

3065 thickness_units, conversion=us%Z_to_m)

3066

3067 if (lfpmix) then

3068 cs%id_uE_h = register_diag_field('ocean_model', 'uE_h' , cs%diag%axesTL, &

3069 time, 'x-zonal Eulerian' , 'm s-1', conversion=us%L_T_to_m_s)

3070 cs%id_vE_h = register_diag_field('ocean_model', 'vE_h' , cs%diag%axesTL, &

3071 time, 'y-merid Eulerian' , 'm s-1', conversion=us%L_T_to_m_s)

3072 cs%id_uInc_h = register_diag_field('ocean_model','uInc_h',cs%diag%axesTL, &

3073 time, 'x-zonal Eulerian' , 'm s-1', conversion=us%L_T_to_m_s)

3074 cs%id_vInc_h = register_diag_field('ocean_model','vInc_h',cs%diag%axesTL, &

3075 time, 'x-zonal Eulerian' , 'm s-1', conversion=us%L_T_to_m_s)

3076 cs%id_uStk = register_diag_field('ocean_model', 'uStk' , cs%diag%axesTL, &

3077 time, 'x-FP du increment' , 'm s-1', conversion=us%L_T_to_m_s)

3078 cs%id_vStk = register_diag_field('ocean_model', 'vStk' , cs%diag%axesTL, &

3079 time, 'y-FP dv increment' , 'm s-1', conversion=us%L_T_to_m_s)

3080

3081 cs%id_FPtau2s = register_diag_field('ocean_model','Omega_tau2s',cs%diag%axesTi, &

3082 time, 'Stress direction from shear','radians')

3083 cs%id_FPtau2w = register_diag_field('ocean_model','Omega_tau2w',cs%diag%axesTi, &

3084 time, 'Stress direction from wind','radians')

3085 cs%id_uStk0 = register_diag_field('ocean_model', 'uStk0' , diag%axesT1, &

3086 time, 'Zonal Surface Stokes', 'm s-1', conversion=us%L_T_to_m_s)

3087 cs%id_vStk0 = register_diag_field('ocean_model', 'vStk0' , diag%axesT1, &

3088 time, 'Merid Surface Stokes', 'm s-1', conversion=us%L_T_to_m_s)

3089 endif

3090

3091 cs%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', diag%axesCuL, time, &

3092 'Zonal Acceleration from Vertical Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)

3093 if (cs%id_du_dt_visc > 0) call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)

3094 cs%id_dv_dt_visc = register_diag_field('ocean_model', 'dv_dt_visc', diag%axesCvL, time, &

3095 'Meridional Acceleration from Vertical Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)

3096 if (cs%id_dv_dt_visc > 0) call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)

3097 cs%id_GLwork = register_diag_field('ocean_model', 'GLwork', diag%axesTL, time, &

3098 'Sign-definite Kinetic Energy Source from GL90 Vertical Viscosity', &

3099 'm3 s-3', conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)

3100 cs%id_du_dt_visc_gl90 = register_diag_field('ocean_model', 'du_dt_visc_gl90', diag%axesCuL, time, &

3101 'Zonal Acceleration from GL90 Vertical Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)

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

3103 call safe_alloc_ptr(adp%du_dt_visc_gl90,isdb,iedb,jsd,jed,nz)

3104 call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)

3105 endif

3106 cs%id_dv_dt_visc_gl90 = register_diag_field('ocean_model', 'dv_dt_visc_gl90', diag%axesCvL, time, &

3107 'Meridional Acceleration from GL90 Vertical Viscosity', 'm s-2', conversion=us%L_T2_to_m_s2)

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

3109 call safe_alloc_ptr(adp%dv_dt_visc_gl90,isd,ied,jsdb,jedb,nz)

3110 call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)

3111 endif

3112 cs%id_du_dt_str = register_diag_field('ocean_model', 'du_dt_str', diag%axesCuL, time, &

3113 'Zonal Acceleration from Surface Wind Stresses', 'm s-2', conversion=us%L_T2_to_m_s2)

3114 if (cs%id_du_dt_str > 0) call safe_alloc_ptr(adp%du_dt_str,isdb,iedb,jsd,jed,nz)

3115 cs%id_dv_dt_str = register_diag_field('ocean_model', 'dv_dt_str', diag%axesCvL, time, &

3116 'Meridional Acceleration from Surface Wind Stresses', 'm s-2', conversion=us%L_T2_to_m_s2)

3117 if (cs%id_dv_dt_str > 0) call safe_alloc_ptr(adp%dv_dt_str,isd,ied,jsdb,jedb,nz)

3118

3119 cs%id_taux_bot = register_diag_field('ocean_model', 'taux_bot', diag%axesCu1, &

3120 time, 'Zonal Bottom Stress from Ocean to Earth', &

3121 'Pa', conversion=us%RZ_to_kg_m2*us%L_T2_to_m_s2)

3122 cs%id_tauy_bot = register_diag_field('ocean_model', 'tauy_bot', diag%axesCv1, &

3123 time, 'Meridional Bottom Stress from Ocean to Earth', &

3124 'Pa', conversion=us%RZ_to_kg_m2*us%L_T2_to_m_s2)

3125

3126 !CS%id_hf_du_dt_visc = register_diag_field('ocean_model', 'hf_du_dt_visc', diag%axesCuL, Time, &

3127 ! 'Fractional Thickness-weighted Zonal Acceleration from Vertical Viscosity', &

3128 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)

3129 !if (CS%id_hf_du_dt_visc > 0) then

3130 ! call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)

3131 ! call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)

3132 !endif

3133

3134 !CS%id_hf_dv_dt_visc = register_diag_field('ocean_model', 'hf_dv_dt_visc', diag%axesCvL, Time, &

3135 ! 'Fractional Thickness-weighted Meridional Acceleration from Vertical Viscosity', &

3136 ! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)

3137 !if (CS%id_hf_dv_dt_visc > 0) then

3138 ! call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)

3139 ! call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)

3140 !endif

3141

3142 cs%id_hf_du_dt_visc_2d = register_diag_field('ocean_model', 'hf_du_dt_visc_2d', diag%axesCu1, time, &

3143 'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Vertical Viscosity', &

3144 'm s-2', conversion=us%L_T2_to_m_s2)

3145 if (cs%id_hf_du_dt_visc_2d > 0) then

3146 call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)

3147 call safe_alloc_ptr(adp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)

3148 endif

3149

3150 cs%id_hf_dv_dt_visc_2d = register_diag_field('ocean_model', 'hf_dv_dt_visc_2d', diag%axesCv1, time, &

3151 'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Vertical Viscosity', &

3152 'm s-2', conversion=us%L_T2_to_m_s2)

3153 if (cs%id_hf_dv_dt_visc_2d > 0) then

3154 call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)

3155 call safe_alloc_ptr(adp%diag_hfrac_v,isd,ied,jsdb,jedb,nz)

3156 endif

3157

3158 cs%id_h_du_dt_visc = register_diag_field('ocean_model', 'h_du_dt_visc', diag%axesCuL, time, &

3159 'Thickness Multiplied Zonal Acceleration from Horizontal Viscosity', &

3160 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

3161 if (cs%id_h_du_dt_visc > 0) then

3162 call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)

3163 call safe_alloc_ptr(adp%diag_hu,isdb,iedb,jsd,jed,nz)

3164 endif

3165

3166 cs%id_h_dv_dt_visc = register_diag_field('ocean_model', 'h_dv_dt_visc', diag%axesCvL, time, &

3167 'Thickness Multiplied Meridional Acceleration from Horizontal Viscosity', &

3168 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

3169 if (cs%id_h_dv_dt_visc > 0) then

3170 call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)

3171 call safe_alloc_ptr(adp%diag_hv,isd,ied,jsdb,jedb,nz)

3172 endif

3173

3174 cs%id_h_du_dt_str = register_diag_field('ocean_model', 'h_du_dt_str', diag%axesCuL, time, &

3175 'Thickness Multiplied Zonal Acceleration from Surface Wind Stresses', &

3176 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

3177 if (cs%id_h_du_dt_str > 0) then

3178 call safe_alloc_ptr(adp%du_dt_str,isdb,iedb,jsd,jed,nz)

3179 call safe_alloc_ptr(adp%diag_hu,isdb,iedb,jsd,jed,nz)

3180 endif

3181

3182 cs%id_h_dv_dt_str = register_diag_field('ocean_model', 'h_dv_dt_str', diag%axesCvL, time, &

3183 'Thickness Multiplied Meridional Acceleration from Surface Wind Stresses', &

3184 'm2 s-2', conversion=gv%H_to_m*us%L_T2_to_m_s2)

3185 if (cs%id_h_dv_dt_str > 0) then

3186 call safe_alloc_ptr(adp%dv_dt_str,isd,ied,jsdb,jedb,nz)

3187 call safe_alloc_ptr(adp%diag_hv,isd,ied,jsdb,jedb,nz)

3188 endif

3189

3190 cs%id_du_dt_str_visc_rem = register_diag_field('ocean_model', 'du_dt_str_visc_rem', diag%axesCuL, time, &

3191 'Zonal Acceleration from Surface Wind Stresses multiplied by viscous remnant', &

3192 'm s-2', conversion=us%L_T2_to_m_s2)

3193 if (cs%id_du_dt_str_visc_rem > 0) then

3194 call safe_alloc_ptr(adp%du_dt_str,isdb,iedb,jsd,jed,nz)

3195 call safe_alloc_ptr(adp%visc_rem_u,isdb,iedb,jsd,jed,nz)

3196 endif

3197

3198 cs%id_dv_dt_str_visc_rem = register_diag_field('ocean_model', 'dv_dt_str_visc_rem', diag%axesCvL, time, &

3199 'Meridional Acceleration from Surface Wind Stresses multiplied by viscous remnant', &

3200 'm s-2', conversion=us%L_T2_to_m_s2)

3201 if (cs%id_dv_dt_str_visc_rem > 0) then

3202 call safe_alloc_ptr(adp%dv_dt_str,isd,ied,jsdb,jedb,nz)

3203 call safe_alloc_ptr(adp%visc_rem_v,isd,ied,jsdb,jedb,nz)

3204 endif

3205

3206 if ((len_trim(cs%u_trunc_file) > 0) .or. (len_trim(cs%v_trunc_file) > 0)) &

3207 call pointaccel_init(mis, time, g, param_file, diag, dirs, cs%PointAccel_CSp)

3208

3209end subroutine vertvisc_init

3210

3211!> Update the CFL truncation value as a function of time.

3212!! If called with the optional argument activate=.true., record the

3213!! value of Time as the beginning of the ramp period.

3214subroutine updatecfltruncationvalue(Time, CS, US, activate)

3215 type(time_type), target, intent(in) :: time !< Current model time

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

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

3218 logical, optional, intent(in) :: activate !< Specify whether to record the value of

3219 !! Time as the beginning of the ramp period

3220

3221 ! Local variables

3222 real :: deltatime ! The time since CS%rampStartTime [T ~> s], which may be negative.

3223 real :: wghta ! The relative weight of the final value [nondim]

3224 character(len=12) :: msg

3225

3226 if (cs%truncRampTime==0.) return ! This indicates to ramping is turned off

3227

3228 ! We use the optional argument to indicate this Time should be recorded as the

3229 ! beginning of the ramp-up period.

3230 if (present(activate)) then

3231 if (activate) then

3232 cs%rampStartTime = time ! Record the current time

3233 cs%CFLrampingIsActivated = .true.

3234 endif

3235 endif

3236 if (.not.cs%CFLrampingIsActivated) return

3237 deltatime = max(0., time_minus_signed(time, cs%rampStartTime, scale=us%s_to_T))

3238 if (deltatime >= cs%truncRampTime) then

3239 cs%CFL_trunc = cs%CFL_truncE

3240 cs%truncRampTime = 0. ! This turns off ramping after this call

3241 else

3242 wghta = min( 1., deltatime / cs%truncRampTime ) ! Linear profile in time

3243 !wghtA = wghtA*wghtA ! Convert linear profile to parabolic profile in time

3244 !wghtA = wghtA*wghtA*(3. - 2.*wghtA) ! Convert linear profile to cosine profile

3245 wghta = 1. - ( (1. - wghta)**2 ) ! Convert linear profile to inverted parabolic profile

3246 cs%CFL_trunc = cs%CFL_truncS + wghta * ( cs%CFL_truncE - cs%CFL_truncS )

3247 endif

3248 write(msg(1:12),'(es12.3)') cs%CFL_trunc

3249 call mom_error(note, "MOM_vert_friction: updateCFLtruncationValue set CFL limit to "//trim(msg))

3250end subroutine updatecfltruncationvalue

3251

3252!> Clean up and deallocate the vertical friction module

3253subroutine vertvisc_end(CS)

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

3255 !! will be deallocated in this subroutine.

3256

3257 if ((len_trim(cs%u_trunc_file) > 0) .or. (len_trim(cs%v_trunc_file) > 0)) &

3258 deallocate(cs%PointAccel_CSp)

3259

3260 dealloc_(cs%a_u) ; dealloc_(cs%h_u)

3261 dealloc_(cs%a_v) ; dealloc_(cs%h_v)

3262 if (associated(cs%a1_shelf_u)) deallocate(cs%a1_shelf_u)

3263 if (associated(cs%a1_shelf_v)) deallocate(cs%a1_shelf_v)

3264 if (allocated(cs%kappa_gl90_2d)) deallocate(cs%kappa_gl90_2d)

3265end subroutine vertvisc_end

3266

3267!> \namespace mom_vert_friction

3268!! \author Robert Hallberg

3269!! \date April 1994 - October 2006

3270!!

3271!! The vertical diffusion of momentum is fully implicit. This is

3272!! necessary to allow for vanishingly small layers. The coupling

3273!! is based on the distance between the centers of adjacent layers,

3274!! except where a layer is close to the bottom compared with a

3275!! bottom boundary layer thickness when a bottom drag law is used.

3276!! A stress top b.c. and a no slip bottom b.c. are used. There

3277!! is no limit on the time step for vertvisc.

3278!!

3279!! Near the bottom, the horizontal thickness interpolation scheme

3280!! changes to an upwind biased estimate to control the effect of

3281!! spurious Montgomery potential gradients at the bottom where

3282!! nearly massless layers layers ride over the topography. Within a

3283!! few boundary layer depths of the bottom, the harmonic mean

3284!! thickness (i.e. (2 h+ h-) / (h+ + h-) ) is used if the velocity

3285!! is from the thinner side and the arithmetic mean thickness

3286!! (i.e. (h+ + h-)/2) is used if the velocity is from the thicker

3287!! side. Both of these thickness estimates are second order

3288!! accurate. Above this the arithmetic mean thickness is used.

3289!!

3290!! In addition, vertvisc truncates any velocity component that exceeds a

3291!! maximum CFL number to a fraction of this value. This basically keeps

3292!! instabilities spatially localized. The number of times the velocity is

3293!! truncated is reported each time the energies are saved, and if

3294!! exceeds CS%Maxtrunc the model will stop itself and change the time

3295!! to a large value. This has proven very useful in (1) diagnosing

3296!! model failures and (2) letting the model settle down to a

3297!! meaningful integration from a poorly specified initial condition.

3298!!

3299!! The same code is used for the two velocity components, by

3300!! indirectly referencing the velocities and defining a handful of

3301!! direction-specific defined variables.

3302!!

3303!! Macros written all in capital letters are defined in MOM_memory.h.

3304!!

3305!! A small fragment of the grid is shown below:

3306!! \verbatim

3307!! j+1 x ^ x ^ x At x: q

3308!! j+1 > o > o > At ^: v, frhatv, tauy

3309!! j x ^ x ^ x At >: u, frhatu, taux

3310!! j > o > o > At o: h

3311!! j-1 x ^ x ^ x

3312!! i-1 i i+1 At x & ^:

3313!! i i+1 At > & o:

3314!! \endverbatim

3315!!

3316!! The boundaries always run through q grid points (x).

3317end module mom_vert_friction