src/parameterizations/vertical/MOM_set_viscosity.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		!> Calculates various values related to the bottom boundary layer, such as the viscosity and
8		!! thickness of the BBL (set_viscous_BBL).
9		module MOM_set_visc
10
11		use MOM_ALE, only : ALE_CS, ALE_remap_velocities, ALE_remap_interface_vals, ALE_remap_vertex_vals
12		use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, CLOCK_ROUTINE
13		use MOM_cvmix_conv, only : cvmix_conv_is_used
14		use MOM_CVMix_ddiff, only : CVMix_ddiff_is_used
15		use MOM_cvmix_shear, only : cvmix_shear_is_used
16		use MOM_debugging, only : uvchksum, hchksum
17		use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
18		use MOM_diag_mediator, only : diag_ctrl, time_type
19		use MOM_domains, only : pass_var, CORNER
20		use MOM_EOS, only : calculate_density, calculate_density_derivs, calculate_specific_vol_derivs
21		use MOM_error_handler, only : MOM_error, FATAL, WARNING
22		use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
23		use MOM_file_parser, only : openParameterBlock, closeParameterBlock
24		use MOM_forcing_type, only : forcing, mech_forcing, find_ustar
25		use MOM_grid, only : ocean_grid_type
26		use MOM_hor_index, only : hor_index_type
27		use MOM_interface_heights, only : thickness_to_dz
28		use MOM_intrinsic_functions, only : cuberoot
29		use MOM_io, only : slasher, MOM_read_data, vardesc, var_desc
30		use MOM_kappa_shear, only : kappa_shear_is_used, kappa_shear_at_vertex
31		use MOM_open_boundary, only : ocean_OBC_type, OBC_segment_type, OBC_NONE, OBC_DIRECTION_E
32		use MOM_open_boundary, only : OBC_DIRECTION_W, OBC_DIRECTION_N, OBC_DIRECTION_S
33		use MOM_restart, only : register_restart_field, query_initialized, MOM_restart_CS
34		use MOM_restart, only : register_restart_field_as_obsolete, register_restart_pair
35		use MOM_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc
36		use MOM_unit_scaling, only : unit_scale_type
37		use MOM_variables, only : thermo_var_ptrs, vertvisc_type, porous_barrier_type
38		use MOM_verticalGrid, only : verticalGrid_type, get_thickness_units
39
40		implicit none ; private
41
42		#include <MOM_memory.h>
43
44		public set_viscous_BBL, set_viscous_ML, set_visc_init, set_visc_end
45		public set_visc_register_restarts, set_u_at_v, set_v_at_u
46		public remap_vertvisc_aux_vars
47
48		! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
49		! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
50		! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
51		! vary with the Boussinesq approximation, the Boussinesq variant is given first.
52
53		!> Control structure for MOM_set_visc
54		type, public :: set_visc_CS ; private
55		logical :: initialized = .false. !< True if this control structure has been initialized.
56		real :: Hbbl !< The static bottom boundary layer thickness [H ~> m or kg m-2].
57		!! Runtime parameter `HBBL`.
58		real :: dz_bbl !< The static bottom boundary layer thickness in height units [Z ~> m].
59		!! Runtime parameter `HBBL`.
60		real :: cdrag !< The quadratic drag coefficient [nondim].
61		!! Runtime parameter `CDRAG`.
62		real :: c_Smag !< The Laplacian Smagorinsky coefficient for
63		!! calculating the drag in channels [nondim].
64		real :: drag_bg_vel !< An assumed unresolved background velocity for
65		!! calculating the bottom drag [L T-1 ~> m s-1].
66		!! Runtime parameter `DRAG_BG_VEL`.
67		!! Should not be used if BBL_USE_TIDAL_BG is True.
68		real :: BBL_thick_min !< The minimum bottom boundary layer thickness [Z ~> m].
69		!! This might be Kv / (cdrag * drag_bg_vel) to give
70		!! Kv as the minimum near-bottom viscosity.
71		real :: Htbl_shelf !< A nominal thickness of the surface boundary layer for use
72		!! in calculating the near-surface velocity [H ~> m or kg m-2].
73		real :: Htbl_shelf_min !< The minimum surface boundary layer thickness [Z ~> m].
74		real :: KV_BBL_min !< The minimum viscosity in the bottom boundary layer [H Z T-1 ~> m2 s-1 or Pa s]
75		real :: KV_TBL_min !< The minimum viscosity in the top boundary layer [H Z T-1 ~> m2 s-1 or Pa s]
76		logical :: bottomdraglaw !< If true, the bottom stress is calculated with a
77		!! drag law c_drag\|u\|u. The velocity magnitude
78		!! may be an assumed value or it may be based on the
79		!! actual velocity in the bottommost `HBBL`, depending
80		!! on whether linear_drag is true.
81		!! Runtime parameter `BOTTOMDRAGLAW`.
82		logical :: bottomdragmap !< If true, apply the spatially varying drag coefficient (cdrag_2d)
83		!! instead of the spatially uniform drag coefficient (cdrag).
84		logical :: body_force_drag !< If true, the bottom stress is imposed as an explicit body force
85		!! applied over a fixed distance from the bottom, rather than as an
86		!! implicit calculation based on an enhanced near-bottom viscosity.
87		logical :: BBL_use_EOS !< If true, use the equation of state in determining
88		!! the properties of the bottom boundary layer.
89		logical :: linear_drag !< If true, the drag law is cdrag`DRAG_BG_VEL`u.
90		!! Runtime parameter `LINEAR_DRAG`.
91		logical :: Channel_drag !< If true, the drag is exerted directly on each layer
92		!! according to what fraction of the bottom they overlie.
93		real :: Chan_drag_max_vol !< The maximum bottom boundary layer volume within which the
94		!! channel drag is applied, normalized by the full cell area,
95		!! or a negative value to apply no maximum [Z ~> m].
96		real :: channel_break_depth !< When CHANNEL_DRAG is true, the bathymetric depth interpolated
97		!! to the vorticity point is a combination of the harmonic mean of the
98		!! adjacent velocity point depths below this depth [Z ~> m] and the
99		!! arithmetic mean of the adjacent depths above it, to roughly mimic a
100		!! continental shelf break profile. The internal version of this depth
101		!! uses the same offset (G%Z_ref) as the bathymetry.
102		logical :: correct_BBL_bounds !< If true, uses the correct bounds on the BBL thickness and
103		!! viscosity so that the bottom layer feels the intended drag.
104		logical :: RiNo_mix !< If true, use Richardson number dependent mixing.
105		logical :: dynamic_viscous_ML !< If true, use a bulk Richardson number criterion to
106		!! determine the mixed layer thickness for viscosity.
107		real :: bulk_Ri_ML !< The bulk mixed layer used to determine the
108		!! thickness of the viscous mixed layer [nondim]
109		real :: omega !< The Earth's rotation rate [T-1 ~> s-1].
110		real :: ustar_min !< A minimum value of ustar to avoid numerical
111		!! problems [H T-1 ~> m s-1 or kg m-2 s-1]. If the value is
112		!! small enough, this should not affect the solution.
113		real :: TKE_decay !< The ratio of the natural Ekman depth to the TKE
114		!! decay scale [nondim]
115		real :: omega_frac !< When setting the decay scale for turbulence, use this
116		!! fraction of the absolute rotation rate blended with the local
117		!! value of f, as sqrt((1-of)f^2 + of4*omega^2) [nondim]
118		real :: tideampfac2 !< A factor to multiply by tideamp to convert to a mean ustar,
119		!! accounts for conversion of amplitude to mean magnitude over
120		!! a time average much longer than the tidal periods and for
121		!! non-commuting conversion of mean tideamp to mean ustar**3 [nondim]
122		logical :: concave_trigonometric_L !< If true, use trigonometric expressions to determine the
123		!! fractional open interface lengths for concave topography.
124		integer :: answer_date !< The vintage of the order of arithmetic and expressions in the set
125		!! viscosity calculations. Values below 20190101 recover the answers
126		!! from the end of 2018, while higher values use updated and more robust
127		!! forms of the same expressions.
128		logical :: debug !< If true, write verbose checksums for debugging purposes.
129		logical :: BBL_use_tidal_bg !< If true, use a tidal background amplitude for the bottom velocity
130		!! when computing the bottom stress.
131		character(len=200) :: inputdir !< The directory for input files.
132		type(ocean_OBC_type), pointer :: OBC => NULL() !< Open boundaries control structure
133		type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
134		!! regulate the timing of diagnostic output.
135		! Allocatable data arrays
136		real, allocatable, dimension(:,:) :: cdrag_u !< The spatially varying quadratic drag coefficient [nondim]
137		real, allocatable, dimension(:,:) :: cdrag_v !< The spatially varying quadratic drag coefficient [nondim]
138		real, allocatable, dimension(:,:) :: tideamp !< RMS tidal amplitude at h points [Z T-1 ~> m s-1]
139		! Diagnostic arrays
140		real, allocatable, dimension(:,:) :: bbl_u !< BBL mean U current [L T-1 ~> m s-1]
141		real, allocatable, dimension(:,:) :: bbl_v !< BBL mean V current [L T-1 ~> m s-1]
142		!>@{ Diagnostics handles
143		integer :: id_bbl_thick_u = -1, id_kv_bbl_u = -1, id_bbl_u = -1
144		integer :: id_bbl_thick_v = -1, id_kv_bbl_v = -1, id_bbl_v = -1
145		integer :: id_Ray_u = -1, id_Ray_v = -1
146		integer :: id_nkml_visc_u = -1, id_nkml_visc_v = -1
147		!>@}
148		end type set_visc_CS
149
150		contains
151
152		!> Calculates the thickness of the bottom boundary layer and the viscosity within that layer.
153	12	subroutine set_viscous_BBL(u, v, h, tv, visc, G, GV, US, CS, pbv)
154		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
155		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
156		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
157		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
158		intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
159		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
160		intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
161		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
162		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
163		type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
164		!! available thermodynamic fields. Absent fields
165		!! have NULL pointers.
166		type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
167		!! related fields.
168		type(set_visc_CS), intent(inout) :: CS !< The control structure returned by a previous
169		!! call to set_visc_init.
170		type(porous_barrier_type),intent(in) :: pbv !< porous barrier fractional cell metrics
171
172		! Local variables
173		real, dimension(SZIB_(G),SZJB_(G)) :: &
174	24	ustar, & ! The bottom friction velocity [H T-1 ~> m s-1 or kg m-2 s-1].
175	24	T_EOS, & ! The temperature used to calculate the partial derivatives
176		! of density with T and S [C ~> degC].
177	24	S_EOS, & ! The salinity used to calculate the partial derivatives
178		! of density with T and S [S ~> ppt].
179	24	dR_dT, & ! Partial derivative of the density in the bottom boundary
180		! layer with temperature [R C-1 ~> kg m-3 degC-1].
181	24	dR_dS, & ! Partial derivative of the density in the bottom boundary
182		! layer with salinity [R S-1 ~> kg m-3 ppt-1].
183	24	press, & ! The pressure at which dR_dT and dR_dS are evaluated [R L2 T-2 ~> Pa].
184	24	umag_avg, & ! The average magnitude of velocities in the bottom boundary layer [L T-1 ~> m s-1].
185	24	h_bbl_drag, & ! The thickness over which to apply drag as a body force [H ~> m or kg m-2].
186	24	dz_bbl_drag ! The vertical height over which to apply drag as a body force [Z ~> m].
187		real :: htot ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
188		real :: dztot ! Distance from the bottom up to some point [Z ~> m].
189		real :: htot_vel ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
190		real :: dztot_vel ! Distance from the bottom up to some point [Z ~> m].
191
192		real :: Rhtot ! Running sum of thicknesses times the layer potential
193		! densities [H R ~> kg m-2 or kg2 m-5].
194		real, dimension(SZIB_(G),SZJ_(G)) :: &
195	12	D_u, & ! Bottom depth linearly interpolated to u points [Z ~> m].
196	24	mask_u ! A mask that disables any contributions from u points that
197		! are land or past open boundary conditions [nondim], 0 or 1.
198		real, dimension(SZI_(G),SZJB_(G)) :: &
199	24	D_v, & ! Bottom depth linearly interpolated to v points [Z ~> m].
200	24	mask_v ! A mask that disables any contributions from v points that
201		! are land or past open boundary conditions [nondim], 0 or 1.
202		real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)) :: &
203	24	h_at_vel, & ! Layer thickness at a velocity point, using an upwind-biased
204		! second order accurate estimate based on the previous velocity
205		! direction [H ~> m or kg m-2].
206	24	h_vel, & ! Arithmetic mean of the layer thicknesses adjacent to a
207		! velocity point [H ~> m or kg m-2].
208	24	dz_at_vel, & ! Vertical extent of a layer, using an upwind-biased
209		! second order accurate estimate based on the previous velocity
210		! direction [Z ~> m].
211	24	dz_vel, & ! Arithmetic mean of the difference in across the layers adjacent
212		! to a velocity point [Z ~> m].
213	24	T_vel, & ! Arithmetic mean of the layer temperatures adjacent to a
214		! velocity point [C ~> degC].
215	24	S_vel, & ! Arithmetic mean of the layer salinities adjacent to a
216		! velocity point [S ~> ppt].
217	24	SpV_vel, & ! Arithmetic mean of the layer averaged specific volumes adjacent to a
218		! velocity point [R-1 ~> m3 kg-1].
219	24	Rml_vel ! Arithmetic mean of the layer coordinate densities adjacent
220		! to a velocity point [R ~> kg m-3].
221	24	real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! Height change across layers [Z ~> m]
222
223		real :: h_vel_pos ! The arithmetic mean thickness at a velocity point
224		! plus H_neglect to avoid 0 values [H ~> m or kg m-2].
225		real :: ustarsq ! 400 times the square of ustar, times
226		! Rho0 divided by G_Earth and the conversion
227		! from m to thickness units [H R ~> kg m-2 or kg2 m-5].
228		real :: cdrag ! The drag coefficient [nondim].
229		real :: cdrag_sqrt ! Square root of the drag coefficient [nondim].
230		real :: cdrag_sqrt_H ! Square root of the drag coefficient, times a unit conversion factor
231		! from lateral lengths to layer thicknesses [H L-1 ~> nondim or kg m-3].
232		real :: cdrag_sqrt_H_RL ! Square root of the drag coefficient, times a unit conversion factor from
233		! density times lateral lengths to layer thicknesses [H L-1 R-1 ~> m3 kg-1 or nondim]
234		real :: cdrag_L_to_H ! The drag coefficient times conversion factors from lateral
235		! distance to thickness units [H L-1 ~> nondim or kg m-3]
236		real :: cdrag_RL_to_H ! The drag coefficient times conversion factors from density times lateral
237		! distance to thickness units [H L-1 R-1 ~> m3 kg-1 or nondim]
238		real :: cdrag_conv ! The drag coefficient times a combination of static conversion factors and in
239		! situ density or Boussinesq reference density [H L-1 ~> nondim or kg m-3]
240		real :: oldfn ! The integrated energy required to
241		! entrain up to the bottom of the layer,
242		! divided by G_Earth [H R ~> kg m-2 or kg2 m-5].
243		real :: Dfn ! The increment in oldfn for entraining
244		! the layer [H R ~> kg m-2 or kg2 m-5].
245		real :: frac_used ! The fraction of the present layer that contributes to Dh and Ddz [nondim]
246		real :: Dh ! The increment in layer thickness from
247		! the present layer [H ~> m or kg m-2].
248		real :: Ddz ! The increment in height change from the present layer [Z ~> m].
249		real :: bbl_thick ! The thickness of the bottom boundary layer [Z ~> m].
250		real :: BBL_thick_max ! A huge upper bound on the boundary layer thickness [Z ~> m].
251		real :: kv_bbl ! The bottom boundary layer viscosity [H Z T-1 ~> m2 s-1 or Pa s]
252		real :: C2f ! C2f = 2*f at velocity points [T-1 ~> s-1].
253	24	real :: u2_bg(SZIB_(G),SZJB_(G)) ! The square of an assumed background velocity, for calculating the mean
254		! magnitude near the bottom for use in the quadratic bottom drag [L2 T-2 ~> m2 s-2].
255		real :: hwtot ! Sum of the thicknesses used to calculate
256		! the near-bottom velocity magnitude [H ~> m or kg m-2].
257		real :: I_hwtot ! The Adcroft reciprocal of hwtot [H-1 ~> m-1 or m2 kg-1].
258		real :: dzwtot ! The vertical extent of the region used to calculate
259		! the near-bottom velocity magnitude [Z ~> m].
260		real :: hutot ! Running sum of thicknesses times the velocity
261		! magnitudes [H L T-1 ~> m2 s-1 or kg m-1 s-1].
262		real :: Thtot ! Running sum of thickness times temperature [C H ~> degC m or degC kg m-2].
263		real :: Shtot ! Running sum of thickness times salinity [S H ~> ppt m or ppt kg m-2].
264		real :: SpV_htot ! Running sum of thickness times specific volume [H R-1 ~> m4 kg-1 or m]
265		real :: hweight ! The thickness of a layer that is within Hbbl
266		! of the bottom [H ~> m or kg m-2].
267		real :: dzweight ! The counterpart of hweight in height units [Z ~> m].
268		real :: v_at_u, u_at_v ! v at a u point or vice versa [L T-1 ~> m s-1].
269		real :: Rho0x400_G ! 400*Rho0/G_Earth, times unit conversion factors
270		! [R T2 H-1 ~> kg s2 m-4 or s2 m-1].
271		! The 400 is a constant proposed by Killworth and Edwards, 1999.
272		real, dimension(SZI_(G),SZJ_(G),max(GV%nk_rho_varies,1)) :: &
273	24	Rml ! The mixed layer coordinate density [R ~> kg m-3].
274	24	real :: p_ref(SZI_(G),SZJ_(G)) ! The pressure used to calculate the coordinate
275		! density [R L2 T-2 ~> Pa] (usually set to 2e7 Pa = 2000 dbar).
276
277		real :: D_vel ! The bottom depth relative to the shelfbreak depth at a velocity point [Z ~> m].
278		real :: Dp, Dm ! The bottom depths at the edges of a velocity cell relative to the
279		! shelfbreak depth [Z ~> m].
280		real :: D_vel_p, D_vel_m ! The bottom depths in adjacent velocity points relative to the
281		! shelfbreak depth [Z ~> m].
282		real :: crv ! crv is the curvature of the bottom depth across a
283		! cell, times the cell width squared [Z ~> m].
284		real :: slope ! The absolute value of the bottom depth slope across
285		! a cell times the cell width [Z ~> m].
286		real :: Vol_bbl_chan ! The volume of the bottom boundary layer as used in the channel
287		! drag parameterization, normalized by the full horizontal area
288		! of the velocity cell [Z ~> m].
289	24	real :: vol_below(SZK_(GV)+1) ! The volume below each interface, normalized by the full
290		! horizontal area of a velocity cell [Z ~> m].
291	24	real :: L(SZK_(GV)+1) ! The fraction of the full cell width that is open at
292		! the depth of each interface [nondim].
293		! The next 9 variables are only used for debugging.
294	24	real :: L_trig(SZK_(GV)+1) ! The fraction of the full cell width that is open at
295		! the depth of each interface from trigonometric expressions [nondim].
296	24	real :: vol_err_trig(SZK_(GV)+1) ! The error in the volume below based on L_trig [Z ~> m]
297	24	real :: vol_err_iter(SZK_(GV)+1) ! The error in the volume below based on L_iter [Z ~> m]
298	24	real :: norm_err_trig(SZK_(GV)+1) ! vol_err_trig normalized by vol_below [nondim]
299	24	real :: norm_err_iter(SZK_(GV)+1) ! vol_err_iter normalized by vol_below [nondim]
300	24	real :: dL_trig_itt(SZK_(GV)+1) ! The difference between estimates of the fraction of the full cell
301		! width that is open at the depth of each interface [nondim].
302		real :: max_dL_trig_itt ! The largest difference between L and L_trig, for debugging [nondim]
303		real :: max_norm_err_trig ! The largest magnitude value of norm_err_trig in a column [nondim]
304		real :: max_norm_err_iter ! The largest magnitude value of norm_err_iter in a column [nondim]
305
306		real :: h_neglect ! A thickness that is so small it is usually lost
307		! in roundoff and can be neglected [H ~> m or kg m-2].
308		real :: dz_neglect ! A vertical distance that is so small it is usually lost
309		! in roundoff and can be neglected [Z ~> m].
310		real :: ustH ! ustar converted to units of H T-1 [H T-1 ~> m s-1 or kg m-2 s-1].
311		real :: root ! A temporary variable [H T-1 ~> m s-1 or kg m-2 s-1].
312
313		real :: Cell_width ! The transverse width of the velocity cell [L ~> m].
314		real :: Rayleigh ! A factor that is multiplied by the layer's velocity magnitude
315		! to give the Rayleigh drag velocity, times a lateral distance to
316		! thickness conversion factor [H L-1 ~> nondim or kg m-3].
317		real :: gam ! The ratio of the change in the open interface width
318		! to the open interface width atop a cell [nondim].
319		real :: BBL_frac ! The fraction of a layer's drag that goes into the
320		! viscous bottom boundary layer [nondim].
321		real :: BBL_visc_frac ! The fraction of all the drag that is expressed as
322		! a viscous bottom boundary layer [nondim].
323		real :: h_bbl_fr ! The fraction of the bottom boundary layer in a layer [nondim].
324		real :: h_sum ! The sum of the thicknesses of the layers below the one being
325		! worked on [H ~> m or kg m-2].
326		real :: tideampfac2_x_0p5 ! tideampfac2 multiplied by the c-grid averaging factor of 0.5
327		real, parameter :: C1_3 = 1.0/3.0, C1_6 = 1.0/6.0, C1_12 = 1.0/12.0 ! Rational constants [nondim]
328		real :: tmp ! A temporary variable, sometimes in [Z ~> m]
329	12	logical :: use_BBL_EOS, do_i(SZIB_(G),SZJB_(G))
330		integer, dimension(2,2) :: EOSdom ! The computational domain for the equation of state
331		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, m, n, K2, nkmb, nkml, jstart
332		integer :: is_OBC, ie_OBC, js_OBC, je_OBC
333		type(ocean_OBC_type), pointer :: OBC => NULL()
334
335	12	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
336	12	Isq = G%isc-1 ; Ieq = G%IecB ; Jsq = G%jsc-1 ; Jeq = G%JecB
337	12	nkmb = GV%nk_rho_varies ; nkml = GV%nkml
338	12	h_neglect = GV%H_subroundoff
339	12	dz_neglect = GV%dZ_subroundoff
340
341	12	Rho0x400_G = 400.0*(GV%H_to_RZ / GV%g_Earth_Z_T2)
342	12	tideampfac2_x_0p5 = CS%tideampfac2*0.5
343
344	12	if (.not.CS%initialized) call MOM_error(FATAL,"MOM_set_viscosity(BBL): "//&
345	0	"Module must be initialized before it is used.")
346
347	12	if (.not.CS%bottomdraglaw) return
348
349	12	if (CS%debug) then
350	0	call uvchksum("Start set_viscous_BBL [uv]", u, v, G%HI, haloshift=1, unscale=US%L_T_to_m_s)
351	0	call hchksum(h,"Start set_viscous_BBL h", G%HI, haloshift=1, unscale=GV%H_to_m)
352	0	if (associated(tv%T)) call hchksum(tv%T, "Start set_viscous_BBL T", G%HI, haloshift=1, unscale=US%C_to_degC)
353	0	if (associated(tv%S)) call hchksum(tv%S, "Start set_viscous_BBL S", G%HI, haloshift=1, unscale=US%S_to_ppt)
354	0	if (allocated(tv%SpV_avg)) &
355	0	call hchksum(tv%SpV_avg, "Start set_viscous_BBL SpV_avg", G%HI, haloshift=1, unscale=US%kg_m3_to_R)
356	0	if (allocated(tv%SpV_avg)) call hchksum(tv%SpV_avg, "Cornerless SpV_avg", G%HI, &
357	0	haloshift=1, omit_corners=.true., unscale=US%kg_m3_to_R)
358	0	if (associated(tv%T)) call hchksum(tv%T, "Cornerless T", G%HI, haloshift=1, &
359	0	omit_corners=.true., unscale=US%C_to_degC)
360	0	if (associated(tv%S)) call hchksum(tv%S, "Cornerless S", G%HI, haloshift=1, &
361	0	omit_corners=.true., unscale=US%S_to_ppt)
362		endif
363
364	12	use_BBL_EOS = associated(tv%eqn_of_state) .and. CS%BBL_use_EOS
365	12	OBC => CS%OBC
366
367	12	if (.not.CS%bottomdragmap) then
368	12	cdrag_sqrt = sqrt(CS%cdrag)
369	12	cdrag_sqrt_H = cdrag_sqrt * US%L_to_m * GV%m_to_H
370	12	cdrag_sqrt_H_RL = cdrag_sqrt * US%L_to_Z * GV%RZ_to_H
371	12	cdrag_L_to_H = CS%cdrag * US%L_to_m * GV%m_to_H
372	12	cdrag_RL_to_H = CS%cdrag * US%L_to_Z * GV%RZ_to_H
373		endif
374	12	BBL_thick_max = G%Rad_Earth_L * US%L_to_Z
375	12	K2 = max(nkmb+1, 2)
376
377		!$omp target enter data map(alloc: dz)
378
379		! Find the vertical distances across layers.
380	12	call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1, do_offload=.true.)
381
382		! With a linear drag law, the friction velocity is already known.
383		! if (CS%linear_drag) ustar(:) = cdrag_sqrt_H*CS%drag_bg_vel
384
385		!$omp target enter data map(to: tv, tv%T, tv%S, tv%p_surf, CS) map(alloc: Rml, p_ref, ustar, &
386		!$omp umag_avg, u2_bg, mask_u, mask_v, h_bbl_drag, dz_bbl_drag, do_i, dR_dS, dR_dT, D_u, D_v, &
387		!$omp press, S_EOS, T_EOS, Rml_vel)
388
389	12	if ((nkml>0) .and. .not.use_BBL_EOS) then
390	0	EOSdom(1,1) = Isq - (G%isd-1) ; EOSdom(1,2) = G%iec+1 - (G%isd-1)
391	0	EOSdom(2,1) = Jsq - (G%jsd-1) ; EOSdom(2,2) = G%jec+1 - (G%jsd-1)
392	0	do concurrent (j=Jsq:Jeq+1, i=Isq:Ieq+1)
393	0	p_ref(i,j) = tv%P_Ref
394		enddo
395	0	do k=1,nkmb
396	0	call calculate_density(tv%T(:,:,k), tv%S(:,:,k), p_ref, Rml(:,:,k), tv%eqn_of_state, EOSdom)
397		enddo
398		endif
399
400	12	do concurrent (J=js-1:je, i=is-1:ie+1)
401	89304	D_v(i,J) = 0.5*(G%bathyT(i,j) + G%bathyT(i,j+1))
402	90780	mask_v(i,J) = G%mask2dCv(i,J)
403		enddo
404	1464	do concurrent (j=js-1:je+1, I=is-1:ie)
405	90024	D_u(I,j) = 0.5*(G%bathyT(i,j) + G%bathyT(i+1,j))
406	91488	mask_u(I,j) = G%mask2dCu(I,j)
407		enddo
408
409	12	if (associated(OBC) .and. CS%Channel_drag) then
410		!$omp target update from(mask_u, mask_v, D_u, D_v)
411		! Use a one-sided projection of bottom depths at OBC points.
412	0	if (OBC%v_N_OBCs_on_PE) then
413	0	Js_OBC = max(js-1, OBC%Js_v_N_obc) ; Je_OBC = min(je, OBC%Je_v_N_obc)
414	0	is_OBC = max(is-1, OBC%is_v_N_obc) ; ie_OBC = min(ie+1, OBC%ie_v_N_obc)
415		!$OMP parallel do default(shared)
416	0	do J=Js_OBC,Je_OBC ; do i=is_OBC,ie_OBC
417	0	if (OBC%segnum_v(i,J) > 0) D_v(i,J) = G%bathyT(i,j) ! OBC_DIRECTION_N
418		enddo ; enddo
419		endif
420	0	if (OBC%v_S_OBCs_on_PE) then
421		!$omp target update from(D_v)
422	0	Js_OBC = max(js-1, OBC%Js_v_S_obc) ; Je_OBC = min(je, OBC%Je_v_S_obc)
423	0	is_OBC = max(is-1, OBC%is_v_S_obc) ; ie_OBC = min(ie+1, OBC%ie_v_S_obc)
424		!$OMP parallel do default(shared)
425	0	do J=Js_OBC,Je_OBC ; do i=is_OBC,ie_OBC
426	0	if (OBC%segnum_v(i,J) < 0) D_v(i,J) = G%bathyT(i,j+1) ! OBC_DIRECTION_S
427		enddo ; enddo
428		!$omp target update to(D_v)
429		endif
430	0	if (OBC%u_E_OBCs_on_PE) then
431		!$omp target update from(D_u)
432	0	js_OBC = max(js-1, OBC%js_u_E_obc) ; je_OBC = min(je+1, OBC%je_u_E_obc)
433	0	Is_OBC = max(is-1, OBC%Is_u_E_obc) ; Ie_OBC = min(ie, OBC%Ie_u_E_obc)
434		!$OMP parallel do default(shared)
435	0	do j=js_OBC,je_OBC ; do I=Is_OBC,Ie_OBC
436	0	if (OBC%segnum_u(I,j) > 0) D_u(I,j) = G%bathyT(i,j) ! OBC_DIRECTION_E
437		enddo ; enddo
438		!$omp target update to(D_u)
439		endif
440	0	if (OBC%u_W_OBCs_on_PE) then
441	0	js_OBC = max(js-1, OBC%js_u_W_obc) ; je_OBC = min(je+1, OBC%je_u_W_obc)
442	0	Is_OBC = max(is-1, OBC%Is_u_W_obc) ; Ie_OBC = min(ie, OBC%Ie_u_W_obc)
443		!$OMP parallel do default(shared)
444	0	do j=js_OBC,je_OBC ; do I=Is_OBC,Ie_OBC
445	0	if (OBC%segnum_u(I,j) < 0) D_u(I,j) = G%bathyT(i+1,j) ! OBC_DIRECTION_W
446		enddo ; enddo
447		endif
448
449	0	do n=1,OBC%number_of_segments
450		! Now project bottom depths across cell-corner points in the OBCs. The two
451		! projections have to occur in sequence and can not be combined easily.
452	0	if (.not. OBC%segment(n)%on_pe) cycle
453		! Use a one-sided projection of bottom depths at OBC points.
454	0	I = OBC%segment(n)%HI%IsdB ; J = OBC%segment(n)%HI%JsdB
455	0	if (OBC%segment(n)%is_N_or_S .and. (J >= js-1) .and. (J <= je)) then
456	0	do I = max(is-1,OBC%segment(n)%HI%IsdB), min(ie,OBC%segment(n)%HI%IedB)
457	0	if (OBC%segment(n)%direction == OBC_DIRECTION_N) then
458	0	D_u(I,j+1) = D_u(I,j) ; mask_u(I,j+1) = 0.0
459	0	elseif (OBC%segment(n)%direction == OBC_DIRECTION_S) then
460	0	D_u(I,j) = D_u(I,j+1) ; mask_u(I,j) = 0.0
461		endif
462		enddo
463	0	elseif (OBC%segment(n)%is_E_or_W .and. (I >= is-1) .and. (I <= ie)) then
464	0	do J = max(js-1,OBC%segment(n)%HI%JsdB), min(je,OBC%segment(n)%HI%JedB)
465	0	if (OBC%segment(n)%direction == OBC_DIRECTION_E) then
466	0	D_v(i+1,J) = D_v(i,J) ; mask_v(i+1,J) = 0.0
467	0	elseif (OBC%segment(n)%direction == OBC_DIRECTION_W) then
468	0	D_v(i,J) = D_v(i+1,J) ; mask_v(i,J) = 0.0
469		endif
470		enddo
471		endif
472		enddo
473		!$omp target update to(mask_u, mask_v, D_u, D_v)
474		endif
475
476	12	if (.not.use_BBL_EOS) then
477	0	do concurrent (k=1:nz, j=G%jsdB:G%Jedb, i=G%isdB:G%iedB)
478	0	Rml_vel(i,j,k) = 0.0
479		enddo
480		endif
481
482		! Resetting Ray_[uv] is required by body force drag.
483	12	if (allocated(visc%Ray_u)) then
484	1560	do concurrent (k=1:nz, j=G%jsd:G%jed, i=G%isdB:G%iedB)
485	8001624	visc%Ray_u(i,j,k) = 0.0
486		enddo
487		endif
488	12	if (allocated(visc%Ray_v)) then
489	1548	do concurrent (k=1:nz, j=G%jsdB:G%jedB, i=G%isd:G%ied)
490	8056332	visc%Ray_v(i,j,k) = 0.0
491		enddo
492		endif
493
494		!$omp target enter data map(alloc: S_vel, T_vel, SpV_vel, h_vel, h_at_vel, dz_vel, &
495		!$omp dz_at_vel)
496
497	36	do m=1,2
498	24	if (m==1) then
499		! m=1 refers to u-points
500	12	is = Isq ; ie = Ieq
501	12	jstart = G%Jsc
502		else
503		! m=2 refers to v-points
504	12	is = G%isc ; ie = G%iec
505	12	jstart = Jsq
506		endif
507
508	24	do concurrent (j=jstart:Jeq)
509	1452	if (m==1) then
510	720	do concurrent (i=is:ie)
511	87840	do_i(i,j) = (G%mask2dCu(I,j) > 0.0)
512		enddo
513		else
514	732	do concurrent (i=is:ie)
515	88572	do_i(i,j) = (G%mask2dCv(i,J) > 0.0)
516		enddo
517		endif
518
519		! Calculate thickness at velocity points (u or v depending on value of m).
520		! Also interpolate the ML density or T/S properties.
521	1452	if (m==1) then ! u-points
522	720	do concurrent (k=1:nz, I=is:ie)
523	6534000	if (do_i(I,j)) then
524	4453200	if (u(I,j,k) * (h(i+1,j,k) - h(i,j,k)) >= 0) then
525		! If the flow is from thin to thick then bias towards the thinner thickness
526		h_at_vel(I,j,k) = 2.0h(i,j,k)h(i+1,j,k) / &
527	2551229	(h(i,j,k) + h(i+1,j,k) + h_neglect)
528		dz_at_vel(I,j,k) = 2.0dz(i,j,k)dz(i+1,j,k) / &
529	2551229	(dz(i,j,k) + dz(i+1,j,k) + dz_neglect)
530		else
531		! If the flow is from thick to thin then use the simple average thickness
532	1901971	h_at_vel(I,j,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
533	1901971	dz_at_vel(I,j,k) = 0.5 * (dz(i,j,k) + dz(i+1,j,k))
534		endif
535		endif
536	6534000	h_vel(I,j,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
537	6621840	dz_vel(I,j,k) = 0.5 * (dz(i,j,k) + dz(i+1,j,k))
538		enddo
539	720	if (use_BBL_EOS) then ; do concurrent (k=1:nz, I=is:ie)
540		! Perhaps these should be thickness weighted.
541	6534000	T_vel(I,j,k) = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
542	6621840	S_vel(I,j,k) = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
543	0	enddo ; else ; do concurrent (k=1:nkmb, I=is:ie)
544	0	Rml_vel(I,j,k) = 0.5 * (Rml(i,j,k) + Rml(i+1,j,k))
545		enddo ; endif
546	720	if (allocated(tv%SpV_avg)) then ; do concurrent (k=1:nz, I=is:ie)
547	0	SpV_vel(I,j,k) = 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i+1,j,k))
548		enddo ; endif
549		else ! v-points
550	732	do concurrent (k=1:nz, i=is:ie)
551	6588000	if (do_i(i,j)) then
552	4404600	if (v(i,J,k) * (h(i,j+1,k) - h(i,j,k)) >= 0) then
553		! If the flow is from thin to thick then bias towards the thinner thickness
554		h_at_vel(i,j,k) = 2.0h(i,j,k)h(i,j+1,k) / &
555	2337718	(h(i,j,k) + h(i,j+1,k) + h_neglect)
556		dz_at_vel(i,j,k) = 2.0dz(i,j,k)dz(i,j+1,k) / &
557	2337718	(dz(i,j,k) + dz(i,j+1,k) + dz_neglect)
558		else
559		! If the flow is from thick to thin then use the simple average thickness
560	2066882	h_at_vel(i,j,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
561	2066882	dz_at_vel(i,j,k) = 0.5 * (dz(i,j,k) + dz(i,j+1,k))
562		endif
563		endif
564	6588000	h_vel(i,j,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
565	6676572	dz_vel(i,j,k) = 0.5 * (dz(i,j,k) + dz(i,j+1,k))
566		enddo
567	732	if (use_BBL_EOS) then ; do concurrent (k=1:nz, i=is:ie)
568		! Perhaps these should be thickness weighted.
569	6588000	T_vel(i,j,k) = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
570	6676572	S_vel(i,j,k) = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
571	0	enddo ; else ; do concurrent (k=1:nkmb, i=is:ie)
572	0	Rml_vel(i,j,k) = 0.5 * (Rml(i,j,k) + Rml(i,j+1,k))
573		enddo ; endif
574	732	if (allocated(tv%SpV_avg)) then ; do concurrent (k=1:nz, i=is:ie)
575	0	SpV_vel(i,j,k) = 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j+1,k))
576		enddo ; endif
577		endif
578
579	1452	if (associated(OBC)) then ; if (OBC%number_of_segments > 0) then
580		! Apply a zero gradient projection of thickness across OBC points.
581	0	if (m==1) then
582		do concurrent (I=is:ie, do_i(I,j) .and. (OBC%segnum_u(I,j) /= 0)) &
583	0	DO_LOCALITY(local(k))
584	0	if (OBC%segnum_u(I,j) > 0) then ! OBC_DIRECTION_E
585	0	do k=1,nz
586	0	h_at_vel(I,j,k) = h(i,j,k) ; h_vel(I,j,k) = h(i,j,k)
587	0	dz_at_vel(I,j,k) = dz(i,j,k) ; dz_vel(I,j,k) = dz(i,j,k)
588		enddo
589	0	if (use_BBL_EOS) then
590	0	do k=1,nz
591	0	T_vel(I,j,k) = tv%T(i,j,k) ; S_vel(I,j,k) = tv%S(i,j,k)
592		enddo
593		else
594	0	do k=1,nkmb
595	0	Rml_vel(I,j,k) = Rml(i,j,k)
596		enddo
597		endif
598	0	if (allocated(tv%SpV_avg)) then ; do k=1,nz
599	0	SpV_vel(I,j,k) = tv%SpV_avg(i,j,k)
600		enddo ; endif
601	0	elseif (OBC%segnum_u(I,j) < 0) then ! OBC_DIRECTION_W
602	0	do k=1,nz
603	0	h_at_vel(I,j,k) = h(i+1,j,k) ; h_vel(I,j,k) = h(i+1,j,k)
604	0	dz_at_vel(I,j,k) = dz(i+1,j,k) ; dz_vel(I,j,k) = dz(i+1,j,k)
605		enddo
606	0	if (use_BBL_EOS) then
607	0	do k=1,nz
608	0	T_vel(I,j,k) = tv%T(i+1,j,k) ; S_vel(I,j,k) = tv%S(i+1,j,k)
609		enddo
610		else
611	0	do k=1,nkmb
612	0	Rml_vel(I,j,k) = Rml(i+1,j,k)
613		enddo
614		endif
615	0	if (allocated(tv%SpV_avg)) then ; do k=1,nz
616	0	SpV_vel(I,j,k) = tv%SpV_avg(i+1,j,k)
617		enddo ; endif
618		endif
619		enddo
620		else
621		do concurrent (i=is:ie, do_i(i,j) .and. (OBC%segnum_v(i,J) /= 0)) &
622	0	DO_LOCALITY(local(k))
623	0	if (OBC%segnum_v(i,J) > 0) then ! OBC_DIRECTION_N
624	0	do k=1,nz
625	0	h_at_vel(i,j,k) = h(i,j,k) ; h_vel(i,j,k) = h(i,j,k)
626	0	dz_at_vel(i,j,k) = dz(i,j,k) ; dz_vel(i,j,k) = dz(i,j,k)
627		enddo
628	0	if (use_BBL_EOS) then
629	0	do k=1,nz
630	0	T_vel(i,j,k) = tv%T(i,j,k) ; S_vel(i,j,k) = tv%S(i,j,k)
631		enddo
632		else
633	0	do k=1,nkmb
634	0	Rml_vel(i,j,k) = Rml(i,j,k)
635		enddo
636		endif
637	0	if (allocated(tv%SpV_avg)) then ; do k=1,nz
638	0	SpV_vel(i,j,k) = tv%SpV_avg(i,j,k)
639		enddo ; endif
640	0	elseif (OBC%segnum_v(i,J) < 0) then ! OBC_DIRECTION_S
641	0	do k=1,nz
642	0	h_at_vel(i,j,k) = h(i,j+1,k) ; h_vel(i,j,k) = h(i,j+1,k)
643	0	dz_at_vel(i,j,k) = dz(i,j+1,k) ; dz_vel(i,j,k) = dz(i,j+1,k)
644		enddo
645	0	if (use_BBL_EOS) then
646	0	do k=1,nz
647	0	T_vel(i,j,k) = tv%T(i,j+1,k) ; S_vel(i,j,k) = tv%S(i,j+1,k)
648		enddo
649		else
650	0	do k=1,nkmb
651	0	Rml_vel(i,j,k) = Rml(i,j+1,k)
652		enddo
653		endif
654	0	if (allocated(tv%SpV_avg)) then ; do k=1,nz
655	0	SpV_vel(i,j,k) = tv%SpV_avg(i,j+1,k)
656		enddo ; endif
657		endif
658		enddo
659		endif
660		endif ; endif
661
662		! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant
663	1452	if (CS%BBL_use_tidal_bg) then
664	0	do concurrent (i=is:ie, do_i(i,j)) ; if (m==1) then
665		u2_bg(I,j) = tideampfac2_x_0p5 * ( G%mask2dT(i,j)(CS%tideamp(i,j)CS%tideamp(i,j))+ &
666	0	G%mask2dT(i+1,j)(CS%tideamp(i+1,j)CS%tideamp(i+1,j)) )
667		else
668		u2_bg(i,j) = tideampfac2_x_0p5 * ( G%mask2dT(i,j)(CS%tideamp(i,j)CS%tideamp(i,j))+ &
669	0	G%mask2dT(i,j+1)(CS%tideamp(i,j+1)CS%tideamp(i,j+1)) )
670		endif ; enddo
671		else
672	176412	do concurrent (i=is:ie, do_i(i,j))
673	176412	u2_bg(i,j) = CS%drag_bg_vel * CS%drag_bg_vel
674		enddo
675		endif
676
677	1452	if (use_BBL_EOS .or. CS%body_force_drag .or. .not.CS%linear_drag) then
678		! Calculate the mean velocity magnitude over the bottommost CS%Hbbl of
679		! the water column for determining the quadratic bottom drag.
680		! Used in ustar(i,j)
681		do concurrent (i=is:ie, do_i(i,j)) DO_LOCALITY(local(k, cdrag_sqrt)) &
682	176412	DO_LOCALITY(local_init(cdrag_sqrt_H, cdrag_sqrt_H_RL))
683	118104	htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
684	118104	dztot_vel = 0.0 ; dzwtot = 0.0
685	118104	Thtot = 0.0 ; Shtot = 0.0 ; SpV_htot = 0.0
686
687	118104	if (CS%bottomdragmap) then
688	0	if (m==1) then
689	0	cdrag_sqrt = sqrt(CS%cdrag_u(i,j))
690		else
691	0	cdrag_sqrt = sqrt(CS%cdrag_v(i,j))
692		endif
693	0	cdrag_sqrt_H = cdrag_sqrt * US%L_to_m * GV%m_to_H
694	0	cdrag_sqrt_H_RL = cdrag_sqrt * US%L_to_Z * GV%RZ_to_H
695		endif
696
697	2962762	do k=nz,1,-1
698	2962744	if (htot_vel>=CS%Hbbl) exit ! terminate the k loop
699
700	2844658	hweight = MIN(CS%Hbbl - htot_vel, h_at_vel(i,j,k))
701	2844658	if (hweight < 1.5*GV%Angstrom_H + h_neglect) cycle
702	2844658	dzweight = MIN(CS%dz_bbl - dztot_vel, dz_at_vel(i,j,k))
703
704	2844658	htot_vel = htot_vel + h_at_vel(i,j,k)
705	2844658	hwtot = hwtot + hweight
706	2844658	dztot_vel = dztot_vel + dz_at_vel(i,j,k)
707	2844658	dzwtot = dzwtot + dzweight
708
709	2844658	if ((.not.CS%linear_drag) .and. (hweight >= 0.0)) then ; if (m==1) then
710	1439521	v_at_u = set_v_at_u(v, h, G, GV, i, j, k, mask_v, OBC)
711	1439521	hutot = hutot + hweight * sqrt(u(I,j,k)u(I,j,k) + v_at_uv_at_u + u2_bg(I,j))
712		else
713	1405137	u_at_v = set_u_at_v(u, h, G, GV, i, j, k, mask_u, OBC)
714	1405137	hutot = hutot + hweight * sqrt(v(i,J,k)v(i,J,k) + u_at_vu_at_v + u2_bg(i,j))
715		endif ; endif
716
717	2844658	if (use_BBL_EOS .and. (hweight >= 0.0)) then
718	2844658	Thtot = Thtot + hweight * T_vel(i,j,k)
719	2844658	Shtot = Shtot + hweight * S_vel(i,j,k)
720		endif
721	2844676	if (allocated(tv%SpV_avg) .and. (hweight >= 0.0)) then
722	0	SpV_htot = SpV_htot + hweight * SpV_vel(i,j,k)
723		endif
724		enddo ! end of k loop
725
726		! Find the Adcroft reciprocal of the total thickness weights
727	118104	I_hwtot = 0.0 ; if (hwtot > 0.0) I_hwtot = 1.0 / hwtot
728
729		! Set u* based on u*^2 = Cdrag u_bbl^2
730	118104	if ((hwtot <= 0.0) .or. (CS%linear_drag .and. .not.allocated(tv%SpV_avg))) then
731	0	ustar(i,j) = cdrag_sqrt_H * CS%drag_bg_vel
732	118104	elseif (CS%linear_drag .and. allocated(tv%SpV_avg)) then
733	0	ustar(i,j) = cdrag_sqrt_H_RL * CS%drag_bg_vel * (hwtot / SpV_htot)
734	118104	elseif (allocated(tv%SpV_avg)) then ! (.not.CS%linear_drag)
735	0	ustar(i,j) = cdrag_sqrt_H_RL * hutot / SpV_htot
736		else ! (.not.CS%linear_drag .and. .not.allocated(tv%SpV_avg))
737	118104	ustar(i,j) = cdrag_sqrt_H * hutot / hwtot
738		endif
739
740	118104	umag_avg(i,j) = hutot * I_hwtot
741	118104	h_bbl_drag(i,j) = hwtot
742	118104	dz_bbl_drag(i,j) = dzwtot
743
744	118104	if (use_BBL_EOS) then ; if (hwtot > 0.0) then
745	118104	T_EOS(i,j) = Thtot/hwtot ; S_EOS(i,j) = Shtot/hwtot
746		else
747	0	T_EOS(i,j) = 0.0 ; S_EOS(i,j) = 0.0
748		endif ; endif
749
750		! Diagnostic BBL flow speed at u- and v-points.
751	294516	if (CS%id_bbl_u>0 .and. m==1) then
752	0	if (hwtot > 0.0) CS%bbl_u(I,j) = hutot/hwtot
753	118104	elseif (CS%id_bbl_v>0 .and. m==2) then
754	0	if (hwtot > 0.0) CS%bbl_v(i,J) = hutot/hwtot
755		endif
756		enddo
757		else
758		do concurrent (i=is:ie) DO_LOCALITY(local(cdrag_sqrt)) &
759	0	DO_LOCALITY(local_init(cdrag_sqrt_H))
760	0	if (CS%bottomdragmap) then
761	0	if (m==1) then
762	0	cdrag_sqrt = sqrt(CS%cdrag_u(i,j))
763		else
764	0	cdrag_sqrt = sqrt(CS%cdrag_v(i,j))
765		endif
766	0	cdrag_sqrt_H = cdrag_sqrt * US%L_to_m * GV%m_to_H
767		endif
768	0	ustar(i,j) = cdrag_sqrt_H * CS%drag_bg_vel
769		enddo
770		endif ! Not linear_drag
771
772	2928	if (use_BBL_EOS) then
773	1452	if (associated(tv%p_surf)) then
774	88572	if (m==1) then ; do concurrent (i=is:ie) ; press(I,j) = 0.5*(tv%p_surf(i,j) + tv%p_surf(i+1,j)) ; enddo
775	88572	else ; do concurrent (i=is:ie) ; press(i,j) = 0.5*(tv%p_surf(i,j) + tv%p_surf(i,j+1)) ; enddo ; endif
776		else
777	0	do concurrent (i=is:ie) ; press(i,j) = 0.0 ; enddo
778		endif
779
780	351372	do concurrent (i=is:ie, .not.do_i(i,j)) ; T_EOS(i,j) = 0.0 ; S_EOS(i,j) = 0.0 ; enddo
781
782	1452	do concurrent (i=is:ie)
783	13298412	do k=1,nz
784	13296960	press(i,j) = press(i,j) + (GV%H_to_RZGV%g_Earth) h_vel(i,j,k)
785		enddo
786		enddo
787		endif
788		enddo ! end of j loop
789
790	24	if (use_BBL_EOS) then
791	24	EOSdom(1,1) = is-G%IsdB+1 ; EOSdom(1,2) = ie-G%IsdB+1
792	24	EOSdom(2,1) = jstart-G%JsdB+1 ; EOSdom(2,2) = Jeq-G%JsdB+1
793	24	call calculate_density_derivs(T_EOS, S_EOS, press, dR_dT, dR_dS, tv%eqn_of_state, EOSdom)
794		endif
795
796		! Find a BBL thickness given by equation 2.20 of Killworth and Edwards, 1999:
797		! ( f h / Cn u* )^2 + ( N h / Ci u* ) = 1
798		! where Cn=0.5 and Ci=20 (constants suggested by Zilitinkevich and Mironov, 1996).
799		! Eq. 2.20 can be expressed in terms of boundary layer thicknesses limited by
800		! rotation (h_f) and stratification (h_N):
801		! ( h / h_f )^2 + ( h / h_N ) = 1
802		! When stratification dominates h_N<<h_f, and vice versa.
803		!$omp target teams loop collapse(2) thread_limit(128) private( &
804		!$omp k, ustarsq, htot, dztot, Thtot, Shtot, oldfn, Dfn, Dh, Ddz, &
805		!$omp frac_used, Rhtot, C2f, ustH, bbl_thick, Vol_bbl_chan, vol_below, &
806		!$omp D_vel, D_vel_p, D_vel_m, Dp, tmp, Dm, crv, slope, &
807		!$omp max_dL_trig_itt, max_norm_err_trig, max_norm_err_iter, &
808		!$omp norm_err_trig, norm_err_iter, dL_trig_itt, vol_err_trig, L_trig, &
809		!$omp L, vol_err_iter, BBL_visc_frac, BBL_frac, cdrag_conv, h_vel_pos, &
810		!$omp Cell_width, gam, Rayleigh, v_at_u, u_at_v, kv_bbl, h_sum, &
811		!$omp I_hwtot, h_bbl_fr, root, cdrag &
812		!$omp ) firstprivate(cdrag_L_to_H, cdrag_RL_to_H)
813	176448	do j=jstart,jeq ; do i=is,ie ; if (do_i(i,j)) then
814		! The 400.0 in this expression is the square of a Ci introduced in KW99, eq. 2.22.
815	118104	ustarsq = Rho0x400_G * ustar(i,j)*2 ! Note not in units of u^2 but [H R ~> kg m-2 or kg2 m-5]
816	118104	htot = 0.0
817	118104	dztot = 0.0
818
819	118104	if (CS%bottomdragmap) then
820	0	if (m==1) then
821	0	cdrag = CS%cdrag_u(i,j)
822		else
823	0	cdrag = CS%cdrag_v(i,j)
824		endif
825	0	cdrag_L_to_H = cdrag * US%L_to_m * GV%m_to_H
826	0	cdrag_RL_to_H = cdrag * US%L_to_Z * GV%RZ_to_H
827		endif
828
829		! Calculate the thickness of a stratification limited BBL ignoring rotation:
830		! h_N = Ci u* / N (limit of KW99 eq. 2.20 for \|f\|->0)
831		! For layer mode, N^2 = g'/h. Since (Ci u*)^2 = (h_N N)^2 = h_N g' then
832		! h_N = (Ci u*)^2 / g' (KW99, eq, 2.22)
833		! Starting from the bottom, integrate the stratification upward until h_N N balances Ci u*
834		! or in layer mode
835		! h_N Delta rho ~ (Ci u*)^2 rho0 / g
836		! where the rhs is stored in variable ustarsq.
837		! The method was described in Stephens and Hallberg 2000 (unpublished and lost manuscript).
838	118104	if (use_BBL_EOS) then
839	118104	Thtot = 0.0 ; Shtot = 0.0 ; oldfn = 0.0
840	3048932	do k=nz,2,-1
841	3045350	if (h_at_vel(i,j,k) <= 0.0) cycle
842
843		! Delta rho * h_bbl assuming everything below is homogenized
844		oldfn = dR_dT(i,j)(Thtot - T_vel(i,j,k)htot) + &
845	3045350	dR_dS(i,j)(Shtot - S_vel(i,j,k)htot)
846	3045350	if (oldfn >= ustarsq) exit
847
848		! Local Delta rho * h_bbl at interface
849		Dfn = (dR_dT(i,j)*(T_vel(i,j,k) - T_vel(i,j,k-1)) + &
850		dR_dS(i,j)(S_vel(i,j,k) - S_vel(i,j,k-1))) &
851	2930828	(h_at_vel(i,j,k) + htot)
852
853	2930828	if ((oldfn + Dfn) <= ustarsq) then
854		! Use whole layer
855	2816266	Dh = h_at_vel(i,j,k)
856	2816266	Ddz = dz_at_vel(i,j,k)
857		else
858		! Use only part of the layer
859	114562	frac_used = sqrt((ustarsq-oldfn) / (Dfn))
860	114562	Dh = h_at_vel(i,j,k) * frac_used
861	114562	Ddz = dz_at_vel(i,j,k) * frac_used
862		endif
863
864		! Increment total BBL thickness and cumulative T and S
865	2930828	htot = htot + Dh
866	2930828	dztot = dztot + Ddz
867	2934410	Thtot = Thtot + T_vel(i,j,k)Dh ; Shtot = Shtot + S_vel(i,j,k)Dh
868		enddo
869	118104	if ((oldfn < ustarsq) .and. h_at_vel(i,j,1) > 0.0) then
870		! Layer 1 might be part of the BBL.
871	3582	if (dR_dT(i,j) * (Thtot - T_vel(i,j,1)*htot) + &
872		dR_dS(i,j) * (Shtot - S_vel(i,j,1)*htot) < ustarsq) then
873	3542	htot = htot + h_at_vel(i,j,1)
874	3542	dztot = dztot + dz_at_vel(i,j,1)
875		endif
876		endif ! Examination of layer 1.
877		else ! Use Rlay and/or the coordinate density as density variables.
878	0	Rhtot = 0.0
879	0	do k=nz,K2,-1
880	0	oldfn = Rhtot - GV%Rlay(k)*htot
881	0	Dfn = (GV%Rlay(k) - GV%Rlay(k-1))*(h_at_vel(i,j,k)+htot)
882
883	0	if (oldfn >= ustarsq) then
884	0	cycle
885	0	elseif ((oldfn + Dfn) <= ustarsq) then
886	0	Dh = h_at_vel(i,j,k)
887	0	Ddz = dz_at_vel(i,j,k)
888		else
889	0	frac_used = sqrt((ustarsq-oldfn) / (Dfn))
890	0	Dh = h_at_vel(i,j,k) * frac_used
891	0	Ddz = dz_at_vel(i,j,k) * frac_used
892		endif
893
894	0	htot = htot + Dh
895	0	dztot = dztot + Ddz
896	0	Rhtot = Rhtot + GV%Rlay(k)*Dh
897		enddo
898	0	if (nkml>0) then
899	0	do k=nkmb,2,-1
900	0	oldfn = Rhtot - Rml_vel(i,j,k)*htot
901	0	Dfn = (Rml_vel(i,j,k) - Rml_vel(i,j,k-1)) * (h_at_vel(i,j,k)+htot)
902
903	0	if (oldfn >= ustarsq) then
904	0	cycle
905	0	elseif ((oldfn + Dfn) <= ustarsq) then
906	0	Dh = h_at_vel(i,j,k)
907	0	Ddz = dz_at_vel(i,j,k)
908		else
909	0	frac_used = sqrt((ustarsq-oldfn) / (Dfn))
910	0	Dh = h_at_vel(i,j,k) * frac_used
911	0	Ddz = dz_at_vel(i,j,k) * frac_used
912		endif
913
914	0	htot = htot + Dh
915	0	dztot = dztot + Ddz
916	0	Rhtot = Rhtot + Rml_vel(i,j,k)*Dh
917		enddo
918	0	if (Rhtot - Rml_vel(i,j,1)*htot < ustarsq) then
919	0	htot = htot + h_at_vel(i,j,1)
920	0	dztot = dztot + dz_at_vel(i,j,1)
921		endif
922		else
923	0	if (Rhtot - GV%Rlay(1)*htot < ustarsq) then
924	0	htot = htot + h_at_vel(i,j,1)
925	0	dztot = dztot + dz_at_vel(i,j,1)
926		endif
927		endif
928		endif ! use_BBL_EOS
929
930		! Value of 2*f at u- or v-points.
931	118104	if (m==1) then ; C2f = G%CoriolisBu(I,J-1) + G%CoriolisBu(I,J)
932	58728	else ; C2f = G%CoriolisBu(I-1,J) + G%CoriolisBu(I,J) ; endif
933
934		! The thickness of a rotation limited BBL ignoring stratification is
935		! h_f ~ Cn u* / f (limit of KW99 eq. 2.20 for N->0).
936		! The buoyancy limit of BBL thickness (h_N) is already in the variable htot from above.
937		! Substituting x = h_N/h into KW99 eq. 2.20 yields the quadratic
938		! x^2 - x = (h_N / h_f)^2
939		! for which the positive root is
940		! xp = 1/2 + sqrt( 1/4 + (h_N/h_f)^2 )
941		! and thus h_bbl = h_N / xp . Since h_f = Cn u*/f and Cn=0.5
942		! xp = 1/2 + sqrt( 1/4 + (2 f h_N/u*)^2 )
943		! To avoid dividing by zero if u*=0 then
944		! xp u* = 1/2 u* + sqrt( 1/4 u*^2 + (2 f h_N)^2 )
945	118104	if (CS%cdrag * u2_bg(i,j) <= 0.0) then
946		! This avoids NaNs and overflows, and could be used in all cases,
947		! but is not bitwise identical to the current code.
948	0	ustH = ustar(i,j) ; root = sqrt(0.25ustH2 + (htotC2f)**2)
949	0	if (dztotustH <= (CS%BBL_thick_min+dz_neglect) (0.5*ustH + root)) then
950	0	bbl_thick = CS%BBL_thick_min
951		else
952		! The following expression reads
953		! h_bbl = h_N u* / ( 1/2 u* + sqrt( 1/4 u*^2 + ( 2 f h_N )^2 ) )
954		! which is h_bbl = h_N u/(xp u) as described above.
955	0	bbl_thick = (dztot * ustH) / (0.5*ustH + root)
956		endif
957		else
958		! The following expression reads
959		! h_bbl = h_N / ( 1/2 + sqrt( 1/4 + ( 2 f h_N / u* )^2 ) )
960		! which is h_bbl = h_N/xp as described above.
961	118104	bbl_thick = dztot / (0.5 + sqrt(0.25 + htothtotC2fC2f / (ustar(i,j)ustar(i,j)) ) )
962
963	118104	if (bbl_thick < CS%BBL_thick_min) bbl_thick = CS%BBL_thick_min
964		endif
965
966		! Store the normalized bottom boundary layer volume.
967	118104	if (CS%Channel_drag) Vol_bbl_chan = bbl_thick
968
969		! If there is Richardson number dependent mixing, that determines
970		! the vertical extent of the bottom boundary layer, and there is no
971		! need to set that scale here. In fact, viscously reducing the
972		! shears over an excessively large region reduces the efficacy of
973		! the Richardson number dependent mixing.
974		! In other words, if using RiNo_mix then CS%dz_bbl acts as an upper bound on
975		! bbl_thick.
976	118104	if ((bbl_thick > 0.5CS%dz_bbl) .and. (CS%RiNo_mix)) bbl_thick = 0.5CS%dz_bbl
977
978		! If drag is a body force, bbl_thick is HBBL
979	118104	if (CS%body_force_drag) bbl_thick = dz_bbl_drag(i,j)
980
981	118104	if (CS%Channel_drag) then
982
983	118104	vol_below(nz+1) = 0.0
984	8975904	do K=nz,1,-1
985	8975904	vol_below(K) = vol_below(K+1) + dz_vel(i,j,k)
986		enddo
987
988		! Find the bathymetry at adjacent points relative to the shelf break. For now this
989		! shelf break depth is set with a global constant, but it could vary in space.
990	118104	if (m==1) then
991	59376	D_vel = D_u(I,j) - CS%channel_break_depth
992	59376	D_vel_p = G%mask2dCu(I,j+1) * (D_u(I,j+1) - CS%channel_break_depth)
993	59376	D_vel_m = G%mask2dCu(I,j-1) * (D_u(I,j-1) - CS%channel_break_depth)
994		else
995	58728	D_vel = D_v(i,J) - CS%channel_break_depth
996	58728	D_vel_p = G%mask2dCv(i+1,J) * (D_v(i+1,J) - CS%channel_break_depth)
997	58728	D_vel_m = G%mask2dCv(i-1,J) * (D_v(i-1,J) - CS%channel_break_depth)
998		endif
999		! This profile uses a harmonic mean bottom depth below some reference value to
1000		! roughly mimic the topographic shape at and beneath a continental shelf break.
1001		! Above this a simple arithmetic mean is used.
1002	118104	if ((D_vel > 0.0) .and. (D_vel_p > 0.0)) then
1003	115644	Dp = 2.0 * D_vel * D_vel_p / (D_vel + D_vel_p)
1004		else ! This is above the shelf-break, noting that D is positive downward.
1005	2460	Dp = 0.5 * (min(D_vel, 0.0) + min(D_vel_p, 0.0))
1006		endif
1007	118104	if ((D_vel > 0.0) .and. (D_vel_m > 0.0)) then
1008	115644	Dm = 2.0 * D_vel * D_vel_m / (D_vel + D_vel_m)
1009		else ! This is above the shelf-break, noting that D is positive downward.
1010	2460	Dm = 0.5 * (min(D_vel, 0.0) + min(D_vel_m, 0.0))
1011		endif
1012	118104	if (Dm > Dp) then ; tmp = Dp ; Dp = Dm ; Dm = tmp ; endif
1013	118104	crv = 3.0(Dp + Dm - 2.0D_vel)
1014	118104	slope = Dp - Dm
1015
1016		! If the curvature is small enough, there is no reason not to assume
1017		! a uniformly sloping or flat bottom.
1018	118104	if (abs(crv) < 1e-2*(slope + CS%BBL_thick_min)) crv = 0.0
1019
1020		! Determine the normalized open length (L) at each interface.
1021	118104	if (crv == 0.0) then
1022	768	call find_L_open_uniform_slope(vol_below, Dp, Dm, L, GV)
1023	117336	elseif (crv > 0.0) then
1024	11664	if (CS%concave_trigonometric_L) then
1025	0	call find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, L, GV)
1026		else
1027	11664	call find_L_open_concave_iterative(vol_below, D_vel, Dp, Dm, L, GV)
1028	11664	if (CS%debug) then
1029		! The tests in this block reveal that the iterative and trigonometric solutions are
1030		! mathematically equivalent, but in some cases the iterative solution is consistent
1031		! at roundoff, but that the trigonmetric solutions have errors that can be several
1032		! orders of magnitude larger in some cases.
1033	0	call find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, L_trig, GV)
1034	0	call test_L_open_concave(vol_below, D_vel, Dp, Dm, L_trig, vol_err_trig, GV)
1035	0	call test_L_open_concave(vol_below, D_vel, Dp, Dm, L, vol_err_iter, GV)
1036	0	max_dL_trig_itt = 0.0 ; max_norm_err_trig = 0.0 ; max_norm_err_iter = 0.0
1037	0	norm_err_trig(:) = 0.0 ; norm_err_iter(:) = 0.0
1038	0	do K=1,nz+1
1039	0	dL_trig_itt(K) = L_trig(K) - L(K)
1040	0	if (abs(dL_trig_itt(K)) > abs(max_dL_trig_itt)) max_dL_trig_itt = dL_trig_itt(K)
1041	0	norm_err_trig(K) = vol_err_trig(K) / (vol_below(K) + dz_neglect)
1042	0	norm_err_iter(K) = vol_err_iter(K) / (vol_below(K) + dz_neglect)
1043	0	if (abs(norm_err_trig(K)) > abs(max_norm_err_trig)) max_norm_err_trig = norm_err_trig(K)
1044	0	if (abs(norm_err_iter(K)) > abs(max_norm_err_iter)) max_norm_err_iter = norm_err_iter(K)
1045		enddo
1046	0	if (abs(max_dL_trig_itt) > 1.0e-13) &
1047	0	K = nz+1 ! This is here only to use as a break point for a debugger.
1048	0	if (abs(max_norm_err_trig) > 1.0e-13) &
1049	0	K = nz+1 ! This is here only to use as a break point for a debugger.
1050	0	if (abs(max_norm_err_iter) > 1.0e-13) &
1051	0	K = nz+1 ! This is here only to use as a break point for a debugger.
1052		endif
1053		endif
1054		else ! crv < 0.0
1055	105672	call find_L_open_convex(vol_below, D_vel, Dp, Dm, L, GV, US, CS)
1056		endif ! end of crv<0 cases.
1057
1058		! Determine the Rayleigh drag contributions.
1059
1060		! The drag within the bottommost Vol_bbl_chan is applied as a part of an enhanced bottom
1061		! viscosity, while above this the drag is applied directly to the layers in question as a
1062		! Rayleigh drag term.
1063
1064		! Restrict the volume over which the channel drag is applied from the previously determined value.
1065	118104	if (CS%Chan_drag_max_vol >= 0.0) Vol_bbl_chan = min(Vol_bbl_chan, CS%Chan_drag_max_vol)
1066
1067	118104	BBL_visc_frac = 0.0
1068	8975904	do K=nz,1,-1
1069		!modify L(K) for porous barrier parameterization
1070	8857800	if (m==1) then ; L(K) = L(K)*pbv%por_layer_widthU(I,j,K)
1071	4404600	else ; L(K) = L(K)*pbv%por_layer_widthV(i,J,K) ; endif
1072
1073		! Determine the drag contributing to the bottom boundary layer
1074		! and the Rayleigh drag that acts on each layer.
1075	8857800	if (L(K) > L(K+1)) then
1076	2944908	if (vol_below(K+1) < Vol_bbl_chan) then
1077	2662188	BBL_frac = (1.0-vol_below(K+1)/Vol_bbl_chan)**2
1078	2662188	BBL_visc_frac = BBL_visc_frac + BBL_frac*(L(K) - L(K+1))
1079		else
1080	282720	BBL_frac = 0.0
1081		endif
1082
1083	2944908	if (allocated(tv%SpV_avg)) then
1084	0	cdrag_conv = cdrag_RL_to_H / SpV_vel(i,j,k)
1085		else
1086	2944908	cdrag_conv = cdrag_L_to_H
1087		endif
1088
1089	2944908	h_vel_pos = h_vel(i,j,k) + h_neglect
1090	2944908	if (m==1) then ; Cell_width = G%dy_Cu(I,j)*pbv%por_face_areaU(I,j,k)
1091	1456992	else ; Cell_width = G%dx_Cv(i,J)*pbv%por_face_areaV(i,J,k) ; endif
1092	2944908	gam = 1.0 - L(K+1)/L(K)
1093		Rayleigh = cdrag_conv * (L(K)-L(K+1)) * (1.0-BBL_frac) * &
1094		(12.0CS%c_Smagh_vel_pos) / (12.0CS%c_Smagh_vel_pos + &
1095	2944908	cdrag_conv * gam(1.0-gam)(1.0-1.5gam) L(K)*2 Cell_width)
1096		else ! This layer feels no drag.
1097	5912892	Rayleigh = 0.0
1098		endif
1099
1100	8975904	if (m==1) then
1101	4453200	if (Rayleigh > 0.0) then
1102	1428468	v_at_u = set_v_at_u(v, h, G, GV, i, j, k, mask_v, OBC)
1103	1428468	visc%Ray_u(I,j,k) = Rayleigh * sqrt(u(I,j,k)u(I,j,k) + v_at_uv_at_u + u2_bg(I,j))
1104	3024732	else ; visc%Ray_u(I,j,k) = 0.0 ; endif
1105		else
1106	4404600	if (Rayleigh > 0.0) then
1107	1398216	u_at_v = set_u_at_v(u, h, G, GV, i, j, k, mask_u, OBC)
1108	1398216	visc%Ray_v(i,J,k) = Rayleigh * sqrt(v(i,J,k)v(i,J,k) + u_at_vu_at_v + u2_bg(i,j))
1109	3006384	else ; visc%Ray_v(i,J,k) = 0.0 ; endif
1110		endif
1111
1112		enddo ! k loop to determine visc%Ray_[uv].
1113
1114		! Set the near-bottom viscosity to a value which will give
1115		! the correct stress when the shear occurs over bbl_thick.
1116		! See next block for explanation.
1117	118104	if (CS%correct_BBL_bounds .and. &
1118		cdrag_sqrtustar(i,j)bbl_thick*BBL_visc_frac <= CS%Kv_BBL_min) then
1119		! If the bottom stress implies less viscosity than Kv_BBL_min then
1120		! set kv_bbl to the bound and recompute bbl_thick to be consistent
1121		! but with a ridiculously large upper bound on thickness (for Cd u*=0)
1122	0	kv_bbl = CS%Kv_BBL_min
1123	0	if ((cdrag_sqrtustar(i,j))BBL_visc_frac*BBL_thick_max > kv_bbl) then
1124	0	bbl_thick = kv_bbl / ( (cdrag_sqrtustar(i,j)) BBL_visc_frac )
1125		else
1126	0	bbl_thick = BBL_thick_max
1127		endif
1128		else
1129	118104	kv_bbl = (cdrag_sqrtustar(i,j)) bbl_thick*BBL_visc_frac
1130		endif
1131
1132		else ! Not Channel_drag.
1133		! Set the near-bottom viscosity to a value which will give
1134		! the correct stress when the shear occurs over bbl_thick.
1135		! - The bottom stress is tau_b = Cdrag * u_bbl^2
1136		! - u_bbl was calculated by averaging flow over CS%Hbbl
1137		! (and includes unresolved tidal components)
1138		! - u_bbl is embedded in u* since u*^2 = Cdrag u_bbl^2
1139		! - The average shear in the BBL is du/dz = 2 * u_bbl / h_bbl
1140		! (which assumes a linear profile, hence the "2")
1141		! - bbl_thick was bounded to <= 0.5 * CS%dz_bbl
1142		! - The viscous stress kv_bbl du/dz should balance tau_b
1143		! Cdrag u_bbl^2 = kv_bbl du/dz
1144		! = 2 kv_bbl u_bbl
1145		! so
1146		! kv_bbl = 0.5 h_bbl Cdrag u_bbl
1147		! = 0.5 h_bbl sqrt(Cdrag) u*
1148	0	if (CS%correct_BBL_bounds .and. &
1149		cdrag_sqrtustar(i,j)bbl_thick <= CS%Kv_BBL_min) then
1150		! If the bottom stress implies less viscosity than Kv_BBL_min then
1151		! set kv_bbl to the bound and recompute bbl_thick to be consistent
1152		! but with a ridiculously large upper bound on thickness (for Cd u*=0)
1153	0	kv_bbl = CS%Kv_BBL_min
1154	0	if ((cdrag_sqrtustar(i,j))BBL_thick_max > kv_bbl) then
1155	0	bbl_thick = kv_bbl / ( cdrag_sqrt*ustar(i,j) )
1156		else
1157	0	bbl_thick = BBL_thick_max
1158		endif
1159		else
1160	0	kv_bbl = (cdrag_sqrtustar(i,j)) bbl_thick
1161		endif
1162		endif
1163
1164	118104	if (CS%body_force_drag) then ; if (h_bbl_drag(i,j) > 0.0) then
1165		! Increment the Rayleigh drag as a way introduce the bottom drag as a body force.
1166	0	h_sum = 0.0
1167	0	I_hwtot = 1.0 / h_bbl_drag(i,j)
1168	0	do k=nz,1,-1
1169	0	h_bbl_fr = min(h_bbl_drag(i,j) - h_sum, h_at_vel(i,j,k)) * I_hwtot
1170	0	if (allocated(tv%SpV_avg)) then
1171	0	cdrag_conv = cdrag_RL_to_H / SpV_vel(i,j,k)
1172		else
1173	0	cdrag_conv = cdrag_L_to_H
1174		endif
1175	0	if (m==1) then
1176	0	visc%Ray_u(I,j,k) = visc%Ray_u(I,j,k) + (cdrag_conv * umag_avg(I,j)) * h_bbl_fr
1177		else
1178	0	visc%Ray_v(i,J,k) = visc%Ray_v(i,J,k) + (cdrag_conv * umag_avg(i,j)) * h_bbl_fr
1179		endif
1180	0	h_sum = h_sum + h_at_vel(i,j,k)
1181	0	if (h_sum >= h_bbl_drag(i,j)) exit ! The top of this layer is above the drag zone.
1182		enddo
1183		! Do not enhance the near-bottom viscosity in this case.
1184	0	Kv_bbl = CS%Kv_BBL_min
1185		endif ; endif
1186
1187	118104	kv_bbl = max(CS%Kv_BBL_min, kv_bbl)
1188	118104	if (m==1) then
1189	59376	visc%bbl_thick_u(I,j) = bbl_thick
1190	59376	if (allocated(visc%Kv_bbl_u)) visc%Kv_bbl_u(I,j) = kv_bbl
1191		else
1192	58728	visc%bbl_thick_v(i,J) = bbl_thick
1193	58728	if (allocated(visc%Kv_bbl_v)) visc%Kv_bbl_v(i,J) = kv_bbl
1194		endif
1195		endif ; enddo ; enddo ! end of i & j loops
1196		enddo ! end of m loop
1197
1198		!$omp target exit data map(release: dz, tv, tv%T, tv%S, S_vel, T_vel, SpV_vel, h_vel, h_at_vel, &
1199		!$omp dz_vel, dz_at_vel, Rml, Rml_vel, p_ref, ustar, umag_avg, u2_bg, mask_u, mask_v, &
1200		!$omp h_bbl_drag, dz_bbl_drag, do_i, dR_dS, dR_dT, D_u, D_v, press, S_EOS, T_EOS, tv%p_surf, CS)
1201
1202		! Offer diagnostics for averaging
1203	12	if (CS%id_bbl_thick_u > 0) &
1204	0	call post_data(CS%id_bbl_thick_u, visc%bbl_thick_u, CS%diag)
1205	12	if (CS%id_kv_bbl_u > 0) &
1206	0	call post_data(CS%id_kv_bbl_u, visc%kv_bbl_u, CS%diag)
1207	12	if (CS%id_bbl_u > 0) &
1208	0	call post_data(CS%id_bbl_u, CS%bbl_u, CS%diag)
1209	12	if (CS%id_bbl_thick_v > 0) &
1210	0	call post_data(CS%id_bbl_thick_v, visc%bbl_thick_v, CS%diag)
1211	12	if (CS%id_kv_bbl_v > 0) &
1212	0	call post_data(CS%id_kv_bbl_v, visc%kv_bbl_v, CS%diag)
1213	12	if (CS%id_bbl_v > 0) &
1214	0	call post_data(CS%id_bbl_v, CS%bbl_v, CS%diag)
1215	12	if (CS%id_Ray_u > 0) &
1216	0	call post_data(CS%id_Ray_u, visc%Ray_u, CS%diag)
1217	12	if (CS%id_Ray_v > 0) &
1218	0	call post_data(CS%id_Ray_v, visc%Ray_v, CS%diag)
1219
1220	12	if (CS%debug) then
1221	0	if (allocated(visc%Ray_u) .and. allocated(visc%Ray_v)) then
1222		!$omp target update from(visc%Ray_u, visc%Ray_v)
1223		call uvchksum("Ray [uv]", visc%Ray_u, visc%Ray_v, G%HI, haloshift=0, &
1224	0	unscale=GV%H_to_m*US%s_to_T, scalar_pair=.true.)
1225		endif
1226	0	if (allocated(visc%kv_bbl_u) .and. allocated(visc%kv_bbl_v)) then
1227		!$omp target update from(visc%Kv_bbl_u, visc%Kv_bbl_v)
1228		call uvchksum("kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, G%HI, &
1229	0	haloshift=0, unscale=GV%HZ_T_to_m2_s, scalar_pair=.true.)
1230		endif
1231	0	if (allocated(visc%bbl_thick_u) .and. allocated(visc%bbl_thick_v)) then
1232		!$omp target update from(visc%bbl_thick_u, visc%bbl_thick_v)
1233		call uvchksum("bbl_thick_[uv]", visc%bbl_thick_u, visc%bbl_thick_v, &
1234	0	G%HI, haloshift=0, unscale=US%Z_to_m, scalar_pair=.true.)
1235		endif
1236		endif
1237		end subroutine set_viscous_BBL
1238
1239		!> Determine the normalized open length of each interface, given the edge depths and normalized
1240		!! volumes below each interface.
1241	768	pure subroutine find_L_open_uniform_slope(vol_below, Dp, Dm, L, GV)
1242		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1243		real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by
1244		!! the full horizontal area of a velocity cell [Z ~> m]
1245		real, intent(in) :: Dp !< The larger of the two depths at the edge
1246		!! of a velocity cell [Z ~> m]
1247		real, intent(in) :: Dm !< The smaller of the two depths at the edge
1248		!! of a velocity cell [Z ~> m]
1249		real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at
1250		!! the depth of each interface [nondim]
1251		!$omp declare target
1252
1253		! Local variables
1254		real :: slope ! The absolute value of the bottom depth slope across a cell times the cell width [Z ~> m].
1255		real :: I_slope ! The inverse of the normalized slope [Z-1 ~> m-1]
1256		real :: Vol_open ! The cell volume above which it is open [Z ~> m].
1257		integer :: K, nz
1258		!$omp declare target
1259
1260	768	nz = GV%ke
1261
1262	768	slope = abs(Dp - Dm)
1263	768	if (slope == 0.0) then
1264	0	L(1:nz) = 1.0 ; L(nz+1) = 0.0
1265		else
1266	768	Vol_open = 0.5*slope
1267	768	I_slope = 1.0 / slope
1268
1269	768	L(nz+1) = 0.0
1270	58368	do K=nz,1,-1
1271	58368	if (vol_below(K) >= Vol_open) then ; L(K) = 1.0
1272		else
1273		! With a uniformly sloping bottom, the calculation of L(K) is the solution of a simple quadratic equation.
1274	15312	L(K) = sqrt(2.0vol_below(K)I_slope)
1275		endif
1276		enddo
1277		endif
1278
1279	768	end subroutine find_L_open_uniform_slope
1280
1281		!> Determine the normalized open length of each interface for concave bathymetry (from the ocean perspective)
1282		!! using trigonometric expressions. In this case there can be two separate open regions.
1283	0	pure subroutine find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, L, GV)
1284		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1285		real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by
1286		!! the full horizontal area of a velocity cell [Z ~> m]
1287		real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m]
1288		real, intent(in) :: Dp !< The larger of the two depths at the edge
1289		!! of a velocity cell [Z ~> m]
1290		real, intent(in) :: Dm !< The smaller of the two depths at the edge
1291		!! of a velocity cell [Z ~> m]
1292		real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at
1293		!! the depth of each interface [nondim]
1294		!$omp declare target
1295
1296		! Local variables
1297		real :: crv ! crv is the curvature of the bottom depth across a
1298		! cell, times the cell width squared [Z ~> m].
1299		real :: crv_3 ! crv/3 [Z ~> m].
1300		real :: slope ! The absolute value of the bottom depth slope across
1301		! a cell times the cell width [Z ~> m].
1302		! The following "volumes" have units of vertical heights because they are normalized
1303		! by the full horizontal area of a velocity cell.
1304		real :: Vol_open ! The cell volume above which the face is fully is open [Z ~> m].
1305		real :: Vol_2_reg ! The cell volume above which there are two separate
1306		! open areas that must be integrated [Z ~> m].
1307		real :: C24_crv ! 24/crv [Z-1 ~> m-1].
1308		real :: apb_4a, ax2_3apb ! Various nondimensional ratios of crv and slope [nondim].
1309		real :: a2x48_apb3, Iapb ! Combinations of crv (a) and slope (b) [Z-1 ~> m-1]
1310		real :: L0 ! A linear estimate of L appropriate for tiny volumes [nondim].
1311		real :: slope_crv ! The slope divided by the curvature [nondim]
1312		real :: tmp_val_m1_to_p1 ! A temporary variable [nondim]
1313		real, parameter :: C1_3 = 1.0/3.0, C1_12 = 1.0/12.0 ! Rational constants [nondim]
1314		real, parameter :: C2pi_3 = 8.0*atan(1.0)/3.0 ! An irrational constant, 2/3 pi. [nondim]
1315		integer :: K, nz
1316		!$omp declare target
1317
1318	0	nz = GV%ke
1319
1320		! Each cell extends from x=-1/2 to 1/2, and has a topography
1321		! given by D(x) = crvx^2 + slopex + D_vel - crv/12.
1322		!crv_3 = (Dp + Dm - 2.0D_vel) ; crv = 3.0crv_3
1323	0	crv_3 = (Dp + Dm - (2.0D_vel)) ; crv = 3.0crv_3
1324	0	slope = Dp - Dm
1325
1326		! Calculate the volume above which the entire cell is open and the volume at which the
1327		! equation that is solved for L changes because there are two separate open regions.
1328	0	if (slope >= crv) then
1329	0	Vol_open = D_vel - Dm ; Vol_2_reg = Vol_open
1330		else
1331	0	slope_crv = slope / crv
1332	0	Vol_open = 0.25slopeslope_crv + C1_12*crv
1333	0	Vol_2_reg = 0.5slope_crv2 (crv - C1_3*slope)
1334		endif
1335		! Define some combinations of crv & slope for later use.
1336	0	C24_crv = 24.0/crv ; Iapb = 1.0/(crv+slope)
1337	0	apb_4a = (slope+crv)/(4.0crv) ; a2x48_apb3 = (48.0(crvcrv))(Iapb**3)
1338	0	ax2_3apb = 2.0C1_3crv*Iapb
1339
1340	0	L(nz+1) = 0.0
1341		! Determine the normalized open length (L) at each interface.
1342	0	do K=nz,1,-1
1343	0	if (vol_below(K) >= Vol_open) then ! The whole cell is open.
1344	0	L(K) = 1.0
1345	0	elseif (vol_below(K) < Vol_2_reg) then
1346		! In this case, there is a contiguous open region and
1347		! vol_below(K) = 0.5L^2(slope + crv/3*(3-4L)).
1348	0	if (a2x48_apb3*vol_below(K) < 1e-8) then ! Could be 1e-7?
1349		! There is a very good approximation here for massless layers.
1350		!L0 = sqrt(2.0vol_below(K)Iapb) ; L(K) = L0(1.0 + ax2_3apbL0)
1351	0	L0 = sqrt(2.0vol_below(K)Iapb) ; L(K) = L0(1.0 + (ax2_3apbL0))
1352		else
1353		!L(K) = apb_4a * (1.0 - &
1354		! 2.0 * cos(C1_3acos(a2x48_apb3vol_below(K) - 1.0) - C2pi_3))
1355		L(K) = apb_4a * (1.0 - &
1356	0	2.0 * cos(C1_3acos((a2x48_apb3vol_below(K)) - 1.0) - C2pi_3))
1357		endif
1358		! To check the answers.
1359		! Vol_err = 0.5(L(K)L(K))(slope + crv_3(3.0-4.0*L(K))) - vol_below(K)
1360		else ! There are two separate open regions.
1361		! vol_below(K) = slope^2/4crv + crv/12 - (crv/12)(1-L)^2(1+2L)
1362		! At the deepest volume, L = slope/crv, at the top L = 1.
1363		! L(K) = 0.5 - cos(C1_3acos(1.0 - C24_crv(Vol_open - vol_below(K))) - C2pi_3)
1364	0	tmp_val_m1_to_p1 = 1.0 - C24_crv*(Vol_open - vol_below(K))
1365	0	tmp_val_m1_to_p1 = max(-1., min(1., tmp_val_m1_to_p1))
1366	0	L(K) = 0.5 - cos(C1_3*acos(tmp_val_m1_to_p1) - C2pi_3)
1367		! To check the answers.
1368		! Vol_err = Vol_open - 0.25crv_3(1.0+2.0L(K)) (1.0-L(K))**2 - vol_below(K)
1369		endif
1370		enddo ! k loop to determine L(K) in the concave case
1371
1372	0	end subroutine find_L_open_concave_trigonometric
1373
1374
1375
1376		!> Determine the normalized open length of each interface for concave bathymetry (from the ocean perspective) using
1377		!! iterative methods to solve the relevant cubic equations. In this case there can be two separate open regions.
1378	11664	pure subroutine find_L_open_concave_iterative(vol_below, D_vel, Dp, Dm, L, GV)
1379		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1380		real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by
1381		!! the full horizontal area of a velocity cell [Z ~> m]
1382		real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m]
1383		real, intent(in) :: Dp !< The larger of the two depths at the edge
1384		!! of a velocity cell [Z ~> m]
1385		real, intent(in) :: Dm !< The smaller of the two depths at the edge
1386		!! of a velocity cell [Z ~> m]
1387		real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at
1388		!! the depth of each interface [nondim]
1389		!$omp declare target
1390
1391		! Local variables
1392		real :: crv ! crv is the curvature of the bottom depth across a
1393		! cell, times the cell width squared [Z ~> m].
1394		real :: crv_3 ! crv/3 [Z ~> m].
1395		real :: slope ! The absolute value of the bottom depth slope across
1396		! a cell times the cell width [Z ~> m].
1397
1398		! The following "volumes" have units of vertical heights because they are normalized
1399		! by the full horizontal area of a velocity cell.
1400		real :: Vol_open ! The cell volume above which the face is fully is open [Z ~> m].
1401		real :: Vol_2_reg ! The cell volume above which there are two separate
1402		! open areas that must be integrated [Z ~> m].
1403		real :: L_2_reg ! The value of L when vol_below is Vol_2_reg [nondim]
1404		real :: vol_inflect_1 ! The volume at which there is an inflection point in the expression
1405		! relating L to vol_err when there is a single open region [Z ~> m]
1406		real :: vol_inflect_2 ! The volume at which there is an inflection point in the expression
1407		! relating L to vol_err when there are two open regions [Z ~> m]
1408
1409		real :: L_inflect_1 ! The value of L that sits at an inflection point in the expression
1410		! relating L to vol_err when there is a single open region [nondim]
1411		real :: L_inflect_2 ! The value of L that sits at an inflection point in the expression
1412		! relating L to vol_err when there is are two open regions [nondim]
1413		real :: L_max, L_min ! Maximum and minimum bounds on the solution for L for an interface [nondim]
1414		real :: vol_err ! The difference between the volume below an interface for a given value
1415		! of L and the target value [Z ~> m]
1416		real :: dVol_dL ! The partial derivative of the volume below with L [Z ~> m]
1417		real :: vol_err_max ! The value of vol_err when L is L_max [Z ~> m]
1418
1419		! The following combinations of slope and crv are reused across layers, and hence are pre-calculated
1420		! for efficiency. All are non-negative.
1421		real :: Icrvpslope ! The inverse of the sum of crv and slope [Z-1 ~> m-1]
1422		real :: slope_crv ! The slope divided by the curvature [nondim]
1423		! These are only used if the slope exceeds or matches the curvature.
1424		real :: smc ! The slope minus the curvature [Z ~> m]
1425		real :: C3c_m_s ! 3 times the curvature minus the slope [Z ~> m]
1426		real :: I_3c_m_s ! The inverse of 3 times the curvature minus the slope [Z-1 ~> m-1]
1427		! These are only used if the curvature exceeds the slope.
1428		real :: C4_crv ! The inverse of a quarter of the curvature [Z-1 ~> m-1]
1429		real :: sxcms_c ! The slope times the difference between the curvature and slope
1430		! divided by the curvature [Z ~> m]
1431		real :: slope2_4crv ! A quarter of the slope squared divided by the curvature [Z ~> m]
1432		real :: I_3s_m_c ! The inverse of 3 times the slope minus the curvature [Z-1 ~> m-1]
1433		real :: C3s_m_c ! 3 times the slope minus the curvature [Z ~> m]
1434
1435		real, parameter :: C1_3 = 1.0 / 3.0, C1_12 = 1.0 / 12.0 ! Rational constants [nondim]
1436		integer :: K, nz, itt
1437		integer, parameter :: max_itt = 10
1438		!$omp declare target
1439
1440	11664	nz = GV%ke
1441
1442		! Each cell extends from x=-1/2 to 1/2, and has a topography
1443		! given by D(x) = crvx^2 + slopex + D_vel - crv/12.
1444
1445	11664	crv_3 = (Dp + Dm - 2.0D_vel) ; crv = 3.0crv_3
1446	11664	slope = Dp - Dm
1447
1448		! Calculate the volume above which the entire cell is open and the volume at which the
1449		! equation that is solved for L changes because there are two separate open regions.
1450	11664	if (slope >= crv) then
1451	9024	Vol_open = D_vel - Dm ; Vol_2_reg = Vol_open
1452	9024	L_2_reg = 1.0
1453	9024	if (crv + slope >= 4.0*crv) then
1454	5544	L_inflect_1 = 1.0 ; Vol_inflect_1 = Vol_open
1455		else
1456	3480	slope_crv = slope / crv
1457	3480	L_inflect_1 = 0.25 + 0.25*slope_crv
1458	3480	vol_inflect_1 = 0.25C1_12 ((slope_crv + 1.0)*2 (slope + crv))
1459		endif
1460		! Precalculate some combinations of crv & slope for later use.
1461	9024	smc = slope - crv
1462	9024	C3c_m_s = 3.0*crv - slope
1463	9024	if (C3c_m_s > 2.0*smc) I_3c_m_s = 1.0 / C3c_m_s
1464		else
1465	2640	slope_crv = slope / crv
1466	2640	Vol_open = 0.25slopeslope_crv + C1_12*crv
1467	2640	Vol_2_reg = 0.5slope_crv2 (crv - C1_3*slope)
1468	2640	L_2_reg = slope_crv
1469
1470		! The inflection point is useful to know because below the inflection point
1471		! Newton's method converges monotonically from above and conversely above it.
1472		! These are the inflection point values of L and vol_below with a single open segment.
1473	2640	vol_inflect_1 = 0.25C1_12 ((slope_crv + 1.0)*2 (slope + crv))
1474	2640	L_inflect_1 = 0.25 + 0.25*slope_crv
1475		! These are the inflection point values of L and vol_below when there are two open segments.
1476		! Vol_inflect_2 = Vol_open - 0.125 * crv_3, which is equivalent to:
1477	2640	vol_inflect_2 = 0.25slopeslope_crv + 0.125*crv_3
1478	2640	L_inflect_2 = 0.5
1479		! Precalculate some combinations of crv & slope for later use.
1480	2640	C4_crv = 4.0 / crv
1481	2640	slope2_4crv = 0.25 * slope * slope_crv
1482	2640	sxcms_c = slope_crv*(crv - slope)
1483	2640	C3s_m_c = 3.0*slope - crv
1484	2640	if (C3s_m_c > 2.0*sxcms_c) I_3s_m_c = 1.0 / C3s_m_c
1485		endif
1486		! Define some combinations of crv & slope for later use.
1487	11664	Icrvpslope = 1.0 / (crv+slope)
1488
1489	11664	L(nz+1) = 0.0
1490		! Determine the normalized open length (L) at each interface.
1491	886464	do K=nz,1,-1
1492	886464	if (vol_below(K) >= Vol_open) then ! The whole cell is open.
1493	553728	L(K) = 1.0
1494	321072	elseif (vol_below(K) < Vol_2_reg) then
1495		! In this case, there is a single contiguous open region from x=1/2-L to 1/2.
1496		! Changing the horizontal variable in the expression from D(x) to D(L) gives:
1497		! x(L) = 1/2 - L
1498		! D(L) = crv(0.5 - L)^2 + slope(0.5 - L) + D_vel - crv/12
1499		! D(L) = crvL^2 - crvL + crv/4 + slope*(1/2 - L) + D_vel - crv/12
1500		! D(L) = crvL^2 - (slope+crv)L + slope/2 + D_vel + crv/6
1501		! D(0) = slope/2 + D_vel + crv/6 = (Dp - Dm)/2 + D_vel + (Dp + Dm - 2*D_vel)/2 = Dp
1502		! D(1) = crv - slope - crv + slope/2 + Dvel + crv/6 = D_vel - slope/2 + crv/6 = Dm
1503		!
1504		! vol_below = integral(y = 0 to L) D(y) dy - L * D(L)
1505		! = crv/3L^3 - (slope+crv)/2L^2 + (slope/2 + D_vel + crv/6)*L -
1506		! (crvL^2 - (slope+crv)L + slope/2 + D_vel + crv/6) * L
1507		! = -2/3 * crv * L^3 + 1/2 * (slope+crv) * L^2
1508		! vol_below(K) = 0.5L(K)2(slope + crv_3(3-4L(K)))
1509		! L(K) is between L(K+1) and slope_crv.
1510	319080	L_max = min(L_2_reg, 1.0)
1511	319080	if (vol_below(K) <= vol_inflect_1) L_max = min(L_max, L_inflect_1)
1512
1513	319080	L_min = L(K+1)
1514	319080	if (vol_below(K) >= vol_inflect_1) L_min = max(L_min, L_inflect_1)
1515
1516		! Ignoring the cubic term gives an under-estimate but is very accurate for near bottom
1517		! layers, so use this as a potential floor.
1518	319080	if (2.0vol_below(K)Icrvpslope > L_min*2) L_min = sqrt(2.0vol_below(K)*Icrvpslope)
1519
1520		! Start with L_min in most cases.
1521	319080	L(k) = L_min
1522
1523	319080	if (vol_below(K) <= vol_inflect_1) then
1524		! Starting with L_min below L_inflect_1, only the first overshooting iteration of Newton's
1525		! method needs bounding.
1526	318384	L(k) = L_min
1527	318384	vol_err = 0.5L(K)2 (slope + crv(1.0 - 4.0C1_3*L(K))) - vol_below(K)
1528		! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L_min is already the best solution.
1529	318384	if (vol_err < 0.0) then
1530	318384	dVol_dL = L(K) * (slope + crv(1.0 - 2.0L(k)))
1531	318384	if (L(K)dVol_dL > vol_err + L_maxdVol_dL) then
1532	24	L(K) = L_max
1533		else
1534	318360	L(K) = L(K) - (vol_err / dVol_dL)
1535		endif
1536
1537		! Subsequent iterations of Newton's method do not need bounds.
1538	937427	do itt=1,max_itt
1539	937427	vol_err = 0.5L(K)2 (slope + crv(1.0 - 4.0C1_3*L(K))) - vol_below(K)
1540	937427	dVol_dL = L(K) * (slope + crv(1.0 - 2.0L(k)))
1541	937427	if (abs(vol_err) < max(1.0e-15L(K), 1.0e-25)dVol_dL) exit
1542	619043	L(K) = L(K) - (vol_err / dVol_dL)
1543		enddo
1544		endif
1545		else ! (vol_below(K) > vol_inflect_1)
1546		! Iteration from below converges monotonically, but we need to deal with the case where we are
1547		! close to the peak of the topography and Newton's method mimics the convergence of bisection.
1548
1549		! Evaluate the error when L(K) = L_min as a possible first guess.
1550	696	L(k) = L_min
1551	696	vol_err = 0.5L(K)2 (slope + crv(1.0 - 4.0C1_3*L(K))) - vol_below(K)
1552		! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L_min is already the best solution.
1553	696	if (vol_err < 0.0) then
1554
1555		! These two upper estimates deal with the possibility that this point may be near
1556		! the upper extrema, where the error term might be approximately parabolic and
1557		! Newton's method would converge slowly like simple bisection.
1558	696	if (slope < crv) then
1559		! if ((L_2_reg - L_min)(3.0slope - crv) > 2.0slope_crv(crv-slope)) then
1560	168	if ((L_2_reg - L_min)C3s_m_c > 2.0sxcms_c) then
1561		! There is a decent upper estimate of L from the approximate quadratic equation found
1562		! by examining the error expressions at L ~= L_2_reg and ignoring the cubic term.
1563		L_max = (slope_crv(2.0slope) - sqrt(sxcms_c**2 + &
1564	48	2.0C3s_m_c(Vol_2_reg - vol_below(K))) ) * I_3s_m_c
1565		! The line above is equivalent to:
1566		! L_max = (slope_crv(2.0slope) - sqrt(slope_crv*2(crv-slope)**2 + &
1567		! 2.0(3.0slope - crv)*(Vol_2_reg - vol_below(K))) ) / &
1568		! (3.0*slope - crv)
1569		else
1570	120	L_max = slope_crv
1571		endif
1572		else ! (slope >= crv)
1573	528	if ((1.0 - L_min)C3c_m_s > 2.0smc) then
1574		! There is a decent upper estimate of L from the approximate quadratic equation found
1575		! by examining the error expressions at L ~= 1 and ignoring the cubic term.
1576	168	L_max = ( 2.0crv - sqrt(smc2 + 2.0C3c_m_s * (Vol_open - vol_below(K))) ) * I_3c_m_s
1577		! The line above is equivalent to:
1578		! L_max = ( 2.0crv - sqrt((slope - crv)2 + 2.0(3.0crv - slope) (Vol_open - vol_below(K))) ) / &
1579		! (3.0*crv - slope)
1580		else
1581	360	L_max = 1.0
1582		endif
1583		endif
1584	696	Vol_err_max = 0.5L_max2 (slope + crv(1.0 - 4.0C1_3*L_max)) - vol_below(K)
1585		! if (Vol_err_max < 0.0) call MOM_error(FATAL, &
1586		! "Vol_err_max should never be negative in find_L_open_concave_iterative.")
1587	696	if ((Vol_err_max < abs(Vol_err)) .and. (L_max < 1.0)) then
1588		! Start with 1 bounded Newton's method step from L_max
1589	264	dVol_dL = L_max * (slope + crv(1.0 - 2.0L_max))
1590	264	L(K) = max(L_min, L_max - (vol_err_max / dVol_dL) )
1591		! else ! Could use the fact that Vol_err is known to take an iteration?
1592		endif
1593
1594		! Subsequent iterations of Newton's method do not need bounds.
1595	2976	do itt=1,max_itt
1596	2976	vol_err = 0.5L(K)2 (slope + crv(1.0 - 4.0C1_3*L(K))) - vol_below(K)
1597	2976	dVol_dL = L(K) * (slope + crv(1.0 - 2.0L(k)))
1598	2976	if (abs(vol_err) < max(1.0e-15L(K), 1.0e-25)dVol_dL) exit
1599	2280	L(K) = L(K) - (vol_err / dVol_dL)
1600		enddo
1601		endif
1602
1603		endif
1604
1605		! To check the answers.
1606		! Vol_err = 0.5(L(K)L(K))(slope + crv_3(3.0-4.0*L(K))) - vol_below(K)
1607		else ! There are two separate open regions.
1608		! vol_below(K) = slope^2/(4crv) + crv/12 - (crv/12)(1-L)^2*(1+2L)
1609		! At the deepest volume, L = slope/crv, at the top L = 1.
1610
1611		! To check the answers.
1612		! Vol_err = Vol_open - 0.25crv_3(1.0+2.0L(K)) (1.0-L(K))**2 - vol_below(K)
1613		! or equivalently:
1614		! Vol_err = Vol_open - 0.25crv_3(3.0-2.0(1.0-L(K))) (1.0-L(K))**2 - vol_below(K)
1615		! ! Note that: Vol_open = 0.25slopeslope_crv + C1_12*crv
1616		! Vol_err = 0.25slopeslope_crv + 0.25crv_3( 1.0 - (1.0 + 2.0L(K)) (1.0-L(K))**2 ) - vol_below(K)
1617		! Vol_err = 0.25crv_3L(K)*2( 3.0 - 2.0L(K) ) + 0.25slope*slope_crv - vol_below(K)
1618
1619		! Derivation of the L_max limit below:
1620		! Vol_open - vol_below(K) = 0.25crv_3(3.0-2.0(1.0-L(K))) (1.0-L(K))**2
1621		! (3.0-2.0(1.0-L(K))) (1.0-L(K))*2 = (Vol_open - vol_below(K)) / (0.25crv_3)
1622		! When 1-L(K) << 1:
1623		! 3.0 * (1.0-L_max)*2 = (Vol_open - vol_below(K)) / (0.25crv_3)
1624		! (1.0-L_max)*2 = (Vol_open - vol_below(K)) / (0.25crv)
1625
1626		! Derivation of the L_min limit below:
1627		! Vol_err = 0.25crv_3L(K)*2( 3.0 - 2.0L(K) ) + 0.25slope*slope_crv - vol_below(K)
1628		! crvL(K)2( 1.0 - 2.0C1_3L(K) ) = 4.0vol_below(K) - slopeslope_crv
1629		! When L(K) << 1:
1630		! crvL_min2 = 4.0vol_below(K) - slope*slope_crv
1631		! L_min = sqrt((4.0vol_below(K) - slopeslope_crv)/crv)
1632		! Noting that L(K) >= slope_crv, when L(K)-slope_crv << 1:
1633		! (crv + 2.0C1_3slope)L_min2 = 4.0vol_below(K) - slope*slope_crv
1634		! L_min = sqrt((4.0vol_below(K) - slopeslope_crv)/(crv + 2.0C1_3slope))
1635
1636	1992	if (vol_below(K) <= Vol_inflect_2) then
1637		! Newton's Method would converge monotonically from above, but overshoot from below.
1638	1776	L_min = max(L(K+1), L_2_reg) ! L_2_reg = slope_crv
1639		! This under-estimate of L(K) is accurate for L ~= slope_crv:
1640	1776	if ((4.0vol_below(K) - slopeslope_crv) > (crv + 2.0C1_3slope)L_min*2) &
1641	960	L_min = max(L_min, sqrt((4.0vol_below(K) - slopeslope_crv) / (crv + 2.0C1_3slope)))
1642	1776	L_max = 0.5 ! = L_inflect_2
1643
1644		! Starting with L_min below L_inflect_2, only the first overshooting iteration of Newton's
1645		! method needs bounding.
1646	1776	L(k) = L_min
1647	1776	Vol_err = crv_3L(K)2( 0.75 - 0.5*L(K) ) + (slope2_4crv - vol_below(K))
1648
1649		! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L_min is already the best solution.
1650	1776	if (vol_err < 0.0) then
1651	1776	dVol_dL = 0.5crv (L(K) * (1.0 - L(K)))
1652	1776	if (L(K)dVol_dL >= vol_err + L_maxdVol_dL) then
1653	0	L(K) = L_max
1654		else
1655	1776	L(K) = L(K) - (vol_err / dVol_dL)
1656		endif
1657		! Subsequent iterations of Newton's method do not need bounds.
1658	5619	do itt=1,max_itt
1659	5619	Vol_err = crv_3 * (L(K)*2 (0.75 - 0.5*L(K))) + (slope2_4crv - vol_below(K))
1660	5619	dVol_dL = 0.5crv (L(K)*(1.0 - L(K)))
1661	5619	if (abs(vol_err) < max(1.0e-15L(K), 1.0e-25)dVol_dL) exit
1662	3843	L(K) = L(K) - (vol_err / dVol_dL)
1663		enddo
1664		endif
1665		else ! (vol_below(K) > Vol_inflect_2)
1666		! Newton's Method would converge monotonically from below, but overshoots from above, and
1667		! we may need to deal with the case where we are close to the peak of the topography.
1668	216	L_min = max(L(K+1), 0.5)
1669	216	L(k) = L_min
1670
1671	216	Vol_err = crv_3 * (L(K)*2 ( 0.75 - 0.5*L(K))) + (slope2_4crv - vol_below(K))
1672		! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L(k) is already the best solution.
1673	216	if (Vol_err < 0.0) then
1674		! This over-estimate of L(K) is accurate for L ~= 1:
1675	216	L_max = 1.0 - sqrt( (Vol_open - vol_below(K)) * C4_crv )
1676	216	Vol_err_max = crv_3 * (L_max*2 ( 0.75 - 0.5*L_max)) + (slope2_4crv - vol_below(K))
1677		! if (Vol_err_max < 0.0) call MOM_error(FATAL, &
1678		! "Vol_err_max should never be negative in find_L_open_concave_iterative.")
1679	216	if ((Vol_err_max < abs(Vol_err)) .and. (L_max < 1.0)) then
1680		! Start with 1 bounded Newton's method step from L_max
1681	216	dVol_dL = 0.5crv (L_max * (1.0 - L_max))
1682	216	L(K) = max(L_min, L_max - (vol_err_max / dVol_dL) )
1683		! else ! Could use the fact that Vol_err is known to take an iteration?
1684		endif
1685
1686		! Subsequent iterations of Newton's method do not need bounds.
1687	816	do itt=1,max_itt
1688	816	Vol_err = crv_3 * (L(K)*2 ( 0.75 - 0.5*L(K))) + (slope2_4crv - vol_below(K))
1689	816	dVol_dL = 0.5crv (L(K) * (1.0 - L(K)))
1690	816	if (abs(vol_err) < max(1.0e-15L(K), 1.0e-25)dVol_dL) exit
1691	600	L(K) = L(K) - (vol_err / dVol_dL)
1692		enddo
1693		endif
1694		endif
1695
1696		endif
1697		enddo ! k loop to determine L(K) in the concave case
1698
1699	11664	end subroutine find_L_open_concave_iterative
1700
1701
1702
1703		!> Test the validity the normalized open lengths of each interface for concave bathymetry (from the ocean perspective)
1704		!! by evaluating and returing the relevant cubic equations.
1705	0	pure subroutine test_L_open_concave(vol_below, D_vel, Dp, Dm, L, vol_err, GV)
1706		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1707		real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by
1708		!! the full horizontal area of a velocity cell [Z ~> m]
1709		real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m]
1710		real, intent(in) :: Dp !< The larger of the two depths at the edge
1711		!! of a velocity cell [Z ~> m]
1712		real, intent(in) :: Dm !< The smaller of the two depths at the edge
1713		!! of a velocity cell [Z ~> m]
1714		real, dimension(SZK_(GV)+1), intent(in) :: L !< The fraction of the full cell width that is open at
1715		!! the depth of each interface [nondim]
1716		real, dimension(SZK_(GV)+1), intent(out) :: vol_err !< The difference between vol_below and the
1717		!! value obtained from using L in the cubic equation [Z ~> m]
1718		!$omp declare target
1719
1720		! Local variables
1721		real :: crv ! crv is the curvature of the bottom depth across a
1722		! cell, times the cell width squared [Z ~> m].
1723		real :: crv_3 ! crv/3 [Z ~> m].
1724		real :: slope ! The absolute value of the bottom depth slope across
1725		! a cell times the cell width [Z ~> m].
1726
1727		! The following "volumes" have units of vertical heights because they are normalized
1728		! by the full horizontal area of a velocity cell.
1729		real :: Vol_open ! The cell volume above which the face is fully is open [Z ~> m].
1730		real :: Vol_2_reg ! The cell volume above which there are two separate
1731		! open areas that must be integrated [Z ~> m].
1732		real :: L_2_reg ! The value of L when vol_below is Vol_2_reg [nondim]
1733
1734		! The following combinations of slope and crv are reused across layers, and hence are pre-calculated
1735		! for efficiency. All are non-negative.
1736		real :: slope_crv ! The slope divided by the curvature [nondim]
1737		! These are only used if the curvature exceeds the slope.
1738		real :: slope2_4crv ! A quarter of the slope squared divided by the curvature [Z ~> m]
1739
1740		real, parameter :: C1_3 = 1.0 / 3.0, C1_12 = 1.0 / 12.0 ! Rational constants [nondim]
1741		integer :: K, nz
1742		!$omp declare target
1743
1744	0	nz = GV%ke
1745
1746		! Each cell extends from x=-1/2 to 1/2, and has a topography
1747		! given by D(x) = crvx^2 + slopex + D_vel - crv/12.
1748
1749	0	crv_3 = (Dp + Dm - 2.0D_vel) ; crv = 3.0crv_3
1750	0	slope = Dp - Dm
1751
1752		! Calculate the volume above which the entire cell is open and the volume at which the
1753		! equation that is solved for L changes because there are two separate open regions.
1754	0	if (slope >= crv) then
1755	0	Vol_open = D_vel - Dm ; Vol_2_reg = Vol_open
1756	0	L_2_reg = 1.0
1757	0	if (crv + slope >= 4.0*crv) then
1758	0	slope_crv = 1.0
1759		else
1760	0	slope_crv = slope / crv
1761		endif
1762		else
1763	0	slope_crv = slope / crv
1764	0	Vol_open = 0.25slopeslope_crv + C1_12*crv
1765	0	Vol_2_reg = 0.5slope_crv2 (crv - C1_3*slope)
1766	0	L_2_reg = slope_crv
1767		endif
1768	0	slope2_4crv = 0.25 * slope * slope_crv
1769
1770		! Determine the volume error based on the normalized open length (L) at each interface.
1771	0	Vol_err(nz+1) = 0.0
1772	0	do K=nz,1,-1
1773	0	if (L(K) >= 1.0) then
1774	0	Vol_err(K) = max(Vol_open - vol_below(K), 0.0)
1775	0	elseif (L(K) <= L_2_reg) then
1776	0	vol_err(K) = 0.5L(K)2 (slope + crv(1.0 - 4.0C1_3*L(K))) - vol_below(K)
1777		else ! There are two separate open regions.
1778	0	Vol_err(K) = crv_3 * (L(K)*2 ( 0.75 - 0.5*L(K))) + (slope2_4crv - vol_below(K))
1779		endif
1780		enddo ! k loop to determine L(K) in the concave case
1781
1782	0	end subroutine test_L_open_concave
1783
1784
1785		!> Determine the normalized open length of each interface for convex bathymetry (from the ocean
1786		!! perspective) using Newton's method iterations. In this case there is a single open region
1787		!! with the minimum depth at one edge of the cell.
1788	105672	pure subroutine find_L_open_convex(vol_below, D_vel, Dp, Dm, L, GV, US, CS)
1789		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1790		real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by
1791		!! the full horizontal area of a velocity cell [Z ~> m]
1792		real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m]
1793		real, intent(in) :: Dp !< The larger of the two depths at the edge
1794		!! of a velocity cell [Z ~> m]
1795		real, intent(in) :: Dm !< The smaller of the two depths at the edge
1796		!! of a velocity cell [Z ~> m]
1797		real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at
1798		!! the depth of each interface [nondim]
1799		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1800		type(set_visc_CS), intent(in) :: CS !< The control structure returned by a previous
1801		!! call to set_visc_init.
1802		!$omp declare target
1803
1804		! Local variables
1805		real :: crv ! crv is the curvature of the bottom depth across a
1806		! cell, times the cell width squared [Z ~> m].
1807		real :: crv_3 ! crv/3 [Z ~> m].
1808		real :: slope ! The absolute value of the bottom depth slope across
1809		! a cell times the cell width [Z ~> m].
1810		! All of the following "volumes" have units of vertical heights because they are normalized
1811		! by the full horizontal area of a velocity cell.
1812		real :: Vol_err ! The error in the volume with the latest estimate of
1813		! L, or the error for the interface below [Z ~> m].
1814		real :: Vol_quit ! The volume error below which to quit iterating [Z ~> m].
1815		real :: Vol_tol ! A volume error tolerance [Z ~> m].
1816		real :: Vol_open ! The cell volume above which the face is fully open [Z ~> m].
1817		real :: Vol_direct ! With less than Vol_direct [Z ~> m], there is a direct
1818		! solution of a cubic equation for L.
1819		real :: Vol_err_max ! The volume error for the upper bound on the correct value for L [Z ~> m]
1820		real :: Vol_err_min ! The volume error for the lower bound on the correct value for L [Z ~> m]
1821		real :: Vol_0 ! A deeper volume with known width L0 [Z ~> m].
1822		real :: dVol ! vol - Vol_0 [Z ~> m].
1823		real :: dV_dL2 ! The partial derivative of volume with L squared
1824		! evaluated at L=L0 [Z ~> m].
1825		real :: L_direct ! The value of L above volume Vol_direct [nondim].
1826		real :: L_max, L_min ! Upper and lower bounds on the correct value for L [nondim].
1827		real :: L0 ! The value of L above volume Vol_0 [nondim].
1828		real :: Iapb, Ibma_2 ! Combinations of crv (a) and slope (b) [Z-1 ~> m-1]
1829		real :: C24_crv ! 24/crv [Z-1 ~> m-1].
1830		real :: curv_tol ! Numerator of curvature cubed, used to estimate
1831		! accuracy of a single L(:) Newton iteration [Z5 ~> m5]
1832		real, parameter :: C1_3 = 1.0/3.0, C1_6 = 1.0/6.0 ! Rational constants [nondim]
1833		logical :: use_L0, do_one_L_iter ! Control flags for L(:) Newton iteration
1834		integer :: K, nz, itt
1835		integer, parameter:: maxitt = 20
1836		!$omp declare target
1837
1838	105672	nz = GV%ke
1839
1840		! Each cell extends from x=-1/2 to 1/2, and has a topography
1841		! given by D(x) = crvx^2 + slopex + D_vel - crv/12.
1842	105672	crv_3 = (Dp + Dm - 2.0D_vel) ; crv = 3.0crv_3
1843	105672	slope = Dp - Dm
1844
1845		! Calculate the volume above which the entire cell is open and the volume at which the
1846		! equation that is solved for L changes because there is a direct solution.
1847	105672	Vol_open = D_vel - Dm
1848	105672	if (slope >= -crv) then
1849	89760	Iapb = 1.0e30*US%Z_to_m ; if (slope+crv /= 0.0) Iapb = 1.0/(crv+slope)
1850	89760	Vol_direct = 0.0 ; L_direct = 0.0 ; C24_crv = 0.0
1851		else
1852	15912	C24_crv = 24.0/crv ; Iapb = 1.0/(crv+slope)
1853	15912	L_direct = 1.0 + slope/crv ! L_direct < 1 because crv < 0
1854	15912	Vol_direct = -C1_6crvL_direct**3
1855		endif
1856	105672	Ibma_2 = 2.0 / (slope - crv)
1857
1858	105672	if (CS%answer_date < 20190101) Vol_quit = (0.9*GV%Angstrom_Z + GV%dZ_subroundoff)
1859
1860	105672	L(nz+1) = 0.0 ; Vol_err = 0.0
1861		! Determine the normalized open length (L) at each interface.
1862	8031072	do K=nz,1,-1
1863	8031072	if (vol_below(K) >= Vol_open) then
1864	5431235	L(K) = 1.0
1865	2494165	elseif (vol_below(K) <= Vol_direct) then
1866	351360	if (CS%answer_date < 20260501) then
1867	0	L(K) = (-0.25C24_crvvol_below(K))**C1_3
1868		else
1869	351360	L(K) = cuberoot(-0.25C24_crvvol_below(K))
1870		endif
1871		else
1872		! x_R is at 1/2 but x_L is in the interior & L is found by iteratively solving
1873		! vol_below(K) = 0.5L^2(slope + crv/3*(3-4L))
1874
1875		! Vol_err = 0.5(L(K+1)L(K+1))(slope + crv_3(3.0-4.0*L(K+1))) - vol_below(K+1)
1876		! Change to ...
1877		! if (min(vol_below(K+1) + Vol_err, vol_below(K)) <= Vol_direct) then ?
1878	2142805	if (vol_below(K+1) + Vol_err <= Vol_direct) then
1879	98556	L0 = L_direct ; Vol_0 = Vol_direct
1880		else
1881	2044249	L0 = L(K+1) ; Vol_0 = vol_below(K+1) + Vol_err
1882		! Change to Vol_0 = min(vol_below(K+1) + Vol_err, vol_below(K)) ?
1883		endif
1884
1885		! Try a relatively simple solution that usually works well
1886		! for massless layers.
1887	2142805	dV_dL2 = 0.5(slope+crv) - crvL0 ; dVol = (vol_below(K)-Vol_0)
1888		! dV_dL2 = 0.5(slope+crv) - crvL0 ; dVol = max(vol_below(K)-Vol_0, 0.0)
1889
1890	2142805	use_L0 = .false.
1891	2142805	do_one_L_iter = .false.
1892	2142805	if (CS%answer_date < 20190101) then
1893		curv_tol = GV%Angstrom_ZdV_dL2*2 &
1894	0	* (0.25 * dV_dL2 * GV%Angstrom_Z - crv * L0 * dVol)
1895	0	do_one_L_iter = (crv * crv * dVol**3) < curv_tol
1896		else
1897		! The following code is more robust when GV%Angstrom_H=0, but
1898		! it changes answers.
1899	2142805	use_L0 = (dVol <= 0.)
1900
1901	2142805	Vol_tol = max(0.5 * GV%Angstrom_Z + GV%dZ_subroundoff, 1e-14 * vol_below(K))
1902	2142805	Vol_quit = max(0.9 * GV%Angstrom_Z + GV%dZ_subroundoff, 1e-14 * vol_below(K))
1903
1904		curv_tol = Vol_tol * dV_dL2**2 &
1905	2142805	* (dV_dL2 * Vol_tol - 2.0 * crv * L0 * dVol)
1906	2142805	do_one_L_iter = (crv * crv * dVol**3) < curv_tol
1907		endif
1908
1909	2142805	if (use_L0) then
1910	0	L(K) = L0
1911	0	Vol_err = 0.5(L(K)L(K))(slope + crv_3(3.0-4.0*L(K))) - vol_below(K)
1912	2142805	elseif (do_one_L_iter) then
1913		! One iteration of Newton's method should give an estimate
1914		! that is accurate to within Vol_tol.
1915	0	L(K) = sqrt(L0*L0 + dVol / dV_dL2)
1916	0	Vol_err = 0.5(L(K)L(K))(slope + crv_3(3.0-4.0*L(K))) - vol_below(K)
1917		else
1918	2142805	if (dV_dL2(1.0-L0L0) < dVol + &
1919		dV_dL2 * (Vol_open - vol_below(K))*Ibma_2) then
1920	101353	L_max = sqrt(1.0 - (Vol_open - vol_below(K))*Ibma_2)
1921		else
1922	2041452	L_max = sqrt(L0*L0 + dVol / dV_dL2)
1923		endif
1924	2142805	L_min = sqrt(L0L0 + dVol / (0.5(slope+crv) - crv*L_max))
1925
1926	2142805	Vol_err_min = 0.5(L_min2)(slope + crv_3(3.0-4.0L_min)) - vol_below(K)
1927	2142805	Vol_err_max = 0.5(L_max2)(slope + crv_3(3.0-4.0L_max)) - vol_below(K)
1928		! if ((abs(Vol_err_min) <= Vol_quit) .or. (Vol_err_min >= Vol_err_max)) then
1929	2142805	if (abs(Vol_err_min) <= Vol_quit) then
1930	0	L(K) = L_min ; Vol_err = Vol_err_min
1931		else
1932		L(K) = sqrt((L_min*2Vol_err_max - L_max*2Vol_err_min) / &
1933	2142805	(Vol_err_max - Vol_err_min))
1934	2807612	do itt=1,maxitt
1935	2807612	Vol_err = 0.5(L(K)L(K))(slope + crv_3(3.0-4.0*L(K))) - vol_below(K)
1936	2807612	if (abs(Vol_err) <= Vol_quit) exit
1937		! Take a Newton's method iteration. This equation has proven
1938		! robust enough not to need bracketing.
1939	664807	L(K) = L(K) - Vol_err / (L(K)* (slope + crv - 2.0crvL(K)))
1940		! This would be a Newton's method iteration for L^2:
1941		! L(K) = sqrt(L(K)L(K) - Vol_err / (0.5(slope+crv) - crv*L(K)))
1942		enddo
1943		endif ! end of iterative solver
1944		endif ! end of 1-boundary alternatives.
1945		endif ! end of 0, 1- and 2- boundary cases.
1946		enddo ! k loop to determine L(K) in the convex case
1947
1948	105672	end subroutine find_L_open_convex
1949
1950		!> This subroutine finds a thickness-weighted value of v at the u-points.
1951	2867989	pure function set_v_at_u(v, h, G, GV, i, j, k, mask2dCv, OBC)
1952		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1953		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1954		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1955		intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
1956		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1957		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1958		integer, intent(in) :: i !< The i-index of the u-location to work on.
1959		integer, intent(in) :: j !< The j-index of the u-location to work on.
1960		integer, intent(in) :: k !< The k-index of the u-location to work on.
1961		real, dimension(SZI_(G),SZJB_(G)),&
1962		intent(in) :: mask2dCv !< A multiplicative mask of the v-points [nondim]
1963		type(ocean_OBC_type), pointer :: OBC !< A pointer to an open boundary condition structure
1964		real :: set_v_at_u !< The return value of v at u points points in the
1965		!! same units as u, i.e. [L T-1 ~> m s-1] or other units.
1966		!$omp declare target
1967
1968		! This subroutine finds a thickness-weighted value of v at the u-points.
1969		real :: hwt(0:1,-1:0) ! Masked weights used to average u onto v [H ~> m or kg m-2].
1970		real :: hwt_tot ! The sum of the masked thicknesses [H ~> m or kg m-2].
1971		integer :: i0, j0, i1, j1
1972
1973	20075923	do j0 = -1,0 ; do i0 = 0,1 ; i1 = i+i0 ; J1 = J+j0
1974	17207934	hwt(i0,j0) = (h(i1,j1,k) + h(i1,j1+1,k)) * mask2dCv(i1,J1)
1975		enddo ; enddo
1976
1977	2867989	if (associated(OBC)) then ; if (OBC%number_of_segments > 0) then
1978	0	do j0 = -1,0 ; do i0 = 0,1 ; if (OBC%segnum_v(i+i0,J+j0) /= 0) then
1979	0	i1 = i+i0 ; J1 = J+j0
1980	0	if (OBC%segnum_v(i1,j1) > 0) then ! OBC_DIRECTION_N
1981	0	hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dCv(i1,J1)
1982	0	elseif (OBC%segnum_v(i1,J1) < 0) then ! OBC_DIRECTION_S
1983	0	hwt(i0,j0) = 2.0 * h(i1,J1+1,k) * mask2dCv(i1,J1)
1984		endif
1985		endif ; enddo ; enddo
1986		endif ; endif
1987
1988	2867989	hwt_tot = (hwt(0,-1) + hwt(1,0)) + (hwt(1,-1) + hwt(0,0))
1989	2867989	set_v_at_u = 0.0
1990	2867989	if (hwt_tot > 0.0) set_v_at_u = &
1991		(((hwt(0,0) * v(i,J,k)) + (hwt(1,-1) * v(i+1,J-1,k))) + &
1992	2867989	((hwt(1,0) * v(i+1,J,k)) + (hwt(0,-1) * v(i,J-1,k)))) / hwt_tot
1993
1994	2867989	end function set_v_at_u
1995
1996		!> This subroutine finds a thickness-weighted value of u at the v-points.
1997	2803353	pure function set_u_at_v(u, h, G, GV, i, j, k, mask2dCu, OBC)
1998		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1999		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
2000		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
2001		intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1] or other units.
2002		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2003		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
2004		integer, intent(in) :: i !< The i-index of the u-location to work on.
2005		integer, intent(in) :: j !< The j-index of the u-location to work on.
2006		integer, intent(in) :: k !< The k-index of the u-location to work on.
2007		real, dimension(SZIB_(G),SZJ_(G)), &
2008		intent(in) :: mask2dCu !< A multiplicative mask of the u-points [nondim]
2009		type(ocean_OBC_type), pointer :: OBC !< A pointer to an open boundary condition structure
2010		real :: set_u_at_v !< The return value of u at v points in the
2011		!! same units as u, i.e. [L T-1 ~> m s-1] or other units.
2012		!$omp declare target
2013
2014		! This subroutine finds a thickness-weighted value of u at the v-points.
2015		real :: hwt(-1:0,0:1) ! Masked weights used to average u onto v [H ~> m or kg m-2].
2016		real :: hwt_tot ! The sum of the masked thicknesses [H ~> m or kg m-2].
2017		integer :: i0, j0, i1, j1
2018
2019	19623471	do j0 = 0,1 ; do i0 = -1,0 ; I1 = I+i0 ; j1 = j+j0
2020	16820118	hwt(i0,j0) = (h(i1,j1,k) + h(i1+1,j1,k)) * mask2dCu(I1,j1)
2021		enddo ; enddo
2022
2023	2803353	if (associated(OBC)) then ; if (OBC%number_of_segments > 0) then
2024	0	do j0 = 0,1 ; do i0 = -1,0 ; if ((OBC%segnum_u(I+i0,j+j0) /= 0)) then
2025	0	I1 = I+i0 ; j1 = j+j0
2026	0	if (OBC%segnum_u(I1,j1) > 0) then ! OBC_DIRECTION_E
2027	0	hwt(i0,j0) = 2.0 * h(I1,j1,k) * mask2dCu(I1,j1)
2028	0	elseif (OBC%segnum_u(I1,j1) < 0) then ! OBC_DIRECTION_W
2029	0	hwt(i0,j0) = 2.0 * h(I1+1,j1,k) * mask2dCu(I1,j1)
2030		endif
2031		endif ; enddo ; enddo
2032		endif ; endif
2033
2034	2803353	hwt_tot = (hwt(-1,0) + hwt(0,1)) + (hwt(0,0) + hwt(-1,1))
2035	2803353	set_u_at_v = 0.0
2036	2803353	if (hwt_tot > 0.0) set_u_at_v = &
2037		(((hwt(0,0) * u(I,j,k)) + (hwt(-1,1) * u(I-1,j+1,k))) + &
2038	2803353	((hwt(-1,0) * u(I-1,j,k)) + (hwt(0,1) * u(I,j+1,k)))) / hwt_tot
2039
2040	2803353	end function set_u_at_v
2041
2042		!> Calculates the thickness of the surface boundary layer for applying an elevated viscosity.
2043		!!
2044		!! A bulk Richardson criterion or the thickness of the topmost NKML layers (with a bulk mixed layer)
2045		!! are currently used. The thicknesses are given in terms of fractional layers, so that this
2046		!! thickness will move as the thickness of the topmost layers change.
2047	24	subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)
2048		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
2049		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2050		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2051		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
2052		intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
2053		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
2054		intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
2055		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2056		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
2057		type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
2058		!! thermodynamic fields. Absent fields have
2059		!! NULL pointers.
2060		type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
2061		type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
2062		!! related fields.
2063		real, intent(in) :: dt !< Time increment [T ~> s].
2064		type(set_visc_CS), intent(inout) :: CS !< The control structure returned by a previous
2065		!! call to set_visc_init.
2066
2067		! Local variables
2068		real, dimension(SZIB_(G)) :: &
2069	48	htot, & ! The total thickness of the layers that are within the
2070		! surface mixed layer [H ~> m or kg m-2].
2071	48	dztot, & ! The distance from the surface to the bottom of the layers that are
2072		! within the surface mixed layer [Z ~> m]
2073	48	Thtot, & ! The integrated temperature of layers that are within the
2074		! surface mixed layer [H C ~> m degC or kg degC m-2].
2075	48	Shtot, & ! The integrated salt of layers that are within the
2076		! surface mixed layer [H S ~> m ppt or kg ppt m-2].
2077	48	SpV_htot, & ! Running sum of thickness times specific volume [H R-1 ~> m4 kg-1 or m]
2078	48	Rhtot, & ! The integrated density of layers that are within the surface mixed layer
2079		! [H R ~> kg m-2 or kg2 m-5]. Rhtot is only used if no
2080		! equation of state is used.
2081	48	uhtot, & ! The depth integrated zonal velocity within the surface
2082		! mixed layer [H L T-1 ~> m2 s-1 or kg m-1 s-1].
2083	48	vhtot, & ! The depth integrated meridional velocity within the surface
2084		! mixed layer [H L T-1 ~> m2 s-1 or kg m-1 s-1].
2085	48	Idecay_len_TKE, & ! The inverse of a turbulence decay length scale [H-1 ~> m-1 or m2 kg-1].
2086	48	dR_dT, & ! Partial derivative of the density at the base of layer nkml
2087		! (roughly the base of the mixed layer) with temperature [R C-1 ~> kg m-3 degC-1].
2088	48	dR_dS, & ! Partial derivative of the density at the base of layer nkml
2089		! (roughly the base of the mixed layer) with salinity [R S-1 ~> kg m-3 ppt-1].
2090	48	dSpV_dT, & ! Partial derivative of the specific volume at the base of layer nkml
2091		! (roughly the base of the mixed layer) with temperature [R-1 C-1 ~> m3 kg-1 degC-1].
2092	48	dSpV_dS, & ! Partial derivative of the specific volume at the base of layer nkml
2093		! (roughly the base of the mixed layer) with salinity [R-1 S-1 ~> m3 kg-1 ppt-1].
2094	48	ustar, & ! The surface friction velocity under ice shelves [H T-1 ~> m s-1 or kg m-2 s-1].
2095	48	press, & ! The pressure at which dR_dT and dR_dS are evaluated [R L2 T-2 ~> Pa].
2096	48	T_EOS, & ! The potential temperature at which dR_dT and dR_dS are evaluated [C ~> degC]
2097	48	S_EOS ! The salinity at which dR_dT and dR_dS are evaluated [S ~> ppt].
2098	48	real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! Height change across layers [Z ~> m]
2099		real, dimension(SZIB_(G),SZJ_(G)) :: &
2100	48	mask_u ! A mask that disables any contributions from u points that
2101		! are land or past open boundary conditions [nondim], 0 or 1.
2102		real, dimension(SZI_(G),SZJB_(G)) :: &
2103	48	mask_v ! A mask that disables any contributions from v points that
2104		! are land or past open boundary conditions [nondim], 0 or 1.
2105	48	real :: U_star_2d(SZI_(G),SZJ_(G)) ! The wind friction velocity in thickness-based units,
2106		! calculated using the Boussinesq reference density or the time-evolving
2107		! surface density in non-Boussinesq mode [H T-1 ~> m s-1 or kg m-2 s-1]
2108	48	real :: h_at_vel(SZIB_(G),SZK_(GV))! Layer thickness at velocity points,
2109		! using an upwind-biased second order accurate estimate based
2110		! on the previous velocity direction [H ~> m or kg m-2].
2111	48	real :: dz_at_vel(SZIB_(G),SZK_(GV)) ! Vertical extent of a layer at velocity points,
2112		! using an upwind-biased second order accurate estimate based
2113		! on the previous velocity direction [Z ~> m].
2114	48	integer :: k_massive(SZIB_(G)) ! The k-index of the deepest layer yet found
2115		! that has more than h_tiny thickness and will be in the
2116		! viscous mixed layer.
2117		real :: Uh2 ! The squared magnitude of the difference between the velocity
2118		! integrated through the mixed layer and the velocity of the
2119		! interior layer layer times the depth of the mixed layer
2120		! [H2 L2 T-2 ~> m4 s-2 or kg2 m-2 s-2].
2121		real :: htot_vel ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
2122		real :: hwtot ! Sum of the thicknesses used to calculate
2123		! the near-bottom velocity magnitude [H ~> m or kg m-2].
2124		real :: hutot ! Running sum of thicknesses times the velocity
2125		! magnitudes [H L T-1 ~> m2 s-1 or kg m-1 s-1].
2126		real :: hweight ! The thickness of a layer that is within Hbbl
2127		! of the bottom [H ~> m or kg m-2].
2128		real :: tbl_thick ! The thickness of the top boundary layer [Z ~> m].
2129
2130		real :: hlay ! The layer thickness at velocity points [H ~> m or kg m-2].
2131		real :: I_2hlay ! 1 / 2*hlay [H-1 ~> m-1 or m2 kg-1].
2132		real :: T_lay ! The layer temperature at velocity points [C ~> degC].
2133		real :: S_lay ! The layer salinity at velocity points [S ~> ppt].
2134		real :: Rlay ! The layer potential density at velocity points [R ~> kg m-3].
2135		real :: Rlb ! The potential density of the layer below [R ~> kg m-3].
2136		real :: v_at_u ! The meridional velocity at a zonal velocity point [L T-1 ~> m s-1].
2137		real :: u_at_v ! The zonal velocity at a meridional velocity point [L T-1 ~> m s-1].
2138		real :: gHprime ! The mixed-layer internal gravity wave speed squared, based
2139		! on the mixed layer thickness and density difference across
2140		! the base of the mixed layer [L2 T-2 ~> m2 s-2].
2141		real :: RiBulk ! The bulk Richardson number below which water is in the
2142		! viscous mixed layer, including reduction for turbulent decay [nondim]
2143		real :: dt_Rho0 ! The time step divided by the conversion from the layer
2144		! thickness to layer mass [T H Z-1 R-1 ~> s m3 kg-1 or s].
2145		real :: g_H_Rho0 ! The gravitational acceleration times the conversion from H to m divided
2146		! by the mean density [L2 T-2 H-1 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2].
2147		real :: ustarsq ! 400 times the square of ustar, times
2148		! Rho0 divided by G_Earth and the conversion
2149		! from m to thickness units [H R ~> kg m-2 or kg2 m-5].
2150		real :: cdrag_sqrt ! Square root of the drag coefficient [nondim].
2151		real :: cdrag_sqrt_H ! Square root of the drag coefficient, times a unit conversion
2152		! factor from lateral lengths to layer thicknesses [H L-1 ~> nondim or kg m-3].
2153		real :: cdrag_sqrt_H_RL ! Square root of the drag coefficient, times a unit conversion factor from
2154		! density times lateral lengths to layer thicknesses [H L-1 R-1 ~> m3 kg-1 or nondim]
2155		real :: oldfn ! The integrated energy required to
2156		! entrain up to the bottom of the layer,
2157		! divided by G_Earth [H R ~> kg m-2 or kg2 m-5].
2158		real :: Dfn ! The increment in oldfn for entraining
2159		! the layer [H R ~> kg m-2 or kg2 m-5].
2160		real :: frac_used ! The fraction of the present layer that contributes to Dh and Ddz [nondim]
2161		real :: Dh ! The increment in layer thickness from the present layer [H ~> m or kg m-2].
2162		real :: Ddz ! The increment in height change from the present layer [Z ~> m].
2163	48	real :: u2_bg(SZIB_(G)) ! The square of an assumed background velocity, for
2164		! calculating the mean magnitude near the top for use in
2165		! the quadratic surface drag [L2 T-2 ~> m2 s-2].
2166		real :: h_tiny ! A very small thickness [H ~> m or kg m-2]. Layers that are less than
2167		! h_tiny can not be the deepest in the viscous mixed layer.
2168		real :: absf ! The absolute value of f averaged to velocity points [T-1 ~> s-1].
2169		real :: U_star ! The friction velocity at velocity points [H T-1 ~> m s-1 or kg m-2 s-1].
2170		real :: h_neglect ! A thickness that is so small it is usually lost
2171		! in roundoff and can be neglected [H ~> m or kg m-2].
2172		real :: dz_neglect ! A vertical distance that is so small it is usually lost
2173		! in roundoff and can be neglected [Z ~> m].
2174		real :: Rho0x400_G ! 400*Rho0/G_Earth, times unit conversion factors
2175		! [R T2 H-1 ~> kg s2 m-4 or s2 m-1].
2176		! The 400 is a constant proposed by Killworth and Edwards, 1999.
2177		real :: ustar1 ! ustar [H T-1 ~> m s-1 or kg m-2 s-1]
2178		real :: h2f2 ! (h2f)^2 [H2 T-2 ~> m2 s-2 or kg2 m-4 s-2]
2179	24	logical :: use_EOS, do_any, do_any_shelf, do_i(SZIB_(G))
2180		logical :: nonBous_ML ! If true, use the non-Boussinesq form of some energy and
2181		! stratification calculations.
2182		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, K2, nkmb, nkml, n
2183		type(ocean_OBC_type), pointer :: OBC => NULL()
2184
2185	24	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
2186	24	Isq = G%isc-1 ; Ieq = G%IecB ; Jsq = G%jsc-1 ; Jeq = G%JecB
2187	24	nkmb = GV%nk_rho_varies ; nkml = GV%nkml
2188
2189	24	if (.not.CS%initialized) call MOM_error(FATAL,"MOM_set_viscosity(visc_ML): "//&
2190	0	"Module must be initialized before it is used.")
2191
2192		! TODO: Remove this check and move it outside of the function call.
2193	24	if (.not.(CS%dynamic_viscous_ML .or. associated(forces%frac_shelf_u) .or. &
2194	24	associated(forces%frac_shelf_v)) ) return
2195
2196		! NOTE: Requried since this is called by the GPU-enabled dycore, but it could
2197		! also be implicitly fixing other functions.
2198		!$omp target update from(u, v)
2199
2200	0	Rho0x400_G = 400.0*(GV%H_to_RZ / GV%g_Earth_Z_T2)
2201	0	cdrag_sqrt = sqrt(CS%cdrag)
2202	0	cdrag_sqrt_H = cdrag_sqrt * US%L_to_m * GV%m_to_H
2203	0	cdrag_sqrt_H_RL = cdrag_sqrt * US%L_to_Z * GV%RZ_to_H
2204
2205	0	OBC => CS%OBC
2206	0	use_EOS = associated(tv%eqn_of_state)
2207	0	nonBous_ML = allocated(tv%SpV_avg)
2208	0	dt_Rho0 = dt / GV%H_to_RZ
2209	0	h_neglect = GV%H_subroundoff
2210	0	h_tiny = 2.0*GV%Angstrom_H + h_neglect
2211	0	dz_neglect = GV%dZ_subroundoff
2212	0	g_H_Rho0 = (GV%g_Earth*GV%H_to_Z) / (GV%Rho0)
2213
2214	0	if (associated(forces%frac_shelf_u) .neqv. associated(forces%frac_shelf_v)) &
2215		call MOM_error(FATAL, "set_viscous_ML: one of forces%frac_shelf_u and "//&
2216	0	"forces%frac_shelf_v is associated, but the other is not.")
2217
2218		! Extract the friction velocity from the forcing type.
2219	0	call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1, H_T_units=.true.)
2220
2221	0	if (associated(forces%frac_shelf_u)) then
2222		! This configuration has ice shelves, and the appropriate variables need to be
2223		! allocated. If the arrays have already been allocated, these calls do nothing.
2224	0	if (.not.allocated(visc%taux_shelf)) &
2225	0	allocate(visc%taux_shelf(G%IsdB:G%IedB, G%jsd:G%jed), source=0.0)
2226	0	if (.not.allocated(visc%tauy_shelf)) &
2227	0	allocate(visc%tauy_shelf(G%isd:G%ied, G%JsdB:G%JedB), source=0.0)
2228	0	if (.not.allocated(visc%tbl_thick_shelf_u)) &
2229	0	allocate(visc%tbl_thick_shelf_u(G%IsdB:G%IedB, G%jsd:G%jed), source=0.0)
2230	0	if (.not.allocated(visc%tbl_thick_shelf_v)) &
2231	0	allocate(visc%tbl_thick_shelf_v(G%isd:G%ied, G%JsdB:G%JedB), source=0.0)
2232	0	if (.not.allocated(visc%kv_tbl_shelf_u)) &
2233	0	allocate(visc%kv_tbl_shelf_u(G%IsdB:G%IedB, G%jsd:G%jed), source=0.0)
2234	0	if (.not.allocated(visc%kv_tbl_shelf_v)) &
2235	0	allocate(visc%kv_tbl_shelf_v(G%isd:G%ied, G%JsdB:G%JedB), source=0.0)
2236
2237		! With a linear drag law under shelves, the friction velocity is already known.
2238		! if (CS%linear_drag) ustar(:) = cdrag_sqrt_H*CS%drag_bg_vel
2239
2240		! Find the vertical distances across layers.
2241	0	call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
2242		endif
2243
2244		!$OMP parallel do default(shared)
2245	0	do J=js-1,je ; do i=is-1,ie+1
2246	0	mask_v(i,J) = G%mask2dCv(i,J)
2247		enddo ; enddo
2248		!$OMP parallel do default(shared)
2249	0	do j=js-1,je+1 ; do I=is-1,ie
2250	0	mask_u(I,j) = G%mask2dCu(I,j)
2251		enddo ; enddo
2252
2253	0	if (associated(OBC)) then ; do n=1,OBC%number_of_segments
2254		! Project bottom depths across cell-corner points in the OBCs.
2255	0	if (.not. OBC%segment(n)%on_pe) cycle
2256		! Use a one-sided projection of bottom depths at OBC points.
2257	0	I = OBC%segment(n)%HI%IsdB ; J = OBC%segment(n)%HI%JsdB
2258	0	if (OBC%segment(n)%is_N_or_S .and. (J >= js-1) .and. (J <= je)) then
2259	0	do I = max(is-1,OBC%segment(n)%HI%IsdB), min(ie,OBC%segment(n)%HI%IedB)
2260	0	if (OBC%segment(n)%direction == OBC_DIRECTION_N) mask_u(I,j+1) = 0.0
2261	0	if (OBC%segment(n)%direction == OBC_DIRECTION_S) mask_u(I,j) = 0.0
2262		enddo
2263	0	elseif (OBC%segment(n)%is_E_or_W .and. (I >= is-1) .and. (I <= ie)) then
2264	0	do J = max(js-1,OBC%segment(n)%HI%JsdB), min(je,OBC%segment(n)%HI%JedB)
2265	0	if (OBC%segment(n)%direction == OBC_DIRECTION_E) mask_v(i+1,J) = 0.0
2266	0	if (OBC%segment(n)%direction == OBC_DIRECTION_W) mask_v(i,J) = 0.0
2267		enddo
2268		endif
2269		enddo ; endif
2270
2271		!$OMP parallel do default(private) shared(u,v,h,dz,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, &
2272		!$OMP nonBous_ML,h_neglect,dz_neglect,h_tiny,g_H_Rho0, &
2273		!$OMP js,je,OBC,Isq,Ieq,nz,nkml,U_star_2d,mask_v, &
2274		!$OMP cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL,Rho0x400_G)
2275	0	do j=js,je ! u-point loop
2276	0	if (CS%dynamic_viscous_ML) then
2277	0	do_any = .false.
2278	0	do I=Isq,Ieq
2279	0	htot(I) = 0.0
2280	0	if (G%mask2dCu(I,j) < 0.5) then
2281	0	do_i(I) = .false. ; visc%nkml_visc_u(I,j) = nkml
2282		else
2283	0	do_i(I) = .true. ; do_any = .true.
2284	0	k_massive(I) = nkml
2285	0	Thtot(I) = 0.0 ; Shtot(I) = 0.0 ; Rhtot(i) = 0.0
2286	0	uhtot(I) = dt_Rho0 * forces%taux(I,j)
2287		vhtot(I) = 0.25 * dt_Rho0 * ((forces%tauy(i,J) + forces%tauy(i+1,J-1)) + &
2288	0	(forces%tauy(i,J-1) + forces%tauy(i+1,J)))
2289
2290	0	if (CS%omega_frac >= 1.0) then ; absf = 2.0*CS%omega ; else
2291	0	absf = 0.5*(abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I,J-1)))
2292	0	if (CS%omega_frac > 0.0) &
2293	0	absf = sqrt(CS%omega_frac4.0CS%omega*2 + (1.0-CS%omega_frac)absf**2)
2294		endif
2295	0	U_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i+1,j)))
2296	0	Idecay_len_TKE(I) = (absf / U_star) * CS%TKE_decay
2297		endif
2298		enddo
2299
2300	0	if (do_any) then ; do k=1,nz
2301	0	if (k > nkml) then
2302	0	do_any = .false.
2303	0	if (use_EOS .and. (k==nkml+1)) then
2304		! Find dRho/dT and dRho_dS.
2305	0	do I=Isq,Ieq
2306	0	press(I) = (GV%H_to_RZGV%g_Earth) htot(I)
2307	0	if (associated(tv%p_surf)) press(I) = press(I) + 0.5*(tv%p_surf(i,j)+tv%p_surf(i+1,j))
2308	0	k2 = max(1,nkml)
2309	0	I_2hlay = 1.0 / (h(i,j,k2) + h(i+1,j,k2) + h_neglect)
2310	0	T_EOS(I) = ((h(i,j,k2)tv%T(i,j,k2)) + (h(i+1,j,k2)tv%T(i+1,j,k2))) * I_2hlay
2311	0	S_EOS(I) = ((h(i,j,k2)tv%S(i,j,k2)) + (h(i+1,j,k2)tv%S(i+1,j,k2))) * I_2hlay
2312		enddo
2313		call calculate_density_derivs(T_EOS, S_EOS, press, dR_dT, dR_dS, tv%eqn_of_state, &
2314	0	(/Isq-G%IsdB+1,Ieq-G%IsdB+1/) )
2315	0	if (nonBous_ML) then
2316		call calculate_specific_vol_derivs(T_EOS, S_EOS, press, dSpV_dT, dSpV_dS, tv%eqn_of_state, &
2317	0	(/Isq-G%IsdB+1,Ieq-G%IsdB+1/) )
2318		endif
2319		endif
2320
2321	0	do I=Isq,Ieq ; if (do_i(I)) then
2322
2323	0	hlay = 0.5*(h(i,j,k) + h(i+1,j,k))
2324	0	if (hlay > h_tiny) then ! Only consider non-vanished layers.
2325	0	I_2hlay = 1.0 / (h(i,j,k) + h(i+1,j,k))
2326		v_at_u = 0.5 * ((h(i,j,k) * (v(i,J,k) + v(i,J-1,k))) + &
2327	0	(h(i+1,j,k) * (v(i+1,J,k) + v(i+1,J-1,k)))) * I_2hlay
2328	0	Uh2 = (uhtot(I) - htot(I)u(I,j,k))2 + (vhtot(I) - htot(I)v_at_u)**2
2329
2330	0	if (use_EOS) then
2331	0	T_lay = ((h(i,j,k)tv%T(i,j,k)) + (h(i+1,j,k)tv%T(i+1,j,k))) * I_2hlay
2332	0	S_lay = ((h(i,j,k)tv%S(i,j,k)) + (h(i+1,j,k)tv%S(i+1,j,k))) * I_2hlay
2333	0	if (nonBous_ML) then
2334		gHprime = (GV%g_Earth * GV%H_to_RZ) * (dSpV_dT(I) * (Thtot(I) - T_lay*htot(I)) + &
2335	0	dSpV_dS(I) * (Shtot(I) - S_lay*htot(I)))
2336		else
2337		gHprime = g_H_Rho0 * (dR_dT(I) * (T_lay*htot(I) - Thtot(I)) + &
2338	0	dR_dS(I) * (S_lay*htot(I) - Shtot(I)))
2339		endif
2340		else
2341	0	gHprime = g_H_Rho0 * (GV%Rlay(k)*htot(I) - Rhtot(I))
2342		endif
2343
2344	0	if (gHprime > 0.0) then
2345	0	RiBulk = CS%bulk_Ri_ML * exp(-htot(I) * Idecay_len_TKE(I))
2346	0	if (RiBulk * Uh2 <= (htot(I)*2) gHprime) then
2347	0	visc%nkml_visc_u(I,j) = real(k_massive(I))
2348	0	do_i(I) = .false.
2349	0	elseif (RiBulk * Uh2 <= (htot(I) + hlay)*2 gHprime) then
2350		visc%nkml_visc_u(I,j) = real(k-1) + &
2351	0	( sqrt(RiBulk * Uh2 / gHprime) - htot(I) ) / hlay
2352	0	do_i(I) = .false.
2353		endif
2354		endif
2355	0	k_massive(I) = k
2356		endif ! hlay > h_tiny
2357
2358	0	if (do_i(I)) do_any = .true.
2359		endif ; enddo
2360
2361	0	if (.not.do_any) exit ! All columns are done.
2362		endif
2363
2364	0	do I=Isq,Ieq ; if (do_i(I)) then
2365	0	htot(I) = htot(I) + 0.5 * (h(i,j,k) + h(i+1,j,k))
2366	0	uhtot(I) = uhtot(I) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * u(I,j,k)
2367		vhtot(I) = vhtot(I) + 0.25 * ((h(i,j,k) * (v(i,J,k) + v(i,J-1,k))) + &
2368	0	(h(i+1,j,k) * (v(i+1,J,k) + v(i+1,J-1,k))))
2369	0	if (use_EOS) then
2370	0	Thtot(I) = Thtot(I) + 0.5 * ((h(i,j,k)tv%T(i,j,k)) + (h(i+1,j,k)tv%T(i+1,j,k)))
2371	0	Shtot(I) = Shtot(I) + 0.5 * ((h(i,j,k)tv%S(i,j,k)) + (h(i+1,j,k)tv%S(i+1,j,k)))
2372		else
2373	0	Rhtot(i) = Rhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * GV%Rlay(k)
2374		endif
2375		endif ; enddo
2376		enddo ; endif
2377
2378	0	if (do_any) then ; do I=Isq,Ieq ; if (do_i(I)) then
2379	0	visc%nkml_visc_u(I,j) = k_massive(I)
2380		endif ; enddo ; endif
2381		endif ! dynamic_viscous_ML
2382
2383	0	do_any_shelf = .false.
2384	0	if (associated(forces%frac_shelf_u)) then
2385	0	do I=Isq,Ieq
2386	0	if (forces%frac_shelf_u(I,j)*G%mask2dCu(I,j) == 0.0) then
2387	0	do_i(I) = .false.
2388	0	visc%tbl_thick_shelf_u(I,j) = 0.0 ; visc%kv_tbl_shelf_u(I,j) = 0.0
2389		else
2390	0	do_i(I) = .true. ; do_any_shelf = .true.
2391		endif
2392		enddo
2393		endif
2394
2395	0	if (do_any_shelf) then
2396	0	do k=1,nz ; do I=Isq,Ieq ; if (do_i(I)) then
2397	0	if (u(I,j,k) * (h(i+1,j,k) - h(i,j,k)) >= 0) then
2398		h_at_vel(i,k) = 2.0h(i,j,k)h(i+1,j,k) / &
2399	0	(h(i,j,k) + h(i+1,j,k) + h_neglect)
2400		dz_at_vel(i,k) = 2.0dz(i,j,k)dz(i+1,j,k) / &
2401	0	(dz(i,j,k) + dz(i+1,j,k) + dz_neglect)
2402		else
2403	0	h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
2404	0	dz_at_vel(i,k) = 0.5 * (dz(i,j,k) + dz(i+1,j,k))
2405		endif
2406		else
2407	0	h_at_vel(I,k) = 0.0
2408	0	dz_at_vel(I,k) = 0.0
2409	0	ustar(I) = 0.0
2410		endif ; enddo ; enddo
2411
2412	0	do I=Isq,Ieq ; if (do_i(I)) then
2413	0	htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
2414	0	Thtot(I) = 0.0 ; Shtot(I) = 0.0 ; SpV_htot(I) = 0.0
2415	0	if (use_EOS .or. .not.CS%linear_drag) then ; do k=1,nz
2416	0	if (htot_vel>=CS%Htbl_shelf) exit ! terminate the k loop
2417	0	hweight = MIN(CS%Htbl_shelf - htot_vel, h_at_vel(i,k))
2418	0	if (hweight <= 1.5*GV%Angstrom_H + h_neglect) cycle
2419
2420	0	htot_vel = htot_vel + h_at_vel(i,k)
2421	0	hwtot = hwtot + hweight
2422
2423	0	if (.not.CS%linear_drag) then
2424	0	v_at_u = set_v_at_u(v, h, G, GV, i, j, k, mask_v, OBC)
2425		! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant
2426	0	if (CS%BBL_use_tidal_bg) then
2427		u2_bg(I) = 0.5( G%mask2dT(i,j)(CS%tideamp(i,j)*CS%tideamp(i,j))+ &
2428	0	G%mask2dT(i+1,j)(CS%tideamp(i+1,j)CS%tideamp(i+1,j)) )
2429		else
2430	0	u2_bg(I) = CS%drag_bg_vel * CS%drag_bg_vel
2431		endif
2432	0	hutot = hutot + hweight * sqrt(u(I,j,k)2 + v_at_u2 + u2_bg(I))
2433		endif
2434	0	if (use_EOS) then
2435	0	Thtot(I) = Thtot(I) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
2436	0	Shtot(I) = Shtot(I) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
2437		endif
2438	0	if (allocated(tv%SpV_avg)) then
2439	0	SpV_htot(I) = SpV_htot(I) + hweight * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i+1,j,k))
2440		endif
2441		enddo ; endif
2442
2443	0	if ((hwtot <= 0.0) .or. (CS%linear_drag .and. .not.allocated(tv%SpV_avg))) then
2444	0	ustar(I) = cdrag_sqrt_H * CS%drag_bg_vel
2445	0	elseif (CS%linear_drag .and. allocated(tv%SpV_avg)) then
2446	0	ustar(I) = cdrag_sqrt_H_RL * CS%drag_bg_vel * (hwtot / SpV_htot(I))
2447	0	elseif (allocated(tv%SpV_avg)) then ! (.not.CS%linear_drag)
2448	0	ustar(I) = cdrag_sqrt_H_RL * hutot / SpV_htot(I)
2449		else ! (.not.CS%linear_drag .and. .not.allocated(tv%SpV_avg))
2450	0	ustar(I) = cdrag_sqrt_H * hutot / hwtot
2451		endif
2452
2453	0	if (use_EOS) then ; if (hwtot > 0.0) then
2454	0	T_EOS(I) = Thtot(I)/hwtot ; S_EOS(I) = Shtot(I)/hwtot
2455		else
2456	0	T_EOS(I) = 0.0 ; S_EOS(I) = 0.0
2457		endif ; endif
2458		! if (allocated(tv%SpV_avg)) SpV_av(I) = SpVhtot(I) / hwtot
2459		endif ; enddo ! I-loop
2460
2461	0	if (use_EOS) then
2462		call calculate_density_derivs(T_EOS, S_EOS, forces%p_surf(:,j), dR_dT, dR_dS, &
2463	0	tv%eqn_of_state, (/Isq-G%IsdB+1,Ieq-G%IsdB+1/) )
2464		endif
2465
2466	0	do I=Isq,Ieq ; if (do_i(I)) then
2467		! The 400.0 in this expression is the square of a constant proposed
2468		! by Killworth and Edwards, 1999, in equation (2.20).
2469	0	ustarsq = Rho0x400_G * ustar(i)**2
2470	0	htot(i) = 0.0 ; dztot(i) = 0.0
2471	0	if (use_EOS) then
2472	0	Thtot(i) = 0.0 ; Shtot(i) = 0.0 ; oldfn = 0.0
2473	0	do k=1,nz-1
2474	0	if (h_at_vel(i,k) <= 0.0) cycle
2475	0	T_Lay = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
2476	0	S_Lay = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
2477	0	oldfn = dR_dT(i)(T_Layhtot(i) - Thtot(i)) + dR_dS(i)(S_Layhtot(i) - Shtot(i))
2478	0	if (oldfn >= ustarsq) exit
2479
2480		Dfn = (dR_dT(i)(0.5(tv%T(i,j,k+1)+tv%T(i+1,j,k+1)) - T_Lay) + &
2481		dR_dS(i)(0.5(tv%S(i,j,k+1)+tv%S(i+1,j,k+1)) - S_Lay)) * &
2482	0	(h_at_vel(i,k)+htot(i))
2483	0	if ((oldfn + Dfn) <= ustarsq) then
2484	0	Dh = h_at_vel(i,k)
2485	0	Ddz = dz_at_vel(i,k)
2486		else
2487	0	frac_used = sqrt((ustarsq-oldfn) / (Dfn))
2488	0	Dh = h_at_vel(i,k) * frac_used
2489	0	Ddz = dz_at_vel(i,k) * frac_used
2490		endif
2491
2492	0	htot(i) = htot(i) + Dh
2493	0	dztot(i) = dztot(i) + Ddz
2494	0	Thtot(i) = Thtot(i) + T_LayDh ; Shtot(i) = Shtot(i) + S_LayDh
2495		enddo
2496	0	if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0)) then
2497	0	T_Lay = 0.5*(tv%T(i,j,nz) + tv%T(i+1,j,nz))
2498	0	S_Lay = 0.5*(tv%S(i,j,nz) + tv%S(i+1,j,nz))
2499	0	if (dR_dT(i)(T_Layhtot(i) - Thtot(i)) + &
2500		dR_dS(i)(S_Layhtot(i) - Shtot(i)) < ustarsq) then
2501	0	htot(i) = htot(i) + h_at_vel(i,nz)
2502	0	dztot(i) = dztot(i) + dz_at_vel(i,nz)
2503		endif
2504		endif ! Examination of layer nz.
2505		else ! Use Rlay as the density variable.
2506	0	Rhtot = 0.0
2507	0	do k=1,nz-1
2508	0	Rlay = GV%Rlay(k) ; Rlb = GV%Rlay(k+1)
2509
2510	0	oldfn = Rlay*htot(i) - Rhtot(i)
2511	0	if (oldfn >= ustarsq) exit
2512
2513	0	Dfn = (Rlb - Rlay)*(h_at_vel(i,k)+htot(i))
2514	0	if ((oldfn + Dfn) <= ustarsq) then
2515	0	Dh = h_at_vel(i,k)
2516	0	Ddz = dz_at_vel(i,k)
2517		else
2518	0	frac_used = sqrt((ustarsq-oldfn) / (Dfn))
2519	0	Dh = h_at_vel(i,k) * frac_used
2520	0	Ddz = dz_at_vel(i,k) * frac_used
2521		endif
2522
2523	0	htot(i) = htot(i) + Dh
2524	0	dztot(i) = dztot(i) + Ddz
2525	0	Rhtot(i) = Rhtot(i) + Rlay*Dh
2526		enddo
2527	0	if (GV%Rlay(nz)*htot(i) - Rhtot(i) < ustarsq) then
2528	0	htot(i) = htot(i) + h_at_vel(i,nz)
2529	0	dztot(i) = dztot(i) + dz_at_vel(i,nz)
2530		endif
2531		endif ! use_EOS
2532
2533		! visc%tbl_thick_shelf_u(I,j) = max(CS%Htbl_shelf_min, &
2534		! dztot(I) / (0.5 + sqrt(0.25 + &
2535		! ((htot(i)(G%CoriolisBu(I,J-1)+G%CoriolisBu(I,J)))*2) / &
2536		! (ustar(i)**2) )) )
2537	0	ustar1 = ustar(i)
2538	0	h2f2 = (htot(i)(G%CoriolisBu(I,J-1)+G%CoriolisBu(I,J)) + h_neglectCS%omega)**2
2539		tbl_thick = max(CS%Htbl_shelf_min, &
2540	0	( dztot(I)ustar(i) ) / ( 0.5ustar1 + sqrt((0.5ustar1)*2 + h2f2 ) ) )
2541	0	visc%tbl_thick_shelf_u(I,j) = tbl_thick
2542	0	visc%Kv_tbl_shelf_u(I,j) = max(CS%Kv_TBL_min, cdrag_sqrtustar1tbl_thick)
2543		endif ; enddo ! I-loop
2544		endif ! do_any_shelf
2545
2546		enddo ! j-loop at u-points
2547
2548		!$OMP parallel do default(private) shared(u,v,h,dz,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, &
2549		!$OMP nonBous_ML,h_neglect,dz_neglect,h_tiny,g_H_Rho0, &
2550		!$OMP is,ie,OBC,Jsq,Jeq,nz,nkml,U_star_2d,mask_u, &
2551		!$OMP cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL,Rho0x400_G)
2552	0	do J=Jsq,Jeq ! v-point loop
2553	0	if (CS%dynamic_viscous_ML) then
2554	0	do_any = .false.
2555	0	do i=is,ie
2556	0	htot(i) = 0.0
2557	0	if (G%mask2dCv(i,J) < 0.5) then
2558	0	do_i(i) = .false. ; visc%nkml_visc_v(i,J) = nkml
2559		else
2560	0	do_i(i) = .true. ; do_any = .true.
2561	0	k_massive(i) = nkml
2562	0	Thtot(i) = 0.0 ; Shtot(i) = 0.0 ; Rhtot(i) = 0.0
2563	0	vhtot(i) = dt_Rho0 * forces%tauy(i,J)
2564		uhtot(i) = 0.25 * dt_Rho0 * ((forces%taux(I,j) + forces%taux(I-1,j+1)) + &
2565	0	(forces%taux(I-1,j) + forces%taux(I,j+1)))
2566
2567	0	if (CS%omega_frac >= 1.0) then ; absf = 2.0*CS%omega ; else
2568	0	absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))
2569	0	if (CS%omega_frac > 0.0) &
2570	0	absf = sqrt(CS%omega_frac4.0CS%omega*2 + (1.0-CS%omega_frac)absf**2)
2571		endif
2572
2573	0	U_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i,j+1)))
2574	0	Idecay_len_TKE(i) = (absf / U_star) * CS%TKE_decay
2575
2576		endif
2577		enddo
2578
2579	0	if (do_any) then ; do k=1,nz
2580	0	if (k > nkml) then
2581	0	do_any = .false.
2582	0	if (use_EOS .and. (k==nkml+1)) then
2583		! Find dRho/dT and dRho_dS.
2584	0	do i=is,ie
2585	0	press(i) = (GV%H_to_RZ * GV%g_Earth) * htot(i)
2586	0	if (associated(tv%p_surf)) press(i) = press(i) + 0.5*(tv%p_surf(i,j)+tv%p_surf(i,j+1))
2587	0	k2 = max(1,nkml)
2588	0	I_2hlay = 1.0 / (h(i,j,k2) + h(i,j+1,k2) + h_neglect)
2589	0	T_EOS(i) = ((h(i,j,k2)tv%T(i,j,k2)) + (h(i,j+1,k2)tv%T(i,j+1,k2))) * I_2hlay
2590	0	S_EOS(i) = ((h(i,j,k2)tv%S(i,j,k2)) + (h(i,j+1,k2)tv%S(i,j+1,k2))) * I_2hlay
2591		enddo
2592		call calculate_density_derivs(T_EOS, S_EOS, press, dR_dT, dR_dS, &
2593	0	tv%eqn_of_state, (/is-G%IsdB+1,ie-G%IsdB+1/) )
2594	0	if (nonBous_ML) then
2595		call calculate_specific_vol_derivs(T_EOS, S_EOS, press, dSpV_dT, dSpV_dS, tv%eqn_of_state, &
2596	0	(/is-G%IsdB+1,ie-G%IsdB+1/) )
2597		endif
2598		endif
2599
2600	0	do i=is,ie ; if (do_i(i)) then
2601
2602	0	hlay = 0.5*(h(i,j,k) + h(i,j+1,k))
2603	0	if (hlay > h_tiny) then ! Only consider non-vanished layers.
2604	0	I_2hlay = 1.0 / (h(i,j,k) + h(i,j+1,k))
2605		u_at_v = 0.5 * ((h(i,j,k) * (u(I-1,j,k) + u(I,j,k))) + &
2606	0	(h(i,j+1,k) * (u(I-1,j+1,k) + u(I,j+1,k)))) * I_2hlay
2607	0	Uh2 = (vhtot(i) - htot(i)v(i,J,k))2 + (uhtot(i) - htot(i)u_at_v)**2
2608
2609	0	if (use_EOS) then
2610	0	T_lay = ((h(i,j,k)tv%T(i,j,k)) + (h(i,j+1,k)tv%T(i,j+1,k))) * I_2hlay
2611	0	S_lay = ((h(i,j,k)tv%S(i,j,k)) + (h(i,j+1,k)tv%S(i,j+1,k))) * I_2hlay
2612	0	if (nonBous_ML) then
2613		gHprime = (GV%g_Earth * GV%H_to_RZ) * (dSpV_dT(i) * (Thtot(i) - T_lay*htot(i)) + &
2614	0	dSpV_dS(i) * (Shtot(i) - S_lay*htot(i)))
2615		else
2616		gHprime = g_H_Rho0 * (dR_dT(i) * (T_lay*htot(i) - Thtot(i)) + &
2617	0	dR_dS(i) * (S_lay*htot(i) - Shtot(i)))
2618		endif
2619		else
2620	0	gHprime = g_H_Rho0 * (GV%Rlay(k)*htot(i) - Rhtot(i))
2621		endif
2622
2623	0	if (gHprime > 0.0) then
2624	0	RiBulk = CS%bulk_Ri_ML * exp(-htot(i) * Idecay_len_TKE(i))
2625	0	if (RiBulk * Uh2 <= htot(i)*2 gHprime) then
2626	0	visc%nkml_visc_v(i,J) = real(k_massive(i))
2627	0	do_i(i) = .false.
2628	0	elseif (RiBulk * Uh2 <= (htot(i) + hlay)*2 gHprime) then
2629		visc%nkml_visc_v(i,J) = real(k-1) + &
2630	0	( sqrt(RiBulk * Uh2 / gHprime) - htot(i) ) / hlay
2631	0	do_i(i) = .false.
2632		endif
2633		endif
2634	0	k_massive(i) = k
2635		endif ! hlay > h_tiny
2636
2637	0	if (do_i(i)) do_any = .true.
2638		endif ; enddo
2639
2640	0	if (.not.do_any) exit ! All columns are done.
2641		endif
2642
2643	0	do i=is,ie ; if (do_i(i)) then
2644	0	htot(i) = htot(i) + 0.5 * (h(i,J,k) + h(i,j+1,k))
2645	0	vhtot(i) = vhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * v(i,J,k)
2646		uhtot(i) = uhtot(i) + 0.25 * ((h(i,j,k) * (u(I-1,j,k) + u(I,j,k))) + &
2647	0	(h(i,j+1,k) * (u(I-1,j+1,k) + u(I,j+1,k))))
2648	0	if (use_EOS) then
2649	0	Thtot(i) = Thtot(i) + 0.5 * ((h(i,j,k)tv%T(i,j,k)) + (h(i,j+1,k)tv%T(i,j+1,k)))
2650	0	Shtot(i) = Shtot(i) + 0.5 * ((h(i,j,k)tv%S(i,j,k)) + (h(i,j+1,k)tv%S(i,j+1,k)))
2651		else
2652	0	Rhtot(i) = Rhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * GV%Rlay(k)
2653		endif
2654		endif ; enddo
2655		enddo ; endif
2656
2657	0	if (do_any) then ; do i=is,ie ; if (do_i(i)) then
2658	0	visc%nkml_visc_v(i,J) = k_massive(i)
2659		endif ; enddo ; endif
2660
2661		endif ! dynamic_viscous_ML
2662
2663	0	do_any_shelf = .false.
2664	0	if (associated(forces%frac_shelf_v)) then
2665	0	do i=is,ie
2666	0	if (forces%frac_shelf_v(i,J)*G%mask2dCv(i,J) == 0.0) then
2667	0	do_i(i) = .false.
2668	0	visc%tbl_thick_shelf_v(i,J) = 0.0 ; visc%kv_tbl_shelf_v(i,J) = 0.0
2669		else
2670	0	do_i(i) = .true. ; do_any_shelf = .true.
2671		endif
2672		enddo
2673		endif
2674
2675	0	if (do_any_shelf) then
2676	0	do k=1,nz ; do i=is,ie ; if (do_i(i)) then
2677	0	if (v(i,J,k) * (h(i,j+1,k) - h(i,j,k)) >= 0) then
2678		h_at_vel(i,k) = 2.0h(i,j,k)h(i,j+1,k) / &
2679	0	(h(i,j,k) + h(i,j+1,k) + h_neglect)
2680		dz_at_vel(i,k) = 2.0dz(i,j,k)dz(i,j+1,k) / &
2681	0	(dz(i,j,k) + dz(i,j+1,k) + dz_neglect)
2682		else
2683	0	h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
2684	0	dz_at_vel(i,k) = 0.5 * (dz(i,j,k) + dz(i,j+1,k))
2685		endif
2686		else
2687	0	h_at_vel(I,k) = 0.0
2688	0	dz_at_vel(I,k) = 0.0
2689	0	ustar(i) = 0.0
2690		endif ; enddo ; enddo
2691
2692	0	do i=is,ie ; if (do_i(i)) then
2693	0	htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
2694	0	Thtot(i) = 0.0 ; Shtot(i) = 0.0 ; SpV_htot(i) = 0.0
2695	0	if (use_EOS .or. .not.CS%linear_drag) then ; do k=1,nz
2696	0	if (htot_vel>=CS%Htbl_shelf) exit ! terminate the k loop
2697	0	hweight = MIN(CS%Htbl_shelf - htot_vel, h_at_vel(i,k))
2698	0	if (hweight <= 1.5*GV%Angstrom_H + h_neglect) cycle
2699
2700	0	htot_vel = htot_vel + h_at_vel(i,k)
2701	0	hwtot = hwtot + hweight
2702
2703	0	if (.not.CS%linear_drag) then
2704	0	u_at_v = set_u_at_v(u, h, G, GV, i, J, k, mask_u, OBC)
2705		! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant
2706	0	if (CS%BBL_use_tidal_bg) then
2707		u2_bg(i) = 0.5( G%mask2dT(i,j)(CS%tideamp(i,j)*CS%tideamp(i,j))+ &
2708	0	G%mask2dT(i,j+1)(CS%tideamp(i,j+1)CS%tideamp(i,j+1)) )
2709		else
2710	0	u2_bg(i) = CS%drag_bg_vel * CS%drag_bg_vel
2711		endif
2712	0	hutot = hutot + hweight * sqrt(v(i,J,k)2 + u_at_v2 + u2_bg(i))
2713		endif
2714	0	if (use_EOS) then
2715	0	Thtot(i) = Thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
2716	0	Shtot(i) = Shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
2717		endif
2718	0	if (allocated(tv%SpV_avg)) then
2719	0	SpV_htot(i) = SpV_htot(i) + hweight * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j+1,k))
2720		endif
2721		enddo ; endif
2722
2723	0	if ((hwtot <= 0.0) .or. (CS%linear_drag .and. .not.allocated(tv%SpV_avg))) then
2724	0	ustar(i) = cdrag_sqrt_H * CS%drag_bg_vel
2725	0	elseif (CS%linear_drag .and. allocated(tv%SpV_avg)) then
2726	0	ustar(i) = cdrag_sqrt_H_RL * CS%drag_bg_vel * (hwtot / SpV_htot(i))
2727	0	elseif (allocated(tv%SpV_avg)) then ! (.not.CS%linear_drag)
2728	0	ustar(i) = cdrag_sqrt_H_RL * hutot / SpV_htot(i)
2729		else ! (.not.CS%linear_drag .and. .not.allocated(tv%SpV_avg))
2730	0	ustar(i) = cdrag_sqrt_H * hutot / hwtot
2731		endif
2732
2733	0	if (use_EOS) then ; if (hwtot > 0.0) then
2734	0	T_EOS(i) = Thtot(i)/hwtot ; S_EOS(i) = Shtot(i)/hwtot
2735		else
2736	0	T_EOS(i) = 0.0 ; S_EOS(i) = 0.0
2737		endif ; endif
2738		endif ; enddo ! I-loop
2739
2740	0	if (use_EOS) then
2741		call calculate_density_derivs(T_EOS, S_EOS, forces%p_surf(:,j), dR_dT, dR_dS, &
2742	0	tv%eqn_of_state, (/is-G%IsdB+1,ie-G%IsdB+1/) )
2743		endif
2744
2745	0	do i=is,ie ; if (do_i(i)) then
2746		! The 400.0 in this expression is the square of a constant proposed
2747		! by Killworth and Edwards, 1999, in equation (2.20).
2748	0	ustarsq = Rho0x400_G * ustar(i)**2
2749	0	htot(i) = 0.0
2750	0	dztot(i) = 0.0
2751	0	if (use_EOS) then
2752	0	Thtot(i) = 0.0 ; Shtot(i) = 0.0 ; oldfn = 0.0
2753	0	do k=1,nz-1
2754	0	if (h_at_vel(i,k) <= 0.0) cycle
2755	0	T_Lay = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
2756	0	S_Lay = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
2757	0	oldfn = dR_dT(i)(T_Layhtot(i) - Thtot(i)) + dR_dS(i)(S_Layhtot(i) - Shtot(i))
2758	0	if (oldfn >= ustarsq) exit
2759
2760		Dfn = (dR_dT(i)(0.5(tv%T(i,j,k+1)+tv%T(i,j+1,k+1)) - T_Lay) + &
2761		dR_dS(i)(0.5(tv%S(i,j,k+1)+tv%S(i,j+1,k+1)) - S_Lay)) * &
2762	0	(h_at_vel(i,k)+htot(i))
2763	0	if ((oldfn + Dfn) <= ustarsq) then
2764	0	Dh = h_at_vel(i,k)
2765	0	Ddz = dz_at_vel(i,k)
2766		else
2767	0	frac_used = sqrt((ustarsq-oldfn) / (Dfn))
2768	0	Dh = h_at_vel(i,k) * frac_used
2769	0	Ddz = dz_at_vel(i,k) * frac_used
2770		endif
2771
2772	0	htot(i) = htot(i) + Dh
2773	0	dztot(i) = dztot(i) + Ddz
2774	0	Thtot(i) = Thtot(i) + T_LayDh ; Shtot(i) = Shtot(i) + S_LayDh
2775		enddo
2776	0	if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0)) then
2777	0	T_Lay = 0.5*(tv%T(i,j,nz) + tv%T(i,j+1,nz))
2778	0	S_Lay = 0.5*(tv%S(i,j,nz) + tv%S(i,j+1,nz))
2779	0	if (dR_dT(i)(T_Layhtot(i) - Thtot(i)) + &
2780		dR_dS(i)(S_Layhtot(i) - Shtot(i)) < ustarsq) then
2781	0	htot(i) = htot(i) + h_at_vel(i,nz)
2782	0	dztot(i) = dztot(i) + dz_at_vel(i,nz)
2783		endif
2784		endif ! Examination of layer nz.
2785		else ! Use Rlay as the density variable.
2786	0	Rhtot = 0.0
2787	0	do k=1,nz-1
2788	0	Rlay = GV%Rlay(k) ; Rlb = GV%Rlay(k+1)
2789
2790	0	oldfn = Rlay*htot(i) - Rhtot(i)
2791	0	if (oldfn >= ustarsq) exit
2792
2793	0	Dfn = (Rlb - Rlay)*(h_at_vel(i,k)+htot(i))
2794	0	if ((oldfn + Dfn) <= ustarsq) then
2795	0	Dh = h_at_vel(i,k)
2796	0	Ddz = dz_at_vel(i,k)
2797		else
2798	0	frac_used = sqrt((ustarsq-oldfn) / (Dfn))
2799	0	Dh = h_at_vel(i,k) * frac_used
2800	0	Ddz = dz_at_vel(i,k) * frac_used
2801		endif
2802
2803	0	htot(i) = htot(i) + Dh
2804	0	dztot(i) = dztot(i) + Ddz
2805	0	Rhtot = Rhtot + Rlay*Dh
2806		enddo
2807	0	if (GV%Rlay(nz)*htot(i) - Rhtot(i) < ustarsq) then
2808	0	htot(i) = htot(i) + h_at_vel(i,nz)
2809	0	dztot(i) = dztot(i) + dz_at_vel(i,nz)
2810		endif
2811		endif ! use_EOS
2812
2813		! visc%tbl_thick_shelf_v(i,J) = max(CS%Htbl_shelf_min, &
2814		! dztot(i) / (0.5 + sqrt(0.25 + &
2815		! (htot(i)(G%CoriolisBu(I-1,J)+G%CoriolisBu(I,J)))*2 / &
2816		! (ustar(i))**2 )) )
2817	0	ustar1 = ustar(i)
2818	0	h2f2 = (htot(i)(G%CoriolisBu(I-1,J)+G%CoriolisBu(I,J)) + h_neglectCS%omega)**2
2819		tbl_thick = max(CS%Htbl_shelf_min, &
2820	0	( dztot(i)ustar(i) ) / ( 0.5ustar1 + sqrt((0.5ustar1)*2 + h2f2 ) ) )
2821	0	visc%tbl_thick_shelf_v(i,J) = tbl_thick
2822	0	visc%Kv_tbl_shelf_v(i,J) = max(CS%Kv_TBL_min, cdrag_sqrtustar1tbl_thick)
2823
2824		endif ; enddo ! i-loop
2825		endif ! do_any_shelf
2826
2827		enddo ! J-loop at v-points
2828
2829		!$omp target update to(visc%nkml_visc_u, visc%nkml_visc_v) if (CS%dynamic_viscous_ML)
2830
2831	0	if (CS%debug) then
2832	0	if (allocated(visc%nkml_visc_u) .and. allocated(visc%nkml_visc_v)) &
2833		call uvchksum("nkml_visc_[uv]", visc%nkml_visc_u, visc%nkml_visc_v, &
2834	0	G%HI, haloshift=0, scalar_pair=.true.)
2835		endif
2836	0	if (CS%id_nkml_visc_u > 0) call post_data(CS%id_nkml_visc_u, visc%nkml_visc_u, CS%diag)
2837	0	if (CS%id_nkml_visc_v > 0) call post_data(CS%id_nkml_visc_v, visc%nkml_visc_v, CS%diag)
2838
2839		end subroutine set_viscous_ML
2840
2841		!> Register any fields associated with the vertvisc_type.
2842	1	subroutine set_visc_register_restarts(HI, G, GV, US, param_file, visc, restart_CS, use_ice_shelf)
2843		type(hor_index_type), intent(in) :: HI !< A horizontal index type structure.
2844		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2845		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2846		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2847		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
2848		!! parameters.
2849		type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical
2850		!! viscosities and related fields.
2851		!! Allocated here.
2852		type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure
2853		logical, intent(in) :: use_ice_shelf !< if true, register tau_shelf restarts
2854		! Local variables
2855		logical :: use_kappa_shear, KS_at_vertex
2856		logical :: adiabatic, useKPP, useEPBL, use_ideal_age
2857		logical :: do_brine_plume, use_hor_bnd_diff, use_neutral_diffusion, use_fpmix
2858		logical :: use_CVMix_shear, MLE_use_PBL_MLD, MLE_use_Bodner, use_CVMix_conv
2859		integer :: isd, ied, jsd, jed, nz
2860		real :: hfreeze !< If hfreeze > 0 [Z ~> m], melt potential will be computed.
2861		character(len=16) :: Kv_units, Kd_units
2862		character(len=40) :: mdl = "MOM_set_visc" ! This module's name.
2863	11	type(vardesc) :: u_desc, v_desc
2864	1	isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke
2865
2866		call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
2867	1	do_not_log=.true.)
2868
2869	1	use_kappa_shear = .false. ; KS_at_vertex = .false. ; use_CVMix_shear = .false.
2870	1	useKPP = .false. ; useEPBL = .false. ; use_CVMix_conv = .false.
2871
2872	1	if (.not.adiabatic) then
2873	1	use_kappa_shear = kappa_shear_is_used(param_file)
2874	1	KS_at_vertex = kappa_shear_at_vertex(param_file)
2875	1	use_CVMix_shear = CVMix_shear_is_used(param_file)
2876	1	use_CVMix_conv = CVMix_conv_is_used(param_file)
2877		call get_param(param_file, mdl, "USE_KPP", useKPP, &
2878		"If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
2879		"to calculate diffusivities and non-local transport in the OBL.", &
2880	1	default=.false., do_not_log=.true.)
2881		call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", useEPBL, &
2882		"If true, use an implied energetics planetary boundary "//&
2883		"layer scheme to determine the diffusivity and viscosity "//&
2884	1	"in the surface boundary layer.", default=.false., do_not_log=.true.)
2885		endif
2886
2887	1	if (GV%Boussinesq) then
2888	1	Kv_units = "m2 s-1" ; Kd_units = "m2 s-1"
2889		else
2890	0	Kv_units = "Pa s" ; Kd_units = "kg m-1 s-1"
2891		endif
2892
2893	1	if (use_kappa_shear .or. useKPP .or. useEPBL .or. use_CVMix_shear .or. use_CVMix_conv) then
2894	1	call safe_alloc_ptr(visc%Kd_shear, isd, ied, jsd, jed, nz+1)
2895		call register_restart_field(visc%Kd_shear, "Kd_shear", .false., restart_CS, &
2896		"Shear-driven turbulent diffusivity at interfaces", &
2897	1	units=Kd_units, conversion=GV%HZ_T_to_MKS, z_grid='i')
2898		endif
2899	1	if (useKPP .or. useEPBL .or. use_CVMix_shear .or. use_CVMix_conv .or. &
2900		(use_kappa_shear .and. .not.KS_at_vertex )) then
2901	1	call safe_alloc_ptr(visc%Kv_shear, isd, ied, jsd, jed, nz+1)
2902		call register_restart_field(visc%Kv_shear, "Kv_shear", .false., restart_CS, &
2903		"Shear-driven turbulent viscosity at interfaces", &
2904	1	units=Kv_units, conversion=GV%HZ_T_to_MKS, z_grid='i')
2905		endif
2906	1	if (use_kappa_shear .and. KS_at_vertex) then
2907	1	call safe_alloc_ptr(visc%TKE_turb, HI%IsdB, HI%IedB, HI%JsdB, HI%JedB, nz+1)
2908	1	call safe_alloc_ptr(visc%Kv_shear_Bu, HI%IsdB, HI%IedB, HI%JsdB, HI%JedB, nz+1)
2909		call register_restart_field(visc%Kv_shear_Bu, "Kv_shear_Bu", .false., restart_CS, &
2910		"Shear-driven turbulent viscosity at vertex interfaces", &
2911	1	units=Kv_units, conversion=GV%HZ_T_to_MKS, hor_grid="Bu", z_grid='i')
2912	0	elseif (use_kappa_shear) then
2913	0	call safe_alloc_ptr(visc%TKE_turb, isd, ied, jsd, jed, nz+1)
2914		endif
2915
2916	1	if (useKPP) then
2917		! MOM_bkgnd_mixing uses Kv_slow when KPP is defined.
2918	0	call safe_alloc_ptr(visc%Kv_slow, isd, ied, jsd, jed, nz+1)
2919		endif
2920
2921		! visc%MLD and visc%h_ML are used to communicate the state of the (e)PBL or KPP to the rest of the model
2922		call get_param(param_file, mdl, "MLE_USE_PBL_MLD", MLE_use_PBL_MLD, &
2923	1	default=.false., do_not_log=.true.)
2924		! visc%h_ML needs to be allocated when melt potential is computed (HFREEZE>0) or one of
2925		! several other parameterizations are in use.
2926		call get_param(param_file, mdl, "HFREEZE", hfreeze, &
2927	1	units="m", default=-1.0, scale=US%m_to_Z, do_not_log=.true.)
2928		call get_param(param_file, mdl, "DO_BRINE_PLUME", do_brine_plume, &
2929		"If true, use a brine plume parameterization from Nguyen et al., 2009.", &
2930	1	default=.false., do_not_log=.true.)
2931		call get_param(param_file, mdl, "USE_HORIZONTAL_BOUNDARY_DIFFUSION", use_hor_bnd_diff, &
2932	1	default=.false., do_not_log=.true.)
2933		call get_param(param_file, mdl, "USE_NEUTRAL_DIFFUSION", use_neutral_diffusion, &
2934	1	default=.false., do_not_log=.true.)
2935	1	if (use_neutral_diffusion) &
2936		call get_param(param_file, mdl, "NDIFF_INTERIOR_ONLY", use_neutral_diffusion, &
2937	0	default=.false., do_not_log=.true.)
2938		call get_param(param_file, mdl, "FPMIX", use_fpmix, &
2939	1	default=.false., do_not_log=.true.)
2940		call get_param(param_file, mdl, "USE_IDEAL_AGE_TRACER", use_ideal_age, &
2941	1	default=.false., do_not_log=.true.)
2942	1	call openParameterBlock(param_file, 'MLE', do_not_log=.true.)
2943		call get_param(param_file, mdl, "USE_BODNER23", MLE_use_Bodner, &
2944	1	default=.false., do_not_log=.true.)
2945	1	call closeParameterBlock(param_file)
2946
2947	1	if (MLE_use_PBL_MLD .or. MLE_use_Bodner) then
2948	1	call safe_alloc_ptr(visc%MLD, isd, ied, jsd, jed)
2949		endif
2950		if ((hfreeze >= 0.0) .or. MLE_use_PBL_MLD .or. do_brine_plume .or. use_fpmix .or. &
2951	1	use_neutral_diffusion .or. use_hor_bnd_diff .or. use_ideal_age) then
2952	1	call safe_alloc_ptr(visc%h_ML, isd, ied, jsd, jed)
2953		endif
2954
2955	1	if (MLE_use_PBL_MLD) then
2956		call register_restart_field(visc%MLD, "MLD", .false., restart_CS, &
2957	1	"Instantaneous active mixing layer depth", units="m", conversion=US%Z_to_m)
2958		endif
2959		if (MLE_use_PBL_MLD .or. do_brine_plume .or. use_fpmix .or. &
2960	1	use_neutral_diffusion .or. use_hor_bnd_diff) then
2961		call register_restart_field(visc%h_ML, "h_ML", .false., restart_CS, &
2962		"Instantaneous active mixing layer thickness", &
2963	1	units=get_thickness_units(GV), conversion=GV%H_to_mks)
2964		endif
2965
2966		! visc%sfc_buoy_flx is used to communicate the state of the (e)PBL or KPP to the rest of the model
2967	1	if (MLE_use_PBL_MLD .or. MLE_use_Bodner) then
2968	1	call safe_alloc_ptr(visc%sfc_buoy_flx, isd, ied, jsd, jed)
2969		call register_restart_field(visc%sfc_buoy_flx, "SFC_BFLX", .false., restart_CS, &
2970		"Instantaneous surface buoyancy flux", "m2 s-3", &
2971	1	conversion=US%Z_to_m*2US%s_to_T**3)
2972		endif
2973
2974	1	if (use_ice_shelf) then
2975	0	if (.not.allocated(visc%taux_shelf)) &
2976	0	allocate(visc%taux_shelf(G%IsdB:G%IedB, G%jsd:G%jed), source=0.0)
2977	0	if (.not.allocated(visc%tauy_shelf)) &
2978	0	allocate(visc%tauy_shelf(G%isd:G%ied, G%JsdB:G%JedB), source=0.0)
2979		u_desc = var_desc("u_taux_shelf", "Pa", "the zonal stress on the ocean under ice shelves", &
2980	0	hor_grid='Cu',z_grid='1')
2981		v_desc = var_desc("v_tauy_shelf", "Pa", "the meridional stress on the ocean under ice shelves", &
2982	0	hor_grid='Cv',z_grid='1')
2983		call register_restart_pair(visc%taux_shelf, visc%tauy_shelf, u_desc, v_desc, &
2984	0	.false., restart_CS, conversion=US%RZ_T_to_kg_m2s*US%L_T_to_m_s)
2985		endif
2986
2987	11	end subroutine set_visc_register_restarts
2988
2989		!> This subroutine does remapping for the auxiliary restart variables in a vertvisc_type
2990		!! that are used across timesteps
2991	0	subroutine remap_vertvisc_aux_vars(G, GV, visc, h_old, h_new, ALE_CSp, OBC)
2992		type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
2993		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2994		type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical
2995		!! viscosities and related fields.
2996		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2997		intent(in) :: h_old !< Thickness of source grid [H ~> m or kg m-2]
2998		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2999		intent(in) :: h_new !< Thickness of destination grid [H ~> m or kg m-2]
3000		type(ALE_CS), pointer :: ALE_CSp !< ALE control structure to use when remapping
3001		type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
3002
3003	0	if (associated(visc%Kd_shear)) then
3004	0	call ALE_remap_interface_vals(ALE_CSp, G, GV, h_old, h_new, visc%Kd_shear)
3005		endif
3006
3007	0	if (associated(visc%Kv_shear)) then
3008	0	call ALE_remap_interface_vals(ALE_CSp, G, GV, h_old, h_new, visc%Kv_shear)
3009		endif
3010
3011	0	if (associated(visc%Kv_shear_Bu)) then
3012	0	call ALE_remap_vertex_vals(ALE_CSp, G, GV, h_old, h_new, visc%Kv_shear_Bu)
3013		endif
3014
3015	0	end subroutine remap_vertvisc_aux_vars
3016
3017		!> Initializes the MOM_set_visc control structure
3018	1	subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS, OBC)
3019		type(time_type), target, intent(in) :: Time !< The current model time.
3020		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
3021		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
3022		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3023		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
3024		!! parameters.
3025		type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
3026		!! output.
3027		type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
3028		!! related fields.
3029		type(set_visc_CS), intent(inout) :: CS !< Vertical viscosity control structure
3030		type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure
3031		type(ocean_OBC_type), pointer :: OBC !< A pointer to an open boundary condition structure
3032
3033		! Local variables
3034		real :: Csmag_chan_dflt ! The default value for SMAG_CONST_CHANNEL [nondim]
3035		real :: smag_const1 ! The default value for the Smagorinsky Laplacian coefficient [nondim]
3036		real :: TKE_decay_dflt ! The default value of a coefficient scaling the vertical decay
3037		! rate of TKE [nondim]
3038		real :: bulk_Ri_ML_dflt ! The default bulk Richardson number for a bulk mixed layer [nondim]
3039		real :: Kv_background ! The background kinematic viscosity in the interior [Z2 T-1 ~> m2 s-1]
3040		real :: omega_frac_dflt ! The default value for the fraction of the absolute rotation rate that
3041		! is used in place of the absolute value of the local Coriolis
3042		! parameter in the denominator of some expressions [nondim]
3043		real :: Chan_max_thick_dflt ! The default value for CHANNEL_DRAG_MAX_THICK [Z ~> m]
3044		real :: tideamp_factor ! A factor to multiply by tideamp when converting to mean tidal magnitude [nondim]
3045		real :: shelfbreak_depth ! When CHANNEL_DRAG is true, the bathymetric depth interpolated
3046		! to the vorticity point is a combination of the harmonic mean of the
3047		! adjacent velocity point depths below this depth [Z ~> m] and the
3048		! arithmetic mean of the adjacent depths above it, to roughly mimic a
3049		! continental shelf break profile.
3050	1	real, allocatable, dimension(:,:) :: cdrag_h !< The spatially varying quadratic drag coefficient [nondim]
3051
3052		integer :: i, j, is, ie, js, je
3053		integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
3054		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
3055		logical :: adiabatic, use_omega, MLE_use_PBL_MLD
3056		logical :: use_KPP
3057		logical :: use_regridding ! If true, use the ALE algorithm rather than layered
3058		! isopycnal or stacked shallow water mode.
3059		logical :: use_temperature ! If true, temperature and salinity are used as state variables.
3060		logical :: use_EOS ! If true, density calculated from T & S using an equation of state.
3061		character(len=200) :: filename, cdrag_file, tideamp_file ! Input file names or paths
3062		character(len=80) :: cdrag_var, tideamp_var ! Input file variable names
3063		! This include declares and sets the variable "version".
3064		# include "version_variable.h"
3065		character(len=40) :: mdl = "MOM_set_visc" ! This module's name.
3066
3067	1	CS%initialized = .true.
3068	1	CS%OBC => OBC
3069
3070	1	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
3071	1	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke
3072	1	IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
3073
3074	1	CS%diag => diag
3075
3076		! Set default, read and log parameters
3077	1	call log_version(param_file, mdl, version, "")
3078	1	CS%RiNo_mix = .false.
3079	1	call get_param(param_file, mdl, "INPUTDIR", CS%inputdir, default=".")
3080	1	CS%inputdir = slasher(CS%inputdir)
3081		call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
3082		"This sets the default value for the various _ANSWER_DATE parameters.", &
3083	1	default=99991231)
3084		call get_param(param_file, mdl, "SET_VISC_ANSWER_DATE", CS%answer_date, &
3085		"The vintage of the order of arithmetic and expressions in the set viscosity "//&
3086		"calculations. Values below 20190101 recover the answers from the end of 2018, "//&
3087		"while higher values use updated and more robust forms of the same expressions.", &
3088	1	default=default_answer_date, do_not_log=.not.GV%Boussinesq)
3089	1	if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701)
3090		call get_param(param_file, mdl, "BOTTOMDRAGLAW", CS%bottomdraglaw, &
3091		"If true, the bottom stress is calculated with a drag "//&
3092		"law of the form c_drag\|u\|u. The velocity magnitude "//&
3093		"may be an assumed value or it may be based on the "//&
3094		"actual velocity in the bottommost HBBL, depending on "//&
3095	1	"LINEAR_DRAG.", default=.true.)
3096		call get_param(param_file, mdl, "DRAG_AS_BODY_FORCE", CS%body_force_drag, &
3097		"If true, the bottom stress is imposed as an explicit body force "//&
3098		"applied over a fixed distance from the bottom, rather than as an "//&
3099		"implicit calculation based on an enhanced near-bottom viscosity. "//&
3100		"The thickness of the bottom boundary layer is HBBL.", &
3101	1	default=.false., do_not_log=.not.CS%bottomdraglaw)
3102		call get_param(param_file, mdl, "CHANNEL_DRAG", CS%Channel_drag, &
3103		"If true, the bottom drag is exerted directly on each "//&
3104		"layer proportional to the fraction of the bottom it overlies.", &
3105	1	default=.false.)
3106		call get_param(param_file, mdl, "CHANNEL_DRAG_SHELFBREAK_DEPTH", shelfbreak_depth, &
3107		"When CHANNEL_DRAG is true, the bathymetric depth interpolated to the "//&
3108		"vorticity point is a combination of the harmonic mean of the adjacent "//&
3109		"velocity point depths below this depth and the arithmetic mean of the "//&
3110		"depths above it, to roughly mimic a continental shelf break profile. "//&
3111		"Setting this to exceed MAXIMUM_DEPTH leads to linear interpolation of "//&
3112		"the topography between velocity points.", &
3113	1	default=0.0, units="m", scale=US%m_to_Z, do_not_log=.not.CS%Channel_drag)
3114	1	CS%channel_break_depth = shelfbreak_depth - G%Z_ref
3115
3116		call get_param(param_file, mdl, "LINEAR_DRAG", CS%linear_drag, &
3117		"If LINEAR_DRAG and BOTTOMDRAGLAW are defined the drag "//&
3118	1	"law is cdragDRAG_BG_VELu.", default=.false.)
3119		call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
3120	1	do_not_log=.true.)
3121	1	if (adiabatic) then
3122		call log_param(param_file, mdl, "ADIABATIC",adiabatic, &
3123		"There are no diapycnal mass fluxes if ADIABATIC is true. "//&
3124		"This assumes that KD = 0.0 and that there is no buoyancy forcing, "//&
3125	0	"but makes the model faster by eliminating subroutine calls.", default=.false.)
3126		endif
3127
3128	1	if (.not.adiabatic) then
3129	1	CS%RiNo_mix = kappa_shear_is_used(param_file)
3130		endif
3131
3132	1	call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.)
3133
3134		call get_param(param_file, mdl, "DYNAMIC_VISCOUS_ML", CS%dynamic_viscous_ML, &
3135		"If true, use a bulk Richardson number criterion to "//&
3136		"determine the mixed layer thickness for viscosity.", &
3137	1	default=.false.)
3138	1	if (CS%dynamic_viscous_ML) then
3139	0	call get_param(param_file, mdl, "BULK_RI_ML", bulk_Ri_ML_dflt, units="nondim", default=0.0)
3140		call get_param(param_file, mdl, "BULK_RI_ML_VISC", CS%bulk_Ri_ML, &
3141		"The efficiency with which mean kinetic energy released by mechanically "//&
3142		"forced entrainment of the mixed layer is converted to turbulent "//&
3143		"kinetic energy. By default, BULK_RI_ML_VISC = BULK_RI_ML or 0.", &
3144	0	units="nondim", default=bulk_Ri_ML_dflt)
3145	0	call get_param(param_file, mdl, "TKE_DECAY", TKE_decay_dflt, units="nondim", default=0.0)
3146		call get_param(param_file, mdl, "TKE_DECAY_VISC", CS%TKE_decay, &
3147		"TKE_DECAY_VISC relates the vertical rate of decay of "//&
3148		"the TKE available for mechanical entrainment to the "//&
3149		"natural Ekman depth for use in calculating the dynamic "//&
3150		"mixed layer viscosity. By default, TKE_DECAY_VISC = TKE_DECAY or 0.", &
3151	0	units="nondim", default=TKE_decay_dflt)
3152		call get_param(param_file, mdl, "ML_USE_OMEGA", use_omega, &
3153		"If true, use the absolute rotation rate instead of the "//&
3154		"vertical component of rotation when setting the decay "//&
3155	0	"scale for turbulence.", default=.false., do_not_log=.true.)
3156	0	omega_frac_dflt = 0.0
3157	0	if (use_omega) then
3158	0	call MOM_error(WARNING, "ML_USE_OMEGA is deprecated; use ML_OMEGA_FRAC=1.0 instead.")
3159	0	omega_frac_dflt = 1.0
3160		endif
3161		call get_param(param_file, mdl, "ML_OMEGA_FRAC", CS%omega_frac, &
3162		"When setting the decay scale for turbulence, use this "//&
3163		"fraction of the absolute rotation rate blended with the "//&
3164		"local value of f, as sqrt((1-of)f^2 + of4*omega^2).", &
3165	0	units="nondim", default=omega_frac_dflt)
3166		call get_param(param_file, mdl, "OMEGA", CS%omega, &
3167		"The rotation rate of the earth.", &
3168	0	units="s-1", default=7.2921e-5, scale=US%T_to_s)
3169		! This give a minimum decay scale that is typically much less than Angstrom.
3170	0	CS%ustar_min = 2e-4CS%omega(GV%Angstrom_H + GV%H_subroundoff)
3171		else
3172		call get_param(param_file, mdl, "OMEGA", CS%omega, &
3173		"The rotation rate of the earth.", &
3174	1	units="s-1", default=7.2921e-5, scale=US%T_to_s)
3175		endif
3176
3177		call get_param(param_file, mdl, "HBBL", CS%dz_bbl, &
3178		"The thickness of a bottom boundary layer with a viscosity increased by "//&
3179		"KV_EXTRA_BBL if BOTTOMDRAGLAW is not defined, or the thickness over which "//&
3180		"near-bottom velocities are averaged for the drag law if BOTTOMDRAGLAW is "//&
3181		"defined but LINEAR_DRAG is not.", &
3182	1	units="m", scale=US%m_to_Z, fail_if_missing=.true.) ! Rescaled later
3183	1	if (CS%bottomdraglaw) then
3184		call get_param(param_file, mdl, "CDRAG", CS%cdrag, &
3185		"CDRAG is the drag coefficient relating the magnitude of "//&
3186		"the velocity field to the bottom stress. CDRAG is only "//&
3187	1	"used if BOTTOMDRAGLAW is defined.", units="nondim", default=0.003)
3188		call get_param(param_file, mdl, "CDRAG_MAP", CS%bottomdragmap, &
3189		"If true, apply a spatially varying scaling factor to CDRAG, "//&
3190	1	"specified by CDRAG_VAR in CDRAG_FILE.", default=.false.)
3191		call get_param(param_file, mdl, "CDRAG_FILE", cdrag_file, &
3192		"The name of the file with the spatially varying bottom drag "//&
3193	1	"scaling factor.", default="", do_not_log=.not.CS%bottomdragmap)
3194		call get_param(param_file, mdl, "CDRAG_VAR", cdrag_var, &
3195		"The name of the variable in CDRAG_FILE with the spatially "//&
3196		"varying bottom drag scaling factor at h points.", &
3197	1	default="", do_not_log=.not.CS%bottomdragmap)
3198		call get_param(param_file, mdl, "BBL_USE_TIDAL_BG", CS%BBL_use_tidal_bg, &
3199		"Flag to use the tidal RMS amplitude in place of constant "//&
3200		"background velocity for computing u* in the BBL. "//&
3201		"This flag is only used when BOTTOMDRAGLAW is true and "//&
3202	1	"LINEAR_DRAG is false.", default=.false.)
3203	1	if (CS%BBL_use_tidal_bg) then
3204		call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, &
3205		"The path to the file containing the spatially varying "//&
3206	0	"tidal amplitudes with INT_TIDE_DISSIPATION.", default="tideamp.nc")
3207		call get_param(param_file, mdl, "TIDEAMP_VARNAME", tideamp_var, &
3208		"The name of the tidal amplitude variable in the input file.", &
3209	0	default="tideamp")
3210		! This value is here only to detect whether it is inadvertently used. CS%drag_bg_vel should
3211		! not be used if CS%BBL_use_tidal_bg is True. For this reason, we do not apply dimensions,
3212		! nor dimensional testing in this mode. If we ever detect a dimensional sensitivity to
3213		! this parameter, in this mode, then it means it is being used inappropriately.
3214	0	CS%drag_bg_vel = 1.e30
3215		call get_param(param_file, mdl, "TIDEAMP_FACTOR", tideamp_factor, &
3216		"A parameter to multiply by tideamp when converting to ustar. "//&
3217		"It accounts for converting the amplitude to a mean magintude (approx 1/sqrt(2)) "//&
3218		"and possibly also for non-commuting averaging operators when converting to ustar**3. "//&
3219		"It is ignored if negative and uncapped so it can be greater than 1 if desired.",&
3220	0	units="nondim", default=-1.0)
3221	0	if (tideamp_factor < 0.0) then
3222	0	CS%tideampfac2 = 1.0
3223		else
3224	0	CS%tideampfac2 = tideamp_factor*tideamp_factor
3225		endif
3226		else
3227		call get_param(param_file, mdl, "DRAG_BG_VEL", CS%drag_bg_vel, &
3228		"DRAG_BG_VEL is either the assumed bottom velocity (with "//&
3229		"LINEAR_DRAG) or an unresolved velocity that is "//&
3230		"combined with the resolved velocity to estimate the "//&
3231		"velocity magnitude. DRAG_BG_VEL is only used when "//&
3232	1	"BOTTOMDRAGLAW is defined.", units="m s-1", default=0.0, scale=US%m_s_to_L_T)
3233		endif
3234		call get_param(param_file, mdl, "USE_REGRIDDING", use_regridding, &
3235	1	do_not_log=.true., default=.false. )
3236		call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
3237	1	default=.true., do_not_log=.true.)
3238		call get_param(param_file, mdl, "USE_EOS", use_EOS, &
3239	1	default=use_temperature, do_not_log=.true.)
3240		call get_param(param_file, mdl, "BBL_USE_EOS", CS%BBL_use_EOS, &
3241		"If true, use the equation of state in determining the properties of the "//&
3242		"bottom boundary layer. Otherwise use the layer target potential densities. "//&
3243		"The default of this parameter is the value of USE_EOS.", &
3244	1	default=use_EOS, do_not_log=.not.use_temperature)
3245	1	if (use_regridding .and. (.not. CS%BBL_use_EOS)) &
3246	0	call MOM_error(FATAL,"When using MOM6 in ALE mode it is required to set BBL_USE_EOS to True.")
3247		endif
3248		call get_param(param_file, mdl, "BBL_THICK_MIN", CS%BBL_thick_min, &
3249		"The minimum bottom boundary layer thickness that can be "//&
3250		"used with BOTTOMDRAGLAW. This might be "//&
3251		"Kv/(cdrag*drag_bg_vel) to give Kv as the minimum "//&
3252	1	"near-bottom viscosity.", units="m", default=0.0, scale=US%m_to_Z)
3253		call get_param(param_file, mdl, "HTBL_SHELF_MIN", CS%Htbl_shelf_min, &
3254		"The minimum top boundary layer thickness that can be "//&
3255		"used with BOTTOMDRAGLAW. This might be "//&
3256		"Kv/(cdrag*drag_bg_vel) to give Kv as the minimum "//&
3257	1	"near-top viscosity.", units="m", default=US%Z_to_m*CS%BBL_thick_min, scale=US%m_to_Z)
3258		call get_param(param_file, mdl, "HTBL_SHELF", CS%Htbl_shelf, &
3259		"The thickness over which near-surface velocities are "//&
3260		"averaged for the drag law under an ice shelf. By "//&
3261		"default this is the same as HBBL", &
3262	1	units="m", default=US%Z_to_m*CS%dz_bbl, scale=GV%m_to_H)
3263
3264		call get_param(param_file, mdl, "KV", Kv_background, &
3265		"The background kinematic viscosity in the interior. "//&
3266		"The molecular value, ~1e-6 m2 s-1, may be used.", &
3267	1	units="m2 s-1", scale=US%m2_s_to_Z2_T, fail_if_missing=.true.)
3268
3269		call get_param(param_file, mdl, "USE_KPP", use_KPP, &
3270		"If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
3271		"to calculate diffusivities and non-local transport in the OBL.", &
3272	1	do_not_log=.true., default=.false.)
3273
3274		call get_param(param_file, mdl, "KV_BBL_MIN", CS%KV_BBL_min, &
3275		"The minimum viscosities in the bottom boundary layer.", &
3276	1	units="m2 s-1", default=US%Z2_T_to_m2_s*Kv_background, scale=GV%m2_s_to_HZ_T)
3277		call get_param(param_file, mdl, "KV_TBL_MIN", CS%KV_TBL_min, &
3278		"The minimum viscosities in the top boundary layer.", &
3279	1	units="m2 s-1", default=US%Z2_T_to_m2_s*Kv_background, scale=GV%m2_s_to_HZ_T)
3280		call get_param(param_file, mdl, "CORRECT_BBL_BOUNDS", CS%correct_BBL_bounds, &
3281		"If true, uses the correct bounds on the BBL thickness and "//&
3282		"viscosity so that the bottom layer feels the intended drag.", &
3283	1	default=.false.)
3284
3285	1	if (CS%Channel_drag) then
3286	1	call get_param(param_file, mdl, "SMAG_LAP_CONST", smag_const1, units="nondim", default=-1.0)
3287
3288	1	cSmag_chan_dflt = 0.15
3289	1	if (smag_const1 >= 0.0) cSmag_chan_dflt = smag_const1
3290
3291		call get_param(param_file, mdl, "SMAG_CONST_CHANNEL", CS%c_Smag, &
3292		"The nondimensional Laplacian Smagorinsky constant used "//&
3293		"in calculating the channel drag if it is enabled. The "//&
3294		"default is to use the same value as SMAG_LAP_CONST if "//&
3295		"it is defined, or 0.15 if it is not. The value used is "//&
3296		"also 0.15 if the specified value is negative.", &
3297	1	units="nondim", default=cSmag_chan_dflt, do_not_log=.not.CS%Channel_drag)
3298	1	if (CS%c_Smag < 0.0) CS%c_Smag = 0.15
3299
3300		call get_param(param_file, mdl, "TRIG_CHANNEL_DRAG_WIDTHS", CS%concave_trigonometric_L, &
3301		"If true, use trigonometric expressions to determine the fractional open "//&
3302		"interface lengths for concave topography.", &
3303	1	default=.true., do_not_log=.not.CS%Channel_drag)
3304		endif
3305
3306	1	Chan_max_thick_dflt = -1.0*US%m_to_Z
3307	1	if (CS%RiNo_mix) Chan_max_thick_dflt = 0.5*CS%dz_bbl
3308	1	if (CS%body_force_drag) Chan_max_thick_dflt = CS%dz_bbl
3309		call get_param(param_file, mdl, "CHANNEL_DRAG_MAX_BBL_THICK", CS%Chan_drag_max_vol, &
3310		"The maximum bottom boundary layer thickness over which the channel drag is "//&
3311		"exerted, or a negative value for no fixed limit, instead basing the BBL "//&
3312		"thickness on the bottom stress, rotation and stratification. The default is "//&
3313		"proportional to HBBL if USE_JACKSON_PARAM or DRAG_AS_BODY_FORCE is true.", &
3314		units="m", default=US%Z_to_m*Chan_max_thick_dflt, scale=US%m_to_Z, &
3315	1	do_not_log=.not.CS%Channel_drag)
3316
3317		call get_param(param_file, mdl, "MLE_USE_PBL_MLD", MLE_use_PBL_MLD, &
3318	1	default=.false., do_not_log=.true.)
3319
3320	1	CS%Hbbl = CS%dz_bbl * (US%Z_to_m * GV%m_to_H) ! Rescaled for use in expressions in thickness units.
3321
3322		! Update scalar data to device
3323		!$omp target update to(visc)
3324		!$omp target update to(CS)
3325
3326		! Now update arrays that were defined in set_visc_register_restarts()
3327		!$omp target enter data map(to: visc%Kv_shear) &
3328		!$omp if (associated(visc%Kv_shear))
3329		!$omp target enter data map(to: visc%Kv_shear_Bu) &
3330		!$omp if (associated(visc%Kv_shear_Bu))
3331
3332	1	if (CS%RiNo_mix .and. kappa_shear_at_vertex(param_file)) then
3333		! This is necessary for reproducibility across restarts in non-symmetric mode.
3334	1	call pass_var(visc%Kv_shear_Bu, G%Domain, position=CORNER, complete=.true.)
3335		endif
3336
3337	1	if (CS%bottomdraglaw) then
3338	8841	allocate(visc%bbl_thick_u(IsdB:IedB,jsd:jed), source=0.0)
3339	8902	allocate(visc%bbl_thick_v(isd:ied,JsdB:JedB), source=0.0)
3340		!$omp target enter data map(to: visc%bbl_thick_u, visc%bbl_thick_v)
3341	8841	allocate(visc%kv_bbl_u(IsdB:IedB,jsd:jed), source=0.0)
3342	8902	allocate(visc%kv_bbl_v(isd:ied,JsdB:JedB), source=0.0)
3343		!$omp target enter data map(to: visc%kv_bbl_u, visc%kv_bbl_v)
3344	8773	allocate(visc%ustar_bbl(isd:ied,jsd:jed), source=0.0)
3345	8773	allocate(visc%BBL_meanKE_loss(isd:ied,jsd:jed), source=0.0)
3346	8773	allocate(visc%BBL_meanKE_loss_sqrtCd(isd:ied,jsd:jed), source=0.0)
3347
3348		CS%id_bbl_thick_u = register_diag_field('ocean_model', 'bbl_thick_u', &
3349	1	diag%axesCu1, Time, 'BBL thickness at u points', 'm', conversion=US%Z_to_m)
3350		CS%id_kv_bbl_u = register_diag_field('ocean_model', 'kv_bbl_u', diag%axesCu1, &
3351	1	Time, 'BBL viscosity at u points', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
3352		CS%id_bbl_u = register_diag_field('ocean_model', 'bbl_u', diag%axesCu1, &
3353	1	Time, 'BBL mean u current', 'm s-1', conversion=US%L_T_to_m_s)
3354	1	if (CS%id_bbl_u>0) then
3355	0	allocate(CS%bbl_u(IsdB:IedB,jsd:jed), source=0.0)
3356		endif
3357		CS%id_bbl_thick_v = register_diag_field('ocean_model', 'bbl_thick_v', &
3358	1	diag%axesCv1, Time, 'BBL thickness at v points', 'm', conversion=US%Z_to_m)
3359		CS%id_kv_bbl_v = register_diag_field('ocean_model', 'kv_bbl_v', diag%axesCv1, &
3360	1	Time, 'BBL viscosity at v points', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
3361		CS%id_bbl_v = register_diag_field('ocean_model', 'bbl_v', diag%axesCv1, &
3362	1	Time, 'BBL mean v current', 'm s-1', conversion=US%L_T_to_m_s)
3363	1	if (CS%id_bbl_v>0) then
3364	0	allocate(CS%bbl_v(isd:ied,JsdB:JedB), source=0.0)
3365		endif
3366	1	if (CS%bottomdragmap) then
3367	0	if (len_trim(cdrag_file)==0 .or. len_trim(cdrag_var)==0) then
3368	0	call MOM_error(FATAL,"CDRAG_FILE and CDRAG_VAR are required when using CDRAG_MAP.")
3369		endif
3370	0	allocate(cdrag_h(isd:ied,jsd:jed), source=0.0)
3371	0	allocate(CS%cdrag_u(IsdB:IedB,jsd:jed), source=0.0)
3372	0	allocate(CS%cdrag_v(isd:ied,JsdB:JedB), source=0.0)
3373	0	filename = trim(CS%inputdir) // trim(cdrag_file)
3374	0	call log_param(param_file, mdl, "INPUTDIR/CDRAG_FILE", filename)
3375	0	call MOM_read_data(filename, cdrag_var, cdrag_h, G%domain, scale=CS%cdrag)
3376	0	call pass_var(cdrag_h, G%domain)
3377	0	do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0) then
3378		CS%cdrag_u(I,j) = (G%mask2dT(i,j) * cdrag_h(i,j) + G%mask2dT(i+1,j) * cdrag_h(i+1,j)) / &
3379	0	(G%mask2dT(i,j) + G%mask2dT(i+1,j))
3380		endif ; enddo ; enddo
3381	0	do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0) then
3382		CS%cdrag_v(i,J) = (G%mask2dT(i,j) * cdrag_h(i,j) + G%mask2dT(i,j+1) * cdrag_h(i,j+1)) / &
3383	0	(G%mask2dT(i,j) + G%mask2dT(i,j+1))
3384		endif ; enddo ; enddo
3385	0	deallocate(cdrag_h)
3386		endif
3387	1	if (CS%BBL_use_tidal_bg) then
3388	0	allocate(CS%tideamp(isd:ied,jsd:jed), source=0.0)
3389	0	filename = trim(CS%inputdir) // trim(tideamp_file)
3390	0	call log_param(param_file, mdl, "INPUTDIR/TIDEAMP_FILE", filename)
3391	0	call MOM_read_data(filename, tideamp_var, CS%tideamp, G%domain, scale=US%m_to_Z*US%T_to_s)
3392	0	call pass_var(CS%tideamp,G%domain)
3393		endif
3394		endif
3395	1	if (CS%Channel_drag .or. CS%body_force_drag) then
3396	663076	allocate(visc%Ray_u(IsdB:IedB,jsd:jed,nz), source=0.0)
3397	667651	allocate(visc%Ray_v(isd:ied,JsdB:JedB,nz), source=0.0)
3398		!$omp target enter data map(to: visc%Ray_u, visc%Ray_v)
3399		CS%id_Ray_u = register_diag_field('ocean_model', 'Rayleigh_u', diag%axesCuL, &
3400	1	Time, 'Rayleigh drag velocity at u points', 'm s-1', conversion=GV%H_to_m*US%s_to_T)
3401		CS%id_Ray_v = register_diag_field('ocean_model', 'Rayleigh_v', diag%axesCvL, &
3402	1	Time, 'Rayleigh drag velocity at v points', 'm s-1', conversion=GV%H_to_m*US%s_to_T)
3403		endif
3404
3405
3406	1	if (CS%dynamic_viscous_ML) then
3407	0	allocate(visc%nkml_visc_u(IsdB:IedB,jsd:jed), source=0.0)
3408	0	allocate(visc%nkml_visc_v(isd:ied,JsdB:JedB), source=0.0)
3409		!$omp target enter data map(to: visc%nkml_visc_u, visc%nkml_visc_v)
3410
3411		CS%id_nkml_visc_u = register_diag_field('ocean_model', 'nkml_visc_u', &
3412	0	diag%axesCu1, Time, 'Number of layers in viscous mixed layer at u points', 'nondim')
3413		CS%id_nkml_visc_v = register_diag_field('ocean_model', 'nkml_visc_v', &
3414	0	diag%axesCv1, Time, 'Number of layers in viscous mixed layer at v points', 'nondim')
3415		endif
3416
3417	1	call register_restart_field_as_obsolete('Kd_turb','Kd_shear', restart_CS)
3418	1	call register_restart_field_as_obsolete('Kv_turb','Kv_shear', restart_CS)
3419
3420	1	end subroutine set_visc_init
3421
3422		!> This subroutine dellocates any memory in the set_visc control structure.
3423	1	subroutine set_visc_end(visc, CS)
3424		type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
3425		!! related fields. Elements are deallocated here.
3426		type(set_visc_CS), intent(inout) :: CS !< The control structure returned by a previous
3427		!! call to set_visc_init.
3428
3429	1	if (allocated(visc%bbl_thick_u)) then
3430		!$omp target exit data map(delete: visc%bbl_thick_u)
3431	1	deallocate(visc%bbl_thick_u)
3432		endif
3433	1	if (allocated(visc%bbl_thick_v)) then
3434		!$omp target exit data map(delete: visc%bbl_thick_v)
3435	1	deallocate(visc%bbl_thick_v)
3436		endif
3437	1	if (allocated(visc%kv_bbl_u)) then
3438		!$omp target exit data map(delete: visc%kv_bbl_u)
3439	1	deallocate(visc%kv_bbl_u)
3440		endif
3441	1	if (allocated(visc%kv_bbl_v)) then
3442		!$omp target exit data map(delete: visc%kv_bbl_v)
3443	1	deallocate(visc%kv_bbl_v)
3444		endif
3445	1	if (allocated(CS%bbl_u)) deallocate(CS%bbl_u)
3446	1	if (allocated(CS%bbl_v)) deallocate(CS%bbl_v)
3447	1	if (allocated(visc%Ray_u)) then
3448		!$omp target exit data map(delete: visc%Ray_u)
3449	1	deallocate(visc%Ray_u)
3450		endif
3451	1	if (allocated(visc%Ray_v)) then
3452		!$omp target exit data map(delete: visc%Ray_v)
3453	1	deallocate(visc%Ray_v)
3454		endif
3455	1	if (allocated(visc%nkml_visc_u)) then
3456		!$omp target exit data map(delete: visc%nkml_visc_u)
3457	0	deallocate(visc%nkml_visc_u)
3458		endif
3459	1	if (allocated(visc%nkml_visc_v)) then
3460		!$omp target exit data map(delete: visc%nkml_visc_v)
3461	0	deallocate(visc%nkml_visc_v)
3462		endif
3463	1	if (associated(visc%Kd_shear)) deallocate(visc%Kd_shear)
3464	1	if (associated(visc%Kv_slow)) deallocate(visc%Kv_slow)
3465	1	if (associated(visc%TKE_turb)) deallocate(visc%TKE_turb)
3466	1	if (associated(visc%Kv_shear)) deallocate(visc%Kv_shear)
3467	1	if (associated(visc%Kv_shear_Bu)) deallocate(visc%Kv_shear_Bu)
3468	1	if (allocated(visc%ustar_bbl)) deallocate(visc%ustar_bbl)
3469	1	if (allocated(visc%BBL_meanKE_loss)) deallocate(visc%BBL_meanKE_loss)
3470	1	if (allocated(visc%BBL_meanKE_loss_sqrtCd)) deallocate(visc%BBL_meanKE_loss_sqrtCd)
3471	1	if (allocated(visc%taux_shelf)) deallocate(visc%taux_shelf)
3472	1	if (allocated(visc%tauy_shelf)) deallocate(visc%tauy_shelf)
3473	1	if (allocated(visc%tbl_thick_shelf_u)) deallocate(visc%tbl_thick_shelf_u)
3474	1	if (allocated(visc%tbl_thick_shelf_v)) deallocate(visc%tbl_thick_shelf_v)
3475	1	if (allocated(visc%kv_tbl_shelf_u)) deallocate(visc%kv_tbl_shelf_u)
3476	1	if (allocated(visc%kv_tbl_shelf_v)) deallocate(visc%kv_tbl_shelf_v)
3477	1	end subroutine set_visc_end
3478
3479		!> \namespace mom_set_visc
3480		!!
3481		!! This would also be the module in which other viscous quantities that are flow-independent might be set.
3482		!! This information is transmitted to other modules via a vertvisc type structure.
3483		!!
3484		!! The same code is used for the two velocity components, by indirectly referencing the velocities and
3485		!! defining a handful of direction-specific defined variables.
3486
3487	0	end module MOM_set_visc