MOM6
MOM_set_viscosity.F90
Go to the documentation of this file.
1 !> Calculates various values related to the bottom boundary layer, such as the viscosity and
2 !! thickness of the BBL (set_viscous_BBL).
3 module mom_set_visc
4 
5 ! This file is part of MOM6. See LICENSE.md for the license.
6 
7 use mom_debugging, only : uvchksum, hchksum
8 use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, clock_routine
9 use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
10 use mom_diag_mediator, only : diag_ctrl, time_type
11 use mom_domains, only : pass_var, corner
12 use mom_error_handler, only : mom_error, fatal, warning
13 use mom_file_parser, only : get_param, log_param, log_version, param_file_type
14 use mom_forcing_type, only : forcing, mech_forcing
15 use mom_grid, only : ocean_grid_type
16 use mom_hor_index, only : hor_index_type
17 use mom_kappa_shear, only : kappa_shear_is_used, kappa_shear_at_vertex
18 use mom_cvmix_shear, only : cvmix_shear_is_used
19 use mom_cvmix_conv, only : cvmix_conv_is_used
20 use mom_cvmix_ddiff, only : cvmix_ddiff_is_used
21 use mom_restart, only : register_restart_field, query_initialized, mom_restart_cs
22 use mom_restart, only : register_restart_field_as_obsolete
23 use mom_safe_alloc, only : safe_alloc_ptr, safe_alloc_alloc
24 use mom_unit_scaling, only : unit_scale_type
25 use mom_variables, only : thermo_var_ptrs, vertvisc_type
26 use mom_verticalgrid, only : verticalgrid_type
27 use mom_eos, only : calculate_density, calculate_density_derivs
28 use mom_open_boundary, only : ocean_obc_type, obc_none, obc_direction_e
29 use mom_open_boundary, only : obc_direction_w, obc_direction_n, obc_direction_s
30 use mom_open_boundary, only : obc_segment_type
31 implicit none ; private
32 
33 #include <MOM_memory.h>
34 
35 public set_viscous_bbl, set_viscous_ml, set_visc_init, set_visc_end
36 public set_visc_register_restarts
37 
38 ! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
39 ! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
40 ! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
41 ! vary with the Boussinesq approximation, the Boussinesq variant is given first.
42 
43 !> Control structure for MOM_set_visc
44 type, public :: set_visc_cs ; private
45  real :: hbbl !< The static bottom boundary layer thickness [H ~> m or kg m-2]
46  real :: cdrag !< The quadratic drag coefficient.
47  real :: c_smag !< The Laplacian Smagorinsky coefficient for
48  !! calculating the drag in channels.
49  real :: drag_bg_vel !< An assumed unresolved background velocity for
50  !! calculating the bottom drag [L T-1 ~> m s-1].
51  real :: bbl_thick_min !< The minimum bottom boundary layer thickness [H ~> m or kg m-2].
52  !! This might be Kv / (cdrag * drag_bg_vel) to give
53  !! Kv as the minimum near-bottom viscosity.
54  real :: htbl_shelf !< A nominal thickness of the surface boundary layer for use
55  !! in calculating the near-surface velocity [H ~> m or kg m-2].
56  real :: htbl_shelf_min !< The minimum surface boundary layer thickness [H ~> m or kg m-2].
57  real :: kv_bbl_min !< The minimum viscosity in the bottom boundary layer [Z2 T-1 ~> m2 s-1].
58  real :: kv_tbl_min !< The minimum viscosity in the top boundary layer [Z2 T-1 ~> m2 s-1].
59  logical :: bottomdraglaw !< If true, the bottom stress is calculated with a
60  !! drag law c_drag*|u|*u. The velocity magnitude
61  !! may be an assumed value or it may be based on the
62  !! actual velocity in the bottommost HBBL, depending
63  !! on whether linear_drag is true.
64  logical :: bbl_use_eos !< If true, use the equation of state in determining
65  !! the properties of the bottom boundary layer.
66  logical :: linear_drag !< If true, the drag law is cdrag*DRAG_BG_VEL*u.
67  logical :: channel_drag !< If true, the drag is exerted directly on each
68  !! layer according to what fraction of the bottom
69  !! they overlie.
70  logical :: rino_mix !< If true, use Richardson number dependent mixing.
71  logical :: dynamic_viscous_ml !< If true, use a bulk Richardson number criterion to
72  !! determine the mixed layer thickness for viscosity.
73  real :: bulk_ri_ml !< The bulk mixed layer used to determine the
74  !! thickness of the viscous mixed layer. Nondim.
75  real :: omega !< The Earth's rotation rate [T-1 ~> s-1].
76  real :: ustar_min !< A minimum value of ustar to avoid numerical
77  !! problems [Z T-1 ~> m s-1]. If the value is small enough,
78  !! this should not affect the solution.
79  real :: tke_decay !< The ratio of the natural Ekman depth to the TKE
80  !! decay scale, nondimensional.
81  real :: omega_frac !< When setting the decay scale for turbulence, use
82  !! this fraction of the absolute rotation rate blended
83  !! with the local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).
84  logical :: answers_2018 !< If true, use the order of arithmetic and expressions that recover the
85  !! answers from the end of 2018. Otherwise, use updated and more robust
86  !! forms of the same expressions.
87  logical :: debug !< If true, write verbose checksums for debugging purposes.
88  type(ocean_obc_type), pointer :: obc => null() !< Open boundaries control structure
89  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
90  !! regulate the timing of diagnostic output.
91  !>@{ Diagnostics handles
92  integer :: id_bbl_thick_u = -1, id_kv_bbl_u = -1
93  integer :: id_bbl_thick_v = -1, id_kv_bbl_v = -1
94  integer :: id_ray_u = -1, id_ray_v = -1
95  integer :: id_nkml_visc_u = -1, id_nkml_visc_v = -1
96  !!@}
97 end type set_visc_cs
98 
99 contains
100 
101 !> Calculates the thickness of the bottom boundary layer and the viscosity within that layer.
102 !! A drag law is used, either linearized about an assumed bottom velocity or using
103 !! the actual near-bottom velocities combined with an assumed
104 !! unresolved velocity. The bottom boundary layer thickness is
105 !! limited by a combination of stratification and rotation, as in the
106 !! paper of Killworth and Edwards, JPO 1999. It is not necessary to
107 !! calculate the thickness and viscosity every time step; instead
108 !! previous values may be used.
109 subroutine set_viscous_bbl(u, v, h, tv, visc, G, GV, US, CS, symmetrize)
110  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
111  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
112  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
113  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
114  intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
115  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
116  intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
117  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
118  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
119  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
120  !! available thermodynamic fields. Absent fields
121  !! have NULL ptrs..
122  type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
123  !! related fields.
124  type(set_visc_cs), pointer :: cs !< The control structure returned by a previous
125  !! call to vertvisc_init.
126  logical, optional, intent(in) :: symmetrize !< If present and true, do extra calculations
127  !! of those values in visc that would be
128  !! calculated with symmetric memory.
129 
130  ! Local variables
131  real, dimension(SZIB_(G)) :: &
132  ustar, & ! The bottom friction velocity [Z T-1 ~> m s-1].
133  t_eos, & ! The temperature used to calculate the partial derivatives
134  ! of density with T and S [degC].
135  s_eos, & ! The salinity used to calculate the partial derivatives
136  ! of density with T and S [ppt].
137  dr_dt, & ! Partial derivative of the density in the bottom boundary
138  ! layer with temperature [R degC-1 ~> kg m-3 degC-1].
139  dr_ds, & ! Partial derivative of the density in the bottom boundary
140  ! layer with salinity [R ppt-1 ~> kg m-3 ppt-1].
141  press ! The pressure at which dR_dT and dR_dS are evaluated [Pa].
142  real :: htot ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
143  real :: htot_vel ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
144 
145  real :: rhtot ! Running sum of thicknesses times the layer potential
146  ! densities [H R ~> kg m-2 or kg2 m-5].
147  real, dimension(SZIB_(G),SZJ_(G)) :: &
148  d_u, & ! Bottom depth interpolated to u points [Z ~> m].
149  mask_u ! A mask that disables any contributions from u points that
150  ! are land or past open boundary conditions [nondim], 0 or 1.
151  real, dimension(SZI_(G),SZJB_(G)) :: &
152  d_v, & ! Bottom depth interpolated to v points [Z ~> m].
153  mask_v ! A mask that disables any contributions from v points that
154  ! are land or past open boundary conditions [nondim], 0 or 1.
155  real, dimension(SZIB_(G),SZK_(G)) :: &
156  h_at_vel, & ! Layer thickness at a velocity point, using an upwind-biased
157  ! second order accurate estimate based on the previous velocity
158  ! direction [H ~> m or kg m-2].
159  h_vel, & ! Arithmetic mean of the layer thicknesses adjacent to a
160  ! velocity point [H ~> m or kg m-2].
161  t_vel, & ! Arithmetic mean of the layer temperatures adjacent to a
162  ! velocity point [degC].
163  s_vel, & ! Arithmetic mean of the layer salinities adjacent to a
164  ! velocity point [ppt].
165  rml_vel ! Arithmetic mean of the layer coordinate densities adjacent
166  ! to a velocity point [R ~> kg m-3].
167 
168  real :: h_vel_pos ! The arithmetic mean thickness at a velocity point
169  ! plus H_neglect to avoid 0 values [H ~> m or kg m-2].
170  real :: ustarsq ! 400 times the square of ustar, times
171  ! Rho0 divided by G_Earth and the conversion
172  ! from m to thickness units [H R ~> kg m-2 or kg2 m-5].
173  real :: cdrag_sqrt_z ! Square root of the drag coefficient, times a unit conversion
174  ! factor from lateral lengths to vertical depths [Z L-1 ~> 1].
175  real :: cdrag_sqrt ! Square root of the drag coefficient [nondim].
176  real :: oldfn ! The integrated energy required to
177  ! entrain up to the bottom of the layer,
178  ! divided by G_Earth [H R ~> kg m-2 or kg2 m-5].
179  real :: dfn ! The increment in oldfn for entraining
180  ! the layer [H R ~> kg m-2 or kg2 m-5].
181  real :: dh ! The increment in layer thickness from
182  ! the present layer [H ~> m or kg m-2].
183  real :: bbl_thick ! The thickness of the bottom boundary layer [H ~> m or kg m-2].
184  real :: bbl_thick_z ! The thickness of the bottom boundary layer [Z ~> m].
185  real :: c2f ! C2f = 2*f at velocity points [T-1 ~> s-1].
186 
187  real :: u_bg_sq ! The square of an assumed background
188  ! velocity, for calculating the mean
189  ! magnitude near the bottom for use in the
190  ! quadratic bottom drag [L2 T-2 ~> m2 s-2].
191  real :: hwtot ! Sum of the thicknesses used to calculate
192  ! the near-bottom velocity magnitude [H ~> m or kg m-2].
193  real :: hutot ! Running sum of thicknesses times the
194  ! velocity magnitudes [H T T-1 ~> m2 s-1 or kg m-1 s-1].
195  real :: thtot ! Running sum of thickness times temperature [degC H ~> degC m or degC kg m-2].
196  real :: shtot ! Running sum of thickness times salinity [ppt H ~> ppt m or ppt kg m-2].
197  real :: hweight ! The thickness of a layer that is within Hbbl
198  ! of the bottom [H ~> m or kg m-2].
199  real :: v_at_u, u_at_v ! v at a u point or vice versa [L T-1 ~> m s-1].
200  real :: rho0x400_g ! 400*Rho0/G_Earth, times unit conversion factors
201  ! [R T2 H Z-2 ~> kg s2 m-4 or kg2 s2 m-7].
202  ! The 400 is a constant proposed by Killworth and Edwards, 1999.
203  real, dimension(SZI_(G),SZJ_(G),max(GV%nk_rho_varies,1)) :: &
204  rml ! The mixed layer coordinate density [R ~> kg m-3].
205  real :: p_ref(szi_(g)) ! The pressure used to calculate the coordinate
206  ! density [Pa] (usually set to 2e7 Pa = 2000 dbar).
207 
208  real :: d_vel ! The bottom depth at a velocity point [H ~> m or kg m-2].
209  real :: dp, dm ! The depths at the edges of a velocity cell [H ~> m or kg m-2].
210  real :: a ! a is the curvature of the bottom depth across a
211  ! cell, times the cell width squared [H ~> m or kg m-2].
212  real :: a_3, a_12 ! a/3 and a/12 [H ~> m or kg m-2].
213  real :: c24_a ! 24/a [H-1 ~> m-1 or m2 kg-1].
214  real :: slope ! The absolute value of the bottom depth slope across
215  ! a cell times the cell width [H ~> m or kg m-2].
216  real :: apb_4a, ax2_3apb ! Various nondimensional ratios of a and slope.
217  real :: a2x48_apb3, iapb, ibma_2 ! Combinations of a and slope [H-1 ~> m-1 or m2 kg-1].
218  ! All of the following "volumes" have units of thickness because they are normalized
219  ! by the full horizontal area of a velocity cell.
220  real :: vol_open ! The cell volume above which it is open [H ~> m or kg m-2].
221  real :: vol_direct ! With less than Vol_direct [H ~> m or kg m-2], there is a direct
222  ! solution of a cubic equation for L.
223  real :: vol_2_reg ! The cell volume above which there are two separate
224  ! open areas that must be integrated [H ~> m or kg m-2].
225  real :: vol ! The volume below the interface whose normalized
226  ! width is being sought [H ~> m or kg m-2].
227  real :: vol_below ! The volume below the interface below the one that
228  ! is currently under consideration [H ~> m or kg m-2].
229  real :: vol_err ! The error in the volume with the latest estimate of
230  ! L, or the error for the interface below [H ~> m or kg m-2].
231  real :: vol_quit ! The volume error below which to quit iterating [H ~> m or kg m-2].
232  real :: vol_tol ! A volume error tolerance [H ~> m or kg m-2].
233  real :: l(szk_(g)+1) ! The fraction of the full cell width that is open at
234  ! the depth of each interface [nondim].
235  real :: l_direct ! The value of L above volume Vol_direct [nondim].
236  real :: l_max, l_min ! Upper and lower bounds on the correct value for L [nondim].
237  real :: vol_err_max ! The volume errors for the upper and lower bounds on
238  real :: vol_err_min ! the correct value for L [H ~> m or kg m-2].
239  real :: vol_0 ! A deeper volume with known width L0 [H ~> m or kg m-2].
240  real :: l0 ! The value of L above volume Vol_0 [nondim].
241  real :: dvol ! vol - Vol_0 [H ~> m or kg m-2].
242  real :: dv_dl2 ! The partial derivative of volume with L squared
243  ! evaluated at L=L0 [H ~> m or kg m-2].
244  real :: h_neglect ! A thickness that is so small it is usually lost
245  ! in roundoff and can be neglected [H ~> m or kg m-2].
246  real :: usth ! ustar converted to units of H T-1 [H T-1 ~> m s-1 or kg m-2 s-1].
247  real :: root ! A temporary variable [H T-1 ~> m s-1 or kg m-2 s-1].
248 
249  real :: cell_width ! The transverse width of the velocity cell [L ~> m].
250  real :: rayleigh ! A nondimensional value that is multiplied by the layer's
251  ! velocity magnitude to give the Rayleigh drag velocity, times
252  ! a lateral to vertical distance conversion factor [Z L-1 ~> 1].
253  real :: gam ! The ratio of the change in the open interface width
254  ! to the open interface width atop a cell [nondim].
255  real :: bbl_frac ! The fraction of a layer's drag that goes into the
256  ! viscous bottom boundary layer [nondim].
257  real :: bbl_visc_frac ! The fraction of all the drag that is expressed as
258  ! a viscous bottom boundary layer [nondim].
259  real, parameter :: c1_3 = 1.0/3.0, c1_6 = 1.0/6.0, c1_12 = 1.0/12.0
260  real :: c2pi_3 ! An irrational constant, 2/3 pi.
261  real :: tmp ! A temporary variable.
262  real :: tmp_val_m1_to_p1
263  real :: curv_tol ! Numerator of curvature cubed, used to estimate
264  ! accuracy of a single L(:) Newton iteration
265  logical :: use_l0, do_one_l_iter ! Control flags for L(:) Newton iteration
266  logical :: use_bbl_eos, do_i(szib_(g))
267  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, m, n, k2, nkmb, nkml
268  integer :: itt, maxitt=20
269  type(ocean_obc_type), pointer :: obc => null()
270 
271  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
272  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
273  nkmb = gv%nk_rho_varies ; nkml = gv%nkml
274  h_neglect = gv%H_subroundoff
275  rho0x400_g = 400.0*(gv%Rho0 / (us%L_to_Z**2 * gv%g_Earth)) * gv%Z_to_H
276  vol_quit = 0.9*gv%Angstrom_H + h_neglect
277  c2pi_3 = 8.0*atan(1.0)/3.0
278 
279  if (.not.associated(cs)) call mom_error(fatal,"MOM_set_viscosity(BBL): "//&
280  "Module must be initialized before it is used.")
281  if (.not.cs%bottomdraglaw) return
282 
283  if (present(symmetrize)) then ; if (symmetrize) then
284  jsq = js-1 ; isq = is-1
285  endif ; endif
286 
287  if (cs%debug) then
288  call uvchksum("Start set_viscous_BBL [uv]", u, v, g%HI, haloshift=1, scale=us%L_T_to_m_s)
289  call hchksum(h,"Start set_viscous_BBL h", g%HI, haloshift=1, scale=gv%H_to_m)
290  if (associated(tv%T)) call hchksum(tv%T, "Start set_viscous_BBL T", g%HI, haloshift=1)
291  if (associated(tv%S)) call hchksum(tv%S, "Start set_viscous_BBL S", g%HI, haloshift=1)
292  endif
293 
294  use_bbl_eos = associated(tv%eqn_of_state) .and. cs%BBL_use_EOS
295  obc => cs%OBC
296 
297  u_bg_sq = cs%drag_bg_vel * cs%drag_bg_vel
298  cdrag_sqrt = sqrt(cs%cdrag)
299  cdrag_sqrt_z = us%L_to_Z * sqrt(cs%cdrag)
300  k2 = max(nkmb+1, 2)
301 
302 ! With a linear drag law, the friction velocity is already known.
303 ! if (CS%linear_drag) ustar(:) = cdrag_sqrt_Z*CS%drag_bg_vel
304 
305  if ((nkml>0) .and. .not.use_bbl_eos) then
306  do i=isq,ieq+1 ; p_ref(i) = tv%P_ref ; enddo
307  !$OMP parallel do default(shared)
308  do j=jsq,jeq+1 ; do k=1,nkmb
309  call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, &
310  rml(:,j,k), isq, ieq-isq+2, tv%eqn_of_state, scale=us%kg_m3_to_R)
311  enddo ; enddo
312  endif
313 
314  !$OMP parallel do default(shared)
315  do j=js-1,je ; do i=is-1,ie+1
316  d_v(i,j) = 0.5*(g%bathyT(i,j) + g%bathyT(i,j+1))
317  mask_v(i,j) = g%mask2dCv(i,j)
318  enddo ; enddo
319  !$OMP parallel do default(shared)
320  do j=js-1,je+1 ; do i=is-1,ie
321  d_u(i,j) = 0.5*(g%bathyT(i,j) + g%bathyT(i+1,j))
322  mask_u(i,j) = g%mask2dCu(i,j)
323  enddo ; enddo
324 
325  if (associated(obc)) then ; do n=1,obc%number_of_segments
326  if (.not. obc%segment(n)%on_pe) cycle
327  ! Use a one-sided projection of bottom depths at OBC points.
328  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
329  if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je)) then
330  do i = max(is-1,obc%segment(n)%HI%isd), min(ie+1,obc%segment(n)%HI%ied)
331  if (obc%segment(n)%direction == obc_direction_n) d_v(i,j) = g%bathyT(i,j)
332  if (obc%segment(n)%direction == obc_direction_s) d_v(i,j) = g%bathyT(i,j+1)
333  enddo
334  elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ie)) then
335  do j = max(js-1,obc%segment(n)%HI%jsd), min(je+1,obc%segment(n)%HI%jed)
336  if (obc%segment(n)%direction == obc_direction_e) d_u(i,j) = g%bathyT(i,j)
337  if (obc%segment(n)%direction == obc_direction_w) d_u(i,j) = g%bathyT(i+1,j)
338  enddo
339  endif
340  enddo ; endif
341  if (associated(obc)) then ; do n=1,obc%number_of_segments
342  ! Now project bottom depths across cell-corner points in the OBCs. The two
343  ! projections have to occur in sequence and can not be combined easily.
344  if (.not. obc%segment(n)%on_pe) cycle
345  ! Use a one-sided projection of bottom depths at OBC points.
346  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
347  if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je)) then
348  do i = max(is-1,obc%segment(n)%HI%IsdB), min(ie,obc%segment(n)%HI%IedB)
349  if (obc%segment(n)%direction == obc_direction_n) then
350  d_u(i,j+1) = d_u(i,j) ; mask_u(i,j+1) = 0.0
351  elseif (obc%segment(n)%direction == obc_direction_s) then
352  d_u(i,j) = d_u(i,j+1) ; mask_u(i,j) = 0.0
353  endif
354  enddo
355  elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ie)) then
356  do j = max(js-1,obc%segment(n)%HI%JsdB), min(je,obc%segment(n)%HI%JedB)
357  if (obc%segment(n)%direction == obc_direction_e) then
358  d_v(i+1,j) = d_v(i,j) ; mask_v(i+1,j) = 0.0
359  elseif (obc%segment(n)%direction == obc_direction_w) then
360  d_v(i,j) = d_v(i+1,j) ; mask_v(i,j) = 0.0
361  endif
362  enddo
363  endif
364  enddo ; endif
365 
366  if (.not.use_bbl_eos) rml_vel(:,:) = 0.0
367 
368  !$OMP parallel do default(private) shared(u,v,h,tv,visc,G,GV,US,CS,Rml,is,ie,js,je,nz,nkmb, &
369  !$OMP nkml,Isq,Ieq,Jsq,Jeq,h_neglect,Rho0x400_G,C2pi_3, &
370  !$OMP U_bg_sq,cdrag_sqrt_Z,cdrag_sqrt,K2,use_BBL_EOS, &
371  !$OMP OBC,maxitt,Vol_quit,D_u,D_v,mask_u,mask_v)
372  do j=jsq,jeq ; do m=1,2
373 
374  if (m==1) then
375  if (j<g%Jsc) cycle
376  is = isq ; ie = ieq
377  do i=is,ie
378  do_i(i) = .false.
379  if (g%mask2dCu(i,j) > 0) do_i(i) = .true.
380  enddo
381  else
382  is = g%isc ; ie = g%iec
383  do i=is,ie
384  do_i(i) = .false.
385  if (g%mask2dCv(i,j) > 0) do_i(i) = .true.
386  enddo
387  endif
388 
389  if (m==1) then
390  do k=1,nz ; do i=is,ie
391  if (do_i(i)) then
392  if (u(i,j,k) *(h(i+1,j,k) - h(i,j,k)) >= 0) then
393  h_at_vel(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / &
394  (h(i,j,k) + h(i+1,j,k) + h_neglect)
395  else
396  h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
397  endif
398  endif
399  h_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
400  enddo ; enddo
401  if (use_bbl_eos) then ; do k=1,nz ; do i=is,ie
402  ! Perhaps these should be thickness weighted.
403  t_vel(i,k) = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
404  s_vel(i,k) = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
405  enddo ; enddo ; else ; do k=1,nkmb ; do i=is,ie
406  rml_vel(i,k) = 0.5 * (rml(i,j,k) + rml(i+1,j,k))
407  enddo ; enddo ; endif
408  else
409  do k=1,nz ; do i=is,ie
410  if (do_i(i)) then
411  if (v(i,j,k) * (h(i,j+1,k) - h(i,j,k)) >= 0) then
412  h_at_vel(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / &
413  (h(i,j,k) + h(i,j+1,k) + h_neglect)
414  else
415  h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
416  endif
417  endif
418  h_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
419  enddo ; enddo
420  if (use_bbl_eos) then ; do k=1,nz ; do i=is,ie
421  t_vel(i,k) = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
422  s_vel(i,k) = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
423  enddo ; enddo ; else ; do k=1,nkmb ; do i=is,ie
424  rml_vel(i,k) = 0.5 * (rml(i,j,k) + rml(i,j+1,k))
425  enddo ; enddo ; endif
426  endif
427 
428  if (associated(obc)) then ; if (obc%number_of_segments > 0) then
429  ! Apply a zero gradient projection of thickness across OBC points.
430  if (m==1) then
431  do i=is,ie ; if (do_i(i) .and. (obc%segnum_u(i,j) /= obc_none)) then
432  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
433  do k=1,nz
434  h_at_vel(i,k) = h(i,j,k) ; h_vel(i,k) = h(i,j,k)
435  enddo
436  if (use_bbl_eos) then
437  do k=1,nz
438  t_vel(i,k) = tv%T(i,j,k) ; s_vel(i,k) = tv%S(i,j,k)
439  enddo
440  else
441  do k=1,nkmb
442  rml_vel(i,k) = rml(i,j,k)
443  enddo
444  endif
445  elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
446  do k=1,nz
447  h_at_vel(i,k) = h(i+1,j,k) ; h_vel(i,k) = h(i+1,j,k)
448  enddo
449  if (use_bbl_eos) then
450  do k=1,nz
451  t_vel(i,k) = tv%T(i+1,j,k) ; s_vel(i,k) = tv%S(i+1,j,k)
452  enddo
453  else
454  do k=1,nkmb
455  rml_vel(i,k) = rml(i+1,j,k)
456  enddo
457  endif
458  endif
459  endif ; enddo
460  else
461  do i=is,ie ; if (do_i(i) .and. (obc%segnum_v(i,j) /= obc_none)) then
462  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
463  do k=1,nz
464  h_at_vel(i,k) = h(i,j,k) ; h_vel(i,k) = h(i,j,k)
465  enddo
466  if (use_bbl_eos) then
467  do k=1,nz
468  t_vel(i,k) = tv%T(i,j,k) ; s_vel(i,k) = tv%S(i,j,k)
469  enddo
470  else
471  do k=1,nkmb
472  rml_vel(i,k) = rml(i,j,k)
473  enddo
474  endif
475  elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
476  do k=1,nz
477  h_at_vel(i,k) = h(i,j+1,k) ; h_vel(i,k) = h(i,j+1,k)
478  enddo
479  if (use_bbl_eos) then
480  do k=1,nz
481  t_vel(i,k) = tv%T(i,j+1,k) ; s_vel(i,k) = tv%S(i,j+1,k)
482  enddo
483  else
484  do k=1,nkmb
485  rml_vel(i,k) = rml(i,j+1,k)
486  enddo
487  endif
488  endif
489  endif ; enddo
490  endif
491  endif ; endif
492 
493  if (use_bbl_eos .or. .not.cs%linear_drag) then
494  do i=is,ie ; if (do_i(i)) then
495 ! This block of code calculates the mean velocity magnitude over
496 ! the bottommost CS%Hbbl of the water column for determining
497 ! the quadratic bottom drag.
498  htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
499  thtot = 0.0 ; shtot = 0.0
500  do k=nz,1,-1
501 
502  if (htot_vel>=cs%Hbbl) exit ! terminate the k loop
503 
504  hweight = min(cs%Hbbl - htot_vel, h_at_vel(i,k))
505  if (hweight < 1.5*gv%Angstrom_H + h_neglect) cycle
506 
507  htot_vel = htot_vel + h_at_vel(i,k)
508  hwtot = hwtot + hweight
509 
510  if ((.not.cs%linear_drag) .and. (hweight >= 0.0)) then ; if (m==1) then
511  v_at_u = set_v_at_u(v, h, g, i, j, k, mask_v, obc)
512  hutot = hutot + hweight * sqrt(u(i,j,k)*u(i,j,k) + &
513  v_at_u*v_at_u + u_bg_sq)
514  else
515  u_at_v = set_u_at_v(u, h, g, i, j, k, mask_u, obc)
516  hutot = hutot + hweight * sqrt(v(i,j,k)*v(i,j,k) + &
517  u_at_v*u_at_v + u_bg_sq)
518  endif ; endif
519 
520  if (use_bbl_eos .and. (hweight >= 0.0)) then
521  thtot = thtot + hweight * t_vel(i,k)
522  shtot = shtot + hweight * s_vel(i,k)
523  endif
524  enddo ! end of k loop
525 
526  if (.not.cs%linear_drag .and. (hwtot > 0.0)) then
527  ustar(i) = cdrag_sqrt_z*hutot/hwtot
528  else
529  ustar(i) = cdrag_sqrt_z*cs%drag_bg_vel
530  endif
531 
532  if (use_bbl_eos) then ; if (hwtot > 0.0) then
533  t_eos(i) = thtot/hwtot ; s_eos(i) = shtot/hwtot
534  else
535  t_eos(i) = 0.0 ; s_eos(i) = 0.0
536  endif ; endif
537  endif ; enddo
538  else
539  do i=is,ie ; ustar(i) = cdrag_sqrt_z*cs%drag_bg_vel ; enddo
540  endif ! Not linear_drag
541 
542  if (use_bbl_eos) then
543  do i=is,ie
544  press(i) = 0.0 ! or = forces%p_surf(i,j)
545  if (.not.do_i(i)) then ; t_eos(i) = 0.0 ; s_eos(i) = 0.0 ; endif
546  enddo
547  do k=1,nz ; do i=is,ie
548  press(i) = press(i) + gv%H_to_Pa * h_vel(i,k)
549  enddo ; enddo
550  call calculate_density_derivs(t_eos, s_eos, press, dr_dt, dr_ds, &
551  is-g%IsdB+1, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
552  endif
553 
554  do i=is,ie ; if (do_i(i)) then
555 ! The 400.0 in this expression is the square of a constant proposed
556 ! by Killworth and Edwards, 1999, in equation (2.20).
557  ustarsq = rho0x400_g * ustar(i)**2
558  htot = 0.0
559 
560 ! This block of code calculates the thickness of a stratification
561 ! limited bottom boundary layer, using the prescription from
562 ! Killworth and Edwards, 1999, as described in Stephens and Hallberg
563 ! 2000 (unpublished and lost manuscript).
564  if (use_bbl_eos) then
565  thtot = 0.0 ; shtot = 0.0 ; oldfn = 0.0
566  do k=nz,2,-1
567  if (h_at_vel(i,k) <= 0.0) cycle
568 
569  oldfn = dr_dt(i)*(thtot - t_vel(i,k)*htot) + &
570  dr_ds(i)*(shtot - s_vel(i,k)*htot)
571  if (oldfn >= ustarsq) exit
572 
573  dfn = (dr_dt(i)*(t_vel(i,k) - t_vel(i,k-1)) + &
574  dr_ds(i)*(s_vel(i,k) - s_vel(i,k-1))) * &
575  (h_at_vel(i,k) + htot)
576 
577  if ((oldfn + dfn) <= ustarsq) then
578  dh = h_at_vel(i,k)
579  else
580  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
581  endif
582 
583  htot = htot + dh
584  thtot = thtot + t_vel(i,k)*dh ; shtot = shtot + s_vel(i,k)*dh
585  enddo
586  if ((oldfn < ustarsq) .and. h_at_vel(i,1) > 0.0) then
587  ! Layer 1 might be part of the BBL.
588  if (dr_dt(i) * (thtot - t_vel(i,1)*htot) + &
589  dr_ds(i) * (shtot - s_vel(i,1)*htot) < ustarsq) &
590  htot = htot + h_at_vel(i,1)
591  endif ! Examination of layer 1.
592  else ! Use Rlay and/or the coordinate density as density variables.
593  rhtot = 0.0
594  do k=nz,k2,-1
595  oldfn = rhtot - gv%Rlay(k)*htot
596  dfn = (gv%Rlay(k) - gv%Rlay(k-1))*(h_at_vel(i,k)+htot)
597 
598  if (oldfn >= ustarsq) then
599  cycle
600  elseif ((oldfn + dfn) <= ustarsq) then
601  dh = h_at_vel(i,k)
602  else
603  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
604  endif
605 
606  htot = htot + dh
607  rhtot = rhtot + gv%Rlay(k)*dh
608  enddo
609  if (nkml>0) then
610  do k=nkmb,2,-1
611  oldfn = rhtot - rml_vel(i,k)*htot
612  dfn = (rml_vel(i,k) - rml_vel(i,k-1)) * (h_at_vel(i,k)+htot)
613 
614  if (oldfn >= ustarsq) then
615  cycle
616  elseif ((oldfn + dfn) <= ustarsq) then
617  dh = h_at_vel(i,k)
618  else
619  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
620  endif
621 
622  htot = htot + dh
623  rhtot = rhtot + rml_vel(i,k)*dh
624  enddo
625  if (rhtot - rml_vel(i,1)*htot < ustarsq) htot = htot + h_at_vel(i,1)
626  else
627  if (rhtot - gv%Rlay(1)*htot < ustarsq) htot = htot + h_at_vel(i,1)
628  endif
629  endif ! use_BBL_EOS
630 
631 ! The Coriolis limit is 0.5*ustar/f. The buoyancy limit here is htot.
632 ! The bottom boundary layer thickness is found by solving the same
633 ! equation as in Killworth and Edwards: (h/h_f)^2 + h/h_N = 1.
634 
635  if (m==1) then ; c2f = g%CoriolisBu(i,j-1) + g%CoriolisBu(i,j)
636  else ; c2f = g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j) ; endif
637 
638  if (cs%cdrag * u_bg_sq <= 0.0) then
639  ! This avoids NaNs and overflows, and could be used in all cases,
640  ! but is not bitwise identical to the current code.
641  usth = ustar(i)*gv%Z_to_H ; root = sqrt(0.25*usth**2 + (htot*c2f)**2)
642  if (htot*usth <= (cs%BBL_thick_min+h_neglect) * (0.5*usth + root)) then
643  bbl_thick = cs%BBL_thick_min
644  else
645  bbl_thick = (htot * usth) / (0.5*usth + root)
646  endif
647  else
648  bbl_thick = htot / (0.5 + sqrt(0.25 + htot*htot*c2f*c2f/ &
649  ((ustar(i)*ustar(i)) * (gv%Z_to_H**2)) ) )
650 
651  if (bbl_thick < cs%BBL_thick_min) bbl_thick = cs%BBL_thick_min
652  endif
653 ! If there is Richardson number dependent mixing, that determines
654 ! the vertical extent of the bottom boundary layer, and there is no
655 ! need to set that scale here. In fact, viscously reducing the
656 ! shears over an excessively large region reduces the efficacy of
657 ! the Richardson number dependent mixing.
658  if ((bbl_thick > 0.5*cs%Hbbl) .and. (cs%RiNo_mix)) bbl_thick = 0.5*cs%Hbbl
659 
660  if (cs%Channel_drag) then
661  ! The drag within the bottommost bbl_thick is applied as a part of
662  ! an enhanced bottom viscosity, while above this the drag is applied
663  ! directly to the layers in question as a Rayleigh drag term.
664  if (m==1) then
665  d_vel = d_u(i,j)
666  tmp = g%mask2dCu(i,j+1) * d_u(i,j+1)
667  dp = 2.0 * d_vel * tmp / (d_vel + tmp)
668  tmp = g%mask2dCu(i,j-1) * d_u(i,j-1)
669  dm = 2.0 * d_vel * tmp / (d_vel + tmp)
670  else
671  d_vel = d_v(i,j)
672  tmp = g%mask2dCv(i+1,j) * d_v(i+1,j)
673  dp = 2.0 * d_vel * tmp / (d_vel + tmp)
674  tmp = g%mask2dCv(i-1,j) * d_v(i-1,j)
675  dm = 2.0 * d_vel * tmp / (d_vel + tmp)
676  endif
677  if (dm > dp) then ; tmp = dp ; dp = dm ; dm = tmp ; endif
678 
679  ! Convert the D's to the units of thickness.
680  dp = gv%Z_to_H*dp ; dm = gv%Z_to_H*dm ; d_vel = gv%Z_to_H*d_vel
681 
682  a_3 = (dp + dm - 2.0*d_vel) ; a = 3.0*a_3 ; a_12 = 0.25*a_3
683  slope = dp - dm
684  ! If the curvature is small enough, there is no reason not to assume
685  ! a uniformly sloping or flat bottom.
686  if (abs(a) < 1e-2*(slope + cs%BBL_thick_min)) a = 0.0
687  ! Each cell extends from x=-1/2 to 1/2, and has a topography
688  ! given by D(x) = a*x^2 + b*x + D - a/12.
689 
690  ! Calculate the volume above which the entire cell is open and the
691  ! other volumes at which the equation that is solved for L changes.
692  if (a > 0.0) then
693  if (slope >= a) then
694  vol_open = d_vel - dm ; vol_2_reg = vol_open
695  else
696  tmp = slope/a
697  vol_open = 0.25*slope*tmp + c1_12*a
698  vol_2_reg = 0.5*tmp**2 * (a - c1_3*slope)
699  endif
700  ! Define some combinations of a & b for later use.
701  c24_a = 24.0/a ; iapb = 1.0/(a+slope)
702  apb_4a = (slope+a)/(4.0*a) ; a2x48_apb3 = (48.0*(a*a))*(iapb**3)
703  ax2_3apb = 2.0*c1_3*a*iapb
704  elseif (a == 0.0) then
705  vol_open = 0.5*slope
706  if (slope > 0) iapb = 1.0/slope
707  else ! a < 0.0
708  vol_open = d_vel - dm
709  if (slope >= -a) then
710  iapb = 1.0e30 ; if (slope+a /= 0.0) iapb = 1.0/(a+slope)
711  vol_direct = 0.0 ; l_direct = 0.0 ; c24_a = 0.0
712  else
713  c24_a = 24.0/a ; iapb = 1.0/(a+slope)
714  l_direct = 1.0 + slope/a ! L_direct < 1 because a < 0
715  vol_direct = -c1_6*a*l_direct**3
716  endif
717  ibma_2 = 2.0 / (slope - a)
718  endif
719 
720  l(nz+1) = 0.0 ; vol = 0.0 ; vol_err = 0.0 ; bbl_visc_frac = 0.0
721  ! Determine the normalized open length at each interface.
722  do k=nz,1,-1
723  vol_below = vol
724 
725  vol = vol + h_vel(i,k)
726  h_vel_pos = h_vel(i,k) + h_neglect
727 
728  if (vol >= vol_open) then ; l(k) = 1.0
729  elseif (a == 0) then ! The bottom has no curvature.
730  l(k) = sqrt(2.0*vol*iapb)
731  elseif (a > 0) then
732  ! There may be a minimum depth, and there are
733  ! analytic expressions for L for all cases.
734  if (vol < vol_2_reg) then
735  ! In this case, there is a contiguous open region and
736  ! vol = 0.5*L^2*(slope + a/3*(3-4L)).
737  if (a2x48_apb3*vol < 1e-8) then ! Could be 1e-7?
738  ! There is a very good approximation here for massless layers.
739  l0 = sqrt(2.0*vol*iapb) ; l(k) = l0*(1.0 + ax2_3apb*l0)
740  else
741  l(k) = apb_4a * (1.0 - &
742  2.0 * cos(c1_3*acos(a2x48_apb3*vol - 1.0) - c2pi_3))
743  endif
744  ! To check the answers.
745  ! Vol_err = 0.5*(L(K)*L(K))*(slope + a_3*(3.0-4.0*L(K))) - vol
746  else ! There are two separate open regions.
747  ! vol = slope^2/4a + a/12 - (a/12)*(1-L)^2*(1+2L)
748  ! At the deepest volume, L = slope/a, at the top L = 1.
749  !L(K) = 0.5 - cos(C1_3*acos(1.0 - C24_a*(Vol_open - vol)) - C2pi_3)
750  tmp_val_m1_to_p1 = 1.0 - c24_a*(vol_open - vol)
751  tmp_val_m1_to_p1 = max(-1., min(1., tmp_val_m1_to_p1))
752  l(k) = 0.5 - cos(c1_3*acos(tmp_val_m1_to_p1) - c2pi_3)
753  ! To check the answers.
754  ! Vol_err = Vol_open - a_12*(1.0+2.0*L(K)) * (1.0-L(K))**2 - vol
755  endif
756  else ! a < 0.
757  if (vol <= vol_direct) then
758  ! Both edges of the cell are bounded by walls.
759  l(k) = (-0.25*c24_a*vol)**c1_3
760  else
761  ! x_R is at 1/2 but x_L is in the interior & L is found by solving
762  ! vol = 0.5*L^2*(slope + a/3*(3-4L))
763 
764  ! Vol_err = 0.5*(L(K+1)*L(K+1))*(slope + a_3*(3.0-4.0*L(K+1))) - vol_below
765  ! Change to ...
766  ! if (min(Vol_below + Vol_err, vol) <= Vol_direct) then ?
767  if (vol_below + vol_err <= vol_direct) then
768  l0 = l_direct ; vol_0 = vol_direct
769  else
770  l0 = l(k+1) ; vol_0 = vol_below + vol_err
771  ! Change to Vol_0 = min(Vol_below + Vol_err, vol) ?
772  endif
773 
774  ! Try a relatively simple solution that usually works well
775  ! for massless layers.
776  dv_dl2 = 0.5*(slope+a) - a*l0 ; dvol = (vol-vol_0)
777  ! dV_dL2 = 0.5*(slope+a) - a*L0 ; dVol = max(vol-Vol_0, 0.0)
778 
779  use_l0 = .false.
780  do_one_l_iter = .false.
781  if (cs%answers_2018) then
782  curv_tol = gv%Angstrom_H*dv_dl2**2 &
783  * (0.25 * dv_dl2 * gv%Angstrom_H - a * l0 * dvol)
784  do_one_l_iter = (a * a * dvol**3) < curv_tol
785  else
786  ! The following code is more robust when GV%Angstrom_H=0, but
787  ! it changes answers.
788  use_l0 = (dvol <= 0.)
789 
790  vol_tol = max(0.5 * gv%Angstrom_H + gv%H_subroundoff, 1e-14 * vol)
791  vol_quit = max(0.9 * gv%Angstrom_H + gv%H_subroundoff, 1e-14 * vol)
792 
793  curv_tol = vol_tol * dv_dl2**2 &
794  * (dv_dl2 * vol_tol - 2.0 * a * l0 * dvol)
795  do_one_l_iter = (a * a * dvol**3) < curv_tol
796  endif
797 
798  if (use_l0) then
799  l(k) = l0
800  vol_err = 0.5*(l(k)*l(k))*(slope + a_3*(3.0-4.0*l(k))) - vol
801  elseif (do_one_l_iter) then
802  ! One iteration of Newton's method should give an estimate
803  ! that is accurate to within Vol_tol.
804  l(k) = sqrt(l0*l0 + dvol / dv_dl2)
805  vol_err = 0.5*(l(k)*l(k))*(slope + a_3*(3.0-4.0*l(k))) - vol
806  else
807  if (dv_dl2*(1.0-l0*l0) < dvol + &
808  dv_dl2 * (vol_open - vol)*ibma_2) then
809  l_max = sqrt(1.0 - (vol_open - vol)*ibma_2)
810  else
811  l_max = sqrt(l0*l0 + dvol / dv_dl2)
812  endif
813  l_min = sqrt(l0*l0 + dvol / (0.5*(slope+a) - a*l_max))
814 
815  vol_err_min = 0.5*(l_min**2)*(slope + a_3*(3.0-4.0*l_min)) - vol
816  vol_err_max = 0.5*(l_max**2)*(slope + a_3*(3.0-4.0*l_max)) - vol
817  ! if ((abs(Vol_err_min) <= Vol_quit) .or. (Vol_err_min >= Vol_err_max)) then
818  if (abs(vol_err_min) <= vol_quit) then
819  l(k) = l_min ; vol_err = vol_err_min
820  else
821  l(k) = sqrt((l_min**2*vol_err_max - l_max**2*vol_err_min) / &
822  (vol_err_max - vol_err_min))
823  do itt=1,maxitt
824  vol_err = 0.5*(l(k)*l(k))*(slope + a_3*(3.0-4.0*l(k))) - vol
825  if (abs(vol_err) <= vol_quit) exit
826  ! Take a Newton's method iteration. This equation has proven
827  ! robust enough not to need bracketing.
828  l(k) = l(k) - vol_err / (l(k)* (slope + a - 2.0*a*l(k)))
829  ! This would be a Newton's method iteration for L^2:
830  ! L(K) = sqrt(L(K)*L(K) - Vol_err / (0.5*(slope+a) - a*L(K)))
831  enddo
832  endif ! end of iterative solver
833  endif ! end of 1-boundary alternatives.
834  endif ! end of a<0 cases.
835  endif
836 
837  ! Determine the drag contributing to the bottom boundary layer
838  ! and the Raleigh drag that acts on each layer.
839  if (l(k) > l(k+1)) then
840  if (vol_below < bbl_thick) then
841  bbl_frac = (1.0-vol_below/bbl_thick)**2
842  bbl_visc_frac = bbl_visc_frac + bbl_frac*(l(k) - l(k+1))
843  else
844  bbl_frac = 0.0
845  endif
846 
847  if (m==1) then ; cell_width = g%dy_Cu(i,j)
848  else ; cell_width = g%dx_Cv(i,j) ; endif
849  gam = 1.0 - l(k+1)/l(k)
850  rayleigh = us%L_to_Z * cs%cdrag * (l(k)-l(k+1)) * (1.0-bbl_frac) * &
851  (12.0*cs%c_Smag*h_vel_pos) / (12.0*cs%c_Smag*h_vel_pos + &
852  us%L_to_Z*gv%Z_to_H * cs%cdrag * gam*(1.0-gam)*(1.0-1.5*gam) * l(k)**2 * cell_width)
853  else ! This layer feels no drag.
854  rayleigh = 0.0
855  endif
856 
857  if (m==1) then
858  if (rayleigh > 0.0) then
859  v_at_u = set_v_at_u(v, h, g, i, j, k, mask_v, obc)
860  visc%Ray_u(i,j,k) = rayleigh*sqrt(u(i,j,k)*u(i,j,k) + &
861  v_at_u*v_at_u + u_bg_sq)
862  else ; visc%Ray_u(i,j,k) = 0.0 ; endif
863  else
864  if (rayleigh > 0.0) then
865  u_at_v = set_u_at_v(u, h, g, i, j, k, mask_u, obc)
866  visc%Ray_v(i,j,k) = rayleigh*sqrt(v(i,j,k)*v(i,j,k) + &
867  u_at_v*u_at_v + u_bg_sq)
868  else ; visc%Ray_v(i,j,k) = 0.0 ; endif
869  endif
870 
871  enddo ! k loop to determine L(K).
872 
873  bbl_thick_z = bbl_thick * gv%H_to_Z
874  if (m==1) then
875  visc%Kv_bbl_u(i,j) = max(cs%Kv_BBL_min, &
876  cdrag_sqrt*ustar(i)*bbl_thick_z*bbl_visc_frac)
877  visc%bbl_thick_u(i,j) = bbl_thick_z
878  else
879  visc%Kv_bbl_v(i,j) = max(cs%Kv_BBL_min, &
880  cdrag_sqrt*ustar(i)*bbl_thick_z*bbl_visc_frac)
881  visc%bbl_thick_v(i,j) = bbl_thick_z
882  endif
883 
884  else ! Not Channel_drag.
885 ! Here the near-bottom viscosity is set to a value which will give
886 ! the correct stress when the shear occurs over bbl_thick.
887  bbl_thick_z = bbl_thick * gv%H_to_Z
888  if (m==1) then
889  visc%Kv_bbl_u(i,j) = max(cs%Kv_BBL_min, cdrag_sqrt*ustar(i)*bbl_thick_z)
890  visc%bbl_thick_u(i,j) = bbl_thick_z
891  else
892  visc%Kv_bbl_v(i,j) = max(cs%Kv_BBL_min, cdrag_sqrt*ustar(i)*bbl_thick_z)
893  visc%bbl_thick_v(i,j) = bbl_thick_z
894  endif
895  endif
896  endif ; enddo ! end of i loop
897  enddo ; enddo ! end of m & j loops
898 
899 ! Offer diagnostics for averaging
900  if (cs%id_bbl_thick_u > 0) &
901  call post_data(cs%id_bbl_thick_u, visc%bbl_thick_u, cs%diag)
902  if (cs%id_kv_bbl_u > 0) &
903  call post_data(cs%id_kv_bbl_u, visc%kv_bbl_u, cs%diag)
904  if (cs%id_bbl_thick_v > 0) &
905  call post_data(cs%id_bbl_thick_v, visc%bbl_thick_v, cs%diag)
906  if (cs%id_kv_bbl_v > 0) &
907  call post_data(cs%id_kv_bbl_v, visc%kv_bbl_v, cs%diag)
908  if (cs%id_Ray_u > 0) &
909  call post_data(cs%id_Ray_u, visc%Ray_u, cs%diag)
910  if (cs%id_Ray_v > 0) &
911  call post_data(cs%id_Ray_v, visc%Ray_v, cs%diag)
912 
913  if (cs%debug) then
914  if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) &
915  call uvchksum("Ray [uv]", visc%Ray_u, visc%Ray_v, g%HI, haloshift=0, scale=us%Z_to_m)
916  if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) &
917  call uvchksum("kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
918  if (associated(visc%bbl_thick_u) .and. associated(visc%bbl_thick_v)) &
919  call uvchksum("bbl_thick_[uv]", visc%bbl_thick_u, &
920  visc%bbl_thick_v, g%HI, haloshift=0, scale=us%Z_to_m)
921  endif
922 
923 end subroutine set_viscous_bbl
924 
925 !> This subroutine finds a thickness-weighted value of v at the u-points.
926 function set_v_at_u(v, h, G, i, j, k, mask2dCv, OBC)
927  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
928  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
929  intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
930  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
931  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
932  integer, intent(in) :: i !< The i-index of the u-location to work on.
933  integer, intent(in) :: j !< The j-index of the u-location to work on.
934  integer, intent(in) :: k !< The k-index of the u-location to work on.
935  real, dimension(SZI_(G),SZJB_(G)),&
936  intent(in) :: mask2dcv !< A multiplicative mask of the v-points
937  type(ocean_obc_type), pointer :: obc !< A pointer to an open boundary condition structure
938  real :: set_v_at_u !< The return value of v at u points points in the
939  !! same units as u, i.e. [L T-1 ~> m s-1] or other units.
940 
941  ! This subroutine finds a thickness-weighted value of v at the u-points.
942  real :: hwt(0:1,-1:0) ! Masked weights used to average u onto v [H ~> m or kg m-2].
943  real :: hwt_tot ! The sum of the masked thicknesses [H ~> m or kg m-2].
944  integer :: i0, j0, i1, j1
945 
946  do j0 = -1,0 ; do i0 = 0,1 ; i1 = i+i0 ; j1 = j+j0
947  hwt(i0,j0) = (h(i1,j1,k) + h(i1,j1+1,k)) * mask2dcv(i1,j1)
948  enddo ; enddo
949 
950  if (associated(obc)) then ; if (obc%number_of_segments > 0) then
951  do j0 = -1,0 ; do i0 = 0,1 ; if ((obc%segnum_v(i+i0,j+j0) /= obc_none)) then
952  i1 = i+i0 ; j1 = j+j0
953  if (obc%segment(obc%segnum_v(i1,j1))%direction == obc_direction_n) then
954  hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dcv(i1,j1)
955  elseif (obc%segment(obc%segnum_v(i1,j1))%direction == obc_direction_s) then
956  hwt(i0,j0) = 2.0 * h(i1,j1+1,k) * mask2dcv(i1,j1)
957  endif
958  endif ; enddo ; enddo
959  endif ; endif
960 
961  hwt_tot = (hwt(0,-1) + hwt(1,0)) + (hwt(1,-1) + hwt(0,0))
962  set_v_at_u = 0.0
963  if (hwt_tot > 0.0) set_v_at_u = &
964  ((hwt(0,0) * v(i,j,k) + hwt(1,-1) * v(i+1,j-1,k)) + &
965  (hwt(1,0) * v(i+1,j,k) + hwt(0,-1) * v(i,j-1,k))) / hwt_tot
966 
967 end function set_v_at_u
968 
969 !> This subroutine finds a thickness-weighted value of u at the v-points.
970 function set_u_at_v(u, h, G, i, j, k, mask2dCu, OBC)
971  type(ocean_grid_type), intent(in) :: g !< The ocean's grid structure
972  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
973  intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1] or other units.
974  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
975  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
976  integer, intent(in) :: i !< The i-index of the u-location to work on.
977  integer, intent(in) :: j !< The j-index of the u-location to work on.
978  integer, intent(in) :: k !< The k-index of the u-location to work on.
979  real, dimension(SZIB_(G),SZJ_(G)), &
980  intent(in) :: mask2dcu !< A multiplicative mask of the u-points
981  type(ocean_obc_type), pointer :: obc !< A pointer to an open boundary condition structure
982  real :: set_u_at_v !< The return value of u at v points in the
983  !! same units as u, i.e. [L T-1 ~> m s-1] or other units.
984 
985  ! This subroutine finds a thickness-weighted value of u at the v-points.
986  real :: hwt(-1:0,0:1) ! Masked weights used to average u onto v [H ~> m or kg m-2].
987  real :: hwt_tot ! The sum of the masked thicknesses [H ~> m or kg m-2].
988  integer :: i0, j0, i1, j1
989 
990  do j0 = 0,1 ; do i0 = -1,0 ; i1 = i+i0 ; j1 = j+j0
991  hwt(i0,j0) = (h(i1,j1,k) + h(i1+1,j1,k)) * mask2dcu(i1,j1)
992  enddo ; enddo
993 
994  if (associated(obc)) then ; if (obc%number_of_segments > 0) then
995  do j0 = 0,1 ; do i0 = -1,0 ; if ((obc%segnum_u(i+i0,j+j0) /= obc_none)) then
996  i1 = i+i0 ; j1 = j+j0
997  if (obc%segment(obc%segnum_u(i1,j1))%direction == obc_direction_e) then
998  hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dcu(i1,j1)
999  elseif (obc%segment(obc%segnum_u(i1,j1))%direction == obc_direction_w) then
1000  hwt(i0,j0) = 2.0 * h(i1+1,j1,k) * mask2dcu(i1,j1)
1001  endif
1002  endif ; enddo ; enddo
1003  endif ; endif
1004 
1005  hwt_tot = (hwt(-1,0) + hwt(0,1)) + (hwt(0,0) + hwt(-1,1))
1006  set_u_at_v = 0.0
1007  if (hwt_tot > 0.0) set_u_at_v = &
1008  ((hwt(0,0) * u(i,j,k) + hwt(-1,1) * u(i-1,j+1,k)) + &
1009  (hwt(-1,0) * u(i-1,j,k) + hwt(0,1) * u(i,j+1,k))) / hwt_tot
1010 
1011 end function set_u_at_v
1012 
1013 !> Calculates the thickness of the surface boundary layer for applying an elevated viscosity.
1014 !!
1015 !! A bulk Richardson criterion or the thickness of the topmost NKML layers (with a bulk mixed layer)
1016 !! are currently used. The thicknesses are given in terms of fractional layers, so that this
1017 !! thickness will move as the thickness of the topmost layers change.
1018 subroutine set_viscous_ml(u, v, h, tv, forces, visc, dt, G, GV, US, CS, symmetrize)
1019  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
1020  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
1021  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1022  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1023  intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
1024  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1025  intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
1026  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1027  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
1028  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
1029  !! thermodynamic fields. Absent fields have
1030  !! NULL ptrs.
1031  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
1032  type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
1033  !! related fields.
1034  real, intent(in) :: dt !< Time increment [T ~> s].
1035  type(set_visc_cs), pointer :: cs !< The control structure returned by a previous
1036  !! call to vertvisc_init.
1037  logical, optional, intent(in) :: symmetrize !< If present and true, do extra calculations
1038  !! of those values in visc that would be
1039  !! calculated with symmetric memory.
1040 
1041  ! Local variables
1042  real, dimension(SZIB_(G)) :: &
1043  htot, & ! The total depth of the layers being that are within the
1044  ! surface mixed layer [H ~> m or kg m-2].
1045  thtot, & ! The integrated temperature of layers that are within the
1046  ! surface mixed layer [H degC ~> m degC or kg degC m-2].
1047  shtot, & ! The integrated salt of layers that are within the
1048  ! surface mixed layer [H ppt ~> m ppt or kg ppt m-2].
1049  rhtot, & ! The integrated density of layers that are within the surface mixed layer
1050  ! [H R ~> kg m-2 or kg2 m-5]. Rhtot is only used if no
1051  ! equation of state is used.
1052  uhtot, & ! The depth integrated zonal and meridional velocities within
1053  vhtot, & ! the surface mixed layer [H L T-1 ~> m2 s-1 or kg m-1 s-1].
1054  idecay_len_tke, & ! The inverse of a turbulence decay length scale [H-1 ~> m-1 or m2 kg-1].
1055  dr_dt, & ! Partial derivative of the density at the base of layer nkml
1056  ! (roughly the base of the mixed layer) with temperature [R degC-1 ~> kg m-3 degC-1].
1057  dr_ds, & ! Partial derivative of the density at the base of layer nkml
1058  ! (roughly the base of the mixed layer) with salinity [R ppt-1 ~> kg m-3 ppt-1].
1059  ustar, & ! The surface friction velocity under ice shelves [Z T-1 ~> m s-1].
1060  press, & ! The pressure at which dR_dT and dR_dS are evaluated [Pa].
1061  t_eos, & ! The potential temperature at which dR_dT and dR_dS are evaluated [degC]
1062  s_eos ! The salinity at which dR_dT and dR_dS are evaluated [ppt].
1063  real, dimension(SZIB_(G),SZJ_(G)) :: &
1064  mask_u ! A mask that disables any contributions from u points that
1065  ! are land or past open boundary conditions [nondim], 0 or 1.
1066  real, dimension(SZI_(G),SZJB_(G)) :: &
1067  mask_v ! A mask that disables any contributions from v points that
1068  ! are land or past open boundary conditions [nondim], 0 or 1.
1069  real :: h_at_vel(szib_(g),szk_(g))! Layer thickness at velocity points,
1070  ! using an upwind-biased second order accurate estimate based
1071  ! on the previous velocity direction [H ~> m or kg m-2].
1072  integer :: k_massive(szib_(g)) ! The k-index of the deepest layer yet found
1073  ! that has more than h_tiny thickness and will be in the
1074  ! viscous mixed layer.
1075  real :: uh2 ! The squared magnitude of the difference between the velocity
1076  ! integrated through the mixed layer and the velocity of the
1077  ! interior layer layer times the depth of the the mixed layer
1078  ! [H2 Z2 T-2 ~> m4 s-2 or kg2 m-2 s-2].
1079  real :: htot_vel ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
1080  real :: hwtot ! Sum of the thicknesses used to calculate
1081  ! the near-bottom velocity magnitude [H ~> m or kg m-2].
1082  real :: hutot ! Running sum of thicknesses times the
1083  ! velocity magnitudes [H L T-1 ~> m2 s-1 or kg m-1 s-1].
1084  real :: hweight ! The thickness of a layer that is within Hbbl
1085  ! of the bottom [H ~> m or kg m-2].
1086  real :: tbl_thick_z ! The thickness of the top boundary layer [Z ~> m].
1087 
1088  real :: hlay ! The layer thickness at velocity points [H ~> m or kg m-2].
1089  real :: i_2hlay ! 1 / 2*hlay [H-1 ~> m-1 or m2 kg-1].
1090  real :: t_lay ! The layer temperature at velocity points [degC].
1091  real :: s_lay ! The layer salinity at velocity points [ppt].
1092  real :: rlay ! The layer potential density at velocity points [R ~> kg m-3].
1093  real :: rlb ! The potential density of the layer below [R ~> kg m-3].
1094  real :: v_at_u ! The meridonal velocity at a zonal velocity point [L T-1 ~> m s-1].
1095  real :: u_at_v ! The zonal velocity at a meridonal velocity point [L T-1 ~> m s-1].
1096  real :: ghprime ! The mixed-layer internal gravity wave speed squared, based
1097  ! on the mixed layer thickness and density difference across
1098  ! the base of the mixed layer [L2 T-2 ~> m2 s-2].
1099  real :: ribulk ! The bulk Richardson number below which water is in the
1100  ! viscous mixed layer, including reduction for turbulent
1101  ! decay. Nondimensional.
1102  real :: dt_rho0 ! The time step divided by the conversion from the layer
1103  ! thickness to layer mass [T H Z-1 R-1 ~> s m3 kg-1 or s].
1104  real :: g_h_rho0 ! The gravitational acceleration times the conversion from H to m divided
1105  ! by the mean density [L2 T-2 H-1 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2].
1106  real :: ustarsq ! 400 times the square of ustar, times
1107  ! Rho0 divided by G_Earth and the conversion
1108  ! from m to thickness units [H R ~> kg m-2 or kg2 m-5].
1109  real :: cdrag_sqrt_z ! Square root of the drag coefficient, times a unit conversion
1110  ! factor from lateral lengths to vertical depths [Z L-1 ~> 1].
1111  real :: cdrag_sqrt ! Square root of the drag coefficient [nondim].
1112  real :: oldfn ! The integrated energy required to
1113  ! entrain up to the bottom of the layer,
1114  ! divided by G_Earth [H R ~> kg m-2 or kg2 m-5].
1115  real :: dfn ! The increment in oldfn for entraining
1116  ! the layer [H R ~> kg m-2 or kg2 m-5].
1117  real :: dh ! The increment in layer thickness from
1118  ! the present layer [H ~> m or kg m-2].
1119  real :: u_bg_sq ! The square of an assumed background velocity, for
1120  ! calculating the mean magnitude near the top for use in
1121  ! the quadratic surface drag [L2 T-2 ~> m2 s-2].
1122  real :: h_tiny ! A very small thickness [H ~> m or kg m-2]. Layers that are less than
1123  ! h_tiny can not be the deepest in the viscous mixed layer.
1124  real :: absf ! The absolute value of f averaged to velocity points [T-1 ~> s-1].
1125  real :: u_star ! The friction velocity at velocity points [Z T-1 ~> m s-1].
1126  real :: h_neglect ! A thickness that is so small it is usually lost
1127  ! in roundoff and can be neglected [H ~> m or kg m-2].
1128  real :: rho0x400_g ! 400*Rho0/G_Earth, times unit conversion factors
1129  ! [R T2 H Z-2 ~> kg s2 m-4 or kg2 s2 m-7].
1130  ! The 400 is a constant proposed by Killworth and Edwards, 1999.
1131  real :: ustar1 ! ustar [H T-1 ~> m s-1 or kg m-2 s-1]
1132  real :: h2f2 ! (h*2*f)^2 [H2 T-2 ~> m2 s-2 or kg2 m-4 s-2]
1133  logical :: use_eos, do_any, do_any_shelf, do_i(szib_(g))
1134  integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, k2, nkmb, nkml, n
1135  type(ocean_obc_type), pointer :: obc => null()
1136 
1137  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1138  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1139  nkmb = gv%nk_rho_varies ; nkml = gv%nkml
1140 
1141  if (.not.associated(cs)) call mom_error(fatal,"MOM_set_viscosity(visc_ML): "//&
1142  "Module must be initialized before it is used.")
1143  if (.not.(cs%dynamic_viscous_ML .or. associated(forces%frac_shelf_u) .or. &
1144  associated(forces%frac_shelf_v)) ) return
1145 
1146  if (present(symmetrize)) then ; if (symmetrize) then
1147  jsq = js-1 ; isq = is-1
1148  endif ; endif
1149 
1150  rho0x400_g = 400.0*(gv%Rho0/(us%L_to_Z**2 * gv%g_Earth)) * gv%Z_to_H
1151  u_bg_sq = cs%drag_bg_vel * cs%drag_bg_vel
1152  cdrag_sqrt = sqrt(cs%cdrag)
1153  cdrag_sqrt_z = us%L_to_Z * sqrt(cs%cdrag)
1154 
1155  obc => cs%OBC
1156  use_eos = associated(tv%eqn_of_state)
1157  dt_rho0 = dt / gv%H_to_RZ
1158  h_neglect = gv%H_subroundoff
1159  h_tiny = 2.0*gv%Angstrom_H + h_neglect
1160  g_h_rho0 = (gv%g_Earth*gv%H_to_Z) / (gv%Rho0)
1161 
1162  if (associated(forces%frac_shelf_u) .neqv. associated(forces%frac_shelf_v)) &
1163  call mom_error(fatal, "set_viscous_ML: one of forces%frac_shelf_u and "//&
1164  "forces%frac_shelf_v is associated, but the other is not.")
1165 
1166  if (associated(forces%frac_shelf_u)) then
1167  ! This configuration has ice shelves, and the appropriate variables need to be
1168  ! allocated. If the arrays have already been allocated, these calls do nothing.
1169  call safe_alloc_ptr(visc%tauy_shelf, g%isd, g%ied, g%JsdB, g%JedB)
1170  call safe_alloc_ptr(visc%tbl_thick_shelf_u, g%IsdB, g%IedB, g%jsd, g%jed)
1171  call safe_alloc_ptr(visc%tbl_thick_shelf_v, g%isd, g%ied, g%JsdB, g%JedB)
1172  call safe_alloc_ptr(visc%kv_tbl_shelf_u, g%IsdB, g%IedB, g%jsd, g%jed)
1173  call safe_alloc_ptr(visc%kv_tbl_shelf_v, g%isd, g%ied, g%JsdB, g%JedB)
1174  call safe_alloc_ptr(visc%taux_shelf, g%IsdB, g%IedB, g%jsd, g%jed)
1175  call safe_alloc_ptr(visc%tauy_shelf, g%isd, g%ied, g%JsdB, g%JedB)
1176 
1177  ! With a linear drag law under shelves, the friction velocity is already known.
1178 ! if (CS%linear_drag) ustar(:) = cdrag_sqrt_Z*CS%drag_bg_vel
1179  endif
1180 
1181  !$OMP parallel do default(shared)
1182  do j=js-1,je ; do i=is-1,ie+1
1183  mask_v(i,j) = g%mask2dCv(i,j)
1184  enddo ; enddo
1185  !$OMP parallel do default(shared)
1186  do j=js-1,je+1 ; do i=is-1,ie
1187  mask_u(i,j) = g%mask2dCu(i,j)
1188  enddo ; enddo
1189 
1190  if (associated(obc)) then ; do n=1,obc%number_of_segments
1191  ! Now project bottom depths across cell-corner points in the OBCs. The two
1192  ! projections have to occur in sequence and can not be combined easily.
1193  if (.not. obc%segment(n)%on_pe) cycle
1194  ! Use a one-sided projection of bottom depths at OBC points.
1195  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
1196  if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je)) then
1197  do i = max(is-1,obc%segment(n)%HI%IsdB), min(ie,obc%segment(n)%HI%IedB)
1198  if (obc%segment(n)%direction == obc_direction_n) mask_u(i,j+1) = 0.0
1199  if (obc%segment(n)%direction == obc_direction_s) mask_u(i,j) = 0.0
1200  enddo
1201  elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= je)) then
1202  do j = max(js-1,obc%segment(n)%HI%JsdB), min(je,obc%segment(n)%HI%JedB)
1203  if (obc%segment(n)%direction == obc_direction_e) mask_v(i+1,j) = 0.0
1204  if (obc%segment(n)%direction == obc_direction_w) mask_v(i,j) = 0.0
1205  enddo
1206  endif
1207  enddo ; endif
1208 
1209  !$OMP parallel do default(private) shared(u,v,h,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, &
1210  !$OMP h_neglect,h_tiny,g_H_Rho0,js,je,OBC,Isq,Ieq,nz, &
1211  !$OMP U_bg_sq,mask_v,cdrag_sqrt,cdrag_sqrt_Z,Rho0x400_G,nkml)
1212  do j=js,je ! u-point loop
1213  if (cs%dynamic_viscous_ML) then
1214  do_any = .false.
1215  do i=isq,ieq
1216  htot(i) = 0.0
1217  if (g%mask2dCu(i,j) < 0.5) then
1218  do_i(i) = .false. ; visc%nkml_visc_u(i,j) = nkml
1219  else
1220  do_i(i) = .true. ; do_any = .true.
1221  k_massive(i) = nkml
1222  thtot(i) = 0.0 ; shtot(i) = 0.0 ; rhtot(i) = 0.0
1223  uhtot(i) = dt_rho0 * forces%taux(i,j)
1224  vhtot(i) = 0.25 * dt_rho0 * ((forces%tauy(i,j) + forces%tauy(i+1,j-1)) + &
1225  (forces%tauy(i,j-1) + forces%tauy(i+1,j)))
1226 
1227  if (cs%omega_frac >= 1.0) then ; absf = 2.0*cs%omega ; else
1228  absf = 0.5*(abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i,j-1)))
1229  if (cs%omega_frac > 0.0) &
1230  absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
1231  endif
1232  u_star = max(cs%ustar_min, 0.5 * (forces%ustar(i,j) + forces%ustar(i+1,j)))
1233  idecay_len_tke(i) = ((absf / u_star) * cs%TKE_decay) * gv%H_to_Z
1234  endif
1235  enddo
1236 
1237  if (do_any) then ; do k=1,nz
1238  if (k > nkml) then
1239  do_any = .false.
1240  if (use_eos .and. (k==nkml+1)) then
1241  ! Find dRho/dT and dRho_dS.
1242  do i=isq,ieq
1243  press(i) = gv%H_to_Pa * htot(i)
1244  k2 = max(1,nkml)
1245  i_2hlay = 1.0 / (h(i,j,k2) + h(i+1,j,k2) + h_neglect)
1246  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
1247  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
1248  enddo
1249  call calculate_density_derivs(t_eos, s_eos, press, dr_dt, dr_ds, &
1250  isq-g%IsdB+1, ieq-isq+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
1251  endif
1252 
1253  do i=isq,ieq ; if (do_i(i)) then
1254 
1255  hlay = 0.5*(h(i,j,k) + h(i+1,j,k))
1256  if (hlay > h_tiny) then ! Only consider non-vanished layers.
1257  i_2hlay = 1.0 / (h(i,j,k) + h(i+1,j,k))
1258  v_at_u = 0.5 * (h(i,j,k) * (v(i,j,k) + v(i,j-1,k)) + &
1259  h(i+1,j,k) * (v(i+1,j,k) + v(i+1,j-1,k))) * i_2hlay
1260  uh2 = ((uhtot(i) - htot(i)*u(i,j,k))**2 + (vhtot(i) - htot(i)*v_at_u)**2)
1261 
1262  if (use_eos) then
1263  t_lay = (h(i,j,k)*tv%T(i,j,k) + h(i+1,j,k)*tv%T(i+1,j,k)) * i_2hlay
1264  s_lay = (h(i,j,k)*tv%S(i,j,k) + h(i+1,j,k)*tv%S(i+1,j,k)) * i_2hlay
1265  ghprime = g_h_rho0 * (dr_dt(i) * (t_lay*htot(i) - thtot(i)) + &
1266  dr_ds(i) * (s_lay*htot(i) - shtot(i)))
1267  else
1268  ghprime = g_h_rho0 * (gv%Rlay(k)*htot(i) - rhtot(i))
1269  endif
1270 
1271  if (ghprime > 0.0) then
1272  ribulk = cs%bulk_Ri_ML * exp(-htot(i) * idecay_len_tke(i))
1273  if (ribulk * uh2 <= (htot(i)**2) * ghprime) then
1274  visc%nkml_visc_u(i,j) = real(k_massive(i))
1275  do_i(i) = .false.
1276  elseif (ribulk * uh2 <= (htot(i) + hlay)**2 * ghprime) then
1277  visc%nkml_visc_u(i,j) = real(k-1) + &
1278  ( sqrt(ribulk * uh2 / ghprime) - htot(i) ) / hlay
1279  do_i(i) = .false.
1280  endif
1281  endif
1282  k_massive(i) = k
1283  endif ! hlay > h_tiny
1284 
1285  if (do_i(i)) do_any = .true.
1286  endif ; enddo
1287 
1288  if (.not.do_any) exit ! All columns are done.
1289  endif
1290 
1291  do i=isq,ieq ; if (do_i(i)) then
1292  htot(i) = htot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k))
1293  uhtot(i) = uhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * u(i,j,k)
1294  vhtot(i) = vhtot(i) + 0.25 * (h(i,j,k) * (v(i,j,k) + v(i,j-1,k)) + &
1295  h(i+1,j,k) * (v(i+1,j,k) + v(i+1,j-1,k)))
1296  if (use_eos) then
1297  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))
1298  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))
1299  else
1300  rhtot(i) = rhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * gv%Rlay(k)
1301  endif
1302  endif ; enddo
1303  enddo ; endif
1304 
1305  if (do_any) then ; do i=isq,ieq ; if (do_i(i)) then
1306  visc%nkml_visc_u(i,j) = k_massive(i)
1307  endif ; enddo ; endif
1308  endif ! dynamic_viscous_ML
1309 
1310  do_any_shelf = .false.
1311  if (associated(forces%frac_shelf_u)) then
1312  do i=isq,ieq
1313  if (forces%frac_shelf_u(i,j)*g%mask2dCu(i,j) == 0.0) then
1314  do_i(i) = .false.
1315  visc%tbl_thick_shelf_u(i,j) = 0.0 ; visc%kv_tbl_shelf_u(i,j) = 0.0
1316  else
1317  do_i(i) = .true. ; do_any_shelf = .true.
1318  endif
1319  enddo
1320  endif
1321 
1322  if (do_any_shelf) then
1323  do k=1,nz ; do i=isq,ieq ; if (do_i(i)) then
1324  if (u(i,j,k) *(h(i+1,j,k) - h(i,j,k)) >= 0) then
1325  h_at_vel(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / &
1326  (h(i,j,k) + h(i+1,j,k) + h_neglect)
1327  else
1328  h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
1329  endif
1330  else
1331  h_at_vel(i,k) = 0.0 ; ustar(i) = 0.0
1332  endif ; enddo ; enddo
1333 
1334  do i=isq,ieq ; if (do_i(i)) then
1335  htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
1336  thtot(i) = 0.0 ; shtot(i) = 0.0
1337  if (use_eos .or. .not.cs%linear_drag) then ; do k=1,nz
1338  if (htot_vel>=cs%Htbl_shelf) exit ! terminate the k loop
1339  hweight = min(cs%Htbl_shelf - htot_vel, h_at_vel(i,k))
1340  if (hweight <= 1.5*gv%Angstrom_H + h_neglect) cycle
1341 
1342  htot_vel = htot_vel + h_at_vel(i,k)
1343  hwtot = hwtot + hweight
1344 
1345  if (.not.cs%linear_drag) then
1346  v_at_u = set_v_at_u(v, h, g, i, j, k, mask_v, obc)
1347  hutot = hutot + hweight * sqrt(u(i,j,k)**2 + &
1348  v_at_u**2 + u_bg_sq)
1349  endif
1350  if (use_eos) then
1351  thtot(i) = thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
1352  shtot(i) = shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
1353  endif
1354  enddo ; endif
1355 
1356  if ((.not.cs%linear_drag) .and. (hwtot > 0.0)) then
1357  ustar(i) = cdrag_sqrt_z * hutot/hwtot
1358  else
1359  ustar(i) = cdrag_sqrt_z * cs%drag_bg_vel
1360  endif
1361 
1362  if (use_eos) then ; if (hwtot > 0.0) then
1363  t_eos(i) = thtot(i)/hwtot ; s_eos(i) = shtot(i)/hwtot
1364  else
1365  t_eos(i) = 0.0 ; s_eos(i) = 0.0
1366  endif ; endif
1367  endif ; enddo ! I-loop
1368 
1369  if (use_eos) then
1370  call calculate_density_derivs(t_eos, s_eos, forces%p_surf(:,j), &
1371  dr_dt, dr_ds, isq-g%IsdB+1, ieq-isq+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
1372  endif
1373 
1374  do i=isq,ieq ; if (do_i(i)) then
1375  ! The 400.0 in this expression is the square of a constant proposed
1376  ! by Killworth and Edwards, 1999, in equation (2.20).
1377  ustarsq = rho0x400_g * ustar(i)**2
1378  htot(i) = 0.0
1379  if (use_eos) then
1380  thtot(i) = 0.0 ; shtot(i) = 0.0
1381  do k=1,nz-1
1382  if (h_at_vel(i,k) <= 0.0) cycle
1383  t_lay = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
1384  s_lay = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
1385  oldfn = dr_dt(i)*(t_lay*htot(i) - thtot(i)) + dr_ds(i)*(s_lay*htot(i) - shtot(i))
1386  if (oldfn >= ustarsq) exit
1387 
1388  dfn = (dr_dt(i)*(0.5*(tv%T(i,j,k+1)+tv%T(i+1,j,k+1)) - t_lay) + &
1389  dr_ds(i)*(0.5*(tv%S(i,j,k+1)+tv%S(i+1,j,k+1)) - s_lay)) * &
1390  (h_at_vel(i,k)+htot(i))
1391  if ((oldfn + dfn) <= ustarsq) then
1392  dh = h_at_vel(i,k)
1393  else
1394  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
1395  endif
1396 
1397  htot(i) = htot(i) + dh
1398  thtot(i) = thtot(i) + t_lay*dh ; shtot(i) = shtot(i) + s_lay*dh
1399  enddo
1400  if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0)) then
1401  t_lay = 0.5*(tv%T(i,j,nz) + tv%T(i+1,j,nz))
1402  s_lay = 0.5*(tv%S(i,j,nz) + tv%S(i+1,j,nz))
1403  if (dr_dt(i)*(t_lay*htot(i) - thtot(i)) + &
1404  dr_ds(i)*(s_lay*htot(i) - shtot(i)) < ustarsq) &
1405  htot(i) = htot(i) + h_at_vel(i,nz)
1406  endif ! Examination of layer nz.
1407  else ! Use Rlay as the density variable.
1408  rhtot = 0.0
1409  do k=1,nz-1
1410  rlay = gv%Rlay(k) ; rlb = gv%Rlay(k+1)
1411 
1412  oldfn = rlay*htot(i) - rhtot(i)
1413  if (oldfn >= ustarsq) exit
1414 
1415  dfn = (rlb - rlay)*(h_at_vel(i,k)+htot(i))
1416  if ((oldfn + dfn) <= ustarsq) then
1417  dh = h_at_vel(i,k)
1418  else
1419  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
1420  endif
1421 
1422  htot(i) = htot(i) + dh
1423  rhtot(i) = rhtot(i) + rlay*dh
1424  enddo
1425  if (gv%Rlay(nz)*htot(i) - rhtot(i) < ustarsq) &
1426  htot(i) = htot(i) + h_at_vel(i,nz)
1427  endif ! use_EOS
1428 
1429  !visc%tbl_thick_shelf_u(I,j) = GV%H_to_Z * max(CS%Htbl_shelf_min, &
1430  ! htot(I) / (0.5 + sqrt(0.25 + &
1431  ! (htot(i)*(G%CoriolisBu(I,J-1)+G%CoriolisBu(I,J)))**2 / &
1432  ! (ustar(i)*GV%Z_to_H)**2 )) )
1433  ustar1 = ustar(i)*gv%Z_to_H
1434  h2f2 = (htot(i)*(g%CoriolisBu(i,j-1)+g%CoriolisBu(i,j)) + h_neglect*cs%omega)**2
1435  tbl_thick_z = gv%H_to_Z * max(cs%Htbl_shelf_min, &
1436  ( htot(i)*ustar1 ) / ( 0.5*ustar1 + sqrt((0.5*ustar1)**2 + h2f2 ) ) )
1437  visc%tbl_thick_shelf_u(i,j) = tbl_thick_z
1438  visc%Kv_tbl_shelf_u(i,j) = max(cs%Kv_TBL_min, cdrag_sqrt*ustar(i)*tbl_thick_z)
1439  endif ; enddo ! I-loop
1440  endif ! do_any_shelf
1441 
1442  enddo ! j-loop at u-points
1443 
1444  !$OMP parallel do default(private) shared(u,v,h,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, &
1445  !$OMP h_neglect,h_tiny,g_H_Rho0,is,ie,OBC,Jsq,Jeq,nz, &
1446  !$OMP U_bg_sq,cdrag_sqrt,cdrag_sqrt_Z,Rho0x400_G,nkml,mask_u)
1447  do j=jsq,jeq ! v-point loop
1448  if (cs%dynamic_viscous_ML) then
1449  do_any = .false.
1450  do i=is,ie
1451  htot(i) = 0.0
1452  if (g%mask2dCv(i,j) < 0.5) then
1453  do_i(i) = .false. ; visc%nkml_visc_v(i,j) = nkml
1454  else
1455  do_i(i) = .true. ; do_any = .true.
1456  k_massive(i) = nkml
1457  thtot(i) = 0.0 ; shtot(i) = 0.0 ; rhtot(i) = 0.0
1458  vhtot(i) = dt_rho0 * forces%tauy(i,j)
1459  uhtot(i) = 0.25 * dt_rho0 * ((forces%taux(i,j) + forces%taux(i-1,j+1)) + &
1460  (forces%taux(i-1,j) + forces%taux(i,j+1)))
1461 
1462  if (cs%omega_frac >= 1.0) then ; absf = 2.0*cs%omega ; else
1463  absf = 0.5*(abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))
1464  if (cs%omega_frac > 0.0) &
1465  absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
1466  endif
1467 
1468  u_star = max(cs%ustar_min, 0.5 * (forces%ustar(i,j) + forces%ustar(i,j+1)))
1469  idecay_len_tke(i) = ((absf / u_star) * cs%TKE_decay) * gv%H_to_Z
1470 
1471  endif
1472  enddo
1473 
1474  if (do_any) then ; do k=1,nz
1475  if (k > nkml) then
1476  do_any = .false.
1477  if (use_eos .and. (k==nkml+1)) then
1478  ! Find dRho/dT and dRho_dS.
1479  do i=is,ie
1480  press(i) = gv%H_to_Pa * htot(i)
1481  k2 = max(1,nkml)
1482  i_2hlay = 1.0 / (h(i,j,k2) + h(i,j+1,k2) + h_neglect)
1483  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
1484  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
1485  enddo
1486  call calculate_density_derivs(t_eos, s_eos, press, dr_dt, dr_ds, &
1487  is-g%IsdB+1, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
1488  endif
1489 
1490  do i=is,ie ; if (do_i(i)) then
1491 
1492  hlay = 0.5*(h(i,j,k) + h(i,j+1,k))
1493  if (hlay > h_tiny) then ! Only consider non-vanished layers.
1494  i_2hlay = 1.0 / (h(i,j,k) + h(i,j+1,k))
1495  u_at_v = 0.5 * (h(i,j,k) * (u(i-1,j,k) + u(i,j,k)) + &
1496  h(i,j+1,k) * (u(i-1,j+1,k) + u(i,j+1,k))) * i_2hlay
1497  uh2 = ((uhtot(i) - htot(i)*u_at_v)**2 + (vhtot(i) - htot(i)*v(i,j,k))**2)
1498 
1499  if (use_eos) then
1500  t_lay = (h(i,j,k)*tv%T(i,j,k) + h(i,j+1,k)*tv%T(i,j+1,k)) * i_2hlay
1501  s_lay = (h(i,j,k)*tv%S(i,j,k) + h(i,j+1,k)*tv%S(i,j+1,k)) * i_2hlay
1502  ghprime = g_h_rho0 * (dr_dt(i) * (t_lay*htot(i) - thtot(i)) + &
1503  dr_ds(i) * (s_lay*htot(i) - shtot(i)))
1504  else
1505  ghprime = g_h_rho0 * (gv%Rlay(k)*htot(i) - rhtot(i))
1506  endif
1507 
1508  if (ghprime > 0.0) then
1509  ribulk = cs%bulk_Ri_ML * exp(-htot(i) * idecay_len_tke(i))
1510  if (ribulk * uh2 <= htot(i)**2 * ghprime) then
1511  visc%nkml_visc_v(i,j) = real(k_massive(i))
1512  do_i(i) = .false.
1513  elseif (ribulk * uh2 <= (htot(i) + hlay)**2 * ghprime) then
1514  visc%nkml_visc_v(i,j) = real(k-1) + &
1515  ( sqrt(ribulk * uh2 / ghprime) - htot(i) ) / hlay
1516  do_i(i) = .false.
1517  endif
1518  endif
1519  k_massive(i) = k
1520  endif ! hlay > h_tiny
1521 
1522  if (do_i(i)) do_any = .true.
1523  endif ; enddo
1524 
1525  if (.not.do_any) exit ! All columns are done.
1526  endif
1527 
1528  do i=is,ie ; if (do_i(i)) then
1529  htot(i) = htot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k))
1530  vhtot(i) = vhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * v(i,j,k)
1531  uhtot(i) = uhtot(i) + 0.25 * (h(i,j,k) * (u(i-1,j,k) + u(i,j,k)) + &
1532  h(i,j+1,k) * (u(i-1,j+1,k) + u(i,j+1,k)))
1533  if (use_eos) then
1534  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))
1535  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))
1536  else
1537  rhtot(i) = rhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * gv%Rlay(k)
1538  endif
1539  endif ; enddo
1540  enddo ; endif
1541 
1542  if (do_any) then ; do i=is,ie ; if (do_i(i)) then
1543  visc%nkml_visc_v(i,j) = k_massive(i)
1544  endif ; enddo ; endif
1545  endif ! dynamic_viscous_ML
1546 
1547  do_any_shelf = .false.
1548  if (associated(forces%frac_shelf_v)) then
1549  do i=is,ie
1550  if (forces%frac_shelf_v(i,j)*g%mask2dCv(i,j) == 0.0) then
1551  do_i(i) = .false.
1552  visc%tbl_thick_shelf_v(i,j) = 0.0 ; visc%kv_tbl_shelf_v(i,j) = 0.0
1553  else
1554  do_i(i) = .true. ; do_any_shelf = .true.
1555  endif
1556  enddo
1557  endif
1558 
1559  if (do_any_shelf) then
1560  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
1561  if (v(i,j,k) * (h(i,j+1,k) - h(i,j,k)) >= 0) then
1562  h_at_vel(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / &
1563  (h(i,j,k) + h(i,j+1,k) + h_neglect)
1564  else
1565  h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
1566  endif
1567  else
1568  h_at_vel(i,k) = 0.0 ; ustar(i) = 0.0
1569  endif ; enddo ; enddo
1570 
1571  do i=is,ie ; if (do_i(i)) then
1572  htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
1573  thtot(i) = 0.0 ; shtot(i) = 0.0
1574  if (use_eos .or. .not.cs%linear_drag) then ; do k=1,nz
1575  if (htot_vel>=cs%Htbl_shelf) exit ! terminate the k loop
1576  hweight = min(cs%Htbl_shelf - htot_vel, h_at_vel(i,k))
1577  if (hweight <= 1.5*gv%Angstrom_H + h_neglect) cycle
1578 
1579  htot_vel = htot_vel + h_at_vel(i,k)
1580  hwtot = hwtot + hweight
1581 
1582  if (.not.cs%linear_drag) then
1583  u_at_v = set_u_at_v(u, h, g, i, j, k, mask_u, obc)
1584  hutot = hutot + hweight * sqrt(v(i,j,k)**2 + &
1585  u_at_v**2 + u_bg_sq)
1586  endif
1587  if (use_eos) then
1588  thtot(i) = thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
1589  shtot(i) = shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
1590  endif
1591  enddo ; endif
1592 
1593  if (.not.cs%linear_drag) then ; if (hwtot > 0.0) then
1594  ustar(i) = cdrag_sqrt_z * hutot/hwtot
1595  else
1596  ustar(i) = cdrag_sqrt_z * cs%drag_bg_vel
1597  endif ; endif
1598 
1599  if (use_eos) then ; if (hwtot > 0.0) then
1600  t_eos(i) = thtot(i)/hwtot ; s_eos(i) = shtot(i)/hwtot
1601  else
1602  t_eos(i) = 0.0 ; s_eos(i) = 0.0
1603  endif ; endif
1604  endif ; enddo ! I-loop
1605 
1606  if (use_eos) then
1607  call calculate_density_derivs(t_eos, s_eos, forces%p_surf(:,j), &
1608  dr_dt, dr_ds, is-g%IsdB+1, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
1609  endif
1610 
1611  do i=is,ie ; if (do_i(i)) then
1612  ! The 400.0 in this expression is the square of a constant proposed
1613  ! by Killworth and Edwards, 1999, in equation (2.20).
1614  ustarsq = rho0x400_g * ustar(i)**2
1615  htot(i) = 0.0
1616  if (use_eos) then
1617  thtot(i) = 0.0 ; shtot(i) = 0.0
1618  do k=1,nz-1
1619  if (h_at_vel(i,k) <= 0.0) cycle
1620  t_lay = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
1621  s_lay = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
1622  oldfn = dr_dt(i)*(t_lay*htot(i) - thtot(i)) + dr_ds(i)*(s_lay*htot(i) - shtot(i))
1623  if (oldfn >= ustarsq) exit
1624 
1625  dfn = (dr_dt(i)*(0.5*(tv%T(i,j,k+1)+tv%T(i,j+1,k+1)) - t_lay) + &
1626  dr_ds(i)*(0.5*(tv%S(i,j,k+1)+tv%S(i,j+1,k+1)) - s_lay)) * &
1627  (h_at_vel(i,k)+htot(i))
1628  if ((oldfn + dfn) <= ustarsq) then
1629  dh = h_at_vel(i,k)
1630  else
1631  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
1632  endif
1633 
1634  htot(i) = htot(i) + dh
1635  thtot(i) = thtot(i) + t_lay*dh ; shtot(i) = shtot(i) + s_lay*dh
1636  enddo
1637  if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0)) then
1638  t_lay = 0.5*(tv%T(i,j,nz) + tv%T(i,j+1,nz))
1639  s_lay = 0.5*(tv%S(i,j,nz) + tv%S(i,j+1,nz))
1640  if (dr_dt(i)*(t_lay*htot(i) - thtot(i)) + &
1641  dr_ds(i)*(s_lay*htot(i) - shtot(i)) < ustarsq) &
1642  htot(i) = htot(i) + h_at_vel(i,nz)
1643  endif ! Examination of layer nz.
1644  else ! Use Rlay as the density variable.
1645  rhtot = 0.0
1646  do k=1,nz-1
1647  rlay = gv%Rlay(k) ; rlb = gv%Rlay(k+1)
1648 
1649  oldfn = rlay*htot(i) - rhtot(i)
1650  if (oldfn >= ustarsq) exit
1651 
1652  dfn = (rlb - rlay)*(h_at_vel(i,k)+htot(i))
1653  if ((oldfn + dfn) <= ustarsq) then
1654  dh = h_at_vel(i,k)
1655  else
1656  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
1657  endif
1658 
1659  htot(i) = htot(i) + dh
1660  rhtot = rhtot + rlay*dh
1661  enddo
1662  if (gv%Rlay(nz)*htot(i) - rhtot(i) < ustarsq) &
1663  htot(i) = htot(i) + h_at_vel(i,nz)
1664  endif ! use_EOS
1665 
1666  !visc%tbl_thick_shelf_v(i,J) = GV%H_to_Z * max(CS%Htbl_shelf_min, &
1667  ! htot(i) / (0.5 + sqrt(0.25 + &
1668  ! (htot(i)*(G%CoriolisBu(I-1,J)+G%CoriolisBu(I,J)))**2 / &
1669  ! (ustar(i)*GV%Z_to_H)**2 )) )
1670  ustar1 = ustar(i)*gv%Z_to_H
1671  h2f2 = (htot(i)*(g%CoriolisBu(i-1,j)+g%CoriolisBu(i,j)) + h_neglect*cs%omega)**2
1672  tbl_thick_z = gv%H_to_Z * max(cs%Htbl_shelf_min, &
1673  ( htot(i)*ustar1 ) / ( 0.5*ustar1 + sqrt((0.5*ustar1)**2 + h2f2 ) ) )
1674  visc%tbl_thick_shelf_v(i,j) = tbl_thick_z
1675  visc%Kv_tbl_shelf_v(i,j) = max(cs%Kv_TBL_min, cdrag_sqrt*ustar(i)*tbl_thick_z)
1676 
1677  endif ; enddo ! i-loop
1678  endif ! do_any_shelf
1679 
1680  enddo ! J-loop at v-points
1681 
1682  if (cs%debug) then
1683  if (associated(visc%nkml_visc_u) .and. associated(visc%nkml_visc_v)) &
1684  call uvchksum("nkml_visc_[uv]", visc%nkml_visc_u, &
1685  visc%nkml_visc_v, g%HI,haloshift=0)
1686  endif
1687  if (cs%id_nkml_visc_u > 0) &
1688  call post_data(cs%id_nkml_visc_u, visc%nkml_visc_u, cs%diag)
1689  if (cs%id_nkml_visc_v > 0) &
1690  call post_data(cs%id_nkml_visc_v, visc%nkml_visc_v, cs%diag)
1691 
1692 end subroutine set_viscous_ml
1693 
1694 !> Register any fields associated with the vertvisc_type.
1695 subroutine set_visc_register_restarts(HI, GV, param_file, visc, restart_CS)
1696  type(hor_index_type), intent(in) :: hi !< A horizontal index type structure.
1697  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
1698  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1699  !! parameters.
1700  type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical
1701  !! viscosities and related fields.
1702  !! Allocated here.
1703  type(mom_restart_cs), pointer :: restart_cs !< A pointer to the restart control structure.
1704  ! Local variables
1705  logical :: use_kappa_shear, ks_at_vertex
1706  logical :: adiabatic, usekpp, useepbl
1707  logical :: use_cvmix_shear, mle_use_pbl_mld, use_cvmix_conv
1708  integer :: isd, ied, jsd, jed, nz
1709  real :: hfreeze !< If hfreeze > 0 [m], melt potential will be computed.
1710  character(len=40) :: mdl = "MOM_set_visc" ! This module's name.
1711  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
1712 
1713  call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
1714  do_not_log=.true.)
1715 
1716  use_kappa_shear = .false. ; ks_at_vertex = .false. ; use_cvmix_shear = .false.
1717  usekpp = .false. ; useepbl = .false. ; use_cvmix_conv = .false.
1718 
1719  if (.not.adiabatic) then
1720  use_kappa_shear = kappa_shear_is_used(param_file)
1721  ks_at_vertex = kappa_shear_at_vertex(param_file)
1722  use_cvmix_shear = cvmix_shear_is_used(param_file)
1723  use_cvmix_conv = cvmix_conv_is_used(param_file)
1724  call get_param(param_file, mdl, "USE_KPP", usekpp, &
1725  "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
1726  "to calculate diffusivities and non-local transport in the OBL.", &
1727  default=.false., do_not_log=.true.)
1728  call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", useepbl, &
1729  "If true, use an implied energetics planetary boundary "//&
1730  "layer scheme to determine the diffusivity and viscosity "//&
1731  "in the surface boundary layer.", default=.false., do_not_log=.true.)
1732  endif
1733 
1734  if (use_kappa_shear .or. usekpp .or. useepbl .or. use_cvmix_shear .or. use_cvmix_conv) then
1735  call safe_alloc_ptr(visc%Kd_shear, isd, ied, jsd, jed, nz+1)
1736  call register_restart_field(visc%Kd_shear, "Kd_shear", .false., restart_cs, &
1737  "Shear-driven turbulent diffusivity at interfaces", "m2 s-1", z_grid='i')
1738  endif
1739  if (usekpp .or. useepbl .or. use_cvmix_shear .or. use_cvmix_conv .or. &
1740  (use_kappa_shear .and. .not.ks_at_vertex )) then
1741  call safe_alloc_ptr(visc%Kv_shear, isd, ied, jsd, jed, nz+1)
1742  call register_restart_field(visc%Kv_shear, "Kv_shear", .false., restart_cs, &
1743  "Shear-driven turbulent viscosity at interfaces", "m2 s-1", z_grid='i')
1744  endif
1745  if (use_kappa_shear .and. ks_at_vertex) then
1746  call safe_alloc_ptr(visc%TKE_turb, hi%IsdB, hi%IedB, hi%JsdB, hi%JedB, nz+1)
1747  call safe_alloc_ptr(visc%Kv_shear_Bu, hi%IsdB, hi%IedB, hi%JsdB, hi%JedB, nz+1)
1748  call register_restart_field(visc%Kv_shear_Bu, "Kv_shear_Bu", .false., restart_cs, &
1749  "Shear-driven turbulent viscosity at vertex interfaces", "m2 s-1", &
1750  hor_grid="Bu", z_grid='i')
1751  elseif (use_kappa_shear) then
1752  call safe_alloc_ptr(visc%TKE_turb, isd, ied, jsd, jed, nz+1)
1753  endif
1754 
1755  if (usekpp) then
1756  ! MOM_bkgnd_mixing uses Kv_slow when KPP is defined.
1757  call safe_alloc_ptr(visc%Kv_slow, isd, ied, jsd, jed, nz+1)
1758  endif
1759 
1760  ! visc%MLD is used to communicate the state of the (e)PBL or KPP to the rest of the model
1761  call get_param(param_file, mdl, "MLE_USE_PBL_MLD", mle_use_pbl_mld, &
1762  default=.false., do_not_log=.true.)
1763  ! visc%MLD needs to be allocated when melt potential is computed (HFREEZE>0)
1764  call get_param(param_file, mdl, "HFREEZE", hfreeze, &
1765  default=-1.0, do_not_log=.true.)
1766 
1767  if (mle_use_pbl_mld) then
1768  call safe_alloc_ptr(visc%MLD, isd, ied, jsd, jed)
1769  call register_restart_field(visc%MLD, "MLD", .false., restart_cs, &
1770  "Instantaneous active mixing layer depth", "m")
1771  endif
1772 
1773  if (hfreeze >= 0.0 .and. .not.mle_use_pbl_mld) then
1774  call safe_alloc_ptr(visc%MLD, isd, ied, jsd, jed)
1775  endif
1776 
1777 
1778 end subroutine set_visc_register_restarts
1779 
1780 !> Initializes the MOM_set_visc control structure
1781 subroutine set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS, OBC)
1782  type(time_type), target, intent(in) :: time !< The current model time.
1783  type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
1784  type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
1785  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
1786  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1787  !! parameters.
1788  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
1789  !! output.
1790  type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
1791  !! related fields. Allocated here.
1792  type(set_visc_cs), pointer :: cs !< A pointer that is set to point to the control
1793  !! structure for this module
1794  type(mom_restart_cs), pointer :: restart_cs !< A pointer to the restart control structure.
1795  type(ocean_obc_type), pointer :: obc !< A pointer to an open boundary condition structure
1796 
1797  ! Local variables
1798  real :: csmag_chan_dflt, smag_const1, tke_decay_dflt, bulk_ri_ml_dflt
1799  real :: kv_background
1800  real :: omega_frac_dflt
1801  real :: z_rescale ! A rescaling factor for heights from the representation in
1802  ! a restart file to the internal representation in this run.
1803  real :: i_t_rescale ! A rescaling factor for time from the internal representation in this run
1804  ! to the representation in a restart file.
1805  real :: z2_t_rescale ! A rescaling factor for vertical diffusivities and viscosities from the
1806  ! representation in a restart file to the internal representation in this run.
1807  integer :: i, j, k, is, ie, js, je, n
1808  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb, nz
1809  logical :: default_2018_answers
1810  logical :: use_kappa_shear, adiabatic, use_omega
1811  logical :: use_cvmix_ddiff, differential_diffusion, use_kpp
1812  type(obc_segment_type), pointer :: segment => null() ! pointer to OBC segment type
1813  ! This include declares and sets the variable "version".
1814 # include "version_variable.h"
1815  character(len=40) :: mdl = "MOM_set_visc" ! This module's name.
1816 
1817  if (associated(cs)) then
1818  call mom_error(warning, "set_visc_init called with an associated "// &
1819  "control structure.")
1820  return
1821  endif
1822  allocate(cs)
1823 
1824  cs%OBC => obc
1825 
1826  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1827  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
1828  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1829 
1830  cs%diag => diag
1831 
1832  ! Set default, read and log parameters
1833  call log_version(param_file, mdl, version, "")
1834  cs%RiNo_mix = .false. ; use_cvmix_ddiff = .false.
1835  differential_diffusion = .false.
1836  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
1837  "This sets the default value for the various _2018_ANSWERS parameters.", &
1838  default=.true.)
1839  call get_param(param_file, mdl, "SET_VISC_2018_ANSWERS", cs%answers_2018, &
1840  "If true, use the order of arithmetic and expressions that recover the "//&
1841  "answers from the end of 2018. Otherwise, use updated and more robust "//&
1842  "forms of the same expressions.", default=default_2018_answers)
1843  call get_param(param_file, mdl, "BOTTOMDRAGLAW", cs%bottomdraglaw, &
1844  "If true, the bottom stress is calculated with a drag "//&
1845  "law of the form c_drag*|u|*u. The velocity magnitude "//&
1846  "may be an assumed value or it may be based on the "//&
1847  "actual velocity in the bottommost HBBL, depending on "//&
1848  "LINEAR_DRAG.", default=.true.)
1849  call get_param(param_file, mdl, "CHANNEL_DRAG", cs%Channel_drag, &
1850  "If true, the bottom drag is exerted directly on each "//&
1851  "layer proportional to the fraction of the bottom it "//&
1852  "overlies.", default=.false.)
1853  call get_param(param_file, mdl, "LINEAR_DRAG", cs%linear_drag, &
1854  "If LINEAR_DRAG and BOTTOMDRAGLAW are defined the drag "//&
1855  "law is cdrag*DRAG_BG_VEL*u.", default=.false.)
1856  call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
1857  do_not_log=.true.)
1858  if (adiabatic) then
1859  call log_param(param_file, mdl, "ADIABATIC",adiabatic, &
1860  "There are no diapycnal mass fluxes if ADIABATIC is "//&
1861  "true. This assumes that KD = KDML = 0.0 and that "//&
1862  "there is no buoyancy forcing, but makes the model "//&
1863  "faster by eliminating subroutine calls.", default=.false.)
1864  endif
1865 
1866  if (.not.adiabatic) then
1867  cs%RiNo_mix = kappa_shear_is_used(param_file)
1868  call get_param(param_file, mdl, "DOUBLE_DIFFUSION", differential_diffusion, &
1869  "If true, increase diffusivites for temperature or salt "//&
1870  "based on double-diffusive parameterization from MOM4/KPP.", &
1871  default=.false.)
1872  use_cvmix_ddiff = cvmix_ddiff_is_used(param_file)
1873  endif
1874 
1875  call get_param(param_file, mdl, "PRANDTL_TURB", visc%Prandtl_turb, &
1876  "The turbulent Prandtl number applied to shear "//&
1877  "instability.", units="nondim", default=1.0)
1878  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
1879 
1880  call get_param(param_file, mdl, "DYNAMIC_VISCOUS_ML", cs%dynamic_viscous_ML, &
1881  "If true, use a bulk Richardson number criterion to "//&
1882  "determine the mixed layer thickness for viscosity.", &
1883  default=.false.)
1884  if (cs%dynamic_viscous_ML) then
1885  call get_param(param_file, mdl, "BULK_RI_ML", bulk_ri_ml_dflt, default=0.0)
1886  call get_param(param_file, mdl, "BULK_RI_ML_VISC", cs%bulk_Ri_ML, &
1887  "The efficiency with which mean kinetic energy released "//&
1888  "by mechanically forced entrainment of the mixed layer "//&
1889  "is converted to turbulent kinetic energy. By default, "//&
1890  "BULK_RI_ML_VISC = BULK_RI_ML or 0.", units="nondim", &
1891  default=bulk_ri_ml_dflt)
1892  call get_param(param_file, mdl, "TKE_DECAY", tke_decay_dflt, default=0.0)
1893  call get_param(param_file, mdl, "TKE_DECAY_VISC", cs%TKE_decay, &
1894  "TKE_DECAY_VISC relates the vertical rate of decay of "//&
1895  "the TKE available for mechanical entrainment to the "//&
1896  "natural Ekman depth for use in calculating the dynamic "//&
1897  "mixed layer viscosity. By default, "//&
1898  "TKE_DECAY_VISC = TKE_DECAY or 0.", units="nondim", &
1899  default=tke_decay_dflt)
1900  call get_param(param_file, mdl, "ML_USE_OMEGA", use_omega, &
1901  "If true, use the absolute rotation rate instead of the "//&
1902  "vertical component of rotation when setting the decay "//&
1903  "scale for turbulence.", default=.false., do_not_log=.true.)
1904  omega_frac_dflt = 0.0
1905  if (use_omega) then
1906  call mom_error(warning, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
1907  omega_frac_dflt = 1.0
1908  endif
1909  call get_param(param_file, mdl, "ML_OMEGA_FRAC", cs%omega_frac, &
1910  "When setting the decay scale for turbulence, use this "//&
1911  "fraction of the absolute rotation rate blended with the "//&
1912  "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
1913  units="nondim", default=omega_frac_dflt)
1914  call get_param(param_file, mdl, "OMEGA", cs%omega, &
1915  "The rotation rate of the earth.", units="s-1", &
1916  default=7.2921e-5, scale=us%T_to_s)
1917  ! This give a minimum decay scale that is typically much less than Angstrom.
1918  cs%ustar_min = 2e-4*cs%omega*(gv%Angstrom_Z + gv%H_to_Z*gv%H_subroundoff)
1919  else
1920  call get_param(param_file, mdl, "OMEGA", cs%omega, &
1921  "The rotation rate of the earth.", units="s-1", &
1922  default=7.2921e-5, scale=us%T_to_s)
1923  endif
1924 
1925  call get_param(param_file, mdl, "HBBL", cs%Hbbl, &
1926  "The thickness of a bottom boundary layer with a "//&
1927  "viscosity of KVBBL if BOTTOMDRAGLAW is not defined, or "//&
1928  "the thickness over which near-bottom velocities are "//&
1929  "averaged for the drag law if BOTTOMDRAGLAW is defined "//&
1930  "but LINEAR_DRAG is not.", units="m", fail_if_missing=.true.) ! Rescaled later
1931  if (cs%bottomdraglaw) then
1932  call get_param(param_file, mdl, "CDRAG", cs%cdrag, &
1933  "CDRAG is the drag coefficient relating the magnitude of "//&
1934  "the velocity field to the bottom stress. CDRAG is only "//&
1935  "used if BOTTOMDRAGLAW is defined.", units="nondim", &
1936  default=0.003)
1937  call get_param(param_file, mdl, "DRAG_BG_VEL", cs%drag_bg_vel, &
1938  "DRAG_BG_VEL is either the assumed bottom velocity (with "//&
1939  "LINEAR_DRAG) or an unresolved velocity that is "//&
1940  "combined with the resolved velocity to estimate the "//&
1941  "velocity magnitude. DRAG_BG_VEL is only used when "//&
1942  "BOTTOMDRAGLAW is defined.", units="m s-1", default=0.0, scale=us%m_s_to_L_T)
1943  call get_param(param_file, mdl, "BBL_USE_EOS", cs%BBL_use_EOS, &
1944  "If true, use the equation of state in determining the "//&
1945  "properties of the bottom boundary layer. Otherwise use "//&
1946  "the layer target potential densities.", default=.false.)
1947  endif
1948  call get_param(param_file, mdl, "BBL_THICK_MIN", cs%BBL_thick_min, &
1949  "The minimum bottom boundary layer thickness that can be "//&
1950  "used with BOTTOMDRAGLAW. This might be "//&
1951  "Kv/(cdrag*drag_bg_vel) to give Kv as the minimum "//&
1952  "near-bottom viscosity.", units="m", default=0.0) ! Rescaled later
1953  call get_param(param_file, mdl, "HTBL_SHELF_MIN", cs%Htbl_shelf_min, &
1954  "The minimum top boundary layer thickness that can be "//&
1955  "used with BOTTOMDRAGLAW. This might be "//&
1956  "Kv/(cdrag*drag_bg_vel) to give Kv as the minimum "//&
1957  "near-top viscosity.", units="m", default=cs%BBL_thick_min, scale=gv%m_to_H)
1958  call get_param(param_file, mdl, "HTBL_SHELF", cs%Htbl_shelf, &
1959  "The thickness over which near-surface velocities are "//&
1960  "averaged for the drag law under an ice shelf. By "//&
1961  "default this is the same as HBBL", units="m", default=cs%Hbbl, scale=gv%m_to_H)
1962  ! These unit conversions are out outside the get_param calls because the are also defaults.
1963  cs%Hbbl = cs%Hbbl * gv%m_to_H ! Rescale
1964  cs%BBL_thick_min = cs%BBL_thick_min * gv%m_to_H ! Rescale
1965 
1966  call get_param(param_file, mdl, "KV", kv_background, &
1967  "The background kinematic viscosity in the interior. "//&
1968  "The molecular value, ~1e-6 m2 s-1, may be used.", &
1969  units="m2 s-1", fail_if_missing=.true.)
1970 
1971  call get_param(param_file, mdl, "USE_KPP", use_kpp, &
1972  "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
1973  "to calculate diffusivities and non-local transport in the OBL.", &
1974  do_not_log=.true., default=.false.)
1975 
1976  call get_param(param_file, mdl, "KV_BBL_MIN", cs%KV_BBL_min, &
1977  "The minimum viscosities in the bottom boundary layer.", &
1978  units="m2 s-1", default=kv_background, scale=us%m2_s_to_Z2_T)
1979  call get_param(param_file, mdl, "KV_TBL_MIN", cs%KV_TBL_min, &
1980  "The minimum viscosities in the top boundary layer.", &
1981  units="m2 s-1", default=kv_background, scale=us%m2_s_to_Z2_T)
1982 
1983  if (cs%Channel_drag) then
1984  call get_param(param_file, mdl, "SMAG_LAP_CONST", smag_const1, default=-1.0)
1985 
1986  csmag_chan_dflt = 0.15
1987  if (smag_const1 >= 0.0) csmag_chan_dflt = smag_const1
1988 
1989  call get_param(param_file, mdl, "SMAG_CONST_CHANNEL", cs%c_Smag, &
1990  "The nondimensional Laplacian Smagorinsky constant used "//&
1991  "in calculating the channel drag if it is enabled. The "//&
1992  "default is to use the same value as SMAG_LAP_CONST if "//&
1993  "it is defined, or 0.15 if it is not. The value used is "//&
1994  "also 0.15 if the specified value is negative.", &
1995  units="nondim", default=csmag_chan_dflt)
1996  if (cs%c_Smag < 0.0) cs%c_Smag = 0.15
1997  endif
1998 
1999  if (cs%RiNo_mix .and. kappa_shear_at_vertex(param_file)) then
2000  ! This is necessary for reproduciblity across restarts in non-symmetric mode.
2001  call pass_var(visc%Kv_shear_Bu, g%Domain, position=corner, complete=.true.)
2002  endif
2003 
2004  if (cs%bottomdraglaw) then
2005  allocate(visc%bbl_thick_u(isdb:iedb,jsd:jed)) ; visc%bbl_thick_u = 0.0
2006  allocate(visc%kv_bbl_u(isdb:iedb,jsd:jed)) ; visc%kv_bbl_u = 0.0
2007  allocate(visc%bbl_thick_v(isd:ied,jsdb:jedb)) ; visc%bbl_thick_v = 0.0
2008  allocate(visc%kv_bbl_v(isd:ied,jsdb:jedb)) ; visc%kv_bbl_v = 0.0
2009  allocate(visc%ustar_bbl(isd:ied,jsd:jed)) ; visc%ustar_bbl = 0.0
2010  allocate(visc%TKE_bbl(isd:ied,jsd:jed)) ; visc%TKE_bbl = 0.0
2011 
2012  cs%id_bbl_thick_u = register_diag_field('ocean_model', 'bbl_thick_u', &
2013  diag%axesCu1, time, 'BBL thickness at u points', 'm', conversion=us%Z_to_m)
2014  cs%id_kv_bbl_u = register_diag_field('ocean_model', 'kv_bbl_u', diag%axesCu1, &
2015  time, 'BBL viscosity at u points', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
2016  cs%id_bbl_thick_v = register_diag_field('ocean_model', 'bbl_thick_v', &
2017  diag%axesCv1, time, 'BBL thickness at v points', 'm', conversion=us%Z_to_m)
2018  cs%id_kv_bbl_v = register_diag_field('ocean_model', 'kv_bbl_v', diag%axesCv1, &
2019  time, 'BBL viscosity at v points', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
2020  endif
2021  if (cs%Channel_drag) then
2022  allocate(visc%Ray_u(isdb:iedb,jsd:jed,nz)) ; visc%Ray_u = 0.0
2023  allocate(visc%Ray_v(isd:ied,jsdb:jedb,nz)) ; visc%Ray_v = 0.0
2024  cs%id_Ray_u = register_diag_field('ocean_model', 'Rayleigh_u', diag%axesCuL, &
2025  time, 'Rayleigh drag velocity at u points', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
2026  cs%id_Ray_v = register_diag_field('ocean_model', 'Rayleigh_v', diag%axesCvL, &
2027  time, 'Rayleigh drag velocity at v points', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
2028  endif
2029 
2030  if (use_cvmix_ddiff .or. differential_diffusion) then
2031  allocate(visc%Kd_extra_T(isd:ied,jsd:jed,nz+1)) ; visc%Kd_extra_T = 0.0
2032  allocate(visc%Kd_extra_S(isd:ied,jsd:jed,nz+1)) ; visc%Kd_extra_S = 0.0
2033  endif
2034 
2035  if (cs%dynamic_viscous_ML) then
2036  allocate(visc%nkml_visc_u(isdb:iedb,jsd:jed)) ; visc%nkml_visc_u = 0.0
2037  allocate(visc%nkml_visc_v(isd:ied,jsdb:jedb)) ; visc%nkml_visc_v = 0.0
2038  cs%id_nkml_visc_u = register_diag_field('ocean_model', 'nkml_visc_u', &
2039  diag%axesCu1, time, 'Number of layers in viscous mixed layer at u points', 'm')
2040  cs%id_nkml_visc_v = register_diag_field('ocean_model', 'nkml_visc_v', &
2041  diag%axesCv1, time, 'Number of layers in viscous mixed layer at v points', 'm')
2042  endif
2043 
2044  call register_restart_field_as_obsolete('Kd_turb','Kd_shear', restart_cs)
2045  call register_restart_field_as_obsolete('Kv_turb','Kv_shear', restart_cs)
2046 
2047  ! Account for possible changes in dimensional scaling for variables that have been
2048  ! read from a restart file.
2049  z_rescale = 1.0
2050  if ((us%m_to_Z_restart /= 0.0) .and. (us%m_to_Z_restart /= us%m_to_Z)) &
2051  z_rescale = us%m_to_Z / us%m_to_Z_restart
2052  i_t_rescale = 1.0
2053  if ((us%s_to_T_restart /= 0.0) .and. (us%s_to_T_restart /= us%s_to_T)) &
2054  i_t_rescale = us%s_to_T_restart / us%s_to_T
2055  z2_t_rescale = z_rescale**2*i_t_rescale
2056 
2057  if (z2_t_rescale /= 1.0) then
2058  if (associated(visc%Kd_shear)) then ; if (query_initialized(visc%Kd_shear, "Kd_shear", restart_cs)) then
2059  do k=1,nz+1 ; do j=js,je ; do i=is,ie
2060  visc%Kd_shear(i,j,k) = z2_t_rescale * visc%Kd_shear(i,j,k)
2061  enddo ; enddo ; enddo
2062  endif ; endif
2063 
2064  if (associated(visc%Kv_shear)) then ; if (query_initialized(visc%Kv_shear, "Kv_shear", restart_cs)) then
2065  do k=1,nz+1 ; do j=js,je ; do i=is,ie
2066  visc%Kv_shear(i,j,k) = z2_t_rescale * visc%Kv_shear(i,j,k)
2067  enddo ; enddo ; enddo
2068  endif ; endif
2069 
2070  if (associated(visc%Kv_shear_Bu)) then ; if (query_initialized(visc%Kv_shear_Bu, "Kv_shear_Bu", restart_cs)) then
2071  do k=1,nz+1 ; do j=js-1,je ; do i=is-1,ie
2072  visc%Kv_shear_Bu(i,j,k) = z2_t_rescale * visc%Kv_shear_Bu(i,j,k)
2073  enddo ; enddo ; enddo
2074  endif ; endif
2075 
2076  endif
2077 
2078 end subroutine set_visc_init
2079 
2080 !> This subroutine dellocates any memory in the set_visc control structure.
2081 subroutine set_visc_end(visc, CS)
2082  type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
2083  !! related fields. Elements are deallocated here.
2084  type(set_visc_cs), pointer :: cs !< The control structure returned by a previous
2085  !! call to vertvisc_init.
2086  if (cs%bottomdraglaw) then
2087  deallocate(visc%bbl_thick_u) ; deallocate(visc%bbl_thick_v)
2088  deallocate(visc%kv_bbl_u) ; deallocate(visc%kv_bbl_v)
2089  endif
2090  if (cs%Channel_drag) then
2091  deallocate(visc%Ray_u) ; deallocate(visc%Ray_v)
2092  endif
2093  if (cs%dynamic_viscous_ML) then
2094  deallocate(visc%nkml_visc_u) ; deallocate(visc%nkml_visc_v)
2095  endif
2096  if (associated(visc%Kd_shear)) deallocate(visc%Kd_shear)
2097  if (associated(visc%Kv_slow)) deallocate(visc%Kv_slow)
2098  if (associated(visc%TKE_turb)) deallocate(visc%TKE_turb)
2099  if (associated(visc%Kv_shear)) deallocate(visc%Kv_shear)
2100  if (associated(visc%Kv_shear_Bu)) deallocate(visc%Kv_shear_Bu)
2101  if (associated(visc%ustar_bbl)) deallocate(visc%ustar_bbl)
2102  if (associated(visc%TKE_bbl)) deallocate(visc%TKE_bbl)
2103  if (associated(visc%taux_shelf)) deallocate(visc%taux_shelf)
2104  if (associated(visc%tauy_shelf)) deallocate(visc%tauy_shelf)
2105  if (associated(visc%tbl_thick_shelf_u)) deallocate(visc%tbl_thick_shelf_u)
2106  if (associated(visc%tbl_thick_shelf_v)) deallocate(visc%tbl_thick_shelf_v)
2107  if (associated(visc%kv_tbl_shelf_u)) deallocate(visc%kv_tbl_shelf_u)
2108  if (associated(visc%kv_tbl_shelf_v)) deallocate(visc%kv_tbl_shelf_v)
2109 
2110  deallocate(cs)
2111 end subroutine set_visc_end
2112 
2113 !> \namespace mom_set_visc
2114 !!
2115 !! This would also be the module in which other viscous quantities that are flow-independent might be set.
2116 !! This information is transmitted to other modules via a vertvisc type structure.
2117 !!
2118 !! The same code is used for the two velocity components, by indirectly referencing the velocities and
2119 !! defining a handful of direction-specific defined variables.
2120 
2121 end module mom_set_visc
mom_cvmix_ddiff::cvmix_ddiff_is_used
logical function, public cvmix_ddiff_is_used(param_file)
Reads the parameter "USE_CVMIX_DDIFF" and returns state. This function allows other modules to know w...
Definition: MOM_CVMix_ddiff.F90:299
mom_forcing_type::mech_forcing
Structure that contains pointers to the mechanical forcing at the surface used to drive the liquid oc...
Definition: MOM_forcing_type.F90:187
mom_open_boundary::obc_direction_n
integer, parameter, public obc_direction_n
Indicates the boundary is an effective northern boundary.
Definition: MOM_open_boundary.F90:63
mom_open_boundary::obc_direction_s
integer, parameter, public obc_direction_s
Indicates the boundary is an effective southern boundary.
Definition: MOM_open_boundary.F90:64
mom_safe_alloc
Convenience functions for safely allocating memory without accidentally reallocating pointer and caus...
Definition: MOM_safe_alloc.F90:3
mom_set_visc::set_viscous_ml
subroutine, public set_viscous_ml(u, v, h, tv, forces, visc, dt, G, GV, US, CS, symmetrize)
Calculates the thickness of the surface boundary layer for applying an elevated viscosity.
Definition: MOM_set_viscosity.F90:1019
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_restart::register_restart_field_as_obsolete
subroutine, public register_restart_field_as_obsolete(field_name, replacement_name, CS)
Definition: MOM_restart.F90:128
mom_cvmix_conv::cvmix_conv_is_used
logical function, public cvmix_conv_is_used(param_file)
Reads the parameter "USE_CVMix_CONVECTION" and returns state. This function allows other modules to k...
Definition: MOM_CVMix_conv.F90:264
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_set_visc::set_v_at_u
real function set_v_at_u(v, h, G, i, j, k, mask2dCv, OBC)
This subroutine finds a thickness-weighted value of v at the u-points.
Definition: MOM_set_viscosity.F90:927
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_variables::thermo_var_ptrs
Pointers to an assortment of thermodynamic fields that may be available, including potential temperat...
Definition: MOM_variables.F90:78
mom_open_boundary::obc_direction_e
integer, parameter, public obc_direction_e
Indicates the boundary is an effective eastern boundary.
Definition: MOM_open_boundary.F90:65
mom_open_boundary::obc_direction_w
integer, parameter, public obc_direction_w
Indicates the boundary is an effective western boundary.
Definition: MOM_open_boundary.F90:66
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes_in, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1878
mom_set_visc::set_viscous_bbl
subroutine, public set_viscous_bbl(u, v, h, tv, visc, G, GV, US, CS, symmetrize)
Calculates the thickness of the bottom boundary layer and the viscosity within that layer....
Definition: MOM_set_viscosity.F90:110
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_domains::pass_var
Do a halo update on an array.
Definition: MOM_domains.F90:49
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_hor_index
Defines the horizontal index type (hor_index_type) used for providing index ranges.
Definition: MOM_hor_index.F90:2
mom_restart::mom_restart_cs
A restart registry and the control structure for restarts.
Definition: MOM_restart.F90:72
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_safe_alloc::safe_alloc_alloc
Allocate a 2-d or 3-d allocatable array.
Definition: MOM_safe_alloc.F90:18
mom_forcing_type
This module implements boundary forcing for MOM6.
Definition: MOM_forcing_type.F90:2
mom_set_visc::set_visc_register_restarts
subroutine, public set_visc_register_restarts(HI, GV, param_file, visc, restart_CS)
Register any fields associated with the vertvisc_type.
Definition: MOM_set_viscosity.F90:1696
mom_variables::vertvisc_type
Vertical viscosities, drag coefficients, and related fields.
Definition: MOM_variables.F90:196
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_restart
The MOM6 facility for reading and writing restart files, and querying what has been read.
Definition: MOM_restart.F90:2
mom_cvmix_shear::cvmix_shear_is_used
logical function, public cvmix_shear_is_used(param_file)
Reads the parameters "LMD94" and "PP81" and returns state. This function allows other modules to know...
Definition: MOM_CVMix_shear.F90:303
mom_kappa_shear::kappa_shear_at_vertex
logical function, public kappa_shear_at_vertex(param_file)
This function indicates to other modules whether the Jackson et al shear mixing parameterization will...
Definition: MOM_kappa_shear.F90:2120
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_variables
Provides transparent structures with groups of MOM6 variables and supporting routines.
Definition: MOM_variables.F90:2
mom_kappa_shear::kappa_shear_is_used
logical function, public kappa_shear_is_used(param_file)
This function indicates to other modules whether the Jackson et al shear mixing parameterization will...
Definition: MOM_kappa_shear.F90:2109
mom_open_boundary
Controls where open boundary conditions are applied.
Definition: MOM_open_boundary.F90:2
mom_set_visc::set_visc_init
subroutine, public set_visc_init(Time, G, GV, US, param_file, diag, visc, CS, restart_CS, OBC)
Initializes the MOM_set_visc control structure.
Definition: MOM_set_viscosity.F90:1782
mom_cpu_clock
Wraps the MPP cpu clock functions.
Definition: MOM_cpu_clock.F90:2
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_eos::calculate_density_derivs
Calculate the derivatives of density with temperature and salinity from T, S, and P.
Definition: MOM_EOS.F90:70
mom_cvmix_ddiff
Interface to CVMix double diffusion scheme.
Definition: MOM_CVMix_ddiff.F90:2
mom_hor_index::hor_index_type
Container for horizontal index ranges for data, computational and global domains.
Definition: MOM_hor_index.F90:15
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_cvmix_shear
Interface to CVMix interior shear schemes.
Definition: MOM_CVMix_shear.F90:2
mom_open_boundary::ocean_obc_type
Open-boundary data.
Definition: MOM_open_boundary.F90:195
mom_restart::register_restart_field
Register fields for restarts.
Definition: MOM_restart.F90:107
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_forcing_type::forcing
Structure that contains pointers to the boundary forcing used to drive the liquid ocean simulated by ...
Definition: MOM_forcing_type.F90:50
mom_set_visc
Calculates various values related to the bottom boundary layer, such as the viscosity and thickness o...
Definition: MOM_set_viscosity.F90:3
mom_debugging
Provides checksumming functions for debugging.
Definition: MOM_debugging.F90:7
mom_set_visc::set_visc_end
subroutine, public set_visc_end(visc, CS)
This subroutine dellocates any memory in the set_visc control structure.
Definition: MOM_set_viscosity.F90:2082
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_open_boundary::obc_none
integer, parameter, public obc_none
Indicates the use of no open boundary.
Definition: MOM_open_boundary.F90:58
mom_file_parser::log_param
An overloaded interface to log the values of various types of parameters.
Definition: MOM_file_parser.F90:96
mom_restart::query_initialized
Indicate whether a field has been read from a restart file.
Definition: MOM_restart.F90:116
mom_safe_alloc::safe_alloc_ptr
Allocate a pointer to a 1-d, 2-d or 3-d array.
Definition: MOM_safe_alloc.F90:12
mom_open_boundary::obc_segment_type
Open boundary segment data structure.
Definition: MOM_open_boundary.F90:103
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_set_visc::set_u_at_v
real function set_u_at_v(u, h, G, i, j, k, mask2dCu, OBC)
This subroutine finds a thickness-weighted value of u at the v-points.
Definition: MOM_set_viscosity.F90:971
mom_kappa_shear
Shear-dependent mixing following Jackson et al. 2008.
Definition: MOM_kappa_shear.F90:2
mom_set_visc::set_visc_cs
Control structure for MOM_set_visc.
Definition: MOM_set_viscosity.F90:44
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239
mom_cvmix_conv
Interface to CVMix convection scheme.
Definition: MOM_CVMix_conv.F90:2
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60