src/parameterizations/vertical/MOM_vert_friction.F90

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

1		! This file is part of MOM6, the Modular Ocean Model version 6.
2		! See the LICENSE file for licensing information.
3		! SPDX-License-Identifier: Apache-2.0
4
5		#include "do_concurrent_compat.h"
6
7		!> Implements vertical viscosity (vertvisc)
8		module MOM_vert_friction
9
10		use MOM_domains, only : pass_var, To_All, Omit_corners
11		use MOM_domains, only : pass_vector, Scalar_Pair
12		use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
13		use MOM_diag_mediator, only : post_product_u, post_product_sum_u
14		use MOM_diag_mediator, only : post_product_v, post_product_sum_v
15		use MOM_diag_mediator, only : diag_ctrl, query_averaging_enabled
16		use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type
17		use MOM_domains, only : To_North, To_East
18		use MOM_debugging, only : uvchksum, hchksum
19		use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE
20		use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
21		use MOM_forcing_type, only : mech_forcing, find_ustar
22		use MOM_get_input, only : directories
23		use MOM_grid, only : ocean_grid_type
24		use MOM_io, only : MOM_read_data, slasher
25		use MOM_open_boundary, only : ocean_OBC_type, OBC_NONE, OBC_DIRECTION_E
26		use MOM_open_boundary, only : OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S
27		use MOM_PointAccel, only : write_u_accel, write_v_accel, PointAccel_init
28		use MOM_PointAccel, only : PointAccel_CS
29		use MOM_time_manager, only : time_type, time_minus_signed
30		use MOM_unit_scaling, only : unit_scale_type
31		use MOM_variables, only : thermo_var_ptrs, vertvisc_type
32		use MOM_variables, only : cont_diag_ptrs, accel_diag_ptrs
33		use MOM_variables, only : ocean_internal_state
34		use MOM_verticalGrid, only : verticalGrid_type
35		use MOM_wave_interface, only : wave_parameters_CS
36		use MOM_set_visc, only : set_v_at_u, set_u_at_v
37		use MOM_lateral_mixing_coeffs, only : VarMix_CS
38
39		use CVMix_kpp, only : cvmix_kpp_composite_Gshape
40
41		implicit none ; private
42
43		#include <MOM_memory.h>
44
45		public vertvisc, vertvisc_remnant, vertvisc_coef
46		public vertvisc_limit_vel, vertvisc_init, vertvisc_end
47		public updateCFLtruncationValue
48		public vertFPmix
49
50		! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
51		! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
52		! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
53		! vary with the Boussinesq approximation, the Boussinesq variant is given first.
54
55		!> The control structure with parameters and memory for the MOM_vert_friction module
56		type, public :: vertvisc_CS ; private
57		logical :: initialized = .false. !< True if this control structure has been initialized.
58		real :: Hmix !< The mixed layer thickness [Z ~> m].
59		real :: Hmix_stress !< The mixed layer thickness over which the wind
60		!! stress is applied with direct_stress [H ~> m or kg m-2].
61		real :: Kvml_invZ2 !< The extra vertical viscosity scale in [H Z T-1 ~> m2 s-1 or Pa s] in a
62		!! surface mixed layer with a characteristic thickness given by Hmix,
63		!! and scaling proportional to (Hmix/z)^2, where z is the distance
64		!! from the surface; this can get very large with thin layers.
65		real :: Kv !< The interior vertical viscosity [H Z T-1 ~> m2 s-1 or Pa s].
66		real :: Hbbl !< The static bottom boundary layer thickness [Z ~> m].
67		real :: Hbbl_gl90 !< The static bottom boundary layer thickness used for GL90 [Z ~> m].
68		real :: Kv_extra_bbl !< An extra vertical viscosity in the bottom boundary layer of thickness
69		!! Hbbl when there is not a bottom drag law in use [H Z T-1 ~> m2 s-1 or Pa s].
70		real :: vonKar !< The von Karman constant as used for mixed layer viscosity [nondim]
71
72		logical :: use_GL90_in_SSW !< If true, use the GL90 parameterization in stacked shallow water mode (SSW).
73		!! The calculation of the GL90 viscosity coefficient uses the fact that in SSW
74		!! we simply have 1/N^2 = h/g^prime, where g^prime is the reduced gravity.
75		!! This identity does not generalize to non-SSW setups.
76		logical :: use_GL90_N2 !< If true, use GL90 vertical viscosity coefficient that is depth-independent;
77		!! this corresponds to a kappa_GM that scales as N^2 with depth.
78		real :: kappa_gl90 !< The scalar diffusivity used in the GL90 vertical viscosity scheme
79		!! [L2 H Z-1 T-1 ~> m2 s-1 or Pa s]
80		logical :: read_kappa_gl90 !< If true, read a file containing the spatially varying kappa_gl90
81		real :: alpha_gl90 !< Coefficient used to compute a depth-independent GL90 vertical
82		!! viscosity via Kv_gl90 = alpha_gl90 * f^2. Note that the implied
83		!! Kv_gl90 corresponds to a kappa_gl90 that scales as N^2 with depth.
84		!! [H Z T ~> m2 s or kg s m-1]
85		real :: vel_underflow !< Velocity components smaller than vel_underflow
86		!! are set to 0 [L T-1 ~> m s-1].
87		real :: CFL_trunc !< Velocity components will be truncated when they
88		!! are large enough that the corresponding CFL number
89		!! exceeds this value [nondim].
90		real :: CFL_report !< The value of the CFL number that will cause the
91		!! accelerations to be reported [nondim]. CFL_report
92		!! will often equal CFL_trunc.
93		real :: truncRampTime !< The time-scale over which to ramp up the value of
94		!! CFL_trunc from CFL_truncS to CFL_truncE [T ~> s]
95		real :: CFL_truncS !< The start value of CFL_trunc [nondim]
96		real :: CFL_truncE !< The end/target value of CFL_trunc [nondim]
97		logical :: CFLrampingIsActivated = .false. !< True if the ramping has been initialized
98		type(time_type) :: rampStartTime !< The time at which the ramping of CFL_trunc starts
99
100		real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NK_INTERFACE_) :: &
101		a_u !< The u-drag coefficient across an interface [H T-1 ~> m s-1 or Pa s m-1]
102		real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NK_INTERFACE_) :: &
103		a_u_gl90 !< The u-drag coefficient associated with GL90 across an interface [H T-1 ~> m s-1 or Pa s m-1]
104		real ALLOCABLE_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
105		h_u !< The effective layer thickness at u-points [H ~> m or kg m-2].
106		real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NK_INTERFACE_) :: &
107		a_v !< The v-drag coefficient across an interface [H T-1 ~> m s-1 or Pa s m-1]
108		real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NK_INTERFACE_) :: &
109		a_v_gl90 !< The v-drag coefficient associated with GL90 across an interface [H T-1 ~> m s-1 or Pa s m-1]
110		real ALLOCABLE_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
111		h_v !< The effective layer thickness at v-points [H ~> m or kg m-2].
112		real, pointer, dimension(:,:) :: a1_shelf_u => NULL() !< The u-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		real, pointer, dimension(:,:) :: a1_shelf_v => NULL() !< The v-momentum coupling coefficient under
115		!! ice shelves [H T-1 ~> m s-1 or Pa s m-1]. Retained to determine stress under shelves.
116
117		logical :: split !< If true, use the split time stepping scheme.
118		logical :: bottomdraglaw !< If true, the bottom stress is calculated with a
119		!! drag law c_drag\|u\|u. The velocity magnitude
120		!! may be an assumed value or it may be based on the
121		!! actual velocity in the bottommost HBBL, depending
122		!! on whether linear_drag is true.
123		logical :: harmonic_visc !< If true, the harmonic mean thicknesses are used
124		!! to calculate the viscous coupling between layers
125		!! except near the bottom. Otherwise the arithmetic
126		!! mean thickness is used except near the bottom.
127		real :: harm_BL_val !< A scale to determine when water is in the boundary
128		!! layers based solely on harmonic mean thicknesses
129		!! for the purpose of determining the extent to which
130		!! the thicknesses used in the viscosities are upwinded [nondim].
131		logical :: direct_stress !< If true, the wind stress is distributed over the topmost Hmix_stress
132		!! of fluid, and an added mixed layer viscosity or a physically based
133		!! boundary layer turbulence parameterization is not needed for stability.
134		logical :: dynamic_viscous_ML !< If true, use the results from a dynamic
135		!! calculation, perhaps based on a bulk Richardson
136		!! number criterion, to determine the mixed layer
137		!! thickness for viscosity.
138		logical :: fixed_LOTW_ML !< If true, use a Law-of-the-wall prescription for the mixed layer
139		!! viscosity within a boundary layer that is the lesser of Hmix and the
140		!! total depth of the ocean in a column.
141		logical :: apply_LOTW_floor !< If true, use a Law-of-the-wall prescription to set a lower bound
142		!! on the viscous coupling between layers within the surface boundary
143		!! layer, based the distance of interfaces from the surface. This only
144		!! acts when there are large changes in the thicknesses of successive
145		!! layers or when the viscosity is set externally and the wind stress
146		!! has subsequently increased.
147		integer :: answer_date !< The vintage of the order of arithmetic and expressions in the viscous
148		!! calculations. Values below 20190101 recover the answers from the end
149		!! of 2018, while higher values use expressions that do not use an
150		!! arbitrary and hard-coded maximum viscous coupling coefficient between
151		!! layers. In non-Boussinesq cases, values below 20230601 recover a
152		!! form of the viscosity within the mixed layer that breaks up the
153		!! magnitude of the wind stress with BULKMIXEDLAYER, DYNAMIC_VISCOUS_ML
154		!! or FIXED_DEPTH_LOTW_ML, but not LOTW_VISCOUS_ML_FLOOR.
155		logical :: debug !< If true, write verbose checksums for debugging purposes.
156		integer :: nkml !< The number of layers in the mixed layer.
157		integer, pointer :: ntrunc !< The number of times the velocity has been
158		!! truncated since the last call to write_energy.
159		character(len=200) :: u_trunc_file !< The complete path to a file in which a column of
160		!! u-accelerations are written if velocity truncations occur.
161		character(len=200) :: v_trunc_file !< The complete path to a file in which a column of
162		!! v-accelerations are written if velocity truncations occur.
163		logical :: StokesMixing !< If true, do Stokes drift mixing via the Lagrangian current
164		!! (Eulerian plus Stokes drift). False by default and set
165		!! via STOKES_MIXING_COMBINED.
166
167		type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
168		!! timing of diagnostic output.
169		real, allocatable, dimension(:,:) :: kappa_gl90_2d !< 2D kappa_gl90 at h-points [L2 H Z-1 T-1 ~> m2 s-1 or Pa s]
170
171		!>@{ Diagnostic identifiers
172		integer :: id_du_dt_visc = -1, id_dv_dt_visc = -1, id_du_dt_visc_gl90 = -1, id_dv_dt_visc_gl90 = -1
173		integer :: id_GLwork = -1
174		integer :: id_au_vv = -1, id_av_vv = -1, id_au_gl90_vv = -1, id_av_gl90_vv = -1
175		integer :: id_du_dt_str = -1, id_dv_dt_str = -1
176		integer :: id_h_u = -1, id_h_v = -1, id_hML_u = -1 , id_hML_v = -1
177		integer :: id_Omega_w2x = -1, id_FPtau2s = -1 , id_FPtau2w = -1
178		integer :: id_uE_h = -1, id_vE_h = -1
179		integer :: id_uStk = -1, id_vStk = -1
180		integer :: id_uStk0 = -1, id_vStk0 = -1
181		integer :: id_uInc_h= -1, id_vInc_h= -1
182		integer :: id_taux_bot = -1, id_tauy_bot = -1
183		integer :: id_Kv_slow = -1, id_Kv_u = -1, id_Kv_v = -1
184		integer :: id_Kv_gl90_u = -1, id_Kv_gl90_v = -1
185		! integer :: id_hf_du_dt_visc = -1, id_hf_dv_dt_visc = -1
186		integer :: id_h_du_dt_visc = -1, id_h_dv_dt_visc = -1
187		integer :: id_hf_du_dt_visc_2d = -1, id_hf_dv_dt_visc_2d = -1
188		integer :: id_h_du_dt_str = -1, id_h_dv_dt_str = -1
189		integer :: id_du_dt_str_visc_rem = -1, id_dv_dt_str_visc_rem = -1
190		!>@}
191
192		type(PointAccel_CS), pointer :: PointAccel_CSp => NULL() !< A pointer to the control structure
193		!! for recording accelerations leading to velocity truncations
194
195		type(group_pass_type) :: pass_KE_uv !< A handle used for group halo passes
196		end type vertvisc_CS
197
198		contains
199
200		!> Add nonlocal stress increments to ui^n and vi^n.
201	0	subroutine vertFPmix(ui, vi, uold, vold, hbl_h, h, forces, dt, lpost, Cemp_NL, G, GV, US, CS, OBC, Waves)
202		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
203		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
204		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
205		intent(inout) :: ui !< Zonal velocity after vertvisc [L T-1 ~> m s-1]
206		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
207		intent(inout) :: vi !< Meridional velocity after vertvisc [L T-1 ~> m s-1]
208		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
209		intent(inout) :: uold !< Old Zonal velocity [L T-1 ~> m s-1]
210		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
211		intent(inout) :: vold !< Old Meridional velocity [L T-1 ~> m s-1]
212		real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: hbl_h !< boundary layer depth [H ~> m]
213		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
214		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
215		type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
216		real, intent(in) :: dt !< Time increment [T ~> s]
217		real, intent(in) :: Cemp_NL !< empirical coefficient of non-local momentum mixing [nondim]
218		logical, intent(in) :: lpost !< Compute and make available FPMix diagnostics
219		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
220		type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
221		type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure
222		type(wave_parameters_CS), &
223		optional, pointer :: Waves !< Container for wave/Stokes information
224
225		! local variables
226	0	real, dimension(SZIB_(G),SZJ_(G)) :: hbl_u !< boundary layer depth (u-pts) [H ~> m]
227	0	real, dimension(SZI_(G),SZJB_(G)) :: hbl_v !< boundary layer depth (v-pts) [H ~> m]
228	0	real, dimension(SZIB_(G),SZJ_(G)) :: taux_u !< kinematic zonal wind stress (u-pts) [L2 T-2 ~> m2 s-2]
229	0	real, dimension(SZI_(G),SZJB_(G)) :: tauy_v !< kinematic merid wind stress (v-pts) [L2 T-2 ~> m2 s-2]
230	0	real, dimension(SZI_(G),SZJ_(G)) :: uS0 !< surface zonal Stokes drift h-pts [L T-1 ~> m s-1]
231	0	real, dimension(SZI_(G),SZJ_(G)) :: vS0 !< surface zonal Stokes drift h-pts [L T-1 ~> m s-1]
232	0	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uE_u !< zonal Eulerian u-pts [L T-1 ~> m s-1]
233	0	real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: uE_h !< zonal Eulerian h-pts [L T-1 ~> m s-1]
234	0	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vE_v !< merid Eulerian v-pts [L T-1 ~> m s-1]
235	0	real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: vE_h !< merid Eulerian h-pts [L T-1 ~> m s-1]
236	0	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: uInc_u !< zonal Eulerian u-pts [L T-1 ~> m s-1]
237	0	real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: uInc_h !< zonal Eulerian h-pts [L T-1 ~> m s-1]
238	0	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: vInc_v !< merid Eulerian v-pts [L T-1 ~> m s-1]
239	0	real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: vInc_h !< merid Eulerian h-pts [L T-1 ~> m s-1]
240	0	real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: uStk !< zonal Stokes Drift (h-pts) [L T-1 ~> m s-1]
241	0	real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)) :: vStk !< merid Stokes Drift (h-pts) [L T-1 ~> m s-1]
242	0	real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)+1) :: omega_tau2s !< angle stress to shear (h-pts) [rad]
243	0	real, dimension(SZI_(G) ,SZJ_(G),SZK_(GV)+1) :: omega_tau2w !< angle stress to wind (h-pts) [rad]
244		real :: omega_tmp, omega_s2x, omega_tau2x !< temporary angle wrt the x axis [rad]
245		real :: Irho0 !< Inverse of the mean density rescaled to [Z L-1 R-1 ~> m3 kg-1]
246		real :: pi !< ! The ratio of the circumference of a circle to its diameter [nondim]
247		real :: tmp_u, tmp_v !< temporary ocean mask weights on u and v points [nondim]
248		real :: fexp !< temporary exponential function [nondim]
249		real :: sigma !< temporary normalize boundary layer coordinate [nondim]
250		real :: Gat1, Gsig, dGdsig !< Shape parameters [nondim]
251		real :: du, dv !< Intermediate velocity differences [L T-1 ~> m s-1]
252		real :: depth !< Cumulative of thicknesses [H ~> m]
253		integer :: b, kp1, k, nz !< band and vertical indices
254		integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq !< horizontal indices
255
256	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
257	0	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke
258
259	0	pi = 4. * atan2(1.,1.)
260	0	Irho0 = 1.0 / GV%Rho0
261
262	0	call pass_var(hbl_h , G%Domain, halo=1)
263
264		! u-points
265	0	do j = js,je
266	0	do I = Isq,Ieq
267	0	taux_u(I,j) = forces%taux(I,j) * Irho0
268	0	if ( (G%mask2dCu(I,j) > 0.5) ) then
269		! h to u-pts
270	0	tmp_u = MAX (1.0 ,(G%mask2dT(i,j) + G%mask2dT(i+1,j) ) )
271	0	hbl_u(I,j) = ((G%mask2dT(i,j) * hbl_h(i,j)) + (G%mask2dT(i+1,j) * hbl_h(i+1,j))) / tmp_u
272	0	depth = 0.
273	0	Gat1 = 0.
274	0	do k=1, nz
275		! cell center
276	0	depth = depth + 0.5*CS%h_u(I,j,k)
277	0	uE_u(I,j,k) = ui(I,j,k) - waves%Us_x(I,j,k)
278	0	if ( depth < hbl_u(I,j) ) then
279	0	sigma = depth / hbl_u(i,j)
280		! cell bottom
281	0	depth = depth + 0.5*CS%h_u(I,j,k)
282	0	call cvmix_kpp_composite_Gshape(sigma,Gat1,Gsig,dGdsig)
283		! nonlocal boundary-layer increment
284	0	uInc_u(I,j,k) = dt * Cemp_NL * taux_u(I,j) * dGdsig / hbl_u(I,j)
285	0	ui(I,j,k) = ui(I,j,k) + uInc_u(I,j,k)
286		else
287	0	uInc_u(I,j,k) = 0.0
288		endif
289		enddo
290		else
291	0	do k=1, nz
292	0	uInc_u(I,j,k) = 0.0
293		enddo
294		endif
295		enddo
296		enddo
297
298		! v-points
299	0	do J = Jsq,Jeq
300	0	do i = is,ie
301	0	tauy_v(i,J) = forces%tauy(i,J) * Irho0
302	0	if ( (G%mask2dCv(i,J) > 0.5) ) then
303		! h to v-pts
304	0	tmp_v = max( 1.0 ,(G%mask2dT(i,j) + G%mask2dT(i,j+1)))
305	0	hbl_v(i,J) = (G%mask2dT(i,j) * hbl_h(i,J) + G%mask2dT(i,j+1) * hbl_h(i,j+1)) / tmp_v
306	0	depth = 0.
307	0	Gat1 = 0.
308	0	do k=1, nz
309		! cell center
310	0	depth = depth + 0.5* CS%h_v(i,J,k)
311	0	vE_v(i,J,k) = vi(i,J,k) - waves%Us_y(i,J,k)
312	0	if ( depth < hbl_v(i,J) ) then
313	0	sigma = depth / hbl_v(i,J)
314		! cell bottom
315	0	depth = depth + 0.5* CS%h_v(i,J,k)
316	0	call cvmix_kpp_composite_Gshape(sigma,Gat1,Gsig,dGdsig)
317		! nonlocal boundary-layer increment
318	0	vInc_v(i,J,k) = dt * Cemp_NL * tauy_v(i,J) * dGdsig / hbl_v(i,J)
319	0	vi(i,J,k) = vi(i,J,k) + vInc_v(i,J,k)
320		else
321	0	vInc_v(i,J,k) = 0.0
322		endif
323		enddo
324		else
325	0	do k=1, nz
326	0	vInc_v(i,J,k) = 0.0
327		enddo
328		endif
329		enddo
330		enddo
331
332		! Compute and store diagnostics, only during the corrector step.
333	0	if (lpost) then
334	0	call pass_vector(uE_u , vE_v , G%Domain, To_All)
335	0	call pass_vector(uInc_u, vInc_v , G%Domain, To_All)
336	0	uStk = 0.0
337	0	vStk = 0.0
338	0	uS0 = 0.0
339	0	vS0 = 0.0
340
341	0	do j = js,je
342	0	do i = is,ie
343	0	if (G%mask2dT(i,j) > 0.5) then
344		! u to h-pts
345	0	tmp_u = max( 1.0 ,(G%mask2dCu(i,j) + G%mask2dCu(i-1,j)))
346		! v to h-pts
347	0	tmp_v = max( 1.0 ,(G%mask2dCv(i,j) + G%mask2dCv(i,j-1)))
348	0	do k = 1,nz
349	0	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
350	0	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
351	0	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
352	0	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
353		enddo
354		! Wind, Stress and Shear align at surface
355	0	Omega_tau2w(i,j,1) = 0.0
356	0	Omega_tau2s(i,j,1) = 0.0
357	0	do k = 1,nz
358	0	kp1 = min( nz , k+1)
359	0	du = uE_h(i,j,k) - uE_h(i,j,kp1)
360	0	dv = vE_h(i,j,k) - vE_h(i,j,kp1)
361	0	omega_s2x = atan2( dv , du )
362
363	0	du = du + uInc_h(i,j,k) - uInc_h(i,j,kp1)
364	0	dv = dv + vInc_h(i,j,k) - vInc_h(i,j,kp1)
365	0	omega_tau2x = atan2( dv , du )
366
367	0	omega_tmp = omega_tau2x - forces%omega_w2x(i,j)
368	0	if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi
369	0	if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi
370	0	Omega_tau2w(i,j,kp1) = omega_tmp
371
372	0	omega_tmp = omega_tau2x - omega_s2x
373	0	if ( (omega_tmp > pi ) ) omega_tmp = omega_tmp - 2.*pi
374	0	if ( (omega_tmp < (0.-pi)) ) omega_tmp = omega_tmp + 2.*pi
375	0	Omega_tau2s(i,j,kp1) = omega_tmp
376
377		enddo
378		endif
379
380		! Stokes drift
381	0	do b=1,waves%NumBands
382	0	uS0(i,j) = uS0(i,j) + waves%UStk_Hb(i,j,b) ! or forces%UStkb(i,j,b)
383	0	vS0(i,j) = vS0(i,j) + waves%VStk_Hb(i,j,b) ! or forces%VStkb(i,j,b)
384		enddo
385	0	depth = 0.0
386	0	do k = 1,nz
387	0	do b = 1, waves%NumBands
388		! cell center
389	0	fexp = exp(-2. * waves%WaveNum_Cen(b) * (depth+0.5*h(i,j,k)) )
390	0	uStk(i,j,k) = uStk(i,j,k) + waves%UStk_Hb(i,j,b) * fexp
391	0	vStk(i,j,k) = vStk(i,j,k) + waves%VStk_Hb(i,j,b) * fexp
392		enddo
393		! cell bottom
394	0	depth = depth + h(i,j,k)
395		enddo
396		enddo
397		enddo
398
399		! post FPmix diagnostics
400	0	if (CS%id_uE_h > 0) call post_data(CS%id_uE_h , uE_h , CS%diag)
401	0	if (CS%id_vE_h > 0) call post_data(CS%id_vE_h , vE_h , CS%diag)
402	0	if (CS%id_uInc_h > 0) call post_data(CS%id_uInc_h , uInc_h , CS%diag)
403	0	if (CS%id_vInc_h > 0) call post_data(CS%id_vInc_h , vInc_h , CS%diag)
404	0	if (CS%id_FPtau2s > 0) call post_data(CS%id_FPtau2s, Omega_tau2s, CS%diag)
405	0	if (CS%id_FPtau2w > 0) call post_data(CS%id_FPtau2w, Omega_tau2w, CS%diag)
406	0	if (CS%id_uStk0 > 0) call post_data(CS%id_uStk0 , uS0 , CS%diag)
407	0	if (CS%id_vStk0 > 0) call post_data(CS%id_vStk0 , vS0 , CS%diag)
408	0	if (CS%id_uStk > 0) call post_data(CS%id_uStk , uStk , CS%diag)
409	0	if (CS%id_vStk > 0) call post_data(CS%id_vStk , vStk , CS%diag)
410	0	if (CS%id_Omega_w2x > 0) call post_data(CS%id_Omega_w2x, forces%omega_w2x, CS%diag)
411
412		endif
413
414	0	end subroutine vertFPmix
415
416
417		!> Compute coupling coefficient associated with vertical viscosity parameterization as in Greatbatch and Lamb
418		!! (1990), hereafter referred to as the GL90 vertical viscosity parameterization. This vertical viscosity scheme
419		!! redistributes momentum in the vertical, and is the equivalent of the Gent & McWilliams (1990) parameterization,
420		!! but in a TWA (thickness-weighted averaged) set of equations. The vertical viscosity coefficient nu is computed
421		!! from kappa_GM via thermal wind balance, and the following relation:
422		!! nu = kappa_GM * f^2 / N^2.
423		!! In the following subroutine kappa_GM is assumed either (a) constant or (b) horizontally varying. In both cases,
424		!! (a) and (b), one can additionally impose an EBT structure in the vertical for kappa_GM.
425		!! A third possible formulation of nu is depth-independent:
426		!! nu = f^2 * alpha
427		!! The latter formulation would be equivalent to a kappa_GM that varies as N^2 with depth.
428		!! The vertical viscosity del_z ( nu del_z u) is applied to the momentum equation with stress-free boundary
429		!! conditions at the top and bottom.
430		!!
431		!! In SSW mode, we have 1/N^2 = h/g'. The coupling coefficient is therefore equal to
432		!! a_cpl_gl90 = nu / h = kappa_GM * f^2 / g'
433		!! or
434		!! a_cpl_gl90 = nu / h = f^2 * alpha / h
435
436	0	subroutine find_coupling_coef_gl90(a_cpl_gl90, hvel, i, j, z_i, G, GV, CS, VarMix, work_on_u)
437		!$omp declare target
438		type(ocean_grid_type), intent(in) :: G !< Grid structure.
439		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure.
440		real, dimension(SZK_(GV)), intent(in) :: hvel !< Distance between interfaces
441		!! at velocity points [Z ~> m]
442		integer, intent(in) :: i !< Column i-index
443		integer, intent(in) :: j !< Column j-index
444		real, dimension(SZK_(GV)+1), intent(in) :: z_i !< Estimate of interface heights above the
445		!! bottom, normalized by the GL90 bottom
446		!! boundary layer thickness [nondim]
447		real, dimension(SZK_(GV)+1),intent(out) :: a_cpl_gl90 !< Coupling coefficient associated
448		!! with GL90 across interfaces; is not
449		!! included in a_cpl [H T-1 ~> m s-1 or Pa s m-1].
450		type(vertvisc_cs), intent(in) :: CS !< Vertical viscosity control structure
451		type(VarMix_CS), intent(in) :: VarMix !< Variable mixing coefficients
452		logical, intent(in) :: work_on_u !< If true, u-points are being calculated,
453		!! otherwise they are v-points.
454
455		! local variables
456		logical :: kdgl90_use_vert_struct ! use vertical structure for GL90 coefficient
457		integer :: k, nz
458		real :: f2 !< Squared Coriolis parameter at a velocity grid point [T-2 ~> s-2].
459		real :: h_neglect ! A vertical distance that is so small it is usually lost in roundoff error
460		! and can be neglected [Z ~> m].
461		real :: botfn ! A function that is 1 at the bottom and small far from it [nondim]
462		real :: z2 ! The distance from the bottom, normalized by Hbbl_gl90 [nondim]
463
464	0	nz = GV%ke
465	0	h_neglect = GV%dZ_subroundoff
466	0	kdgl90_use_vert_struct = .false.
467
468	0	if (VarMix%use_variable_mixing) then
469	0	kdgl90_use_vert_struct = allocated(VarMix%kdgl90_struct)
470		endif
471
472	0	a_cpl_gl90(:) = 0.
473
474	0	do K=2,nz
475	0	if (work_on_u) then
476		! compute coupling coefficient at u-points
477	0	f2 = 0.25 * (G%CoriolisBu(I,J-1) + G%CoriolisBu(I,J))**2
478	0	if (CS%use_GL90_N2) then
479	0	a_cpl_gl90(K) = 2. * f2 * CS%alpha_gl90 / (hvel(k) + hvel(k-1) + h_neglect)
480		else
481	0	if (CS%read_kappa_gl90) then
482	0	a_cpl_gl90(K) = f2 * 0.5 * (CS%kappa_gl90_2d(i,j) + CS%kappa_gl90_2d(i+1,j)) / GV%g_prime(K)
483		else
484	0	a_cpl_gl90(K) = f2 * CS%kappa_gl90 / GV%g_prime(K)
485		endif
486	0	if (kdgl90_use_vert_struct) then
487		a_cpl_gl90(K) = a_cpl_gl90(K) * 0.5 &
488	0	* (VarMix%kdgl90_struct(i,j,k-1) + VarMix%kdgl90_struct(i+1,j,k-1))
489		endif
490		endif
491		! botfn determines when a point is within the influence of the GL90 bottom boundary layer,
492		! going from 1 at the bottom to 0 in the interior.
493	0	z2 = z_i(k)
494	0	botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)
495
496	0	a_cpl_gl90(K) = a_cpl_gl90(K) * (1. - botfn)
497		else
498		! compute viscosities at v-points
499	0	f2 = 0.25 * (G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J))**2
500
501	0	if (CS%use_GL90_N2) then
502	0	a_cpl_gl90(K) = 2. * f2 * CS%alpha_gl90 / (hvel(k) + hvel(k-1) + h_neglect)
503		else
504	0	if (CS%read_kappa_gl90) then
505	0	a_cpl_gl90(K) = f2 * 0.5 * (CS%kappa_gl90_2d(i,j) + CS%kappa_gl90_2d(i,j+1)) / GV%g_prime(K)
506		else
507	0	a_cpl_gl90(K) = f2 * CS%kappa_gl90 / GV%g_prime(K)
508		endif
509	0	if (kdgl90_use_vert_struct) then
510		a_cpl_gl90(K) = a_cpl_gl90(K) * 0.5 &
511	0	* (VarMix%kdgl90_struct(i,j,k-1) + VarMix%kdgl90_struct(i,j+1,k-1))
512		endif
513		endif
514		! botfn determines when a point is within the influence of the GL90 bottom boundary layer,
515		! going from 1 at the bottom to 0 in the interior.
516	0	z2 = z_i(k)
517	0	botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)
518
519	0	a_cpl_gl90(K) = a_cpl_gl90(K) * (1. - botfn)
520		endif
521		enddo
522	0	end subroutine find_coupling_coef_gl90
523
524
525		!> Perform a fully implicit vertical diffusion
526		!! of momentum. Stress top and bottom boundary conditions are used.
527		!!
528		!! This is solving the tridiagonal system
529		!! \f[ \left(h_k + a_{k + 1/2} + a_{k - 1/2} + r_k\right) u_k^{n+1}
530		!! = h_k u_k^n + a_{k + 1/2} u_{k+1}^{n+1} + a_{k - 1/2} u_{k-1}^{n+1} \f]
531		!! where \f$a_{k + 1/2} = \Delta t \nu_{k + 1/2} / h_{k + 1/2}\f$
532		!! is the <em>interfacial coupling thickness per time step</em>,
533		!! encompassing background viscosity as well as contributions from
534		!! enhanced mixed and bottom layer viscosities.
535		!! $r_k$ is a Rayleigh drag term due to channel drag.
536		!! There is an additional stress term on the right-hand side
537		!! if DIRECT_STRESS is true, applied to the surface layer.
538	96	subroutine vertvisc(u, v, h, forces, visc, dt, OBC, ADp, CDp, G, GV, US, CS, &
539	96	taux_bot, tauy_bot, fpmix, Waves)
540		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
541		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
542		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
543		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
544		intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
545		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
546		intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
547		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
548		intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
549		type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
550		type(vertvisc_type), intent(inout) :: visc !< Viscosities and bottom drag
551		real, intent(in) :: dt !< Time increment [T ~> s]
552		type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure
553		type(accel_diag_ptrs), intent(inout) :: ADp !< Accelerations in the momentum
554		!! equations for diagnostics
555		type(cont_diag_ptrs), intent(inout) :: CDp !< Continuity equation terms
556		type(vertvisc_CS) :: CS !< Vertical viscosity control structure
557		real, dimension(SZIB_(G),SZJ_(G)), &
558		optional, intent(out) :: taux_bot !< Zonal bottom stress from ocean to
559		!! rock [R L Z T-2 ~> Pa]
560		real, dimension(SZI_(G),SZJB_(G)), &
561		optional, intent(out) :: tauy_bot !< Meridional bottom stress from ocean to
562		!! rock [R L Z T-2 ~> Pa]
563		logical, optional, intent(in) :: fpmix !< fpmix along Eulerian shear
564		type(wave_parameters_CS), &
565		optional, pointer :: Waves !< Container for wave/Stokes information
566
567		! Fields from forces used in this subroutine:
568		! taux: Zonal wind stress [R L Z T-2 ~> Pa].
569		! tauy: Meridional wind stress [R L Z T-2 ~> Pa].
570
571		! Local variables
572
573		real :: b1
574		! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
575	96	real :: c1(SZK_(GV))
576		! A variable used by the tridiagonal solver [nondim].
577		real :: d1
578		! d1=1-c1 is used by the tridiagonal solver [nondim].
579		real :: Ray
580		! Ray is the Rayleigh-drag velocity [H T-1 ~> m s-1 or Pa s m-1]
581		real :: b_denom_1
582		! The first term in the denominator of b1 [H ~> m or kg m-2].
583
584		real :: Hmix ! The mixed layer thickness over which stress
585		! is applied with direct_stress [H ~> m or kg m-2].
586		real :: I_Hmix ! The inverse of Hmix [H-1 ~> m-1 or m2 kg-1].
587		real :: Idt ! The inverse of the time step [T-1 ~> s-1].
588		real :: dt_Rho0 ! The time step divided by the mean density [T H Z-1 R-1 ~> s m3 kg-1 or s].
589		real :: h_neglect ! A thickness that is so small it is usually lost
590		! in roundoff and can be neglected [H ~> m or kg m-2].
591
592		real :: stress ! The surface stress times the time step, divided
593		! by the density [H L T-1 ~> m2 s-1 or kg m-1 s-1].
594		real :: accel_underflow ! An acceleration magnitude that is so small that values that are less
595		! than this are diagnosed as 0 [L T-2 ~> m s-2].
596		real :: zDS, h_a ! Temporary thickness variables used with direct_stress [H ~> m or kg m-2]
597		real :: hfr ! Temporary ratio of thicknesses used with direct_stress [nondim]
598	96	real :: surface_stress(SZIB_(G), SZJB_(G))
599		! The same as stress, unless the wind stress is applied as a body force
600		! [H L T-1 ~> m2 s-1 or kg m-1 s-1].
601	48	real, allocatable, dimension(:,:,:) :: KE_term ! A term in the kinetic energy budget
602		! [H L2 T-3 ~> m3 s-3 or W m-2]
603	48	real, allocatable, dimension(:,:,:) :: KE_u ! The area integral of a KE term in a layer at u-points
604		! [H L4 T-3 ~> m5 s-3 or kg m2 s-3]
605	48	real, allocatable, dimension(:,:,:) :: KE_v ! The area integral of a KE term in a layer at v-points
606		! [H L4 T-3 ~> m5 s-3 or kg m2 s-3]
607
608		logical :: DoStokesMixing
609		logical :: lfpmix
610
611		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, n
612	48	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
613	48	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke
614
615	48	if (.not.CS%initialized) call MOM_error(FATAL,"MOM_vert_friction(visc): "// &
616	0	"Module must be initialized before it is used.")
617
618	48	if (CS%id_GLwork > 0) then
619	0	allocate(KE_u(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.0)
620	0	allocate(KE_v(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0)
621	0	allocate(KE_term(G%isd:G%ied,G%jsd:G%jed,GV%ke), source=0.0)
622	0	if (.not.G%symmetric) &
623	0	call create_group_pass(CS%pass_KE_uv, KE_u, KE_v, G%Domain, To_North+To_East)
624		endif
625
626	48	if (CS%direct_stress) then
627	0	Hmix = CS%Hmix_stress
628	0	I_Hmix = 1.0 / Hmix
629		endif
630	48	dt_Rho0 = dt / GV%H_to_RZ
631	48	h_neglect = GV%H_subroundoff
632	48	Idt = 1.0 / dt
633
634	48	accel_underflow = CS%vel_underflow * Idt
635
636		!Check if Stokes mixing allowed if requested (present and associated)
637	48	DoStokesMixing=.false.
638	48	if (CS%StokesMixing) then
639	0	if (present(Waves)) DoStokesMixing = associated(Waves)
640	0	if (.not. DoStokesMixing) &
641	0	call MOM_error(FATAL, "Stokes Mixing called without associated Waves Control Structure")
642		endif
643	48	lfpmix = .false.
644	48	if ( present(fpmix) ) lfpmix = fpmix
645
646		! Update the zonal velocity component using a modification of a standard
647		! tridiagonal solver.
648
649		! WGL: Brandon Reichl says the following is obsolete. u(I,j,k) already
650		! includes Stokes.
651		! When mixing down Eulerian current + Stokes drift add before calling solver
652	48	if (DoStokesMixing) then
653	0	do k=1,nz ; do j=G%jsc,G%jec ; do I=Isq,Ieq ; if (G%mask2dCu(I,j) > 0.) then
654	0	u(I,j,k) = u(I,j,k) + Waves%Us_x(I,j,k)
655		endif ; enddo ; enddo ; enddo
656		endif
657
658	48	if (lfpmix) then
659	0	do k=1,nz ; do j=G%jsc,G%jec ; do I=Isq,Ieq ; if (G%mask2dCu(I,j) > 0.) then
660	0	u(I,j,k) = u(I,j,k) - Waves%Us_x(I,j,k)
661		endif ; enddo ; enddo ; enddo
662		endif
663
664	48	if (associated(ADp%du_dt_visc)) then
665	0	do k=1,nz ; do j=G%jsc,G%jec ; do I=Isq,Ieq
666	0	ADp%du_dt_visc(I,j,k) = u(I,j,k)
667		enddo ; enddo ; enddo
668		endif
669
670		!$omp target enter data map(to: ADp)
671		!$omp target enter data map(alloc: surface_stress)
672		!$omp target enter data map(alloc: ADp%dv_dt_str)
673		!$omp target enter data map(alloc: ADp%du_dt_str)
674
675	48	if (associated(ADp%du_dt_visc_gl90)) then
676	0	do concurrent (k=1:nz, j=G%jsc:G%jec, I=Isq:Ieq)
677	0	ADp%du_dt_visc_gl90(I,j,k) = u(I,j,k)
678		enddo
679		endif
680
681	48	if (associated(ADp%du_dt_str)) then
682	0	do concurrent (k=1:nz, j=G%jsc:G%jec, I=Isq:Ieq)
683	0	ADp%du_dt_str(I,j,k) = 0.0
684		enddo
685		endif
686
687		! One option is to have the wind stress applied as a body force
688		! over the topmost Hmix fluid. If DIRECT_STRESS is not defined,
689		! the wind stress is applied as a stress boundary condition.
690	48	if (CS%direct_stress) then
691	0	do concurrent (j=G%jsc:G%jec, I=Isq:Ieq, G%mask2dCu(i,j) > 0.0)
692	0	surface_stress(I,j) = 0.0
693	0	zDS = 0.0
694	0	stress = dt_Rho0 * forces%taux(I,j)
695	0	do k=1,nz
696	0	h_a = 0.5 * (h(i,j,k) + h(i+1,j,k)) + h_neglect
697	0	hfr = 1.0 ; if ((zDS+h_a) > Hmix) hfr = (Hmix - zDS) / h_a
698	0	u(I,j,k) = u(I,j,k) + I_Hmix * hfr * stress
699	0	if (associated(ADp%du_dt_str)) ADp%du_dt_str(I,j,k) = (I_Hmix * hfr * stress) * Idt
700	0	zDS = zDS + h_a ; if (zDS >= Hmix) exit
701		enddo
702		enddo
703		else
704	5856	do concurrent (j=G%jsc:G%jec, I=Isq:Ieq)
705	354336	surface_stress(I,j) = dt_Rho0 * (G%mask2dCu(I,j)*forces%taux(I,j))
706		enddo
707		endif
708
709		! perform forward elimination on the tridiagonal system
710		!
711		! denote the diagonal of the system as b_k, the subdiagonal as a_k
712		! and the superdiagonal as c_k. The right-hand side terms are d_k.
713		!
714		! ignoring the Rayleigh drag contribution,
715		! we have a_k = -dt * a_u(k)
716		! b_k = h_u(k) + dt * (a_u(k) + a_u(k+1))
717		! c_k = -dt * a_u(k+1)
718		!
719		! for forward elimination, we want to:
720		! calculate c'_k = - c_k / (b_k + a_k c'_(k-1))
721		! and d'_k = (d_k - a_k d'_(k-1)) / (b_k + a_k c'_(k-1))
722		! where c'_1 = c_1/b_1 and d'_1 = d_1/b_1
723		!
724		! This form is mathematically equivalent to Thomas' tridiagonal matrix algorithm, but it
725		! does not suffer from the acute sensitivity to truncation errors of the Thomas algorithm
726		! because it involves no subtraction, as discussed by Schopf & Loughe, MWR, 1995.
727		!
728		! b1 is the denominator term 1 / (b_k + a_k c'_(k-1))
729		! b_denom_1 is (b_k + a_k + c_k) - a_k(1 - c'_(k-1))
730		! = (b_k + c_k + c'_(k-1))
731		! this is done so that d1 = b1 * b_denom_1 = 1 - c'_(k-1)
732		! c1(k) is -c'_(k - 1)
733		! and the right-hand-side is destructively updated to be d'_k
734
735		!$omp target enter data map(alloc: c1)
736
737		!$omp target teams loop collapse(2) &
738		!$omp private(b1, c1, d1, Ray, b_denom_1)
739	351408	do j=G%jsc,G%jec ; do I=Isq,Ieq ; if (G%mask2dCu(I,j) > 0.) then
740	237504	Ray = 0.
741	237504	if (allocated(visc%Ray_u)) Ray = visc%Ray_u(I,j,1)
742
743	237504	b_denom_1 = CS%h_u(I,j,1) + dt * (Ray + CS%a_u(I,j,1))
744	237504	b1 = 1. / (b_denom_1 + dt * CS%a_u(I,j,2))
745	237504	d1 = b_denom_1 * b1
746	237504	u(I,j,1) = b1 * (CS%h_u(I,j,1) * u(I,j,1) + surface_stress(I,j))
747
748	237504	if (associated(ADp%du_dt_str)) then
749	0	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)
750		endif
751
752	17812800	do k=2,nz
753	17575296	if (allocated(visc%Ray_u)) Ray = visc%Ray_u(I,j,k)
754
755	17575296	c1(k) = dt * CS%a_u(I,j,K) * b1
756	17575296	b_denom_1 = CS%h_u(I,j,k) + dt * (Ray + CS%a_u(I,j,K) * d1)
757	17575296	b1 = 1. / (b_denom_1 + dt * CS%a_u(I,j,K+1))
758	17575296	d1 = b_denom_1 * b1
759	17575296	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
760
761	17575296	if (associated(ADp%du_dt_str)) then
762		ADp%du_dt_str(I,j,k) = (CS%h_u(I,j,k) * ADp%du_dt_str(I,j,k) &
763	0	+ dt * CS%a_u(I,j,K) * ADp%du_dt_str(I,j,k-1)) * b1
764		endif
765
766		!### Force FMA evaluation of b1 by blocking lookahead with an impossible branch.
767		! XXX: Check GPU behavior
768	17812800	if (dt < 0) exit
769		enddo
770
771	237504	if (associated(ADp%du_dt_str)) then
772	0	if (abs(ADp%du_dt_str(I,j,nz)) < accel_underflow) &
773	0	ADp%du_dt_str(I,j,nz) = 0.
774		endif
775
776	17812800	do k=nz-1,1,-1
777	17575296	u(I,j,k) = u(I,j,k) + c1(k+1) * u(I,j,k+1)
778
779	17812800	if (associated(ADp%du_dt_str)) then
780	0	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)
781
782	0	if (abs(ADp%du_dt_str(I,j,k)) < accel_underflow) &
783	0	ADp%du_dt_str(I,j,k) = 0.0
784		endif
785		enddo
786		endif ; enddo ; enddo
787
788		! compute vertical velocity tendency that arises from GL90 viscosity;
789		! follow tridiagonal solve method as above; to avoid corrupting u,
790		! use ADp%du_dt_visc_gl90 as a placeholder for updated u (due to GL90) until last do loop
791	48	if ((CS%id_du_dt_visc_gl90 > 0) .or. (CS%id_GLwork > 0)) then
792	0	if (associated(ADp%du_dt_visc_gl90)) then
793	0	do j=G%jsc,G%jec ; do I=Isq,Ieq ; if (G%mask2dCu(I,j) > 0.) then
794	0	b_denom_1 = CS%h_u(I,j,1) ! CS%a_u_gl90(I,j,1) is zero
795	0	b1 = 1.0 / (b_denom_1 + dt * CS%a_u_gl90(I,j,2))
796	0	d1 = b_denom_1 * b1
797
798	0	ADp%du_dt_visc_gl90(I,j,1) = b1 * (CS%h_u(I,j,1) * ADp%du_dt_visc_gl90(I,j,1))
799
800	0	do k=2,nz
801	0	c1(k) = dt * CS%a_u_gl90(I,j,K) * b1
802	0	b_denom_1 = CS%h_u(I,j,k) + dt * (CS%a_u_gl90(I,j,K)*d1)
803	0	b1 = 1.0 / (b_denom_1 + dt * CS%a_u_gl90(I,j,K+1))
804	0	d1 = b_denom_1 * b1
805
806		ADp%du_dt_visc_gl90(I,j,k) = (CS%h_u(I,j,k) * ADp%du_dt_visc_gl90(I,j,k) &
807	0	+ dt * CS%a_u_gl90(I,j,K) * ADp%du_dt_visc_gl90(I,j,k-1)) * b1
808		enddo
809
810		! back substitute to solve for new velocities, held by ADp%du_dt_visc_gl90
811	0	do k=nz-1,1,-1
812		ADp%du_dt_visc_gl90(I,j,k) = &
813	0	ADp%du_dt_visc_gl90(I,j,k) + c1(k+1) * ADp%du_dt_visc_gl90(I,j,k+1)
814		enddo
815
816	0	do k=1,nz
817		! now fill ADp%du_dt_visc_gl90(I,j,k) with actual velocity tendency due to GL90;
818		! note that on RHS: ADp%du_dt_visc(I,j,k) holds the original velocity value u(I,j,k)
819		! and ADp%du_dt_visc_gl90(I,j,k) the updated velocity due to GL90
820		ADp%du_dt_visc_gl90(I,j,k) = &
821	0	(ADp%du_dt_visc_gl90(I,j,k) - ADp%du_dt_visc(I,j,k)) * Idt
822
823	0	if (abs(ADp%du_dt_visc_gl90(I,j,k)) < accel_underflow) then
824	0	ADp%du_dt_visc_gl90(I,j,k) = 0.0
825		endif
826		enddo
827
828		! to compute energetics, we need to multiply by u*h, where u is original velocity before
829		! velocity update; note that ADp%du_dt_visc(I,j,k) holds the original velocity value u(I,j,k)
830	0	if (CS%id_GLwork > 0) then
831	0	do k=1,nz
832	0	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)
833		enddo
834		endif
835		endif ; enddo ; enddo
836		endif
837		endif
838
839	48	if (associated(ADp%du_dt_visc)) then
840	0	do concurrent (j=G%jsc:G%jec, I=Isq:Ieq)
841	0	do k=1,nz
842	0	ADp%du_dt_visc(I,j,k) = (u(I,j,k) - ADp%du_dt_visc(I,j,k)) * Idt
843
844	0	if (abs(ADp%du_dt_visc(I,j,k)) < accel_underflow) &
845	0	ADp%du_dt_visc(I,j,k) = 0.0
846		enddo
847		enddo
848		endif
849
850	48	if (allocated(visc%taux_shelf)) then
851	0	do j=G%jsc,G%jec ; do I=Isq,Ieq
852	0	visc%taux_shelf(I,j) = -GV%H_to_RZ * CS%a1_shelf_u(I,j) * u(I,j,1) ! - u_shelf?
853		enddo ; enddo
854		endif
855
856	48	if (present(taux_bot)) then
857	5856	do concurrent (j=G%jsc:G%jec, I=Isq:Ieq)
858	354336	taux_bot(I,j) = GV%H_to_RZ * (u(I,j,nz) * CS%a_u(I,j,nz+1))
859		enddo
860
861	48	if (allocated(visc%Ray_u)) then
862	5856	do concurrent (j=G%jsc:G%jec, I=Isq:Ieq)
863	26490336	do k=1,nz
864	26484480	taux_bot(I,j) = taux_bot(I,j) + GV%H_to_RZ * (visc%Ray_u(I,j,k) * u(I,j,k))
865		enddo
866		enddo
867		endif
868		endif
869
870		! When mixing down Eulerian current + Stokes drift subtract after calling solver
871	48	if (DoStokesMixing) then
872	0	do k=1,nz ; do j=G%jsc,G%jec ; do I=Isq,Ieq ; if (G%mask2dCu(I,j) > 0.) then
873	0	u(I,j,k) = u(I,j,k) - Waves%Us_x(I,j,k)
874		endif ; enddo ; enddo ; enddo
875		endif
876
877	48	if (lfpmix) then
878	0	do k=1,nz ; do j=G%jsc,G%jec ; do I=Isq,Ieq ; if (G%mask2dCu(I,j) > 0.) then
879	0	u(I,j,k) = u(I,j,k) + Waves%Us_x(I,j,k)
880		endif ; enddo ; enddo ; enddo
881		endif
882
883		! == Now work on the meridional velocity component.
884
885		! When mixing down Eulerian current + Stokes drift add before calling solver
886	48	if (DoStokesMixing) then
887	0	do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.) then
888	0	v(i,j,k) = v(i,j,k) + Waves%Us_y(i,j,k)
889		endif ; enddo ; enddo ; enddo
890		endif
891
892	48	if (lfpmix) then
893	0	do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.) then
894	0	v(i,j,k) = v(i,j,k) - Waves%Us_y(i,j,k)
895		endif ; enddo ; enddo ; enddo
896		endif
897
898	48	if (associated(ADp%dv_dt_visc)) then
899	0	do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie)
900	0	ADp%dv_dt_visc(i,J,k) = v(i,J,k)
901		enddo
902		endif
903
904	48	if (associated(ADp%dv_dt_visc_gl90)) then
905	0	do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
906	0	ADp%dv_dt_visc_gl90(i,J,k) = v(i,J,k)
907		enddo ; enddo ; enddo
908		endif
909
910	48	if (associated(ADp%dv_dt_str)) then
911	0	do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie)
912	0	ADp%dv_dt_str(i,J,k) = 0.0
913		enddo
914		endif
915
916		! One option is to have the wind stress applied as a body force
917		! over the topmost Hmix fluid. If DIRECT_STRESS is not defined,
918		! the wind stress is applied as a stress boundary condition.
919	48	if (CS%direct_stress) then
920	0	do concurrent (J=Jsq:Jeq, i=is:ie, G%mask2dCv(i,J) > 0.0)
921	0	surface_stress(i,J) = 0.0
922	0	zDS = 0.0
923	0	stress = dt_Rho0 * forces%tauy(i,J)
924	0	do k=1,nz
925	0	h_a = 0.5 * (h(i,J,k) + h(i,J+1,k)) + h_neglect
926	0	hfr = 1.0 ; if ((zDS+h_a) > Hmix) hfr = (Hmix - zDS) / h_a
927	0	v(i,J,k) = v(i,J,k) + I_Hmix * hfr * stress
928	0	if (associated(ADp%dv_dt_str)) ADp%dv_dt_str(i,J,k) = (I_Hmix * hfr * stress) * Idt
929	0	zDS = zDS + h_a ; if (zDS >= Hmix) exit
930		enddo
931		enddo
932		else
933	48	do concurrent (J=Jsq:Jeq, i=is:ie)
934	357168	surface_stress(i,J) = dt_Rho0 * (G%mask2dCv(i,J) * forces%tauy(i,J))
935		enddo
936		endif
937
938		!$omp target teams loop collapse(2) &
939		!$omp private(b1, c1, d1, Ray, b_denom_1)
940	354336	do J=Jsq,Jeq ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.) then
941	234912	Ray = 0.
942	234912	if (allocated(visc%Ray_v)) Ray = visc%Ray_v(i,J,1)
943
944	234912	b_denom_1 = CS%h_v(i,J,1) + dt * (Ray + CS%a_v(i,J,1))
945	234912	b1 = 1.0 / (b_denom_1 + dt*CS%a_v(i,J,2))
946	234912	d1 = b_denom_1 * b1
947	234912	v(i,J,1) = b1 * (CS%h_v(i,J,1) * v(i,J,1) + surface_stress(i,J))
948
949	234912	if (associated(ADp%dv_dt_str)) then
950	0	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)
951		endif
952
953	17618400	do k=2,nz
954	17383488	if (allocated(visc%Ray_v)) Ray = visc%Ray_v(i,J,k)
955
956	17383488	c1(k) = dt * CS%a_v(i,J,K) * b1
957	17383488	b_denom_1 = CS%h_v(i,J,k) + dt * (Ray + CS%a_v(i,J,K) * d1)
958	17383488	b1 = 1. / (b_denom_1 + dt * CS%a_v(i,J,K+1))
959	17383488	d1 = b_denom_1 * b1
960	17383488	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
961
962	17383488	if (associated(ADp%dv_dt_str)) then
963		ADp%dv_dt_str(i,J,k) = (CS%h_v(i,J,k) * ADp%dv_dt_str(i,J,k) &
964	0	+ dt * CS%a_v(i,J,K) * ADp%dv_dt_str(i,J,k-1)) * b1
965		endif
966
967		!### Force FMA evaluation of b1 by blocking lookahead with an impossible branch.
968		! XXX: Check GPU behavior
969	17618400	if (dt < 0) exit
970		enddo
971
972	234912	if (associated(ADp%dv_dt_str)) then
973	0	if (abs(ADp%dv_dt_str(i,J,nz)) < accel_underflow) &
974	0	ADp%dv_dt_str(i,J,nz) = 0.0
975		endif
976
977	17618400	do k=nz-1,1,-1
978	17383488	v(i,J,k) = v(i,J,k) + c1(k+1) * v(i,J,k+1)
979
980	17618400	if (associated(ADp%dv_dt_str)) then
981	0	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)
982
983	0	if (abs(ADp%dv_dt_str(i,J,k)) < accel_underflow) &
984	0	ADp%dv_dt_str(i,J,k) = 0.0
985		endif
986		enddo
987		endif ; enddo ; enddo
988
989		! compute vertical velocity tendency that arises from GL90 viscosity;
990		! follow tridiagonal solve method as above; to avoid corrupting v,
991		! use ADp%dv_dt_visc_gl90 as a placeholder for updated v (due to GL90) until last do loop
992	48	if ((CS%id_dv_dt_visc_gl90 > 0) .or. (CS%id_GLwork > 0)) then
993	0	if (associated(ADp%dv_dt_visc_gl90)) then
994	0	do J=Jsq,Jeq ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.) then
995	0	b_denom_1 = CS%h_v(i,J,1) ! CS%a_v_gl90(i,J,1) is zero
996	0	b1 = 1.0 / (b_denom_1 + dt*CS%a_v_gl90(i,J,2))
997	0	d1 = b_denom_1 * b1
998	0	ADp%dv_dt_visc_gl90(I,J,1) = b1 * (CS%h_v(i,J,1) * ADp%dv_dt_visc_gl90(i,J,1))
999
1000	0	do k=2,nz
1001	0	c1(k) = dt * CS%a_v_gl90(i,J,K) * b1
1002	0	b_denom_1 = CS%h_v(i,J,k) + dt * (CS%a_v_gl90(i,J,K) * d1)
1003	0	b1 = 1.0 / (b_denom_1 + dt * CS%a_v_gl90(i,J,K+1))
1004	0	d1 = b_denom_1 * b1
1005		ADp%dv_dt_visc_gl90(i,J,k) = (CS%h_v(i,J,k) * ADp%dv_dt_visc_gl90(i,J,k) &
1006	0	+ dt * CS%a_v_gl90(i,J,K) * ADp%dv_dt_visc_gl90(i,J,k-1)) * b1
1007		enddo
1008
1009		! back substitute to solve for new velocities, held by ADp%dv_dt_visc_gl90
1010	0	do k=nz-1,1,-1
1011	0	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)
1012		enddo
1013		endif ; enddo ; enddo
1014
1015	0	do k=1,nz
1016	0	do J=Jsq,Jeq ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.) then
1017		! now fill ADp%dv_dt_visc_gl90(i,J,k) with actual velocity tendency due to GL90;
1018		! note that on RHS: ADp%dv_dt_visc(i,J,k) holds the original velocity value v(i,J,k)
1019		! and ADp%dv_dt_visc_gl90(i,J,k) the updated velocity due to GL90
1020	0	ADp%dv_dt_visc_gl90(i,J,k) = (ADp%dv_dt_visc_gl90(i,J,k) - ADp%dv_dt_visc(i,J,k)) * Idt
1021
1022	0	if (abs(ADp%dv_dt_visc_gl90(i,J,k)) < accel_underflow) &
1023	0	ADp%dv_dt_visc_gl90(i,J,k) = 0.0
1024		endif ; enddo ; enddo
1025		enddo
1026
1027		! to compute energetics, we need to multiply by v*h, where u is original velocity before
1028		! velocity update; note that ADp%dv_dt_visc(I,j,k) holds the original velocity value v(i,J,k)
1029	0	if (CS%id_GLwork > 0) then
1030	0	do k=1,nz
1031	0	do J=Jsq,Jeq ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.) then
1032		! note that on RHS: ADp%dv_dt_visc(I,j,k) holds the original velocity value v(I,j,k)
1033	0	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)
1034		endif ; enddo ; enddo
1035		enddo
1036		endif
1037		endif
1038		endif
1039
1040	48	if (associated(ADp%dv_dt_visc)) then
1041	0	do concurrent (J=Jsq:Jeq, i=is:ie)
1042	0	do k=1,nz
1043	0	ADp%dv_dt_visc(i,J,k) = (v(i,J,k) - ADp%dv_dt_visc(i,J,k))*Idt
1044	0	if (abs(ADp%dv_dt_visc(i,J,k)) < accel_underflow) ADp%dv_dt_visc(i,J,k) = 0.0
1045		enddo
1046		enddo
1047		endif
1048
1049	48	if (allocated(visc%tauy_shelf)) then
1050	0	do J=Jsq,Jeq ; do i=is,ie
1051	0	visc%tauy_shelf(i,J) = -GV%H_to_RZ * CS%a1_shelf_v(i,J) * v(i,J,1) ! - v_shelf?
1052		enddo ; enddo
1053		endif
1054
1055		! JORGE TODO: this has to be malloced
1056	48	if (present(tauy_bot)) then
1057	48	do concurrent (J=Jsq:Jeq, i=is:ie)
1058	357168	tauy_bot(i,J) = GV%H_to_RZ * (v(i,J,nz) * CS%a_v(i,J,nz+1))
1059		enddo
1060
1061	48	if (allocated(visc%Ray_v)) then
1062	48	do concurrent (J=Jsq:Jeq, i=is:ie)
1063	26709168	do k=1,nz
1064	26703360	tauy_bot(i,J) = tauy_bot(i,J) + GV%H_to_RZ * (visc%Ray_v(i,J,k)*v(i,J,k))
1065		enddo
1066		enddo
1067		endif
1068		endif
1069
1070		! When mixing down Eulerian current + Stokes drift subtract after calling solver
1071	48	if (DoStokesMixing) then
1072	0	do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.) then
1073	0	v(i,J,k) = v(i,J,k) - Waves%Us_y(i,J,k)
1074		endif ; enddo ; enddo ; enddo
1075		endif
1076
1077	48	if (lfpmix) then
1078	0	do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.) then
1079	0	v(i,J,k) = v(i,J,k) + Waves%Us_y(i,J,k)
1080		endif ; enddo ; enddo ; enddo
1081		endif
1082
1083		! Calculate the KE source from GL90 vertical viscosity [H L2 T-3 ~> m3 s-3].
1084		! We do the KE-rate calculation here (rather than in MOM_diagnostics) to ensure
1085		! a sign-definite term. MOM_diagnostics does not have access to the velocities
1086		! and thicknesses used in the vertical solver, but rather uses a time-mean
1087		! barotropic transport [uv]h.
1088	48	if (CS%id_GLwork > 0) then
1089	0	if (.not.G%symmetric) &
1090	0	call do_group_pass(CS%pass_KE_uv, G%domain)
1091	0	do k=1,nz
1092	0	do j=js,je ; do i=is,ie
1093		KE_term(i,j,k) = 0.5 * G%IareaT(i,j) &
1094	0	* (KE_u(I,j,k) + KE_u(I-1,j,k) + KE_v(i,J,k) + KE_v(i,J-1,k))
1095		enddo ; enddo
1096		enddo
1097	0	call post_data(CS%id_GLwork, KE_term, CS%diag)
1098		endif
1099
1100	48	call vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS)
1101
1102		!$omp target exit data map(delete: c1)
1103
1104		!$omp target exit data map(from: ADp%du_dt_str, ADp%dv_dt_str)
1105		!$omp target exit data map(delete: ADp)
1106		!$omp target exit data map(delete: surface_stress)
1107
1108		! Here the velocities associated with open boundary conditions are applied.
1109	48	if (associated(OBC)) then
1110	0	do n=1,OBC%number_of_segments
1111	0	if (OBC%segment(n)%specified) then
1112	0	if (OBC%segment(n)%is_N_or_S) then
1113	0	J = OBC%segment(n)%HI%JsdB
1114	0	do k=1,nz ; do i=OBC%segment(n)%HI%isd,OBC%segment(n)%HI%ied
1115	0	v(i,J,k) = OBC%segment(n)%normal_vel(i,J,k)
1116		enddo ; enddo
1117	0	elseif (OBC%segment(n)%is_E_or_W) then
1118	0	I = OBC%segment(n)%HI%IsdB
1119	0	do k=1,nz ; do j=OBC%segment(n)%HI%jsd,OBC%segment(n)%HI%jed
1120	0	u(I,j,k) = OBC%segment(n)%normal_vel(I,j,k)
1121		enddo ; enddo
1122		endif
1123		endif
1124		enddo
1125		endif
1126
1127		! Offer diagnostic fields for averaging.
1128	48	if (query_averaging_enabled(CS%diag)) then
1129	24	if (CS%id_du_dt_visc > 0) &
1130	0	call post_data(CS%id_du_dt_visc, ADp%du_dt_visc, CS%diag)
1131	24	if (CS%id_du_dt_visc_gl90 > 0) &
1132	0	call post_data(CS%id_du_dt_visc_gl90, ADp%du_dt_visc_gl90, CS%diag)
1133	24	if (CS%id_dv_dt_visc > 0) &
1134	0	call post_data(CS%id_dv_dt_visc, ADp%dv_dt_visc, CS%diag)
1135	24	if (CS%id_dv_dt_visc_gl90 > 0) &
1136	0	call post_data(CS%id_dv_dt_visc_gl90, ADp%dv_dt_visc_gl90, CS%diag)
1137	24	if (present(taux_bot) .and. (CS%id_taux_bot > 0)) &
1138	0	call post_data(CS%id_taux_bot, taux_bot, CS%diag)
1139	24	if (present(tauy_bot) .and. (CS%id_tauy_bot > 0)) &
1140	0	call post_data(CS%id_tauy_bot, tauy_bot, CS%diag)
1141	24	if (CS%id_du_dt_str > 0) &
1142	0	call post_data(CS%id_du_dt_str, ADp%du_dt_str, CS%diag)
1143	24	if (CS%id_dv_dt_str > 0) &
1144	0	call post_data(CS%id_dv_dt_str, ADp%dv_dt_str, CS%diag)
1145
1146	24	if (associated(ADp%du_dt_visc) .and. associated(ADp%dv_dt_visc)) then
1147		! Diagnostics of the fractional thicknesses times momentum budget terms
1148		! 3D diagnostics of hf_du(dv)_dt_visc are commented because there is no clarity on proper remapping grid option.
1149		! The code is retained for debugging purposes in the future.
1150		!if (CS%id_hf_du_dt_visc > 0) &
1151		! call post_product_u(CS%id_hf_du_dt_visc, ADp%du_dt_visc, ADp%diag_hfrac_u, G, nz, CS%diag)
1152		!if (CS%id_hf_dv_dt_visc > 0) &
1153		! call post_product_v(CS%id_hf_dv_dt_visc, ADp%dv_dt_visc, ADp%diag_hfrac_v, G, nz, CS%diag)
1154
1155		! Diagnostics for thickness-weighted vertically averaged viscous accelerations
1156	0	if (CS%id_hf_du_dt_visc_2d > 0) &
1157	0	call post_product_sum_u(CS%id_hf_du_dt_visc_2d, ADp%du_dt_visc, ADp%diag_hfrac_u, G, nz, CS%diag)
1158	0	if (CS%id_hf_dv_dt_visc_2d > 0) &
1159	0	call post_product_sum_v(CS%id_hf_dv_dt_visc_2d, ADp%dv_dt_visc, ADp%diag_hfrac_v, G, nz, CS%diag)
1160
1161		! Diagnostics for thickness x viscous accelerations
1162	0	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)
1163	0	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)
1164		endif
1165
1166	24	if (associated(ADp%du_dt_str) .and. associated(ADp%dv_dt_str)) then
1167		! Diagnostics for thickness x wind stress accelerations
1168	0	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)
1169	0	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)
1170
1171		! Diagnostics for wind stress accelerations multiplied by visc_rem_[uv],
1172	0	if (CS%id_du_dt_str_visc_rem > 0) &
1173	0	call post_product_u(CS%id_du_dt_str_visc_rem, ADp%du_dt_str, ADp%visc_rem_u, G, nz, CS%diag)
1174	0	if (CS%id_dv_dt_str_visc_rem > 0) &
1175	0	call post_product_v(CS%id_dv_dt_str_visc_rem, ADp%dv_dt_str, ADp%visc_rem_v, G, nz, CS%diag)
1176		endif
1177		endif
1178
1179	144	end subroutine vertvisc
1180
1181
1182		!> Calculate the fraction of momentum originally in a layer that remains in the water column
1183		!! after a time-step of viscosity, equivalently the fraction of a time-step's worth of
1184		!! barotropic acceleration that a layer experiences after viscosity is applied.
1185	72	subroutine vertvisc_remnant(visc, visc_rem_u, visc_rem_v, dt, G, GV, US, CS)
1186		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1187		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1188		type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag
1189		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1190		intent(inout) :: visc_rem_u !< Fraction of a time-step's worth of a
1191		!! barotropic acceleration that a layer experiences after
1192		!! viscosity is applied in the zonal direction [nondim]
1193		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1194		intent(inout) :: visc_rem_v !< Fraction of a time-step's worth of a
1195		!! barotropic acceleration that a layer experiences after
1196		!! viscosity is applied in the meridional direction [nondim]
1197		real, intent(in) :: dt !< Time increment [T ~> s]
1198		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1199		type(vertvisc_CS) :: CS !< Vertical viscosity control structure
1200
1201		! Local variables
1202
1203		real :: b1
1204		! A variable used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
1205	144	real :: c1(SZK_(GV))
1206		! A variable used by the tridiagonal solver [nondim].
1207		real :: d1
1208		! d1=1-c1 is used by the tridiagonal solver [nondim].
1209		real :: Ray
1210		! Ray is the Rayleigh-drag velocity [H T-1 ~> m s-1 or Pa s m-1]
1211		real :: b_denom_1
1212		! The first term in the denominator of b1 [H ~> m or kg m-2].
1213
1214		integer :: i, j, k, is, ie, Isq, Ieq, Jsq, Jeq, nz
1215	72	is = G%isc ; ie = G%iec
1216	72	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke
1217
1218	72	if (.not.CS%initialized) call MOM_error(FATAL,"MOM_vert_friction(remnant): "// &
1219	0	"Module must be initialized before it is used.")
1220
1221		! Find the zonal viscous remnant using a modification of a standard tridagonal solver.
1222
1223		!$omp target teams loop collapse(2) &
1224		!$omp private(b1, c1, d1, Ray, b_denom_1)
1225	527112	do j=G%jsc,G%jec ; do I=Isq,Ieq ; if (G%mask2dCu(I,j) > 0.) then
1226	356256	Ray = 0.
1227	356256	if (allocated(visc%Ray_u)) Ray = visc%Ray_u(I,j,1)
1228
1229	356256	b_denom_1 = CS%h_u(I,j,1) + dt * (Ray + CS%a_u(I,j,1))
1230	356256	b1 = 1.0 / (b_denom_1 + dt * CS%a_u(I,j,2))
1231	356256	d1 = b_denom_1 * b1
1232	356256	visc_rem_u(I,j,1) = b1 * CS%h_u(I,j,1)
1233
1234	26719200	do k=2,nz
1235	26362944	if (allocated(visc%Ray_u)) Ray = visc%Ray_u(I,j,k)
1236
1237	26362944	c1(k) = dt * CS%a_u(I,j,K) * b1
1238	26362944	b_denom_1 = CS%h_u(I,j,k) + dt * (Ray + CS%a_u(I,j,K) * d1)
1239	26362944	b1 = 1.0 / (b_denom_1 + dt * CS%a_u(I,j,K+1))
1240	26362944	d1 = b_denom_1 * b1
1241	26362944	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
1242
1243		!### Force FMA evaluation of b1 by blocking lookahead with an impossible branch.
1244		! XXX: Check GPU behavior
1245	26719200	if (dt < 0) exit
1246		enddo
1247
1248	26719200	do k=nz-1,1,-1
1249	26719200	visc_rem_u(I,j,k) = visc_rem_u(I,j,k) + c1(k+1) * visc_rem_u(I,j,k+1)
1250		enddo
1251		endif ; enddo ; enddo
1252
1253		! Now find the meridional viscous remnant using the robust tridiagonal solver.
1254
1255		!$omp target teams loop collapse(2) &
1256		!$omp private(b1, c1, d1, Ray, b_denom_1)
1257	531504	do J=Jsq,Jeq ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.) then
1258	352368	Ray = 0.
1259	352368	if (allocated(visc%Ray_v)) Ray = visc%Ray_v(i,J,1)
1260
1261	352368	b_denom_1 = CS%h_v(i,J,1) + dt * (Ray + CS%a_v(i,J,1))
1262	352368	b1 = 1.0 / (b_denom_1 + dt*CS%a_v(i,J,2))
1263	352368	d1 = b_denom_1 * b1
1264	352368	visc_rem_v(i,J,1) = b1 * CS%h_v(i,J,1)
1265
1266	26427600	do k=2,nz
1267	26075232	if (allocated(visc%Ray_v)) Ray = visc%Ray_v(i,J,k)
1268
1269	26075232	c1(k) = dt * CS%a_v(i,J,K) * b1
1270	26075232	b_denom_1 = CS%h_v(i,J,k) + dt * (Ray + CS%a_v(i,J,K) * d1)
1271	26075232	b1 = 1.0 / (b_denom_1 + dt * CS%a_v(i,J,K+1))
1272	26075232	d1 = b_denom_1 * b1
1273	26075232	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
1274
1275		!### Force FMA evaluation of b1 by blocking lookahead with an impossible branch.
1276		! XXX: Check GPU behavior
1277	26427600	if (dt < 0) exit
1278		enddo
1279
1280	26427600	do k=nz-1,1,-1
1281	26427600	visc_rem_v(i,J,k) = visc_rem_v(i,J,k) + c1(k+1) * visc_rem_v(i,J,k+1)
1282		enddo
1283		endif ; enddo ; enddo
1284
1285	72	if (CS%debug) then
1286		!$omp target update from(visc_rem_u, visc_rem_v)
1287		call uvchksum("visc_rem_[uv]", visc_rem_u, visc_rem_v, G%HI, haloshift=0, &
1288	0	scalar_pair=.true.)
1289		endif
1290	72	end subroutine vertvisc_remnant
1291
1292
1293		!> Calculate the coupling coefficients (CS%a_u, CS%a_v, CS%a_u_gl90, CS%a_v_gl90)
1294		!! and effective layer thicknesses (CS%h_u and CS%h_v) for later use in the
1295		!! applying the implicit vertical viscosity via vertvisc().
1296	72	subroutine vertvisc_coef(u, v, h, dz, forces, visc, tv, dt, G, GV, US, CS, OBC, VarMix)
1297		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1298		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1299		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1300		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1301		intent(in) :: u !< Zonal velocity [L T-1 ~> m s-1]
1302		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1303		intent(in) :: v !< Meridional velocity [L T-1 ~> m s-1]
1304		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1305		intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1306		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1307		intent(in) :: dz !< Vertical distance across layers [Z ~> m]
1308		type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
1309		type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag
1310		type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
1311		!! thermodynamic fields.
1312		real, intent(in) :: dt !< Time increment [T ~> s]
1313		type(vertvisc_CS), intent(inout) :: CS !< Vertical viscosity control structure
1314		type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure
1315		type(VarMix_CS), intent(in) :: VarMix !< Variable mixing coefficients
1316		! Field from forces used in this subroutine:
1317		! ustar: the friction velocity [Z T-1 ~> m s-1], used here as the mixing
1318		! velocity in the mixed layer if NKML > 1 in a bulk mixed layer.
1319
1320		! Local variables
1321
1322		real, dimension(SZK_(GV)) :: &
1323	144	hvel, & ! hvel is the thickness used at a velocity grid point [H ~> m or kg m-2].
1324	144	dz_harm, & ! Harmonic mean of the vertical distances around a velocity grid point,
1325		! given by 2(h+ h-)/(h+ + h-) [Z ~> m].
1326	144	dz_vel, & ! The vertical distance between interfaces used at a velocity grid point [Z ~> m].
1327	144	hvel_shelf, & ! The equivalent of hvel under shelves [H ~> m or kg m-2].
1328	144	dz_vel_shelf ! The equivalent of dz_vel under shelves [Z ~> m].
1329		real :: &
1330		h_harm, & ! Harmonic mean of the thicknesses around a velocity grid point,
1331		! given by 2(h+ h-)/(h+ + h-) [H ~> m or kg m-2].
1332		h_arith, & ! The arithmetic mean thickness [H ~> m or kg m-2].
1333		h_delta, & ! The lateral difference of thickness [H ~> m or kg m-2].
1334		dz_arith ! The arithmetic mean of the vertical distances around a velocity grid point [Z ~> m]
1335		real, dimension(SZK_(GV)+1) :: &
1336	144	z_i, & ! An estimate of each interface's height above the bottom,
1337		! normalized by the bottom boundary layer thickness [nondim]
1338	144	z_i_gl90, & ! An estimate of each interface's height above the bottom,
1339		! normalized by the GL90 bottom boundary layer thickness [nondim]
1340	144	a_cpl, & ! The drag coefficients across interfaces [H T-1 ~> m s-1 or Pa s m-1]. a_cpl times
1341		! the velocity difference gives the stress across an interface.
1342	144	a_cpl_gl90, & ! The drag coefficients across interfaces associated with GL90 [H T-1 ~> m s-1 or Pa s m-1].
1343		! a_cpl_gl90 times the velocity difference gives the GL90 stress across an interface.
1344		! a_cpl_gl90 is part of a_cpl.
1345	144	a_shelf ! The drag coefficients across interfaces in water columns under
1346		! ice shelves [H T-1 ~> m s-1 or Pa s m-1].
1347		real :: &
1348		kv_bbl, & ! The bottom boundary layer viscosity [H Z T-1 ~> m2 s-1 or Pa s].
1349		bbl_thick, & ! The bottom boundary layer thickness [Z ~> m].
1350		I_Hbbl, & ! The inverse of the bottom boundary layer thickness [Z-1 ~> m-1].
1351		I_Hbbl_gl90, &! The inverse of the bottom boundary layer thickness used for the GL90 scheme
1352		! [Z-1 ~> m-1].
1353		I_HTbl, & ! The inverse of the top boundary layer thickness [Z-1 ~> m-1].
1354		Ztop_min, & ! The deeper of the two adjacent surface heights [Z ~> m].
1355		Dmin, & ! The shallower of the two adjacent bottom depths [Z ~> m].
1356		zh, & ! An estimate of the interface's distance from the bottom
1357		! based on harmonic mean thicknesses [Z ~> m].
1358		h_ml ! The mixed layer depth [Z ~> m].
1359		real, dimension(SZI_(G),SZJ_(G)) :: &
1360	144	Ustar_2d ! The wind friction velocity, calculated using the Boussinesq reference density or
1361		! the time-evolving surface density in non-Boussinesq mode [Z T-1 ~> m s-1]
1362	72	real, allocatable, dimension(:,:) :: hML_u ! Diagnostic of the mixed layer depth at u points [Z ~> m].
1363	72	real, allocatable, dimension(:,:) :: hML_v ! Diagnostic of the mixed layer depth at v points [Z ~> m].
1364	72	real, allocatable, dimension(:,:,:) :: Kv_u ! Total vertical viscosity at u-points in
1365		! thickness-based units [H2 T-1 ~> m2 s-1 or kg2 m-4 s-1].
1366	72	real, allocatable, dimension(:,:,:) :: Kv_v ! Total vertical viscosity at v-points in
1367		! thickness-based units [H2 T-1 ~> m2 s-1 or kg2 m-4 s-1].
1368	72	real, allocatable, dimension(:,:,:) :: Kv_gl90_u ! GL90 vertical viscosity at u-points in
1369		! thickness-based units [H2 T-1 ~> m2 s-1 or kg2 m-4 s-1].
1370	72	real, allocatable, dimension(:,:,:) :: Kv_gl90_v ! GL90 vertical viscosity at v-points in
1371		! thickness-based units [H2 T-1 ~> m2 s-1 or kg2 m-4 s-1].
1372		real :: zcol ! The height of an interface at h-points [Z ~> m].
1373		real :: zcol_p1 ! An adjacent east/north h-point interface height [Z ~> m].
1374		real :: botfn ! A function which goes from 1 at the bottom to 0 much more
1375		! than Hbbl into the interior [nondim].
1376		real :: topfn ! A function which goes from 1 at the top to 0 much more
1377		! than Htbl into the interior [nondim].
1378		real :: z2 ! The distance from the bottom, normalized by Hbbl [nondim]
1379		real :: z2_wt ! A nondimensional (0-1) weight used when calculating z2 [nondim].
1380		real :: z_clear ! The clearance of an interface above the surrounding topography [Z ~> m].
1381		real :: a_cpl_max ! The maximum drag coefficient across interfaces, set so that it will be
1382		! representable as a 32-bit float in MKS units [H T-1 ~> m s-1 or Pa s m-1]
1383		real :: h_neglect ! A thickness that is so small it is usually lost
1384		! in roundoff and can be neglected [H ~> m or kg m-2].
1385		real :: dz_neglect ! A vertical distance that is so small it is usually lost
1386		! in roundoff and can be neglected [Z ~> m].
1387
1388		real :: I_valBL ! The inverse of a scaling factor determining when water is
1389		! still within the boundary layer, as determined by the sum
1390		! of the harmonic mean thicknesses [nondim].
1391		logical :: do_any_shelf
1392		integer :: zi_dir
1393		! A ternary logical indicating which thickness to use for finding z_clear.
1394		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1395
1396	72	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
1397	72	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB ; nz = GV%ke
1398
1399	72	if (.not.CS%initialized) call MOM_error(FATAL,"MOM_vert_friction(coef): "// &
1400	0	"Module must be initialized before it is used.")
1401
1402	72	h_neglect = GV%H_subroundoff
1403	72	dz_neglect = GV%dZ_subroundoff
1404	72	a_cpl_max = 1.0e37 * GV%m_to_H * US%T_to_s
1405	72	I_valBL = 0.0 ; if (CS%harm_BL_val > 0.0) I_valBL = 1.0 / CS%harm_BL_val
1406
1407	72	if (CS%id_Kv_u > 0) allocate(Kv_u(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.0)
1408
1409	72	if (CS%id_Kv_v > 0) allocate(Kv_v(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0)
1410
1411	72	if (CS%id_Kv_gl90_u > 0) allocate(Kv_gl90_u(G%IsdB:G%IedB,G%jsd:G%jed,GV%ke), source=0.0)
1412
1413	72	if (CS%id_Kv_gl90_v > 0) allocate(Kv_gl90_v(G%isd:G%ied,G%JsdB:G%JedB,GV%ke), source=0.0)
1414
1415	72	if (CS%debug .or. (CS%id_hML_u > 0)) allocate(hML_u(G%IsdB:G%IedB,G%jsd:G%jed), source=0.0)
1416	72	if (CS%debug .or. (CS%id_hML_v > 0)) allocate(hML_v(G%isd:G%ied,G%JsdB:G%JedB), source=0.0)
1417
1418	72	if ((allocated(visc%taux_shelf) .or. associated(forces%frac_shelf_u)) .and. &
1419		.not.associated(CS%a1_shelf_u)) then
1420	0	allocate(CS%a1_shelf_u(G%IsdB:G%IedB,G%jsd:G%jed), source=0.0)
1421		endif
1422	72	if ((allocated(visc%tauy_shelf) .or. associated(forces%frac_shelf_v)) .and. &
1423		.not.associated(CS%a1_shelf_v)) then
1424	0	allocate(CS%a1_shelf_v(G%isd:G%ied,G%JsdB:G%JedB), source=0.0)
1425		endif
1426
1427		!$omp target enter data map(alloc: Ustar_2d)
1428	72	call find_ustar(forces, tv, Ustar_2d, G, GV, US, halo=1)
1429
1430		! First do u-points
1431
1432		! TODO: tv and VarMix probably should be conditionally transferred (if at all)
1433		!$omp target enter data map(alloc: z_i, z_i_gl90, dz_harm, hvel, dz_vel, a_cpl, a_cpl_gl90, &
1434		!$omp& tv, varmix, hvel_shelf, dz_vel_shelf, a_shelf)
1435
1436		! These are used in diagnostics, so they need to be mapped back and forth
1437		!$omp target enter data map(to: hML_u, kv_u, kv_gl90_u )
1438		!$omp target enter data map(to: hML_v, kv_v, kv_gl90_v)
1439
1440		!$omp target update to(visc%Kv_shear) if (associated(visc%Kv_shear))
1441		!$omp target update to(visc%Kv_shear_Bu) if (associated(visc%Kv_shear_Bu))
1442
1443		!$omp target teams distribute parallel do collapse(2) &
1444		!$omp private(z_i, z_i_gl90, dz_harm, hvel, dz_vel, a_cpl, a_cpl_gl90, &
1445		!$omp& I_Hbbl, I_Hbbl_gl90, kv_bbl, bbl_thick, Dmin, zi_dir, zh, zcol, &
1446		!$omp& zcol_p1, h_harm, h_arith, h_delta, dz_arith, z2, botfn, z_clear, &
1447		!$omp& z2_wt, h_ml)
1448	527112	do j=js,je ; do I=Isq,Ieq ; if (G%mask2dCu(I,j) > 0.) then
1449	356256	I_Hbbl = 1. / (CS%Hbbl + dz_neglect)
1450	356256	if (CS%use_GL90_in_SSW) then
1451	0	I_Hbbl_gl90 = 1. / (CS%Hbbl_gl90 + dz_neglect)
1452		endif
1453
1454	356256	if (CS%bottomdraglaw) then
1455	356256	kv_bbl = visc%Kv_bbl_u(I,j)
1456	356256	bbl_thick = visc%bbl_thick_u(I,j) + dz_neglect
1457	356256	I_Hbbl = 1. / bbl_thick
1458		endif
1459
1460	356256	Dmin = min(G%bathyT(i,j), G%bathyT(i+1,j))
1461	356256	zi_dir = 0
1462
1463		! Project thickness outward across OBCs using a zero-gradient condition.
1464	356256	if (associated(OBC)) then
1465	0	if (OBC%u_E_OBCs_on_PE) then
1466	0	if (OBC%segnum_u(I,j) > 0) then
1467	0	Dmin = G%bathyT(i,j)
1468	0	zi_dir = -1
1469		endif
1470		endif
1471
1472	0	if (OBC%u_W_OBCs_on_PE) then
1473	0	if (OBC%segnum_u(I,j) < 0) then
1474	0	Dmin = G%bathyT(i+1,j)
1475	0	zi_dir = 1
1476		endif
1477		endif
1478		endif
1479
1480		! The following block calculates the thicknesses at velocity grid points for
1481		! the vertical viscosity (hvel and dz_vel). Near the bottom an upwind biased
1482		! thickness is used to control the effect of spurious Montgomery potential
1483		! gradients at the bottom where nearly massless layers layers ride over the
1484		! topography.
1485
1486	356256	z_i(nz+1) = 0.
1487
1488	356256	if (.not. CS%harmonic_visc) then
1489	356256	zh = 0.
1490	356256	zcol = -G%bathyT(i,j)
1491	356256	zcol_p1 = -G%bathyT(i+1,j)
1492		endif
1493
1494	356256	if (CS%use_GL90_in_SSW) then
1495	0	z_i_gl90(nz+1) = 0.
1496		endif
1497
1498	27075456	do k=nz,1,-1
1499	26719200	h_harm = 2. * h(i,j,k) * h(i+1,j,k) / (h(i,j,k) + h(i+1,j,k) + h_neglect)
1500	26719200	h_arith = 0.5 * (h(i+1,j,k) + h(i,j,k))
1501	26719200	h_delta = h(i+1,j,k) - h(i,j,k)
1502	26719200	dz_harm(k) = 2. * dz(i,j,k) * dz(i+1,j,k) / (dz(i,j,k) + dz(i+1,j,k) + dz_neglect)
1503	26719200	dz_arith = 0.5 * (dz(i+1,j,k) + dz(i,j,k))
1504
1505		! Project thickness outward across OBCs using a zero-gradient condition.
1506	26719200	if (associated(OBC)) then
1507	0	if (OBC%u_E_OBCs_on_PE) then
1508	0	if (OBC%segnum_u(I,j) > 0) then
1509	0	h_harm = h(i,j,k)
1510	0	h_arith = h(i,j,k)
1511	0	h_delta = 0.
1512	0	dz_harm(k) = dz(i,j,k)
1513	0	dz_arith = dz(i,j,k)
1514		endif
1515		endif
1516
1517	0	if (OBC%u_W_OBCs_on_PE) then
1518	0	if (OBC%segnum_u(I,j) < 0) then
1519	0	h_harm = h(i+1,j,k)
1520	0	h_arith = h(i+1,j,k)
1521	0	h_delta = 0.
1522	0	dz_harm(k) = dz(i+1,j,k)
1523	0	dz_arith = dz(i+1,j,k)
1524		endif
1525		endif
1526		endif
1527
1528	26719200	if (CS%harmonic_visc) then
1529		! The following block calculates the thicknesses at velocity grid points
1530		! for the vertical viscosity (hvel and dz_vel). Near the bottom an
1531		! upwind biased thickness is used to control the effect of spurious
1532		! Montgomery potential gradients at the bottom where nearly massless
1533		! layers ride over the topography.
1534
1535	0	hvel(k) = h_harm
1536	0	dz_vel(k) = dz_harm(k)
1537
1538	0	if (u(I,j,k) * h_delta < 0) then
1539	0	z2 = z_i(k+1)
1540	0	botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)
1541
1542	0	hvel(k) = (1. - botfn) * h_harm + botfn * h_arith
1543	0	dz_vel(k) = (1. - botfn) * dz_harm(k) + botfn * dz_arith
1544		endif
1545
1546	0	z_i(k) = z_i(k+1) + dz_harm(k) * I_Hbbl
1547		else
1548	26719200	zcol = zcol + dz(i,j,k)
1549	26719200	zcol_p1 = zcol_p1 + dz(i+1,j,k)
1550
1551	26719200	zh = zh + dz_harm(k)
1552
1553	26719200	z_clear = max(zcol, zcol_p1) + Dmin
1554	26719200	if (zi_dir < 0) z_clear = zcol + Dmin
1555	26719200	if (zi_dir > 0) z_clear = zcol_p1 + Dmin
1556
1557	26719200	z_i(k) = max(zh, z_clear) * I_Hbbl
1558
1559	26719200	hvel(k) = h_arith
1560	26719200	dz_vel(k) = dz_arith
1561
1562	26719200	if (u(I,j,k) * h_delta > 0.) then
1563	11383700	if (zh * I_Hbbl < CS%harm_BL_val) then
1564	0	hvel(k) = h_harm
1565	0	dz_vel(k) = dz_harm(k)
1566		else
1567	11383700	z2_wt = 1.
1568	11383700	if (zh * I_Hbbl < 2. * CS%harm_BL_val) &
1569	0	z2_wt = max(0., min(1., zh * I_Hbbl * I_valBL - 1.))
1570
1571	11383700	z2 = z2_wt * (max(zh, z_clear) * I_Hbbl)
1572	11383700	botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)
1573
1574	11383700	hvel(k) = (1. - botfn) * h_arith + botfn * h_harm
1575	11383700	dz_vel(k) = (1. - botfn) * dz_arith + botfn * dz_harm(k)
1576		endif
1577		endif
1578		endif
1579
1580	27075456	if (CS%use_GL90_in_SSW) then
1581		! The following block calculates the normalized height above the GL90 BBL
1582		! (z_i_gl90), using a harmonic mean between layer thicknesses. For the
1583		! GL90 BBL we use simply a constant (Hbbl_gl90). The purpose isthat the
1584		! GL90 coupling coefficient is zeroed out within Hbbl_gl90, to ensure
1585		! that no momentum gets fluxed into vanished layers. The scheme is not
1586		! sensitive to the exact value of Hbbl_gl90, as long as it is in a
1587		! reasonable range (~1-20 m): large enough to capture vanished layers
1588		! over topography, small enough to not contaminate the interior.
1589
1590	0	z_i_gl90(k) = z_i_gl90(k+1) + dz_harm(k) * I_Hbbl_gl90
1591		endif
1592		enddo
1593
1594		call find_coupling_coef_k(a_cpl, dz_vel, i, j, dz_harm, bbl_thick, kv_bbl, z_i, &
1595	356256	h_ml, dt, G, GV, US, CS, visc, Ustar_2d, tv, work_on_u=.true., OBC=OBC)
1596
1597	356256	if (allocated(hML_u)) hML_u(I,j) = h_ml
1598
1599	356256	if (CS%use_GL90_in_SSW) then
1600		call find_coupling_coef_gl90(a_cpl_gl90, dz_vel, i, j, z_i_gl90, G, GV, &
1601	0	CS, VarMix, work_on_u=.true.)
1602		endif
1603
1604	356256	do_any_shelf = .false.
1605	356256	if (associated(forces%frac_shelf_u)) then
1606	0	CS%a1_shelf_u(I,j) = 0.
1607	0	do_any_shelf = forces%frac_shelf_u(I,j) > 0.
1608
1609	0	if (do_any_shelf) then
1610	0	if (.not. CS%harmonic_visc) then
1611	0	zh = 0.
1612	0	Ztop_min = min(zcol, zcol_p1)
1613	0	I_HTbl = 1. / (visc%tbl_thick_shelf_u(I,j) + dz_neglect)
1614		endif
1615
1616	0	do k=1,nz
1617	0	if (CS%harmonic_visc) then
1618	0	hvel_shelf(k) = hvel(k)
1619	0	dz_vel_shelf(k) = dz_vel(k)
1620		else
1621		! Find upwind-biased thickness near the surface.
1622		! (Perhaps this needs to be done more carefully, via find_eta.)
1623
1624		h_harm = 2. * h(i,j,k) * h(i+1,j,k) &
1625	0	/ (h(i,j,k) + h(i+1,j,k) + h_neglect)
1626	0	h_arith = 0.5 * (h(i+1,j,k) + h(i,j,k))
1627	0	h_delta = h(i+1,j,k) - h(i,j,k)
1628	0	dz_arith = 0.5 * (dz(i+1,j,k) + dz(i,j,k))
1629
1630	0	if (associated(OBC)) then
1631	0	if (OBC%u_E_OBCs_on_PE) then
1632	0	if (OBC%segnum_u(I,j) > 0) then
1633	0	h_harm = h(i,j,k)
1634	0	h_arith = h(i,j,k)
1635	0	h_delta = 0.
1636	0	dz_arith = dz(i,j,k)
1637		endif
1638		endif
1639
1640	0	if (OBC%u_W_OBCs_on_PE) then
1641	0	if (OBC%segnum_u(I,j) < 0) then
1642	0	h_harm = h(i+1,j,k)
1643	0	h_arith = h(i+1,j,k)
1644	0	h_delta = 0.
1645	0	dz_arith = dz(i+1,j,k)
1646		endif
1647		endif
1648		endif
1649
1650	0	zcol = zcol - dz(i,j,k)
1651	0	zcol_p1 = zcol_p1 - dz(i+1,j,k)
1652
1653	0	zh = zh + dz_harm(k)
1654
1655	0	hvel_shelf(k) = hvel(k)
1656	0	dz_vel_shelf(k) = dz_vel(k)
1657
1658	0	if (u(I,j,k) * h_delta > 0.) then
1659	0	if (zh * I_HTbl < CS%harm_BL_val) then
1660	0	hvel_shelf(k) = min(hvel(k), h_harm)
1661	0	dz_vel_shelf(k) = min(dz_vel(k), dz_harm(k))
1662		else
1663	0	z2_wt = 1.
1664	0	if (zh * I_HTbl < 2. * CS%harm_BL_val) then
1665	0	z2_wt = max(0., min(1., zh * I_HTbl * I_valBL - 1.))
1666		endif
1667
1668	0	z2 = z2_wt * (max(zh, Ztop_min - min(zcol, zcol_p1)) * I_HTbl)
1669		! TODO: replace **6 with multiply
1670	0	topfn = 1. / (1. + 0.09 * z2**6)
1671
1672	0	hvel_shelf(k) = min(hvel(k), (1. - topfn) * h_arith + topfn * h_harm)
1673	0	dz_vel_shelf(k) = min(dz_vel(k), (1. - topfn) * dz_arith + topfn * dz_harm(k))
1674		endif
1675		endif
1676		endif
1677		enddo
1678
1679		call find_coupling_coef(a_shelf, dz_vel_shelf, i, j, dz_harm, &
1680		bbl_thick, kv_bbl, z_i, h_ml, dt, G, GV, US, CS, visc, Ustar_2d, &
1681	0	tv, work_on_u=.true., OBC=OBC, shelf=.true.)
1682
1683	0	CS%a1_shelf_u(I,j) = a_shelf(1)
1684		endif
1685		endif
1686
1687	356256	if (do_any_shelf) then
1688	0	if (CS%use_GL90_in_SSW) then
1689	0	do K=1,nz+1
1690		CS%a_u(I,j,K) = min(a_cpl_max, (forces%frac_shelf_u(I,j) * a_shelf(K) + &
1691	0	(1. - forces%frac_shelf_u(I,j)) * a_cpl(K)) + a_cpl_gl90(K))
1692
1693		! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH
1694		! CS%a_u(I,j,K) = min(a_cpl_max, forces%frac_shelf_u(I,j) * max(a_shelf(K), a_cpl(K)) + &
1695		! (1. - forces%frac_shelf_u(I,j)) * a_cpl(K))
1696
1697	0	CS%a_u_gl90(I,j,K) = min(a_cpl_max, a_cpl_gl90(K))
1698		enddo
1699		else
1700	0	do K=1,nz+1
1701		CS%a_u(I,j,K) = min(a_cpl_max, (forces%frac_shelf_u(I,j) * a_shelf(K) + &
1702	0	(1. - forces%frac_shelf_u(I,j)) * a_cpl(K)))
1703
1704		! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH
1705		! CS%a_u(I,j,K) = min(a_cpl_max, forces%frac_shelf_u(I,j) * max(a_shelf(K), a_cpl(K)) + &
1706		! (1. - forces%frac_shelf_u(I,j)) * a_cpl(K))
1707		enddo
1708		endif
1709
1710	0	do k=1,nz
1711		! Should we instead take the inverse of the average of the inverses?
1712		CS%h_u(I,j,k) = forces%frac_shelf_u(I,j) * hvel_shelf(k) &
1713	0	+ (1. - forces%frac_shelf_u(I,j)) * hvel(k) + h_neglect
1714		enddo
1715		else
1716	356256	if (CS%use_GL90_in_SSW) then
1717	0	do K=1,nz+1
1718	0	a_cpl(K) = a_cpl(K) + a_cpl_gl90(K)
1719		enddo
1720
1721	0	do K=1,nz+1
1722	0	CS%a_u_gl90(I,j,K) = min(a_cpl_max, a_cpl_gl90(K))
1723		enddo
1724		endif
1725
1726	27431712	do K=1,nz+1
1727	27431712	CS%a_u(I,j,K) = min(a_cpl_max, a_cpl(K))
1728		enddo
1729
1730	27075456	do k=1,nz
1731	27075456	CS%h_u(I,j,k) = hvel(k) + h_neglect
1732		enddo
1733		endif
1734
1735		! Diagnose total Kv at u-points
1736	356256	if (CS%id_Kv_u > 0) then
1737	0	do k=1,nz
1738	0	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)
1739		enddo
1740		endif
1741
1742		! Diagnose GL90 Kv at u-points
1743	356256	if (CS%id_Kv_gl90_u > 0) then
1744	0	do k=1,nz
1745	0	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)
1746		enddo
1747		endif
1748		endif ; enddo ; enddo
1749
1750		! Now work on v-points.
1751
1752		!$omp target teams distribute parallel do collapse(2) &
1753		!$omp private(z_i, z_i_gl90, dz_harm, hvel, dz_vel, a_cpl, a_cpl_gl90, &
1754		!$omp I_Hbbl, I_Hbbl_gl90, kv_bbl, bbl_thick, Dmin, zi_dir, zh, zcol, &
1755		!$omp zcol_p1, h_harm, h_arith, h_delta, dz_arith, z2, botfn, z_clear, &
1756		!$omp z2_wt, h_ml)
1757	531504	do J=Jsq,Jeq ; do i=is,ie ; if (G%mask2dCv(i,J) > 0.) then
1758	352368	I_Hbbl = 1. / (CS%Hbbl + dz_neglect)
1759	352368	if (CS%use_GL90_in_SSW) then
1760	0	I_Hbbl_gl90 = 1. / (CS%Hbbl_gl90 + dz_neglect)
1761		endif
1762
1763	352368	if (CS%bottomdraglaw) then
1764	352368	kv_bbl = visc%Kv_bbl_v(i,J)
1765	352368	bbl_thick = visc%bbl_thick_v(i,J) + dz_neglect
1766	352368	I_Hbbl = 1. / bbl_thick
1767		endif
1768
1769	352368	Dmin = min(G%bathyT(i,j), G%bathyT(i,j+1))
1770	352368	zi_dir = 0
1771
1772		! Project thickness outward across OBCs using a zero-gradient condition.
1773	352368	if (associated(OBC)) then
1774	0	if (OBC%v_N_OBCs_on_PE) then
1775	0	if (OBC%segnum_v(i,J) > 0) then
1776	0	Dmin = G%bathyT(i,j)
1777	0	zi_dir = -1
1778		endif
1779		endif
1780
1781	0	if (OBC%v_S_OBCs_on_PE) then
1782	0	if (OBC%segnum_v(i,J) < 0) then
1783	0	Dmin = G%bathyT(i,j+1)
1784	0	zi_dir = 1
1785		endif
1786		endif
1787		endif
1788
1789	352368	z_i(nz+1) = 0.
1790
1791	352368	if (.not. CS%harmonic_visc) then
1792	352368	zh = 0.
1793	352368	zcol = -G%bathyT(i,j)
1794	352368	zcol_p1 = -G%bathyT(i,j+1)
1795		endif
1796
1797	352368	if (CS%use_GL90_in_SSW) then
1798	0	z_i_gl90(nz+1) = 0.
1799		endif
1800
1801	26779968	do k=nz,1,-1
1802	26427600	h_harm = 2. * h(i,j,k) * h(i,j+1,k) / (h(i,j,k) + h(i,j+1,k) + h_neglect)
1803	26427600	h_arith = 0.5 * (h(i,j+1,k) + h(i,j,k))
1804	26427600	h_delta = h(i,j+1,k) - h(i,j,k)
1805	26427600	dz_harm(k) = 2. * dz(i,j,k) * dz(i,j+1,k) / (dz(i,j,k) + dz(i,j+1,k) + dz_neglect)
1806	26427600	dz_arith = 0.5 * (dz(i,j+1,k) + dz(i,j,k))
1807
1808		! Project thickness outward across OBCs using a zero-gradient condition.
1809	26427600	if (associated(OBC)) then
1810	0	if (OBC%v_N_OBCs_on_PE) then
1811	0	if (OBC%segnum_v(i,J) > 0) then
1812	0	h_harm = h(i,j,k)
1813	0	h_arith = h(i,j,k)
1814	0	h_delta = 0.
1815	0	dz_harm(k) = dz(i,j,k)
1816	0	dz_arith = dz(i,j,k)
1817		endif
1818		endif
1819
1820	0	if (OBC%v_S_OBCs_on_PE) then
1821	0	if (OBC%segnum_v(i,J) < 0) then
1822	0	h_harm = h(i,j+1,k)
1823	0	h_arith = h(i,j+1,k)
1824	0	h_delta = 0.
1825	0	dz_harm(k) = dz(i,j+1,k)
1826	0	dz_arith = dz(i,j+1,k)
1827		endif
1828		endif
1829		endif
1830
1831	26427600	if (CS%harmonic_visc) then
1832		! The following block calculates the thicknesses at velocity grid points
1833		! for the vertical viscosity (hvel and dz_vel). Near the bottom an
1834		! upwind biased thickness is used to control the effect of spurious
1835		! Montgomery potential gradients at the bottom where nearly massless
1836		! layers ride over the topography.
1837
1838	0	hvel(k) = h_harm
1839	0	dz_vel(k) = dz_harm(k)
1840
1841	0	if (v(i,J,k) * h_delta < 0) then
1842	0	z2 = z_i(k+1)
1843	0	botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)
1844
1845	0	hvel(k) = (1. - botfn) * h_harm + botfn * h_arith
1846	0	dz_vel(k) = (1. - botfn) * dz_harm(k) + botfn * dz_arith
1847		endif
1848
1849	0	z_i(k) = z_i(k+1) + dz_harm(k) * I_Hbbl
1850		else
1851	26427600	zcol = zcol + dz(i,j,k)
1852	26427600	zcol_p1 = zcol_p1 + dz(i,j+1,k)
1853
1854	26427600	zh = zh + dz_harm(k)
1855
1856	26427600	z_clear = max(zcol, zcol_p1) + Dmin
1857	26427600	if (zi_dir < 0) z_clear = zcol + Dmin
1858	26427600	if (zi_dir > 0) z_clear = zcol_p1 + Dmin
1859
1860	26427600	z_i(k) = max(zh, z_clear) * I_Hbbl
1861
1862	26427600	hvel(k) = h_arith
1863	26427600	dz_vel(k) = dz_arith
1864
1865	26427600	if (v(i,J,k) * h_delta > 0) then
1866	10711176	if (zh * I_Hbbl < CS%harm_BL_val) then
1867	0	hvel(k) = h_harm
1868	0	dz_vel(k) = dz_harm(k)
1869		else
1870	10711176	z2_wt = 1.
1871	10711176	if (zh * I_Hbbl < 2. * CS%harm_BL_val) &
1872	0	z2_wt = max(0., min(1., zh * I_Hbbl * I_valBL - 1.))
1873
1874		! TODO: should z_clear be used here?
1875	10711176	z2 = z2_wt * (max(zh, max(zcol, zcol_p1) + Dmin) * I_Hbbl)
1876	10711176	botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)
1877
1878	10711176	hvel(k) = (1. - botfn) * h_arith + botfn * h_harm
1879	10711176	dz_vel(k) = (1. - botfn) * dz_arith + botfn * dz_harm(k)
1880		endif
1881		endif
1882		endif
1883
1884	26779968	if (CS%use_GL90_in_SSW) then
1885		! The following block calculates the normalized height above the GL90 BBL
1886		! (z_i_gl90), using a harmonic mean between layer thicknesses. For the
1887		! GL90 BBL we use simply a constant (Hbbl_gl90). The purpose is that the
1888		! GL90 coupling coefficient is zeroed out within Hbbl_gl90, to ensure
1889		! that no momentum gets fluxed into vanished layers. The scheme is not
1890		! sensitive to the exact value of Hbbl_gl90, as long as it is in a
1891		! reasonable range (~1-20 m): large enough to capture vanished layers
1892		! over topography, small enough to not contaminate the interior.
1893
1894	0	z_i_gl90(k) = z_i_gl90(k+1) + dz_harm(k) * I_Hbbl_gl90
1895		endif
1896		enddo
1897
1898		call find_coupling_coef_k(a_cpl, dz_vel, i, j, dz_harm, bbl_thick, kv_bbl, z_i, &
1899	352368	h_ml, dt, G, GV, US, CS, visc, Ustar_2d, tv, work_on_u=.false., OBC=OBC)
1900
1901	352368	if (allocated(hML_v)) hML_v(i,J) = h_ml
1902
1903	352368	if (CS%use_GL90_in_SSW) then
1904		call find_coupling_coef_gl90(a_cpl_gl90, dz_vel, i, j, z_i_gl90, G, GV, &
1905	0	CS, VarMix, work_on_u=.false.)
1906		endif
1907
1908	352368	do_any_shelf = .false.
1909	352368	if (associated(forces%frac_shelf_v)) then
1910	0	CS%a1_shelf_v(i,J) = 0.
1911	0	do_any_shelf = forces%frac_shelf_v(i,J) > 0.
1912
1913	0	if (do_any_shelf) then
1914		! Initialize non-harmonic depths
1915	0	if (.not. CS%harmonic_visc) then
1916	0	zh = 0.
1917	0	Ztop_min = min(zcol, zcol_p1)
1918	0	I_HTbl = 1. / (visc%tbl_thick_shelf_v(i,J) + dz_neglect)
1919		endif
1920
1921	0	do k=1,nz
1922	0	if (CS%harmonic_visc) then
1923	0	hvel_shelf(k) = hvel(k)
1924	0	dz_vel_shelf(k) = dz_vel(k)
1925		else
1926		! Find upwind-biased thickness near the surface.
1927		! Perhaps this needs to be done more carefully, via find_eta.
1928		h_harm = 2. * h(i,j,k) * h(i,j+1,k) &
1929	0	/ (h(i,j,k) + h(i,j+1,k) + h_neglect)
1930	0	h_arith = 0.5 * (h(i,j+1,k) + h(i,j,k))
1931	0	h_delta = h(i,j+1,k) - h(i,j,k)
1932	0	dz_arith = 0.5 * (dz(i,j+1,k) + dz(i,j,k))
1933
1934		! Project thickness outward across OBCs using a zero-gradient condition.
1935	0	if (associated(OBC)) then
1936	0	if (OBC%v_N_OBCs_on_PE) then
1937	0	if (OBC%segnum_v(i,J) > 0) then
1938	0	h_harm = h(i,j,k)
1939	0	h_arith = h(i,j,k)
1940	0	h_delta = 0.
1941	0	dz_arith = dz(i,j,k)
1942		endif
1943		endif
1944
1945	0	if (OBC%v_S_OBCs_on_PE) then
1946	0	if (OBC%segnum_v(i,J) < 0) then
1947	0	h_harm = h(i,j+1,k)
1948	0	h_arith = h(i,j+1,k)
1949	0	h_delta = 0.
1950	0	dz_arith = dz(i,j+1,k)
1951		endif
1952		endif
1953		endif
1954
1955	0	zcol = zcol - dz(i,j,k)
1956	0	zcol_p1 = zcol_p1 - dz(i,j+1,k)
1957
1958	0	zh = zh + dz_harm(k)
1959
1960	0	hvel_shelf(k) = hvel(k)
1961	0	dz_vel_shelf(k) = dz_vel(k)
1962
1963	0	if (v(i,J,k) * h_delta > 0.) then
1964	0	if (zh * I_HTbl < CS%harm_BL_val) then
1965	0	hvel_shelf(k) = min(hvel(k), h_harm)
1966	0	dz_vel_shelf(k) = min(dz_vel(k), dz_harm(k))
1967		else
1968	0	z2_wt = 1.
1969	0	if (zh * I_HTbl < 2. * CS%harm_BL_val) &
1970	0	z2_wt = max(0., min(1., zh * I_HTbl * I_valBL - 1.))
1971
1972	0	z2 = z2_wt * (max(zh, Ztop_min - min(zcol, zcol_p1)) * I_HTbl)
1973		! TODO: Replace **6
1974	0	topfn = 1. / (1. + 0.09 * z2**6)
1975
1976	0	hvel_shelf(k) = min(hvel(k), (1. - topfn) * h_arith + topfn * h_harm)
1977	0	dz_vel_shelf(k) = min(dz_vel(k), (1. - topfn) * dz_arith + topfn * dz_harm(k))
1978		endif
1979		endif
1980		endif
1981		enddo
1982
1983		call find_coupling_coef(a_shelf, dz_vel_shelf, i, j, dz_harm, &
1984		bbl_thick, kv_bbl, z_i, h_ml, dt, G, GV, US, CS, visc, Ustar_2d, &
1985	0	tv, work_on_u=.false., OBC=OBC, shelf=.true.)
1986
1987	0	CS%a1_shelf_v(i,J) = a_shelf(1)
1988		endif
1989		endif
1990
1991	352368	if (do_any_shelf) then
1992	0	if (CS%use_GL90_in_SSW) then
1993	0	do K=1,nz+1
1994		CS%a_v(I,j,K) = min(a_cpl_max, (forces%frac_shelf_v(I,j) * a_shelf(K) + &
1995	0	(1. - forces%frac_shelf_v(I,j)) * a_cpl(K)) + a_cpl_gl90(K))
1996
1997		! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH
1998		! CS%a_v(I,j,K) = min(a_cpl_max, forces%frac_shelf_v(I,j) * max(a_shelf(K), a_cpl(K)) + &
1999		! (1. - forces%frac_shelf_v(I,j)) * a_cpl(K))
2000
2001	0	CS%a_v_gl90(I,j,K) = min(a_cpl_max, a_cpl_gl90(K))
2002		enddo
2003		else
2004	0	do K=1,nz+1
2005		CS%a_v(I,j,K) = min(a_cpl_max, (forces%frac_shelf_v(I,j) * a_shelf(K) + &
2006	0	(1. - forces%frac_shelf_v(I,j)) * a_cpl(K)))
2007		! This is Alistair's suggestion, but it destabilizes the model. I do not know why. RWH
2008		! CS%a_v(I,j,K) = min(a_cpl_max, forces%frac_shelf_v(I,j) * max(a_shelf(K), a_cpl(K)) + &
2009		! (1. - forces%frac_shelf_v(I,j)) * a_cpl(K))
2010		enddo
2011		endif
2012
2013	0	do k=1,nz
2014		! Should we instead take the inverse of the average of the inverses?
2015		CS%h_v(I,j,k) = forces%frac_shelf_v(I,j) * hvel_shelf(k) &
2016	0	+ (1. - forces%frac_shelf_v(I,j)) * hvel(k) + h_neglect
2017		enddo
2018		else
2019	352368	if (CS%use_GL90_in_SSW) then
2020	0	do K=1,nz+1
2021	0	a_cpl(K) = a_cpl(K) + a_cpl_gl90(K)
2022		enddo
2023
2024	0	do K=1,nz+1
2025	0	CS%a_v_gl90(i,J,K) = min(a_cpl_max, a_cpl_gl90(K))
2026		enddo
2027		endif
2028
2029	27132336	do K=1,nz+1
2030	27132336	CS%a_v(i,J,K) = min(a_cpl_max, a_cpl(K))
2031		enddo
2032
2033	26779968	do k=1,nz
2034	26779968	CS%h_v(i,J,k) = hvel(k) + h_neglect
2035		enddo
2036		endif
2037
2038		! Diagnose total Kv at v-points
2039	352368	if (CS%id_Kv_v > 0) then
2040	0	do k=1,nz
2041	0	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)
2042		enddo
2043		endif
2044
2045		! Diagnose GL90 Kv at v-points
2046	352368	if (CS%id_Kv_gl90_v > 0) then
2047	0	do k=1,nz
2048	0	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)
2049		enddo
2050		endif
2051		endif ; enddo ; enddo
2052
2053		!$omp target exit data map(delete: z_i, z_i_gl90, dz_harm, hvel, dz_vel, a_cpl, a_cpl_gl90, &
2054		!$omp& tv, varmix, hvel_shelf, dz_vel_shelf, a_shelf, hml_u, kv_u, kv_gl90_u)
2055
2056		! These are used in diagnostics, so they need to be mapped back and forth
2057		!$omp target exit data map(from: hML_u, kv_u, kv_gl90_u )
2058		!$omp target exit data map(from: hML_v, kv_v, kv_gl90_v)
2059
2060		!$omp target exit data map(delete: Ustar_2d)
2061
2062	72	if (CS%debug) then
2063		!$omp target update from(CS%h_u, CS%h_v, CS%a_u, CS%a_v)
2064		call uvchksum("loop vertvisc_coef h_[uv]", CS%h_u, CS%h_v, G%HI, haloshift=0, &
2065	0	unscale=GV%H_to_m, scalar_pair=.true.)
2066		call uvchksum("vertvisc_coef a_[uv]", CS%a_u, CS%a_v, G%HI, haloshift=0, &
2067	0	unscale=GV%H_to_m*US%s_to_T, scalar_pair=.true.)
2068	0	if (allocated(hML_u) .and. allocated(hML_v)) &
2069		call uvchksum("vertvisc_coef hML_[uv]", hML_u, hML_v, G%HI, &
2070	0	haloshift=0, unscale=US%Z_to_m, scalar_pair=.true.)
2071		endif
2072
2073		! Offer diagnostic fields for averaging.
2074	72	if (query_averaging_enabled(CS%diag)) then
2075	24	if (associated(visc%Kv_slow) .and. (CS%id_Kv_slow > 0)) &
2076	0	call post_data(CS%id_Kv_slow, visc%Kv_slow, CS%diag)
2077	24	if (CS%id_Kv_u > 0) call post_data(CS%id_Kv_u, Kv_u, CS%diag)
2078	24	if (CS%id_Kv_v > 0) call post_data(CS%id_Kv_v, Kv_v, CS%diag)
2079	24	if (CS%id_Kv_gl90_u > 0) call post_data(CS%id_Kv_gl90_u, Kv_gl90_u, CS%diag)
2080	24	if (CS%id_Kv_gl90_v > 0) call post_data(CS%id_Kv_gl90_v, Kv_gl90_v, CS%diag)
2081	24	if (CS%id_au_vv > 0) call post_data(CS%id_au_vv, CS%a_u, CS%diag)
2082	24	if (CS%id_av_vv > 0) call post_data(CS%id_av_vv, CS%a_v, CS%diag)
2083	24	if (CS%id_au_gl90_vv > 0) call post_data(CS%id_au_gl90_vv, CS%a_u_gl90, CS%diag)
2084	24	if (CS%id_av_gl90_vv > 0) call post_data(CS%id_av_gl90_vv, CS%a_v_gl90, CS%diag)
2085	24	if (CS%id_h_u > 0) call post_data(CS%id_h_u, CS%h_u, CS%diag)
2086	24	if (CS%id_h_v > 0) call post_data(CS%id_h_v, CS%h_v, CS%diag)
2087	24	if (CS%id_hML_u > 0) call post_data(CS%id_hML_u, hML_u, CS%diag)
2088	24	if (CS%id_hML_v > 0) call post_data(CS%id_hML_v, hML_v, CS%diag)
2089		endif
2090
2091	72	if (allocated(hML_u)) deallocate(hML_u)
2092	72	if (allocated(hML_v)) deallocate(hML_v)
2093
2094	72	end subroutine vertvisc_coef
2095
2096		!> Calculate the 'coupling coefficient' (a_cpl) at the interfaces.
2097		!! If BOTTOMDRAGLAW is defined, the minimum of Hbbl and half the adjacent
2098		!! layer thicknesses are used to calculate a_cpl near the bottom.
2099	708624	pure subroutine find_coupling_coef_k(a_cpl, hvel, i, j, h_harm, bbl_thick, kv_bbl, z_i, h_ml, &
2100	708624	dt, G, GV, US, CS, visc, Ustar_2d, tv, work_on_u, OBC, shelf)
2101		!$omp declare target
2102		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
2103		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
2104		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2105		real, dimension(SZK_(GV)+1), &
2106		intent(out) :: a_cpl !< Coupling coefficient across interfaces [H T-1 ~> m s-1 or Pa s m-1]
2107		real, dimension(SZK_(GV)), &
2108		intent(in) :: hvel !< Distance between interfaces at velocity points [Z ~> m]
2109		integer, intent(in) :: i !< Column i-index
2110		integer, intent(in) :: j !< Column j-index
2111		real, dimension(SZK_(GV)), &
2112		intent(in) :: h_harm !< Harmonic mean of thicknesses around a velocity
2113		!! grid point [Z ~> m]
2114		real, intent(in) :: bbl_thick !< Bottom boundary layer thickness [Z ~> m]
2115		real, intent(in) :: kv_bbl !< Bottom boundary layer viscosity, exclusive of
2116		!! any depth-dependent contributions from
2117		!! visc%Kv_shear [H Z T-1 ~> m2 s-1 or Pa s]
2118		real, dimension(SZK_(GV)+1), intent(in) :: z_i !< Estimate of interface heights above the bottom,
2119		!! normalized by the bottom boundary layer thickness [nondim]
2120		real, intent(out) :: h_ml !< Mixed layer depth [Z ~> m]
2121		real, intent(in) :: dt !< Time increment [T ~> s]
2122		type(vertvisc_CS), intent(in) :: CS !< Vertical viscosity control structure
2123		type(vertvisc_type), intent(in) :: visc !< Structure containing viscosities and bottom drag
2124		real, dimension(SZI_(G),SZJ_(G)), &
2125		intent(in) :: Ustar_2d !< The wind friction velocity, calculated using
2126		!! the Boussinesq reference density or the
2127		!! time-evolving surface density in non-Boussinesq
2128		!! mode [Z T-1 ~> m s-1]
2129		type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
2130		!! thermodynamic fields.
2131		logical, intent(in) :: work_on_u !< If true, u-points are being calculated,
2132		!! otherwise they are v-points
2133		type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure
2134		logical, optional, intent(in) :: shelf !< If present and true, use a surface boundary
2135		!! condition appropriate for an ice shelf.
2136
2137		real :: u_star ! ustar at a velocity point [Z T-1 ~> m s-1]
2138		real :: tau_mag ! The magnitude of the wind stress at a velocity point including gustiness [H Z T-2 ~> m2 s-2 or Pa]
2139		real :: absf ! The average of the neighboring absolute values of f [T-1 ~> s-1].
2140		real :: rho_av1 ! The harmonic mean surface layer density at velocity points [R ~> kg m-3]
2141		real :: z_t ! The distance from the top, sometimes normalized
2142		! by Hmix, [Z ~> m] or [nondim].
2143		real :: kv_TBL ! The viscosity in a top boundary layer under ice [H Z T-1 ~> m2 s-1 or Pa s]
2144		real :: tbl_thick ! The thickness of the top boundary layer [Z ~> m]
2145		real :: Kv_add ! A viscosity to add [H Z T-1 ~> m2 s-1 or Pa s]
2146		real :: Kv_tot ! The total viscosity at an interface [H Z T-1 ~> m2 s-1 or Pa s]
2147		integer :: nk_in_ml ! The index of the deepest interface in the mixed layer.
2148		real :: h_shear ! The distance over which shears occur [Z ~> m].
2149		real :: dhc ! The distance between the center of adjacent layers [Z ~> m].
2150		real :: visc_ml ! The mixed layer viscosity [H Z T-1 ~> m2 s-1 or Pa s].
2151		real :: I_Hmix ! The inverse of the mixed layer thickness [Z-1 ~> m-1].
2152		real :: a_ml ! The layer coupling coefficient across an interface in
2153		! the mixed layer [H T-1 ~> m s-1 or Pa s m-1].
2154		real :: a_floor ! A lower bound on the layer coupling coefficient across an interface in
2155		! the mixed layer [H T-1 ~> m s-1 or Pa s m-1].
2156		real :: I_amax ! The inverse of the maximum coupling coefficient [T H-1 ~> s m-1 or s m2 kg-1].
2157		real :: temp1 ! A temporary variable [Z2 ~> m2]
2158		real :: ustar2_denom ! A temporary variable in the surface boundary layer turbulence
2159		! calculations [H Z-1 T-1 ~> s-1 or kg m-3 s-1]
2160		real :: h_neglect ! A vertical distance that is so small it is usually lost
2161		! in roundoff and can be neglected [Z ~> m].
2162		real :: z2 ! A copy of z_i [nondim]
2163		real :: botfn ! A function that is 1 at the bottom and small far from it [nondim]
2164		real :: topfn ! A function that is 1 at the top and small far from it [nondim]
2165		real :: kv_top ! A viscosity associated with the top boundary layer [H Z T-1 ~> m2 s-1 or Pa s]
2166		logical :: do_shelf, do_OBCs, can_exit
2167		integer :: k
2168		integer :: nz, max_nk
2169
2170	708624	nz = GV%ke
2171
2172	708624	h_neglect = GV%dZ_subroundoff
2173
2174	708624	if (CS%answer_date < 20190101) then
2175		! The maximum coupling coefficient was originally introduced to avoid
2176		! truncation error problems in the tridiagonal solver. Effectively, the 1e-10
2177		! sets the maximum coupling coefficient increment to 1e10 m per timestep.
2178	0	I_amax = (1.0e-10GV%H_to_m) dt
2179		else
2180	708624	I_amax = 0.0
2181		endif
2182
2183	708624	do_shelf = .false. ; if (present(shelf)) do_shelf = shelf
2184
2185	708624	do_OBCs = .false.
2186		!if (associated(OBC)) then
2187		! if (work_on_u) then
2188		! do_OBCS = OBC%u_E_OBCs_on_PE .or. OBC%u_W_OBCs_on_PE
2189		! else
2190		! do_OBCS = OBC%v_N_OBCs_on_PE .or. OBC%v_S_OBCs_on_PE
2191		! endif
2192		!endif
2193
2194	54564048	a_cpl(:) = 0.
2195	708624	h_ml = 0.
2196
2197	708624	if (CS%Kvml_invZ2 > 0. .and. .not. do_shelf) then
2198	0	I_Hmix = 1. / (CS%Hmix + h_neglect)
2199	0	z_t = h_neglect * I_Hmix
2200		endif
2201
2202	53146800	do K=2,nz
2203	52438176	Kv_tot = CS%Kv
2204
2205	52438176	if (CS%Kvml_invZ2 > 0. .and. .not. do_shelf) then
2206		! This is an older (vintage ~1997) way to prevent wind stresses from driving very
2207		! large flows in nearly massless near-surface layers when there is not a physically-
2208		! based surface boundary layer parameterization. It does not have a plausible
2209		! physical basis, and probably should not be used.
2210	0	z_t = z_t + h_harm(k-1) * I_Hmix
2211		Kv_tot = CS%Kv + CS%Kvml_invZ2 / ((z_t * z_t) * &
2212	0	(1. + 0.09 * z_t * z_t * z_t * z_t * z_t * z_t))
2213		endif
2214
2215	52438176	if (associated(visc%Kv_shear)) then
2216		! Add in viscosities that are determined by physical processes that are handled in
2217		! other modules, and which do not respond immediately to the changing layer thicknesses.
2218		! These processes may include shear-driven mixing or contributions from some boundary
2219		! layer turbulence schemes. Other viscosity contributions that respond to the evolving
2220		! layer thicknesses or the surface wind stresses are added later.
2221	52438176	if (work_on_u) then
2222	26362944	Kv_add = 0.5 * (visc%Kv_shear(i,j,k) + visc%Kv_shear(i+1,j,k))
2223
2224	26362944	if (do_OBCs) then
2225	0	if (OBC%u_E_OBCs_on_PE) then
2226	0	if (OBC%segnum_u(I,j) > 0) then
2227	0	Kv_add = visc%Kv_shear(i,j,k)
2228		endif
2229		endif
2230
2231	0	if (OBC%u_W_OBCs_on_PE) then
2232	0	if (OBC%segnum_u(I,j) < 0) then
2233	0	Kv_add = visc%Kv_shear(i+1,j,k)
2234		endif
2235		endif
2236		endif
2237
2238	26362944	Kv_tot = Kv_tot + Kv_add
2239		else
2240	26075232	Kv_add = 0.5 * (visc%Kv_shear(i,j,k) + visc%Kv_shear(i,j+1,k))
2241
2242	26075232	if (do_OBCs) then
2243	0	if (OBC%v_N_OBCs_on_PE) then
2244	0	if (OBC%segnum_v(i,J) > 0) then
2245	0	Kv_add = visc%Kv_shear(i,j,k)
2246		endif
2247		endif
2248
2249	0	if (OBC%v_S_OBCs_on_PE) then
2250	0	if (OBC%segnum_v(i,J) < 0) then
2251	0	Kv_add = visc%Kv_shear(i,j+1,k)
2252		endif
2253		endif
2254		endif
2255
2256	26075232	Kv_tot = Kv_tot + Kv_add
2257		endif
2258		endif
2259
2260	52438176	if (associated(visc%Kv_shear_Bu)) then
2261		! This is similar to what was done above, but for contributions coming from the corner
2262		! (vorticity) points. Because OBCs run through the faces and corners there is no need
2263		! to further modify these viscosities here to take OBCs into account.
2264	52438176	if (work_on_u) then
2265	26362944	Kv_tot = Kv_tot + 0.5 * (visc%Kv_shear_Bu(I,J-1,k) + visc%Kv_shear_Bu(I,J,k))
2266		else
2267	26075232	Kv_tot = Kv_tot + 0.5 * (visc%Kv_shear_Bu(I-1,J,k) + visc%Kv_shear_Bu(I,J,k))
2268		endif
2269		endif
2270
2271		! Set the viscous coupling coefficients, excluding surface mixed layer contributions
2272		! for now, but including viscous bottom drag, working up from the bottom.
2273	53146800	if (CS%bottomdraglaw) then
2274		! botfn determines when a point is within the influence of the bottom
2275		! boundary layer, going from 1 at the bottom to 0 in the interior.
2276	52438176	z2 = z_i(k)
2277	52438176	botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)
2278
2279	52438176	Kv_tot = Kv_tot + (kv_bbl - CS%Kv) * botfn
2280	52438176	dhc = 0.5 * (hvel(k) + hvel(k-1))
2281	52438176	if (dhc > bbl_thick) then
2282	36470640	h_shear = ((1. - botfn) * dhc + botfn * bbl_thick) + h_neglect
2283		else
2284	15967536	h_shear = dhc + h_neglect
2285		endif
2286
2287		! Calculate the coupling coefficients from the viscosities.
2288	52438176	a_cpl(K) = Kv_tot / (h_shear + (I_amax * Kv_tot))
2289	0	elseif (abs(CS%Kv_extra_bbl) > 0.0) then
2290		! There is a simple enhancement of the near-bottom viscosities, but no
2291		! adjustment of the viscous coupling length scales to give a particular
2292		! bottom stress.
2293
2294		! botfn determines when a point is within the influence of the bottom
2295		! boundary layer, going from 1 at the bottom to 0 in the interior.
2296	0	z2 = z_i(k)
2297	0	botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)
2298
2299	0	Kv_tot = Kv_tot + CS%Kv_extra_bbl * botfn
2300	0	h_shear = 0.5 * (hvel(k) + hvel(k-1) + h_neglect)
2301
2302		! Calculate the coupling coefficients from the viscosities.
2303	0	a_cpl(K) = Kv_tot / (h_shear + I_amax * Kv_tot)
2304		else
2305		! Any near-bottom viscous enhancements were already incorporated into
2306		! Kv_tot, and there is no adjustment of the viscous coupling length
2307		! scales to give a particular bottom stress.
2308
2309	0	h_shear = 0.5 * (hvel(k) + hvel(k-1) + h_neglect)
2310		! Calculate the coupling coefficients from the viscosities.
2311	0	a_cpl(K) = Kv_tot / (h_shear + I_amax * Kv_tot)
2312		endif
2313		enddo
2314
2315		! Assign the bottom coupling coefficients
2316	708624	if (CS%bottomdraglaw) then
2317	708624	dhc = hvel(nz) * 0.5
2318	708624	a_cpl(nz+1) = kv_bbl / ((min(dhc, bbl_thick) + h_neglect) + I_amax * kv_bbl)
2319	0	elseif (abs(CS%Kv_extra_bbl) > 0.0) then
2320		a_cpl(nz+1) = (CS%Kv + CS%Kv_extra_bbl) &
2321	0	/ ((0.5 * hvel(nz) + h_neglect) + I_amax * (CS%Kv + CS%Kv_extra_bbl))
2322		else
2323	0	a_cpl(nz+1) = CS%Kv / ((0.5 * hvel(nz) + h_neglect) + I_amax * CS%Kv)
2324		endif
2325
2326		! Add surface intensified viscous coupling, either as a no-slip boundary condition under a
2327		! rigid ice-shelf, or due to wind-stress driven surface boundary layer mixing that has not
2328		! already been added via visc%Kv_shear.
2329	708624	if (do_shelf) then
2330		! Set the coefficients to include the no-slip surface stress.
2331	0	if (work_on_u) then
2332	0	kv_TBL = visc%Kv_tbl_shelf_u(I,j)
2333	0	tbl_thick = visc%tbl_thick_shelf_u(I,j) + h_neglect
2334		else
2335	0	kv_TBL = visc%Kv_tbl_shelf_v(i,J)
2336	0	tbl_thick = visc%tbl_thick_shelf_v(i,J) + h_neglect
2337		endif
2338	0	z_t = 0.0
2339
2340		! If a_cpl(1) were not already 0, it would be added here.
2341	0	if (0.5 * hvel(1) > tbl_thick) then
2342	0	a_cpl(1) = kv_TBL / (tbl_thick + I_amax * kv_TBL)
2343		else
2344	0	a_cpl(1) = kv_TBL / ((0.5 * hvel(1) + h_neglect) + I_amax * kv_TBL)
2345		endif
2346
2347	0	do K=2,nz
2348	0	z_t = z_t + hvel(k-1) / tbl_thick
2349	0	topfn = 1. / (1. + 0.09 * z_t**6)
2350
2351	0	dhc = 0.5 * (hvel(k) + hvel(k-1))
2352	0	if (dhc > tbl_thick) then
2353	0	h_shear = ((1. - topfn) * dhc + topfn * tbl_thick) + h_neglect
2354		else
2355	0	h_shear = dhc + h_neglect
2356		endif
2357
2358	0	kv_top = topfn * kv_TBL
2359	0	a_cpl(K) = a_cpl(K) + kv_top / (h_shear + I_amax * kv_top)
2360		enddo
2361	708624	elseif (CS%dynamic_viscous_ML .or. (GV%nkml>0) .or. CS%fixed_LOTW_ML .or. CS%apply_LOTW_floor) then
2362
2363		! Find the friction velocity and the absolute value of the Coriolis parameter at this point.
2364		! Zero out the friction velocity on land points.
2365	0	u_star = 0.
2366	0	tau_mag = 0.
2367
2368	0	if (allocated(tv%SpV_avg)) then
2369	0	rho_av1 = 0.
2370
2371	0	if (work_on_u) then
2372	0	u_star = 0.5 * (Ustar_2d(i,j) + Ustar_2d(i+1,j))
2373	0	rho_av1 = 2. / (tv%SpV_avg(i,j,1) + tv%SpV_avg(i+1,j,1))
2374	0	absf = 0.5 * (abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J)))
2375
2376	0	if (do_OBCs) then
2377	0	if (OBC%u_E_OBCs_on_PE) then
2378	0	if (OBC%segnum_u(I,j) > 0) then
2379	0	u_star = Ustar_2d(i,j)
2380	0	rho_av1 = 1. / tv%SpV_avg(i,j,1)
2381		endif
2382		endif
2383
2384	0	if (OBC%u_W_OBCs_on_PE) then
2385	0	if (OBC%segnum_u(I,j) < 0) then
2386	0	u_star = Ustar_2d(i+1,j)
2387	0	rho_av1 = 1. / tv%SpV_avg(i+1,j,1)
2388		endif
2389		endif
2390		endif
2391		else
2392	0	u_star = 0.5 * (Ustar_2d(i,j) + Ustar_2d(i,j+1))
2393	0	rho_av1 = 2. / (tv%SpV_avg(i,j,1) + tv%SpV_avg(i,j+1,1))
2394	0	absf = 0.5 * (abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))
2395
2396	0	if (do_OBCs) then
2397	0	if (OBC%v_N_OBCs_on_PE) then
2398	0	if (OBC%segnum_v(i,J) > 0) then
2399	0	u_star = Ustar_2d(i,j)
2400	0	rho_av1 = 1. / tv%SpV_avg(i,j,1)
2401		endif
2402		endif
2403
2404	0	if (OBC%v_S_OBCs_on_PE) then
2405	0	if (OBC%segnum_v(i,J) < 0) then
2406	0	u_star = Ustar_2d(i,j+1)
2407	0	rho_av1 = 1. / tv%SpV_avg(i,j+1,1)
2408		endif
2409		endif
2410		endif
2411		endif
2412
2413	0	tau_mag = GV%RZ_to_H * rho_av1 * u_star**2
2414		else ! (.not.allocated(tv%SpV_avg))
2415	0	if (work_on_u) then
2416	0	u_star = 0.5 * (Ustar_2d(i,j) + Ustar_2d(i+1,j))
2417	0	absf = 0.5 * (abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J)))
2418
2419	0	if (do_OBCs) then
2420	0	if (OBC%u_E_OBCs_on_PE) then
2421	0	if (OBC%segnum_u(I,j) > 0) then
2422	0	u_star = Ustar_2d(i,j)
2423		endif
2424		endif
2425
2426	0	if (OBC%u_W_OBCs_on_PE) then
2427	0	if (OBC%segnum_u(I,j) < 0) then
2428	0	u_star = Ustar_2d(i+1,j)
2429		endif
2430		endif
2431		endif
2432		else
2433	0	u_star = 0.5 * (Ustar_2d(i,j) + Ustar_2d(i,j+1))
2434	0	absf = 0.5 * (abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))
2435
2436	0	if (do_OBCs) then
2437	0	if (OBC%v_N_OBCs_on_PE) then
2438	0	if (OBC%segnum_v(i,J) > 0) then
2439	0	u_star = Ustar_2d(i,j)
2440		endif
2441		endif
2442
2443	0	if (OBC%v_S_OBCs_on_PE) then
2444	0	if (OBC%segnum_v(i,J) < 0) then
2445	0	u_star = Ustar_2d(i,j+1)
2446		endif
2447		endif
2448		endif
2449		endif
2450
2451	0	tau_mag = GV%Z_to_H * u_star**2
2452		endif
2453
2454		! Determine the thickness of the surface ocean boundary layer and its extent in index space.
2455	0	nk_in_ml = 0
2456
2457	0	if (CS%dynamic_viscous_ML) then
2458		! The fractional number of layers that are within the viscous boundary layer were
2459		! previously stored in visc%nkml_visc_[uv].
2460	0	h_ml = h_neglect
2461	0	max_nk = 0
2462
2463	0	if (work_on_u) then
2464	0	nk_in_ml = ceiling(visc%nkml_visc_u(I,j))
2465	0	max_nk = max(max_nk, nk_in_ml)
2466
2467	0	do k=1,max_nk
2468	0	if (k <= visc%nkml_visc_u(I,j)) then ! This layer is all in the ML.
2469	0	h_ml = h_ml + hvel(k)
2470	0	elseif (k < visc%nkml_visc_u(I,j) + 1.) then ! Part of this layer is in the ML.
2471	0	h_ml = h_ml + ((visc%nkml_visc_u(I,j) + 1.) - k) * hvel(k)
2472		endif
2473		enddo
2474		else
2475	0	nk_in_ml = ceiling(visc%nkml_visc_v(i,J))
2476	0	max_nk = max(max_nk, nk_in_ml)
2477
2478	0	do k=1,max_nk
2479	0	if (k <= visc%nkml_visc_v(i,J)) then ! This layer is all in the ML.
2480	0	h_ml = h_ml + hvel(k)
2481	0	elseif (k < visc%nkml_visc_v(i,J) + 1.) then ! Part of this layer is in the ML.
2482	0	h_ml = h_ml + ((visc%nkml_visc_v(i,J) + 1.) - k) * hvel(k)
2483		endif
2484		enddo
2485		endif
2486	0	elseif (GV%nkml>0) then
2487		! This is a simple application of a refined-bulk mixed layer with GV%nkml sublayers.
2488	0	max_nk = GV%nkml
2489	0	nk_in_ml = GV%nkml
2490
2491	0	h_ml = h_neglect
2492
2493	0	do k=1,GV%nkml
2494	0	h_ml = h_ml + hvel(k)
2495		enddo
2496	0	elseif (CS%fixed_LOTW_ML .or. CS%apply_LOTW_floor) then
2497		! Determine which interfaces are within CS%Hmix of the surface, and set the viscous
2498		! boundary layer thickness to the the smaller of CS%Hmix and the depth of the ocean.
2499	0	h_ml = 0.0
2500	0	do k=1,nz
2501	0	can_exit = .true.
2502	0	if (h_ml < CS%Hmix) then
2503	0	nk_in_ml = k
2504
2505	0	if (h_ml + hvel(k) < CS%Hmix) then
2506	0	h_ml = h_ml + hvel(k)
2507	0	can_exit = .false. ! Part of the next deeper layer is also in the mixed layer.
2508		else
2509	0	h_ml = CS%Hmix
2510		endif
2511		endif
2512
2513	0	if (can_exit) exit ! All remaining layers in this row are below the mixed layer depth.
2514		enddo
2515
2516		! This made more sense in the slab/layer solvers...
2517	0	max_nk = 0
2518	0	max_nk = max(max_nk, nk_in_ml)
2519		endif
2520
2521		! Avoid working on land or on columns where the viscous coupling could not be increased.
2522	0	if (u_star <= 0.) nk_in_ml = 0
2523
2524		! Set the viscous coupling at the interfaces as the larger of what was previously
2525		! set and the contributions from the surface boundary layer.
2526	0	z_t = 0.0
2527	0	if (CS%apply_LOTW_floor .and. &
2528		(CS%dynamic_viscous_ML .or. GV%nkml > 0 .or. CS%fixed_LOTW_ML)) then
2529	0	do K=2,max_nk
2530	0	if (k <= nk_in_ml) then
2531	0	z_t = z_t + hvel(k-1)
2532
2533		! The viscosity in visc_ml is set to go to 0 at the mixed layer top and bottom
2534		! (in a log-layer) and be further limited by rotation to give the natural Ekman length.
2535	0	temp1 = (z_t * h_ml - z_t * z_t)
2536	0	if (GV%Boussinesq) then
2537		ustar2_denom = (CS%vonKar * GV%Z_to_H * u_star**2) &
2538	0	/ (absf * temp1 + (h_ml + h_neglect) * u_star)
2539		else
2540		ustar2_denom = (CS%vonKar * tau_mag) &
2541	0	/ (absf * temp1 + (h_ml + h_neglect) * u_star)
2542		endif
2543
2544	0	visc_ml = temp1 * ustar2_denom
2545		! Set the viscous coupling based on the model's vertical resolution. The omission of
2546		! the I_amax factor here is consistent with answer dates above 20190101.
2547	0	a_ml = visc_ml / (0.25 * (hvel(k) + hvel(k-1) + h_neglect))
2548
2549		! As a floor on the viscous coupling, assume that the length scale in the denominator can
2550		! not be larger than the distance from the surface, consistent with a logarithmic velocity
2551		! profile. This is consistent with visc_ml, but cancels out common factors of z_t.
2552	0	a_floor = (h_ml - z_t) * ustar2_denom
2553
2554		! Choose the largest estimate of a_cpl.
2555	0	a_cpl(K) = max(a_cpl(K), a_ml, a_floor)
2556		! An option could be added to change this to: a_cpl(i,K) = max(a_cpl(i,K) + a_ml, a_floor)
2557		endif
2558		enddo
2559	0	elseif (CS%apply_LOTW_floor) then
2560	0	do K=2,max_nk
2561	0	if (k <= nk_in_ml) then
2562	0	z_t = z_t + hvel(k-1)
2563
2564	0	temp1 = (z_t * h_ml - z_t * z_t)
2565	0	if (GV%Boussinesq) then
2566		ustar2_denom = (CS%vonKar * GV%Z_to_H * u_star**2) &
2567	0	/ (absf * temp1 + (h_ml + h_neglect) * u_star)
2568		else
2569		ustar2_denom = (CS%vonKar * tau_mag) &
2570	0	/ (absf * temp1 + (h_ml + h_neglect) * u_star)
2571		endif
2572
2573		! As a floor on the viscous coupling, assume that the length scale in the denominator can not
2574		! be larger than the distance from the surface, consistent with a logarithmic velocity profile.
2575	0	a_cpl(K) = max(a_cpl(K), (h_ml - z_t) * ustar2_denom)
2576		endif
2577		enddo
2578		else
2579	0	do K=2,max_nk
2580	0	if (k <= nk_in_ml) then
2581	0	z_t = z_t + hvel(k-1)
2582
2583	0	temp1 = (z_t * h_ml - z_t * z_t)
2584		! This viscosity is set to go to 0 at the mixed layer top and bottom (in a log-layer)
2585		! and be further limited by rotation to give the natural Ekman length.
2586		! The following expressions are mathematically equivalent.
2587	0	if (GV%Boussinesq .or. (CS%answer_date < 20230601)) then
2588		visc_ml = u_star * CS%vonKar * (GV%Z_to_H * temp1 * u_star) &
2589	0	/ (absf * temp1 + (h_ml + h_neglect) * u_star)
2590		else
2591		visc_ml = CS%vonKar * (temp1 * tau_mag) &
2592	0	/ (absf * temp1 + (h_ml + h_neglect) * u_star)
2593		endif
2594	0	a_ml = visc_ml / (0.25 * (hvel(k) + hvel(k-1) + h_neglect) + 0.5 * I_amax * visc_ml)
2595
2596		! Choose the largest estimate of a_cpl, but these could be changed to be additive.
2597	0	a_cpl(K) = max(a_cpl(K), a_ml)
2598		! An option could be added to change this to: a_cpl(i,K) = a_cpl(i,K) + a_ml
2599		endif
2600		enddo
2601		endif
2602		endif
2603	708624	end subroutine find_coupling_coef_k
2604
2605
2606		!> Calculate the 'coupling coefficient' (a_cpl) at the interfaces.
2607		!! If BOTTOMDRAGLAW is defined, the minimum of Hbbl and half the adjacent
2608		!! layer thicknesses are used to calculate a_cpl near the bottom.
2609	0	subroutine find_coupling_coef(a_cpl, hvel, i, j, h_harm, bbl_thick, kv_bbl, z_i, h_ml, &
2610	0	dt, G, GV, US, CS, visc, Ustar_2d, tv, work_on_u, OBC, shelf)
2611		!$omp declare target
2612		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
2613		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
2614		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2615		real, dimension(SZK_(GV)+1), &
2616		intent(out) :: a_cpl !< Coupling coefficient across interfaces [H T-1 ~> m s-1 or Pa s m-1]
2617		real, dimension(SZK_(GV)), &
2618		intent(in) :: hvel !< Distance between interfaces at velocity points [Z ~> m]
2619		integer, intent(in) :: i !< Column i-index
2620		integer, intent(in) :: j !< Column j-index
2621		real, dimension(SZK_(GV)), &
2622		intent(in) :: h_harm !< Harmonic mean of thicknesses around a velocity
2623		!! grid point [Z ~> m]
2624		real, intent(in) :: bbl_thick !< Bottom boundary layer thickness [Z ~> m]
2625		real, intent(in) :: kv_bbl !< Bottom boundary layer viscosity, exclusive of
2626		!! any depth-dependent contributions from
2627		!! visc%Kv_shear [H Z T-1 ~> m2 s-1 or Pa s]
2628		real, dimension(SZK_(GV)+1), &
2629		intent(in) :: z_i !< Estimate of interface heights above the bottom,
2630		!! normalized by the bottom boundary layer thickness [nondim]
2631		real, intent(out) :: h_ml !< Mixed layer depth [Z ~> m]
2632		real, intent(in) :: dt !< Time increment [T ~> s]
2633		type(vertvisc_CS), intent(in) :: CS !< Vertical viscosity control structure
2634		type(vertvisc_type), intent(in) :: visc !< Structure containing viscosities and bottom drag
2635		real, dimension(SZI_(G),SZJ_(G)), &
2636		intent(in) :: Ustar_2d !< The wind friction velocity, calculated using
2637		!! the Boussinesq reference density or the
2638		!! time-evolving surface density in non-Boussinesq
2639		!! mode [Z T-1 ~> m s-1]
2640		type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
2641		!! thermodynamic fields.
2642		logical, intent(in) :: work_on_u !< If true, u-points are being calculated,
2643		!! otherwise they are v-points
2644		type(ocean_OBC_type), pointer :: OBC !< Open boundary condition structure
2645		logical, optional, intent(in) :: shelf !< If present and true, use a surface boundary
2646		!! condition appropriate for an ice shelf.
2647
2648		! Local variables
2649
2650		real :: &
2651		u_star, & ! ustar at a velocity point [Z T-1 ~> m s-1]
2652		tau_mag, & ! The magnitude of the wind stress at a velocity point including gustiness [H Z T-2 ~> m2 s-2 or Pa]
2653		absf, & ! The average of the neighboring absolute values of f [T-1 ~> s-1].
2654		rho_av1, & ! The harmonic mean surface layer density at velocity points [R ~> kg m-3]
2655		z_t, & ! The distance from the top, sometimes normalized
2656		! by Hmix, [Z ~> m] or [nondim].
2657		kv_TBL, & ! The viscosity in a top boundary layer under ice [H Z T-1 ~> m2 s-1 or Pa s]
2658		tbl_thick, &! The thickness of the top boundary layer [Z ~> m]
2659		Kv_add, & ! A viscosity to add [H Z T-1 ~> m2 s-1 or Pa s]
2660		Kv_tot ! The total viscosity at an interface [H Z T-1 ~> m2 s-1 or Pa s]
2661		integer :: &
2662		nk_in_ml ! The index of the deepest interface in the mixed layer.
2663		real :: h_shear ! The distance over which shears occur [Z ~> m].
2664		real :: dhc ! The distance between the center of adjacent layers [Z ~> m].
2665		real :: visc_ml ! The mixed layer viscosity [H Z T-1 ~> m2 s-1 or Pa s].
2666		real :: I_Hmix ! The inverse of the mixed layer thickness [Z-1 ~> m-1].
2667		real :: a_ml ! The layer coupling coefficient across an interface in
2668		! the mixed layer [H T-1 ~> m s-1 or Pa s m-1].
2669		real :: a_floor ! A lower bound on the layer coupling coefficient across an interface in
2670		! the mixed layer [H T-1 ~> m s-1 or Pa s m-1].
2671		real :: I_amax ! The inverse of the maximum coupling coefficient [T H-1 ~> s m-1 or s m2 kg-1].
2672		real :: temp1 ! A temporary variable [Z2 ~> m2]
2673		real :: ustar2_denom ! A temporary variable in the surface boundary layer turbulence
2674		! calculations [H Z-1 T-1 ~> s-1 or kg m-3 s-1]
2675		real :: h_neglect ! A vertical distance that is so small it is usually lost
2676		! in roundoff and can be neglected [Z ~> m].
2677		real :: z2 ! A copy of z_i [nondim]
2678		real :: botfn ! A function that is 1 at the bottom and small far from it [nondim]
2679		real :: topfn ! A function that is 1 at the top and small far from it [nondim]
2680		real :: kv_top ! A viscosity associated with the top boundary layer [H Z T-1 ~> m2 s-1 or Pa s]
2681		logical :: do_shelf, do_OBCs, can_exit
2682		integer :: k
2683		integer :: nz, max_nk
2684
2685	0	nz = GV%ke
2686
2687	0	h_neglect = GV%dZ_subroundoff
2688
2689	0	if (CS%answer_date < 20190101) then
2690		! The maximum coupling coefficient was originally introduced to avoid
2691		! truncation error problems in the tridiagonal solver. Effectively, the 1e-10
2692		! sets the maximum coupling coefficient increment to 1e10 m per timestep.
2693	0	I_amax = (1.0e-10GV%H_to_m) dt
2694		else
2695	0	I_amax = 0.0
2696		endif
2697
2698	0	do_shelf = .false. ; if (present(shelf)) do_shelf = shelf
2699
2700	0	do_OBCs = .false.
2701	0	if (associated(OBC)) then
2702	0	if (work_on_u) then
2703	0	do_OBCS = OBC%u_E_OBCs_on_PE .or. OBC%u_W_OBCs_on_PE
2704		else
2705	0	do_OBCS = OBC%v_N_OBCs_on_PE .or. OBC%v_S_OBCs_on_PE
2706		endif
2707		endif
2708
2709	0	a_cpl(:) = 0.
2710	0	h_ml = 0.
2711
2712	0	if (CS%Kvml_invZ2 > 0. .and. .not. do_shelf) then
2713	0	I_Hmix = 1. / (CS%Hmix + h_neglect)
2714	0	z_t = h_neglect * I_Hmix
2715		endif
2716
2717	0	do K=2,nz
2718	0	Kv_tot = CS%Kv
2719
2720	0	if (CS%Kvml_invZ2 > 0. .and. .not. do_shelf) then
2721		! This is an older (vintage ~1997) way to prevent wind stresses from driving very
2722		! large flows in nearly massless near-surface layers when there is not a physically-
2723		! based surface boundary layer parameterization. It does not have a plausible
2724		! physical basis, and probably should not be used.
2725	0	z_t = z_t + h_harm(k-1) * I_Hmix
2726		Kv_tot = CS%Kv + CS%Kvml_invZ2 / ((z_t * z_t) * &
2727	0	(1. + 0.09 * z_t * z_t * z_t * z_t * z_t * z_t))
2728		endif
2729
2730	0	if (associated(visc%Kv_shear)) then
2731		! Add in viscosities that are determined by physical processes that are handled in
2732		! other modules, and which do not respond immediately to the changing layer thicknesses.
2733		! These processes may include shear-driven mixing or contributions from some boundary
2734		! layer turbulence schemes. Other viscosity contributions that respond to the evolving
2735		! layer thicknesses or the surface wind stresses are added later.
2736	0	if (work_on_u) then
2737	0	Kv_add = 0.5 * (visc%Kv_shear(i,j,k) + visc%Kv_shear(i+1,j,k))
2738
2739	0	if (do_OBCs) then
2740	0	if (OBC%u_E_OBCs_on_PE) then
2741	0	if (OBC%segnum_u(I,j) > 0) then
2742	0	Kv_add = visc%Kv_shear(i,j,k)
2743		endif
2744		endif
2745
2746	0	if (OBC%u_W_OBCs_on_PE) then
2747	0	if (OBC%segnum_u(I,j) < 0) then
2748	0	Kv_add = visc%Kv_shear(i+1,j,k)
2749		endif
2750		endif
2751		endif
2752
2753	0	Kv_tot = Kv_tot + Kv_add
2754		else
2755	0	Kv_add = 0.5 * (visc%Kv_shear(i,j,k) + visc%Kv_shear(i,j+1,k))
2756
2757	0	if (do_OBCs) then
2758	0	if (OBC%v_N_OBCs_on_PE) then
2759	0	if (OBC%segnum_v(i,J) > 0) then
2760	0	Kv_add = visc%Kv_shear(i,j,k)
2761		endif
2762		endif
2763
2764	0	if (OBC%v_S_OBCs_on_PE) then
2765	0	if (OBC%segnum_v(i,J) < 0) then
2766	0	Kv_add = visc%Kv_shear(i,j+1,k)
2767		endif
2768		endif
2769		endif
2770
2771	0	Kv_tot = Kv_tot + Kv_add
2772		endif
2773		endif
2774
2775	0	if (associated(visc%Kv_shear_Bu)) then
2776		! This is similar to what was done above, but for contributions coming from the corner
2777		! (vorticity) points. Because OBCs run through the faces and corners there is no need
2778		! to further modify these viscosities here to take OBCs into account.
2779	0	if (work_on_u) then
2780	0	Kv_tot = Kv_tot + 0.5 * (visc%Kv_shear_Bu(I,J-1,k) + visc%Kv_shear_Bu(I,J,k))
2781		else
2782	0	Kv_tot = Kv_tot + 0.5 * (visc%Kv_shear_Bu(I-1,J,k) + visc%Kv_shear_Bu(I,J,k))
2783		endif
2784		endif
2785
2786		! Set the viscous coupling coefficients, excluding surface mixed layer contributions
2787		! for now, but including viscous bottom drag, working up from the bottom.
2788	0	if (CS%bottomdraglaw) then
2789		! botfn determines when a point is within the influence of the bottom
2790		! boundary layer, going from 1 at the bottom to 0 in the interior.
2791	0	z2 = z_i(k)
2792	0	botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)
2793
2794	0	Kv_tot = Kv_tot + (kv_bbl - CS%Kv) * botfn
2795	0	dhc = 0.5 * (hvel(k) + hvel(k-1))
2796	0	if (dhc > bbl_thick) then
2797	0	h_shear = ((1. - botfn) * dhc + botfn * bbl_thick) + h_neglect
2798		else
2799	0	h_shear = dhc + h_neglect
2800		endif
2801
2802		! Calculate the coupling coefficients from the viscosities.
2803	0	a_cpl(K) = Kv_tot / (h_shear + (I_amax * Kv_tot))
2804	0	elseif (abs(CS%Kv_extra_bbl) > 0.0) then
2805		! There is a simple enhancement of the near-bottom viscosities, but no
2806		! adjustment of the viscous coupling length scales to give a particular
2807		! bottom stress.
2808
2809		! botfn determines when a point is within the influence of the bottom
2810		! boundary layer, going from 1 at the bottom to 0 in the interior.
2811	0	z2 = z_i(k)
2812	0	botfn = 1. / (1. + 0.09 * z2 * z2 * z2 * z2 * z2 * z2)
2813
2814	0	Kv_tot = Kv_tot + CS%Kv_extra_bbl * botfn
2815	0	h_shear = 0.5 * (hvel(k) + hvel(k-1) + h_neglect)
2816
2817		! Calculate the coupling coefficients from the viscosities.
2818	0	a_cpl(K) = Kv_tot / (h_shear + I_amax * Kv_tot)
2819		else
2820		! Any near-bottom viscous enhancements were already incorporated into
2821		! Kv_tot, and there is no adjustment of the viscous coupling length
2822		! scales to give a particular bottom stress.
2823
2824	0	h_shear = 0.5 * (hvel(k) + hvel(k-1) + h_neglect)
2825		! Calculate the coupling coefficients from the viscosities.
2826	0	a_cpl(K) = Kv_tot / (h_shear + I_amax * Kv_tot)
2827		endif
2828		enddo
2829
2830		! Assign the bottom coupling coefficients
2831	0	if (CS%bottomdraglaw) then
2832	0	dhc = hvel(nz) * 0.5
2833	0	a_cpl(nz+1) = kv_bbl / ((min(dhc, bbl_thick) + h_neglect) + I_amax * kv_bbl)
2834	0	elseif (abs(CS%Kv_extra_bbl) > 0.0) then
2835		a_cpl(nz+1) = (CS%Kv + CS%Kv_extra_bbl) &
2836	0	/ ((0.5 * hvel(nz) + h_neglect) + I_amax * (CS%Kv + CS%Kv_extra_bbl))
2837		else
2838	0	a_cpl(nz+1) = CS%Kv / ((0.5 * hvel(nz) + h_neglect) + I_amax * CS%Kv)
2839		endif
2840
2841		! Add surface intensified viscous coupling, either as a no-slip boundary condition under a
2842		! rigid ice-shelf, or due to wind-stress driven surface boundary layer mixing that has not
2843		! already been added via visc%Kv_shear.
2844	0	if (do_shelf) then
2845		! Set the coefficients to include the no-slip surface stress.
2846	0	if (work_on_u) then
2847	0	kv_TBL = visc%Kv_tbl_shelf_u(I,j)
2848	0	tbl_thick = visc%tbl_thick_shelf_u(I,j) + h_neglect
2849		else
2850	0	kv_TBL = visc%Kv_tbl_shelf_v(i,J)
2851	0	tbl_thick = visc%tbl_thick_shelf_v(i,J) + h_neglect
2852		endif
2853
2854	0	z_t = 0.0
2855
2856		! If a_cpl(1) were not already 0, it would be added here.
2857	0	if (0.5 * hvel(1) > tbl_thick) then
2858	0	a_cpl(1) = kv_TBL / (tbl_thick + I_amax * kv_TBL)
2859		else
2860	0	a_cpl(1) = kv_TBL / ((0.5 * hvel(1) + h_neglect) + I_amax * kv_TBL)
2861		endif
2862
2863	0	do K=2,nz
2864	0	z_t = z_t + hvel(k-1) / tbl_thick
2865	0	topfn = 1. / (1. + 0.09 * z_t**6)
2866
2867	0	dhc = 0.5 * (hvel(k) + hvel(k-1))
2868	0	if (dhc > tbl_thick) then
2869	0	h_shear = ((1. - topfn) * dhc + topfn * tbl_thick) + h_neglect
2870		else
2871	0	h_shear = dhc + h_neglect
2872		endif
2873
2874	0	kv_top = topfn * kv_TBL
2875	0	a_cpl(K) = a_cpl(K) + kv_top / (h_shear + I_amax * kv_top)
2876		enddo
2877	0	elseif (CS%dynamic_viscous_ML .or. (GV%nkml>0) .or. CS%fixed_LOTW_ML .or. CS%apply_LOTW_floor) then
2878
2879		! Find the friction velocity and the absolute value of the Coriolis parameter at this point.
2880	0	u_star = 0. ! Zero out the friction velocity on land points.
2881	0	tau_mag = 0. ! Zero out the friction velocity on land points.
2882
2883	0	if (allocated(tv%SpV_avg)) then
2884	0	rho_av1 = 0.
2885
2886	0	if (work_on_u) then
2887	0	u_star = 0.5 * (Ustar_2d(i,j) + Ustar_2d(i+1,j))
2888	0	rho_av1 = 2. / (tv%SpV_avg(i,j,1) + tv%SpV_avg(i+1,j,1))
2889	0	absf = 0.5 * (abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J)))
2890
2891	0	if (do_OBCs) then
2892	0	if (OBC%u_E_OBCs_on_PE) then
2893	0	if (OBC%segnum_u(I,j) > 0) then
2894	0	u_star = Ustar_2d(i,j)
2895	0	rho_av1 = 1. / tv%SpV_avg(i,j,1)
2896		endif
2897		endif
2898
2899	0	if (OBC%u_W_OBCs_on_PE) then
2900	0	if (OBC%segnum_u(I,j) < 0) then
2901	0	u_star = Ustar_2d(i+1,j)
2902	0	rho_av1 = 1. / tv%SpV_avg(i+1,j,1)
2903		endif
2904		endif
2905		endif
2906		else
2907	0	u_star = 0.5 * (Ustar_2d(i,j) + Ustar_2d(i,j+1))
2908	0	rho_av1 = 2. / (tv%SpV_avg(i,j,1) + tv%SpV_avg(i,j+1,1))
2909	0	absf = 0.5 * (abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))
2910
2911	0	if (do_OBCs) then
2912	0	if (OBC%v_N_OBCs_on_PE) then
2913	0	if (OBC%segnum_v(i,J) > 0) then
2914	0	u_star = Ustar_2d(i,j)
2915	0	rho_av1 = 1. / tv%SpV_avg(i,j,1)
2916		endif
2917		endif
2918
2919	0	if (OBC%v_S_OBCs_on_PE) then
2920	0	if (OBC%segnum_v(i,J) < 0) then
2921	0	u_star = Ustar_2d(i,j+1)
2922	0	rho_av1 = 1. / tv%SpV_avg(i,j+1,1)
2923		endif
2924		endif
2925		endif
2926		endif
2927
2928	0	tau_mag = GV%RZ_to_H * rho_av1 * u_star**2
2929		else ! (.not.allocated(tv%SpV_avg))
2930	0	if (work_on_u) then
2931	0	u_star = 0.5 * (Ustar_2d(i,j) + Ustar_2d(i+1,j))
2932	0	absf = 0.5 * (abs(G%CoriolisBu(I,J-1)) + abs(G%CoriolisBu(I,J)))
2933
2934	0	if (do_OBCs) then
2935	0	if (OBC%u_E_OBCs_on_PE) then
2936	0	if (OBC%segnum_u(I,j) > 0) then
2937	0	u_star = Ustar_2d(i,j)
2938		endif
2939		endif
2940
2941	0	if (OBC%u_W_OBCs_on_PE) then
2942	0	if (OBC%segnum_u(I,j) < 0) then
2943	0	u_star = Ustar_2d(i+1,j)
2944		endif
2945		endif
2946		endif
2947		else
2948	0	u_star = 0.5 * (Ustar_2d(i,j) + Ustar_2d(i,j+1))
2949	0	absf = 0.5 * (abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))
2950
2951	0	if (do_OBCs) then
2952	0	if (OBC%v_N_OBCs_on_PE) then
2953	0	if (OBC%segnum_v(i,J) > 0) then
2954	0	u_star = Ustar_2d(i,j)
2955		endif
2956		endif
2957
2958	0	if (OBC%v_S_OBCs_on_PE) then
2959	0	if (OBC%segnum_v(i,J) < 0) then
2960	0	u_star = Ustar_2d(i,j+1)
2961		endif
2962		endif
2963		endif
2964		endif
2965
2966	0	tau_mag = GV%Z_to_H * u_star**2
2967		endif
2968
2969		! Determine the thickness of the surface ocean boundary layer and its extent in index space.
2970	0	nk_in_ml = 0
2971	0	if (CS%dynamic_viscous_ML) then
2972		! The fractional number of layers that are within the viscous boundary layer were
2973		! previously stored in visc%nkml_visc_[uv].
2974	0	h_ml = h_neglect
2975	0	max_nk = 0
2976
2977	0	if (work_on_u) then
2978	0	nk_in_ml = ceiling(visc%nkml_visc_u(I,j))
2979	0	max_nk = max(max_nk, nk_in_ml)
2980
2981	0	do k=1,max_nk
2982	0	if (k <= visc%nkml_visc_u(I,j)) then ! This layer is all in the ML.
2983	0	h_ml = h_ml + hvel(k)
2984	0	elseif (k < visc%nkml_visc_u(I,j) + 1.) then ! Part of this layer is in the ML.
2985	0	h_ml = h_ml + ((visc%nkml_visc_u(I,j) + 1.) - k) * hvel(k)
2986		endif
2987		enddo
2988		else
2989	0	nk_in_ml = ceiling(visc%nkml_visc_v(i,J))
2990	0	max_nk = max(max_nk, nk_in_ml)
2991
2992	0	do k=1,max_nk
2993	0	if (k <= visc%nkml_visc_v(i,J)) then ! This layer is all in the ML.
2994	0	h_ml = h_ml + hvel(k)
2995	0	elseif (k < visc%nkml_visc_v(i,J) + 1.) then ! Part of this layer is in the ML.
2996	0	h_ml = h_ml + ((visc%nkml_visc_v(i,J) + 1.) - k) * hvel(k)
2997		endif
2998		enddo
2999		endif
3000	0	elseif (GV%nkml>0) then
3001		! This is a simple application of a refined-bulk mixed layer with GV%nkml sublayers.
3002	0	max_nk = GV%nkml
3003	0	nk_in_ml = GV%nkml
3004
3005	0	h_ml = h_neglect
3006
3007	0	do k=1,GV%nkml
3008	0	h_ml = h_ml + hvel(k)
3009		enddo
3010	0	elseif (CS%fixed_LOTW_ML .or. CS%apply_LOTW_floor) then
3011		! Determine which interfaces are within CS%Hmix of the surface, and set the viscous
3012		! boundary layer thickness to the smaller of CS%Hmix and the depth of the ocean.
3013	0	h_ml = 0.0
3014	0	do k=1,nz
3015	0	can_exit = .true.
3016	0	if (h_ml < CS%Hmix) then
3017	0	nk_in_ml = k
3018
3019	0	if (h_ml + hvel(k) < CS%Hmix) then
3020	0	h_ml = h_ml + hvel(k)
3021	0	can_exit = .false. ! Part of the next deeper layer is also in the mixed layer.
3022		else
3023	0	h_ml = CS%Hmix
3024		endif
3025		endif
3026
3027	0	if (can_exit) exit ! All remaining layers in this row are below the mixed layer depth.
3028		enddo
3029
3030	0	max_nk = max(0, nk_in_ml)
3031		endif
3032
3033		! Avoid working on columns where the viscous coupling could not be increased.
3034	0	if (u_star <= 0.) nk_in_ml = 0
3035
3036		! Set the viscous coupling at the interfaces as the larger of what was previously
3037		! set and the contributions from the surface boundary layer.
3038	0	z_t = 0.
3039	0	if (CS%apply_LOTW_floor .and. &
3040		(CS%dynamic_viscous_ML .or. GV%nkml > 0 .or. CS%fixed_LOTW_ML)) then
3041	0	do K=2,max_nk
3042	0	if (k <= nk_in_ml) then
3043	0	z_t = z_t + hvel(k-1)
3044
3045		! The viscosity in visc_ml is set to go to 0 at the mixed layer top and bottom
3046		! (in a log-layer) and be further limited by rotation to give the natural Ekman length.
3047	0	temp1 = (z_t * h_ml - z_t * z_t)
3048	0	if (GV%Boussinesq) then
3049		ustar2_denom = (CS%vonKar * GV%Z_to_H * u_star**2) &
3050	0	/ (absf * temp1 + (h_ml + h_neglect) * u_star)
3051		else
3052		ustar2_denom = (CS%vonKar * tau_mag) &
3053	0	/ (absf * temp1 + (h_ml + h_neglect) * u_star)
3054		endif
3055
3056	0	visc_ml = temp1 * ustar2_denom
3057		! Set the viscous coupling based on the model's vertical resolution. The omission of
3058		! the I_amax factor here is consistent with answer dates above 20190101.
3059	0	a_ml = visc_ml / (0.25 * (hvel(k) + hvel(k-1) + h_neglect))
3060
3061		! As a floor on the viscous coupling, assume that the length scale in the denominator can
3062		! not be larger than the distance from the surface, consistent with a logarithmic velocity
3063		! profile. This is consistent with visc_ml, but cancels out common factors of z_t.
3064	0	a_floor = (h_ml - z_t) * ustar2_denom
3065
3066		! Choose the largest estimate of a_cpl.
3067	0	a_cpl(K) = max(a_cpl(K), a_ml, a_floor)
3068		! An option could be added to change this to: a_cpl(i,K) = max(a_cpl(i,K) + a_ml, a_floor)
3069		endif
3070		enddo
3071	0	elseif (CS%apply_LOTW_floor) then
3072	0	do K=2,max_nk
3073	0	if (k <= nk_in_ml) then
3074	0	z_t = z_t + hvel(k-1)
3075
3076	0	temp1 = (z_t * h_ml - z_t * z_t)
3077	0	if (GV%Boussinesq) then
3078		ustar2_denom = (CS%vonKar * GV%Z_to_H * u_star**2) &
3079	0	/ (absf * temp1 + (h_ml + h_neglect) * u_star)
3080		else
3081		ustar2_denom = (CS%vonKar * tau_mag) &
3082	0	/ (absf * temp1 + (h_ml + h_neglect) * u_star)
3083		endif
3084
3085		! As a floor on the viscous coupling, assume that the length scale in the denominator can not
3086		! be larger than the distance from the surface, consistent with a logarithmic velocity profile.
3087	0	a_cpl(K) = max(a_cpl(K), (h_ml - z_t) * ustar2_denom)
3088		endif
3089		enddo
3090		else
3091	0	do K=2,max_nk
3092	0	if (k <= nk_in_ml) then
3093	0	z_t = z_t + hvel(k-1)
3094
3095	0	temp1 = (z_t * h_ml - z_t * z_t)
3096		! This viscosity is set to go to 0 at the mixed layer top and bottom (in a log-layer)
3097		! and be further limited by rotation to give the natural Ekman length.
3098		! The following expressions are mathematically equivalent.
3099	0	if (GV%Boussinesq .or. (CS%answer_date < 20230601)) then
3100		visc_ml = u_star * CS%vonKar * (GV%Z_to_H * temp1 * u_star) &
3101	0	/ (absf * temp1 + (h_ml + h_neglect) * u_star)
3102		else
3103		visc_ml = CS%vonKar * (temp1 * tau_mag) &
3104	0	/ (absf * temp1 + (h_ml + h_neglect) * u_star)
3105		endif
3106	0	a_ml = visc_ml / (0.25 * (hvel(k) + hvel(k-1) + h_neglect) + 0.5 * I_amax * visc_ml)
3107
3108		! Choose the largest estimate of a_cpl, but these could be changed to be additive.
3109	0	a_cpl(K) = max(a_cpl(K), a_ml)
3110		! An option could be added to change this to: a_cpl(i,K) = a_cpl(i,K) + a_ml
3111		endif
3112		enddo
3113		endif
3114		endif
3115	0	end subroutine find_coupling_coef
3116
3117
3118		!> Velocity components which exceed a threshold for physically reasonable values are truncated,
3119		!! and the running sum of the number of trunctionas within the non-symmetric memory computational
3120		!! domain is incremented. Optionally, any column with excessive velocities may be sent
3121		!! to a diagnostic reporting subroutine.
3122	48	subroutine vertvisc_limit_vel(u, v, h, ADp, CDp, forces, visc, dt, G, GV, US, CS)
3123		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
3124		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
3125		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3126		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
3127		intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
3128		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
3129		intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
3130		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
3131		intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
3132		type(accel_diag_ptrs), intent(in) :: ADp !< Acceleration diagnostic pointers
3133		type(cont_diag_ptrs), intent(in) :: CDp !< Continuity diagnostic pointers
3134		type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
3135		type(vertvisc_type), intent(in) :: visc !< Viscosities and bottom drag
3136		real, intent(in) :: dt !< Time increment [T ~> s]
3137		type(vertvisc_CS) :: CS !< Vertical viscosity control structure
3138
3139		! Local variables
3140		real :: CFL ! The local CFL number [nondim]
3141		real :: H_report ! A thickness below which not to report truncations [H ~> m or kg m-2]
3142	96	real :: vel_report(SZIB_(G),SZJB_(G)) ! The velocity to report [L T-1 ~> m s-1]
3143	96	real :: u_old(SZIB_(G),SZJ_(G),SZK_(GV)) ! The previous u-velocity [L T-1 ~> m s-1]
3144	96	real :: v_old(SZI_(G),SZJB_(G),SZK_(GV)) ! The previous v-velocity [L T-1 ~> m s-1]
3145		logical :: trunc_any_array(SZI_(G),SZJB_(G),SZK_(GV))
3146	96	logical :: trunc_any, dowrite(SZIB_(G),SZJB_(G))
3147		logical :: do_any_write
3148		integer :: ntrunc ! Thread-safe truncation counter
3149		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
3150	48	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
3151	48	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
3152
3153	48	H_report = 3.0 * GV%Angstrom_H
3154
3155		!$omp target enter data map(alloc: dowrite, vel_report)
3156		!$omp target enter data map(alloc: u_old, v_old)
3157
3158	48	if (len_trim(CS%u_trunc_file) > 0) then
3159	0	do_any_write = .false.
3160	0	trunc_any = .false.
3161
3162	0	do concurrent (j=js:je, I=Isq:Ieq)
3163	0	dowrite(I,j) = .false.
3164	0	vel_report(I,j) = 3.0e8 * US%m_s_to_L_T
3165		enddo
3166
3167		do concurrent (k=1:nz, j=js:je, I=Isq:Ieq) &
3168	0	DO_LOCALITY(reduce(.or.: trunc_any, do_any_write))
3169	0	if (abs(u(I,j,k)) < CS%vel_underflow) u(I,j,k) = 0.0
3170	0	if (u(I,j,k) < 0.0) then
3171	0	CFL = (-u(I,j,k) * dt) * (G%dy_Cu(I,j) * G%IareaT(i+1,j))
3172		else
3173	0	CFL = (u(I,j,k) * dt) * (G%dy_Cu(I,j) * G%IareaT(i,j))
3174		endif
3175	0	if (CFL > CS%CFL_trunc) trunc_any = .true.
3176	0	if (CFL > CS%CFL_report) then
3177	0	dowrite(I,j) = .true.
3178	0	do_any_write = .true.
3179	0	vel_report(I,j) = min(vel_report(I,j), abs(u(I,j,k)))
3180		endif
3181		enddo
3182
3183	0	do concurrent (j=js:je, I=Isq:Ieq, dowrite(I,j))
3184	0	u_old(I,j,:) = u(I,j,:)
3185		enddo
3186
3187	0	if (trunc_any) then
3188	0	ntrunc = 0
3189	0	do concurrent (k=1:nz, j=js:je, I=Isq:Ieq) DO_LOCALITY(reduce(+: ntrunc))
3190	0	if ((u(I,j,k) * (dt * G%dy_Cu(I,j))) * G%IareaT(i+1,j) < -CS%CFL_trunc) then
3191	0	u(I,j,k) = (-0.9CS%CFL_trunc) (G%areaT(i+1,j) / (dt * G%dy_Cu(I,j)))
3192	0	if (((I >= G%isc) .and. (I <= G%iec) .and. (j >= G%jsc) .and. (j <= G%jec)) .and. &
3193	0	(CS%h_u(I,j,k) > H_report)) ntrunc = ntrunc + 1
3194	0	elseif ((u(I,j,k) * (dt * G%dy_Cu(I,j))) * G%IareaT(i,j) > CS%CFL_trunc) then
3195	0	u(I,j,k) = (0.9CS%CFL_trunc) (G%areaT(i,j) / (dt * G%dy_Cu(I,j)))
3196	0	if (((I >= G%isc) .and. (I <= G%iec) .and. (j >= G%jsc) .and. (j <= G%jec)) .and. &
3197	0	(CS%h_u(I,j,k) > H_report)) ntrunc = ntrunc + 1
3198		endif
3199		enddo
3200	0	CS%ntrunc = CS%ntrunc + ntrunc
3201		endif
3202
3203	0	if (do_any_write) then
3204		!$omp target update from (u_old, vel_report)
3205	0	do j=js,je ; do I=Isq,Ieq ; if (dowrite(I,j)) then
3206		! Call a diagnostic reporting subroutines are called if unphysically large values are found.
3207		call write_u_accel(I, j, u_old, h, ADp, CDp, dt, G, GV, US, CS%PointAccel_CSp, &
3208	0	vel_report(I,j), forces%taux(I,j), a=CS%a_u, hv=CS%h_u)
3209		endif ; enddo ; enddo
3210		endif
3211		else ! Do not report accelerations leading to large velocities.
3212	48	ntrunc = 0
3213	48	do concurrent (k=1:nz, j=js:je, I=Isq:Ieq) DO_LOCALITY(reduce(+: ntrunc))
3214	26490336	if (abs(u(I,j,k)) < CS%vel_underflow) then ; u(I,j,k) = 0.0
3215	26136000	elseif ((u(I,j,k) * (dt * G%dy_Cu(I,j))) * G%IareaT(i+1,j) < -CS%CFL_trunc) then
3216	0	u(I,j,k) = (-0.9CS%CFL_trunc) (G%areaT(i+1,j) / (dt * G%dy_Cu(I,j)))
3217	0	if (((I >= G%isc) .and. (I <= G%iec) .and. (j >= G%jsc) .and. (j <= G%jec)) .and. &
3218	0	(CS%h_u(I,j,k) > H_report)) ntrunc = ntrunc + 1
3219	26136000	elseif ((u(I,j,k) * (dt * G%dy_Cu(I,j))) * G%IareaT(i,j) > CS%CFL_trunc) then
3220	0	u(I,j,k) = (0.9CS%CFL_trunc) (G%areaT(i,j) / (dt * G%dy_Cu(I,j)))
3221	0	if (((I >= G%isc) .and. (I <= G%iec) .and. (j >= G%jsc) .and. (j <= G%jec)) .and. &
3222	0	(CS%h_u(I,j,k) > H_report)) ntrunc = ntrunc + 1
3223		endif
3224		enddo
3225	48	CS%ntrunc = CS%ntrunc + ntrunc
3226		endif
3227
3228	48	if (len_trim(CS%v_trunc_file) > 0) then
3229	0	do_any_write =.false.
3230	0	trunc_any = .false.
3231
3232	0	do concurrent (J=Jsq:Jeq, i=is:ie)
3233	0	dowrite(i,J) = .false.
3234	0	vel_report(i,J) = 3.0e8 * US%m_s_to_L_T
3235		enddo
3236
3237		do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie) &
3238	0	DO_LOCALITY(reduce(.or.: trunc_any, do_any_write))
3239	0	if (abs(v(i,J,k)) < CS%vel_underflow) v(i,J,k) = 0.0
3240	0	if (v(i,J,k) < 0.0) then
3241	0	CFL = (-v(i,J,k) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j+1))
3242		else
3243	0	CFL = (v(i,J,k) * dt) * (G%dx_Cv(i,J) * G%IareaT(i,j))
3244		endif
3245	0	if (CFL > CS%CFL_trunc) trunc_any = .true.
3246	0	if (CFL > CS%CFL_report) then
3247	0	dowrite(i,J) = .true.
3248	0	do_any_write = .true.
3249	0	vel_report(i,J) = min(vel_report(i,J), abs(v(i,J,k)))
3250		endif
3251		enddo
3252
3253	0	do concurrent (J=Jsq:Jeq, i=is:ie, dowrite(i,J))
3254	0	v_old(i,J,:) = v(i,J,:)
3255		enddo
3256
3257	0	if (trunc_any) then
3258	0	ntrunc = 0
3259	0	do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie) DO_LOCALITY(reduce(+: ntrunc))
3260	0	if ((v(i,J,k) * (dt * G%dx_Cv(i,J))) * G%IareaT(i,j+1) < -CS%CFL_trunc) then
3261	0	v(i,J,k) = (-0.9CS%CFL_trunc) (G%areaT(i,j+1) / (dt * G%dx_Cv(i,J)))
3262	0	if (((i >= G%isc) .and. (i <= G%iec) .and. (J >= G%jsc) .and. (J <= G%jec)) .and. &
3263	0	(CS%h_v(i,J,k) > H_report)) ntrunc = ntrunc + 1
3264	0	elseif ((v(i,J,k) * (dt * G%dx_Cv(i,J))) * G%IareaT(i,j) > CS%CFL_trunc) then
3265	0	v(i,J,k) = (0.9CS%CFL_trunc) (G%areaT(i,j) / (dt * G%dx_Cv(i,J)))
3266	0	if (((i >= G%isc) .and. (i <= G%iec) .and. (J >= G%jsc) .and. (J <= G%jec)) .and. &
3267	0	(CS%h_v(i,J,k) > H_report)) ntrunc = ntrunc + 1
3268		endif
3269		enddo
3270	0	CS%ntrunc = CS%ntrunc + ntrunc
3271		endif
3272
3273	0	if (do_any_write) then
3274		!$omp target update from(v_old, vel_report)
3275	0	do J=Jsq,Jeq ; do i=is,ie ; if (dowrite(i,J)) then
3276		! Call a diagnostic reporting subroutines are called if unphysically large values are found.
3277		call write_v_accel(i, J, v_old, h, ADp, CDp, dt, G, GV, US, CS%PointAccel_CSp, &
3278	0	vel_report(i,J), forces%tauy(i,J), a=CS%a_v, hv=CS%h_v)
3279		endif ; enddo ; enddo
3280		endif
3281		else ! Do not report accelerations leading to large velocities.
3282	48	ntrunc = 0
3283	48	do concurrent (k=1:nz, J=Jsq:Jeq, i=is:ie) DO_LOCALITY(reduce(+: ntrunc))
3284	26709168	if (abs(v(i,J,k)) < CS%vel_underflow) then ; v(i,J,k) = 0.0
3285	26352000	elseif ((v(i,J,k) * (dt * G%dx_Cv(i,J))) * G%IareaT(i,j+1) < -CS%CFL_trunc) then
3286	0	v(i,J,k) = (-0.9CS%CFL_trunc) (G%areaT(i,j+1) / (dt * G%dx_Cv(i,J)))
3287	0	if (((i >= G%isc) .and. (i <= G%iec) .and. (J >= G%jsc) .and. (J <= G%jec)) .and. &
3288	0	(CS%h_v(i,J,k) > H_report)) ntrunc = ntrunc + 1
3289	26352000	elseif ((v(i,J,k) * (dt * G%dx_Cv(i,J))) * G%IareaT(i,j) > CS%CFL_trunc) then
3290	0	v(i,J,k) = (0.9CS%CFL_trunc) (G%areaT(i,j) / (dt * G%dx_Cv(i,J)))
3291	0	if (((i >= G%isc) .and. (i <= G%iec) .and. (J >= G%jsc) .and. (J <= G%jec)) .and. &
3292	0	(CS%h_v(i,J,k) > H_report)) ntrunc = ntrunc + 1
3293		endif
3294		enddo
3295	48	CS%ntrunc = CS%ntrunc + ntrunc
3296		endif
3297
3298		!$omp target exit data map(release: u_old, v_old, dowrite, vel_report)
3299	48	end subroutine vertvisc_limit_vel
3300
3301
3302		!> Initialize the vertical friction module
3303	1	subroutine vertvisc_init(MIS, Time, G, GV, US, param_file, diag, ADp, dirs, &
3304		ntrunc, CS, fpmix)
3305		type(ocean_internal_state), &
3306		target, intent(in) :: MIS !< The "MOM Internal State", a set of pointers
3307		!! to the fields and accelerations that make
3308		!! up the ocean's physical state
3309		type(time_type), target, intent(in) :: Time !< Current model time
3310		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
3311		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
3312		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3313		type(param_file_type), intent(in) :: param_file !< File to parse for parameters
3314		type(diag_ctrl), target, intent(inout) :: diag !< Diagnostic control structure
3315		type(accel_diag_ptrs), intent(inout) :: ADp !< Acceleration diagnostic pointers
3316		type(directories), intent(in) :: dirs !< Relevant directory paths
3317		integer, target, intent(inout) :: ntrunc !< Number of velocity truncations
3318		type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
3319		logical, optional, intent(in) :: fpmix !< Nonlocal momentum mixing
3320
3321		! Local variables
3322
3323		real :: Kv_BBL ! A viscosity in the bottom boundary layer with a simple scheme [H Z T-1 ~> m2 s-1 or Pa s]
3324		real :: Kv_back_z ! A background kinematic viscosity [Z2 T-1 ~> m2 s-1]
3325		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
3326		integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
3327		logical :: lfpmix
3328		character(len=200) :: kappa_gl90_file, inputdir, kdgl90_varname
3329		! This include declares and sets the variable "version".
3330		# include "version_variable.h"
3331		character(len=40) :: mdl = "MOM_vert_friction" ! This module's name.
3332		character(len=40) :: thickness_units
3333		real :: Kv_mks ! KVML in MKS [m2 s-1]
3334
3335	1	CS%initialized = .true.
3336
3337	1	if (GV%Boussinesq) then ; thickness_units = "m"
3338	0	else ; thickness_units = "kg m-2" ; endif
3339
3340	1	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke
3341	1	IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
3342
3343	1	CS%diag => diag ; CS%ntrunc => ntrunc ; ntrunc = 0
3344
3345	1	lfpmix = .false.
3346	1	if (present(fpmix)) lfpmix = fpmix
3347
3348		! Default, read and log parameters
3349	1	call log_version(param_file, mdl, version, "", log_to_all=.true., debugging=.true.)
3350		call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
3351		"This sets the default value for the various _ANSWER_DATE parameters.", &
3352	1	default=99991231)
3353		call get_param(param_file, mdl, "VERT_FRICTION_ANSWER_DATE", CS%answer_date, &
3354		"The vintage of the order of arithmetic and expressions in the viscous "//&
3355		"calculations. Values below 20190101 recover the answers from the end of 2018, "//&
3356		"while higher values use expressions that do not use an arbitrary hard-coded "//&
3357		"maximum viscous coupling coefficient between layers. Values below 20230601 "//&
3358		"recover a form of the viscosity within the mixed layer that breaks up the "//&
3359		"magnitude of the wind stress in some non-Boussinesq cases.", &
3360	1	default=default_answer_date, do_not_log=.not.GV%Boussinesq)
3361	1	if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701)
3362
3363		call get_param(param_file, mdl, "BOTTOMDRAGLAW", CS%bottomdraglaw, &
3364		"If true, the bottom stress is calculated with a drag "//&
3365		"law of the form c_drag\|u\|u. The velocity magnitude "//&
3366		"may be an assumed value or it may be based on the "//&
3367		"actual velocity in the bottommost HBBL, depending on "//&
3368	1	"LINEAR_DRAG.", default=.true.)
3369		call get_param(param_file, mdl, "DIRECT_STRESS", CS%direct_stress, &
3370		"If true, the wind stress is distributed over the topmost HMIX_STRESS of fluid "//&
3371		"(like in HYCOM), and an added mixed layer viscosity or a physically based "//&
3372		"boundary layer turbulence parameterization is not needed for stability.", &
3373	1	default=.false.)
3374		call get_param(param_file, mdl, "DYNAMIC_VISCOUS_ML", CS%dynamic_viscous_ML, &
3375		"If true, use a bulk Richardson number criterion to "//&
3376		"determine the mixed layer thickness for viscosity.", &
3377	1	default=.false.)
3378		call get_param(param_file, mdl, "FIXED_DEPTH_LOTW_ML", CS%fixed_LOTW_ML, &
3379		"If true, use a Law-of-the-wall prescription for the mixed layer viscosity "//&
3380		"within a boundary layer that is the lesser of HMIX_FIXED and the total "//&
3381	1	"depth of the ocean in a column.", default=.false.)
3382		call get_param(param_file, mdl, "LOTW_VISCOUS_ML_FLOOR", CS%apply_LOTW_floor, &
3383		"If true, use a Law-of-the-wall prescription to set a lower bound on the "//&
3384		"viscous coupling between layers within the surface boundary layer, based "//&
3385		"the distance of interfaces from the surface. This only acts when there "//&
3386		"are large changes in the thicknesses of successive layers or when the "//&
3387		"viscosity is set externally and the wind stress has subsequently increased.", &
3388	1	default=.false.)
3389		call get_param(param_file, mdl, 'VON_KARMAN_CONST', CS%vonKar, &
3390		'The value the von Karman constant as used for mixed layer viscosity.', &
3391	1	units='nondim', default=0.41)
3392		call get_param(param_file, mdl, "U_TRUNC_FILE", CS%u_trunc_file, &
3393		"The absolute path to a file into which the accelerations "//&
3394		"leading to zonal velocity truncations are written. "//&
3395		"Undefine this for efficiency if this diagnostic is not needed.", &
3396	1	default=" ", debuggingParam=.true.)
3397		call get_param(param_file, mdl, "V_TRUNC_FILE", CS%v_trunc_file, &
3398		"The absolute path to a file into which the accelerations "//&
3399		"leading to meridional velocity truncations are written. "//&
3400		"Undefine this for efficiency if this diagnostic is not needed.", &
3401	1	default=" ", debuggingParam=.true.)
3402		call get_param(param_file, mdl, "HARMONIC_VISC", CS%harmonic_visc, &
3403		"If true, use the harmonic mean thicknesses for "//&
3404	1	"calculating the vertical viscosity.", default=.false.)
3405		call get_param(param_file, mdl, "HARMONIC_BL_SCALE", CS%harm_BL_val, &
3406		"A scale to determine when water is in the boundary "//&
3407		"layers based solely on harmonic mean thicknesses for "//&
3408		"the purpose of determining the extent to which the "//&
3409		"thicknesses used in the viscosities are upwinded.", &
3410	1	default=0.0, units="nondim")
3411	1	call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.)
3412
3413	1	if (GV%nkml < 1) then
3414		call get_param(param_file, mdl, "HMIX_FIXED", CS%Hmix, &
3415		"The prescribed depth over which the near-surface viscosity and "//&
3416		"diffusivity are elevated when the bulk mixed layer is not used.", &
3417	1	units="m", scale=US%m_to_Z, fail_if_missing=.true.)
3418		endif
3419	1	if (CS%direct_stress) then
3420	0	if (GV%nkml < 1) then
3421		call get_param(param_file, mdl, "HMIX_STRESS", CS%Hmix_stress, &
3422		"The depth over which the wind stress is applied if DIRECT_STRESS is true.", &
3423	0	units="m", default=US%Z_to_m*CS%Hmix, scale=GV%m_to_H)
3424		else
3425		call get_param(param_file, mdl, "HMIX_STRESS", CS%Hmix_stress, &
3426		"The depth over which the wind stress is applied if DIRECT_STRESS is true.", &
3427	0	units="m", fail_if_missing=.true., scale=GV%m_to_H)
3428		endif
3429	0	if (CS%Hmix_stress <= 0.0) call MOM_error(FATAL, "vertvisc_init: " // &
3430	0	"HMIX_STRESS must be set to a positive value if DIRECT_STRESS is true.")
3431		endif
3432		call get_param(param_file, mdl, "KV", Kv_back_z, &
3433		"The background kinematic viscosity in the interior. "//&
3434		"The molecular value, ~1e-6 m2 s-1, may be used.", &
3435	1	units="m2 s-1", fail_if_missing=.true., scale=US%m2_s_to_Z2_T)
3436		! Convert input kinematic viscosity to dynamic viscosity when non-Boussinesq.
3437	1	CS%Kv = (US%Z2_T_to_m2_sGV%m2_s_to_HZ_T) Kv_back_z
3438
3439		call get_param(param_file, mdl, "USE_GL90_IN_SSW", CS%use_GL90_in_SSW, &
3440		"If true, use simpler method to calculate 1/N^2 in GL90 vertical "// &
3441		"viscosity coefficient. This method is valid in stacked shallow water mode.", &
3442	1	default=.false.)
3443		call get_param(param_file, mdl, "KD_GL90", CS%kappa_gl90, &
3444		"The scalar diffusivity used in GL90 vertical viscosity scheme.", &
3445		units="m2 s-1", default=0.0, scale=US%m_to_LUS%Z_to_LGV%m_to_H*US%T_to_s, &
3446	1	do_not_log=.not.CS%use_GL90_in_SSW)
3447		call get_param(param_file, mdl, "READ_KD_GL90", CS%read_kappa_gl90, &
3448		"If true, read a file (given by KD_GL90_FILE) containing the "//&
3449		"spatially varying diffusivity KD_GL90 used in the GL90 scheme.", default=.false., &
3450	1	do_not_log=.not.CS%use_GL90_in_SSW)
3451	1	if (CS%read_kappa_gl90) then
3452	0	if (CS%kappa_gl90 > 0) then
3453		call MOM_error(FATAL, "MOM_vert_friction.F90, vertvisc_init: KD_GL90 > 0 "// &
3454	0	"is not compatible with READ_KD_GL90 = .TRUE. ")
3455		endif
3456		call get_param(param_file, mdl, "INPUTDIR", inputdir, &
3457		"The directory in which all input files are found.", &
3458	0	default=".", do_not_log=.true.)
3459	0	inputdir = slasher(inputdir)
3460		call get_param(param_file, mdl, "KD_GL90_FILE", kappa_gl90_file, &
3461		"The file containing the spatially varying diffusivity used in the "// &
3462	0	"GL90 scheme.", default="kd_gl90.nc", do_not_log=.not.CS%use_GL90_in_SSW)
3463		call get_param(param_file, mdl, "KD_GL90_VARIABLE", kdgl90_varname, &
3464		"The name of the GL90 diffusivity variable to read "//&
3465	0	"from KD_GL90_FILE.", default="kd_gl90", do_not_log=.not.CS%use_GL90_in_SSW)
3466	0	kappa_gl90_file = trim(inputdir) // trim(kappa_gl90_file)
3467
3468	0	allocate(CS%kappa_gl90_2d(G%isd:G%ied, G%jsd:G%jed), source=0.0)
3469		call MOM_read_data(kappa_gl90_file, kdgl90_varname, CS%kappa_gl90_2d(:,:), G%domain, &
3470	0	scale=US%m_to_LUS%Z_to_LGV%m_to_H*US%T_to_s)
3471	0	call pass_var(CS%kappa_gl90_2d, G%domain)
3472		endif
3473		call get_param(param_file, mdl, "USE_GL90_N2", CS%use_GL90_N2, &
3474		"If true, use GL90 vertical viscosity coefficient that is depth-independent; "// &
3475		"this corresponds to a kappa_GM that scales as N^2 with depth.", &
3476	1	default=.false., do_not_log=.not.CS%use_GL90_in_SSW)
3477	1	if (CS%use_GL90_N2) then
3478	0	if (.not. CS%use_GL90_in_SSW) call MOM_error(FATAL, &
3479		"MOM_vert_friction.F90, vertvisc_init: "//&
3480	0	"When USE_GL90_N2=True, USE_GL90_in_SSW must also be True.")
3481	0	if (CS%kappa_gl90 > 0) then
3482		call MOM_error(FATAL, "MOM_vert_friction.F90, vertvisc_init: KD_GL90 > 0 "// &
3483	0	"is not compatible with USE_GL90_N2 = .TRUE. ")
3484		endif
3485	0	if (CS%read_kappa_gl90) call MOM_error(FATAL, &
3486		"MOM_vert_friction.F90, vertvisc_init: "//&
3487	0	"READ_KD_GL90 = .TRUE. is not compatible with USE_GL90_N2 = .TRUE.")
3488		call get_param(param_file, mdl, "alpha_GL90", CS%alpha_gl90, &
3489		"Coefficient used to compute a depth-independent GL90 vertical "//&
3490		"viscosity via Kv_GL90 = alpha_GL90 * f2. Is only used "// &
3491		"if USE_GL90_N2 is true. Note that the implied Kv_GL90 "// &
3492		"corresponds to a KD_GL90 that scales as N^2 with depth.", &
3493		units="m2 s", default=0.0, scale=GV%m_to_HUS%m_to_ZUS%s_to_T, &
3494	0	do_not_log=.not.CS%use_GL90_in_SSW)
3495		endif
3496		call get_param(param_file, mdl, "HBBL_GL90", CS%Hbbl_gl90, &
3497		"The thickness of the GL90 bottom boundary layer, "//&
3498		"which defines the range over which the GL90 coupling "//&
3499		"coefficient is zeroed out, in order to avoid fluxing "//&
3500		"momentum into vanished layers over steep topography.", &
3501	1	units="m", default=5.0, scale=US%m_to_Z, do_not_log=.not.CS%use_GL90_in_SSW)
3502
3503	1	CS%Kvml_invZ2 = 0.0
3504	1	if (GV%nkml < 1) then
3505		call get_param(param_file, mdl, "KVML", Kv_mks, &
3506		"The scale for an extra kinematic viscosity in the mixed layer", &
3507	1	units="m2 s-1", default=-1.0, do_not_log=.true.)
3508	1	if (Kv_mks >= 0.0) then
3509	0	call MOM_error(WARNING, "KVML is a deprecated parameter. Use KV_ML_INVZ2 instead.")
3510		else
3511	1	Kv_mks = 0.0
3512		endif
3513		call get_param(param_file, mdl, "KV_ML_INVZ2", CS%Kvml_invZ2, &
3514		"An extra kinematic viscosity in a mixed layer of thickness HMIX_FIXED, "//&
3515		"with the actual viscosity scaling as 1/(z*HMIX_FIXED)^2, where z is the "//&
3516		"distance from the surface, to allow for finite wind stresses to be "//&
3517		"transmitted through infinitesimally thin surface layers. This is an "//&
3518		"older option for numerical convenience without a strong physical basis, "//&
3519		"and its use is now discouraged.", &
3520	1	units="m2 s-1", default=Kv_mks, scale=GV%m2_s_to_HZ_T)
3521		endif
3522
3523	1	if (.not.CS%bottomdraglaw) then
3524		call get_param(param_file, mdl, "KV_EXTRA_BBL", CS%Kv_extra_bbl, &
3525		"An extra kinematic viscosity in the benthic boundary layer. "//&
3526		"KV_EXTRA_BBL is not used if BOTTOMDRAGLAW is true.", &
3527	0	units="m2 s-1", default=0.0, scale=GV%m2_s_to_HZ_T, do_not_log=.true.)
3528	0	if (CS%Kv_extra_bbl == 0.0) then
3529		call get_param(param_file, mdl, "KVBBL", Kv_BBL, &
3530		"An extra kinematic viscosity in the benthic boundary layer. "//&
3531		"KV_EXTRA_BBL is not used if BOTTOMDRAGLAW is true.", &
3532		units="m2 s-1", default=US%Z2_T_to_m2_s*Kv_back_z, scale=GV%m2_s_to_HZ_T, &
3533	0	do_not_log=.true.)
3534	0	if (abs(Kv_BBL - CS%Kv) > 1.0e-15*abs(CS%Kv)) then
3535	0	call MOM_error(WARNING, "KVBBL is a deprecated parameter. Use KV_EXTRA_BBL instead.")
3536	0	CS%Kv_extra_bbl = Kv_BBL - CS%Kv
3537		endif
3538		endif
3539		call log_param(param_file, mdl, "KV_EXTRA_BBL", CS%Kv_extra_bbl, &
3540		"An extra kinematic viscosity in the benthic boundary layer. "//&
3541		"KV_EXTRA_BBL is not used if BOTTOMDRAGLAW is true.", &
3542	0	units="m2 s-1", default=0.0, unscale=GV%HZ_T_to_m2_s)
3543		endif
3544		call get_param(param_file, mdl, "HBBL", CS%Hbbl, &
3545		"The thickness of a bottom boundary layer with a viscosity increased by "//&
3546		"KV_EXTRA_BBL if BOTTOMDRAGLAW is not defined, or the thickness over which "//&
3547		"near-bottom velocities are averaged for the drag law if BOTTOMDRAGLAW is "//&
3548		"defined but LINEAR_DRAG is not.", &
3549	1	units="m", fail_if_missing=.true., scale=US%m_to_Z)
3550		call get_param(param_file, mdl, "CFL_TRUNCATE", CS%CFL_trunc, &
3551		"The value of the CFL number that will cause velocity "//&
3552		"components to be truncated; instability can occur past 0.5.", &
3553	1	units="nondim", default=0.5)
3554		call get_param(param_file, mdl, "CFL_REPORT", CS%CFL_report, &
3555		"The value of the CFL number that causes accelerations "//&
3556		"to be reported; the default is CFL_TRUNCATE.", &
3557	1	units="nondim", default=CS%CFL_trunc)
3558		call get_param(param_file, mdl, "CFL_TRUNCATE_RAMP_TIME", CS%truncRampTime, &
3559		"The time over which the CFL truncation value is ramped "//&
3560		"up at the beginning of the run.", &
3561	1	units="s", default=0., scale=US%s_to_T)
3562	1	CS%CFL_truncE = CS%CFL_trunc
3563		call get_param(param_file, mdl, "CFL_TRUNCATE_START", CS%CFL_truncS, &
3564		"The start value of the truncation CFL number used when "//&
3565		"ramping up CFL_TRUNC.", &
3566	1	units="nondim", default=0.)
3567		call get_param(param_file, mdl, "STOKES_MIXING_COMBINED", CS%StokesMixing, &
3568		"Flag to use Stokes drift Mixing via the Lagrangian "//&
3569		" current (Eulerian plus Stokes drift). "//&
3570		" Still needs work and testing, so not recommended for use.",&
3571	1	default=.false.)
3572		!BGR 04/04/2018{
3573		! StokesMixing is required for MOM6 for some Langmuir mixing parameterization.
3574		! The code used here has not been developed for vanishing layers or in
3575		! conjunction with any bottom friction. Therefore, the following line is
3576		! added so this functionality cannot be used without user intervention in
3577		! the code. This will prevent general use of this functionality until proper
3578		! care is given to the previously mentioned issues. Comment out the following
3579		! MOM_error to use, but do so at your own risk and with these points in mind.
3580		!}
3581	1	if (CS%StokesMixing) then
3582		call MOM_error(FATAL, "Stokes mixing requires user intervention in the code.\n"//&
3583		" Model now exiting. See MOM_vert_friction.F90 for \n"//&
3584	0	" details (search 'BGR 04/04/2018' to locate comment).")
3585		endif
3586		call get_param(param_file, mdl, "VEL_UNDERFLOW", CS%vel_underflow, &
3587		"A negligibly small velocity magnitude below which velocity "//&
3588		"components are set to 0. A reasonable value might be "//&
3589		"1e-30 m/s, which is less than an Angstrom divided by "//&
3590	1	"the age of the universe.", units="m s-1", default=0.0, scale=US%m_s_to_L_T)
3591
3592		!$omp target update to(CS)
3593
3594	671917	ALLOC_(CS%a_u(IsdB:IedB,jsd:jed,nz+1)) ; CS%a_u(:,:,:) = 0.0
3595	671917	ALLOC_(CS%a_u_gl90(IsdB:IedB,jsd:jed,nz+1)) ; CS%a_u_gl90(:,:,:) = 0.0
3596	663076	ALLOC_(CS%h_u(IsdB:IedB,jsd:jed,nz)) ; CS%h_u(:,:,:) = 0.0
3597	676553	ALLOC_(CS%a_v(isd:ied,JsdB:JedB,nz+1)) ; CS%a_v(:,:,:) = 0.0
3598	676553	ALLOC_(CS%a_v_gl90(isd:ied,JsdB:JedB,nz+1)) ; CS%a_v_gl90(:,:,:) = 0.0
3599	667651	ALLOC_(CS%h_v(isd:ied,JsdB:JedB,nz)) ; CS%h_v(:,:,:) = 0.0
3600
3601		!$omp target enter data map(to: CS%a_u, CS%a_v)
3602		!$omp target enter data map(to: CS%h_u, CS%h_v)
3603		! TODO: Conditional?
3604		!$omp target enter data map(to: CS%a_u_gl90, CS%a_v_gl90)
3605
3606		CS%id_Kv_slow = register_diag_field('ocean_model', 'Kv_slow', diag%axesTi, Time, &
3607	1	'Slow varying vertical viscosity', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
3608
3609		CS%id_Kv_u = register_diag_field('ocean_model', 'Kv_u', diag%axesCuL, Time, &
3610	1	'Total vertical viscosity at u-points', 'm2 s-1', conversion=GV%H_to_m*2US%s_to_T)
3611
3612		CS%id_Kv_v = register_diag_field('ocean_model', 'Kv_v', diag%axesCvL, Time, &
3613	1	'Total vertical viscosity at v-points', 'm2 s-1', conversion=GV%H_to_m*2US%s_to_T)
3614
3615		CS%id_Kv_gl90_u = register_diag_field('ocean_model', 'Kv_gl90_u', diag%axesCuL, Time, &
3616	1	'GL90 vertical viscosity at u-points', 'm2 s-1', conversion=GV%H_to_m*2US%s_to_T)
3617
3618		CS%id_Kv_gl90_v = register_diag_field('ocean_model', 'Kv_gl90_v', diag%axesCvL, Time, &
3619	1	'GL90 vertical viscosity at v-points', 'm2 s-1', conversion=GV%H_to_m*2US%s_to_T)
3620
3621		CS%id_au_vv = register_diag_field('ocean_model', 'au_visc', diag%axesCui, Time, &
3622	1	'Zonal Viscous Vertical Coupling Coefficient', 'm s-1', conversion=GV%H_to_m*US%s_to_T)
3623
3624		CS%id_av_vv = register_diag_field('ocean_model', 'av_visc', diag%axesCvi, Time, &
3625	1	'Meridional Viscous Vertical Coupling Coefficient', 'm s-1', conversion=GV%H_to_m*US%s_to_T)
3626
3627		CS%id_au_gl90_vv = register_diag_field('ocean_model', 'au_gl90_visc', diag%axesCui, Time, &
3628	1	'Zonal Viscous Vertical GL90 Coupling Coefficient', 'm s-1', conversion=GV%H_to_m*US%s_to_T)
3629
3630		CS%id_av_gl90_vv = register_diag_field('ocean_model', 'av_gl90_visc', diag%axesCvi, Time, &
3631	1	'Meridional Viscous Vertical GL90 Coupling Coefficient', 'm s-1', conversion=GV%H_to_m*US%s_to_T)
3632
3633		CS%id_h_u = register_diag_field('ocean_model', 'Hu_visc', diag%axesCuL, Time, &
3634		'Thickness at Zonal Velocity Points for Viscosity', &
3635	1	thickness_units, conversion=GV%H_to_MKS)
3636		! Alternately, to always give this variable in 'm' use the following line instead:
3637		! 'm', conversion=GV%H_to_m)
3638
3639		CS%id_h_v = register_diag_field('ocean_model', 'Hv_visc', diag%axesCvL, Time, &
3640		'Thickness at Meridional Velocity Points for Viscosity', &
3641	1	thickness_units, conversion=GV%H_to_MKS)
3642
3643		CS%id_hML_u = register_diag_field('ocean_model', 'HMLu_visc', diag%axesCu1, Time, &
3644		'Mixed Layer Thickness at Zonal Velocity Points for Viscosity', &
3645	1	thickness_units, conversion=US%Z_to_m)
3646
3647		CS%id_hML_v = register_diag_field('ocean_model', 'HMLv_visc', diag%axesCv1, Time, &
3648		'Mixed Layer Thickness at Meridional Velocity Points for Viscosity', &
3649	1	thickness_units, conversion=US%Z_to_m)
3650
3651	1	if (lfpmix) then
3652		CS%id_uE_h = register_diag_field('ocean_model', 'uE_h' , CS%diag%axesTL, &
3653	0	Time, 'x-zonal Eulerian' , 'm s-1', conversion=US%L_T_to_m_s)
3654		CS%id_vE_h = register_diag_field('ocean_model', 'vE_h' , CS%diag%axesTL, &
3655	0	Time, 'y-merid Eulerian' , 'm s-1', conversion=US%L_T_to_m_s)
3656		CS%id_uInc_h = register_diag_field('ocean_model','uInc_h',CS%diag%axesTL, &
3657	0	Time, 'x-zonal Eulerian' , 'm s-1', conversion=US%L_T_to_m_s)
3658		CS%id_vInc_h = register_diag_field('ocean_model','vInc_h',CS%diag%axesTL, &
3659	0	Time, 'x-zonal Eulerian' , 'm s-1', conversion=US%L_T_to_m_s)
3660		CS%id_uStk = register_diag_field('ocean_model', 'uStk' , CS%diag%axesTL, &
3661	0	Time, 'x-FP du increment' , 'm s-1', conversion=US%L_T_to_m_s)
3662		CS%id_vStk = register_diag_field('ocean_model', 'vStk' , CS%diag%axesTL, &
3663	0	Time, 'y-FP dv increment' , 'm s-1', conversion=US%L_T_to_m_s)
3664
3665		CS%id_FPtau2s = register_diag_field('ocean_model','Omega_tau2s',CS%diag%axesTi, &
3666	0	Time, 'Stress direction from shear','radians')
3667		CS%id_FPtau2w = register_diag_field('ocean_model','Omega_tau2w',CS%diag%axesTi, &
3668	0	Time, 'Stress direction from wind','radians')
3669		CS%id_uStk0 = register_diag_field('ocean_model', 'uStk0' , diag%axesT1, &
3670	0	Time, 'Zonal Surface Stokes', 'm s-1', conversion=US%L_T_to_m_s)
3671		CS%id_vStk0 = register_diag_field('ocean_model', 'vStk0' , diag%axesT1, &
3672	0	Time, 'Merid Surface Stokes', 'm s-1', conversion=US%L_T_to_m_s)
3673		endif
3674
3675		CS%id_du_dt_visc = register_diag_field('ocean_model', 'du_dt_visc', diag%axesCuL, Time, &
3676	1	'Zonal Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2)
3677	1	if (CS%id_du_dt_visc > 0) call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
3678		CS%id_dv_dt_visc = register_diag_field('ocean_model', 'dv_dt_visc', diag%axesCvL, Time, &
3679	1	'Meridional Acceleration from Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2)
3680	1	if (CS%id_dv_dt_visc > 0) call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
3681		CS%id_GLwork = register_diag_field('ocean_model', 'GLwork', diag%axesTL, Time, &
3682		'Sign-definite Kinetic Energy Source from GL90 Vertical Viscosity', &
3683	1	'm3 s-3', conversion=GV%H_to_m(US%L_T_to_m_s2)US%s_to_T)
3684		CS%id_du_dt_visc_gl90 = register_diag_field('ocean_model', 'du_dt_visc_gl90', diag%axesCuL, Time, &
3685	1	'Zonal Acceleration from GL90 Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2)
3686	1	if ((CS%id_du_dt_visc_gl90 > 0) .or. (CS%id_GLwork > 0)) then
3687	0	call safe_alloc_ptr(ADp%du_dt_visc_gl90,IsdB,IedB,jsd,jed,nz)
3688	0	call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
3689		endif
3690		CS%id_dv_dt_visc_gl90 = register_diag_field('ocean_model', 'dv_dt_visc_gl90', diag%axesCvL, Time, &
3691	1	'Meridional Acceleration from GL90 Vertical Viscosity', 'm s-2', conversion=US%L_T2_to_m_s2)
3692	1	if ((CS%id_dv_dt_visc_gl90 > 0) .or. (CS%id_GLwork > 0)) then
3693	0	call safe_alloc_ptr(ADp%dv_dt_visc_gl90,isd,ied,JsdB,JedB,nz)
3694	0	call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
3695		endif
3696		CS%id_du_dt_str = register_diag_field('ocean_model', 'du_dt_str', diag%axesCuL, Time, &
3697	1	'Zonal Acceleration from Surface Wind Stresses', 'm s-2', conversion=US%L_T2_to_m_s2)
3698	1	if (CS%id_du_dt_str > 0) call safe_alloc_ptr(ADp%du_dt_str,IsdB,IedB,jsd,jed,nz)
3699		CS%id_dv_dt_str = register_diag_field('ocean_model', 'dv_dt_str', diag%axesCvL, Time, &
3700	1	'Meridional Acceleration from Surface Wind Stresses', 'm s-2', conversion=US%L_T2_to_m_s2)
3701	1	if (CS%id_dv_dt_str > 0) call safe_alloc_ptr(ADp%dv_dt_str,isd,ied,JsdB,JedB,nz)
3702
3703		CS%id_taux_bot = register_diag_field('ocean_model', 'taux_bot', diag%axesCu1, &
3704		Time, 'Zonal Bottom Stress from Ocean to Earth', &
3705	1	'Pa', conversion=US%RZ_to_kg_m2*US%L_T2_to_m_s2)
3706		CS%id_tauy_bot = register_diag_field('ocean_model', 'tauy_bot', diag%axesCv1, &
3707		Time, 'Meridional Bottom Stress from Ocean to Earth', &
3708	1	'Pa', conversion=US%RZ_to_kg_m2*US%L_T2_to_m_s2)
3709
3710		!CS%id_hf_du_dt_visc = register_diag_field('ocean_model', 'hf_du_dt_visc', diag%axesCuL, Time, &
3711		! 'Fractional Thickness-weighted Zonal Acceleration from Vertical Viscosity', &
3712		! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
3713		!if (CS%id_hf_du_dt_visc > 0) then
3714		! call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
3715		! call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
3716		!endif
3717
3718		!CS%id_hf_dv_dt_visc = register_diag_field('ocean_model', 'hf_dv_dt_visc', diag%axesCvL, Time, &
3719		! 'Fractional Thickness-weighted Meridional Acceleration from Vertical Viscosity', &
3720		! 'm s-2', v_extensive=.true., conversion=US%L_T2_to_m_s2)
3721		!if (CS%id_hf_dv_dt_visc > 0) then
3722		! call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
3723		! call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
3724		!endif
3725
3726		CS%id_hf_du_dt_visc_2d = register_diag_field('ocean_model', 'hf_du_dt_visc_2d', diag%axesCu1, Time, &
3727		'Depth-sum Fractional Thickness-weighted Zonal Acceleration from Vertical Viscosity', &
3728	1	'm s-2', conversion=US%L_T2_to_m_s2)
3729	1	if (CS%id_hf_du_dt_visc_2d > 0) then
3730	0	call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
3731	0	call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
3732		endif
3733
3734		CS%id_hf_dv_dt_visc_2d = register_diag_field('ocean_model', 'hf_dv_dt_visc_2d', diag%axesCv1, Time, &
3735		'Depth-sum Fractional Thickness-weighted Meridional Acceleration from Vertical Viscosity', &
3736	1	'm s-2', conversion=US%L_T2_to_m_s2)
3737	1	if (CS%id_hf_dv_dt_visc_2d > 0) then
3738	0	call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
3739	0	call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,JsdB,JedB,nz)
3740		endif
3741
3742		CS%id_h_du_dt_visc = register_diag_field('ocean_model', 'h_du_dt_visc', diag%axesCuL, Time, &
3743		'Thickness Multiplied Zonal Acceleration from Horizontal Viscosity', &
3744	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
3745	1	if (CS%id_h_du_dt_visc > 0) then
3746	0	call safe_alloc_ptr(ADp%du_dt_visc,IsdB,IedB,jsd,jed,nz)
3747	0	call safe_alloc_ptr(ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
3748		endif
3749
3750		CS%id_h_dv_dt_visc = register_diag_field('ocean_model', 'h_dv_dt_visc', diag%axesCvL, Time, &
3751		'Thickness Multiplied Meridional Acceleration from Horizontal Viscosity', &
3752	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
3753	1	if (CS%id_h_dv_dt_visc > 0) then
3754	0	call safe_alloc_ptr(ADp%dv_dt_visc,isd,ied,JsdB,JedB,nz)
3755	0	call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,JedB,nz)
3756		endif
3757
3758		CS%id_h_du_dt_str = register_diag_field('ocean_model', 'h_du_dt_str', diag%axesCuL, Time, &
3759		'Thickness Multiplied Zonal Acceleration from Surface Wind Stresses', &
3760	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
3761	1	if (CS%id_h_du_dt_str > 0) then
3762	0	call safe_alloc_ptr(ADp%du_dt_str,IsdB,IedB,jsd,jed,nz)
3763	0	call safe_alloc_ptr(ADp%diag_hu,IsdB,IedB,jsd,jed,nz)
3764		endif
3765
3766		CS%id_h_dv_dt_str = register_diag_field('ocean_model', 'h_dv_dt_str', diag%axesCvL, Time, &
3767		'Thickness Multiplied Meridional Acceleration from Surface Wind Stresses', &
3768	1	'm2 s-2', conversion=GV%H_to_m*US%L_T2_to_m_s2)
3769	1	if (CS%id_h_dv_dt_str > 0) then
3770	0	call safe_alloc_ptr(ADp%dv_dt_str,isd,ied,JsdB,JedB,nz)
3771	0	call safe_alloc_ptr(ADp%diag_hv,isd,ied,JsdB,JedB,nz)
3772		endif
3773
3774		CS%id_du_dt_str_visc_rem = register_diag_field('ocean_model', 'du_dt_str_visc_rem', diag%axesCuL, Time, &
3775		'Zonal Acceleration from Surface Wind Stresses multiplied by viscous remnant', &
3776	1	'm s-2', conversion=US%L_T2_to_m_s2)
3777	1	if (CS%id_du_dt_str_visc_rem > 0) then
3778	0	call safe_alloc_ptr(ADp%du_dt_str,IsdB,IedB,jsd,jed,nz)
3779	0	call safe_alloc_ptr(ADp%visc_rem_u,IsdB,IedB,jsd,jed,nz)
3780		endif
3781
3782		CS%id_dv_dt_str_visc_rem = register_diag_field('ocean_model', 'dv_dt_str_visc_rem', diag%axesCvL, Time, &
3783		'Meridional Acceleration from Surface Wind Stresses multiplied by viscous remnant', &
3784	1	'm s-2', conversion=US%L_T2_to_m_s2)
3785	1	if (CS%id_dv_dt_str_visc_rem > 0) then
3786	0	call safe_alloc_ptr(ADp%dv_dt_str,isd,ied,JsdB,JedB,nz)
3787	0	call safe_alloc_ptr(ADp%visc_rem_v,isd,ied,JsdB,JedB,nz)
3788		endif
3789
3790	1	if ((len_trim(CS%u_trunc_file) > 0) .or. (len_trim(CS%v_trunc_file) > 0)) &
3791	0	call PointAccel_init(MIS, Time, G, param_file, diag, dirs, CS%PointAccel_CSp)
3792
3793	1	end subroutine vertvisc_init
3794
3795		!> Update the CFL truncation value as a function of time.
3796		!! If called with the optional argument activate=.true., record the
3797		!! value of Time as the beginning of the ramp period.
3798	25	subroutine updateCFLtruncationValue(Time, CS, US, activate)
3799		type(time_type), target, intent(in) :: Time !< Current model time
3800		type(vertvisc_CS), pointer :: CS !< Vertical viscosity control structure
3801		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3802		logical, optional, intent(in) :: activate !< Specify whether to record the value of
3803		!! Time as the beginning of the ramp period
3804
3805		! Local variables
3806		real :: deltaTime ! The time since CS%rampStartTime [T ~> s], which may be negative.
3807		real :: wghtA ! The relative weight of the final value [nondim]
3808		character(len=12) :: msg
3809
3810	25	if (CS%truncRampTime==0.) return ! This indicates to ramping is turned off
3811
3812		! We use the optional argument to indicate this Time should be recorded as the
3813		! beginning of the ramp-up period.
3814	0	if (present(activate)) then
3815	0	if (activate) then
3816	0	CS%rampStartTime = Time ! Record the current time
3817	0	CS%CFLrampingIsActivated = .true.
3818		endif
3819		endif
3820	0	if (.not.CS%CFLrampingIsActivated) return
3821	0	deltaTime = max(0., time_minus_signed(Time, CS%rampStartTime, scale=US%s_to_T))
3822	0	if (deltaTime >= CS%truncRampTime) then
3823	0	CS%CFL_trunc = CS%CFL_truncE
3824	0	CS%truncRampTime = 0. ! This turns off ramping after this call
3825		else
3826	0	wghtA = min( 1., deltaTime / CS%truncRampTime ) ! Linear profile in time
3827		!wghtA = wghtA*wghtA ! Convert linear profile to parabolic profile in time
3828		!wghtA = wghtAwghtA(3. - 2.*wghtA) ! Convert linear profile to cosine profile
3829	0	wghtA = 1. - ( (1. - wghtA)**2 ) ! Convert linear profile to inverted parabolic profile
3830	0	CS%CFL_trunc = CS%CFL_truncS + wghtA * ( CS%CFL_truncE - CS%CFL_truncS )
3831		endif
3832	0	write(msg(1:12),'(es12.3)') CS%CFL_trunc
3833	0	call MOM_error(NOTE, "MOM_vert_friction: updateCFLtruncationValue set CFL limit to "//trim(msg))
3834		end subroutine updateCFLtruncationValue
3835
3836		!> Clean up and deallocate the vertical friction module
3837	1	subroutine vertvisc_end(CS)
3838		type(vertvisc_CS), intent(inout) :: CS !< Vertical viscosity control structure that
3839		!! will be deallocated in this subroutine.
3840
3841	1	if ((len_trim(CS%u_trunc_file) > 0) .or. (len_trim(CS%v_trunc_file) > 0)) &
3842	0	deallocate(CS%PointAccel_CSp)
3843
3844	1	DEALLOC_(CS%a_u) ; DEALLOC_(CS%h_u)
3845	1	DEALLOC_(CS%a_v) ; DEALLOC_(CS%h_v)
3846	1	if (associated(CS%a1_shelf_u)) deallocate(CS%a1_shelf_u)
3847	1	if (associated(CS%a1_shelf_v)) deallocate(CS%a1_shelf_v)
3848	1	if (allocated(CS%kappa_gl90_2d)) deallocate(CS%kappa_gl90_2d)
3849	1	end subroutine vertvisc_end
3850
3851		!> \namespace mom_vert_friction
3852		!! \author Robert Hallberg
3853		!! \date April 1994 - October 2006
3854		!!
3855		!! The vertical diffusion of momentum is fully implicit. This is
3856		!! necessary to allow for vanishingly small layers. The coupling
3857		!! is based on the distance between the centers of adjacent layers,
3858		!! except where a layer is close to the bottom compared with a
3859		!! bottom boundary layer thickness when a bottom drag law is used.
3860		!! A stress top b.c. and a no slip bottom b.c. are used. There
3861		!! is no limit on the time step for vertvisc.
3862		!!
3863		!! Near the bottom, the horizontal thickness interpolation scheme
3864		!! changes to an upwind biased estimate to control the effect of
3865		!! spurious Montgomery potential gradients at the bottom where
3866		!! nearly massless layers layers ride over the topography. Within a
3867		!! few boundary layer depths of the bottom, the harmonic mean
3868		!! thickness (i.e. (2 h+ h-) / (h+ + h-) ) is used if the velocity
3869		!! is from the thinner side and the arithmetic mean thickness
3870		!! (i.e. (h+ + h-)/2) is used if the velocity is from the thicker
3871		!! side. Both of these thickness estimates are second order
3872		!! accurate. Above this the arithmetic mean thickness is used.
3873		!!
3874		!! In addition, vertvisc truncates any velocity component that exceeds a
3875		!! maximum CFL number to a fraction of this value. This basically keeps
3876		!! instabilities spatially localized. The number of times the velocity is
3877		!! truncated is reported each time the energies are saved, and if
3878		!! exceeds CS%Maxtrunc the model will stop itself and change the time
3879		!! to a large value. This has proven very useful in (1) diagnosing
3880		!! model failures and (2) letting the model settle down to a
3881		!! meaningful integration from a poorly specified initial condition.
3882		!!
3883		!! The same code is used for the two velocity components, by
3884		!! indirectly referencing the velocities and defining a handful of
3885		!! direction-specific defined variables.
3886		!!
3887		!! Macros written all in capital letters are defined in MOM_memory.h.
3888		!!
3889		!! A small fragment of the grid is shown below:
3890		!! \verbatim
3891		!! j+1 x ^ x ^ x At x: q
3892		!! j+1 > o > o > At ^: v, frhatv, tauy
3893		!! j x ^ x ^ x At >: u, frhatu, taux
3894		!! j > o > o > At o: h
3895		!! j-1 x ^ x ^ x
3896		!! i-1 i i+1 At x & ^:
3897		!! i i+1 At > & o:
3898		!! \endverbatim
3899		!!
3900		!! The boundaries always run through q grid points (x).
3901	0	end module MOM_vert_friction