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		!$omp declare target
1243		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1244		real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by
1245		!! the full horizontal area of a velocity cell [Z ~> m]
1246		real, intent(in) :: Dp !< The larger of the two depths at the edge
1247		!! of a velocity cell [Z ~> m]
1248		real, intent(in) :: Dm !< The smaller of the two depths at the edge
1249		!! of a velocity cell [Z ~> m]
1250		real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at
1251		!! the depth of each interface [nondim]
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
1259	768	nz = GV%ke
1260
1261	768	slope = abs(Dp - Dm)
1262	768	if (slope == 0.0) then
1263	0	L(1:nz) = 1.0 ; L(nz+1) = 0.0
1264		else
1265	768	Vol_open = 0.5*slope
1266	768	I_slope = 1.0 / slope
1267
1268	768	L(nz+1) = 0.0
1269	58368	do K=nz,1,-1
1270	58368	if (vol_below(K) >= Vol_open) then ; L(K) = 1.0
1271		else
1272		! With a uniformly sloping bottom, the calculation of L(K) is the solution of a simple quadratic equation.
1273	15312	L(K) = sqrt(2.0vol_below(K)I_slope)
1274		endif
1275		enddo
1276		endif
1277
1278	768	end subroutine find_L_open_uniform_slope
1279
1280		!> Determine the normalized open length of each interface for concave bathymetry (from the ocean perspective)
1281		!! using trigonometric expressions. In this case there can be two separate open regions.
1282	0	pure subroutine find_L_open_concave_trigonometric(vol_below, D_vel, Dp, Dm, L, GV)
1283		!$omp declare target
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
1295		! Local variables
1296		real :: crv ! crv is the curvature of the bottom depth across a
1297		! cell, times the cell width squared [Z ~> m].
1298		real :: crv_3 ! crv/3 [Z ~> m].
1299		real :: slope ! The absolute value of the bottom depth slope across
1300		! a cell times the cell width [Z ~> m].
1301		! The following "volumes" have units of vertical heights because they are normalized
1302		! by the full horizontal area of a velocity cell.
1303		real :: Vol_open ! The cell volume above which the face is fully is open [Z ~> m].
1304		real :: Vol_2_reg ! The cell volume above which there are two separate
1305		! open areas that must be integrated [Z ~> m].
1306		real :: C24_crv ! 24/crv [Z-1 ~> m-1].
1307		real :: apb_4a, ax2_3apb ! Various nondimensional ratios of crv and slope [nondim].
1308		real :: a2x48_apb3, Iapb ! Combinations of crv (a) and slope (b) [Z-1 ~> m-1]
1309		real :: L0 ! A linear estimate of L appropriate for tiny volumes [nondim].
1310		real :: slope_crv ! The slope divided by the curvature [nondim]
1311		real :: tmp_val_m1_to_p1 ! A temporary variable [nondim]
1312		real, parameter :: C1_3 = 1.0/3.0, C1_12 = 1.0/12.0 ! Rational constants [nondim]
1313		real, parameter :: C2pi_3 = 8.0*atan(1.0)/3.0 ! An irrational constant, 2/3 pi. [nondim]
1314		integer :: K, nz
1315
1316	0	nz = GV%ke
1317
1318		! Each cell extends from x=-1/2 to 1/2, and has a topography
1319		! given by D(x) = crvx^2 + slopex + D_vel - crv/12.
1320		!crv_3 = (Dp + Dm - 2.0D_vel) ; crv = 3.0crv_3
1321	0	crv_3 = (Dp + Dm - (2.0D_vel)) ; crv = 3.0crv_3
1322	0	slope = Dp - Dm
1323
1324		! Calculate the volume above which the entire cell is open and the volume at which the
1325		! equation that is solved for L changes because there are two separate open regions.
1326	0	if (slope >= crv) then
1327	0	Vol_open = D_vel - Dm ; Vol_2_reg = Vol_open
1328		else
1329	0	slope_crv = slope / crv
1330	0	Vol_open = 0.25slopeslope_crv + C1_12*crv
1331	0	Vol_2_reg = 0.5slope_crv2 (crv - C1_3*slope)
1332		endif
1333		! Define some combinations of crv & slope for later use.
1334	0	C24_crv = 24.0/crv ; Iapb = 1.0/(crv+slope)
1335	0	apb_4a = (slope+crv)/(4.0crv) ; a2x48_apb3 = (48.0(crvcrv))(Iapb**3)
1336	0	ax2_3apb = 2.0C1_3crv*Iapb
1337
1338	0	L(nz+1) = 0.0
1339		! Determine the normalized open length (L) at each interface.
1340	0	do K=nz,1,-1
1341	0	if (vol_below(K) >= Vol_open) then ! The whole cell is open.
1342	0	L(K) = 1.0
1343	0	elseif (vol_below(K) < Vol_2_reg) then
1344		! In this case, there is a contiguous open region and
1345		! vol_below(K) = 0.5L^2(slope + crv/3*(3-4L)).
1346	0	if (a2x48_apb3*vol_below(K) < 1e-8) then ! Could be 1e-7?
1347		! There is a very good approximation here for massless layers.
1348		!L0 = sqrt(2.0vol_below(K)Iapb) ; L(K) = L0(1.0 + ax2_3apbL0)
1349	0	L0 = sqrt(2.0vol_below(K)Iapb) ; L(K) = L0(1.0 + (ax2_3apbL0))
1350		else
1351		!L(K) = apb_4a * (1.0 - &
1352		! 2.0 * cos(C1_3acos(a2x48_apb3vol_below(K) - 1.0) - C2pi_3))
1353		L(K) = apb_4a * (1.0 - &
1354	0	2.0 * cos(C1_3acos((a2x48_apb3vol_below(K)) - 1.0) - C2pi_3))
1355		endif
1356		! To check the answers.
1357		! Vol_err = 0.5(L(K)L(K))(slope + crv_3(3.0-4.0*L(K))) - vol_below(K)
1358		else ! There are two separate open regions.
1359		! vol_below(K) = slope^2/4crv + crv/12 - (crv/12)(1-L)^2(1+2L)
1360		! At the deepest volume, L = slope/crv, at the top L = 1.
1361		! L(K) = 0.5 - cos(C1_3acos(1.0 - C24_crv(Vol_open - vol_below(K))) - C2pi_3)
1362	0	tmp_val_m1_to_p1 = 1.0 - C24_crv*(Vol_open - vol_below(K))
1363	0	tmp_val_m1_to_p1 = max(-1., min(1., tmp_val_m1_to_p1))
1364	0	L(K) = 0.5 - cos(C1_3*acos(tmp_val_m1_to_p1) - C2pi_3)
1365		! To check the answers.
1366		! Vol_err = Vol_open - 0.25crv_3(1.0+2.0L(K)) (1.0-L(K))**2 - vol_below(K)
1367		endif
1368		enddo ! k loop to determine L(K) in the concave case
1369
1370	0	end subroutine find_L_open_concave_trigonometric
1371
1372
1373
1374		!> Determine the normalized open length of each interface for concave bathymetry (from the ocean perspective) using
1375		!! iterative methods to solve the relevant cubic equations. In this case there can be two separate open regions.
1376	11664	pure subroutine find_L_open_concave_iterative(vol_below, D_vel, Dp, Dm, L, GV)
1377		!$omp declare target
1378		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1379		real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by
1380		!! the full horizontal area of a velocity cell [Z ~> m]
1381		real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m]
1382		real, intent(in) :: Dp !< The larger of the two depths at the edge
1383		!! of a velocity cell [Z ~> m]
1384		real, intent(in) :: Dm !< The smaller of the two depths at the edge
1385		!! of a velocity cell [Z ~> m]
1386		real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at
1387		!! the depth of each interface [nondim]
1388
1389		! Local variables
1390		real :: crv ! crv is the curvature of the bottom depth across a
1391		! cell, times the cell width squared [Z ~> m].
1392		real :: crv_3 ! crv/3 [Z ~> m].
1393		real :: slope ! The absolute value of the bottom depth slope across
1394		! a cell times the cell width [Z ~> m].
1395
1396		! The following "volumes" have units of vertical heights because they are normalized
1397		! by the full horizontal area of a velocity cell.
1398		real :: Vol_open ! The cell volume above which the face is fully is open [Z ~> m].
1399		real :: Vol_2_reg ! The cell volume above which there are two separate
1400		! open areas that must be integrated [Z ~> m].
1401		real :: L_2_reg ! The value of L when vol_below is Vol_2_reg [nondim]
1402		real :: vol_inflect_1 ! The volume at which there is an inflection point in the expression
1403		! relating L to vol_err when there is a single open region [Z ~> m]
1404		real :: vol_inflect_2 ! The volume at which there is an inflection point in the expression
1405		! relating L to vol_err when there are two open regions [Z ~> m]
1406
1407		real :: L_inflect_1 ! The value of L that sits at an inflection point in the expression
1408		! relating L to vol_err when there is a single open region [nondim]
1409		real :: L_inflect_2 ! The value of L that sits at an inflection point in the expression
1410		! relating L to vol_err when there is are two open regions [nondim]
1411		real :: L_max, L_min ! Maximum and minimum bounds on the solution for L for an interface [nondim]
1412		real :: vol_err ! The difference between the volume below an interface for a given value
1413		! of L and the target value [Z ~> m]
1414		real :: dVol_dL ! The partial derivative of the volume below with L [Z ~> m]
1415		real :: vol_err_max ! The value of vol_err when L is L_max [Z ~> m]
1416
1417		! The following combinations of slope and crv are reused across layers, and hence are pre-calculated
1418		! for efficiency. All are non-negative.
1419		real :: Icrvpslope ! The inverse of the sum of crv and slope [Z-1 ~> m-1]
1420		real :: slope_crv ! The slope divided by the curvature [nondim]
1421		! These are only used if the slope exceeds or matches the curvature.
1422		real :: smc ! The slope minus the curvature [Z ~> m]
1423		real :: C3c_m_s ! 3 times the curvature minus the slope [Z ~> m]
1424		real :: I_3c_m_s ! The inverse of 3 times the curvature minus the slope [Z-1 ~> m-1]
1425		! These are only used if the curvature exceeds the slope.
1426		real :: C4_crv ! The inverse of a quarter of the curvature [Z-1 ~> m-1]
1427		real :: sxcms_c ! The slope times the difference between the curvature and slope
1428		! divided by the curvature [Z ~> m]
1429		real :: slope2_4crv ! A quarter of the slope squared divided by the curvature [Z ~> m]
1430		real :: I_3s_m_c ! The inverse of 3 times the slope minus the curvature [Z-1 ~> m-1]
1431		real :: C3s_m_c ! 3 times the slope minus the curvature [Z ~> m]
1432
1433		real, parameter :: C1_3 = 1.0 / 3.0, C1_12 = 1.0 / 12.0 ! Rational constants [nondim]
1434		integer :: K, nz, itt
1435		integer, parameter :: max_itt = 10
1436
1437	11664	nz = GV%ke
1438
1439		! Each cell extends from x=-1/2 to 1/2, and has a topography
1440		! given by D(x) = crvx^2 + slopex + D_vel - crv/12.
1441
1442	11664	crv_3 = (Dp + Dm - 2.0D_vel) ; crv = 3.0crv_3
1443	11664	slope = Dp - Dm
1444
1445		! Calculate the volume above which the entire cell is open and the volume at which the
1446		! equation that is solved for L changes because there are two separate open regions.
1447	11664	if (slope >= crv) then
1448	9024	Vol_open = D_vel - Dm ; Vol_2_reg = Vol_open
1449	9024	L_2_reg = 1.0
1450	9024	if (crv + slope >= 4.0*crv) then
1451	5544	L_inflect_1 = 1.0 ; Vol_inflect_1 = Vol_open
1452		else
1453	3480	slope_crv = slope / crv
1454	3480	L_inflect_1 = 0.25 + 0.25*slope_crv
1455	3480	vol_inflect_1 = 0.25C1_12 ((slope_crv + 1.0)*2 (slope + crv))
1456		endif
1457		! Precalculate some combinations of crv & slope for later use.
1458	9024	smc = slope - crv
1459	9024	C3c_m_s = 3.0*crv - slope
1460	9024	if (C3c_m_s > 2.0*smc) I_3c_m_s = 1.0 / C3c_m_s
1461		else
1462	2640	slope_crv = slope / crv
1463	2640	Vol_open = 0.25slopeslope_crv + C1_12*crv
1464	2640	Vol_2_reg = 0.5slope_crv2 (crv - C1_3*slope)
1465	2640	L_2_reg = slope_crv
1466
1467		! The inflection point is useful to know because below the inflection point
1468		! Newton's method converges monotonically from above and conversely above it.
1469		! These are the inflection point values of L and vol_below with a single open segment.
1470	2640	vol_inflect_1 = 0.25C1_12 ((slope_crv + 1.0)*2 (slope + crv))
1471	2640	L_inflect_1 = 0.25 + 0.25*slope_crv
1472		! These are the inflection point values of L and vol_below when there are two open segments.
1473		! Vol_inflect_2 = Vol_open - 0.125 * crv_3, which is equivalent to:
1474	2640	vol_inflect_2 = 0.25slopeslope_crv + 0.125*crv_3
1475	2640	L_inflect_2 = 0.5
1476		! Precalculate some combinations of crv & slope for later use.
1477	2640	C4_crv = 4.0 / crv
1478	2640	slope2_4crv = 0.25 * slope * slope_crv
1479	2640	sxcms_c = slope_crv*(crv - slope)
1480	2640	C3s_m_c = 3.0*slope - crv
1481	2640	if (C3s_m_c > 2.0*sxcms_c) I_3s_m_c = 1.0 / C3s_m_c
1482		endif
1483		! Define some combinations of crv & slope for later use.
1484	11664	Icrvpslope = 1.0 / (crv+slope)
1485
1486	11664	L(nz+1) = 0.0
1487		! Determine the normalized open length (L) at each interface.
1488	886464	do K=nz,1,-1
1489	886464	if (vol_below(K) >= Vol_open) then ! The whole cell is open.
1490	553728	L(K) = 1.0
1491	321072	elseif (vol_below(K) < Vol_2_reg) then
1492		! In this case, there is a single contiguous open region from x=1/2-L to 1/2.
1493		! Changing the horizontal variable in the expression from D(x) to D(L) gives:
1494		! x(L) = 1/2 - L
1495		! D(L) = crv(0.5 - L)^2 + slope(0.5 - L) + D_vel - crv/12
1496		! D(L) = crvL^2 - crvL + crv/4 + slope*(1/2 - L) + D_vel - crv/12
1497		! D(L) = crvL^2 - (slope+crv)L + slope/2 + D_vel + crv/6
1498		! D(0) = slope/2 + D_vel + crv/6 = (Dp - Dm)/2 + D_vel + (Dp + Dm - 2*D_vel)/2 = Dp
1499		! D(1) = crv - slope - crv + slope/2 + Dvel + crv/6 = D_vel - slope/2 + crv/6 = Dm
1500		!
1501		! vol_below = integral(y = 0 to L) D(y) dy - L * D(L)
1502		! = crv/3L^3 - (slope+crv)/2L^2 + (slope/2 + D_vel + crv/6)*L -
1503		! (crvL^2 - (slope+crv)L + slope/2 + D_vel + crv/6) * L
1504		! = -2/3 * crv * L^3 + 1/2 * (slope+crv) * L^2
1505		! vol_below(K) = 0.5L(K)2(slope + crv_3(3-4L(K)))
1506		! L(K) is between L(K+1) and slope_crv.
1507	319080	L_max = min(L_2_reg, 1.0)
1508	319080	if (vol_below(K) <= vol_inflect_1) L_max = min(L_max, L_inflect_1)
1509
1510	319080	L_min = L(K+1)
1511	319080	if (vol_below(K) >= vol_inflect_1) L_min = max(L_min, L_inflect_1)
1512
1513		! Ignoring the cubic term gives an under-estimate but is very accurate for near bottom
1514		! layers, so use this as a potential floor.
1515	319080	if (2.0vol_below(K)Icrvpslope > L_min*2) L_min = sqrt(2.0vol_below(K)*Icrvpslope)
1516
1517		! Start with L_min in most cases.
1518	319080	L(k) = L_min
1519
1520	319080	if (vol_below(K) <= vol_inflect_1) then
1521		! Starting with L_min below L_inflect_1, only the first overshooting iteration of Newton's
1522		! method needs bounding.
1523	318384	L(k) = L_min
1524	318384	vol_err = 0.5L(K)2 (slope + crv(1.0 - 4.0C1_3*L(K))) - vol_below(K)
1525		! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L_min is already the best solution.
1526	318384	if (vol_err < 0.0) then
1527	318384	dVol_dL = L(K) * (slope + crv(1.0 - 2.0L(k)))
1528	318384	if (L(K)dVol_dL > vol_err + L_maxdVol_dL) then
1529	24	L(K) = L_max
1530		else
1531	318360	L(K) = L(K) - (vol_err / dVol_dL)
1532		endif
1533
1534		! Subsequent iterations of Newton's method do not need bounds.
1535	937427	do itt=1,max_itt
1536	937427	vol_err = 0.5L(K)2 (slope + crv(1.0 - 4.0C1_3*L(K))) - vol_below(K)
1537	937427	dVol_dL = L(K) * (slope + crv(1.0 - 2.0L(k)))
1538	937427	if (abs(vol_err) < max(1.0e-15L(K), 1.0e-25)dVol_dL) exit
1539	619043	L(K) = L(K) - (vol_err / dVol_dL)
1540		enddo
1541		endif
1542		else ! (vol_below(K) > vol_inflect_1)
1543		! Iteration from below converges monotonically, but we need to deal with the case where we are
1544		! close to the peak of the topography and Newton's method mimics the convergence of bisection.
1545
1546		! Evaluate the error when L(K) = L_min as a possible first guess.
1547	696	L(k) = L_min
1548	696	vol_err = 0.5L(K)2 (slope + crv(1.0 - 4.0C1_3*L(K))) - vol_below(K)
1549		! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L_min is already the best solution.
1550	696	if (vol_err < 0.0) then
1551
1552		! These two upper estimates deal with the possibility that this point may be near
1553		! the upper extrema, where the error term might be approximately parabolic and
1554		! Newton's method would converge slowly like simple bisection.
1555	696	if (slope < crv) then
1556		! if ((L_2_reg - L_min)(3.0slope - crv) > 2.0slope_crv(crv-slope)) then
1557	168	if ((L_2_reg - L_min)C3s_m_c > 2.0sxcms_c) then
1558		! There is a decent upper estimate of L from the approximate quadratic equation found
1559		! by examining the error expressions at L ~= L_2_reg and ignoring the cubic term.
1560		L_max = (slope_crv(2.0slope) - sqrt(sxcms_c**2 + &
1561	48	2.0C3s_m_c(Vol_2_reg - vol_below(K))) ) * I_3s_m_c
1562		! The line above is equivalent to:
1563		! L_max = (slope_crv(2.0slope) - sqrt(slope_crv*2(crv-slope)**2 + &
1564		! 2.0(3.0slope - crv)*(Vol_2_reg - vol_below(K))) ) / &
1565		! (3.0*slope - crv)
1566		else
1567	120	L_max = slope_crv
1568		endif
1569		else ! (slope >= crv)
1570	528	if ((1.0 - L_min)C3c_m_s > 2.0smc) then
1571		! There is a decent upper estimate of L from the approximate quadratic equation found
1572		! by examining the error expressions at L ~= 1 and ignoring the cubic term.
1573	168	L_max = ( 2.0crv - sqrt(smc2 + 2.0C3c_m_s * (Vol_open - vol_below(K))) ) * I_3c_m_s
1574		! The line above is equivalent to:
1575		! L_max = ( 2.0crv - sqrt((slope - crv)2 + 2.0(3.0crv - slope) (Vol_open - vol_below(K))) ) / &
1576		! (3.0*crv - slope)
1577		else
1578	360	L_max = 1.0
1579		endif
1580		endif
1581	696	Vol_err_max = 0.5L_max2 (slope + crv(1.0 - 4.0C1_3*L_max)) - vol_below(K)
1582		! if (Vol_err_max < 0.0) call MOM_error(FATAL, &
1583		! "Vol_err_max should never be negative in find_L_open_concave_iterative.")
1584	696	if ((Vol_err_max < abs(Vol_err)) .and. (L_max < 1.0)) then
1585		! Start with 1 bounded Newton's method step from L_max
1586	264	dVol_dL = L_max * (slope + crv(1.0 - 2.0L_max))
1587	264	L(K) = max(L_min, L_max - (vol_err_max / dVol_dL) )
1588		! else ! Could use the fact that Vol_err is known to take an iteration?
1589		endif
1590
1591		! Subsequent iterations of Newton's method do not need bounds.
1592	2976	do itt=1,max_itt
1593	2976	vol_err = 0.5L(K)2 (slope + crv(1.0 - 4.0C1_3*L(K))) - vol_below(K)
1594	2976	dVol_dL = L(K) * (slope + crv(1.0 - 2.0L(k)))
1595	2976	if (abs(vol_err) < max(1.0e-15L(K), 1.0e-25)dVol_dL) exit
1596	2280	L(K) = L(K) - (vol_err / dVol_dL)
1597		enddo
1598		endif
1599
1600		endif
1601
1602		! To check the answers.
1603		! Vol_err = 0.5(L(K)L(K))(slope + crv_3(3.0-4.0*L(K))) - vol_below(K)
1604		else ! There are two separate open regions.
1605		! vol_below(K) = slope^2/(4crv) + crv/12 - (crv/12)(1-L)^2*(1+2L)
1606		! At the deepest volume, L = slope/crv, at the top L = 1.
1607
1608		! To check the answers.
1609		! Vol_err = Vol_open - 0.25crv_3(1.0+2.0L(K)) (1.0-L(K))**2 - vol_below(K)
1610		! or equivalently:
1611		! Vol_err = Vol_open - 0.25crv_3(3.0-2.0(1.0-L(K))) (1.0-L(K))**2 - vol_below(K)
1612		! ! Note that: Vol_open = 0.25slopeslope_crv + C1_12*crv
1613		! Vol_err = 0.25slopeslope_crv + 0.25crv_3( 1.0 - (1.0 + 2.0L(K)) (1.0-L(K))**2 ) - vol_below(K)
1614		! Vol_err = 0.25crv_3L(K)*2( 3.0 - 2.0L(K) ) + 0.25slope*slope_crv - vol_below(K)
1615
1616		! Derivation of the L_max limit below:
1617		! Vol_open - vol_below(K) = 0.25crv_3(3.0-2.0(1.0-L(K))) (1.0-L(K))**2
1618		! (3.0-2.0(1.0-L(K))) (1.0-L(K))*2 = (Vol_open - vol_below(K)) / (0.25crv_3)
1619		! When 1-L(K) << 1:
1620		! 3.0 * (1.0-L_max)*2 = (Vol_open - vol_below(K)) / (0.25crv_3)
1621		! (1.0-L_max)*2 = (Vol_open - vol_below(K)) / (0.25crv)
1622
1623		! Derivation of the L_min limit below:
1624		! Vol_err = 0.25crv_3L(K)*2( 3.0 - 2.0L(K) ) + 0.25slope*slope_crv - vol_below(K)
1625		! crvL(K)2( 1.0 - 2.0C1_3L(K) ) = 4.0vol_below(K) - slopeslope_crv
1626		! When L(K) << 1:
1627		! crvL_min2 = 4.0vol_below(K) - slope*slope_crv
1628		! L_min = sqrt((4.0vol_below(K) - slopeslope_crv)/crv)
1629		! Noting that L(K) >= slope_crv, when L(K)-slope_crv << 1:
1630		! (crv + 2.0C1_3slope)L_min2 = 4.0vol_below(K) - slope*slope_crv
1631		! L_min = sqrt((4.0vol_below(K) - slopeslope_crv)/(crv + 2.0C1_3slope))
1632
1633	1992	if (vol_below(K) <= Vol_inflect_2) then
1634		! Newton's Method would converge monotonically from above, but overshoot from below.
1635	1776	L_min = max(L(K+1), L_2_reg) ! L_2_reg = slope_crv
1636		! This under-estimate of L(K) is accurate for L ~= slope_crv:
1637	1776	if ((4.0vol_below(K) - slopeslope_crv) > (crv + 2.0C1_3slope)L_min*2) &
1638	960	L_min = max(L_min, sqrt((4.0vol_below(K) - slopeslope_crv) / (crv + 2.0C1_3slope)))
1639	1776	L_max = 0.5 ! = L_inflect_2
1640
1641		! Starting with L_min below L_inflect_2, only the first overshooting iteration of Newton's
1642		! method needs bounding.
1643	1776	L(k) = L_min
1644	1776	Vol_err = crv_3L(K)2( 0.75 - 0.5*L(K) ) + (slope2_4crv - vol_below(K))
1645
1646		! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L_min is already the best solution.
1647	1776	if (vol_err < 0.0) then
1648	1776	dVol_dL = 0.5crv (L(K) * (1.0 - L(K)))
1649	1776	if (L(K)dVol_dL >= vol_err + L_maxdVol_dL) then
1650	0	L(K) = L_max
1651		else
1652	1776	L(K) = L(K) - (vol_err / dVol_dL)
1653		endif
1654		! Subsequent iterations of Newton's method do not need bounds.
1655	5619	do itt=1,max_itt
1656	5619	Vol_err = crv_3 * (L(K)*2 (0.75 - 0.5*L(K))) + (slope2_4crv - vol_below(K))
1657	5619	dVol_dL = 0.5crv (L(K)*(1.0 - L(K)))
1658	5619	if (abs(vol_err) < max(1.0e-15L(K), 1.0e-25)dVol_dL) exit
1659	3843	L(K) = L(K) - (vol_err / dVol_dL)
1660		enddo
1661		endif
1662		else ! (vol_below(K) > Vol_inflect_2)
1663		! Newton's Method would converge monotonically from below, but overshoots from above, and
1664		! we may need to deal with the case where we are close to the peak of the topography.
1665	216	L_min = max(L(K+1), 0.5)
1666	216	L(k) = L_min
1667
1668	216	Vol_err = crv_3 * (L(K)*2 ( 0.75 - 0.5*L(K))) + (slope2_4crv - vol_below(K))
1669		! If vol_err is 0 or positive (perhaps due to roundoff in L(K+1)), L(k) is already the best solution.
1670	216	if (Vol_err < 0.0) then
1671		! This over-estimate of L(K) is accurate for L ~= 1:
1672	216	L_max = 1.0 - sqrt( (Vol_open - vol_below(K)) * C4_crv )
1673	216	Vol_err_max = crv_3 * (L_max*2 ( 0.75 - 0.5*L_max)) + (slope2_4crv - vol_below(K))
1674		! if (Vol_err_max < 0.0) call MOM_error(FATAL, &
1675		! "Vol_err_max should never be negative in find_L_open_concave_iterative.")
1676	216	if ((Vol_err_max < abs(Vol_err)) .and. (L_max < 1.0)) then
1677		! Start with 1 bounded Newton's method step from L_max
1678	216	dVol_dL = 0.5crv (L_max * (1.0 - L_max))
1679	216	L(K) = max(L_min, L_max - (vol_err_max / dVol_dL) )
1680		! else ! Could use the fact that Vol_err is known to take an iteration?
1681		endif
1682
1683		! Subsequent iterations of Newton's method do not need bounds.
1684	816	do itt=1,max_itt
1685	816	Vol_err = crv_3 * (L(K)*2 ( 0.75 - 0.5*L(K))) + (slope2_4crv - vol_below(K))
1686	816	dVol_dL = 0.5crv (L(K) * (1.0 - L(K)))
1687	816	if (abs(vol_err) < max(1.0e-15L(K), 1.0e-25)dVol_dL) exit
1688	600	L(K) = L(K) - (vol_err / dVol_dL)
1689		enddo
1690		endif
1691		endif
1692
1693		endif
1694		enddo ! k loop to determine L(K) in the concave case
1695
1696	11664	end subroutine find_L_open_concave_iterative
1697
1698
1699
1700		!> Test the validity the normalized open lengths of each interface for concave bathymetry (from the ocean perspective)
1701		!! by evaluating and returing the relevant cubic equations.
1702	0	pure subroutine test_L_open_concave(vol_below, D_vel, Dp, Dm, L, vol_err, GV)
1703		!$omp declare target
1704		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1705		real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by
1706		!! the full horizontal area of a velocity cell [Z ~> m]
1707		real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m]
1708		real, intent(in) :: Dp !< The larger of the two depths at the edge
1709		!! of a velocity cell [Z ~> m]
1710		real, intent(in) :: Dm !< The smaller of the two depths at the edge
1711		!! of a velocity cell [Z ~> m]
1712		real, dimension(SZK_(GV)+1), intent(in) :: L !< The fraction of the full cell width that is open at
1713		!! the depth of each interface [nondim]
1714		real, dimension(SZK_(GV)+1), intent(out) :: vol_err !< The difference between vol_below and the
1715		!! value obtained from using L in the cubic equation [Z ~> m]
1716
1717		! Local variables
1718		real :: crv ! crv is the curvature of the bottom depth across a
1719		! cell, times the cell width squared [Z ~> m].
1720		real :: crv_3 ! crv/3 [Z ~> m].
1721		real :: slope ! The absolute value of the bottom depth slope across
1722		! a cell times the cell width [Z ~> m].
1723
1724		! The following "volumes" have units of vertical heights because they are normalized
1725		! by the full horizontal area of a velocity cell.
1726		real :: Vol_open ! The cell volume above which the face is fully is open [Z ~> m].
1727		real :: Vol_2_reg ! The cell volume above which there are two separate
1728		! open areas that must be integrated [Z ~> m].
1729		real :: L_2_reg ! The value of L when vol_below is Vol_2_reg [nondim]
1730
1731		! The following combinations of slope and crv are reused across layers, and hence are pre-calculated
1732		! for efficiency. All are non-negative.
1733		real :: slope_crv ! The slope divided by the curvature [nondim]
1734		! These are only used if the curvature exceeds the slope.
1735		real :: slope2_4crv ! A quarter of the slope squared divided by the curvature [Z ~> m]
1736
1737		real, parameter :: C1_3 = 1.0 / 3.0, C1_12 = 1.0 / 12.0 ! Rational constants [nondim]
1738		integer :: K, nz
1739
1740	0	nz = GV%ke
1741
1742		! Each cell extends from x=-1/2 to 1/2, and has a topography
1743		! given by D(x) = crvx^2 + slopex + D_vel - crv/12.
1744
1745	0	crv_3 = (Dp + Dm - 2.0D_vel) ; crv = 3.0crv_3
1746	0	slope = Dp - Dm
1747
1748		! Calculate the volume above which the entire cell is open and the volume at which the
1749		! equation that is solved for L changes because there are two separate open regions.
1750	0	if (slope >= crv) then
1751	0	Vol_open = D_vel - Dm ; Vol_2_reg = Vol_open
1752	0	L_2_reg = 1.0
1753	0	if (crv + slope >= 4.0*crv) then
1754	0	slope_crv = 1.0
1755		else
1756	0	slope_crv = slope / crv
1757		endif
1758		else
1759	0	slope_crv = slope / crv
1760	0	Vol_open = 0.25slopeslope_crv + C1_12*crv
1761	0	Vol_2_reg = 0.5slope_crv2 (crv - C1_3*slope)
1762	0	L_2_reg = slope_crv
1763		endif
1764	0	slope2_4crv = 0.25 * slope * slope_crv
1765
1766		! Determine the volume error based on the normalized open length (L) at each interface.
1767	0	Vol_err(nz+1) = 0.0
1768	0	do K=nz,1,-1
1769	0	if (L(K) >= 1.0) then
1770	0	Vol_err(K) = max(Vol_open - vol_below(K), 0.0)
1771	0	elseif (L(K) <= L_2_reg) then
1772	0	vol_err(K) = 0.5L(K)2 (slope + crv(1.0 - 4.0C1_3*L(K))) - vol_below(K)
1773		else ! There are two separate open regions.
1774	0	Vol_err(K) = crv_3 * (L(K)*2 ( 0.75 - 0.5*L(K))) + (slope2_4crv - vol_below(K))
1775		endif
1776		enddo ! k loop to determine L(K) in the concave case
1777
1778	0	end subroutine test_L_open_concave
1779
1780
1781		!> Determine the normalized open length of each interface for convex bathymetry (from the ocean
1782		!! perspective) using Newton's method iterations. In this case there is a single open region
1783		!! with the minimum depth at one edge of the cell.
1784	105672	pure subroutine find_L_open_convex(vol_below, D_vel, Dp, Dm, L, GV, US, CS)
1785		!$omp declare target
1786		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1787		real, dimension(SZK_(GV)+1), intent(in) :: vol_below !< The volume below each interface, normalized by
1788		!! the full horizontal area of a velocity cell [Z ~> m]
1789		real, intent(in) :: D_vel !< The average bottom depth at a velocity point [Z ~> m]
1790		real, intent(in) :: Dp !< The larger of the two depths at the edge
1791		!! of a velocity cell [Z ~> m]
1792		real, intent(in) :: Dm !< The smaller of the two depths at the edge
1793		!! of a velocity cell [Z ~> m]
1794		real, dimension(SZK_(GV)+1), intent(out) :: L !< The fraction of the full cell width that is open at
1795		!! the depth of each interface [nondim]
1796		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1797		type(set_visc_CS), intent(in) :: CS !< The control structure returned by a previous
1798		!! call to set_visc_init.
1799
1800		! Local variables
1801		real :: crv ! crv is the curvature of the bottom depth across a
1802		! cell, times the cell width squared [Z ~> m].
1803		real :: crv_3 ! crv/3 [Z ~> m].
1804		real :: slope ! The absolute value of the bottom depth slope across
1805		! a cell times the cell width [Z ~> m].
1806		! All of the following "volumes" have units of vertical heights because they are normalized
1807		! by the full horizontal area of a velocity cell.
1808		real :: Vol_err ! The error in the volume with the latest estimate of
1809		! L, or the error for the interface below [Z ~> m].
1810		real :: Vol_quit ! The volume error below which to quit iterating [Z ~> m].
1811		real :: Vol_tol ! A volume error tolerance [Z ~> m].
1812		real :: Vol_open ! The cell volume above which the face is fully open [Z ~> m].
1813		real :: Vol_direct ! With less than Vol_direct [Z ~> m], there is a direct
1814		! solution of a cubic equation for L.
1815		real :: Vol_err_max ! The volume error for the upper bound on the correct value for L [Z ~> m]
1816		real :: Vol_err_min ! The volume error for the lower bound on the correct value for L [Z ~> m]
1817		real :: Vol_0 ! A deeper volume with known width L0 [Z ~> m].
1818		real :: dVol ! vol - Vol_0 [Z ~> m].
1819		real :: dV_dL2 ! The partial derivative of volume with L squared
1820		! evaluated at L=L0 [Z ~> m].
1821		real :: L_direct ! The value of L above volume Vol_direct [nondim].
1822		real :: L_max, L_min ! Upper and lower bounds on the correct value for L [nondim].
1823		real :: L0 ! The value of L above volume Vol_0 [nondim].
1824		real :: Iapb, Ibma_2 ! Combinations of crv (a) and slope (b) [Z-1 ~> m-1]
1825		real :: C24_crv ! 24/crv [Z-1 ~> m-1].
1826		real :: curv_tol ! Numerator of curvature cubed, used to estimate
1827		! accuracy of a single L(:) Newton iteration [Z5 ~> m5]
1828		real, parameter :: C1_3 = 1.0/3.0, C1_6 = 1.0/6.0 ! Rational constants [nondim]
1829		logical :: use_L0, do_one_L_iter ! Control flags for L(:) Newton iteration
1830		integer :: K, nz, itt
1831		integer, parameter:: maxitt = 20
1832
1833	105672	nz = GV%ke
1834
1835		! Each cell extends from x=-1/2 to 1/2, and has a topography
1836		! given by D(x) = crvx^2 + slopex + D_vel - crv/12.
1837	105672	crv_3 = (Dp + Dm - 2.0D_vel) ; crv = 3.0crv_3
1838	105672	slope = Dp - Dm
1839
1840		! Calculate the volume above which the entire cell is open and the volume at which the
1841		! equation that is solved for L changes because there is a direct solution.
1842	105672	Vol_open = D_vel - Dm
1843	105672	if (slope >= -crv) then
1844	89760	Iapb = 1.0e30*US%Z_to_m ; if (slope+crv /= 0.0) Iapb = 1.0/(crv+slope)
1845	89760	Vol_direct = 0.0 ; L_direct = 0.0 ; C24_crv = 0.0
1846		else
1847	15912	C24_crv = 24.0/crv ; Iapb = 1.0/(crv+slope)
1848	15912	L_direct = 1.0 + slope/crv ! L_direct < 1 because crv < 0
1849	15912	Vol_direct = -C1_6crvL_direct**3
1850		endif
1851	105672	Ibma_2 = 2.0 / (slope - crv)
1852
1853	105672	if (CS%answer_date < 20190101) Vol_quit = (0.9*GV%Angstrom_Z + GV%dZ_subroundoff)
1854
1855	105672	L(nz+1) = 0.0 ; Vol_err = 0.0
1856		! Determine the normalized open length (L) at each interface.
1857	8031072	do K=nz,1,-1
1858	8031072	if (vol_below(K) >= Vol_open) then
1859	5431235	L(K) = 1.0
1860	2494165	elseif (vol_below(K) <= Vol_direct) then
1861	351360	if (CS%answer_date < 20260501) then
1862	0	L(K) = (-0.25C24_crvvol_below(K))**C1_3
1863		else
1864	351360	L(K) = cuberoot(-0.25C24_crvvol_below(K))
1865		endif
1866		else
1867		! x_R is at 1/2 but x_L is in the interior & L is found by iteratively solving
1868		! vol_below(K) = 0.5L^2(slope + crv/3*(3-4L))
1869
1870		! Vol_err = 0.5(L(K+1)L(K+1))(slope + crv_3(3.0-4.0*L(K+1))) - vol_below(K+1)
1871		! Change to ...
1872		! if (min(vol_below(K+1) + Vol_err, vol_below(K)) <= Vol_direct) then ?
1873	2142805	if (vol_below(K+1) + Vol_err <= Vol_direct) then
1874	98556	L0 = L_direct ; Vol_0 = Vol_direct
1875		else
1876	2044249	L0 = L(K+1) ; Vol_0 = vol_below(K+1) + Vol_err
1877		! Change to Vol_0 = min(vol_below(K+1) + Vol_err, vol_below(K)) ?
1878		endif
1879
1880		! Try a relatively simple solution that usually works well
1881		! for massless layers.
1882	2142805	dV_dL2 = 0.5(slope+crv) - crvL0 ; dVol = (vol_below(K)-Vol_0)
1883		! dV_dL2 = 0.5(slope+crv) - crvL0 ; dVol = max(vol_below(K)-Vol_0, 0.0)
1884
1885	2142805	use_L0 = .false.
1886	2142805	do_one_L_iter = .false.
1887	2142805	if (CS%answer_date < 20190101) then
1888		curv_tol = GV%Angstrom_ZdV_dL2*2 &
1889	0	* (0.25 * dV_dL2 * GV%Angstrom_Z - crv * L0 * dVol)
1890	0	do_one_L_iter = (crv * crv * dVol**3) < curv_tol
1891		else
1892		! The following code is more robust when GV%Angstrom_H=0, but
1893		! it changes answers.
1894	2142805	use_L0 = (dVol <= 0.)
1895
1896	2142805	Vol_tol = max(0.5 * GV%Angstrom_Z + GV%dZ_subroundoff, 1e-14 * vol_below(K))
1897	2142805	Vol_quit = max(0.9 * GV%Angstrom_Z + GV%dZ_subroundoff, 1e-14 * vol_below(K))
1898
1899		curv_tol = Vol_tol * dV_dL2**2 &
1900	2142805	* (dV_dL2 * Vol_tol - 2.0 * crv * L0 * dVol)
1901	2142805	do_one_L_iter = (crv * crv * dVol**3) < curv_tol
1902		endif
1903
1904	2142805	if (use_L0) then
1905	0	L(K) = L0
1906	0	Vol_err = 0.5(L(K)L(K))(slope + crv_3(3.0-4.0*L(K))) - vol_below(K)
1907	2142805	elseif (do_one_L_iter) then
1908		! One iteration of Newton's method should give an estimate
1909		! that is accurate to within Vol_tol.
1910	0	L(K) = sqrt(L0*L0 + dVol / dV_dL2)
1911	0	Vol_err = 0.5(L(K)L(K))(slope + crv_3(3.0-4.0*L(K))) - vol_below(K)
1912		else
1913	2142805	if (dV_dL2(1.0-L0L0) < dVol + &
1914		dV_dL2 * (Vol_open - vol_below(K))*Ibma_2) then
1915	101353	L_max = sqrt(1.0 - (Vol_open - vol_below(K))*Ibma_2)
1916		else
1917	2041452	L_max = sqrt(L0*L0 + dVol / dV_dL2)
1918		endif
1919	2142805	L_min = sqrt(L0L0 + dVol / (0.5(slope+crv) - crv*L_max))
1920
1921	2142805	Vol_err_min = 0.5(L_min2)(slope + crv_3(3.0-4.0L_min)) - vol_below(K)
1922	2142805	Vol_err_max = 0.5(L_max2)(slope + crv_3(3.0-4.0L_max)) - vol_below(K)
1923		! if ((abs(Vol_err_min) <= Vol_quit) .or. (Vol_err_min >= Vol_err_max)) then
1924	2142805	if (abs(Vol_err_min) <= Vol_quit) then
1925	0	L(K) = L_min ; Vol_err = Vol_err_min
1926		else
1927		L(K) = sqrt((L_min*2Vol_err_max - L_max*2Vol_err_min) / &
1928	2142805	(Vol_err_max - Vol_err_min))
1929	2807612	do itt=1,maxitt
1930	2807612	Vol_err = 0.5(L(K)L(K))(slope + crv_3(3.0-4.0*L(K))) - vol_below(K)
1931	2807612	if (abs(Vol_err) <= Vol_quit) exit
1932		! Take a Newton's method iteration. This equation has proven
1933		! robust enough not to need bracketing.
1934	664807	L(K) = L(K) - Vol_err / (L(K)* (slope + crv - 2.0crvL(K)))
1935		! This would be a Newton's method iteration for L^2:
1936		! L(K) = sqrt(L(K)L(K) - Vol_err / (0.5(slope+crv) - crv*L(K)))
1937		enddo
1938		endif ! end of iterative solver
1939		endif ! end of 1-boundary alternatives.
1940		endif ! end of 0, 1- and 2- boundary cases.
1941		enddo ! k loop to determine L(K) in the convex case
1942
1943	105672	end subroutine find_L_open_convex
1944
1945		!> This subroutine finds a thickness-weighted value of v at the u-points.
1946	2867989	pure function set_v_at_u(v, h, G, GV, i, j, k, mask2dCv, OBC)
1947		!$omp declare target
1948		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1949		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1950		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1951		intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
1952		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1953		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1954		integer, intent(in) :: i !< The i-index of the u-location to work on.
1955		integer, intent(in) :: j !< The j-index of the u-location to work on.
1956		integer, intent(in) :: k !< The k-index of the u-location to work on.
1957		real, dimension(SZI_(G),SZJB_(G)),&
1958		intent(in) :: mask2dCv !< A multiplicative mask of the v-points [nondim]
1959		type(ocean_OBC_type), pointer :: OBC !< A pointer to an open boundary condition structure
1960		real :: set_v_at_u !< The return value of v at u points points in the
1961		!! same units as u, i.e. [L T-1 ~> m s-1] or other units.
1962
1963		! This subroutine finds a thickness-weighted value of v at the u-points.
1964		real :: hwt(0:1,-1:0) ! Masked weights used to average u onto v [H ~> m or kg m-2].
1965		real :: hwt_tot ! The sum of the masked thicknesses [H ~> m or kg m-2].
1966		integer :: i0, j0, i1, j1
1967
1968	20075923	do j0 = -1,0 ; do i0 = 0,1 ; i1 = i+i0 ; J1 = J+j0
1969	17207934	hwt(i0,j0) = (h(i1,j1,k) + h(i1,j1+1,k)) * mask2dCv(i1,J1)
1970		enddo ; enddo
1971
1972	2867989	if (associated(OBC)) then ; if (OBC%number_of_segments > 0) then
1973	0	do j0 = -1,0 ; do i0 = 0,1 ; if (OBC%segnum_v(i+i0,J+j0) /= 0) then
1974	0	i1 = i+i0 ; J1 = J+j0
1975	0	if (OBC%segnum_v(i1,j1) > 0) then ! OBC_DIRECTION_N
1976	0	hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dCv(i1,J1)
1977	0	elseif (OBC%segnum_v(i1,J1) < 0) then ! OBC_DIRECTION_S
1978	0	hwt(i0,j0) = 2.0 * h(i1,J1+1,k) * mask2dCv(i1,J1)
1979		endif
1980		endif ; enddo ; enddo
1981		endif ; endif
1982
1983	2867989	hwt_tot = (hwt(0,-1) + hwt(1,0)) + (hwt(1,-1) + hwt(0,0))
1984	2867989	set_v_at_u = 0.0
1985	2867989	if (hwt_tot > 0.0) set_v_at_u = &
1986		(((hwt(0,0) * v(i,J,k)) + (hwt(1,-1) * v(i+1,J-1,k))) + &
1987	2867989	((hwt(1,0) * v(i+1,J,k)) + (hwt(0,-1) * v(i,J-1,k)))) / hwt_tot
1988
1989	2867989	end function set_v_at_u
1990
1991		!> This subroutine finds a thickness-weighted value of u at the v-points.
1992	2803353	pure function set_u_at_v(u, h, G, GV, i, j, k, mask2dCu, OBC)
1993		!$omp declare target
1994		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1995		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1996		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1997		intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1] or other units.
1998		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1999		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
2000		integer, intent(in) :: i !< The i-index of the u-location to work on.
2001		integer, intent(in) :: j !< The j-index of the u-location to work on.
2002		integer, intent(in) :: k !< The k-index of the u-location to work on.
2003		real, dimension(SZIB_(G),SZJ_(G)), &
2004		intent(in) :: mask2dCu !< A multiplicative mask of the u-points [nondim]
2005		type(ocean_OBC_type), pointer :: OBC !< A pointer to an open boundary condition structure
2006		real :: set_u_at_v !< The return value of u at v points in the
2007		!! same units as u, i.e. [L T-1 ~> m s-1] or other units.
2008
2009		! This subroutine finds a thickness-weighted value of u at the v-points.
2010		real :: hwt(-1:0,0:1) ! Masked weights used to average u onto v [H ~> m or kg m-2].
2011		real :: hwt_tot ! The sum of the masked thicknesses [H ~> m or kg m-2].
2012		integer :: i0, j0, i1, j1
2013
2014	19623471	do j0 = 0,1 ; do i0 = -1,0 ; I1 = I+i0 ; j1 = j+j0
2015	16820118	hwt(i0,j0) = (h(i1,j1,k) + h(i1+1,j1,k)) * mask2dCu(I1,j1)
2016		enddo ; enddo
2017
2018	2803353	if (associated(OBC)) then ; if (OBC%number_of_segments > 0) then
2019	0	do j0 = 0,1 ; do i0 = -1,0 ; if ((OBC%segnum_u(I+i0,j+j0) /= 0)) then
2020	0	I1 = I+i0 ; j1 = j+j0
2021	0	if (OBC%segnum_u(I1,j1) > 0) then ! OBC_DIRECTION_E
2022	0	hwt(i0,j0) = 2.0 * h(I1,j1,k) * mask2dCu(I1,j1)
2023	0	elseif (OBC%segnum_u(I1,j1) < 0) then ! OBC_DIRECTION_W
2024	0	hwt(i0,j0) = 2.0 * h(I1+1,j1,k) * mask2dCu(I1,j1)
2025		endif
2026		endif ; enddo ; enddo
2027		endif ; endif
2028
2029	2803353	hwt_tot = (hwt(-1,0) + hwt(0,1)) + (hwt(0,0) + hwt(-1,1))
2030	2803353	set_u_at_v = 0.0
2031	2803353	if (hwt_tot > 0.0) set_u_at_v = &
2032		(((hwt(0,0) * u(I,j,k)) + (hwt(-1,1) * u(I-1,j+1,k))) + &
2033	2803353	((hwt(-1,0) * u(I-1,j,k)) + (hwt(0,1) * u(I,j+1,k)))) / hwt_tot
2034
2035	2803353	end function set_u_at_v
2036
2037		!> Calculates the thickness of the surface boundary layer for applying an elevated viscosity.
2038		!!
2039		!! A bulk Richardson criterion or the thickness of the topmost NKML layers (with a bulk mixed layer)
2040		!! are currently used. The thicknesses are given in terms of fractional layers, so that this
2041		!! thickness will move as the thickness of the topmost layers change.
2042	24	subroutine set_viscous_ML(u, v, h, tv, forces, visc, dt, G, GV, US, CS)
2043		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
2044		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2045		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2046		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
2047		intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
2048		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
2049		intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
2050		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2051		intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
2052		type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
2053		!! thermodynamic fields. Absent fields have
2054		!! NULL pointers.
2055		type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
2056		type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
2057		!! related fields.
2058		real, intent(in) :: dt !< Time increment [T ~> s].
2059		type(set_visc_CS), intent(inout) :: CS !< The control structure returned by a previous
2060		!! call to set_visc_init.
2061
2062		! Local variables
2063		real, dimension(SZIB_(G)) :: &
2064	48	htot, & ! The total thickness of the layers that are within the
2065		! surface mixed layer [H ~> m or kg m-2].
2066	48	dztot, & ! The distance from the surface to the bottom of the layers that are
2067		! within the surface mixed layer [Z ~> m]
2068	48	Thtot, & ! The integrated temperature of layers that are within the
2069		! surface mixed layer [H C ~> m degC or kg degC m-2].
2070	48	Shtot, & ! The integrated salt of layers that are within the
2071		! surface mixed layer [H S ~> m ppt or kg ppt m-2].
2072	48	SpV_htot, & ! Running sum of thickness times specific volume [H R-1 ~> m4 kg-1 or m]
2073	48	Rhtot, & ! The integrated density of layers that are within the surface mixed layer
2074		! [H R ~> kg m-2 or kg2 m-5]. Rhtot is only used if no
2075		! equation of state is used.
2076	48	uhtot, & ! The depth integrated zonal velocity within the surface
2077		! mixed layer [H L T-1 ~> m2 s-1 or kg m-1 s-1].
2078	48	vhtot, & ! The depth integrated meridional velocity within the surface
2079		! mixed layer [H L T-1 ~> m2 s-1 or kg m-1 s-1].
2080	48	Idecay_len_TKE, & ! The inverse of a turbulence decay length scale [H-1 ~> m-1 or m2 kg-1].
2081	48	dR_dT, & ! Partial derivative of the density at the base of layer nkml
2082		! (roughly the base of the mixed layer) with temperature [R C-1 ~> kg m-3 degC-1].
2083	48	dR_dS, & ! Partial derivative of the density at the base of layer nkml
2084		! (roughly the base of the mixed layer) with salinity [R S-1 ~> kg m-3 ppt-1].
2085	48	dSpV_dT, & ! Partial derivative of the specific volume at the base of layer nkml
2086		! (roughly the base of the mixed layer) with temperature [R-1 C-1 ~> m3 kg-1 degC-1].
2087	48	dSpV_dS, & ! Partial derivative of the specific volume at the base of layer nkml
2088		! (roughly the base of the mixed layer) with salinity [R-1 S-1 ~> m3 kg-1 ppt-1].
2089	48	ustar, & ! The surface friction velocity under ice shelves [H T-1 ~> m s-1 or kg m-2 s-1].
2090	48	press, & ! The pressure at which dR_dT and dR_dS are evaluated [R L2 T-2 ~> Pa].
2091	48	T_EOS, & ! The potential temperature at which dR_dT and dR_dS are evaluated [C ~> degC]
2092	48	S_EOS ! The salinity at which dR_dT and dR_dS are evaluated [S ~> ppt].
2093	48	real :: dz(SZI_(G),SZJ_(G),SZK_(GV)) ! Height change across layers [Z ~> m]
2094		real, dimension(SZIB_(G),SZJ_(G)) :: &
2095	48	mask_u ! A mask that disables any contributions from u points that
2096		! are land or past open boundary conditions [nondim], 0 or 1.
2097		real, dimension(SZI_(G),SZJB_(G)) :: &
2098	48	mask_v ! A mask that disables any contributions from v points that
2099		! are land or past open boundary conditions [nondim], 0 or 1.
2100	48	real :: U_star_2d(SZI_(G),SZJ_(G)) ! The wind friction velocity in thickness-based units,
2101		! calculated using the Boussinesq reference density or the time-evolving
2102		! surface density in non-Boussinesq mode [H T-1 ~> m s-1 or kg m-2 s-1]
2103	48	real :: h_at_vel(SZIB_(G),SZK_(GV))! Layer thickness at velocity points,
2104		! using an upwind-biased second order accurate estimate based
2105		! on the previous velocity direction [H ~> m or kg m-2].
2106	48	real :: dz_at_vel(SZIB_(G),SZK_(GV)) ! Vertical extent of a layer at velocity points,
2107		! using an upwind-biased second order accurate estimate based
2108		! on the previous velocity direction [Z ~> m].
2109	48	integer :: k_massive(SZIB_(G)) ! The k-index of the deepest layer yet found
2110		! that has more than h_tiny thickness and will be in the
2111		! viscous mixed layer.
2112		real :: Uh2 ! The squared magnitude of the difference between the velocity
2113		! integrated through the mixed layer and the velocity of the
2114		! interior layer layer times the depth of the mixed layer
2115		! [H2 L2 T-2 ~> m4 s-2 or kg2 m-2 s-2].
2116		real :: htot_vel ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
2117		real :: hwtot ! Sum of the thicknesses used to calculate
2118		! the near-bottom velocity magnitude [H ~> m or kg m-2].
2119		real :: hutot ! Running sum of thicknesses times the velocity
2120		! magnitudes [H L T-1 ~> m2 s-1 or kg m-1 s-1].
2121		real :: hweight ! The thickness of a layer that is within Hbbl
2122		! of the bottom [H ~> m or kg m-2].
2123		real :: tbl_thick ! The thickness of the top boundary layer [Z ~> m].
2124
2125		real :: hlay ! The layer thickness at velocity points [H ~> m or kg m-2].
2126		real :: I_2hlay ! 1 / 2*hlay [H-1 ~> m-1 or m2 kg-1].
2127		real :: T_lay ! The layer temperature at velocity points [C ~> degC].
2128		real :: S_lay ! The layer salinity at velocity points [S ~> ppt].
2129		real :: Rlay ! The layer potential density at velocity points [R ~> kg m-3].
2130		real :: Rlb ! The potential density of the layer below [R ~> kg m-3].
2131		real :: v_at_u ! The meridional velocity at a zonal velocity point [L T-1 ~> m s-1].
2132		real :: u_at_v ! The zonal velocity at a meridional velocity point [L T-1 ~> m s-1].
2133		real :: gHprime ! The mixed-layer internal gravity wave speed squared, based
2134		! on the mixed layer thickness and density difference across
2135		! the base of the mixed layer [L2 T-2 ~> m2 s-2].
2136		real :: RiBulk ! The bulk Richardson number below which water is in the
2137		! viscous mixed layer, including reduction for turbulent decay [nondim]
2138		real :: dt_Rho0 ! The time step divided by the conversion from the layer
2139		! thickness to layer mass [T H Z-1 R-1 ~> s m3 kg-1 or s].
2140		real :: g_H_Rho0 ! The gravitational acceleration times the conversion from H to m divided
2141		! by the mean density [L2 T-2 H-1 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2].
2142		real :: ustarsq ! 400 times the square of ustar, times
2143		! Rho0 divided by G_Earth and the conversion
2144		! from m to thickness units [H R ~> kg m-2 or kg2 m-5].
2145		real :: cdrag_sqrt ! Square root of the drag coefficient [nondim].
2146		real :: cdrag_sqrt_H ! Square root of the drag coefficient, times a unit conversion
2147		! factor from lateral lengths to layer thicknesses [H L-1 ~> nondim or kg m-3].
2148		real :: cdrag_sqrt_H_RL ! Square root of the drag coefficient, times a unit conversion factor from
2149		! density times lateral lengths to layer thicknesses [H L-1 R-1 ~> m3 kg-1 or nondim]
2150		real :: oldfn ! The integrated energy required to
2151		! entrain up to the bottom of the layer,
2152		! divided by G_Earth [H R ~> kg m-2 or kg2 m-5].
2153		real :: Dfn ! The increment in oldfn for entraining
2154		! the layer [H R ~> kg m-2 or kg2 m-5].
2155		real :: frac_used ! The fraction of the present layer that contributes to Dh and Ddz [nondim]
2156		real :: Dh ! The increment in layer thickness from the present layer [H ~> m or kg m-2].
2157		real :: Ddz ! The increment in height change from the present layer [Z ~> m].
2158	48	real :: u2_bg(SZIB_(G)) ! The square of an assumed background velocity, for
2159		! calculating the mean magnitude near the top for use in
2160		! the quadratic surface drag [L2 T-2 ~> m2 s-2].
2161		real :: h_tiny ! A very small thickness [H ~> m or kg m-2]. Layers that are less than
2162		! h_tiny can not be the deepest in the viscous mixed layer.
2163		real :: absf ! The absolute value of f averaged to velocity points [T-1 ~> s-1].
2164		real :: U_star ! The friction velocity at velocity points [H T-1 ~> m s-1 or kg m-2 s-1].
2165		real :: h_neglect ! A thickness that is so small it is usually lost
2166		! in roundoff and can be neglected [H ~> m or kg m-2].
2167		real :: dz_neglect ! A vertical distance that is so small it is usually lost
2168		! in roundoff and can be neglected [Z ~> m].
2169		real :: Rho0x400_G ! 400*Rho0/G_Earth, times unit conversion factors
2170		! [R T2 H-1 ~> kg s2 m-4 or s2 m-1].
2171		! The 400 is a constant proposed by Killworth and Edwards, 1999.
2172		real :: ustar1 ! ustar [H T-1 ~> m s-1 or kg m-2 s-1]
2173		real :: h2f2 ! (h2f)^2 [H2 T-2 ~> m2 s-2 or kg2 m-4 s-2]
2174	24	logical :: use_EOS, do_any, do_any_shelf, do_i(SZIB_(G))
2175		logical :: nonBous_ML ! If true, use the non-Boussinesq form of some energy and
2176		! stratification calculations.
2177		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, K2, nkmb, nkml, n
2178		type(ocean_OBC_type), pointer :: OBC => NULL()
2179
2180	24	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
2181	24	Isq = G%isc-1 ; Ieq = G%IecB ; Jsq = G%jsc-1 ; Jeq = G%JecB
2182	24	nkmb = GV%nk_rho_varies ; nkml = GV%nkml
2183
2184	24	if (.not.CS%initialized) call MOM_error(FATAL,"MOM_set_viscosity(visc_ML): "//&
2185	0	"Module must be initialized before it is used.")
2186
2187		! TODO: Remove this check and move it outside of the function call.
2188	24	if (.not.(CS%dynamic_viscous_ML .or. associated(forces%frac_shelf_u) .or. &
2189	24	associated(forces%frac_shelf_v)) ) return
2190
2191		! NOTE: Requried since this is called by the GPU-enabled dycore, but it could
2192		! also be implicitly fixing other functions.
2193		!$omp target update from(u, v)
2194
2195	0	Rho0x400_G = 400.0*(GV%H_to_RZ / GV%g_Earth_Z_T2)
2196	0	cdrag_sqrt = sqrt(CS%cdrag)
2197	0	cdrag_sqrt_H = cdrag_sqrt * US%L_to_m * GV%m_to_H
2198	0	cdrag_sqrt_H_RL = cdrag_sqrt * US%L_to_Z * GV%RZ_to_H
2199
2200	0	OBC => CS%OBC
2201	0	use_EOS = associated(tv%eqn_of_state)
2202	0	nonBous_ML = allocated(tv%SpV_avg)
2203	0	dt_Rho0 = dt / GV%H_to_RZ
2204	0	h_neglect = GV%H_subroundoff
2205	0	h_tiny = 2.0*GV%Angstrom_H + h_neglect
2206	0	dz_neglect = GV%dZ_subroundoff
2207	0	g_H_Rho0 = (GV%g_Earth*GV%H_to_Z) / (GV%Rho0)
2208
2209	0	if (associated(forces%frac_shelf_u) .neqv. associated(forces%frac_shelf_v)) &
2210		call MOM_error(FATAL, "set_viscous_ML: one of forces%frac_shelf_u and "//&
2211	0	"forces%frac_shelf_v is associated, but the other is not.")
2212
2213		! Extract the friction velocity from the forcing type.
2214	0	call find_ustar(forces, tv, U_star_2d, G, GV, US, halo=1, H_T_units=.true.)
2215
2216	0	if (associated(forces%frac_shelf_u)) then
2217		! This configuration has ice shelves, and the appropriate variables need to be
2218		! allocated. If the arrays have already been allocated, these calls do nothing.
2219	0	if (.not.allocated(visc%taux_shelf)) &
2220	0	allocate(visc%taux_shelf(G%IsdB:G%IedB, G%jsd:G%jed), source=0.0)
2221	0	if (.not.allocated(visc%tauy_shelf)) &
2222	0	allocate(visc%tauy_shelf(G%isd:G%ied, G%JsdB:G%JedB), source=0.0)
2223	0	if (.not.allocated(visc%tbl_thick_shelf_u)) &
2224	0	allocate(visc%tbl_thick_shelf_u(G%IsdB:G%IedB, G%jsd:G%jed), source=0.0)
2225	0	if (.not.allocated(visc%tbl_thick_shelf_v)) &
2226	0	allocate(visc%tbl_thick_shelf_v(G%isd:G%ied, G%JsdB:G%JedB), source=0.0)
2227	0	if (.not.allocated(visc%kv_tbl_shelf_u)) &
2228	0	allocate(visc%kv_tbl_shelf_u(G%IsdB:G%IedB, G%jsd:G%jed), source=0.0)
2229	0	if (.not.allocated(visc%kv_tbl_shelf_v)) &
2230	0	allocate(visc%kv_tbl_shelf_v(G%isd:G%ied, G%JsdB:G%JedB), source=0.0)
2231
2232		! With a linear drag law under shelves, the friction velocity is already known.
2233		! if (CS%linear_drag) ustar(:) = cdrag_sqrt_H*CS%drag_bg_vel
2234
2235		! Find the vertical distances across layers.
2236	0	call thickness_to_dz(h, tv, dz, G, GV, US, halo_size=1)
2237		endif
2238
2239		!$OMP parallel do default(shared)
2240	0	do J=js-1,je ; do i=is-1,ie+1
2241	0	mask_v(i,J) = G%mask2dCv(i,J)
2242		enddo ; enddo
2243		!$OMP parallel do default(shared)
2244	0	do j=js-1,je+1 ; do I=is-1,ie
2245	0	mask_u(I,j) = G%mask2dCu(I,j)
2246		enddo ; enddo
2247
2248	0	if (associated(OBC)) then ; do n=1,OBC%number_of_segments
2249		! Project bottom depths across cell-corner points in the OBCs.
2250	0	if (.not. OBC%segment(n)%on_pe) cycle
2251		! Use a one-sided projection of bottom depths at OBC points.
2252	0	I = OBC%segment(n)%HI%IsdB ; J = OBC%segment(n)%HI%JsdB
2253	0	if (OBC%segment(n)%is_N_or_S .and. (J >= js-1) .and. (J <= je)) then
2254	0	do I = max(is-1,OBC%segment(n)%HI%IsdB), min(ie,OBC%segment(n)%HI%IedB)
2255	0	if (OBC%segment(n)%direction == OBC_DIRECTION_N) mask_u(I,j+1) = 0.0
2256	0	if (OBC%segment(n)%direction == OBC_DIRECTION_S) mask_u(I,j) = 0.0
2257		enddo
2258	0	elseif (OBC%segment(n)%is_E_or_W .and. (I >= is-1) .and. (I <= ie)) then
2259	0	do J = max(js-1,OBC%segment(n)%HI%JsdB), min(je,OBC%segment(n)%HI%JedB)
2260	0	if (OBC%segment(n)%direction == OBC_DIRECTION_E) mask_v(i+1,J) = 0.0
2261	0	if (OBC%segment(n)%direction == OBC_DIRECTION_W) mask_v(i,J) = 0.0
2262		enddo
2263		endif
2264		enddo ; endif
2265
2266		!$OMP parallel do default(private) shared(u,v,h,dz,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, &
2267		!$OMP nonBous_ML,h_neglect,dz_neglect,h_tiny,g_H_Rho0, &
2268		!$OMP js,je,OBC,Isq,Ieq,nz,nkml,U_star_2d,mask_v, &
2269		!$OMP cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL,Rho0x400_G)
2270	0	do j=js,je ! u-point loop
2271	0	if (CS%dynamic_viscous_ML) then
2272	0	do_any = .false.
2273	0	do I=Isq,Ieq
2274	0	htot(I) = 0.0
2275	0	if (G%mask2dCu(I,j) < 0.5) then
2276	0	do_i(I) = .false. ; visc%nkml_visc_u(I,j) = nkml
2277		else
2278	0	do_i(I) = .true. ; do_any = .true.
2279	0	k_massive(I) = nkml
2280	0	Thtot(I) = 0.0 ; Shtot(I) = 0.0 ; Rhtot(i) = 0.0
2281	0	uhtot(I) = dt_Rho0 * forces%taux(I,j)
2282		vhtot(I) = 0.25 * dt_Rho0 * ((forces%tauy(i,J) + forces%tauy(i+1,J-1)) + &
2283	0	(forces%tauy(i,J-1) + forces%tauy(i+1,J)))
2284
2285	0	if (CS%omega_frac >= 1.0) then ; absf = 2.0*CS%omega ; else
2286	0	absf = 0.5*(abs(G%CoriolisBu(I,J)) + abs(G%CoriolisBu(I,J-1)))
2287	0	if (CS%omega_frac > 0.0) &
2288	0	absf = sqrt(CS%omega_frac4.0CS%omega*2 + (1.0-CS%omega_frac)absf**2)
2289		endif
2290	0	U_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i+1,j)))
2291	0	Idecay_len_TKE(I) = (absf / U_star) * CS%TKE_decay
2292		endif
2293		enddo
2294
2295	0	if (do_any) then ; do k=1,nz
2296	0	if (k > nkml) then
2297	0	do_any = .false.
2298	0	if (use_EOS .and. (k==nkml+1)) then
2299		! Find dRho/dT and dRho_dS.
2300	0	do I=Isq,Ieq
2301	0	press(I) = (GV%H_to_RZGV%g_Earth) htot(I)
2302	0	if (associated(tv%p_surf)) press(I) = press(I) + 0.5*(tv%p_surf(i,j)+tv%p_surf(i+1,j))
2303	0	k2 = max(1,nkml)
2304	0	I_2hlay = 1.0 / (h(i,j,k2) + h(i+1,j,k2) + h_neglect)
2305	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
2306	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
2307		enddo
2308		call calculate_density_derivs(T_EOS, S_EOS, press, dR_dT, dR_dS, tv%eqn_of_state, &
2309	0	(/Isq-G%IsdB+1,Ieq-G%IsdB+1/) )
2310	0	if (nonBous_ML) then
2311		call calculate_specific_vol_derivs(T_EOS, S_EOS, press, dSpV_dT, dSpV_dS, tv%eqn_of_state, &
2312	0	(/Isq-G%IsdB+1,Ieq-G%IsdB+1/) )
2313		endif
2314		endif
2315
2316	0	do I=Isq,Ieq ; if (do_i(I)) then
2317
2318	0	hlay = 0.5*(h(i,j,k) + h(i+1,j,k))
2319	0	if (hlay > h_tiny) then ! Only consider non-vanished layers.
2320	0	I_2hlay = 1.0 / (h(i,j,k) + h(i+1,j,k))
2321		v_at_u = 0.5 * ((h(i,j,k) * (v(i,J,k) + v(i,J-1,k))) + &
2322	0	(h(i+1,j,k) * (v(i+1,J,k) + v(i+1,J-1,k)))) * I_2hlay
2323	0	Uh2 = (uhtot(I) - htot(I)u(I,j,k))2 + (vhtot(I) - htot(I)v_at_u)**2
2324
2325	0	if (use_EOS) then
2326	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
2327	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
2328	0	if (nonBous_ML) then
2329		gHprime = (GV%g_Earth * GV%H_to_RZ) * (dSpV_dT(I) * (Thtot(I) - T_lay*htot(I)) + &
2330	0	dSpV_dS(I) * (Shtot(I) - S_lay*htot(I)))
2331		else
2332		gHprime = g_H_Rho0 * (dR_dT(I) * (T_lay*htot(I) - Thtot(I)) + &
2333	0	dR_dS(I) * (S_lay*htot(I) - Shtot(I)))
2334		endif
2335		else
2336	0	gHprime = g_H_Rho0 * (GV%Rlay(k)*htot(I) - Rhtot(I))
2337		endif
2338
2339	0	if (gHprime > 0.0) then
2340	0	RiBulk = CS%bulk_Ri_ML * exp(-htot(I) * Idecay_len_TKE(I))
2341	0	if (RiBulk * Uh2 <= (htot(I)*2) gHprime) then
2342	0	visc%nkml_visc_u(I,j) = real(k_massive(I))
2343	0	do_i(I) = .false.
2344	0	elseif (RiBulk * Uh2 <= (htot(I) + hlay)*2 gHprime) then
2345		visc%nkml_visc_u(I,j) = real(k-1) + &
2346	0	( sqrt(RiBulk * Uh2 / gHprime) - htot(I) ) / hlay
2347	0	do_i(I) = .false.
2348		endif
2349		endif
2350	0	k_massive(I) = k
2351		endif ! hlay > h_tiny
2352
2353	0	if (do_i(I)) do_any = .true.
2354		endif ; enddo
2355
2356	0	if (.not.do_any) exit ! All columns are done.
2357		endif
2358
2359	0	do I=Isq,Ieq ; if (do_i(I)) then
2360	0	htot(I) = htot(I) + 0.5 * (h(i,j,k) + h(i+1,j,k))
2361	0	uhtot(I) = uhtot(I) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * u(I,j,k)
2362		vhtot(I) = vhtot(I) + 0.25 * ((h(i,j,k) * (v(i,J,k) + v(i,J-1,k))) + &
2363	0	(h(i+1,j,k) * (v(i+1,J,k) + v(i+1,J-1,k))))
2364	0	if (use_EOS) then
2365	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)))
2366	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)))
2367		else
2368	0	Rhtot(i) = Rhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * GV%Rlay(k)
2369		endif
2370		endif ; enddo
2371		enddo ; endif
2372
2373	0	if (do_any) then ; do I=Isq,Ieq ; if (do_i(I)) then
2374	0	visc%nkml_visc_u(I,j) = k_massive(I)
2375		endif ; enddo ; endif
2376		endif ! dynamic_viscous_ML
2377
2378	0	do_any_shelf = .false.
2379	0	if (associated(forces%frac_shelf_u)) then
2380	0	do I=Isq,Ieq
2381	0	if (forces%frac_shelf_u(I,j)*G%mask2dCu(I,j) == 0.0) then
2382	0	do_i(I) = .false.
2383	0	visc%tbl_thick_shelf_u(I,j) = 0.0 ; visc%kv_tbl_shelf_u(I,j) = 0.0
2384		else
2385	0	do_i(I) = .true. ; do_any_shelf = .true.
2386		endif
2387		enddo
2388		endif
2389
2390	0	if (do_any_shelf) then
2391	0	do k=1,nz ; do I=Isq,Ieq ; if (do_i(I)) then
2392	0	if (u(I,j,k) * (h(i+1,j,k) - h(i,j,k)) >= 0) then
2393		h_at_vel(i,k) = 2.0h(i,j,k)h(i+1,j,k) / &
2394	0	(h(i,j,k) + h(i+1,j,k) + h_neglect)
2395		dz_at_vel(i,k) = 2.0dz(i,j,k)dz(i+1,j,k) / &
2396	0	(dz(i,j,k) + dz(i+1,j,k) + dz_neglect)
2397		else
2398	0	h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
2399	0	dz_at_vel(i,k) = 0.5 * (dz(i,j,k) + dz(i+1,j,k))
2400		endif
2401		else
2402	0	h_at_vel(I,k) = 0.0
2403	0	dz_at_vel(I,k) = 0.0
2404	0	ustar(I) = 0.0
2405		endif ; enddo ; enddo
2406
2407	0	do I=Isq,Ieq ; if (do_i(I)) then
2408	0	htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
2409	0	Thtot(I) = 0.0 ; Shtot(I) = 0.0 ; SpV_htot(I) = 0.0
2410	0	if (use_EOS .or. .not.CS%linear_drag) then ; do k=1,nz
2411	0	if (htot_vel>=CS%Htbl_shelf) exit ! terminate the k loop
2412	0	hweight = MIN(CS%Htbl_shelf - htot_vel, h_at_vel(i,k))
2413	0	if (hweight <= 1.5*GV%Angstrom_H + h_neglect) cycle
2414
2415	0	htot_vel = htot_vel + h_at_vel(i,k)
2416	0	hwtot = hwtot + hweight
2417
2418	0	if (.not.CS%linear_drag) then
2419	0	v_at_u = set_v_at_u(v, h, G, GV, i, j, k, mask_v, OBC)
2420		! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant
2421	0	if (CS%BBL_use_tidal_bg) then
2422		u2_bg(I) = 0.5( G%mask2dT(i,j)(CS%tideamp(i,j)*CS%tideamp(i,j))+ &
2423	0	G%mask2dT(i+1,j)(CS%tideamp(i+1,j)CS%tideamp(i+1,j)) )
2424		else
2425	0	u2_bg(I) = CS%drag_bg_vel * CS%drag_bg_vel
2426		endif
2427	0	hutot = hutot + hweight * sqrt(u(I,j,k)2 + v_at_u2 + u2_bg(I))
2428		endif
2429	0	if (use_EOS) then
2430	0	Thtot(I) = Thtot(I) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
2431	0	Shtot(I) = Shtot(I) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
2432		endif
2433	0	if (allocated(tv%SpV_avg)) then
2434	0	SpV_htot(I) = SpV_htot(I) + hweight * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i+1,j,k))
2435		endif
2436		enddo ; endif
2437
2438	0	if ((hwtot <= 0.0) .or. (CS%linear_drag .and. .not.allocated(tv%SpV_avg))) then
2439	0	ustar(I) = cdrag_sqrt_H * CS%drag_bg_vel
2440	0	elseif (CS%linear_drag .and. allocated(tv%SpV_avg)) then
2441	0	ustar(I) = cdrag_sqrt_H_RL * CS%drag_bg_vel * (hwtot / SpV_htot(I))
2442	0	elseif (allocated(tv%SpV_avg)) then ! (.not.CS%linear_drag)
2443	0	ustar(I) = cdrag_sqrt_H_RL * hutot / SpV_htot(I)
2444		else ! (.not.CS%linear_drag .and. .not.allocated(tv%SpV_avg))
2445	0	ustar(I) = cdrag_sqrt_H * hutot / hwtot
2446		endif
2447
2448	0	if (use_EOS) then ; if (hwtot > 0.0) then
2449	0	T_EOS(I) = Thtot(I)/hwtot ; S_EOS(I) = Shtot(I)/hwtot
2450		else
2451	0	T_EOS(I) = 0.0 ; S_EOS(I) = 0.0
2452		endif ; endif
2453		! if (allocated(tv%SpV_avg)) SpV_av(I) = SpVhtot(I) / hwtot
2454		endif ; enddo ! I-loop
2455
2456	0	if (use_EOS) then
2457		call calculate_density_derivs(T_EOS, S_EOS, forces%p_surf(:,j), dR_dT, dR_dS, &
2458	0	tv%eqn_of_state, (/Isq-G%IsdB+1,Ieq-G%IsdB+1/) )
2459		endif
2460
2461	0	do I=Isq,Ieq ; if (do_i(I)) then
2462		! The 400.0 in this expression is the square of a constant proposed
2463		! by Killworth and Edwards, 1999, in equation (2.20).
2464	0	ustarsq = Rho0x400_G * ustar(i)**2
2465	0	htot(i) = 0.0 ; dztot(i) = 0.0
2466	0	if (use_EOS) then
2467	0	Thtot(i) = 0.0 ; Shtot(i) = 0.0 ; oldfn = 0.0
2468	0	do k=1,nz-1
2469	0	if (h_at_vel(i,k) <= 0.0) cycle
2470	0	T_Lay = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
2471	0	S_Lay = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
2472	0	oldfn = dR_dT(i)(T_Layhtot(i) - Thtot(i)) + dR_dS(i)(S_Layhtot(i) - Shtot(i))
2473	0	if (oldfn >= ustarsq) exit
2474
2475		Dfn = (dR_dT(i)(0.5(tv%T(i,j,k+1)+tv%T(i+1,j,k+1)) - T_Lay) + &
2476		dR_dS(i)(0.5(tv%S(i,j,k+1)+tv%S(i+1,j,k+1)) - S_Lay)) * &
2477	0	(h_at_vel(i,k)+htot(i))
2478	0	if ((oldfn + Dfn) <= ustarsq) then
2479	0	Dh = h_at_vel(i,k)
2480	0	Ddz = dz_at_vel(i,k)
2481		else
2482	0	frac_used = sqrt((ustarsq-oldfn) / (Dfn))
2483	0	Dh = h_at_vel(i,k) * frac_used
2484	0	Ddz = dz_at_vel(i,k) * frac_used
2485		endif
2486
2487	0	htot(i) = htot(i) + Dh
2488	0	dztot(i) = dztot(i) + Ddz
2489	0	Thtot(i) = Thtot(i) + T_LayDh ; Shtot(i) = Shtot(i) + S_LayDh
2490		enddo
2491	0	if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0)) then
2492	0	T_Lay = 0.5*(tv%T(i,j,nz) + tv%T(i+1,j,nz))
2493	0	S_Lay = 0.5*(tv%S(i,j,nz) + tv%S(i+1,j,nz))
2494	0	if (dR_dT(i)(T_Layhtot(i) - Thtot(i)) + &
2495		dR_dS(i)(S_Layhtot(i) - Shtot(i)) < ustarsq) then
2496	0	htot(i) = htot(i) + h_at_vel(i,nz)
2497	0	dztot(i) = dztot(i) + dz_at_vel(i,nz)
2498		endif
2499		endif ! Examination of layer nz.
2500		else ! Use Rlay as the density variable.
2501	0	Rhtot = 0.0
2502	0	do k=1,nz-1
2503	0	Rlay = GV%Rlay(k) ; Rlb = GV%Rlay(k+1)
2504
2505	0	oldfn = Rlay*htot(i) - Rhtot(i)
2506	0	if (oldfn >= ustarsq) exit
2507
2508	0	Dfn = (Rlb - Rlay)*(h_at_vel(i,k)+htot(i))
2509	0	if ((oldfn + Dfn) <= ustarsq) then
2510	0	Dh = h_at_vel(i,k)
2511	0	Ddz = dz_at_vel(i,k)
2512		else
2513	0	frac_used = sqrt((ustarsq-oldfn) / (Dfn))
2514	0	Dh = h_at_vel(i,k) * frac_used
2515	0	Ddz = dz_at_vel(i,k) * frac_used
2516		endif
2517
2518	0	htot(i) = htot(i) + Dh
2519	0	dztot(i) = dztot(i) + Ddz
2520	0	Rhtot(i) = Rhtot(i) + Rlay*Dh
2521		enddo
2522	0	if (GV%Rlay(nz)*htot(i) - Rhtot(i) < ustarsq) then
2523	0	htot(i) = htot(i) + h_at_vel(i,nz)
2524	0	dztot(i) = dztot(i) + dz_at_vel(i,nz)
2525		endif
2526		endif ! use_EOS
2527
2528		! visc%tbl_thick_shelf_u(I,j) = max(CS%Htbl_shelf_min, &
2529		! dztot(I) / (0.5 + sqrt(0.25 + &
2530		! ((htot(i)(G%CoriolisBu(I,J-1)+G%CoriolisBu(I,J)))*2) / &
2531		! (ustar(i)**2) )) )
2532	0	ustar1 = ustar(i)
2533	0	h2f2 = (htot(i)(G%CoriolisBu(I,J-1)+G%CoriolisBu(I,J)) + h_neglectCS%omega)**2
2534		tbl_thick = max(CS%Htbl_shelf_min, &
2535	0	( dztot(I)ustar(i) ) / ( 0.5ustar1 + sqrt((0.5ustar1)*2 + h2f2 ) ) )
2536	0	visc%tbl_thick_shelf_u(I,j) = tbl_thick
2537	0	visc%Kv_tbl_shelf_u(I,j) = max(CS%Kv_TBL_min, cdrag_sqrtustar1tbl_thick)
2538		endif ; enddo ! I-loop
2539		endif ! do_any_shelf
2540
2541		enddo ! j-loop at u-points
2542
2543		!$OMP parallel do default(private) shared(u,v,h,dz,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, &
2544		!$OMP nonBous_ML,h_neglect,dz_neglect,h_tiny,g_H_Rho0, &
2545		!$OMP is,ie,OBC,Jsq,Jeq,nz,nkml,U_star_2d,mask_u, &
2546		!$OMP cdrag_sqrt,cdrag_sqrt_H,cdrag_sqrt_H_RL,Rho0x400_G)
2547	0	do J=Jsq,Jeq ! v-point loop
2548	0	if (CS%dynamic_viscous_ML) then
2549	0	do_any = .false.
2550	0	do i=is,ie
2551	0	htot(i) = 0.0
2552	0	if (G%mask2dCv(i,J) < 0.5) then
2553	0	do_i(i) = .false. ; visc%nkml_visc_v(i,J) = nkml
2554		else
2555	0	do_i(i) = .true. ; do_any = .true.
2556	0	k_massive(i) = nkml
2557	0	Thtot(i) = 0.0 ; Shtot(i) = 0.0 ; Rhtot(i) = 0.0
2558	0	vhtot(i) = dt_Rho0 * forces%tauy(i,J)
2559		uhtot(i) = 0.25 * dt_Rho0 * ((forces%taux(I,j) + forces%taux(I-1,j+1)) + &
2560	0	(forces%taux(I-1,j) + forces%taux(I,j+1)))
2561
2562	0	if (CS%omega_frac >= 1.0) then ; absf = 2.0*CS%omega ; else
2563	0	absf = 0.5*(abs(G%CoriolisBu(I-1,J)) + abs(G%CoriolisBu(I,J)))
2564	0	if (CS%omega_frac > 0.0) &
2565	0	absf = sqrt(CS%omega_frac4.0CS%omega*2 + (1.0-CS%omega_frac)absf**2)
2566		endif
2567
2568	0	U_star = max(CS%ustar_min, 0.5*(U_star_2d(i,j) + U_star_2d(i,j+1)))
2569	0	Idecay_len_TKE(i) = (absf / U_star) * CS%TKE_decay
2570
2571		endif
2572		enddo
2573
2574	0	if (do_any) then ; do k=1,nz
2575	0	if (k > nkml) then
2576	0	do_any = .false.
2577	0	if (use_EOS .and. (k==nkml+1)) then
2578		! Find dRho/dT and dRho_dS.
2579	0	do i=is,ie
2580	0	press(i) = (GV%H_to_RZ * GV%g_Earth) * htot(i)
2581	0	if (associated(tv%p_surf)) press(i) = press(i) + 0.5*(tv%p_surf(i,j)+tv%p_surf(i,j+1))
2582	0	k2 = max(1,nkml)
2583	0	I_2hlay = 1.0 / (h(i,j,k2) + h(i,j+1,k2) + h_neglect)
2584	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
2585	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
2586		enddo
2587		call calculate_density_derivs(T_EOS, S_EOS, press, dR_dT, dR_dS, &
2588	0	tv%eqn_of_state, (/is-G%IsdB+1,ie-G%IsdB+1/) )
2589	0	if (nonBous_ML) then
2590		call calculate_specific_vol_derivs(T_EOS, S_EOS, press, dSpV_dT, dSpV_dS, tv%eqn_of_state, &
2591	0	(/is-G%IsdB+1,ie-G%IsdB+1/) )
2592		endif
2593		endif
2594
2595	0	do i=is,ie ; if (do_i(i)) then
2596
2597	0	hlay = 0.5*(h(i,j,k) + h(i,j+1,k))
2598	0	if (hlay > h_tiny) then ! Only consider non-vanished layers.
2599	0	I_2hlay = 1.0 / (h(i,j,k) + h(i,j+1,k))
2600		u_at_v = 0.5 * ((h(i,j,k) * (u(I-1,j,k) + u(I,j,k))) + &
2601	0	(h(i,j+1,k) * (u(I-1,j+1,k) + u(I,j+1,k)))) * I_2hlay
2602	0	Uh2 = (vhtot(i) - htot(i)v(i,J,k))2 + (uhtot(i) - htot(i)u_at_v)**2
2603
2604	0	if (use_EOS) then
2605	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
2606	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
2607	0	if (nonBous_ML) then
2608		gHprime = (GV%g_Earth * GV%H_to_RZ) * (dSpV_dT(i) * (Thtot(i) - T_lay*htot(i)) + &
2609	0	dSpV_dS(i) * (Shtot(i) - S_lay*htot(i)))
2610		else
2611		gHprime = g_H_Rho0 * (dR_dT(i) * (T_lay*htot(i) - Thtot(i)) + &
2612	0	dR_dS(i) * (S_lay*htot(i) - Shtot(i)))
2613		endif
2614		else
2615	0	gHprime = g_H_Rho0 * (GV%Rlay(k)*htot(i) - Rhtot(i))
2616		endif
2617
2618	0	if (gHprime > 0.0) then
2619	0	RiBulk = CS%bulk_Ri_ML * exp(-htot(i) * Idecay_len_TKE(i))
2620	0	if (RiBulk * Uh2 <= htot(i)*2 gHprime) then
2621	0	visc%nkml_visc_v(i,J) = real(k_massive(i))
2622	0	do_i(i) = .false.
2623	0	elseif (RiBulk * Uh2 <= (htot(i) + hlay)*2 gHprime) then
2624		visc%nkml_visc_v(i,J) = real(k-1) + &
2625	0	( sqrt(RiBulk * Uh2 / gHprime) - htot(i) ) / hlay
2626	0	do_i(i) = .false.
2627		endif
2628		endif
2629	0	k_massive(i) = k
2630		endif ! hlay > h_tiny
2631
2632	0	if (do_i(i)) do_any = .true.
2633		endif ; enddo
2634
2635	0	if (.not.do_any) exit ! All columns are done.
2636		endif
2637
2638	0	do i=is,ie ; if (do_i(i)) then
2639	0	htot(i) = htot(i) + 0.5 * (h(i,J,k) + h(i,j+1,k))
2640	0	vhtot(i) = vhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * v(i,J,k)
2641		uhtot(i) = uhtot(i) + 0.25 * ((h(i,j,k) * (u(I-1,j,k) + u(I,j,k))) + &
2642	0	(h(i,j+1,k) * (u(I-1,j+1,k) + u(I,j+1,k))))
2643	0	if (use_EOS) then
2644	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)))
2645	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)))
2646		else
2647	0	Rhtot(i) = Rhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * GV%Rlay(k)
2648		endif
2649		endif ; enddo
2650		enddo ; endif
2651
2652	0	if (do_any) then ; do i=is,ie ; if (do_i(i)) then
2653	0	visc%nkml_visc_v(i,J) = k_massive(i)
2654		endif ; enddo ; endif
2655
2656		endif ! dynamic_viscous_ML
2657
2658	0	do_any_shelf = .false.
2659	0	if (associated(forces%frac_shelf_v)) then
2660	0	do i=is,ie
2661	0	if (forces%frac_shelf_v(i,J)*G%mask2dCv(i,J) == 0.0) then
2662	0	do_i(i) = .false.
2663	0	visc%tbl_thick_shelf_v(i,J) = 0.0 ; visc%kv_tbl_shelf_v(i,J) = 0.0
2664		else
2665	0	do_i(i) = .true. ; do_any_shelf = .true.
2666		endif
2667		enddo
2668		endif
2669
2670	0	if (do_any_shelf) then
2671	0	do k=1,nz ; do i=is,ie ; if (do_i(i)) then
2672	0	if (v(i,J,k) * (h(i,j+1,k) - h(i,j,k)) >= 0) then
2673		h_at_vel(i,k) = 2.0h(i,j,k)h(i,j+1,k) / &
2674	0	(h(i,j,k) + h(i,j+1,k) + h_neglect)
2675		dz_at_vel(i,k) = 2.0dz(i,j,k)dz(i,j+1,k) / &
2676	0	(dz(i,j,k) + dz(i,j+1,k) + dz_neglect)
2677		else
2678	0	h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
2679	0	dz_at_vel(i,k) = 0.5 * (dz(i,j,k) + dz(i,j+1,k))
2680		endif
2681		else
2682	0	h_at_vel(I,k) = 0.0
2683	0	dz_at_vel(I,k) = 0.0
2684	0	ustar(i) = 0.0
2685		endif ; enddo ; enddo
2686
2687	0	do i=is,ie ; if (do_i(i)) then
2688	0	htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
2689	0	Thtot(i) = 0.0 ; Shtot(i) = 0.0 ; SpV_htot(i) = 0.0
2690	0	if (use_EOS .or. .not.CS%linear_drag) then ; do k=1,nz
2691	0	if (htot_vel>=CS%Htbl_shelf) exit ! terminate the k loop
2692	0	hweight = MIN(CS%Htbl_shelf - htot_vel, h_at_vel(i,k))
2693	0	if (hweight <= 1.5*GV%Angstrom_H + h_neglect) cycle
2694
2695	0	htot_vel = htot_vel + h_at_vel(i,k)
2696	0	hwtot = hwtot + hweight
2697
2698	0	if (.not.CS%linear_drag) then
2699	0	u_at_v = set_u_at_v(u, h, G, GV, i, J, k, mask_u, OBC)
2700		! Set the "back ground" friction velocity scale to either the tidal amplitude or place-holder constant
2701	0	if (CS%BBL_use_tidal_bg) then
2702		u2_bg(i) = 0.5( G%mask2dT(i,j)(CS%tideamp(i,j)*CS%tideamp(i,j))+ &
2703	0	G%mask2dT(i,j+1)(CS%tideamp(i,j+1)CS%tideamp(i,j+1)) )
2704		else
2705	0	u2_bg(i) = CS%drag_bg_vel * CS%drag_bg_vel
2706		endif
2707	0	hutot = hutot + hweight * sqrt(v(i,J,k)2 + u_at_v2 + u2_bg(i))
2708		endif
2709	0	if (use_EOS) then
2710	0	Thtot(i) = Thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
2711	0	Shtot(i) = Shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
2712		endif
2713	0	if (allocated(tv%SpV_avg)) then
2714	0	SpV_htot(i) = SpV_htot(i) + hweight * 0.5 * (tv%SpV_avg(i,j,k) + tv%SpV_avg(i,j+1,k))
2715		endif
2716		enddo ; endif
2717
2718	0	if ((hwtot <= 0.0) .or. (CS%linear_drag .and. .not.allocated(tv%SpV_avg))) then
2719	0	ustar(i) = cdrag_sqrt_H * CS%drag_bg_vel
2720	0	elseif (CS%linear_drag .and. allocated(tv%SpV_avg)) then
2721	0	ustar(i) = cdrag_sqrt_H_RL * CS%drag_bg_vel * (hwtot / SpV_htot(i))
2722	0	elseif (allocated(tv%SpV_avg)) then ! (.not.CS%linear_drag)
2723	0	ustar(i) = cdrag_sqrt_H_RL * hutot / SpV_htot(i)
2724		else ! (.not.CS%linear_drag .and. .not.allocated(tv%SpV_avg))
2725	0	ustar(i) = cdrag_sqrt_H * hutot / hwtot
2726		endif
2727
2728	0	if (use_EOS) then ; if (hwtot > 0.0) then
2729	0	T_EOS(i) = Thtot(i)/hwtot ; S_EOS(i) = Shtot(i)/hwtot
2730		else
2731	0	T_EOS(i) = 0.0 ; S_EOS(i) = 0.0
2732		endif ; endif
2733		endif ; enddo ! I-loop
2734
2735	0	if (use_EOS) then
2736		call calculate_density_derivs(T_EOS, S_EOS, forces%p_surf(:,j), dR_dT, dR_dS, &
2737	0	tv%eqn_of_state, (/is-G%IsdB+1,ie-G%IsdB+1/) )
2738		endif
2739
2740	0	do i=is,ie ; if (do_i(i)) then
2741		! The 400.0 in this expression is the square of a constant proposed
2742		! by Killworth and Edwards, 1999, in equation (2.20).
2743	0	ustarsq = Rho0x400_G * ustar(i)**2
2744	0	htot(i) = 0.0
2745	0	dztot(i) = 0.0
2746	0	if (use_EOS) then
2747	0	Thtot(i) = 0.0 ; Shtot(i) = 0.0 ; oldfn = 0.0
2748	0	do k=1,nz-1
2749	0	if (h_at_vel(i,k) <= 0.0) cycle
2750	0	T_Lay = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
2751	0	S_Lay = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
2752	0	oldfn = dR_dT(i)(T_Layhtot(i) - Thtot(i)) + dR_dS(i)(S_Layhtot(i) - Shtot(i))
2753	0	if (oldfn >= ustarsq) exit
2754
2755		Dfn = (dR_dT(i)(0.5(tv%T(i,j,k+1)+tv%T(i,j+1,k+1)) - T_Lay) + &
2756		dR_dS(i)(0.5(tv%S(i,j,k+1)+tv%S(i,j+1,k+1)) - S_Lay)) * &
2757	0	(h_at_vel(i,k)+htot(i))
2758	0	if ((oldfn + Dfn) <= ustarsq) then
2759	0	Dh = h_at_vel(i,k)
2760	0	Ddz = dz_at_vel(i,k)
2761		else
2762	0	frac_used = sqrt((ustarsq-oldfn) / (Dfn))
2763	0	Dh = h_at_vel(i,k) * frac_used
2764	0	Ddz = dz_at_vel(i,k) * frac_used
2765		endif
2766
2767	0	htot(i) = htot(i) + Dh
2768	0	dztot(i) = dztot(i) + Ddz
2769	0	Thtot(i) = Thtot(i) + T_LayDh ; Shtot(i) = Shtot(i) + S_LayDh
2770		enddo
2771	0	if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0)) then
2772	0	T_Lay = 0.5*(tv%T(i,j,nz) + tv%T(i,j+1,nz))
2773	0	S_Lay = 0.5*(tv%S(i,j,nz) + tv%S(i,j+1,nz))
2774	0	if (dR_dT(i)(T_Layhtot(i) - Thtot(i)) + &
2775		dR_dS(i)(S_Layhtot(i) - Shtot(i)) < ustarsq) then
2776	0	htot(i) = htot(i) + h_at_vel(i,nz)
2777	0	dztot(i) = dztot(i) + dz_at_vel(i,nz)
2778		endif
2779		endif ! Examination of layer nz.
2780		else ! Use Rlay as the density variable.
2781	0	Rhtot = 0.0
2782	0	do k=1,nz-1
2783	0	Rlay = GV%Rlay(k) ; Rlb = GV%Rlay(k+1)
2784
2785	0	oldfn = Rlay*htot(i) - Rhtot(i)
2786	0	if (oldfn >= ustarsq) exit
2787
2788	0	Dfn = (Rlb - Rlay)*(h_at_vel(i,k)+htot(i))
2789	0	if ((oldfn + Dfn) <= ustarsq) then
2790	0	Dh = h_at_vel(i,k)
2791	0	Ddz = dz_at_vel(i,k)
2792		else
2793	0	frac_used = sqrt((ustarsq-oldfn) / (Dfn))
2794	0	Dh = h_at_vel(i,k) * frac_used
2795	0	Ddz = dz_at_vel(i,k) * frac_used
2796		endif
2797
2798	0	htot(i) = htot(i) + Dh
2799	0	dztot(i) = dztot(i) + Ddz
2800	0	Rhtot = Rhtot + Rlay*Dh
2801		enddo
2802	0	if (GV%Rlay(nz)*htot(i) - Rhtot(i) < ustarsq) then
2803	0	htot(i) = htot(i) + h_at_vel(i,nz)
2804	0	dztot(i) = dztot(i) + dz_at_vel(i,nz)
2805		endif
2806		endif ! use_EOS
2807
2808		! visc%tbl_thick_shelf_v(i,J) = max(CS%Htbl_shelf_min, &
2809		! dztot(i) / (0.5 + sqrt(0.25 + &
2810		! (htot(i)(G%CoriolisBu(I-1,J)+G%CoriolisBu(I,J)))*2 / &
2811		! (ustar(i))**2 )) )
2812	0	ustar1 = ustar(i)
2813	0	h2f2 = (htot(i)(G%CoriolisBu(I-1,J)+G%CoriolisBu(I,J)) + h_neglectCS%omega)**2
2814		tbl_thick = max(CS%Htbl_shelf_min, &
2815	0	( dztot(i)ustar(i) ) / ( 0.5ustar1 + sqrt((0.5ustar1)*2 + h2f2 ) ) )
2816	0	visc%tbl_thick_shelf_v(i,J) = tbl_thick
2817	0	visc%Kv_tbl_shelf_v(i,J) = max(CS%Kv_TBL_min, cdrag_sqrtustar1tbl_thick)
2818
2819		endif ; enddo ! i-loop
2820		endif ! do_any_shelf
2821
2822		enddo ! J-loop at v-points
2823
2824		!$omp target update to(visc%nkml_visc_u, visc%nkml_visc_v) if (CS%dynamic_viscous_ML)
2825
2826	0	if (CS%debug) then
2827	0	if (allocated(visc%nkml_visc_u) .and. allocated(visc%nkml_visc_v)) &
2828		call uvchksum("nkml_visc_[uv]", visc%nkml_visc_u, visc%nkml_visc_v, &
2829	0	G%HI, haloshift=0, scalar_pair=.true.)
2830		endif
2831	0	if (CS%id_nkml_visc_u > 0) call post_data(CS%id_nkml_visc_u, visc%nkml_visc_u, CS%diag)
2832	0	if (CS%id_nkml_visc_v > 0) call post_data(CS%id_nkml_visc_v, visc%nkml_visc_v, CS%diag)
2833
2834		end subroutine set_viscous_ML
2835
2836		!> Register any fields associated with the vertvisc_type.
2837	1	subroutine set_visc_register_restarts(HI, G, GV, US, param_file, visc, restart_CS, use_ice_shelf)
2838		type(hor_index_type), intent(in) :: HI !< A horizontal index type structure.
2839		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2840		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2841		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2842		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
2843		!! parameters.
2844		type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical
2845		!! viscosities and related fields.
2846		!! Allocated here.
2847		type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure
2848		logical, intent(in) :: use_ice_shelf !< if true, register tau_shelf restarts
2849		! Local variables
2850		logical :: use_kappa_shear, KS_at_vertex
2851		logical :: adiabatic, useKPP, useEPBL, use_ideal_age
2852		logical :: do_brine_plume, use_hor_bnd_diff, use_neutral_diffusion, use_fpmix
2853		logical :: use_CVMix_shear, MLE_use_PBL_MLD, MLE_use_Bodner, use_CVMix_conv
2854		integer :: isd, ied, jsd, jed, nz
2855		real :: hfreeze !< If hfreeze > 0 [Z ~> m], melt potential will be computed.
2856		character(len=16) :: Kv_units, Kd_units
2857		character(len=40) :: mdl = "MOM_set_visc" ! This module's name.
2858	11	type(vardesc) :: u_desc, v_desc
2859	1	isd = HI%isd ; ied = HI%ied ; jsd = HI%jsd ; jed = HI%jed ; nz = GV%ke
2860
2861		call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
2862	1	do_not_log=.true.)
2863
2864	1	use_kappa_shear = .false. ; KS_at_vertex = .false. ; use_CVMix_shear = .false.
2865	1	useKPP = .false. ; useEPBL = .false. ; use_CVMix_conv = .false.
2866
2867	1	if (.not.adiabatic) then
2868	1	use_kappa_shear = kappa_shear_is_used(param_file)
2869	1	KS_at_vertex = kappa_shear_at_vertex(param_file)
2870	1	use_CVMix_shear = CVMix_shear_is_used(param_file)
2871	1	use_CVMix_conv = CVMix_conv_is_used(param_file)
2872		call get_param(param_file, mdl, "USE_KPP", useKPP, &
2873		"If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
2874		"to calculate diffusivities and non-local transport in the OBL.", &
2875	1	default=.false., do_not_log=.true.)
2876		call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", useEPBL, &
2877		"If true, use an implied energetics planetary boundary "//&
2878		"layer scheme to determine the diffusivity and viscosity "//&
2879	1	"in the surface boundary layer.", default=.false., do_not_log=.true.)
2880		endif
2881
2882	1	if (GV%Boussinesq) then
2883	1	Kv_units = "m2 s-1" ; Kd_units = "m2 s-1"
2884		else
2885	0	Kv_units = "Pa s" ; Kd_units = "kg m-1 s-1"
2886		endif
2887
2888	1	if (use_kappa_shear .or. useKPP .or. useEPBL .or. use_CVMix_shear .or. use_CVMix_conv) then
2889	1	call safe_alloc_ptr(visc%Kd_shear, isd, ied, jsd, jed, nz+1)
2890		call register_restart_field(visc%Kd_shear, "Kd_shear", .false., restart_CS, &
2891		"Shear-driven turbulent diffusivity at interfaces", &
2892	1	units=Kd_units, conversion=GV%HZ_T_to_MKS, z_grid='i')
2893		endif
2894	1	if (useKPP .or. useEPBL .or. use_CVMix_shear .or. use_CVMix_conv .or. &
2895		(use_kappa_shear .and. .not.KS_at_vertex )) then
2896	1	call safe_alloc_ptr(visc%Kv_shear, isd, ied, jsd, jed, nz+1)
2897		call register_restart_field(visc%Kv_shear, "Kv_shear", .false., restart_CS, &
2898		"Shear-driven turbulent viscosity at interfaces", &
2899	1	units=Kv_units, conversion=GV%HZ_T_to_MKS, z_grid='i')
2900		endif
2901	1	if (use_kappa_shear .and. KS_at_vertex) then
2902	1	call safe_alloc_ptr(visc%TKE_turb, HI%IsdB, HI%IedB, HI%JsdB, HI%JedB, nz+1)
2903	1	call safe_alloc_ptr(visc%Kv_shear_Bu, HI%IsdB, HI%IedB, HI%JsdB, HI%JedB, nz+1)
2904		call register_restart_field(visc%Kv_shear_Bu, "Kv_shear_Bu", .false., restart_CS, &
2905		"Shear-driven turbulent viscosity at vertex interfaces", &
2906	1	units=Kv_units, conversion=GV%HZ_T_to_MKS, hor_grid="Bu", z_grid='i')
2907	0	elseif (use_kappa_shear) then
2908	0	call safe_alloc_ptr(visc%TKE_turb, isd, ied, jsd, jed, nz+1)
2909		endif
2910
2911	1	if (useKPP) then
2912		! MOM_bkgnd_mixing uses Kv_slow when KPP is defined.
2913	0	call safe_alloc_ptr(visc%Kv_slow, isd, ied, jsd, jed, nz+1)
2914		endif
2915
2916		! visc%MLD and visc%h_ML are used to communicate the state of the (e)PBL or KPP to the rest of the model
2917		call get_param(param_file, mdl, "MLE_USE_PBL_MLD", MLE_use_PBL_MLD, &
2918	1	default=.false., do_not_log=.true.)
2919		! visc%h_ML needs to be allocated when melt potential is computed (HFREEZE>0) or one of
2920		! several other parameterizations are in use.
2921		call get_param(param_file, mdl, "HFREEZE", hfreeze, &
2922	1	units="m", default=-1.0, scale=US%m_to_Z, do_not_log=.true.)
2923		call get_param(param_file, mdl, "DO_BRINE_PLUME", do_brine_plume, &
2924		"If true, use a brine plume parameterization from Nguyen et al., 2009.", &
2925	1	default=.false., do_not_log=.true.)
2926		call get_param(param_file, mdl, "USE_HORIZONTAL_BOUNDARY_DIFFUSION", use_hor_bnd_diff, &
2927	1	default=.false., do_not_log=.true.)
2928		call get_param(param_file, mdl, "USE_NEUTRAL_DIFFUSION", use_neutral_diffusion, &
2929	1	default=.false., do_not_log=.true.)
2930	1	if (use_neutral_diffusion) &
2931		call get_param(param_file, mdl, "NDIFF_INTERIOR_ONLY", use_neutral_diffusion, &
2932	0	default=.false., do_not_log=.true.)
2933		call get_param(param_file, mdl, "FPMIX", use_fpmix, &
2934	1	default=.false., do_not_log=.true.)
2935		call get_param(param_file, mdl, "USE_IDEAL_AGE_TRACER", use_ideal_age, &
2936	1	default=.false., do_not_log=.true.)
2937	1	call openParameterBlock(param_file, 'MLE', do_not_log=.true.)
2938		call get_param(param_file, mdl, "USE_BODNER23", MLE_use_Bodner, &
2939	1	default=.false., do_not_log=.true.)
2940	1	call closeParameterBlock(param_file)
2941
2942	1	if (MLE_use_PBL_MLD .or. MLE_use_Bodner) then
2943	1	call safe_alloc_ptr(visc%MLD, isd, ied, jsd, jed)
2944		endif
2945		if ((hfreeze >= 0.0) .or. MLE_use_PBL_MLD .or. do_brine_plume .or. use_fpmix .or. &
2946	1	use_neutral_diffusion .or. use_hor_bnd_diff .or. use_ideal_age) then
2947	1	call safe_alloc_ptr(visc%h_ML, isd, ied, jsd, jed)
2948		endif
2949
2950	1	if (MLE_use_PBL_MLD) then
2951		call register_restart_field(visc%MLD, "MLD", .false., restart_CS, &
2952	1	"Instantaneous active mixing layer depth", units="m", conversion=US%Z_to_m)
2953		endif
2954		if (MLE_use_PBL_MLD .or. do_brine_plume .or. use_fpmix .or. &
2955	1	use_neutral_diffusion .or. use_hor_bnd_diff) then
2956		call register_restart_field(visc%h_ML, "h_ML", .false., restart_CS, &
2957		"Instantaneous active mixing layer thickness", &
2958	1	units=get_thickness_units(GV), conversion=GV%H_to_mks)
2959		endif
2960
2961		! visc%sfc_buoy_flx is used to communicate the state of the (e)PBL or KPP to the rest of the model
2962	1	if (MLE_use_PBL_MLD .or. MLE_use_Bodner) then
2963	1	call safe_alloc_ptr(visc%sfc_buoy_flx, isd, ied, jsd, jed)
2964		call register_restart_field(visc%sfc_buoy_flx, "SFC_BFLX", .false., restart_CS, &
2965		"Instantaneous surface buoyancy flux", "m2 s-3", &
2966	1	conversion=US%Z_to_m*2US%s_to_T**3)
2967		endif
2968
2969	1	if (use_ice_shelf) then
2970	0	if (.not.allocated(visc%taux_shelf)) &
2971	0	allocate(visc%taux_shelf(G%IsdB:G%IedB, G%jsd:G%jed), source=0.0)
2972	0	if (.not.allocated(visc%tauy_shelf)) &
2973	0	allocate(visc%tauy_shelf(G%isd:G%ied, G%JsdB:G%JedB), source=0.0)
2974		u_desc = var_desc("u_taux_shelf", "Pa", "the zonal stress on the ocean under ice shelves", &
2975	0	hor_grid='Cu',z_grid='1')
2976		v_desc = var_desc("v_tauy_shelf", "Pa", "the meridional stress on the ocean under ice shelves", &
2977	0	hor_grid='Cv',z_grid='1')
2978		call register_restart_pair(visc%taux_shelf, visc%tauy_shelf, u_desc, v_desc, &
2979	0	.false., restart_CS, conversion=US%RZ_T_to_kg_m2s*US%L_T_to_m_s)
2980		endif
2981
2982	11	end subroutine set_visc_register_restarts
2983
2984		!> This subroutine does remapping for the auxiliary restart variables in a vertvisc_type
2985		!! that are used across timesteps
2986	0	subroutine remap_vertvisc_aux_vars(G, GV, visc, h_old, h_new, ALE_CSp, OBC)
2987		type(ocean_grid_type), intent(inout) :: G !< ocean grid structure
2988		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2989		type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical
2990		!! viscosities and related fields.
2991		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2992		intent(in) :: h_old !< Thickness of source grid [H ~> m or kg m-2]
2993		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2994		intent(in) :: h_new !< Thickness of destination grid [H ~> m or kg m-2]
2995		type(ALE_CS), pointer :: ALE_CSp !< ALE control structure to use when remapping
2996		type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
2997
2998	0	if (associated(visc%Kd_shear)) then
2999	0	call ALE_remap_interface_vals(ALE_CSp, G, GV, h_old, h_new, visc%Kd_shear)
3000		endif
3001
3002	0	if (associated(visc%Kv_shear)) then
3003	0	call ALE_remap_interface_vals(ALE_CSp, G, GV, h_old, h_new, visc%Kv_shear)
3004		endif
3005
3006	0	if (associated(visc%Kv_shear_Bu)) then
3007	0	call ALE_remap_vertex_vals(ALE_CSp, G, GV, h_old, h_new, visc%Kv_shear_Bu)
3008		endif
3009
3010	0	end subroutine remap_vertvisc_aux_vars
3011
3012		!> Initializes the MOM_set_visc control structure
3013	1	subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS, OBC)
3014		type(time_type), target, intent(in) :: Time !< The current model time.
3015		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
3016		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
3017		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3018		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
3019		!! parameters.
3020		type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
3021		!! output.
3022		type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
3023		!! related fields.
3024		type(set_visc_CS), intent(inout) :: CS !< Vertical viscosity control structure
3025		type(MOM_restart_CS), intent(inout) :: restart_CS !< MOM restart control structure
3026		type(ocean_OBC_type), pointer :: OBC !< A pointer to an open boundary condition structure
3027
3028		! Local variables
3029		real :: Csmag_chan_dflt ! The default value for SMAG_CONST_CHANNEL [nondim]
3030		real :: smag_const1 ! The default value for the Smagorinsky Laplacian coefficient [nondim]
3031		real :: TKE_decay_dflt ! The default value of a coefficient scaling the vertical decay
3032		! rate of TKE [nondim]
3033		real :: bulk_Ri_ML_dflt ! The default bulk Richardson number for a bulk mixed layer [nondim]
3034		real :: Kv_background ! The background kinematic viscosity in the interior [Z2 T-1 ~> m2 s-1]
3035		real :: omega_frac_dflt ! The default value for the fraction of the absolute rotation rate that
3036		! is used in place of the absolute value of the local Coriolis
3037		! parameter in the denominator of some expressions [nondim]
3038		real :: Chan_max_thick_dflt ! The default value for CHANNEL_DRAG_MAX_THICK [Z ~> m]
3039		real :: tideamp_factor ! A factor to multiply by tideamp when converting to mean tidal magnitude [nondim]
3040		real :: shelfbreak_depth ! When CHANNEL_DRAG is true, the bathymetric depth interpolated
3041		! to the vorticity point is a combination of the harmonic mean of the
3042		! adjacent velocity point depths below this depth [Z ~> m] and the
3043		! arithmetic mean of the adjacent depths above it, to roughly mimic a
3044		! continental shelf break profile.
3045	1	real, allocatable, dimension(:,:) :: cdrag_h !< The spatially varying quadratic drag coefficient [nondim]
3046
3047		integer :: i, j, is, ie, js, je
3048		integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
3049		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
3050		logical :: adiabatic, use_omega, MLE_use_PBL_MLD
3051		logical :: use_KPP
3052		logical :: use_regridding ! If true, use the ALE algorithm rather than layered
3053		! isopycnal or stacked shallow water mode.
3054		logical :: use_temperature ! If true, temperature and salinity are used as state variables.
3055		logical :: use_EOS ! If true, density calculated from T & S using an equation of state.
3056		character(len=200) :: filename, cdrag_file, tideamp_file ! Input file names or paths
3057		character(len=80) :: cdrag_var, tideamp_var ! Input file variable names
3058		! This include declares and sets the variable "version".
3059		# include "version_variable.h"
3060		character(len=40) :: mdl = "MOM_set_visc" ! This module's name.
3061
3062	1	CS%initialized = .true.
3063	1	CS%OBC => OBC
3064
3065	1	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
3066	1	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed ; nz = GV%ke
3067	1	IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
3068
3069	1	CS%diag => diag
3070
3071		! Set default, read and log parameters
3072	1	call log_version(param_file, mdl, version, "")
3073	1	CS%RiNo_mix = .false.
3074	1	call get_param(param_file, mdl, "INPUTDIR", CS%inputdir, default=".")
3075	1	CS%inputdir = slasher(CS%inputdir)
3076		call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
3077		"This sets the default value for the various _ANSWER_DATE parameters.", &
3078	1	default=99991231)
3079		call get_param(param_file, mdl, "SET_VISC_ANSWER_DATE", CS%answer_date, &
3080		"The vintage of the order of arithmetic and expressions in the set viscosity "//&
3081		"calculations. Values below 20190101 recover the answers from the end of 2018, "//&
3082		"while higher values use updated and more robust forms of the same expressions.", &
3083	1	default=default_answer_date, do_not_log=.not.GV%Boussinesq)
3084	1	if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701)
3085		call get_param(param_file, mdl, "BOTTOMDRAGLAW", CS%bottomdraglaw, &
3086		"If true, the bottom stress is calculated with a drag "//&
3087		"law of the form c_drag\|u\|u. The velocity magnitude "//&
3088		"may be an assumed value or it may be based on the "//&
3089		"actual velocity in the bottommost HBBL, depending on "//&
3090	1	"LINEAR_DRAG.", default=.true.)
3091		call get_param(param_file, mdl, "DRAG_AS_BODY_FORCE", CS%body_force_drag, &
3092		"If true, the bottom stress is imposed as an explicit body force "//&
3093		"applied over a fixed distance from the bottom, rather than as an "//&
3094		"implicit calculation based on an enhanced near-bottom viscosity. "//&
3095		"The thickness of the bottom boundary layer is HBBL.", &
3096	1	default=.false., do_not_log=.not.CS%bottomdraglaw)
3097		call get_param(param_file, mdl, "CHANNEL_DRAG", CS%Channel_drag, &
3098		"If true, the bottom drag is exerted directly on each "//&
3099		"layer proportional to the fraction of the bottom it overlies.", &
3100	1	default=.false.)
3101		call get_param(param_file, mdl, "CHANNEL_DRAG_SHELFBREAK_DEPTH", shelfbreak_depth, &
3102		"When CHANNEL_DRAG is true, the bathymetric depth interpolated to the "//&
3103		"vorticity point is a combination of the harmonic mean of the adjacent "//&
3104		"velocity point depths below this depth and the arithmetic mean of the "//&
3105		"depths above it, to roughly mimic a continental shelf break profile. "//&
3106		"Setting this to exceed MAXIMUM_DEPTH leads to linear interpolation of "//&
3107		"the topography between velocity points.", &
3108	1	default=0.0, units="m", scale=US%m_to_Z, do_not_log=.not.CS%Channel_drag)
3109	1	CS%channel_break_depth = shelfbreak_depth - G%Z_ref
3110
3111		call get_param(param_file, mdl, "LINEAR_DRAG", CS%linear_drag, &
3112		"If LINEAR_DRAG and BOTTOMDRAGLAW are defined the drag "//&
3113	1	"law is cdragDRAG_BG_VELu.", default=.false.)
3114		call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
3115	1	do_not_log=.true.)
3116	1	if (adiabatic) then
3117		call log_param(param_file, mdl, "ADIABATIC",adiabatic, &
3118		"There are no diapycnal mass fluxes if ADIABATIC is true. "//&
3119		"This assumes that KD = 0.0 and that there is no buoyancy forcing, "//&
3120	0	"but makes the model faster by eliminating subroutine calls.", default=.false.)
3121		endif
3122
3123	1	if (.not.adiabatic) then
3124	1	CS%RiNo_mix = kappa_shear_is_used(param_file)
3125		endif
3126
3127	1	call get_param(param_file, mdl, "DEBUG", CS%debug, default=.false.)
3128
3129		call get_param(param_file, mdl, "DYNAMIC_VISCOUS_ML", CS%dynamic_viscous_ML, &
3130		"If true, use a bulk Richardson number criterion to "//&
3131		"determine the mixed layer thickness for viscosity.", &
3132	1	default=.false.)
3133	1	if (CS%dynamic_viscous_ML) then
3134	0	call get_param(param_file, mdl, "BULK_RI_ML", bulk_Ri_ML_dflt, units="nondim", default=0.0)
3135		call get_param(param_file, mdl, "BULK_RI_ML_VISC", CS%bulk_Ri_ML, &
3136		"The efficiency with which mean kinetic energy released by mechanically "//&
3137		"forced entrainment of the mixed layer is converted to turbulent "//&
3138		"kinetic energy. By default, BULK_RI_ML_VISC = BULK_RI_ML or 0.", &
3139	0	units="nondim", default=bulk_Ri_ML_dflt)
3140	0	call get_param(param_file, mdl, "TKE_DECAY", TKE_decay_dflt, units="nondim", default=0.0)
3141		call get_param(param_file, mdl, "TKE_DECAY_VISC", CS%TKE_decay, &
3142		"TKE_DECAY_VISC relates the vertical rate of decay of "//&
3143		"the TKE available for mechanical entrainment to the "//&
3144		"natural Ekman depth for use in calculating the dynamic "//&
3145		"mixed layer viscosity. By default, TKE_DECAY_VISC = TKE_DECAY or 0.", &
3146	0	units="nondim", default=TKE_decay_dflt)
3147		call get_param(param_file, mdl, "ML_USE_OMEGA", use_omega, &
3148		"If true, use the absolute rotation rate instead of the "//&
3149		"vertical component of rotation when setting the decay "//&
3150	0	"scale for turbulence.", default=.false., do_not_log=.true.)
3151	0	omega_frac_dflt = 0.0
3152	0	if (use_omega) then
3153	0	call MOM_error(WARNING, "ML_USE_OMEGA is deprecated; use ML_OMEGA_FRAC=1.0 instead.")
3154	0	omega_frac_dflt = 1.0
3155		endif
3156		call get_param(param_file, mdl, "ML_OMEGA_FRAC", CS%omega_frac, &
3157		"When setting the decay scale for turbulence, use this "//&
3158		"fraction of the absolute rotation rate blended with the "//&
3159		"local value of f, as sqrt((1-of)f^2 + of4*omega^2).", &
3160	0	units="nondim", default=omega_frac_dflt)
3161		call get_param(param_file, mdl, "OMEGA", CS%omega, &
3162		"The rotation rate of the earth.", &
3163	0	units="s-1", default=7.2921e-5, scale=US%T_to_s)
3164		! This give a minimum decay scale that is typically much less than Angstrom.
3165	0	CS%ustar_min = 2e-4CS%omega(GV%Angstrom_H + GV%H_subroundoff)
3166		else
3167		call get_param(param_file, mdl, "OMEGA", CS%omega, &
3168		"The rotation rate of the earth.", &
3169	1	units="s-1", default=7.2921e-5, scale=US%T_to_s)
3170		endif
3171
3172		call get_param(param_file, mdl, "HBBL", CS%dz_bbl, &
3173		"The thickness of a bottom boundary layer with a viscosity increased by "//&
3174		"KV_EXTRA_BBL if BOTTOMDRAGLAW is not defined, or the thickness over which "//&
3175		"near-bottom velocities are averaged for the drag law if BOTTOMDRAGLAW is "//&
3176		"defined but LINEAR_DRAG is not.", &
3177	1	units="m", scale=US%m_to_Z, fail_if_missing=.true.) ! Rescaled later
3178	1	if (CS%bottomdraglaw) then
3179		call get_param(param_file, mdl, "CDRAG", CS%cdrag, &
3180		"CDRAG is the drag coefficient relating the magnitude of "//&
3181		"the velocity field to the bottom stress. CDRAG is only "//&
3182	1	"used if BOTTOMDRAGLAW is defined.", units="nondim", default=0.003)
3183		call get_param(param_file, mdl, "CDRAG_MAP", CS%bottomdragmap, &
3184		"If true, apply a spatially varying scaling factor to CDRAG, "//&
3185	1	"specified by CDRAG_VAR in CDRAG_FILE.", default=.false.)
3186		call get_param(param_file, mdl, "CDRAG_FILE", cdrag_file, &
3187		"The name of the file with the spatially varying bottom drag "//&
3188	1	"scaling factor.", default="", do_not_log=.not.CS%bottomdragmap)
3189		call get_param(param_file, mdl, "CDRAG_VAR", cdrag_var, &
3190		"The name of the variable in CDRAG_FILE with the spatially "//&
3191		"varying bottom drag scaling factor at h points.", &
3192	1	default="", do_not_log=.not.CS%bottomdragmap)
3193		call get_param(param_file, mdl, "BBL_USE_TIDAL_BG", CS%BBL_use_tidal_bg, &
3194		"Flag to use the tidal RMS amplitude in place of constant "//&
3195		"background velocity for computing u* in the BBL. "//&
3196		"This flag is only used when BOTTOMDRAGLAW is true and "//&
3197	1	"LINEAR_DRAG is false.", default=.false.)
3198	1	if (CS%BBL_use_tidal_bg) then
3199		call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, &
3200		"The path to the file containing the spatially varying "//&
3201	0	"tidal amplitudes with INT_TIDE_DISSIPATION.", default="tideamp.nc")
3202		call get_param(param_file, mdl, "TIDEAMP_VARNAME", tideamp_var, &
3203		"The name of the tidal amplitude variable in the input file.", &
3204	0	default="tideamp")
3205		! This value is here only to detect whether it is inadvertently used. CS%drag_bg_vel should
3206		! not be used if CS%BBL_use_tidal_bg is True. For this reason, we do not apply dimensions,
3207		! nor dimensional testing in this mode. If we ever detect a dimensional sensitivity to
3208		! this parameter, in this mode, then it means it is being used inappropriately.
3209	0	CS%drag_bg_vel = 1.e30
3210		call get_param(param_file, mdl, "TIDEAMP_FACTOR", tideamp_factor, &
3211		"A parameter to multiply by tideamp when converting to ustar. "//&
3212		"It accounts for converting the amplitude to a mean magintude (approx 1/sqrt(2)) "//&
3213		"and possibly also for non-commuting averaging operators when converting to ustar**3. "//&
3214		"It is ignored if negative and uncapped so it can be greater than 1 if desired.",&
3215	0	units="nondim", default=-1.0)
3216	0	if (tideamp_factor < 0.0) then
3217	0	CS%tideampfac2 = 1.0
3218		else
3219	0	CS%tideampfac2 = tideamp_factor*tideamp_factor
3220		endif
3221		else
3222		call get_param(param_file, mdl, "DRAG_BG_VEL", CS%drag_bg_vel, &
3223		"DRAG_BG_VEL is either the assumed bottom velocity (with "//&
3224		"LINEAR_DRAG) or an unresolved velocity that is "//&
3225		"combined with the resolved velocity to estimate the "//&
3226		"velocity magnitude. DRAG_BG_VEL is only used when "//&
3227	1	"BOTTOMDRAGLAW is defined.", units="m s-1", default=0.0, scale=US%m_s_to_L_T)
3228		endif
3229		call get_param(param_file, mdl, "USE_REGRIDDING", use_regridding, &
3230	1	do_not_log=.true., default=.false. )
3231		call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", use_temperature, &
3232	1	default=.true., do_not_log=.true.)
3233		call get_param(param_file, mdl, "USE_EOS", use_EOS, &
3234	1	default=use_temperature, do_not_log=.true.)
3235		call get_param(param_file, mdl, "BBL_USE_EOS", CS%BBL_use_EOS, &
3236		"If true, use the equation of state in determining the properties of the "//&
3237		"bottom boundary layer. Otherwise use the layer target potential densities. "//&
3238		"The default of this parameter is the value of USE_EOS.", &
3239	1	default=use_EOS, do_not_log=.not.use_temperature)
3240	1	if (use_regridding .and. (.not. CS%BBL_use_EOS)) &
3241	0	call MOM_error(FATAL,"When using MOM6 in ALE mode it is required to set BBL_USE_EOS to True.")
3242		endif
3243		call get_param(param_file, mdl, "BBL_THICK_MIN", CS%BBL_thick_min, &
3244		"The minimum bottom boundary layer thickness that can be "//&
3245		"used with BOTTOMDRAGLAW. This might be "//&
3246		"Kv/(cdrag*drag_bg_vel) to give Kv as the minimum "//&
3247	1	"near-bottom viscosity.", units="m", default=0.0, scale=US%m_to_Z)
3248		call get_param(param_file, mdl, "HTBL_SHELF_MIN", CS%Htbl_shelf_min, &
3249		"The minimum top 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-top viscosity.", units="m", default=US%Z_to_m*CS%BBL_thick_min, scale=US%m_to_Z)
3253		call get_param(param_file, mdl, "HTBL_SHELF", CS%Htbl_shelf, &
3254		"The thickness over which near-surface velocities are "//&
3255		"averaged for the drag law under an ice shelf. By "//&
3256		"default this is the same as HBBL", &
3257	1	units="m", default=US%Z_to_m*CS%dz_bbl, scale=GV%m_to_H)
3258
3259		call get_param(param_file, mdl, "KV", Kv_background, &
3260		"The background kinematic viscosity in the interior. "//&
3261		"The molecular value, ~1e-6 m2 s-1, may be used.", &
3262	1	units="m2 s-1", scale=US%m2_s_to_Z2_T, fail_if_missing=.true.)
3263
3264		call get_param(param_file, mdl, "USE_KPP", use_KPP, &
3265		"If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
3266		"to calculate diffusivities and non-local transport in the OBL.", &
3267	1	do_not_log=.true., default=.false.)
3268
3269		call get_param(param_file, mdl, "KV_BBL_MIN", CS%KV_BBL_min, &
3270		"The minimum viscosities in the bottom boundary layer.", &
3271	1	units="m2 s-1", default=US%Z2_T_to_m2_s*Kv_background, scale=GV%m2_s_to_HZ_T)
3272		call get_param(param_file, mdl, "KV_TBL_MIN", CS%KV_TBL_min, &
3273		"The minimum viscosities in the top boundary layer.", &
3274	1	units="m2 s-1", default=US%Z2_T_to_m2_s*Kv_background, scale=GV%m2_s_to_HZ_T)
3275		call get_param(param_file, mdl, "CORRECT_BBL_BOUNDS", CS%correct_BBL_bounds, &
3276		"If true, uses the correct bounds on the BBL thickness and "//&
3277		"viscosity so that the bottom layer feels the intended drag.", &
3278	1	default=.false.)
3279
3280	1	if (CS%Channel_drag) then
3281	1	call get_param(param_file, mdl, "SMAG_LAP_CONST", smag_const1, units="nondim", default=-1.0)
3282
3283	1	cSmag_chan_dflt = 0.15
3284	1	if (smag_const1 >= 0.0) cSmag_chan_dflt = smag_const1
3285
3286		call get_param(param_file, mdl, "SMAG_CONST_CHANNEL", CS%c_Smag, &
3287		"The nondimensional Laplacian Smagorinsky constant used "//&
3288		"in calculating the channel drag if it is enabled. The "//&
3289		"default is to use the same value as SMAG_LAP_CONST if "//&
3290		"it is defined, or 0.15 if it is not. The value used is "//&
3291		"also 0.15 if the specified value is negative.", &
3292	1	units="nondim", default=cSmag_chan_dflt, do_not_log=.not.CS%Channel_drag)
3293	1	if (CS%c_Smag < 0.0) CS%c_Smag = 0.15
3294
3295		call get_param(param_file, mdl, "TRIG_CHANNEL_DRAG_WIDTHS", CS%concave_trigonometric_L, &
3296		"If true, use trigonometric expressions to determine the fractional open "//&
3297		"interface lengths for concave topography.", &
3298	1	default=.true., do_not_log=.not.CS%Channel_drag)
3299		endif
3300
3301	1	Chan_max_thick_dflt = -1.0*US%m_to_Z
3302	1	if (CS%RiNo_mix) Chan_max_thick_dflt = 0.5*CS%dz_bbl
3303	1	if (CS%body_force_drag) Chan_max_thick_dflt = CS%dz_bbl
3304		call get_param(param_file, mdl, "CHANNEL_DRAG_MAX_BBL_THICK", CS%Chan_drag_max_vol, &
3305		"The maximum bottom boundary layer thickness over which the channel drag is "//&
3306		"exerted, or a negative value for no fixed limit, instead basing the BBL "//&
3307		"thickness on the bottom stress, rotation and stratification. The default is "//&
3308		"proportional to HBBL if USE_JACKSON_PARAM or DRAG_AS_BODY_FORCE is true.", &
3309		units="m", default=US%Z_to_m*Chan_max_thick_dflt, scale=US%m_to_Z, &
3310	1	do_not_log=.not.CS%Channel_drag)
3311
3312		call get_param(param_file, mdl, "MLE_USE_PBL_MLD", MLE_use_PBL_MLD, &
3313	1	default=.false., do_not_log=.true.)
3314
3315	1	CS%Hbbl = CS%dz_bbl * (US%Z_to_m * GV%m_to_H) ! Rescaled for use in expressions in thickness units.
3316
3317		! Update scalar data to device
3318		!$omp target update to(visc)
3319		!$omp target update to(CS)
3320
3321		! Now update arrays that were defined in set_visc_register_restarts()
3322		!$omp target enter data map(to: visc%Kv_shear) &
3323		!$omp if (associated(visc%Kv_shear))
3324		!$omp target enter data map(to: visc%Kv_shear_Bu) &
3325		!$omp if (associated(visc%Kv_shear_Bu))
3326
3327	1	if (CS%RiNo_mix .and. kappa_shear_at_vertex(param_file)) then
3328		! This is necessary for reproducibility across restarts in non-symmetric mode.
3329	1	call pass_var(visc%Kv_shear_Bu, G%Domain, position=CORNER, complete=.true.)
3330		endif
3331
3332	1	if (CS%bottomdraglaw) then
3333	8841	allocate(visc%bbl_thick_u(IsdB:IedB,jsd:jed), source=0.0)
3334	8902	allocate(visc%bbl_thick_v(isd:ied,JsdB:JedB), source=0.0)
3335		!$omp target enter data map(to: visc%bbl_thick_u, visc%bbl_thick_v)
3336	8841	allocate(visc%kv_bbl_u(IsdB:IedB,jsd:jed), source=0.0)
3337	8902	allocate(visc%kv_bbl_v(isd:ied,JsdB:JedB), source=0.0)
3338		!$omp target enter data map(to: visc%kv_bbl_u, visc%kv_bbl_v)
3339	8773	allocate(visc%ustar_bbl(isd:ied,jsd:jed), source=0.0)
3340	8773	allocate(visc%BBL_meanKE_loss(isd:ied,jsd:jed), source=0.0)
3341	8773	allocate(visc%BBL_meanKE_loss_sqrtCd(isd:ied,jsd:jed), source=0.0)
3342
3343		CS%id_bbl_thick_u = register_diag_field('ocean_model', 'bbl_thick_u', &
3344	1	diag%axesCu1, Time, 'BBL thickness at u points', 'm', conversion=US%Z_to_m)
3345		CS%id_kv_bbl_u = register_diag_field('ocean_model', 'kv_bbl_u', diag%axesCu1, &
3346	1	Time, 'BBL viscosity at u points', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
3347		CS%id_bbl_u = register_diag_field('ocean_model', 'bbl_u', diag%axesCu1, &
3348	1	Time, 'BBL mean u current', 'm s-1', conversion=US%L_T_to_m_s)
3349	1	if (CS%id_bbl_u>0) then
3350	0	allocate(CS%bbl_u(IsdB:IedB,jsd:jed), source=0.0)
3351		endif
3352		CS%id_bbl_thick_v = register_diag_field('ocean_model', 'bbl_thick_v', &
3353	1	diag%axesCv1, Time, 'BBL thickness at v points', 'm', conversion=US%Z_to_m)
3354		CS%id_kv_bbl_v = register_diag_field('ocean_model', 'kv_bbl_v', diag%axesCv1, &
3355	1	Time, 'BBL viscosity at v points', 'm2 s-1', conversion=GV%HZ_T_to_m2_s)
3356		CS%id_bbl_v = register_diag_field('ocean_model', 'bbl_v', diag%axesCv1, &
3357	1	Time, 'BBL mean v current', 'm s-1', conversion=US%L_T_to_m_s)
3358	1	if (CS%id_bbl_v>0) then
3359	0	allocate(CS%bbl_v(isd:ied,JsdB:JedB), source=0.0)
3360		endif
3361	1	if (CS%bottomdragmap) then
3362	0	if (len_trim(cdrag_file)==0 .or. len_trim(cdrag_var)==0) then
3363	0	call MOM_error(FATAL,"CDRAG_FILE and CDRAG_VAR are required when using CDRAG_MAP.")
3364		endif
3365	0	allocate(cdrag_h(isd:ied,jsd:jed), source=0.0)
3366	0	allocate(CS%cdrag_u(IsdB:IedB,jsd:jed), source=0.0)
3367	0	allocate(CS%cdrag_v(isd:ied,JsdB:JedB), source=0.0)
3368	0	filename = trim(CS%inputdir) // trim(cdrag_file)
3369	0	call log_param(param_file, mdl, "INPUTDIR/CDRAG_FILE", filename)
3370	0	call MOM_read_data(filename, cdrag_var, cdrag_h, G%domain, scale=CS%cdrag)
3371	0	call pass_var(cdrag_h, G%domain)
3372	0	do j=js,je ; do I=is-1,ie ; if (G%mask2dCu(I,j) > 0) then
3373		CS%cdrag_u(I,j) = (G%mask2dT(i,j) * cdrag_h(i,j) + G%mask2dT(i+1,j) * cdrag_h(i+1,j)) / &
3374	0	(G%mask2dT(i,j) + G%mask2dT(i+1,j))
3375		endif ; enddo ; enddo
3376	0	do J=js-1,je ; do i=is,ie ; if (G%mask2dCv(i,J) > 0) then
3377		CS%cdrag_v(i,J) = (G%mask2dT(i,j) * cdrag_h(i,j) + G%mask2dT(i,j+1) * cdrag_h(i,j+1)) / &
3378	0	(G%mask2dT(i,j) + G%mask2dT(i,j+1))
3379		endif ; enddo ; enddo
3380	0	deallocate(cdrag_h)
3381		endif
3382	1	if (CS%BBL_use_tidal_bg) then
3383	0	allocate(CS%tideamp(isd:ied,jsd:jed), source=0.0)
3384	0	filename = trim(CS%inputdir) // trim(tideamp_file)
3385	0	call log_param(param_file, mdl, "INPUTDIR/TIDEAMP_FILE", filename)
3386	0	call MOM_read_data(filename, tideamp_var, CS%tideamp, G%domain, scale=US%m_to_Z*US%T_to_s)
3387	0	call pass_var(CS%tideamp,G%domain)
3388		endif
3389		endif
3390	1	if (CS%Channel_drag .or. CS%body_force_drag) then
3391	663076	allocate(visc%Ray_u(IsdB:IedB,jsd:jed,nz), source=0.0)
3392	667651	allocate(visc%Ray_v(isd:ied,JsdB:JedB,nz), source=0.0)
3393		!$omp target enter data map(to: visc%Ray_u, visc%Ray_v)
3394		CS%id_Ray_u = register_diag_field('ocean_model', 'Rayleigh_u', diag%axesCuL, &
3395	1	Time, 'Rayleigh drag velocity at u points', 'm s-1', conversion=GV%H_to_m*US%s_to_T)
3396		CS%id_Ray_v = register_diag_field('ocean_model', 'Rayleigh_v', diag%axesCvL, &
3397	1	Time, 'Rayleigh drag velocity at v points', 'm s-1', conversion=GV%H_to_m*US%s_to_T)
3398		endif
3399
3400
3401	1	if (CS%dynamic_viscous_ML) then
3402	0	allocate(visc%nkml_visc_u(IsdB:IedB,jsd:jed), source=0.0)
3403	0	allocate(visc%nkml_visc_v(isd:ied,JsdB:JedB), source=0.0)
3404		!$omp target enter data map(to: visc%nkml_visc_u, visc%nkml_visc_v)
3405
3406		CS%id_nkml_visc_u = register_diag_field('ocean_model', 'nkml_visc_u', &
3407	0	diag%axesCu1, Time, 'Number of layers in viscous mixed layer at u points', 'nondim')
3408		CS%id_nkml_visc_v = register_diag_field('ocean_model', 'nkml_visc_v', &
3409	0	diag%axesCv1, Time, 'Number of layers in viscous mixed layer at v points', 'nondim')
3410		endif
3411
3412	1	call register_restart_field_as_obsolete('Kd_turb','Kd_shear', restart_CS)
3413	1	call register_restart_field_as_obsolete('Kv_turb','Kv_shear', restart_CS)
3414
3415	1	end subroutine set_visc_init
3416
3417		!> This subroutine dellocates any memory in the set_visc control structure.
3418	1	subroutine set_visc_end(visc, CS)
3419		type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
3420		!! related fields. Elements are deallocated here.
3421		type(set_visc_CS), intent(inout) :: CS !< The control structure returned by a previous
3422		!! call to set_visc_init.
3423
3424	1	if (allocated(visc%bbl_thick_u)) then
3425		!$omp target exit data map(delete: visc%bbl_thick_u)
3426	1	deallocate(visc%bbl_thick_u)
3427		endif
3428	1	if (allocated(visc%bbl_thick_v)) then
3429		!$omp target exit data map(delete: visc%bbl_thick_v)
3430	1	deallocate(visc%bbl_thick_v)
3431		endif
3432	1	if (allocated(visc%kv_bbl_u)) then
3433		!$omp target exit data map(delete: visc%kv_bbl_u)
3434	1	deallocate(visc%kv_bbl_u)
3435		endif
3436	1	if (allocated(visc%kv_bbl_v)) then
3437		!$omp target exit data map(delete: visc%kv_bbl_v)
3438	1	deallocate(visc%kv_bbl_v)
3439		endif
3440	1	if (allocated(CS%bbl_u)) deallocate(CS%bbl_u)
3441	1	if (allocated(CS%bbl_v)) deallocate(CS%bbl_v)
3442	1	if (allocated(visc%Ray_u)) then
3443		!$omp target exit data map(delete: visc%Ray_u)
3444	1	deallocate(visc%Ray_u)
3445		endif
3446	1	if (allocated(visc%Ray_v)) then
3447		!$omp target exit data map(delete: visc%Ray_v)
3448	1	deallocate(visc%Ray_v)
3449		endif
3450	1	if (allocated(visc%nkml_visc_u)) then
3451		!$omp target exit data map(delete: visc%nkml_visc_u)
3452	0	deallocate(visc%nkml_visc_u)
3453		endif
3454	1	if (allocated(visc%nkml_visc_v)) then
3455		!$omp target exit data map(delete: visc%nkml_visc_v)
3456	0	deallocate(visc%nkml_visc_v)
3457		endif
3458	1	if (associated(visc%Kd_shear)) deallocate(visc%Kd_shear)
3459	1	if (associated(visc%Kv_slow)) deallocate(visc%Kv_slow)
3460	1	if (associated(visc%TKE_turb)) deallocate(visc%TKE_turb)
3461	1	if (associated(visc%Kv_shear)) deallocate(visc%Kv_shear)
3462	1	if (associated(visc%Kv_shear_Bu)) deallocate(visc%Kv_shear_Bu)
3463	1	if (allocated(visc%ustar_bbl)) deallocate(visc%ustar_bbl)
3464	1	if (allocated(visc%BBL_meanKE_loss)) deallocate(visc%BBL_meanKE_loss)
3465	1	if (allocated(visc%BBL_meanKE_loss_sqrtCd)) deallocate(visc%BBL_meanKE_loss_sqrtCd)
3466	1	if (allocated(visc%taux_shelf)) deallocate(visc%taux_shelf)
3467	1	if (allocated(visc%tauy_shelf)) deallocate(visc%tauy_shelf)
3468	1	if (allocated(visc%tbl_thick_shelf_u)) deallocate(visc%tbl_thick_shelf_u)
3469	1	if (allocated(visc%tbl_thick_shelf_v)) deallocate(visc%tbl_thick_shelf_v)
3470	1	if (allocated(visc%kv_tbl_shelf_u)) deallocate(visc%kv_tbl_shelf_u)
3471	1	if (allocated(visc%kv_tbl_shelf_v)) deallocate(visc%kv_tbl_shelf_v)
3472	1	end subroutine set_visc_end
3473
3474		!> \namespace mom_set_visc
3475		!!
3476		!! This would also be the module in which other viscous quantities that are flow-independent might be set.
3477		!! This information is transmitted to other modules via a vertvisc type structure.
3478		!!
3479		!! The same code is used for the two velocity components, by indirectly referencing the velocities and
3480		!! defining a handful of direction-specific defined variables.
3481
3482	0	end module MOM_set_visc