MOM6
Data Types | Functions/Subroutines
mom_set_visc Module Reference

Detailed Description

Calculates various values related to the bottom boundary layer, such as the viscosity and thickness of the BBL (set_viscous_BBL).

This would also be the module in which other viscous quantities that are flow-independent might be set. This information is transmitted to other modules via a vertvisc type structure.

The same code is used for the two velocity components, by indirectly referencing the velocities and defining a handful of direction-specific defined variables.

Data Types

type  set_visc_cs
 Control structure for MOM_set_visc. More...
 

Functions/Subroutines

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. More...
 
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. More...
 
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. More...
 
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. More...
 
subroutine, public set_visc_register_restarts (HI, GV, param_file, visc, restart_CS)
 Register any fields associated with the vertvisc_type. More...
 
subroutine, public set_visc_init (Time, G, GV, US, param_file, diag, visc, CS, restart_CS, OBC)
 Initializes the MOM_set_visc control structure. More...
 
subroutine, public set_visc_end (visc, CS)
 This subroutine dellocates any memory in the set_visc control structure. More...
 

Function/Subroutine Documentation

◆ set_u_at_v()

real function mom_set_visc::set_u_at_v ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
type(ocean_grid_type), intent(in)  G,
integer, intent(in)  i,
integer, intent(in)  j,
integer, intent(in)  k,
real, dimension( g %isdb: g %iedb, g %jsd: g %jed), intent(in)  mask2dCu,
type(ocean_obc_type), pointer  OBC 
)
private

This subroutine finds a thickness-weighted value of u at the v-points.

Parameters
[in]gThe ocean's grid structure
[in]uThe zonal velocity [L T-1 ~> m s-1] or other units.
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]iThe i-index of the u-location to work on.
[in]jThe j-index of the u-location to work on.
[in]kThe k-index of the u-location to work on.
[in]mask2dcuA multiplicative mask of the u-points
obcA pointer to an open boundary condition structure
Returns
The return value of u at v points in the same units as u, i.e. [L T-1 ~> m s-1] or other units.

Definition at line 1160 of file MOM_set_viscosity.F90.

1160  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1161  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1162  intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1] or other units.
1163  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1164  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1165  integer, intent(in) :: i !< The i-index of the u-location to work on.
1166  integer, intent(in) :: j !< The j-index of the u-location to work on.
1167  integer, intent(in) :: k !< The k-index of the u-location to work on.
1168  real, dimension(SZIB_(G),SZJ_(G)), &
1169  intent(in) :: mask2dCu !< A multiplicative mask of the u-points
1170  type(ocean_OBC_type), pointer :: OBC !< A pointer to an open boundary condition structure
1171  real :: set_u_at_v !< The return value of u at v points in the
1172  !! same units as u, i.e. [L T-1 ~> m s-1] or other units.
1173 
1174  ! This subroutine finds a thickness-weighted value of u at the v-points.
1175  real :: hwt(-1:0,0:1) ! Masked weights used to average u onto v [H ~> m or kg m-2].
1176  real :: hwt_tot ! The sum of the masked thicknesses [H ~> m or kg m-2].
1177  integer :: i0, j0, i1, j1
1178 
1179  do j0 = 0,1 ; do i0 = -1,0 ; i1 = i+i0 ; j1 = j+j0
1180  hwt(i0,j0) = (h(i1,j1,k) + h(i1+1,j1,k)) * mask2dcu(i1,j1)
1181  enddo ; enddo
1182 
1183  if (associated(obc)) then ; if (obc%number_of_segments > 0) then
1184  do j0 = 0,1 ; do i0 = -1,0 ; if ((obc%segnum_u(i+i0,j+j0) /= obc_none)) then
1185  i1 = i+i0 ; j1 = j+j0
1186  if (obc%segment(obc%segnum_u(i1,j1))%direction == obc_direction_e) then
1187  hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dcu(i1,j1)
1188  elseif (obc%segment(obc%segnum_u(i1,j1))%direction == obc_direction_w) then
1189  hwt(i0,j0) = 2.0 * h(i1+1,j1,k) * mask2dcu(i1,j1)
1190  endif
1191  endif ; enddo ; enddo
1192  endif ; endif
1193 
1194  hwt_tot = (hwt(-1,0) + hwt(0,1)) + (hwt(0,0) + hwt(-1,1))
1195  set_u_at_v = 0.0
1196  if (hwt_tot > 0.0) set_u_at_v = &
1197  ((hwt(0,0) * u(i,j,k) + hwt(-1,1) * u(i-1,j+1,k)) + &
1198  (hwt(-1,0) * u(i-1,j,k) + hwt(0,1) * u(i,j+1,k))) / hwt_tot
1199 

◆ set_v_at_u()

real function mom_set_visc::set_v_at_u ( real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
type(ocean_grid_type), intent(in)  G,
integer, intent(in)  i,
integer, intent(in)  j,
integer, intent(in)  k,
real, dimension(szi_(g),szjb_(g)), intent(in)  mask2dCv,
type(ocean_obc_type), pointer  OBC 
)
private

This subroutine finds a thickness-weighted value of v at the u-points.

Parameters
[in]gThe ocean's grid structure
[in]vThe meridional velocity [L T-1 ~> m s-1]
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]iThe i-index of the u-location to work on.
[in]jThe j-index of the u-location to work on.
[in]kThe k-index of the u-location to work on.
[in]mask2dcvA multiplicative mask of the v-points
obcA pointer to an open boundary condition structure
Returns
The return value of v at u points points in the same units as u, i.e. [L T-1 ~> m s-1] or other units.

Definition at line 1116 of file MOM_set_viscosity.F90.

1116  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1117  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1118  intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
1119  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1120  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1121  integer, intent(in) :: i !< The i-index of the u-location to work on.
1122  integer, intent(in) :: j !< The j-index of the u-location to work on.
1123  integer, intent(in) :: k !< The k-index of the u-location to work on.
1124  real, dimension(SZI_(G),SZJB_(G)),&
1125  intent(in) :: mask2dCv !< A multiplicative mask of the v-points
1126  type(ocean_OBC_type), pointer :: OBC !< A pointer to an open boundary condition structure
1127  real :: set_v_at_u !< The return value of v at u points points in the
1128  !! same units as u, i.e. [L T-1 ~> m s-1] or other units.
1129 
1130  ! This subroutine finds a thickness-weighted value of v at the u-points.
1131  real :: hwt(0:1,-1:0) ! Masked weights used to average u onto v [H ~> m or kg m-2].
1132  real :: hwt_tot ! The sum of the masked thicknesses [H ~> m or kg m-2].
1133  integer :: i0, j0, i1, j1
1134 
1135  do j0 = -1,0 ; do i0 = 0,1 ; i1 = i+i0 ; j1 = j+j0
1136  hwt(i0,j0) = (h(i1,j1,k) + h(i1,j1+1,k)) * mask2dcv(i1,j1)
1137  enddo ; enddo
1138 
1139  if (associated(obc)) then ; if (obc%number_of_segments > 0) then
1140  do j0 = -1,0 ; do i0 = 0,1 ; if ((obc%segnum_v(i+i0,j+j0) /= obc_none)) then
1141  i1 = i+i0 ; j1 = j+j0
1142  if (obc%segment(obc%segnum_v(i1,j1))%direction == obc_direction_n) then
1143  hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dcv(i1,j1)
1144  elseif (obc%segment(obc%segnum_v(i1,j1))%direction == obc_direction_s) then
1145  hwt(i0,j0) = 2.0 * h(i1,j1+1,k) * mask2dcv(i1,j1)
1146  endif
1147  endif ; enddo ; enddo
1148  endif ; endif
1149 
1150  hwt_tot = (hwt(0,-1) + hwt(1,0)) + (hwt(1,-1) + hwt(0,0))
1151  set_v_at_u = 0.0
1152  if (hwt_tot > 0.0) set_v_at_u = &
1153  ((hwt(0,0) * v(i,j,k) + hwt(1,-1) * v(i+1,j-1,k)) + &
1154  (hwt(1,0) * v(i+1,j,k) + hwt(0,-1) * v(i,j-1,k))) / hwt_tot
1155 

◆ set_visc_end()

subroutine, public mom_set_visc::set_visc_end ( type(vertvisc_type), intent(inout)  visc,
type(set_visc_cs), pointer  CS 
)

This subroutine dellocates any memory in the set_visc control structure.

Parameters
[in,out]viscA structure containing vertical viscosities and related fields. Elements are deallocated here.
csThe control structure returned by a previous call to set_visc_init.

Definition at line 2315 of file MOM_set_viscosity.F90.

2315  type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
2316  !! related fields. Elements are deallocated here.
2317  type(set_visc_CS), pointer :: CS !< The control structure returned by a previous
2318  !! call to set_visc_init.
2319  if (cs%bottomdraglaw) then
2320  deallocate(visc%bbl_thick_u) ; deallocate(visc%bbl_thick_v)
2321  deallocate(visc%kv_bbl_u) ; deallocate(visc%kv_bbl_v)
2322  if (allocated(cs%bbl_u)) deallocate(cs%bbl_u)
2323  if (allocated(cs%bbl_v)) deallocate(cs%bbl_v)
2324  endif
2325  if (cs%Channel_drag) then
2326  deallocate(visc%Ray_u) ; deallocate(visc%Ray_v)
2327  endif
2328  if (cs%dynamic_viscous_ML) then
2329  deallocate(visc%nkml_visc_u) ; deallocate(visc%nkml_visc_v)
2330  endif
2331  if (associated(visc%Kd_shear)) deallocate(visc%Kd_shear)
2332  if (associated(visc%Kv_slow)) deallocate(visc%Kv_slow)
2333  if (associated(visc%TKE_turb)) deallocate(visc%TKE_turb)
2334  if (associated(visc%Kv_shear)) deallocate(visc%Kv_shear)
2335  if (associated(visc%Kv_shear_Bu)) deallocate(visc%Kv_shear_Bu)
2336  if (associated(visc%ustar_bbl)) deallocate(visc%ustar_bbl)
2337  if (associated(visc%TKE_bbl)) deallocate(visc%TKE_bbl)
2338  if (associated(visc%taux_shelf)) deallocate(visc%taux_shelf)
2339  if (associated(visc%tauy_shelf)) deallocate(visc%tauy_shelf)
2340  if (associated(visc%tbl_thick_shelf_u)) deallocate(visc%tbl_thick_shelf_u)
2341  if (associated(visc%tbl_thick_shelf_v)) deallocate(visc%tbl_thick_shelf_v)
2342  if (associated(visc%kv_tbl_shelf_u)) deallocate(visc%kv_tbl_shelf_u)
2343  if (associated(visc%kv_tbl_shelf_v)) deallocate(visc%kv_tbl_shelf_v)
2344 
2345  deallocate(cs)

◆ set_visc_init()

subroutine, public mom_set_visc::set_visc_init ( type(time_type), intent(in), target  Time,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(diag_ctrl), intent(inout), target  diag,
type(vertvisc_type), intent(inout)  visc,
type(set_visc_cs), pointer  CS,
type(mom_restart_cs), pointer  restart_CS,
type(ocean_obc_type), pointer  OBC 
)

Initializes the MOM_set_visc control structure.

Parameters
[in]timeThe current model time.
[in,out]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]param_fileA structure to parse for run-time parameters.
[in,out]diagA structure that is used to regulate diagnostic output.
[in,out]viscA structure containing vertical viscosities and related fields. Allocated here.
csA pointer that is set to point to the control structure for this module
restart_csA pointer to the restart control structure.
obcA pointer to an open boundary condition structure

Definition at line 1973 of file MOM_set_viscosity.F90.

1973  type(time_type), target, intent(in) :: Time !< The current model time.
1974  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
1975  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1976  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1977  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1978  !! parameters.
1979  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
1980  !! output.
1981  type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
1982  !! related fields. Allocated here.
1983  type(set_visc_CS), pointer :: CS !< A pointer that is set to point to the control
1984  !! structure for this module
1985  type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure.
1986  type(ocean_OBC_type), pointer :: OBC !< A pointer to an open boundary condition structure
1987 
1988  ! Local variables
1989  real :: Csmag_chan_dflt, smag_const1, TKE_decay_dflt, bulk_Ri_ML_dflt
1990  real :: Kv_background
1991  real :: omega_frac_dflt
1992  real :: Z_rescale ! A rescaling factor for heights from the representation in
1993  ! a restart file to the internal representation in this run.
1994  real :: I_T_rescale ! A rescaling factor for time from the internal representation in this run
1995  ! to the representation in a restart file.
1996  real :: Z2_T_rescale ! A rescaling factor for vertical diffusivities and viscosities from the
1997  ! representation in a restart file to the internal representation in this run.
1998  integer :: i, j, k, is, ie, js, je, n
1999  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
2000  logical :: default_2018_answers
2001  logical :: use_kappa_shear, adiabatic, use_omega, MLE_use_PBL_MLD
2002  logical :: use_KPP
2003  logical :: use_regridding
2004  character(len=200) :: filename, tideamp_file
2005  type(OBC_segment_type), pointer :: segment => null() ! pointer to OBC segment type
2006  ! This include declares and sets the variable "version".
2007 # include "version_variable.h"
2008  character(len=40) :: mdl = "MOM_set_visc" ! This module's name.
2009 
2010  if (associated(cs)) then
2011  call mom_error(warning, "set_visc_init called with an associated "// &
2012  "control structure.")
2013  return
2014  endif
2015  allocate(cs)
2016 
2017  cs%OBC => obc
2018 
2019  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
2020  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
2021  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
2022 
2023  cs%diag => diag
2024 
2025  ! Set default, read and log parameters
2026  call log_version(param_file, mdl, version, "")
2027  cs%RiNo_mix = .false.
2028  call get_param(param_file, mdl, "INPUTDIR", cs%inputdir, default=".")
2029  cs%inputdir = slasher(cs%inputdir)
2030  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
2031  "This sets the default value for the various _2018_ANSWERS parameters.", &
2032  default=.false.)
2033  call get_param(param_file, mdl, "SET_VISC_2018_ANSWERS", cs%answers_2018, &
2034  "If true, use the order of arithmetic and expressions that recover the "//&
2035  "answers from the end of 2018. Otherwise, use updated and more robust "//&
2036  "forms of the same expressions.", default=default_2018_answers)
2037  call get_param(param_file, mdl, "BOTTOMDRAGLAW", cs%bottomdraglaw, &
2038  "If true, the bottom stress is calculated with a drag "//&
2039  "law of the form c_drag*|u|*u. The velocity magnitude "//&
2040  "may be an assumed value or it may be based on the "//&
2041  "actual velocity in the bottommost HBBL, depending on "//&
2042  "LINEAR_DRAG.", default=.true.)
2043  call get_param(param_file, mdl, "CHANNEL_DRAG", cs%Channel_drag, &
2044  "If true, the bottom drag is exerted directly on each "//&
2045  "layer proportional to the fraction of the bottom it "//&
2046  "overlies.", default=.false.)
2047  call get_param(param_file, mdl, "LINEAR_DRAG", cs%linear_drag, &
2048  "If LINEAR_DRAG and BOTTOMDRAGLAW are defined the drag "//&
2049  "law is cdrag*DRAG_BG_VEL*u.", default=.false.)
2050  call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
2051  do_not_log=.true.)
2052  if (adiabatic) then
2053  call log_param(param_file, mdl, "ADIABATIC",adiabatic, &
2054  "There are no diapycnal mass fluxes if ADIABATIC is "//&
2055  "true. This assumes that KD = KDML = 0.0 and that "//&
2056  "there is no buoyancy forcing, but makes the model "//&
2057  "faster by eliminating subroutine calls.", default=.false.)
2058  endif
2059 
2060  if (.not.adiabatic) then
2061  cs%RiNo_mix = kappa_shear_is_used(param_file)
2062  endif
2063 
2064  call get_param(param_file, mdl, "PRANDTL_TURB", visc%Prandtl_turb, &
2065  "The turbulent Prandtl number applied to shear "//&
2066  "instability.", units="nondim", default=1.0)
2067  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
2068 
2069  call get_param(param_file, mdl, "DYNAMIC_VISCOUS_ML", cs%dynamic_viscous_ML, &
2070  "If true, use a bulk Richardson number criterion to "//&
2071  "determine the mixed layer thickness for viscosity.", &
2072  default=.false.)
2073  if (cs%dynamic_viscous_ML) then
2074  call get_param(param_file, mdl, "BULK_RI_ML", bulk_ri_ml_dflt, default=0.0)
2075  call get_param(param_file, mdl, "BULK_RI_ML_VISC", cs%bulk_Ri_ML, &
2076  "The efficiency with which mean kinetic energy released "//&
2077  "by mechanically forced entrainment of the mixed layer "//&
2078  "is converted to turbulent kinetic energy. By default, "//&
2079  "BULK_RI_ML_VISC = BULK_RI_ML or 0.", units="nondim", &
2080  default=bulk_ri_ml_dflt)
2081  call get_param(param_file, mdl, "TKE_DECAY", tke_decay_dflt, default=0.0)
2082  call get_param(param_file, mdl, "TKE_DECAY_VISC", cs%TKE_decay, &
2083  "TKE_DECAY_VISC relates the vertical rate of decay of "//&
2084  "the TKE available for mechanical entrainment to the "//&
2085  "natural Ekman depth for use in calculating the dynamic "//&
2086  "mixed layer viscosity. By default, "//&
2087  "TKE_DECAY_VISC = TKE_DECAY or 0.", units="nondim", &
2088  default=tke_decay_dflt)
2089  call get_param(param_file, mdl, "ML_USE_OMEGA", use_omega, &
2090  "If true, use the absolute rotation rate instead of the "//&
2091  "vertical component of rotation when setting the decay "//&
2092  "scale for turbulence.", default=.false., do_not_log=.true.)
2093  omega_frac_dflt = 0.0
2094  if (use_omega) then
2095  call mom_error(warning, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
2096  omega_frac_dflt = 1.0
2097  endif
2098  call get_param(param_file, mdl, "ML_OMEGA_FRAC", cs%omega_frac, &
2099  "When setting the decay scale for turbulence, use this "//&
2100  "fraction of the absolute rotation rate blended with the "//&
2101  "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
2102  units="nondim", default=omega_frac_dflt)
2103  call get_param(param_file, mdl, "OMEGA", cs%omega, &
2104  "The rotation rate of the earth.", units="s-1", &
2105  default=7.2921e-5, scale=us%T_to_s)
2106  ! This give a minimum decay scale that is typically much less than Angstrom.
2107  cs%ustar_min = 2e-4*cs%omega*(gv%Angstrom_Z + gv%H_to_Z*gv%H_subroundoff)
2108  else
2109  call get_param(param_file, mdl, "OMEGA", cs%omega, &
2110  "The rotation rate of the earth.", units="s-1", &
2111  default=7.2921e-5, scale=us%T_to_s)
2112  endif
2113 
2114  call get_param(param_file, mdl, "HBBL", cs%Hbbl, &
2115  "The thickness of a bottom boundary layer with a "//&
2116  "viscosity of KVBBL if BOTTOMDRAGLAW is not defined, or "//&
2117  "the thickness over which near-bottom velocities are "//&
2118  "averaged for the drag law if BOTTOMDRAGLAW is defined "//&
2119  "but LINEAR_DRAG is not.", units="m", fail_if_missing=.true.) ! Rescaled later
2120  if (cs%bottomdraglaw) then
2121  call get_param(param_file, mdl, "CDRAG", cs%cdrag, &
2122  "CDRAG is the drag coefficient relating the magnitude of "//&
2123  "the velocity field to the bottom stress. CDRAG is only "//&
2124  "used if BOTTOMDRAGLAW is defined.", units="nondim", &
2125  default=0.003)
2126  call get_param(param_file, mdl, "BBL_USE_TIDAL_BG", cs%BBL_use_tidal_bg, &
2127  "Flag to use the tidal RMS amplitude in place of constant "//&
2128  "background velocity for computing u* in the BBL. "//&
2129  "This flag is only used when BOTTOMDRAGLAW is true and "//&
2130  "LINEAR_DRAG is false.", default=.false.)
2131  if (cs%BBL_use_tidal_bg) then
2132  call get_param(param_file, mdl, "TIDEAMP_FILE", tideamp_file, &
2133  "The path to the file containing the spatially varying "//&
2134  "tidal amplitudes with INT_TIDE_DISSIPATION.", default="tideamp.nc")
2135  else
2136  call get_param(param_file, mdl, "DRAG_BG_VEL", cs%drag_bg_vel, &
2137  "DRAG_BG_VEL is either the assumed bottom velocity (with "//&
2138  "LINEAR_DRAG) or an unresolved velocity that is "//&
2139  "combined with the resolved velocity to estimate the "//&
2140  "velocity magnitude. DRAG_BG_VEL is only used when "//&
2141  "BOTTOMDRAGLAW is defined.", units="m s-1", default=0.0, scale=us%m_s_to_L_T)
2142  endif
2143  call get_param(param_file, mdl, "USE_REGRIDDING", use_regridding, &
2144  do_not_log = .true., default = .false. )
2145  call get_param(param_file, mdl, "BBL_USE_EOS", cs%BBL_use_EOS, &
2146  "If true, use the equation of state in determining the "//&
2147  "properties of the bottom boundary layer. Otherwise use "//&
2148  "the layer target potential densities. The default of "//&
2149  "this is determined by USE_REGRIDDING.", default=use_regridding)
2150  if (use_regridding .and. (.not. cs%BBL_use_EOS)) &
2151  call mom_error(fatal,"When using MOM6 in ALE mode it is required to "//&
2152  "set BBL_USE_EOS to True")
2153  endif
2154  call get_param(param_file, mdl, "BBL_THICK_MIN", cs%BBL_thick_min, &
2155  "The minimum bottom boundary layer thickness that can be "//&
2156  "used with BOTTOMDRAGLAW. This might be "//&
2157  "Kv/(cdrag*drag_bg_vel) to give Kv as the minimum "//&
2158  "near-bottom viscosity.", units="m", default=0.0) ! Rescaled later
2159  call get_param(param_file, mdl, "HTBL_SHELF_MIN", cs%Htbl_shelf_min, &
2160  "The minimum top boundary layer thickness that can be "//&
2161  "used with BOTTOMDRAGLAW. This might be "//&
2162  "Kv/(cdrag*drag_bg_vel) to give Kv as the minimum "//&
2163  "near-top viscosity.", units="m", default=cs%BBL_thick_min, scale=gv%m_to_H)
2164  call get_param(param_file, mdl, "HTBL_SHELF", cs%Htbl_shelf, &
2165  "The thickness over which near-surface velocities are "//&
2166  "averaged for the drag law under an ice shelf. By "//&
2167  "default this is the same as HBBL", units="m", default=cs%Hbbl, scale=gv%m_to_H)
2168  ! These unit conversions are out outside the get_param calls because the are also defaults.
2169  cs%Hbbl = cs%Hbbl * gv%m_to_H ! Rescale
2170  cs%BBL_thick_min = cs%BBL_thick_min * gv%m_to_H ! Rescale
2171 
2172  call get_param(param_file, mdl, "KV", kv_background, &
2173  "The background kinematic viscosity in the interior. "//&
2174  "The molecular value, ~1e-6 m2 s-1, may be used.", &
2175  units="m2 s-1", fail_if_missing=.true.)
2176 
2177  call get_param(param_file, mdl, "USE_KPP", use_kpp, &
2178  "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
2179  "to calculate diffusivities and non-local transport in the OBL.", &
2180  do_not_log=.true., default=.false.)
2181 
2182  call get_param(param_file, mdl, "KV_BBL_MIN", cs%KV_BBL_min, &
2183  "The minimum viscosities in the bottom boundary layer.", &
2184  units="m2 s-1", default=kv_background, scale=us%m2_s_to_Z2_T)
2185  call get_param(param_file, mdl, "KV_TBL_MIN", cs%KV_TBL_min, &
2186  "The minimum viscosities in the top boundary layer.", &
2187  units="m2 s-1", default=kv_background, scale=us%m2_s_to_Z2_T)
2188  call get_param(param_file, mdl, "CORRECT_BBL_BOUNDS", cs%correct_BBL_bounds, &
2189  "If true, uses the correct bounds on the BBL thickness and "//&
2190  "viscosity so that the bottom layer feels the intended drag.", &
2191  default=.false.)
2192 
2193  if (cs%Channel_drag) then
2194  call get_param(param_file, mdl, "SMAG_LAP_CONST", smag_const1, default=-1.0)
2195 
2196  csmag_chan_dflt = 0.15
2197  if (smag_const1 >= 0.0) csmag_chan_dflt = smag_const1
2198 
2199  call get_param(param_file, mdl, "SMAG_CONST_CHANNEL", cs%c_Smag, &
2200  "The nondimensional Laplacian Smagorinsky constant used "//&
2201  "in calculating the channel drag if it is enabled. The "//&
2202  "default is to use the same value as SMAG_LAP_CONST if "//&
2203  "it is defined, or 0.15 if it is not. The value used is "//&
2204  "also 0.15 if the specified value is negative.", &
2205  units="nondim", default=csmag_chan_dflt)
2206  if (cs%c_Smag < 0.0) cs%c_Smag = 0.15
2207  endif
2208 
2209  call get_param(param_file, mdl, "MLE_USE_PBL_MLD", mle_use_pbl_mld, &
2210  default=.false., do_not_log=.true.)
2211 
2212  if (cs%RiNo_mix .and. kappa_shear_at_vertex(param_file)) then
2213  ! This is necessary for reproduciblity across restarts in non-symmetric mode.
2214  call pass_var(visc%Kv_shear_Bu, g%Domain, position=corner, complete=.true.)
2215  endif
2216 
2217  if (cs%bottomdraglaw) then
2218  allocate(visc%bbl_thick_u(isdb:iedb,jsd:jed)) ; visc%bbl_thick_u(:,:) = 0.0
2219  allocate(visc%kv_bbl_u(isdb:iedb,jsd:jed)) ; visc%kv_bbl_u(:,:) = 0.0
2220  allocate(visc%bbl_thick_v(isd:ied,jsdb:jedb)) ; visc%bbl_thick_v(:,:) = 0.0
2221  allocate(visc%kv_bbl_v(isd:ied,jsdb:jedb)) ; visc%kv_bbl_v(:,:) = 0.0
2222  allocate(visc%ustar_bbl(isd:ied,jsd:jed)) ; visc%ustar_bbl(:,:) = 0.0
2223  allocate(visc%TKE_bbl(isd:ied,jsd:jed)) ; visc%TKE_bbl(:,:) = 0.0
2224 
2225  cs%id_bbl_thick_u = register_diag_field('ocean_model', 'bbl_thick_u', &
2226  diag%axesCu1, time, 'BBL thickness at u points', 'm', conversion=us%Z_to_m)
2227  cs%id_kv_bbl_u = register_diag_field('ocean_model', 'kv_bbl_u', diag%axesCu1, &
2228  time, 'BBL viscosity at u points', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
2229  cs%id_bbl_u = register_diag_field('ocean_model', 'bbl_u', diag%axesCu1, &
2230  time, 'BBL mean u current', 'm s-1', conversion=us%L_T_to_m_s)
2231  if (cs%id_bbl_u>0) then
2232  allocate(cs%bbl_u(isdb:iedb,jsd:jed)) ; cs%bbl_u(:,:) = 0.0
2233  endif
2234  cs%id_bbl_thick_v = register_diag_field('ocean_model', 'bbl_thick_v', &
2235  diag%axesCv1, time, 'BBL thickness at v points', 'm', conversion=us%Z_to_m)
2236  cs%id_kv_bbl_v = register_diag_field('ocean_model', 'kv_bbl_v', diag%axesCv1, &
2237  time, 'BBL viscosity at v points', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
2238  cs%id_bbl_v = register_diag_field('ocean_model', 'bbl_v', diag%axesCv1, &
2239  time, 'BBL mean v current', 'm s-1', conversion=us%L_T_to_m_s)
2240  if (cs%id_bbl_v>0) then
2241  allocate(cs%bbl_v(isd:ied,jsdb:jedb)) ; cs%bbl_v(:,:) = 0.0
2242  endif
2243  if (cs%BBL_use_tidal_bg) then
2244  allocate(cs%tideamp(isd:ied,jsd:jed)) ; cs%tideamp(:,:) = 0.0
2245  filename = trim(cs%inputdir) // trim(tideamp_file)
2246  call log_param(param_file, mdl, "INPUTDIR/TIDEAMP_FILE", filename)
2247  call mom_read_data(filename, 'tideamp', cs%tideamp, g%domain, timelevel=1, scale=us%m_to_Z*us%T_to_s)
2248  call pass_var(cs%tideamp,g%domain)
2249  endif
2250  endif
2251  if (cs%Channel_drag) then
2252  allocate(visc%Ray_u(isdb:iedb,jsd:jed,nz)) ; visc%Ray_u(:,:,:) = 0.0
2253  allocate(visc%Ray_v(isd:ied,jsdb:jedb,nz)) ; visc%Ray_v(:,:,:) = 0.0
2254  cs%id_Ray_u = register_diag_field('ocean_model', 'Rayleigh_u', diag%axesCuL, &
2255  time, 'Rayleigh drag velocity at u points', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
2256  cs%id_Ray_v = register_diag_field('ocean_model', 'Rayleigh_v', diag%axesCvL, &
2257  time, 'Rayleigh drag velocity at v points', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
2258  endif
2259 
2260 
2261  if (cs%dynamic_viscous_ML) then
2262  allocate(visc%nkml_visc_u(isdb:iedb,jsd:jed)) ; visc%nkml_visc_u(:,:) = 0.0
2263  allocate(visc%nkml_visc_v(isd:ied,jsdb:jedb)) ; visc%nkml_visc_v(:,:) = 0.0
2264  cs%id_nkml_visc_u = register_diag_field('ocean_model', 'nkml_visc_u', &
2265  diag%axesCu1, time, 'Number of layers in viscous mixed layer at u points', 'm')
2266  cs%id_nkml_visc_v = register_diag_field('ocean_model', 'nkml_visc_v', &
2267  diag%axesCv1, time, 'Number of layers in viscous mixed layer at v points', 'm')
2268  endif
2269 
2270  call register_restart_field_as_obsolete('Kd_turb','Kd_shear', restart_cs)
2271  call register_restart_field_as_obsolete('Kv_turb','Kv_shear', restart_cs)
2272 
2273  ! Account for possible changes in dimensional scaling for variables that have been
2274  ! read from a restart file.
2275  z_rescale = 1.0
2276  if ((us%m_to_Z_restart /= 0.0) .and. (us%m_to_Z_restart /= us%m_to_Z)) &
2277  z_rescale = us%m_to_Z / us%m_to_Z_restart
2278  i_t_rescale = 1.0
2279  if ((us%s_to_T_restart /= 0.0) .and. (us%s_to_T_restart /= us%s_to_T)) &
2280  i_t_rescale = us%s_to_T_restart / us%s_to_T
2281  z2_t_rescale = z_rescale**2*i_t_rescale
2282 
2283  if (z2_t_rescale /= 1.0) then
2284  if (associated(visc%Kd_shear)) then ; if (query_initialized(visc%Kd_shear, "Kd_shear", restart_cs)) then
2285  do k=1,nz+1 ; do j=js,je ; do i=is,ie
2286  visc%Kd_shear(i,j,k) = z2_t_rescale * visc%Kd_shear(i,j,k)
2287  enddo ; enddo ; enddo
2288  endif ; endif
2289 
2290  if (associated(visc%Kv_shear)) then ; if (query_initialized(visc%Kv_shear, "Kv_shear", restart_cs)) then
2291  do k=1,nz+1 ; do j=js,je ; do i=is,ie
2292  visc%Kv_shear(i,j,k) = z2_t_rescale * visc%Kv_shear(i,j,k)
2293  enddo ; enddo ; enddo
2294  endif ; endif
2295 
2296  if (associated(visc%Kv_shear_Bu)) then ; if (query_initialized(visc%Kv_shear_Bu, "Kv_shear_Bu", restart_cs)) then
2297  do k=1,nz+1 ; do j=js-1,je ; do i=is-1,ie
2298  visc%Kv_shear_Bu(i,j,k) = z2_t_rescale * visc%Kv_shear_Bu(i,j,k)
2299  enddo ; enddo ; enddo
2300  endif ; endif
2301  endif
2302 
2303  if (mle_use_pbl_mld .and. (z_rescale /= 1.0)) then
2304  if (associated(visc%MLD)) then ; if (query_initialized(visc%MLD, "MLD", restart_cs)) then
2305  do j=js,je ; do i=is,ie
2306  visc%MLD(i,j) = z_rescale * visc%MLD(i,j)
2307  enddo ; enddo
2308  endif ; endif
2309  endif
2310 

◆ set_visc_register_restarts()

subroutine, public mom_set_visc::set_visc_register_restarts ( type(hor_index_type), intent(in)  HI,
type(verticalgrid_type), intent(in)  GV,
type(param_file_type), intent(in)  param_file,
type(vertvisc_type), intent(inout)  visc,
type(mom_restart_cs), pointer  restart_CS 
)

Register any fields associated with the vertvisc_type.

Parameters
[in]hiA horizontal index type structure.
[in]gvThe ocean's vertical grid structure.
[in]param_fileA structure to parse for run-time parameters.
[in,out]viscA structure containing vertical viscosities and related fields. Allocated here.
restart_csA pointer to the restart control structure.

Definition at line 1887 of file MOM_set_viscosity.F90.

1887  type(hor_index_type), intent(in) :: HI !< A horizontal index type structure.
1888  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1889  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1890  !! parameters.
1891  type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical
1892  !! viscosities and related fields.
1893  !! Allocated here.
1894  type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure.
1895  ! Local variables
1896  logical :: use_kappa_shear, KS_at_vertex
1897  logical :: adiabatic, useKPP, useEPBL
1898  logical :: use_CVMix_shear, MLE_use_PBL_MLD, use_CVMix_conv
1899  integer :: isd, ied, jsd, jed, nz
1900  real :: hfreeze !< If hfreeze > 0 [m], melt potential will be computed.
1901  character(len=40) :: mdl = "MOM_set_visc" ! This module's name.
1902  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
1903 
1904  call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
1905  do_not_log=.true.)
1906 
1907  use_kappa_shear = .false. ; ks_at_vertex = .false. ; use_cvmix_shear = .false.
1908  usekpp = .false. ; useepbl = .false. ; use_cvmix_conv = .false.
1909 
1910  if (.not.adiabatic) then
1911  use_kappa_shear = kappa_shear_is_used(param_file)
1912  ks_at_vertex = kappa_shear_at_vertex(param_file)
1913  use_cvmix_shear = cvmix_shear_is_used(param_file)
1914  use_cvmix_conv = cvmix_conv_is_used(param_file)
1915  call get_param(param_file, mdl, "USE_KPP", usekpp, &
1916  "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
1917  "to calculate diffusivities and non-local transport in the OBL.", &
1918  default=.false., do_not_log=.true.)
1919  call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", useepbl, &
1920  "If true, use an implied energetics planetary boundary "//&
1921  "layer scheme to determine the diffusivity and viscosity "//&
1922  "in the surface boundary layer.", default=.false., do_not_log=.true.)
1923  endif
1924 
1925  if (use_kappa_shear .or. usekpp .or. useepbl .or. use_cvmix_shear .or. use_cvmix_conv) then
1926  call safe_alloc_ptr(visc%Kd_shear, isd, ied, jsd, jed, nz+1)
1927  call register_restart_field(visc%Kd_shear, "Kd_shear", .false., restart_cs, &
1928  "Shear-driven turbulent diffusivity at interfaces", "m2 s-1", z_grid='i')
1929  endif
1930  if (usekpp .or. useepbl .or. use_cvmix_shear .or. use_cvmix_conv .or. &
1931  (use_kappa_shear .and. .not.ks_at_vertex )) then
1932  call safe_alloc_ptr(visc%Kv_shear, isd, ied, jsd, jed, nz+1)
1933  call register_restart_field(visc%Kv_shear, "Kv_shear", .false., restart_cs, &
1934  "Shear-driven turbulent viscosity at interfaces", "m2 s-1", z_grid='i')
1935  endif
1936  if (use_kappa_shear .and. ks_at_vertex) then
1937  call safe_alloc_ptr(visc%TKE_turb, hi%IsdB, hi%IedB, hi%JsdB, hi%JedB, nz+1)
1938  call safe_alloc_ptr(visc%Kv_shear_Bu, hi%IsdB, hi%IedB, hi%JsdB, hi%JedB, nz+1)
1939  call register_restart_field(visc%Kv_shear_Bu, "Kv_shear_Bu", .false., restart_cs, &
1940  "Shear-driven turbulent viscosity at vertex interfaces", "m2 s-1", &
1941  hor_grid="Bu", z_grid='i')
1942  elseif (use_kappa_shear) then
1943  call safe_alloc_ptr(visc%TKE_turb, isd, ied, jsd, jed, nz+1)
1944  endif
1945 
1946  if (usekpp) then
1947  ! MOM_bkgnd_mixing uses Kv_slow when KPP is defined.
1948  call safe_alloc_ptr(visc%Kv_slow, isd, ied, jsd, jed, nz+1)
1949  endif
1950 
1951  ! visc%MLD is used to communicate the state of the (e)PBL or KPP to the rest of the model
1952  call get_param(param_file, mdl, "MLE_USE_PBL_MLD", mle_use_pbl_mld, &
1953  default=.false., do_not_log=.true.)
1954  ! visc%MLD needs to be allocated when melt potential is computed (HFREEZE>0)
1955  call get_param(param_file, mdl, "HFREEZE", hfreeze, &
1956  default=-1.0, do_not_log=.true.)
1957 
1958  if (mle_use_pbl_mld) then
1959  call safe_alloc_ptr(visc%MLD, isd, ied, jsd, jed)
1960  call register_restart_field(visc%MLD, "MLD", .false., restart_cs, &
1961  "Instantaneous active mixing layer depth", "m")
1962  endif
1963 
1964  if (hfreeze >= 0.0 .and. .not.mle_use_pbl_mld) then
1965  call safe_alloc_ptr(visc%MLD, isd, ied, jsd, jed)
1966  endif
1967 
1968 

◆ set_viscous_bbl()

subroutine, public mom_set_visc::set_viscous_bbl ( real, dimension(szib_(g),szj_(g),szk_(g)), intent(in)  u,
real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)  v,
real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
type(vertvisc_type), intent(inout)  visc,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(set_visc_cs), pointer  CS,
logical, intent(in), optional  symmetrize 
)

Calculates the thickness of the bottom boundary layer and the viscosity within that layer.

A drag law is used, either linearized about an assumed bottom velocity or using the actual near-bottom velocities combined with an assumed unresolved velocity. The bottom boundary layer thickness is limited by a combination of stratification and rotation, as in the paper of Killworth and Edwards, JPO 1999. It is not necessary to calculate the thickness and viscosity every time step; instead previous values may be used.

Viscous Bottom Boundary Layer

If set_visc_cs.bottomdraglaw is True then a bottom boundary layer viscosity and thickness are calculated so that the bottom stress is

\[ \mathbf{\tau}_b = C_d | U_{bbl} | \mathbf{u}_{bbl} \]

If set_visc_cs.bottomdraglaw is True then the term \(|U_{bbl}|\) is set equal to the value in set_visc_cs.drag_bg_vel so that \(C_d |U_{bbl}|\) becomes a Rayleigh bottom drag. Otherwise \(|U_{bbl}|\) is found by averaging the flow over the bottom set_visc_cs.hbbl of the model, adding the amplitude of tides set_visc_cs.tideamp and a constant set_visc_cs.drag_bg_vel. For these calculations the vertical grid at the velocity component locations is found by

\[ \begin{array}{ll} \frac{2 h^- h^+}{h^- + h^+} & u \left( h^+ - h^-\right) >= 0 \\ \frac{1}{2} \left( h^- + h^+ \right) & u \left( h^+ - h^-\right) < 0 \end{array} \]

which biases towards the thin cell if the thin cell is upwind. Biasing the grid toward thin upwind cells helps increase the effect of viscosity and inhibits flow out of these thin cells.

After diagnosing \(|U_{bbl}|\) over a fixed depth an active viscous boundary layer thickness is found using the ideas of Killworth and Edwards, 1999 (hereafter KW99). KW99 solve the equation

\[ \left( \frac{h_{bbl}}{h_f} \right)^2 + \frac{h_{bbl}}{h_N} = 1 \]

for the boundary layer depth \(h_{bbl}\). Here

\[ h_f = \frac{C_n u_*}{f} \]

is the rotation controlled boundary layer depth in the absence of stratification. \(u_*\) is the surface friction speed given by

\[ u_*^2 = C_d |U_{bbl}|^2 \]

and is a function of near bottom model flow.

\[ h_N = \frac{C_i u_*}{N} = \frac{ (C_i u_* )^2 }{g^\prime} \]

is the stratification controlled boundary layer depth. The non-dimensional parameters \(C_n=0.5\) and \(C_i=20\) are suggested by Zilitinkevich and Mironov, 1996.

If a Richardson number dependent mixing scheme is being used, as indicated by set_visc_cs.rino_mix, then the boundary layer thickness is bounded to be no larger than a half of set_visc_cs.hbbl .

Todo:
Channel drag needs to be explained

A BBL viscosity is calculated so that the no-slip boundary condition in the vertical viscosity solver implies the stress \(\mathbf{\tau}_b\).

References

  • Killworth, P. D., and N. R. Edwards, 1999: A Turbulent Bottom Boundary Layer Code for Use in Numerical Ocean Models. J. Phys. Oceanogr., 29, 1221-1238, doi:10.1175/1520-0485(1999)029<1221:ATBBLC>2.0.CO;2
  • Zilitinkevich, S., Mironov, D.V., 1996: A multi-limit formulation for the equilibrium depth of a stably stratified boundary layer. Boundary-Layer Meteorology 81, 325-351. doi:10.1007/BF02430334
Parameters
[in,out]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]uThe zonal velocity [L T-1 ~> m s-1].
[in]vThe meridional velocity [L T-1 ~> m s-1].
[in]hLayer thicknesses [H ~> m or kg m-2].
[in]tvA structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs..
[in,out]viscA structure containing vertical viscosities and related fields.
csThe control structure returned by a previous call to set_visc_init.
[in]symmetrizeIf present and true, do extra calculations of those values in visc that would be calculated with symmetric memory.

Definition at line 193 of file MOM_set_viscosity.F90.

193  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
194  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
195  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
196  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
197  intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
198  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
199  intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
200  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
201  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
202  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
203  !! available thermodynamic fields. Absent fields
204  !! have NULL ptrs..
205  type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
206  !! related fields.
207  type(set_visc_CS), pointer :: CS !< The control structure returned by a previous
208  !! call to set_visc_init.
209  logical, optional, intent(in) :: symmetrize !< If present and true, do extra calculations
210  !! of those values in visc that would be
211  !! calculated with symmetric memory.
212 
213  ! Local variables
214  real, dimension(SZIB_(G)) :: &
215  ustar, & ! The bottom friction velocity [Z T-1 ~> m s-1].
216  T_EOS, & ! The temperature used to calculate the partial derivatives
217  ! of density with T and S [degC].
218  s_eos, & ! The salinity used to calculate the partial derivatives
219  ! of density with T and S [ppt].
220  dr_dt, & ! Partial derivative of the density in the bottom boundary
221  ! layer with temperature [R degC-1 ~> kg m-3 degC-1].
222  dr_ds, & ! Partial derivative of the density in the bottom boundary
223  ! layer with salinity [R ppt-1 ~> kg m-3 ppt-1].
224  press ! The pressure at which dR_dT and dR_dS are evaluated [R L2 T-2 ~> Pa].
225  real :: htot ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
226  real :: htot_vel ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
227 
228  real :: Rhtot ! Running sum of thicknesses times the layer potential
229  ! densities [H R ~> kg m-2 or kg2 m-5].
230  real, dimension(SZIB_(G),SZJ_(G)) :: &
231  D_u, & ! Bottom depth interpolated to u points [Z ~> m].
232  mask_u ! A mask that disables any contributions from u points that
233  ! are land or past open boundary conditions [nondim], 0 or 1.
234  real, dimension(SZI_(G),SZJB_(G)) :: &
235  D_v, & ! Bottom depth interpolated to v points [Z ~> m].
236  mask_v ! A mask that disables any contributions from v points that
237  ! are land or past open boundary conditions [nondim], 0 or 1.
238  real, dimension(SZIB_(G),SZK_(G)) :: &
239  h_at_vel, & ! Layer thickness at a velocity point, using an upwind-biased
240  ! second order accurate estimate based on the previous velocity
241  ! direction [H ~> m or kg m-2].
242  h_vel, & ! Arithmetic mean of the layer thicknesses adjacent to a
243  ! velocity point [H ~> m or kg m-2].
244  t_vel, & ! Arithmetic mean of the layer temperatures adjacent to a
245  ! velocity point [degC].
246  s_vel, & ! Arithmetic mean of the layer salinities adjacent to a
247  ! velocity point [ppt].
248  rml_vel ! Arithmetic mean of the layer coordinate densities adjacent
249  ! to a velocity point [R ~> kg m-3].
250 
251  real :: h_vel_pos ! The arithmetic mean thickness at a velocity point
252  ! plus H_neglect to avoid 0 values [H ~> m or kg m-2].
253  real :: ustarsq ! 400 times the square of ustar, times
254  ! Rho0 divided by G_Earth and the conversion
255  ! from m to thickness units [H R ~> kg m-2 or kg2 m-5].
256  real :: cdrag_sqrt_Z ! Square root of the drag coefficient, times a unit conversion
257  ! factor from lateral lengths to vertical depths [Z L-1 ~> 1].
258  real :: cdrag_sqrt ! Square root of the drag coefficient [nondim].
259  real :: oldfn ! The integrated energy required to
260  ! entrain up to the bottom of the layer,
261  ! divided by G_Earth [H R ~> kg m-2 or kg2 m-5].
262  real :: Dfn ! The increment in oldfn for entraining
263  ! the layer [H R ~> kg m-2 or kg2 m-5].
264  real :: Dh ! The increment in layer thickness from
265  ! the present layer [H ~> m or kg m-2].
266  real :: bbl_thick ! The thickness of the bottom boundary layer [H ~> m or kg m-2].
267  real :: bbl_thick_Z ! The thickness of the bottom boundary layer [Z ~> m].
268  real :: kv_bbl ! The bottom boundary layer viscosity [Z2 T-1 ~> m2 s-1].
269  real :: C2f ! C2f = 2*f at velocity points [T-1 ~> s-1].
270 
271  real :: U_bg_sq ! The square of an assumed background
272  ! velocity, for calculating the mean
273  ! magnitude near the bottom for use in the
274  ! quadratic bottom drag [L2 T-2 ~> m2 s-2].
275  real :: hwtot ! Sum of the thicknesses used to calculate
276  ! the near-bottom velocity magnitude [H ~> m or kg m-2].
277  real :: hutot ! Running sum of thicknesses times the
278  ! velocity magnitudes [H T T-1 ~> m2 s-1 or kg m-1 s-1].
279  real :: Thtot ! Running sum of thickness times temperature [degC H ~> degC m or degC kg m-2].
280  real :: Shtot ! Running sum of thickness times salinity [ppt H ~> ppt m or ppt kg m-2].
281  real :: hweight ! The thickness of a layer that is within Hbbl
282  ! of the bottom [H ~> m or kg m-2].
283  real :: v_at_u, u_at_v ! v at a u point or vice versa [L T-1 ~> m s-1].
284  real :: Rho0x400_G ! 400*Rho0/G_Earth, times unit conversion factors
285  ! [R T2 H Z-2 ~> kg s2 m-4 or kg2 s2 m-7].
286  ! The 400 is a constant proposed by Killworth and Edwards, 1999.
287  real, dimension(SZI_(G),SZJ_(G),max(GV%nk_rho_varies,1)) :: &
288  Rml ! The mixed layer coordinate density [R ~> kg m-3].
289  real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate
290  ! density [R L2 T-2 ~> Pa] (usually set to 2e7 Pa = 2000 dbar).
291 
292  real :: D_vel ! The bottom depth at a velocity point [H ~> m or kg m-2].
293  real :: Dp, Dm ! The depths at the edges of a velocity cell [H ~> m or kg m-2].
294  real :: a ! a is the curvature of the bottom depth across a
295  ! cell, times the cell width squared [H ~> m or kg m-2].
296  real :: a_3, a_12 ! a/3 and a/12 [H ~> m or kg m-2].
297  real :: C24_a ! 24/a [H-1 ~> m-1 or m2 kg-1].
298  real :: slope ! The absolute value of the bottom depth slope across
299  ! a cell times the cell width [H ~> m or kg m-2].
300  real :: apb_4a, ax2_3apb ! Various nondimensional ratios of a and slope.
301  real :: a2x48_apb3, Iapb, Ibma_2 ! Combinations of a and slope [H-1 ~> m-1 or m2 kg-1].
302  ! All of the following "volumes" have units of thickness because they are normalized
303  ! by the full horizontal area of a velocity cell.
304  real :: Vol_open ! The cell volume above which it is open [H ~> m or kg m-2].
305  real :: Vol_direct ! With less than Vol_direct [H ~> m or kg m-2], there is a direct
306  ! solution of a cubic equation for L.
307  real :: Vol_2_reg ! The cell volume above which there are two separate
308  ! open areas that must be integrated [H ~> m or kg m-2].
309  real :: vol ! The volume below the interface whose normalized
310  ! width is being sought [H ~> m or kg m-2].
311  real :: vol_below ! The volume below the interface below the one that
312  ! is currently under consideration [H ~> m or kg m-2].
313  real :: Vol_err ! The error in the volume with the latest estimate of
314  ! L, or the error for the interface below [H ~> m or kg m-2].
315  real :: Vol_quit ! The volume error below which to quit iterating [H ~> m or kg m-2].
316  real :: Vol_tol ! A volume error tolerance [H ~> m or kg m-2].
317  real :: L(SZK_(G)+1) ! The fraction of the full cell width that is open at
318  ! the depth of each interface [nondim].
319  real :: L_direct ! The value of L above volume Vol_direct [nondim].
320  real :: L_max, L_min ! Upper and lower bounds on the correct value for L [nondim].
321  real :: Vol_err_max ! The volume errors for the upper and lower bounds on
322  real :: Vol_err_min ! the correct value for L [H ~> m or kg m-2].
323  real :: Vol_0 ! A deeper volume with known width L0 [H ~> m or kg m-2].
324  real :: L0 ! The value of L above volume Vol_0 [nondim].
325  real :: dVol ! vol - Vol_0 [H ~> m or kg m-2].
326  real :: dV_dL2 ! The partial derivative of volume with L squared
327  ! evaluated at L=L0 [H ~> m or kg m-2].
328  real :: h_neglect ! A thickness that is so small it is usually lost
329  ! in roundoff and can be neglected [H ~> m or kg m-2].
330  real :: ustH ! ustar converted to units of H T-1 [H T-1 ~> m s-1 or kg m-2 s-1].
331  real :: root ! A temporary variable [H T-1 ~> m s-1 or kg m-2 s-1].
332 
333  real :: Cell_width ! The transverse width of the velocity cell [L ~> m].
334  real :: Rayleigh ! A nondimensional value that is multiplied by the layer's
335  ! velocity magnitude to give the Rayleigh drag velocity, times
336  ! a lateral to vertical distance conversion factor [Z L-1 ~> 1].
337  real :: gam ! The ratio of the change in the open interface width
338  ! to the open interface width atop a cell [nondim].
339  real :: BBL_frac ! The fraction of a layer's drag that goes into the
340  ! viscous bottom boundary layer [nondim].
341  real :: BBL_visc_frac ! The fraction of all the drag that is expressed as
342  ! a viscous bottom boundary layer [nondim].
343  real, parameter :: C1_3 = 1.0/3.0, c1_6 = 1.0/6.0, c1_12 = 1.0/12.0
344  real :: C2pi_3 ! An irrational constant, 2/3 pi.
345  real :: tmp ! A temporary variable.
346  real :: tmp_val_m1_to_p1
347  real :: curv_tol ! Numerator of curvature cubed, used to estimate
348  ! accuracy of a single L(:) Newton iteration
349  logical :: use_L0, do_one_L_iter ! Control flags for L(:) Newton iteration
350  logical :: use_BBL_EOS, do_i(SZIB_(G))
351  integer, dimension(2) :: EOSdom ! The computational domain for the equation of state
352  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, m, n, K2, nkmb, nkml
353  integer :: itt, maxitt=20
354  type(ocean_OBC_type), pointer :: OBC => null()
355 
356  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
357  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
358  nkmb = gv%nk_rho_varies ; nkml = gv%nkml
359  h_neglect = gv%H_subroundoff
360  rho0x400_g = 400.0*(gv%Rho0 / (us%L_to_Z**2 * gv%g_Earth)) * gv%Z_to_H
361  vol_quit = 0.9*gv%Angstrom_H + h_neglect
362  c2pi_3 = 8.0*atan(1.0)/3.0
363 
364  if (.not.associated(cs)) call mom_error(fatal,"MOM_set_viscosity(BBL): "//&
365  "Module must be initialized before it is used.")
366  if (.not.cs%bottomdraglaw) return
367 
368  if (present(symmetrize)) then ; if (symmetrize) then
369  jsq = js-1 ; isq = is-1
370  endif ; endif
371 
372  if (cs%debug) then
373  call uvchksum("Start set_viscous_BBL [uv]", u, v, g%HI, haloshift=1, scale=us%L_T_to_m_s)
374  call hchksum(h,"Start set_viscous_BBL h", g%HI, haloshift=1, scale=gv%H_to_m)
375  if (associated(tv%T)) call hchksum(tv%T, "Start set_viscous_BBL T", g%HI, haloshift=1)
376  if (associated(tv%S)) call hchksum(tv%S, "Start set_viscous_BBL S", g%HI, haloshift=1)
377  endif
378 
379  use_bbl_eos = associated(tv%eqn_of_state) .and. cs%BBL_use_EOS
380  obc => cs%OBC
381 
382  u_bg_sq = cs%drag_bg_vel * cs%drag_bg_vel
383  cdrag_sqrt = sqrt(cs%cdrag)
384  cdrag_sqrt_z = us%L_to_Z * sqrt(cs%cdrag)
385  k2 = max(nkmb+1, 2)
386 
387 ! With a linear drag law, the friction velocity is already known.
388 ! if (CS%linear_drag) ustar(:) = cdrag_sqrt_Z*CS%drag_bg_vel
389 
390  if ((nkml>0) .and. .not.use_bbl_eos) then
391  eosdom(1) = isq - (g%isd-1) ; eosdom(2) = g%iec+1 - (g%isd-1)
392  do i=isq,ieq+1 ; p_ref(i) = tv%P_Ref ; enddo
393  !$OMP parallel do default(shared)
394  do k=1,nkmb ; do j=jsq,jeq+1
395  call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, rml(:,j,k), tv%eqn_of_state, &
396  eosdom)
397  enddo ; enddo
398  endif
399 
400  !$OMP parallel do default(shared)
401  do j=js-1,je ; do i=is-1,ie+1
402  d_v(i,j) = 0.5*(g%bathyT(i,j) + g%bathyT(i,j+1))
403  mask_v(i,j) = g%mask2dCv(i,j)
404  enddo ; enddo
405  !$OMP parallel do default(shared)
406  do j=js-1,je+1 ; do i=is-1,ie
407  d_u(i,j) = 0.5*(g%bathyT(i,j) + g%bathyT(i+1,j))
408  mask_u(i,j) = g%mask2dCu(i,j)
409  enddo ; enddo
410 
411  if (associated(obc)) then ; do n=1,obc%number_of_segments
412  if (.not. obc%segment(n)%on_pe) cycle
413  ! Use a one-sided projection of bottom depths at OBC points.
414  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
415  if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je)) then
416  do i = max(is-1,obc%segment(n)%HI%isd), min(ie+1,obc%segment(n)%HI%ied)
417  if (obc%segment(n)%direction == obc_direction_n) d_v(i,j) = g%bathyT(i,j)
418  if (obc%segment(n)%direction == obc_direction_s) d_v(i,j) = g%bathyT(i,j+1)
419  enddo
420  elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ie)) then
421  do j = max(js-1,obc%segment(n)%HI%jsd), min(je+1,obc%segment(n)%HI%jed)
422  if (obc%segment(n)%direction == obc_direction_e) d_u(i,j) = g%bathyT(i,j)
423  if (obc%segment(n)%direction == obc_direction_w) d_u(i,j) = g%bathyT(i+1,j)
424  enddo
425  endif
426  enddo ; endif
427  if (associated(obc)) then ; do n=1,obc%number_of_segments
428  ! Now project bottom depths across cell-corner points in the OBCs. The two
429  ! projections have to occur in sequence and can not be combined easily.
430  if (.not. obc%segment(n)%on_pe) cycle
431  ! Use a one-sided projection of bottom depths at OBC points.
432  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
433  if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je)) then
434  do i = max(is-1,obc%segment(n)%HI%IsdB), min(ie,obc%segment(n)%HI%IedB)
435  if (obc%segment(n)%direction == obc_direction_n) then
436  d_u(i,j+1) = d_u(i,j) ; mask_u(i,j+1) = 0.0
437  elseif (obc%segment(n)%direction == obc_direction_s) then
438  d_u(i,j) = d_u(i,j+1) ; mask_u(i,j) = 0.0
439  endif
440  enddo
441  elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ie)) then
442  do j = max(js-1,obc%segment(n)%HI%JsdB), min(je,obc%segment(n)%HI%JedB)
443  if (obc%segment(n)%direction == obc_direction_e) then
444  d_v(i+1,j) = d_v(i,j) ; mask_v(i+1,j) = 0.0
445  elseif (obc%segment(n)%direction == obc_direction_w) then
446  d_v(i,j) = d_v(i+1,j) ; mask_v(i,j) = 0.0
447  endif
448  enddo
449  endif
450  enddo ; endif
451 
452  if (.not.use_bbl_eos) rml_vel(:,:) = 0.0
453 
454  !$OMP parallel do default(private) shared(u,v,h,tv,visc,G,GV,US,CS,Rml,nz,nkmb, &
455  !$OMP nkml,Isq,Ieq,Jsq,Jeq,h_neglect,Rho0x400_G,C2pi_3, &
456  !$OMP U_bg_sq,cdrag_sqrt_Z,cdrag_sqrt,K2,use_BBL_EOS, &
457  !$OMP OBC,maxitt,D_u,D_v,mask_u,mask_v) &
458  !$OMP firstprivate(Vol_quit)
459  do j=jsq,jeq ; do m=1,2
460 
461  if (m==1) then
462  ! m=1 refers to u-points
463  if (j<g%Jsc) cycle
464  is = isq ; ie = ieq
465  do i=is,ie
466  do_i(i) = .false.
467  if (g%mask2dCu(i,j) > 0) do_i(i) = .true.
468  enddo
469  else
470  ! m=2 refers to v-points
471  is = g%isc ; ie = g%iec
472  do i=is,ie
473  do_i(i) = .false.
474  if (g%mask2dCv(i,j) > 0) do_i(i) = .true.
475  enddo
476  endif
477 
478  ! Calculate thickness at velocity points (u or v depending on value of m).
479  ! Also interpolate the ML density or T/S properties.
480  if (m==1) then ! u-points
481  do k=1,nz ; do i=is,ie
482  if (do_i(i)) then
483  if (u(i,j,k) *(h(i+1,j,k) - h(i,j,k)) >= 0) then
484  ! If the flow is from thin to thick then bias towards the thinner thickness
485  h_at_vel(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / &
486  (h(i,j,k) + h(i+1,j,k) + h_neglect)
487  else
488  ! If the flow is from thick to thin then use the simple average thickness
489  h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
490  endif
491  endif
492  h_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
493  enddo ; enddo
494  if (use_bbl_eos) then ; do k=1,nz ; do i=is,ie
495  ! Perhaps these should be thickness weighted.
496  t_vel(i,k) = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
497  s_vel(i,k) = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
498  enddo ; enddo ; else ; do k=1,nkmb ; do i=is,ie
499  rml_vel(i,k) = 0.5 * (rml(i,j,k) + rml(i+1,j,k))
500  enddo ; enddo ; endif
501  else ! v-points
502  do k=1,nz ; do i=is,ie
503  if (do_i(i)) then
504  if (v(i,j,k) * (h(i,j+1,k) - h(i,j,k)) >= 0) then
505  ! If the flow is from thin to thick then bias towards the thinner thickness
506  h_at_vel(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / &
507  (h(i,j,k) + h(i,j+1,k) + h_neglect)
508  else
509  ! If the flow is from thick to thin then use the simple average thickness
510  h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
511  endif
512  endif
513  h_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
514  enddo ; enddo
515  if (use_bbl_eos) then ; do k=1,nz ; do i=is,ie
516  t_vel(i,k) = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
517  s_vel(i,k) = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
518  enddo ; enddo ; else ; do k=1,nkmb ; do i=is,ie
519  rml_vel(i,k) = 0.5 * (rml(i,j,k) + rml(i,j+1,k))
520  enddo ; enddo ; endif
521  endif
522 
523  if (associated(obc)) then ; if (obc%number_of_segments > 0) then
524  ! Apply a zero gradient projection of thickness across OBC points.
525  if (m==1) then
526  do i=is,ie ; if (do_i(i) .and. (obc%segnum_u(i,j) /= obc_none)) then
527  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
528  do k=1,nz
529  h_at_vel(i,k) = h(i,j,k) ; h_vel(i,k) = h(i,j,k)
530  enddo
531  if (use_bbl_eos) then
532  do k=1,nz
533  t_vel(i,k) = tv%T(i,j,k) ; s_vel(i,k) = tv%S(i,j,k)
534  enddo
535  else
536  do k=1,nkmb
537  rml_vel(i,k) = rml(i,j,k)
538  enddo
539  endif
540  elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
541  do k=1,nz
542  h_at_vel(i,k) = h(i+1,j,k) ; h_vel(i,k) = h(i+1,j,k)
543  enddo
544  if (use_bbl_eos) then
545  do k=1,nz
546  t_vel(i,k) = tv%T(i+1,j,k) ; s_vel(i,k) = tv%S(i+1,j,k)
547  enddo
548  else
549  do k=1,nkmb
550  rml_vel(i,k) = rml(i+1,j,k)
551  enddo
552  endif
553  endif
554  endif ; enddo
555  else
556  do i=is,ie ; if (do_i(i) .and. (obc%segnum_v(i,j) /= obc_none)) then
557  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
558  do k=1,nz
559  h_at_vel(i,k) = h(i,j,k) ; h_vel(i,k) = h(i,j,k)
560  enddo
561  if (use_bbl_eos) then
562  do k=1,nz
563  t_vel(i,k) = tv%T(i,j,k) ; s_vel(i,k) = tv%S(i,j,k)
564  enddo
565  else
566  do k=1,nkmb
567  rml_vel(i,k) = rml(i,j,k)
568  enddo
569  endif
570  elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
571  do k=1,nz
572  h_at_vel(i,k) = h(i,j+1,k) ; h_vel(i,k) = h(i,j+1,k)
573  enddo
574  if (use_bbl_eos) then
575  do k=1,nz
576  t_vel(i,k) = tv%T(i,j+1,k) ; s_vel(i,k) = tv%S(i,j+1,k)
577  enddo
578  else
579  do k=1,nkmb
580  rml_vel(i,k) = rml(i,j+1,k)
581  enddo
582  endif
583  endif
584  endif ; enddo
585  endif
586  endif ; endif
587 
588  if (use_bbl_eos .or. .not.cs%linear_drag) then
589  ! Calculate the mean velocity magnitude over the bottommost CS%Hbbl of
590  ! the water column for determining the quadratic bottom drag.
591  ! Used in ustar(i)
592  do i=is,ie ; if (do_i(i)) then
593  htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
594  thtot = 0.0 ; shtot = 0.0
595  do k=nz,1,-1
596 
597  if (htot_vel>=cs%Hbbl) exit ! terminate the k loop
598 
599  hweight = min(cs%Hbbl - htot_vel, h_at_vel(i,k))
600  if (hweight < 1.5*gv%Angstrom_H + h_neglect) cycle
601 
602  htot_vel = htot_vel + h_at_vel(i,k)
603  hwtot = hwtot + hweight
604 
605  if ((.not.cs%linear_drag) .and. (hweight >= 0.0)) then ; if (m==1) then
606  v_at_u = set_v_at_u(v, h, g, i, j, k, mask_v, obc)
607  if (cs%BBL_use_tidal_bg) then
608  u_bg_sq = 0.5*( g%mask2dT(i,j)*(cs%tideamp(i,j)*cs%tideamp(i,j))+ &
609  g%mask2dT(i+1,j)*(cs%tideamp(i+1,j)*cs%tideamp(i+1,j)) )
610  endif
611  hutot = hutot + hweight * sqrt(u(i,j,k)*u(i,j,k) + &
612  v_at_u*v_at_u + u_bg_sq)
613  else
614  u_at_v = set_u_at_v(u, h, g, i, j, k, mask_u, obc)
615  if (cs%BBL_use_tidal_bg) then
616  u_bg_sq = 0.5*( g%mask2dT(i,j)*(cs%tideamp(i,j)*cs%tideamp(i,j))+ &
617  g%mask2dT(i,j+1)*(cs%tideamp(i,j+1)*cs%tideamp(i,j+1)) )
618  endif
619  hutot = hutot + hweight * sqrt(v(i,j,k)*v(i,j,k) + &
620  u_at_v*u_at_v + u_bg_sq)
621  endif ; endif
622 
623  if (use_bbl_eos .and. (hweight >= 0.0)) then
624  thtot = thtot + hweight * t_vel(i,k)
625  shtot = shtot + hweight * s_vel(i,k)
626  endif
627  enddo ! end of k loop
628 
629  ! Set u* based on u*^2 = Cdrag u_bbl^2
630  if (.not.cs%linear_drag .and. (hwtot > 0.0)) then
631  ustar(i) = cdrag_sqrt_z*hutot/hwtot
632  else
633  ustar(i) = cdrag_sqrt_z*cs%drag_bg_vel
634  endif
635 
636  if (use_bbl_eos) then ; if (hwtot > 0.0) then
637  t_eos(i) = thtot/hwtot ; s_eos(i) = shtot/hwtot
638  else
639  t_eos(i) = 0.0 ; s_eos(i) = 0.0
640  endif ; endif
641 
642  ! Diagnostic BBL flow speed at u- and v-points.
643  if (cs%id_bbl_u>0 .and. m==1) then
644  if (hwtot > 0.0) cs%bbl_u(i,j) = hutot/hwtot
645  elseif (cs%id_bbl_v>0 .and. m==2) then
646  if (hwtot > 0.0) cs%bbl_v(i,j) = hutot/hwtot
647  endif
648 
649  endif ; enddo
650  else
651  do i=is,ie ; ustar(i) = cdrag_sqrt_z*cs%drag_bg_vel ; enddo
652  endif ! Not linear_drag
653 
654  if (use_bbl_eos) then
655  if (associated(tv%p_surf)) then
656  if (m==1) then ; do i=is,ie ; press(i) = 0.5*(tv%p_surf(i,j) + tv%p_surf(i+1,j)) ; enddo
657  else ; do i=is,ie ; press(i) = 0.5*(tv%p_surf(i,j) + tv%p_surf(i,j+1)) ; enddo ; endif
658  else
659  do i=is,ie ; press(i) = 0.0 ; enddo
660  endif
661  do i=is,ie ; if (.not.do_i(i)) then ; t_eos(i) = 0.0 ; s_eos(i) = 0.0 ; endif ; enddo
662  do k=1,nz ; do i=is,ie
663  press(i) = press(i) + (gv%H_to_RZ*gv%g_Earth) * h_vel(i,k)
664  enddo ; enddo
665  call calculate_density_derivs(t_eos, s_eos, press, dr_dt, dr_ds, tv%eqn_of_state, &
666  (/is-g%IsdB+1,ie-g%IsdB+1/) )
667  endif
668 
669  ! Find a BBL thickness given by equation 2.20 of Killworth and Edwards, 1999:
670  ! ( f h / Cn u* )^2 + ( N h / Ci u* ) = 1
671  ! where Cn=0.5 and Ci=20 (constants suggested by Zilitinkevich and Mironov, 1996).
672  ! Eq. 2.20 can be expressed in terms of boundary layer thicknesses limited by
673  ! rotation (h_f) and stratification (h_N):
674  ! ( h / h_f )^2 + ( h / h_N ) = 1
675  ! When stratification dominates h_N<<h_f, and vice versa.
676  do i=is,ie ; if (do_i(i)) then
677  ! The 400.0 in this expression is the square of a Ci introduced in KW99, eq. 2.22.
678  ustarsq = rho0x400_g * ustar(i)**2 ! Note not in units of u*^2 but [H R]
679  htot = 0.0
680 
681  ! Calculate the thickness of a stratification limited BBL ignoring rotation:
682  ! h_N = Ci u* / N (limit of KW99 eq. 2.20 for |f|->0)
683  ! For layer mode, N^2 = g'/h. Since (Ci u*)^2 = (h_N N)^2 = h_N g' then
684  ! h_N = (Ci u*)^2 / g' (KW99, eq, 2.22)
685  ! Starting from the bottom, integrate the stratification upward until h_N N balances Ci u*
686  ! or in layer mode
687  ! h_N Delta rho ~ (Ci u*)^2 rho0 / g
688  ! where the rhs is stored in variable ustarsq.
689  ! The method was described in Stephens and Hallberg 2000 (unpublished and lost manuscript).
690  if (use_bbl_eos) then
691  thtot = 0.0 ; shtot = 0.0 ; oldfn = 0.0
692  do k=nz,2,-1
693  if (h_at_vel(i,k) <= 0.0) cycle
694 
695  ! Delta rho * h_bbl assuming everything below is homogenized
696  oldfn = dr_dt(i)*(thtot - t_vel(i,k)*htot) + &
697  dr_ds(i)*(shtot - s_vel(i,k)*htot)
698  if (oldfn >= ustarsq) exit
699 
700  ! Local Delta rho * h_bbl at interface
701  dfn = (dr_dt(i)*(t_vel(i,k) - t_vel(i,k-1)) + &
702  dr_ds(i)*(s_vel(i,k) - s_vel(i,k-1))) * &
703  (h_at_vel(i,k) + htot)
704 
705  if ((oldfn + dfn) <= ustarsq) then
706  ! Use whole layer
707  dh = h_at_vel(i,k)
708  else
709  ! Use only part of the layer
710  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
711  endif
712 
713  ! Increment total BBL thickness and cumulative T and S
714  htot = htot + dh
715  thtot = thtot + t_vel(i,k)*dh ; shtot = shtot + s_vel(i,k)*dh
716  enddo
717  if ((oldfn < ustarsq) .and. h_at_vel(i,1) > 0.0) then
718  ! Layer 1 might be part of the BBL.
719  if (dr_dt(i) * (thtot - t_vel(i,1)*htot) + &
720  dr_ds(i) * (shtot - s_vel(i,1)*htot) < ustarsq) &
721  htot = htot + h_at_vel(i,1)
722  endif ! Examination of layer 1.
723  else ! Use Rlay and/or the coordinate density as density variables.
724  rhtot = 0.0
725  do k=nz,k2,-1
726  oldfn = rhtot - gv%Rlay(k)*htot
727  dfn = (gv%Rlay(k) - gv%Rlay(k-1))*(h_at_vel(i,k)+htot)
728 
729  if (oldfn >= ustarsq) then
730  cycle
731  elseif ((oldfn + dfn) <= ustarsq) then
732  dh = h_at_vel(i,k)
733  else
734  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
735  endif
736 
737  htot = htot + dh
738  rhtot = rhtot + gv%Rlay(k)*dh
739  enddo
740  if (nkml>0) then
741  do k=nkmb,2,-1
742  oldfn = rhtot - rml_vel(i,k)*htot
743  dfn = (rml_vel(i,k) - rml_vel(i,k-1)) * (h_at_vel(i,k)+htot)
744 
745  if (oldfn >= ustarsq) then
746  cycle
747  elseif ((oldfn + dfn) <= ustarsq) then
748  dh = h_at_vel(i,k)
749  else
750  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
751  endif
752 
753  htot = htot + dh
754  rhtot = rhtot + rml_vel(i,k)*dh
755  enddo
756  if (rhtot - rml_vel(i,1)*htot < ustarsq) htot = htot + h_at_vel(i,1)
757  else
758  if (rhtot - gv%Rlay(1)*htot < ustarsq) htot = htot + h_at_vel(i,1)
759  endif
760  endif ! use_BBL_EOS
761 
762  ! Value of 2*f at u- or v-points.
763  if (m==1) then ; c2f = g%CoriolisBu(i,j-1) + g%CoriolisBu(i,j)
764  else ; c2f = g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j) ; endif
765 
766  ! The thickness of a rotation limited BBL ignoring stratification is
767  ! h_f ~ Cn u* / f (limit of KW99 eq. 2.20 for N->0).
768  ! The buoyancy limit of BBL thickness (h_N) is already in the variable htot from above.
769  ! Substituting x = h_N/h into KW99 eq. 2.20 yields the quadratic
770  ! x^2 - x = (h_N / h_f)^2
771  ! for which the positive root is
772  ! xp = 1/2 + sqrt( 1/4 + (h_N/h_f)^2 )
773  ! and thus h_bbl = h_N / xp . Since h_f = Cn u*/f and Cn=0.5
774  ! xp = 1/2 + sqrt( 1/4 + (2 f h_N/u*)^2 )
775  ! To avoid dividing by zero if u*=0 then
776  ! xp u* = 1/2 u* + sqrt( 1/4 u*^2 + (2 f h_N)^2 )
777  if (cs%cdrag * u_bg_sq <= 0.0) then
778  ! This avoids NaNs and overflows, and could be used in all cases,
779  ! but is not bitwise identical to the current code.
780  usth = ustar(i)*gv%Z_to_H ; root = sqrt(0.25*usth**2 + (htot*c2f)**2)
781  if (htot*usth <= (cs%BBL_thick_min+h_neglect) * (0.5*usth + root)) then
782  bbl_thick = cs%BBL_thick_min
783  else
784  ! The following expression reads
785  ! h_bbl = h_N u* / ( 1/2 u* + sqrt( 1/4 u*^2 + ( 2 f h_N )^2 ) )
786  ! which is h_bbl = h_N u*/(xp u*) as described above.
787  bbl_thick = (htot * usth) / (0.5*usth + root)
788  endif
789  else
790  ! The following expression reads
791  ! h_bbl = h_N / ( 1/2 + sqrt( 1/4 + ( 2 f h_N / u* )^2 ) )
792  ! which is h_bbl = h_N/xp as described above.
793  bbl_thick = htot / (0.5 + sqrt(0.25 + htot*htot*c2f*c2f/ &
794  ((ustar(i)*ustar(i)) * (gv%Z_to_H**2)) ) )
795 
796  if (bbl_thick < cs%BBL_thick_min) bbl_thick = cs%BBL_thick_min
797  endif
798 
799  ! If there is Richardson number dependent mixing, that determines
800  ! the vertical extent of the bottom boundary layer, and there is no
801  ! need to set that scale here. In fact, viscously reducing the
802  ! shears over an excessively large region reduces the efficacy of
803  ! the Richardson number dependent mixing.
804  ! In other words, if using RiNo_mix then CS%Hbbl acts as an upper bound on
805  ! bbl_thick.
806  if ((bbl_thick > 0.5*cs%Hbbl) .and. (cs%RiNo_mix)) bbl_thick = 0.5*cs%Hbbl
807 
808  if (cs%Channel_drag) then
809  ! The drag within the bottommost bbl_thick is applied as a part of
810  ! an enhanced bottom viscosity, while above this the drag is applied
811  ! directly to the layers in question as a Rayleigh drag term.
812  if (m==1) then
813  d_vel = d_u(i,j)
814  tmp = g%mask2dCu(i,j+1) * d_u(i,j+1)
815  dp = 2.0 * d_vel * tmp / (d_vel + tmp)
816  tmp = g%mask2dCu(i,j-1) * d_u(i,j-1)
817  dm = 2.0 * d_vel * tmp / (d_vel + tmp)
818  else
819  d_vel = d_v(i,j)
820  tmp = g%mask2dCv(i+1,j) * d_v(i+1,j)
821  dp = 2.0 * d_vel * tmp / (d_vel + tmp)
822  tmp = g%mask2dCv(i-1,j) * d_v(i-1,j)
823  dm = 2.0 * d_vel * tmp / (d_vel + tmp)
824  endif
825  if (dm > dp) then ; tmp = dp ; dp = dm ; dm = tmp ; endif
826 
827  ! Convert the D's to the units of thickness.
828  dp = gv%Z_to_H*dp ; dm = gv%Z_to_H*dm ; d_vel = gv%Z_to_H*d_vel
829 
830  a_3 = (dp + dm - 2.0*d_vel) ; a = 3.0*a_3 ; a_12 = 0.25*a_3
831  slope = dp - dm
832  ! If the curvature is small enough, there is no reason not to assume
833  ! a uniformly sloping or flat bottom.
834  if (abs(a) < 1e-2*(slope + cs%BBL_thick_min)) a = 0.0
835  ! Each cell extends from x=-1/2 to 1/2, and has a topography
836  ! given by D(x) = a*x^2 + b*x + D - a/12.
837 
838  ! Calculate the volume above which the entire cell is open and the
839  ! other volumes at which the equation that is solved for L changes.
840  if (a > 0.0) then
841  if (slope >= a) then
842  vol_open = d_vel - dm ; vol_2_reg = vol_open
843  else
844  tmp = slope/a
845  vol_open = 0.25*slope*tmp + c1_12*a
846  vol_2_reg = 0.5*tmp**2 * (a - c1_3*slope)
847  endif
848  ! Define some combinations of a & b for later use.
849  c24_a = 24.0/a ; iapb = 1.0/(a+slope)
850  apb_4a = (slope+a)/(4.0*a) ; a2x48_apb3 = (48.0*(a*a))*(iapb**3)
851  ax2_3apb = 2.0*c1_3*a*iapb
852  elseif (a == 0.0) then
853  vol_open = 0.5*slope
854  if (slope > 0) iapb = 1.0/slope
855  else ! a < 0.0
856  vol_open = d_vel - dm
857  if (slope >= -a) then
858  iapb = 1.0e30 ; if (slope+a /= 0.0) iapb = 1.0/(a+slope)
859  vol_direct = 0.0 ; l_direct = 0.0 ; c24_a = 0.0
860  else
861  c24_a = 24.0/a ; iapb = 1.0/(a+slope)
862  l_direct = 1.0 + slope/a ! L_direct < 1 because a < 0
863  vol_direct = -c1_6*a*l_direct**3
864  endif
865  ibma_2 = 2.0 / (slope - a)
866  endif
867 
868  l(nz+1) = 0.0 ; vol = 0.0 ; vol_err = 0.0 ; bbl_visc_frac = 0.0
869  ! Determine the normalized open length at each interface.
870  do k=nz,1,-1
871  vol_below = vol
872 
873  vol = vol + h_vel(i,k)
874  h_vel_pos = h_vel(i,k) + h_neglect
875 
876  if (vol >= vol_open) then ; l(k) = 1.0
877  elseif (a == 0) then ! The bottom has no curvature.
878  l(k) = sqrt(2.0*vol*iapb)
879  elseif (a > 0) then
880  ! There may be a minimum depth, and there are
881  ! analytic expressions for L for all cases.
882  if (vol < vol_2_reg) then
883  ! In this case, there is a contiguous open region and
884  ! vol = 0.5*L^2*(slope + a/3*(3-4L)).
885  if (a2x48_apb3*vol < 1e-8) then ! Could be 1e-7?
886  ! There is a very good approximation here for massless layers.
887  l0 = sqrt(2.0*vol*iapb) ; l(k) = l0*(1.0 + ax2_3apb*l0)
888  else
889  l(k) = apb_4a * (1.0 - &
890  2.0 * cos(c1_3*acos(a2x48_apb3*vol - 1.0) - c2pi_3))
891  endif
892  ! To check the answers.
893  ! Vol_err = 0.5*(L(K)*L(K))*(slope + a_3*(3.0-4.0*L(K))) - vol
894  else ! There are two separate open regions.
895  ! vol = slope^2/4a + a/12 - (a/12)*(1-L)^2*(1+2L)
896  ! At the deepest volume, L = slope/a, at the top L = 1.
897  !L(K) = 0.5 - cos(C1_3*acos(1.0 - C24_a*(Vol_open - vol)) - C2pi_3)
898  tmp_val_m1_to_p1 = 1.0 - c24_a*(vol_open - vol)
899  tmp_val_m1_to_p1 = max(-1., min(1., tmp_val_m1_to_p1))
900  l(k) = 0.5 - cos(c1_3*acos(tmp_val_m1_to_p1) - c2pi_3)
901  ! To check the answers.
902  ! Vol_err = Vol_open - a_12*(1.0+2.0*L(K)) * (1.0-L(K))**2 - vol
903  endif
904  else ! a < 0.
905  if (vol <= vol_direct) then
906  ! Both edges of the cell are bounded by walls.
907  l(k) = (-0.25*c24_a*vol)**c1_3
908  else
909  ! x_R is at 1/2 but x_L is in the interior & L is found by solving
910  ! vol = 0.5*L^2*(slope + a/3*(3-4L))
911 
912  ! Vol_err = 0.5*(L(K+1)*L(K+1))*(slope + a_3*(3.0-4.0*L(K+1))) - vol_below
913  ! Change to ...
914  ! if (min(Vol_below + Vol_err, vol) <= Vol_direct) then ?
915  if (vol_below + vol_err <= vol_direct) then
916  l0 = l_direct ; vol_0 = vol_direct
917  else
918  l0 = l(k+1) ; vol_0 = vol_below + vol_err
919  ! Change to Vol_0 = min(Vol_below + Vol_err, vol) ?
920  endif
921 
922  ! Try a relatively simple solution that usually works well
923  ! for massless layers.
924  dv_dl2 = 0.5*(slope+a) - a*l0 ; dvol = (vol-vol_0)
925  ! dV_dL2 = 0.5*(slope+a) - a*L0 ; dVol = max(vol-Vol_0, 0.0)
926 
927  use_l0 = .false.
928  do_one_l_iter = .false.
929  if (cs%answers_2018) then
930  curv_tol = gv%Angstrom_H*dv_dl2**2 &
931  * (0.25 * dv_dl2 * gv%Angstrom_H - a * l0 * dvol)
932  do_one_l_iter = (a * a * dvol**3) < curv_tol
933  else
934  ! The following code is more robust when GV%Angstrom_H=0, but
935  ! it changes answers.
936  use_l0 = (dvol <= 0.)
937 
938  vol_tol = max(0.5 * gv%Angstrom_H + gv%H_subroundoff, 1e-14 * vol)
939  vol_quit = max(0.9 * gv%Angstrom_H + gv%H_subroundoff, 1e-14 * vol)
940 
941  curv_tol = vol_tol * dv_dl2**2 &
942  * (dv_dl2 * vol_tol - 2.0 * a * l0 * dvol)
943  do_one_l_iter = (a * a * dvol**3) < curv_tol
944  endif
945 
946  if (use_l0) then
947  l(k) = l0
948  vol_err = 0.5*(l(k)*l(k))*(slope + a_3*(3.0-4.0*l(k))) - vol
949  elseif (do_one_l_iter) then
950  ! One iteration of Newton's method should give an estimate
951  ! that is accurate to within Vol_tol.
952  l(k) = sqrt(l0*l0 + dvol / dv_dl2)
953  vol_err = 0.5*(l(k)*l(k))*(slope + a_3*(3.0-4.0*l(k))) - vol
954  else
955  if (dv_dl2*(1.0-l0*l0) < dvol + &
956  dv_dl2 * (vol_open - vol)*ibma_2) then
957  l_max = sqrt(1.0 - (vol_open - vol)*ibma_2)
958  else
959  l_max = sqrt(l0*l0 + dvol / dv_dl2)
960  endif
961  l_min = sqrt(l0*l0 + dvol / (0.5*(slope+a) - a*l_max))
962 
963  vol_err_min = 0.5*(l_min**2)*(slope + a_3*(3.0-4.0*l_min)) - vol
964  vol_err_max = 0.5*(l_max**2)*(slope + a_3*(3.0-4.0*l_max)) - vol
965  ! if ((abs(Vol_err_min) <= Vol_quit) .or. (Vol_err_min >= Vol_err_max)) then
966  if (abs(vol_err_min) <= vol_quit) then
967  l(k) = l_min ; vol_err = vol_err_min
968  else
969  l(k) = sqrt((l_min**2*vol_err_max - l_max**2*vol_err_min) / &
970  (vol_err_max - vol_err_min))
971  do itt=1,maxitt
972  vol_err = 0.5*(l(k)*l(k))*(slope + a_3*(3.0-4.0*l(k))) - vol
973  if (abs(vol_err) <= vol_quit) exit
974  ! Take a Newton's method iteration. This equation has proven
975  ! robust enough not to need bracketing.
976  l(k) = l(k) - vol_err / (l(k)* (slope + a - 2.0*a*l(k)))
977  ! This would be a Newton's method iteration for L^2:
978  ! L(K) = sqrt(L(K)*L(K) - Vol_err / (0.5*(slope+a) - a*L(K)))
979  enddo
980  endif ! end of iterative solver
981  endif ! end of 1-boundary alternatives.
982  endif ! end of a<0 cases.
983  endif
984 
985  ! Determine the drag contributing to the bottom boundary layer
986  ! and the Raleigh drag that acts on each layer.
987  if (l(k) > l(k+1)) then
988  if (vol_below < bbl_thick) then
989  bbl_frac = (1.0-vol_below/bbl_thick)**2
990  bbl_visc_frac = bbl_visc_frac + bbl_frac*(l(k) - l(k+1))
991  else
992  bbl_frac = 0.0
993  endif
994 
995  if (m==1) then ; cell_width = g%dy_Cu(i,j)
996  else ; cell_width = g%dx_Cv(i,j) ; endif
997  gam = 1.0 - l(k+1)/l(k)
998  rayleigh = us%L_to_Z * cs%cdrag * (l(k)-l(k+1)) * (1.0-bbl_frac) * &
999  (12.0*cs%c_Smag*h_vel_pos) / (12.0*cs%c_Smag*h_vel_pos + &
1000  us%L_to_Z*gv%Z_to_H * cs%cdrag * gam*(1.0-gam)*(1.0-1.5*gam) * l(k)**2 * cell_width)
1001  else ! This layer feels no drag.
1002  rayleigh = 0.0
1003  endif
1004 
1005  if (m==1) then
1006  if (rayleigh > 0.0) then
1007  v_at_u = set_v_at_u(v, h, g, i, j, k, mask_v, obc)
1008  visc%Ray_u(i,j,k) = rayleigh*sqrt(u(i,j,k)*u(i,j,k) + &
1009  v_at_u*v_at_u + u_bg_sq)
1010  else ; visc%Ray_u(i,j,k) = 0.0 ; endif
1011  else
1012  if (rayleigh > 0.0) then
1013  u_at_v = set_u_at_v(u, h, g, i, j, k, mask_u, obc)
1014  visc%Ray_v(i,j,k) = rayleigh*sqrt(v(i,j,k)*v(i,j,k) + &
1015  u_at_v*u_at_v + u_bg_sq)
1016  else ; visc%Ray_v(i,j,k) = 0.0 ; endif
1017  endif
1018 
1019  enddo ! k loop to determine L(K).
1020 
1021  ! Set the near-bottom viscosity to a value which will give
1022  ! the correct stress when the shear occurs over bbl_thick.
1023  ! See next block for explanation.
1024  bbl_thick_z = bbl_thick * gv%H_to_Z
1025  if (cs%correct_BBL_bounds .and. &
1026  cdrag_sqrt*ustar(i)*bbl_thick_z*bbl_visc_frac <= cs%Kv_BBL_min) then
1027  ! If the bottom stress implies less viscosity than Kv_BBL_min then
1028  ! set kv_bbl to the bound and recompute bbl_thick to be consistent
1029  ! but with a ridiculously large upper bound on thickness (for Cd u*=0)
1030  kv_bbl = cs%Kv_BBL_min
1031  if (cdrag_sqrt*ustar(i)*bbl_visc_frac*g%Rad_Earth*us%m_to_Z > kv_bbl) then
1032  bbl_thick_z = kv_bbl / ( cdrag_sqrt*ustar(i)*bbl_visc_frac )
1033  else
1034  bbl_thick_z = g%Rad_Earth * us%m_to_Z
1035  endif
1036  else
1037  kv_bbl = cdrag_sqrt*ustar(i)*bbl_thick_z*bbl_visc_frac
1038  endif
1039 
1040  else ! Not Channel_drag.
1041  ! Set the near-bottom viscosity to a value which will give
1042  ! the correct stress when the shear occurs over bbl_thick.
1043  ! - The bottom stress is tau_b = Cdrag * u_bbl^2
1044  ! - u_bbl was calculated by averaging flow over CS%Hbbl
1045  ! (and includes unresolved tidal components)
1046  ! - u_bbl is embedded in u* since u*^2 = Cdrag u_bbl^2
1047  ! - The average shear in the BBL is du/dz = 2 * u_bbl / h_bbl
1048  ! (which assumes a linear profile, hence the "2")
1049  ! - bbl_thick was bounded to <= 0.5 * CS%Hbbl
1050  ! - The viscous stress kv_bbl du/dz should balance tau_b
1051  ! Cdrag u_bbl^2 = kv_bbl du/dz
1052  ! = 2 kv_bbl u_bbl
1053  ! so
1054  ! kv_bbl = 0.5 h_bbl Cdrag u_bbl
1055  ! = 0.5 h_bbl sqrt(Cdrag) u*
1056  bbl_thick_z = bbl_thick * gv%H_to_Z
1057  if (cs%correct_BBL_bounds .and. &
1058  cdrag_sqrt*ustar(i)*bbl_thick_z <= cs%Kv_BBL_min) then
1059  ! If the bottom stress implies less viscosity than Kv_BBL_min then
1060  ! set kv_bbl to the bound and recompute bbl_thick to be consistent
1061  ! but with a ridiculously large upper bound on thickness (for Cd u*=0)
1062  kv_bbl = cs%Kv_BBL_min
1063  if (cdrag_sqrt*ustar(i)*g%Rad_Earth*us%m_to_Z > kv_bbl) then
1064  bbl_thick_z = kv_bbl / ( cdrag_sqrt*ustar(i) )
1065  else
1066  bbl_thick_z = g%Rad_Earth * us%m_to_Z
1067  endif
1068  else
1069  kv_bbl = cdrag_sqrt*ustar(i)*bbl_thick_z
1070  endif
1071  endif
1072  kv_bbl = max(cs%Kv_BBL_min, kv_bbl)
1073  if (m==1) then
1074  visc%Kv_bbl_u(i,j) = kv_bbl
1075  visc%bbl_thick_u(i,j) = bbl_thick_z
1076  else
1077  visc%Kv_bbl_v(i,j) = kv_bbl
1078  visc%bbl_thick_v(i,j) = bbl_thick_z
1079  endif
1080  endif ; enddo ! end of i loop
1081  enddo ; enddo ! end of m & j loops
1082 
1083 ! Offer diagnostics for averaging
1084  if (cs%id_bbl_thick_u > 0) &
1085  call post_data(cs%id_bbl_thick_u, visc%bbl_thick_u, cs%diag)
1086  if (cs%id_kv_bbl_u > 0) &
1087  call post_data(cs%id_kv_bbl_u, visc%kv_bbl_u, cs%diag)
1088  if (cs%id_bbl_u > 0) &
1089  call post_data(cs%id_bbl_u, cs%bbl_u, cs%diag)
1090  if (cs%id_bbl_thick_v > 0) &
1091  call post_data(cs%id_bbl_thick_v, visc%bbl_thick_v, cs%diag)
1092  if (cs%id_kv_bbl_v > 0) &
1093  call post_data(cs%id_kv_bbl_v, visc%kv_bbl_v, cs%diag)
1094  if (cs%id_bbl_v > 0) &
1095  call post_data(cs%id_bbl_v, cs%bbl_v, cs%diag)
1096  if (cs%id_Ray_u > 0) &
1097  call post_data(cs%id_Ray_u, visc%Ray_u, cs%diag)
1098  if (cs%id_Ray_v > 0) &
1099  call post_data(cs%id_Ray_v, visc%Ray_v, cs%diag)
1100 
1101  if (cs%debug) then
1102  if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) &
1103  call uvchksum("Ray [uv]", visc%Ray_u, visc%Ray_v, g%HI, haloshift=0, scale=us%Z_to_m*us%s_to_T)
1104  if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) &
1105  call uvchksum("kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, g%HI, &
1106  haloshift=0, scale=us%Z2_T_to_m2_s, scalar_pair=.true.)
1107  if (associated(visc%bbl_thick_u) .and. associated(visc%bbl_thick_v)) &
1108  call uvchksum("bbl_thick_[uv]", visc%bbl_thick_u, visc%bbl_thick_v, &
1109  g%HI, haloshift=0, scale=us%Z_to_m, scalar_pair=.true.)
1110  endif
1111 

◆ set_viscous_ml()

subroutine, public mom_set_visc::set_viscous_ml ( real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)  u,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in)  v,
real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)  h,
type(thermo_var_ptrs), intent(in)  tv,
type(mech_forcing), intent(in)  forces,
type(vertvisc_type), intent(inout)  visc,
real, intent(in)  dt,
type(ocean_grid_type), intent(inout)  G,
type(verticalgrid_type), intent(in)  GV,
type(unit_scale_type), intent(in)  US,
type(set_visc_cs), pointer  CS,
logical, intent(in), optional  symmetrize 
)

Calculates the thickness of the surface boundary layer for applying an elevated viscosity.

A bulk Richardson criterion or the thickness of the topmost NKML layers (with a bulk mixed layer) are currently used. The thicknesses are given in terms of fractional layers, so that this thickness will move as the thickness of the topmost layers change.

Parameters
[in,out]gThe ocean's grid structure.
[in]gvThe ocean's vertical grid structure.
[in]usA dimensional unit scaling type
[in]uThe zonal velocity [L T-1 ~> m s-1].
[in]vThe meridional velocity [L T-1 ~> m s-1].
[in]hLayer thicknesses [H ~> m or kg m-2].
[in]tvA structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
[in]forcesA structure with the driving mechanical forces
[in,out]viscA structure containing vertical viscosities and related fields.
[in]dtTime increment [T ~> s].
csThe control structure returned by a previous call to set_visc_init.
[in]symmetrizeIf present and true, do extra calculations of those values in visc that would be calculated with symmetric memory.

Definition at line 1208 of file MOM_set_viscosity.F90.

1208  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
1209  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1210  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1211  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1212  intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1].
1213  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1214  intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1].
1215  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1216  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
1217  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
1218  !! thermodynamic fields. Absent fields have
1219  !! NULL ptrs.
1220  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
1221  type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
1222  !! related fields.
1223  real, intent(in) :: dt !< Time increment [T ~> s].
1224  type(set_visc_CS), pointer :: CS !< The control structure returned by a previous
1225  !! call to set_visc_init.
1226  logical, optional, intent(in) :: symmetrize !< If present and true, do extra calculations
1227  !! of those values in visc that would be
1228  !! calculated with symmetric memory.
1229 
1230  ! Local variables
1231  real, dimension(SZIB_(G)) :: &
1232  htot, & ! The total depth of the layers being that are within the
1233  ! surface mixed layer [H ~> m or kg m-2].
1234  thtot, & ! The integrated temperature of layers that are within the
1235  ! surface mixed layer [H degC ~> m degC or kg degC m-2].
1236  shtot, & ! The integrated salt of layers that are within the
1237  ! surface mixed layer [H ppt ~> m ppt or kg ppt m-2].
1238  rhtot, & ! The integrated density of layers that are within the surface mixed layer
1239  ! [H R ~> kg m-2 or kg2 m-5]. Rhtot is only used if no
1240  ! equation of state is used.
1241  uhtot, & ! The depth integrated zonal and meridional velocities within
1242  vhtot, & ! the surface mixed layer [H L T-1 ~> m2 s-1 or kg m-1 s-1].
1243  idecay_len_tke, & ! The inverse of a turbulence decay length scale [H-1 ~> m-1 or m2 kg-1].
1244  dr_dt, & ! Partial derivative of the density at the base of layer nkml
1245  ! (roughly the base of the mixed layer) with temperature [R degC-1 ~> kg m-3 degC-1].
1246  dr_ds, & ! Partial derivative of the density at the base of layer nkml
1247  ! (roughly the base of the mixed layer) with salinity [R ppt-1 ~> kg m-3 ppt-1].
1248  ustar, & ! The surface friction velocity under ice shelves [Z T-1 ~> m s-1].
1249  press, & ! The pressure at which dR_dT and dR_dS are evaluated [R L2 T-2 ~> Pa].
1250  t_eos, & ! The potential temperature at which dR_dT and dR_dS are evaluated [degC]
1251  s_eos ! The salinity at which dR_dT and dR_dS are evaluated [ppt].
1252  real, dimension(SZIB_(G),SZJ_(G)) :: &
1253  mask_u ! A mask that disables any contributions from u points that
1254  ! are land or past open boundary conditions [nondim], 0 or 1.
1255  real, dimension(SZI_(G),SZJB_(G)) :: &
1256  mask_v ! A mask that disables any contributions from v points that
1257  ! are land or past open boundary conditions [nondim], 0 or 1.
1258  real :: h_at_vel(SZIB_(G),SZK_(G))! Layer thickness at velocity points,
1259  ! using an upwind-biased second order accurate estimate based
1260  ! on the previous velocity direction [H ~> m or kg m-2].
1261  integer :: k_massive(SZIB_(G)) ! The k-index of the deepest layer yet found
1262  ! that has more than h_tiny thickness and will be in the
1263  ! viscous mixed layer.
1264  real :: Uh2 ! The squared magnitude of the difference between the velocity
1265  ! integrated through the mixed layer and the velocity of the
1266  ! interior layer layer times the depth of the the mixed layer
1267  ! [H2 Z2 T-2 ~> m4 s-2 or kg2 m-2 s-2].
1268  real :: htot_vel ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
1269  real :: hwtot ! Sum of the thicknesses used to calculate
1270  ! the near-bottom velocity magnitude [H ~> m or kg m-2].
1271  real :: hutot ! Running sum of thicknesses times the
1272  ! velocity magnitudes [H L T-1 ~> m2 s-1 or kg m-1 s-1].
1273  real :: hweight ! The thickness of a layer that is within Hbbl
1274  ! of the bottom [H ~> m or kg m-2].
1275  real :: tbl_thick_Z ! The thickness of the top boundary layer [Z ~> m].
1276 
1277  real :: hlay ! The layer thickness at velocity points [H ~> m or kg m-2].
1278  real :: I_2hlay ! 1 / 2*hlay [H-1 ~> m-1 or m2 kg-1].
1279  real :: T_lay ! The layer temperature at velocity points [degC].
1280  real :: S_lay ! The layer salinity at velocity points [ppt].
1281  real :: Rlay ! The layer potential density at velocity points [R ~> kg m-3].
1282  real :: Rlb ! The potential density of the layer below [R ~> kg m-3].
1283  real :: v_at_u ! The meridonal velocity at a zonal velocity point [L T-1 ~> m s-1].
1284  real :: u_at_v ! The zonal velocity at a meridonal velocity point [L T-1 ~> m s-1].
1285  real :: gHprime ! The mixed-layer internal gravity wave speed squared, based
1286  ! on the mixed layer thickness and density difference across
1287  ! the base of the mixed layer [L2 T-2 ~> m2 s-2].
1288  real :: RiBulk ! The bulk Richardson number below which water is in the
1289  ! viscous mixed layer, including reduction for turbulent
1290  ! decay. Nondimensional.
1291  real :: dt_Rho0 ! The time step divided by the conversion from the layer
1292  ! thickness to layer mass [T H Z-1 R-1 ~> s m3 kg-1 or s].
1293  real :: g_H_Rho0 ! The gravitational acceleration times the conversion from H to m divided
1294  ! by the mean density [L2 T-2 H-1 R-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2].
1295  real :: ustarsq ! 400 times the square of ustar, times
1296  ! Rho0 divided by G_Earth and the conversion
1297  ! from m to thickness units [H R ~> kg m-2 or kg2 m-5].
1298  real :: cdrag_sqrt_Z ! Square root of the drag coefficient, times a unit conversion
1299  ! factor from lateral lengths to vertical depths [Z L-1 ~> 1].
1300  real :: cdrag_sqrt ! Square root of the drag coefficient [nondim].
1301  real :: oldfn ! The integrated energy required to
1302  ! entrain up to the bottom of the layer,
1303  ! divided by G_Earth [H R ~> kg m-2 or kg2 m-5].
1304  real :: Dfn ! The increment in oldfn for entraining
1305  ! the layer [H R ~> kg m-2 or kg2 m-5].
1306  real :: Dh ! The increment in layer thickness from
1307  ! the present layer [H ~> m or kg m-2].
1308  real :: U_bg_sq ! The square of an assumed background velocity, for
1309  ! calculating the mean magnitude near the top for use in
1310  ! the quadratic surface drag [L2 T-2 ~> m2 s-2].
1311  real :: h_tiny ! A very small thickness [H ~> m or kg m-2]. Layers that are less than
1312  ! h_tiny can not be the deepest in the viscous mixed layer.
1313  real :: absf ! The absolute value of f averaged to velocity points [T-1 ~> s-1].
1314  real :: U_star ! The friction velocity at velocity points [Z T-1 ~> m s-1].
1315  real :: h_neglect ! A thickness that is so small it is usually lost
1316  ! in roundoff and can be neglected [H ~> m or kg m-2].
1317  real :: Rho0x400_G ! 400*Rho0/G_Earth, times unit conversion factors
1318  ! [R T2 H Z-2 ~> kg s2 m-4 or kg2 s2 m-7].
1319  ! The 400 is a constant proposed by Killworth and Edwards, 1999.
1320  real :: ustar1 ! ustar [H T-1 ~> m s-1 or kg m-2 s-1]
1321  real :: h2f2 ! (h*2*f)^2 [H2 T-2 ~> m2 s-2 or kg2 m-4 s-2]
1322  logical :: use_EOS, do_any, do_any_shelf, do_i(SZIB_(G))
1323  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, K2, nkmb, nkml, n
1324  type(ocean_OBC_type), pointer :: OBC => null()
1325 
1326  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1327  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1328  nkmb = gv%nk_rho_varies ; nkml = gv%nkml
1329 
1330  if (.not.associated(cs)) call mom_error(fatal,"MOM_set_viscosity(visc_ML): "//&
1331  "Module must be initialized before it is used.")
1332  if (.not.(cs%dynamic_viscous_ML .or. associated(forces%frac_shelf_u) .or. &
1333  associated(forces%frac_shelf_v)) ) return
1334 
1335  if (present(symmetrize)) then ; if (symmetrize) then
1336  jsq = js-1 ; isq = is-1
1337  endif ; endif
1338 
1339  rho0x400_g = 400.0*(gv%Rho0/(us%L_to_Z**2 * gv%g_Earth)) * gv%Z_to_H
1340  u_bg_sq = cs%drag_bg_vel * cs%drag_bg_vel
1341  cdrag_sqrt = sqrt(cs%cdrag)
1342  cdrag_sqrt_z = us%L_to_Z * sqrt(cs%cdrag)
1343 
1344  obc => cs%OBC
1345  use_eos = associated(tv%eqn_of_state)
1346  dt_rho0 = dt / gv%H_to_RZ
1347  h_neglect = gv%H_subroundoff
1348  h_tiny = 2.0*gv%Angstrom_H + h_neglect
1349  g_h_rho0 = (gv%g_Earth*gv%H_to_Z) / (gv%Rho0)
1350 
1351  if (associated(forces%frac_shelf_u) .neqv. associated(forces%frac_shelf_v)) &
1352  call mom_error(fatal, "set_viscous_ML: one of forces%frac_shelf_u and "//&
1353  "forces%frac_shelf_v is associated, but the other is not.")
1354 
1355  if (associated(forces%frac_shelf_u)) then
1356  ! This configuration has ice shelves, and the appropriate variables need to be
1357  ! allocated. If the arrays have already been allocated, these calls do nothing.
1358  call safe_alloc_ptr(visc%tauy_shelf, g%isd, g%ied, g%JsdB, g%JedB)
1359  call safe_alloc_ptr(visc%tbl_thick_shelf_u, g%IsdB, g%IedB, g%jsd, g%jed)
1360  call safe_alloc_ptr(visc%tbl_thick_shelf_v, g%isd, g%ied, g%JsdB, g%JedB)
1361  call safe_alloc_ptr(visc%kv_tbl_shelf_u, g%IsdB, g%IedB, g%jsd, g%jed)
1362  call safe_alloc_ptr(visc%kv_tbl_shelf_v, g%isd, g%ied, g%JsdB, g%JedB)
1363  call safe_alloc_ptr(visc%taux_shelf, g%IsdB, g%IedB, g%jsd, g%jed)
1364  call safe_alloc_ptr(visc%tauy_shelf, g%isd, g%ied, g%JsdB, g%JedB)
1365 
1366  ! With a linear drag law under shelves, the friction velocity is already known.
1367 ! if (CS%linear_drag) ustar(:) = cdrag_sqrt_Z*CS%drag_bg_vel
1368  endif
1369 
1370  !$OMP parallel do default(shared)
1371  do j=js-1,je ; do i=is-1,ie+1
1372  mask_v(i,j) = g%mask2dCv(i,j)
1373  enddo ; enddo
1374  !$OMP parallel do default(shared)
1375  do j=js-1,je+1 ; do i=is-1,ie
1376  mask_u(i,j) = g%mask2dCu(i,j)
1377  enddo ; enddo
1378 
1379  if (associated(obc)) then ; do n=1,obc%number_of_segments
1380  ! Now project bottom depths across cell-corner points in the OBCs. The two
1381  ! projections have to occur in sequence and can not be combined easily.
1382  if (.not. obc%segment(n)%on_pe) cycle
1383  ! Use a one-sided projection of bottom depths at OBC points.
1384  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
1385  if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je)) then
1386  do i = max(is-1,obc%segment(n)%HI%IsdB), min(ie,obc%segment(n)%HI%IedB)
1387  if (obc%segment(n)%direction == obc_direction_n) mask_u(i,j+1) = 0.0
1388  if (obc%segment(n)%direction == obc_direction_s) mask_u(i,j) = 0.0
1389  enddo
1390  elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= je)) then
1391  do j = max(js-1,obc%segment(n)%HI%JsdB), min(je,obc%segment(n)%HI%JedB)
1392  if (obc%segment(n)%direction == obc_direction_e) mask_v(i+1,j) = 0.0
1393  if (obc%segment(n)%direction == obc_direction_w) mask_v(i,j) = 0.0
1394  enddo
1395  endif
1396  enddo ; endif
1397 
1398  !$OMP parallel do default(private) shared(u,v,h,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, &
1399  !$OMP h_neglect,h_tiny,g_H_Rho0,js,je,OBC,Isq,Ieq,nz, &
1400  !$OMP U_bg_sq,mask_v,cdrag_sqrt,cdrag_sqrt_Z,Rho0x400_G,nkml)
1401  do j=js,je ! u-point loop
1402  if (cs%dynamic_viscous_ML) then
1403  do_any = .false.
1404  do i=isq,ieq
1405  htot(i) = 0.0
1406  if (g%mask2dCu(i,j) < 0.5) then
1407  do_i(i) = .false. ; visc%nkml_visc_u(i,j) = nkml
1408  else
1409  do_i(i) = .true. ; do_any = .true.
1410  k_massive(i) = nkml
1411  thtot(i) = 0.0 ; shtot(i) = 0.0 ; rhtot(i) = 0.0
1412  uhtot(i) = dt_rho0 * forces%taux(i,j)
1413  vhtot(i) = 0.25 * dt_rho0 * ((forces%tauy(i,j) + forces%tauy(i+1,j-1)) + &
1414  (forces%tauy(i,j-1) + forces%tauy(i+1,j)))
1415 
1416  if (cs%omega_frac >= 1.0) then ; absf = 2.0*cs%omega ; else
1417  absf = 0.5*(abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i,j-1)))
1418  if (cs%omega_frac > 0.0) &
1419  absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
1420  endif
1421  u_star = max(cs%ustar_min, 0.5 * (forces%ustar(i,j) + forces%ustar(i+1,j)))
1422  idecay_len_tke(i) = ((absf / u_star) * cs%TKE_decay) * gv%H_to_Z
1423  endif
1424  enddo
1425 
1426  if (do_any) then ; do k=1,nz
1427  if (k > nkml) then
1428  do_any = .false.
1429  if (use_eos .and. (k==nkml+1)) then
1430  ! Find dRho/dT and dRho_dS.
1431  do i=isq,ieq
1432  press(i) = (gv%H_to_RZ*gv%g_Earth) * htot(i)
1433  if (associated(tv%p_surf)) press(i) = press(i) + 0.5*(tv%p_surf(i,j)+tv%p_surf(i+1,j))
1434  k2 = max(1,nkml)
1435  i_2hlay = 1.0 / (h(i,j,k2) + h(i+1,j,k2) + h_neglect)
1436  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
1437  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
1438  enddo
1439  call calculate_density_derivs(t_eos, s_eos, press, dr_dt, dr_ds, tv%eqn_of_state, &
1440  (/isq-g%IsdB+1,ieq-g%IsdB+1/) )
1441  endif
1442 
1443  do i=isq,ieq ; if (do_i(i)) then
1444 
1445  hlay = 0.5*(h(i,j,k) + h(i+1,j,k))
1446  if (hlay > h_tiny) then ! Only consider non-vanished layers.
1447  i_2hlay = 1.0 / (h(i,j,k) + h(i+1,j,k))
1448  v_at_u = 0.5 * (h(i,j,k) * (v(i,j,k) + v(i,j-1,k)) + &
1449  h(i+1,j,k) * (v(i+1,j,k) + v(i+1,j-1,k))) * i_2hlay
1450  uh2 = ((uhtot(i) - htot(i)*u(i,j,k))**2 + (vhtot(i) - htot(i)*v_at_u)**2)
1451 
1452  if (use_eos) then
1453  t_lay = (h(i,j,k)*tv%T(i,j,k) + h(i+1,j,k)*tv%T(i+1,j,k)) * i_2hlay
1454  s_lay = (h(i,j,k)*tv%S(i,j,k) + h(i+1,j,k)*tv%S(i+1,j,k)) * i_2hlay
1455  ghprime = g_h_rho0 * (dr_dt(i) * (t_lay*htot(i) - thtot(i)) + &
1456  dr_ds(i) * (s_lay*htot(i) - shtot(i)))
1457  else
1458  ghprime = g_h_rho0 * (gv%Rlay(k)*htot(i) - rhtot(i))
1459  endif
1460 
1461  if (ghprime > 0.0) then
1462  ribulk = cs%bulk_Ri_ML * exp(-htot(i) * idecay_len_tke(i))
1463  if (ribulk * uh2 <= (htot(i)**2) * ghprime) then
1464  visc%nkml_visc_u(i,j) = real(k_massive(I))
1465  do_i(i) = .false.
1466  elseif (ribulk * uh2 <= (htot(i) + hlay)**2 * ghprime) then
1467  visc%nkml_visc_u(i,j) = real(k-1) + &
1468  ( sqrt(RiBulk * Uh2 / gHprime) - htot(I) ) / hlay
1469  do_i(i) = .false.
1470  endif
1471  endif
1472  k_massive(i) = k
1473  endif ! hlay > h_tiny
1474 
1475  if (do_i(i)) do_any = .true.
1476  endif ; enddo
1477 
1478  if (.not.do_any) exit ! All columns are done.
1479  endif
1480 
1481  do i=isq,ieq ; if (do_i(i)) then
1482  htot(i) = htot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k))
1483  uhtot(i) = uhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * u(i,j,k)
1484  vhtot(i) = vhtot(i) + 0.25 * (h(i,j,k) * (v(i,j,k) + v(i,j-1,k)) + &
1485  h(i+1,j,k) * (v(i+1,j,k) + v(i+1,j-1,k)))
1486  if (use_eos) then
1487  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))
1488  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))
1489  else
1490  rhtot(i) = rhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * gv%Rlay(k)
1491  endif
1492  endif ; enddo
1493  enddo ; endif
1494 
1495  if (do_any) then ; do i=isq,ieq ; if (do_i(i)) then
1496  visc%nkml_visc_u(i,j) = k_massive(i)
1497  endif ; enddo ; endif
1498  endif ! dynamic_viscous_ML
1499 
1500  do_any_shelf = .false.
1501  if (associated(forces%frac_shelf_u)) then
1502  do i=isq,ieq
1503  if (forces%frac_shelf_u(i,j)*g%mask2dCu(i,j) == 0.0) then
1504  do_i(i) = .false.
1505  visc%tbl_thick_shelf_u(i,j) = 0.0 ; visc%kv_tbl_shelf_u(i,j) = 0.0
1506  else
1507  do_i(i) = .true. ; do_any_shelf = .true.
1508  endif
1509  enddo
1510  endif
1511 
1512  if (do_any_shelf) then
1513  do k=1,nz ; do i=isq,ieq ; if (do_i(i)) then
1514  if (u(i,j,k) *(h(i+1,j,k) - h(i,j,k)) >= 0) then
1515  h_at_vel(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / &
1516  (h(i,j,k) + h(i+1,j,k) + h_neglect)
1517  else
1518  h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
1519  endif
1520  else
1521  h_at_vel(i,k) = 0.0 ; ustar(i) = 0.0
1522  endif ; enddo ; enddo
1523 
1524  do i=isq,ieq ; if (do_i(i)) then
1525  htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
1526  thtot(i) = 0.0 ; shtot(i) = 0.0
1527  if (use_eos .or. .not.cs%linear_drag) then ; do k=1,nz
1528  if (htot_vel>=cs%Htbl_shelf) exit ! terminate the k loop
1529  hweight = min(cs%Htbl_shelf - htot_vel, h_at_vel(i,k))
1530  if (hweight <= 1.5*gv%Angstrom_H + h_neglect) cycle
1531 
1532  htot_vel = htot_vel + h_at_vel(i,k)
1533  hwtot = hwtot + hweight
1534 
1535  if (.not.cs%linear_drag) then
1536  v_at_u = set_v_at_u(v, h, g, i, j, k, mask_v, obc)
1537  hutot = hutot + hweight * sqrt(u(i,j,k)**2 + &
1538  v_at_u**2 + u_bg_sq)
1539  endif
1540  if (use_eos) then
1541  thtot(i) = thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
1542  shtot(i) = shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
1543  endif
1544  enddo ; endif
1545 
1546  if ((.not.cs%linear_drag) .and. (hwtot > 0.0)) then
1547  ustar(i) = cdrag_sqrt_z * hutot/hwtot
1548  else
1549  ustar(i) = cdrag_sqrt_z * cs%drag_bg_vel
1550  endif
1551 
1552  if (use_eos) then ; if (hwtot > 0.0) then
1553  t_eos(i) = thtot(i)/hwtot ; s_eos(i) = shtot(i)/hwtot
1554  else
1555  t_eos(i) = 0.0 ; s_eos(i) = 0.0
1556  endif ; endif
1557  endif ; enddo ! I-loop
1558 
1559  if (use_eos) then
1560  call calculate_density_derivs(t_eos, s_eos, forces%p_surf(:,j), dr_dt, dr_ds, &
1561  tv%eqn_of_state, (/isq-g%IsdB+1,ieq-g%IsdB+1/) )
1562  endif
1563 
1564  do i=isq,ieq ; if (do_i(i)) then
1565  ! The 400.0 in this expression is the square of a constant proposed
1566  ! by Killworth and Edwards, 1999, in equation (2.20).
1567  ustarsq = rho0x400_g * ustar(i)**2
1568  htot(i) = 0.0
1569  if (use_eos) then
1570  thtot(i) = 0.0 ; shtot(i) = 0.0
1571  do k=1,nz-1
1572  if (h_at_vel(i,k) <= 0.0) cycle
1573  t_lay = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
1574  s_lay = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
1575  oldfn = dr_dt(i)*(t_lay*htot(i) - thtot(i)) + dr_ds(i)*(s_lay*htot(i) - shtot(i))
1576  if (oldfn >= ustarsq) exit
1577 
1578  dfn = (dr_dt(i)*(0.5*(tv%T(i,j,k+1)+tv%T(i+1,j,k+1)) - t_lay) + &
1579  dr_ds(i)*(0.5*(tv%S(i,j,k+1)+tv%S(i+1,j,k+1)) - s_lay)) * &
1580  (h_at_vel(i,k)+htot(i))
1581  if ((oldfn + dfn) <= ustarsq) then
1582  dh = h_at_vel(i,k)
1583  else
1584  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
1585  endif
1586 
1587  htot(i) = htot(i) + dh
1588  thtot(i) = thtot(i) + t_lay*dh ; shtot(i) = shtot(i) + s_lay*dh
1589  enddo
1590  if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0)) then
1591  t_lay = 0.5*(tv%T(i,j,nz) + tv%T(i+1,j,nz))
1592  s_lay = 0.5*(tv%S(i,j,nz) + tv%S(i+1,j,nz))
1593  if (dr_dt(i)*(t_lay*htot(i) - thtot(i)) + &
1594  dr_ds(i)*(s_lay*htot(i) - shtot(i)) < ustarsq) &
1595  htot(i) = htot(i) + h_at_vel(i,nz)
1596  endif ! Examination of layer nz.
1597  else ! Use Rlay as the density variable.
1598  rhtot = 0.0
1599  do k=1,nz-1
1600  rlay = gv%Rlay(k) ; rlb = gv%Rlay(k+1)
1601 
1602  oldfn = rlay*htot(i) - rhtot(i)
1603  if (oldfn >= ustarsq) exit
1604 
1605  dfn = (rlb - rlay)*(h_at_vel(i,k)+htot(i))
1606  if ((oldfn + dfn) <= ustarsq) then
1607  dh = h_at_vel(i,k)
1608  else
1609  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
1610  endif
1611 
1612  htot(i) = htot(i) + dh
1613  rhtot(i) = rhtot(i) + rlay*dh
1614  enddo
1615  if (gv%Rlay(nz)*htot(i) - rhtot(i) < ustarsq) &
1616  htot(i) = htot(i) + h_at_vel(i,nz)
1617  endif ! use_EOS
1618 
1619  !visc%tbl_thick_shelf_u(I,j) = GV%H_to_Z * max(CS%Htbl_shelf_min, &
1620  ! htot(I) / (0.5 + sqrt(0.25 + &
1621  ! (htot(i)*(G%CoriolisBu(I,J-1)+G%CoriolisBu(I,J)))**2 / &
1622  ! (ustar(i)*GV%Z_to_H)**2 )) )
1623  ustar1 = ustar(i)*gv%Z_to_H
1624  h2f2 = (htot(i)*(g%CoriolisBu(i,j-1)+g%CoriolisBu(i,j)) + h_neglect*cs%omega)**2
1625  tbl_thick_z = gv%H_to_Z * max(cs%Htbl_shelf_min, &
1626  ( htot(i)*ustar1 ) / ( 0.5*ustar1 + sqrt((0.5*ustar1)**2 + h2f2 ) ) )
1627  visc%tbl_thick_shelf_u(i,j) = tbl_thick_z
1628  visc%Kv_tbl_shelf_u(i,j) = max(cs%Kv_TBL_min, cdrag_sqrt*ustar(i)*tbl_thick_z)
1629  endif ; enddo ! I-loop
1630  endif ! do_any_shelf
1631 
1632  enddo ! j-loop at u-points
1633 
1634  !$OMP parallel do default(private) shared(u,v,h,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, &
1635  !$OMP h_neglect,h_tiny,g_H_Rho0,is,ie,OBC,Jsq,Jeq,nz, &
1636  !$OMP U_bg_sq,cdrag_sqrt,cdrag_sqrt_Z,Rho0x400_G,nkml,mask_u)
1637  do j=jsq,jeq ! v-point loop
1638  if (cs%dynamic_viscous_ML) then
1639  do_any = .false.
1640  do i=is,ie
1641  htot(i) = 0.0
1642  if (g%mask2dCv(i,j) < 0.5) then
1643  do_i(i) = .false. ; visc%nkml_visc_v(i,j) = nkml
1644  else
1645  do_i(i) = .true. ; do_any = .true.
1646  k_massive(i) = nkml
1647  thtot(i) = 0.0 ; shtot(i) = 0.0 ; rhtot(i) = 0.0
1648  vhtot(i) = dt_rho0 * forces%tauy(i,j)
1649  uhtot(i) = 0.25 * dt_rho0 * ((forces%taux(i,j) + forces%taux(i-1,j+1)) + &
1650  (forces%taux(i-1,j) + forces%taux(i,j+1)))
1651 
1652  if (cs%omega_frac >= 1.0) then ; absf = 2.0*cs%omega ; else
1653  absf = 0.5*(abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))
1654  if (cs%omega_frac > 0.0) &
1655  absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
1656  endif
1657 
1658  u_star = max(cs%ustar_min, 0.5 * (forces%ustar(i,j) + forces%ustar(i,j+1)))
1659  idecay_len_tke(i) = ((absf / u_star) * cs%TKE_decay) * gv%H_to_Z
1660 
1661  endif
1662  enddo
1663 
1664  if (do_any) then ; do k=1,nz
1665  if (k > nkml) then
1666  do_any = .false.
1667  if (use_eos .and. (k==nkml+1)) then
1668  ! Find dRho/dT and dRho_dS.
1669  do i=is,ie
1670  press(i) = (gv%H_to_RZ * gv%g_Earth) * htot(i)
1671  if (associated(tv%p_surf)) press(i) = press(i) + 0.5*(tv%p_surf(i,j)+tv%p_surf(i,j+1))
1672  k2 = max(1,nkml)
1673  i_2hlay = 1.0 / (h(i,j,k2) + h(i,j+1,k2) + h_neglect)
1674  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
1675  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
1676  enddo
1677  call calculate_density_derivs(t_eos, s_eos, press, dr_dt, dr_ds, &
1678  tv%eqn_of_state, (/is-g%IsdB+1,ie-g%IsdB+1/) )
1679  endif
1680 
1681  do i=is,ie ; if (do_i(i)) then
1682 
1683  hlay = 0.5*(h(i,j,k) + h(i,j+1,k))
1684  if (hlay > h_tiny) then ! Only consider non-vanished layers.
1685  i_2hlay = 1.0 / (h(i,j,k) + h(i,j+1,k))
1686  u_at_v = 0.5 * (h(i,j,k) * (u(i-1,j,k) + u(i,j,k)) + &
1687  h(i,j+1,k) * (u(i-1,j+1,k) + u(i,j+1,k))) * i_2hlay
1688  uh2 = ((uhtot(i) - htot(i)*u_at_v)**2 + (vhtot(i) - htot(i)*v(i,j,k))**2)
1689 
1690  if (use_eos) then
1691  t_lay = (h(i,j,k)*tv%T(i,j,k) + h(i,j+1,k)*tv%T(i,j+1,k)) * i_2hlay
1692  s_lay = (h(i,j,k)*tv%S(i,j,k) + h(i,j+1,k)*tv%S(i,j+1,k)) * i_2hlay
1693  ghprime = g_h_rho0 * (dr_dt(i) * (t_lay*htot(i) - thtot(i)) + &
1694  dr_ds(i) * (s_lay*htot(i) - shtot(i)))
1695  else
1696  ghprime = g_h_rho0 * (gv%Rlay(k)*htot(i) - rhtot(i))
1697  endif
1698 
1699  if (ghprime > 0.0) then
1700  ribulk = cs%bulk_Ri_ML * exp(-htot(i) * idecay_len_tke(i))
1701  if (ribulk * uh2 <= htot(i)**2 * ghprime) then
1702  visc%nkml_visc_v(i,j) = real(k_massive(i))
1703  do_i(i) = .false.
1704  elseif (ribulk * uh2 <= (htot(i) + hlay)**2 * ghprime) then
1705  visc%nkml_visc_v(i,j) = real(k-1) + &
1706  ( sqrt(RiBulk * Uh2 / gHprime) - htot(i) ) / hlay
1707  do_i(i) = .false.
1708  endif
1709  endif
1710  k_massive(i) = k
1711  endif ! hlay > h_tiny
1712 
1713  if (do_i(i)) do_any = .true.
1714  endif ; enddo
1715 
1716  if (.not.do_any) exit ! All columns are done.
1717  endif
1718 
1719  do i=is,ie ; if (do_i(i)) then
1720  htot(i) = htot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k))
1721  vhtot(i) = vhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * v(i,j,k)
1722  uhtot(i) = uhtot(i) + 0.25 * (h(i,j,k) * (u(i-1,j,k) + u(i,j,k)) + &
1723  h(i,j+1,k) * (u(i-1,j+1,k) + u(i,j+1,k)))
1724  if (use_eos) then
1725  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))
1726  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))
1727  else
1728  rhtot(i) = rhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * gv%Rlay(k)
1729  endif
1730  endif ; enddo
1731  enddo ; endif
1732 
1733  if (do_any) then ; do i=is,ie ; if (do_i(i)) then
1734  visc%nkml_visc_v(i,j) = k_massive(i)
1735  endif ; enddo ; endif
1736  endif ! dynamic_viscous_ML
1737 
1738  do_any_shelf = .false.
1739  if (associated(forces%frac_shelf_v)) then
1740  do i=is,ie
1741  if (forces%frac_shelf_v(i,j)*g%mask2dCv(i,j) == 0.0) then
1742  do_i(i) = .false.
1743  visc%tbl_thick_shelf_v(i,j) = 0.0 ; visc%kv_tbl_shelf_v(i,j) = 0.0
1744  else
1745  do_i(i) = .true. ; do_any_shelf = .true.
1746  endif
1747  enddo
1748  endif
1749 
1750  if (do_any_shelf) then
1751  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
1752  if (v(i,j,k) * (h(i,j+1,k) - h(i,j,k)) >= 0) then
1753  h_at_vel(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / &
1754  (h(i,j,k) + h(i,j+1,k) + h_neglect)
1755  else
1756  h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
1757  endif
1758  else
1759  h_at_vel(i,k) = 0.0 ; ustar(i) = 0.0
1760  endif ; enddo ; enddo
1761 
1762  do i=is,ie ; if (do_i(i)) then
1763  htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
1764  thtot(i) = 0.0 ; shtot(i) = 0.0
1765  if (use_eos .or. .not.cs%linear_drag) then ; do k=1,nz
1766  if (htot_vel>=cs%Htbl_shelf) exit ! terminate the k loop
1767  hweight = min(cs%Htbl_shelf - htot_vel, h_at_vel(i,k))
1768  if (hweight <= 1.5*gv%Angstrom_H + h_neglect) cycle
1769 
1770  htot_vel = htot_vel + h_at_vel(i,k)
1771  hwtot = hwtot + hweight
1772 
1773  if (.not.cs%linear_drag) then
1774  u_at_v = set_u_at_v(u, h, g, i, j, k, mask_u, obc)
1775  hutot = hutot + hweight * sqrt(v(i,j,k)**2 + &
1776  u_at_v**2 + u_bg_sq)
1777  endif
1778  if (use_eos) then
1779  thtot(i) = thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
1780  shtot(i) = shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
1781  endif
1782  enddo ; endif
1783 
1784  if (.not.cs%linear_drag) then ; if (hwtot > 0.0) then
1785  ustar(i) = cdrag_sqrt_z * hutot/hwtot
1786  else
1787  ustar(i) = cdrag_sqrt_z * cs%drag_bg_vel
1788  endif ; endif
1789 
1790  if (use_eos) then ; if (hwtot > 0.0) then
1791  t_eos(i) = thtot(i)/hwtot ; s_eos(i) = shtot(i)/hwtot
1792  else
1793  t_eos(i) = 0.0 ; s_eos(i) = 0.0
1794  endif ; endif
1795  endif ; enddo ! I-loop
1796 
1797  if (use_eos) then
1798  call calculate_density_derivs(t_eos, s_eos, forces%p_surf(:,j), dr_dt, dr_ds, &
1799  tv%eqn_of_state, (/is-g%IsdB+1,ie-g%IsdB+1/) )
1800  endif
1801 
1802  do i=is,ie ; if (do_i(i)) then
1803  ! The 400.0 in this expression is the square of a constant proposed
1804  ! by Killworth and Edwards, 1999, in equation (2.20).
1805  ustarsq = rho0x400_g * ustar(i)**2
1806  htot(i) = 0.0
1807  if (use_eos) then
1808  thtot(i) = 0.0 ; shtot(i) = 0.0
1809  do k=1,nz-1
1810  if (h_at_vel(i,k) <= 0.0) cycle
1811  t_lay = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
1812  s_lay = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
1813  oldfn = dr_dt(i)*(t_lay*htot(i) - thtot(i)) + dr_ds(i)*(s_lay*htot(i) - shtot(i))
1814  if (oldfn >= ustarsq) exit
1815 
1816  dfn = (dr_dt(i)*(0.5*(tv%T(i,j,k+1)+tv%T(i,j+1,k+1)) - t_lay) + &
1817  dr_ds(i)*(0.5*(tv%S(i,j,k+1)+tv%S(i,j+1,k+1)) - s_lay)) * &
1818  (h_at_vel(i,k)+htot(i))
1819  if ((oldfn + dfn) <= ustarsq) then
1820  dh = h_at_vel(i,k)
1821  else
1822  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
1823  endif
1824 
1825  htot(i) = htot(i) + dh
1826  thtot(i) = thtot(i) + t_lay*dh ; shtot(i) = shtot(i) + s_lay*dh
1827  enddo
1828  if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0)) then
1829  t_lay = 0.5*(tv%T(i,j,nz) + tv%T(i,j+1,nz))
1830  s_lay = 0.5*(tv%S(i,j,nz) + tv%S(i,j+1,nz))
1831  if (dr_dt(i)*(t_lay*htot(i) - thtot(i)) + &
1832  dr_ds(i)*(s_lay*htot(i) - shtot(i)) < ustarsq) &
1833  htot(i) = htot(i) + h_at_vel(i,nz)
1834  endif ! Examination of layer nz.
1835  else ! Use Rlay as the density variable.
1836  rhtot = 0.0
1837  do k=1,nz-1
1838  rlay = gv%Rlay(k) ; rlb = gv%Rlay(k+1)
1839 
1840  oldfn = rlay*htot(i) - rhtot(i)
1841  if (oldfn >= ustarsq) exit
1842 
1843  dfn = (rlb - rlay)*(h_at_vel(i,k)+htot(i))
1844  if ((oldfn + dfn) <= ustarsq) then
1845  dh = h_at_vel(i,k)
1846  else
1847  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn) / (dfn))
1848  endif
1849 
1850  htot(i) = htot(i) + dh
1851  rhtot = rhtot + rlay*dh
1852  enddo
1853  if (gv%Rlay(nz)*htot(i) - rhtot(i) < ustarsq) &
1854  htot(i) = htot(i) + h_at_vel(i,nz)
1855  endif ! use_EOS
1856 
1857  !visc%tbl_thick_shelf_v(i,J) = GV%H_to_Z * max(CS%Htbl_shelf_min, &
1858  ! htot(i) / (0.5 + sqrt(0.25 + &
1859  ! (htot(i)*(G%CoriolisBu(I-1,J)+G%CoriolisBu(I,J)))**2 / &
1860  ! (ustar(i)*GV%Z_to_H)**2 )) )
1861  ustar1 = ustar(i)*gv%Z_to_H
1862  h2f2 = (htot(i)*(g%CoriolisBu(i-1,j)+g%CoriolisBu(i,j)) + h_neglect*cs%omega)**2
1863  tbl_thick_z = gv%H_to_Z * max(cs%Htbl_shelf_min, &
1864  ( htot(i)*ustar1 ) / ( 0.5*ustar1 + sqrt((0.5*ustar1)**2 + h2f2 ) ) )
1865  visc%tbl_thick_shelf_v(i,j) = tbl_thick_z
1866  visc%Kv_tbl_shelf_v(i,j) = max(cs%Kv_TBL_min, cdrag_sqrt*ustar(i)*tbl_thick_z)
1867 
1868  endif ; enddo ! i-loop
1869  endif ! do_any_shelf
1870 
1871  enddo ! J-loop at v-points
1872 
1873  if (cs%debug) then
1874  if (associated(visc%nkml_visc_u) .and. associated(visc%nkml_visc_v)) &
1875  call uvchksum("nkml_visc_[uv]", visc%nkml_visc_u, visc%nkml_visc_v, &
1876  g%HI, haloshift=0, scalar_pair=.true.)
1877  endif
1878  if (cs%id_nkml_visc_u > 0) &
1879  call post_data(cs%id_nkml_visc_u, visc%nkml_visc_u, cs%diag)
1880  if (cs%id_nkml_visc_v > 0) &
1881  call post_data(cs%id_nkml_visc_v, visc%nkml_visc_v, cs%diag)
1882