MOM6
Data Types
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. 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. 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(szib_(g),szj_(g),szk_(g)), intent(in)  u,
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(szib_(g),szj_(g)), 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 [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]mask2dcuA multiplicative mask of the u-points
obcA pointer to an open boundary condition structure
Returns
The return value of u at v points [m s-1].

Definition at line 957 of file MOM_set_viscosity.F90.

957  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
958  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
959  intent(in) :: u !< The zonal velocity [m s-1]
960  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
961  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
962  integer, intent(in) :: i !< The i-index of the u-location to work on.
963  integer, intent(in) :: j !< The j-index of the u-location to work on.
964  integer, intent(in) :: k !< The k-index of the u-location to work on.
965  real, dimension(SZIB_(G),SZJ_(G)), &
966  intent(in) :: mask2dCu !< A multiplicative mask of the u-points
967  type(ocean_OBC_type), pointer :: OBC !< A pointer to an open boundary condition structure
968  real :: set_u_at_v !< The return value of u at v points [m s-1].
969 
970  ! This subroutine finds a thickness-weighted value of u at the v-points.
971  real :: hwt(-1:0,0:1) ! Masked weights used to average u onto v [H ~> m or kg m-2].
972  real :: hwt_tot ! The sum of the masked thicknesses [H ~> m or kg m-2].
973  integer :: i0, j0, i1, j1
974 
975  do j0 = 0,1 ; do i0 = -1,0 ; i1 = i+i0 ; j1 = j+j0
976  hwt(i0,j0) = (h(i1,j1,k) + h(i1+1,j1,k)) * mask2dcu(i1,j1)
977  enddo ; enddo
978 
979  if (associated(obc)) then ; if (obc%number_of_segments > 0) then
980  do j0 = 0,1 ; do i0 = -1,0 ; if ((obc%segnum_u(i+i0,j+j0) /= obc_none)) then
981  i1 = i+i0 ; j1 = j+j0
982  if (obc%segment(obc%segnum_u(i1,j1))%direction == obc_direction_e) then
983  hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dcu(i1,j1)
984  elseif (obc%segment(obc%segnum_u(i1,j1))%direction == obc_direction_w) then
985  hwt(i0,j0) = 2.0 * h(i1+1,j1,k) * mask2dcu(i1,j1)
986  endif
987  endif ; enddo ; enddo
988  endif ; endif
989 
990  hwt_tot = (hwt(-1,0) + hwt(0,1)) + (hwt(0,0) + hwt(-1,1))
991  set_u_at_v = 0.0
992  if (hwt_tot > 0.0) set_u_at_v = &
993  ((hwt(0,0) * u(i,j,k) + hwt(-1,1) * u(i-1,j+1,k)) + &
994  (hwt(-1,0) * u(i-1,j,k) + hwt(0,1) * u(i,j+1,k))) / hwt_tot
995 

◆ set_v_at_u()

real function mom_set_visc::set_v_at_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(ocean_grid_type), intent(in)  G,
integer, intent(in)  i,
integer, intent(in)  j,
integer, intent(in)  k,
real, dimension( g %isd: g %ied, g %jsdb: g %jedb), 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 [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 retur value of v at u points [m s-1].

Definition at line 914 of file MOM_set_viscosity.F90.

914  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
915  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
916  intent(in) :: v !< The meridional velocity [m s-1]
917  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
918  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
919  integer, intent(in) :: i !< The i-index of the u-location to work on.
920  integer, intent(in) :: j !< The j-index of the u-location to work on.
921  integer, intent(in) :: k !< The k-index of the u-location to work on.
922  real, dimension(SZI_(G),SZJB_(G)),&
923  intent(in) :: mask2dCv !< A multiplicative mask of the v-points
924  type(ocean_OBC_type), pointer :: OBC !< A pointer to an open boundary condition structure
925  real :: set_v_at_u !< The retur value of v at u points [m s-1].
926 
927  ! This subroutine finds a thickness-weighted value of v at the u-points.
928  real :: hwt(0:1,-1:0) ! Masked weights used to average u onto v [H ~> m or kg m-2].
929  real :: hwt_tot ! The sum of the masked thicknesses [H ~> m or kg m-2].
930  integer :: i0, j0, i1, j1
931 
932  do j0 = -1,0 ; do i0 = 0,1 ; i1 = i+i0 ; j1 = j+j0
933  hwt(i0,j0) = (h(i1,j1,k) + h(i1,j1+1,k)) * mask2dcv(i1,j1)
934  enddo ; enddo
935 
936  if (associated(obc)) then ; if (obc%number_of_segments > 0) then
937  do j0 = -1,0 ; do i0 = 0,1 ; if ((obc%segnum_v(i+i0,j+j0) /= obc_none)) then
938  i1 = i+i0 ; j1 = j+j0
939  if (obc%segment(obc%segnum_v(i1,j1))%direction == obc_direction_n) then
940  hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dcv(i1,j1)
941  elseif (obc%segment(obc%segnum_v(i1,j1))%direction == obc_direction_s) then
942  hwt(i0,j0) = 2.0 * h(i1,j1+1,k) * mask2dcv(i1,j1)
943  endif
944  endif ; enddo ; enddo
945  endif ; endif
946 
947  hwt_tot = (hwt(0,-1) + hwt(1,0)) + (hwt(1,-1) + hwt(0,0))
948  set_v_at_u = 0.0
949  if (hwt_tot > 0.0) set_v_at_u = &
950  ((hwt(0,0) * v(i,j,k) + hwt(1,-1) * v(i+1,j-1,k)) + &
951  (hwt(1,0) * v(i+1,j,k) + hwt(0,-1) * v(i,j-1,k))) / hwt_tot
952 

◆ 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 vertvisc_init.

Definition at line 2083 of file MOM_set_viscosity.F90.

2083  type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
2084  !! related fields. Elements are deallocated here.
2085  type(set_visc_CS), pointer :: CS !< The control structure returned by a previous
2086  !! call to vertvisc_init.
2087  if (cs%bottomdraglaw) then
2088  deallocate(visc%bbl_thick_u) ; deallocate(visc%bbl_thick_v)
2089  deallocate(visc%kv_bbl_u) ; deallocate(visc%kv_bbl_v)
2090  endif
2091  if (cs%Channel_drag) then
2092  deallocate(visc%Ray_u) ; deallocate(visc%Ray_v)
2093  endif
2094  if (cs%dynamic_viscous_ML) then
2095  deallocate(visc%nkml_visc_u) ; deallocate(visc%nkml_visc_v)
2096  endif
2097  if (associated(visc%Kd_shear)) deallocate(visc%Kd_shear)
2098  if (associated(visc%Kv_slow)) deallocate(visc%Kv_slow)
2099  if (associated(visc%TKE_turb)) deallocate(visc%TKE_turb)
2100  if (associated(visc%Kv_shear)) deallocate(visc%Kv_shear)
2101  if (associated(visc%Kv_shear_Bu)) deallocate(visc%Kv_shear_Bu)
2102  if (associated(visc%ustar_bbl)) deallocate(visc%ustar_bbl)
2103  if (associated(visc%TKE_bbl)) deallocate(visc%TKE_bbl)
2104  if (associated(visc%taux_shelf)) deallocate(visc%taux_shelf)
2105  if (associated(visc%tauy_shelf)) deallocate(visc%tauy_shelf)
2106  if (associated(visc%tbl_thick_shelf_u)) deallocate(visc%tbl_thick_shelf_u)
2107  if (associated(visc%tbl_thick_shelf_v)) deallocate(visc%tbl_thick_shelf_v)
2108  if (associated(visc%kv_tbl_shelf_u)) deallocate(visc%kv_tbl_shelf_u)
2109  if (associated(visc%kv_tbl_shelf_v)) deallocate(visc%kv_tbl_shelf_v)
2110 
2111  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 1765 of file MOM_set_viscosity.F90.

1765  type(time_type), target, intent(in) :: Time !< The current model time.
1766  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
1767  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1768  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1769  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1770  !! parameters.
1771  type(diag_ctrl), target, intent(inout) :: diag !< A structure that is used to regulate diagnostic
1772  !! output.
1773  type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
1774  !! related fields. Allocated here.
1775  type(set_visc_CS), pointer :: CS !< A pointer that is set to point to the control
1776  !! structure for this module
1777  type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure.
1778  type(ocean_OBC_type), pointer :: OBC !< A pointer to an open boundary condition structure
1779 
1780  ! Local variables
1781  real :: Csmag_chan_dflt, smag_const1, TKE_decay_dflt, bulk_Ri_ML_dflt
1782  real :: Kv_background
1783  real :: omega_frac_dflt
1784  real :: Z_rescale ! A rescaling factor for heights from the representation in
1785  ! a restart file to the internal representation in this run.
1786  real :: I_T_rescale ! A rescaling factor for time from the internal representation in this run
1787  ! to the representation in a restart file.
1788  real :: Z2_T_rescale ! A rescaling factor for vertical diffusivities and viscosities from the
1789  ! representation in a restart file to the internal representation in this run.
1790  integer :: i, j, k, is, ie, js, je, n
1791  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
1792  logical :: default_2018_answers
1793  logical :: use_kappa_shear, adiabatic, use_omega
1794  logical :: use_CVMix_ddiff, differential_diffusion, use_KPP
1795  type(OBC_segment_type), pointer :: segment => null() ! pointer to OBC segment type
1796  ! This include declares and sets the variable "version".
1797 # include "version_variable.h"
1798  character(len=40) :: mdl = "MOM_set_visc" ! This module's name.
1799 
1800  if (associated(cs)) then
1801  call mom_error(warning, "set_visc_init called with an associated "// &
1802  "control structure.")
1803  return
1804  endif
1805  allocate(cs)
1806 
1807  cs%OBC => obc
1808 
1809  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
1810  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nz = gv%ke
1811  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
1812 
1813  cs%diag => diag
1814 
1815  ! Set default, read and log parameters
1816  call log_version(param_file, mdl, version, "")
1817  cs%RiNo_mix = .false. ; use_cvmix_ddiff = .false.
1818  differential_diffusion = .false.
1819  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
1820  "This sets the default value for the various _2018_ANSWERS parameters.", &
1821  default=.true.)
1822  call get_param(param_file, mdl, "SET_VISC_2018_ANSWERS", cs%answers_2018, &
1823  "If true, use the order of arithmetic and expressions that recover the "//&
1824  "answers from the end of 2018. Otherwise, use updated and more robust "//&
1825  "forms of the same expressions.", default=default_2018_answers)
1826  call get_param(param_file, mdl, "BOTTOMDRAGLAW", cs%bottomdraglaw, &
1827  "If true, the bottom stress is calculated with a drag "//&
1828  "law of the form c_drag*|u|*u. The velocity magnitude "//&
1829  "may be an assumed value or it may be based on the "//&
1830  "actual velocity in the bottommost HBBL, depending on "//&
1831  "LINEAR_DRAG.", default=.true.)
1832  call get_param(param_file, mdl, "CHANNEL_DRAG", cs%Channel_drag, &
1833  "If true, the bottom drag is exerted directly on each "//&
1834  "layer proportional to the fraction of the bottom it "//&
1835  "overlies.", default=.false.)
1836  call get_param(param_file, mdl, "LINEAR_DRAG", cs%linear_drag, &
1837  "If LINEAR_DRAG and BOTTOMDRAGLAW are defined the drag "//&
1838  "law is cdrag*DRAG_BG_VEL*u.", default=.false.)
1839  call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
1840  do_not_log=.true.)
1841  if (adiabatic) then
1842  call log_param(param_file, mdl, "ADIABATIC",adiabatic, &
1843  "There are no diapycnal mass fluxes if ADIABATIC is "//&
1844  "true. This assumes that KD = KDML = 0.0 and that "//&
1845  "there is no buoyancy forcing, but makes the model "//&
1846  "faster by eliminating subroutine calls.", default=.false.)
1847  endif
1848 
1849  if (.not.adiabatic) then
1850  cs%RiNo_mix = kappa_shear_is_used(param_file)
1851  call get_param(param_file, mdl, "DOUBLE_DIFFUSION", differential_diffusion, &
1852  "If true, increase diffusivites for temperature or salt "//&
1853  "based on double-diffusive parameterization from MOM4/KPP.", &
1854  default=.false.)
1855  use_cvmix_ddiff = cvmix_ddiff_is_used(param_file)
1856  endif
1857 
1858  call get_param(param_file, mdl, "PRANDTL_TURB", visc%Prandtl_turb, &
1859  "The turbulent Prandtl number applied to shear "//&
1860  "instability.", units="nondim", default=1.0)
1861  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
1862 
1863  call get_param(param_file, mdl, "DYNAMIC_VISCOUS_ML", cs%dynamic_viscous_ML, &
1864  "If true, use a bulk Richardson number criterion to "//&
1865  "determine the mixed layer thickness for viscosity.", &
1866  default=.false.)
1867  if (cs%dynamic_viscous_ML) then
1868  call get_param(param_file, mdl, "BULK_RI_ML", bulk_ri_ml_dflt, default=0.0)
1869  call get_param(param_file, mdl, "BULK_RI_ML_VISC", cs%bulk_Ri_ML, &
1870  "The efficiency with which mean kinetic energy released "//&
1871  "by mechanically forced entrainment of the mixed layer "//&
1872  "is converted to turbulent kinetic energy. By default, "//&
1873  "BULK_RI_ML_VISC = BULK_RI_ML or 0.", units="nondim", &
1874  default=bulk_ri_ml_dflt)
1875  call get_param(param_file, mdl, "TKE_DECAY", tke_decay_dflt, default=0.0)
1876  call get_param(param_file, mdl, "TKE_DECAY_VISC", cs%TKE_decay, &
1877  "TKE_DECAY_VISC relates the vertical rate of decay of "//&
1878  "the TKE available for mechanical entrainment to the "//&
1879  "natural Ekman depth for use in calculating the dynamic "//&
1880  "mixed layer viscosity. By default, "//&
1881  "TKE_DECAY_VISC = TKE_DECAY or 0.", units="nondim", &
1882  default=tke_decay_dflt)
1883  call get_param(param_file, mdl, "ML_USE_OMEGA", use_omega, &
1884  "If true, use the absolute rotation rate instead of the "//&
1885  "vertical component of rotation when setting the decay "//&
1886  "scale for turbulence.", default=.false., do_not_log=.true.)
1887  omega_frac_dflt = 0.0
1888  if (use_omega) then
1889  call mom_error(warning, "ML_USE_OMEGA is depricated; use ML_OMEGA_FRAC=1.0 instead.")
1890  omega_frac_dflt = 1.0
1891  endif
1892  call get_param(param_file, mdl, "ML_OMEGA_FRAC", cs%omega_frac, &
1893  "When setting the decay scale for turbulence, use this "//&
1894  "fraction of the absolute rotation rate blended with the "//&
1895  "local value of f, as sqrt((1-of)*f^2 + of*4*omega^2).", &
1896  units="nondim", default=omega_frac_dflt)
1897  call get_param(param_file, mdl, "OMEGA", cs%omega, &
1898  "The rotation rate of the earth.", units="s-1", &
1899  default=7.2921e-5, scale=us%T_to_s)
1900  ! This give a minimum decay scale that is typically much less than Angstrom.
1901  cs%ustar_min = 2e-4*cs%omega*(gv%Angstrom_Z + gv%H_to_Z*gv%H_subroundoff)
1902  else
1903  call get_param(param_file, mdl, "OMEGA", cs%omega, &
1904  "The rotation rate of the earth.", units="s-1", &
1905  default=7.2921e-5, scale=us%T_to_s)
1906  endif
1907 
1908  call get_param(param_file, mdl, "HBBL", cs%Hbbl, &
1909  "The thickness of a bottom boundary layer with a "//&
1910  "viscosity of KVBBL if BOTTOMDRAGLAW is not defined, or "//&
1911  "the thickness over which near-bottom velocities are "//&
1912  "averaged for the drag law if BOTTOMDRAGLAW is defined "//&
1913  "but LINEAR_DRAG is not.", units="m", fail_if_missing=.true.) ! Rescaled later
1914  if (cs%bottomdraglaw) then
1915  call get_param(param_file, mdl, "CDRAG", cs%cdrag, &
1916  "CDRAG is the drag coefficient relating the magnitude of "//&
1917  "the velocity field to the bottom stress. CDRAG is only "//&
1918  "used if BOTTOMDRAGLAW is defined.", units="nondim", &
1919  default=0.003)
1920  call get_param(param_file, mdl, "DRAG_BG_VEL", cs%drag_bg_vel, &
1921  "DRAG_BG_VEL is either the assumed bottom velocity (with "//&
1922  "LINEAR_DRAG) or an unresolved velocity that is "//&
1923  "combined with the resolved velocity to estimate the "//&
1924  "velocity magnitude. DRAG_BG_VEL is only used when "//&
1925  "BOTTOMDRAGLAW is defined.", units="m s-1", default=0.0)
1926  call get_param(param_file, mdl, "BBL_USE_EOS", cs%BBL_use_EOS, &
1927  "If true, use the equation of state in determining the "//&
1928  "properties of the bottom boundary layer. Otherwise use "//&
1929  "the layer target potential densities.", default=.false.)
1930  endif
1931  call get_param(param_file, mdl, "BBL_THICK_MIN", cs%BBL_thick_min, &
1932  "The minimum bottom boundary layer thickness that can be "//&
1933  "used with BOTTOMDRAGLAW. This might be "//&
1934  "Kv/(cdrag*drag_bg_vel) to give Kv as the minimum "//&
1935  "near-bottom viscosity.", units="m", default=0.0) ! Rescaled later
1936  call get_param(param_file, mdl, "HTBL_SHELF_MIN", cs%Htbl_shelf_min, &
1937  "The minimum top boundary layer thickness that can be "//&
1938  "used with BOTTOMDRAGLAW. This might be "//&
1939  "Kv/(cdrag*drag_bg_vel) to give Kv as the minimum "//&
1940  "near-top viscosity.", units="m", default=cs%BBL_thick_min, scale=gv%m_to_H)
1941  call get_param(param_file, mdl, "HTBL_SHELF", cs%Htbl_shelf, &
1942  "The thickness over which near-surface velocities are "//&
1943  "averaged for the drag law under an ice shelf. By "//&
1944  "default this is the same as HBBL", units="m", default=cs%Hbbl, scale=gv%m_to_H)
1945  ! These unit conversions are out outside the get_param calls because the are also defaults.
1946  cs%Hbbl = cs%Hbbl * gv%m_to_H ! Rescale
1947  cs%BBL_thick_min = cs%BBL_thick_min * gv%m_to_H ! Rescale
1948 
1949  call get_param(param_file, mdl, "KV", kv_background, &
1950  "The background kinematic viscosity in the interior. "//&
1951  "The molecular value, ~1e-6 m2 s-1, may be used.", &
1952  units="m2 s-1", fail_if_missing=.true.)
1953 
1954  call get_param(param_file, mdl, "ADD_KV_SLOW", visc%add_Kv_slow, &
1955  "If true, the background vertical viscosity in the interior "//&
1956  "(i.e., tidal + background + shear + convection) is added "//&
1957  "when computing the coupling coefficient. The purpose of this "//&
1958  "flag is to be able to recover previous answers and it will likely "//&
1959  "be removed in the future since this option should always be true.", &
1960  default=.false.)
1961 
1962  call get_param(param_file, mdl, "USE_KPP", use_kpp, &
1963  "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
1964  "to calculate diffusivities and non-local transport in the OBL.", &
1965  do_not_log=.true., default=.false.)
1966 
1967  if (use_kpp .and. visc%add_Kv_slow) call mom_error(fatal,"set_visc_init: "//&
1968  "When USE_KPP=True, ADD_KV_SLOW must be false. Otherwise vertical "//&
1969  "viscosity due to slow processes will be double counted. Please set "//&
1970  "ADD_KV_SLOW=False.")
1971 
1972  call get_param(param_file, mdl, "KV_BBL_MIN", cs%KV_BBL_min, &
1973  "The minimum viscosities in the bottom boundary layer.", &
1974  units="m2 s-1", default=kv_background, scale=us%m2_s_to_Z2_T)
1975  call get_param(param_file, mdl, "KV_TBL_MIN", cs%KV_TBL_min, &
1976  "The minimum viscosities in the top boundary layer.", &
1977  units="m2 s-1", default=kv_background, scale=us%m2_s_to_Z2_T)
1978 
1979  if (cs%Channel_drag) then
1980  call get_param(param_file, mdl, "SMAG_LAP_CONST", smag_const1, default=-1.0)
1981 
1982  csmag_chan_dflt = 0.15
1983  if (smag_const1 >= 0.0) csmag_chan_dflt = smag_const1
1984 
1985  call get_param(param_file, mdl, "SMAG_CONST_CHANNEL", cs%c_Smag, &
1986  "The nondimensional Laplacian Smagorinsky constant used "//&
1987  "in calculating the channel drag if it is enabled. The "//&
1988  "default is to use the same value as SMAG_LAP_CONST if "//&
1989  "it is defined, or 0.15 if it is not. The value used is "//&
1990  "also 0.15 if the specified value is negative.", &
1991  units="nondim", default=csmag_chan_dflt)
1992  if (cs%c_Smag < 0.0) cs%c_Smag = 0.15
1993  endif
1994 
1995  if (cs%RiNo_mix .and. kappa_shear_at_vertex(param_file)) then
1996  ! This is necessary for reproduciblity across restarts in non-symmetric mode.
1997  call pass_var(visc%Kv_shear_Bu, g%Domain, position=corner, complete=.true.)
1998  endif
1999 
2000  if (cs%bottomdraglaw) then
2001  allocate(visc%bbl_thick_u(isdb:iedb,jsd:jed)) ; visc%bbl_thick_u = 0.0
2002  allocate(visc%kv_bbl_u(isdb:iedb,jsd:jed)) ; visc%kv_bbl_u = 0.0
2003  allocate(visc%bbl_thick_v(isd:ied,jsdb:jedb)) ; visc%bbl_thick_v = 0.0
2004  allocate(visc%kv_bbl_v(isd:ied,jsdb:jedb)) ; visc%kv_bbl_v = 0.0
2005  allocate(visc%ustar_bbl(isd:ied,jsd:jed)) ; visc%ustar_bbl = 0.0
2006  allocate(visc%TKE_bbl(isd:ied,jsd:jed)) ; visc%TKE_bbl = 0.0
2007 
2008  cs%id_bbl_thick_u = register_diag_field('ocean_model', 'bbl_thick_u', &
2009  diag%axesCu1, time, 'BBL thickness at u points', 'm', conversion=us%Z_to_m)
2010  cs%id_kv_bbl_u = register_diag_field('ocean_model', 'kv_bbl_u', diag%axesCu1, &
2011  time, 'BBL viscosity at u points', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
2012  cs%id_bbl_thick_v = register_diag_field('ocean_model', 'bbl_thick_v', &
2013  diag%axesCv1, time, 'BBL thickness at v points', 'm', conversion=us%Z_to_m)
2014  cs%id_kv_bbl_v = register_diag_field('ocean_model', 'kv_bbl_v', diag%axesCv1, &
2015  time, 'BBL viscosity at v points', 'm2 s-1', conversion=us%Z2_T_to_m2_s)
2016  endif
2017  if (cs%Channel_drag) then
2018  allocate(visc%Ray_u(isdb:iedb,jsd:jed,nz)) ; visc%Ray_u = 0.0
2019  allocate(visc%Ray_v(isd:ied,jsdb:jedb,nz)) ; visc%Ray_v = 0.0
2020  cs%id_Ray_u = register_diag_field('ocean_model', 'Rayleigh_u', diag%axesCuL, &
2021  time, 'Rayleigh drag velocity at u points', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
2022  cs%id_Ray_v = register_diag_field('ocean_model', 'Rayleigh_v', diag%axesCvL, &
2023  time, 'Rayleigh drag velocity at v points', 'm s-1', conversion=us%Z_to_m*us%s_to_T)
2024  endif
2025 
2026  if (use_cvmix_ddiff .or. differential_diffusion) then
2027  allocate(visc%Kd_extra_T(isd:ied,jsd:jed,nz+1)) ; visc%Kd_extra_T = 0.0
2028  allocate(visc%Kd_extra_S(isd:ied,jsd:jed,nz+1)) ; visc%Kd_extra_S = 0.0
2029  endif
2030 
2031  if (cs%dynamic_viscous_ML) then
2032  allocate(visc%nkml_visc_u(isdb:iedb,jsd:jed)) ; visc%nkml_visc_u = 0.0
2033  allocate(visc%nkml_visc_v(isd:ied,jsdb:jedb)) ; visc%nkml_visc_v = 0.0
2034  cs%id_nkml_visc_u = register_diag_field('ocean_model', 'nkml_visc_u', &
2035  diag%axesCu1, time, 'Number of layers in viscous mixed layer at u points', 'm')
2036  cs%id_nkml_visc_v = register_diag_field('ocean_model', 'nkml_visc_v', &
2037  diag%axesCv1, time, 'Number of layers in viscous mixed layer at v points', 'm')
2038  endif
2039 
2040  call register_restart_field_as_obsolete('Kd_turb','Kd_shear', restart_cs)
2041  call register_restart_field_as_obsolete('Kv_turb','Kv_shear', restart_cs)
2042 
2043  ! Account for possible changes in dimensional scaling for variables that have been
2044  ! read from a restart file.
2045  z_rescale = 1.0
2046  if ((us%m_to_Z_restart /= 0.0) .and. (us%m_to_Z_restart /= us%m_to_Z)) &
2047  z_rescale = us%m_to_Z / us%m_to_Z_restart
2048  i_t_rescale = 1.0
2049  if ((us%s_to_T_restart /= 0.0) .and. (us%s_to_T_restart /= us%s_to_T)) &
2050  i_t_rescale = us%s_to_T_restart / us%s_to_T
2051  z2_t_rescale = z_rescale**2*i_t_rescale
2052 
2053  if (z2_t_rescale /= 1.0) then
2054  if (associated(visc%Kd_shear)) then ; if (query_initialized(visc%Kd_shear, "Kd_shear", restart_cs)) then
2055  do k=1,nz+1 ; do j=js,je ; do i=is,ie
2056  visc%Kd_shear(i,j,k) = z2_t_rescale * visc%Kd_shear(i,j,k)
2057  enddo ; enddo ; enddo
2058  endif ; endif
2059 
2060  if (associated(visc%Kv_shear)) then ; if (query_initialized(visc%Kv_shear, "Kv_shear", restart_cs)) then
2061  do k=1,nz+1 ; do j=js,je ; do i=is,ie
2062  visc%Kv_shear(i,j,k) = z2_t_rescale * visc%Kv_shear(i,j,k)
2063  enddo ; enddo ; enddo
2064  endif ; endif
2065 
2066  if (associated(visc%Kv_shear_Bu)) then ; if (query_initialized(visc%Kv_shear_Bu, "Kv_shear_Bu", restart_cs)) then
2067  do k=1,nz+1 ; do j=js-1,je ; do i=is-1,ie
2068  visc%Kv_shear_Bu(i,j,k) = z2_t_rescale * visc%Kv_shear_Bu(i,j,k)
2069  enddo ; enddo ; enddo
2070  endif ; endif
2071 
2072  if (associated(visc%Kv_slow)) then ; if (query_initialized(visc%Kv_slow, "Kv_slow", restart_cs)) then
2073  do k=1,nz+1 ; do j=js,je ; do i=is,ie
2074  visc%Kv_slow(i,j,k) = z2_t_rescale * visc%Kv_slow(i,j,k)
2075  enddo ; enddo ; enddo
2076  endif ; endif
2077  endif
2078 

◆ 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 1678 of file MOM_set_viscosity.F90.

1678  type(hor_index_type), intent(in) :: HI !< A horizontal index type structure.
1679  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1680  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1681  !! parameters.
1682  type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical
1683  !! viscosities and related fields.
1684  !! Allocated here.
1685  type(MOM_restart_CS), pointer :: restart_CS !< A pointer to the restart control structure.
1686  ! Local variables
1687  logical :: use_kappa_shear, KS_at_vertex
1688  logical :: adiabatic, useKPP, useEPBL
1689  logical :: use_CVMix_shear, MLE_use_PBL_MLD, use_CVMix_conv
1690  integer :: isd, ied, jsd, jed, nz
1691  real :: hfreeze !< If hfreeze > 0 [m], melt potential will be computed.
1692  character(len=40) :: mdl = "MOM_set_visc" ! This module's name.
1693  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke
1694 
1695  call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
1696  do_not_log=.true.)
1697 
1698  use_kappa_shear = .false. ; ks_at_vertex = .false. ; use_cvmix_shear = .false.
1699  usekpp = .false. ; useepbl = .false. ; use_cvmix_conv = .false.
1700 
1701  if (.not.adiabatic) then
1702  use_kappa_shear = kappa_shear_is_used(param_file)
1703  ks_at_vertex = kappa_shear_at_vertex(param_file)
1704  use_cvmix_shear = cvmix_shear_is_used(param_file)
1705  use_cvmix_conv = cvmix_conv_is_used(param_file)
1706  call get_param(param_file, mdl, "USE_KPP", usekpp, &
1707  "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "//&
1708  "to calculate diffusivities and non-local transport in the OBL.", &
1709  default=.false., do_not_log=.true.)
1710  call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", useepbl, &
1711  "If true, use an implied energetics planetary boundary "//&
1712  "layer scheme to determine the diffusivity and viscosity "//&
1713  "in the surface boundary layer.", default=.false., do_not_log=.true.)
1714  endif
1715 
1716  if (use_kappa_shear .or. usekpp .or. useepbl .or. use_cvmix_shear .or. use_cvmix_conv) then
1717  call safe_alloc_ptr(visc%Kd_shear, isd, ied, jsd, jed, nz+1)
1718  call register_restart_field(visc%Kd_shear, "Kd_shear", .false., restart_cs, &
1719  "Shear-driven turbulent diffusivity at interfaces", "m2 s-1", z_grid='i')
1720  endif
1721  if (usekpp .or. useepbl .or. use_cvmix_shear .or. use_cvmix_conv .or. &
1722  (use_kappa_shear .and. .not.ks_at_vertex )) then
1723  call safe_alloc_ptr(visc%Kv_shear, isd, ied, jsd, jed, nz+1)
1724  call register_restart_field(visc%Kv_shear, "Kv_shear", .false., restart_cs, &
1725  "Shear-driven turbulent viscosity at interfaces", "m2 s-1", z_grid='i')
1726  endif
1727  if (use_kappa_shear .and. ks_at_vertex) then
1728  call safe_alloc_ptr(visc%TKE_turb, hi%IsdB, hi%IedB, hi%JsdB, hi%JedB, nz+1)
1729  call safe_alloc_ptr(visc%Kv_shear_Bu, hi%IsdB, hi%IedB, hi%JsdB, hi%JedB, nz+1)
1730  call register_restart_field(visc%Kv_shear_Bu, "Kv_shear_Bu", .false., restart_cs, &
1731  "Shear-driven turbulent viscosity at vertex interfaces", "m2 s-1", &
1732  hor_grid="Bu", z_grid='i')
1733  elseif (use_kappa_shear) then
1734  call safe_alloc_ptr(visc%TKE_turb, isd, ied, jsd, jed, nz+1)
1735  endif
1736 
1737  ! MOM_bkgnd_mixing is always used, so always allocate visc%Kv_slow. GMM
1738  call safe_alloc_ptr(visc%Kv_slow, isd, ied, jsd, jed, nz+1)
1739  call register_restart_field(visc%Kv_slow, "Kv_slow", .false., restart_cs, &
1740  "Vertical turbulent viscosity at interfaces due to slow processes", &
1741  "m2 s-1", z_grid='i')
1742 
1743  ! visc%MLD is used to communicate the state of the (e)PBL or KPP to the rest of the model
1744  call get_param(param_file, mdl, "MLE_USE_PBL_MLD", mle_use_pbl_mld, &
1745  default=.false., do_not_log=.true.)
1746  ! visc%MLD needs to be allocated when melt potential is computed (HFREEZE>0)
1747  call get_param(param_file, mdl, "HFREEZE", hfreeze, &
1748  default=-1.0, do_not_log=.true.)
1749 
1750  if (mle_use_pbl_mld) then
1751  call safe_alloc_ptr(visc%MLD, isd, ied, jsd, jed)
1752  call register_restart_field(visc%MLD, "MLD", .false., restart_cs, &
1753  "Instantaneous active mixing layer depth", "m")
1754  endif
1755 
1756  if (hfreeze >= 0.0 .and. .not.mle_use_pbl_mld) then
1757  call safe_alloc_ptr(visc%MLD, isd, ied, jsd, jed)
1758  endif
1759 
1760 

◆ 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.

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 [m s-1].
[in]vThe meridional velocity [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 vertvisc_init.
[in]symmetrizeIf present and true, do extra calculations of those values in visc that would be calculated with symmetric memory.

Definition at line 110 of file MOM_set_viscosity.F90.

110  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
111  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
112  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
113  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
114  intent(in) :: u !< The zonal velocity [m s-1].
115  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
116  intent(in) :: v !< The meridional velocity [m s-1].
117  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
118  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
119  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any
120  !! available thermodynamic fields. Absent fields
121  !! have NULL ptrs..
122  type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
123  !! related fields.
124  type(set_visc_CS), pointer :: CS !< The control structure returned by a previous
125  !! call to vertvisc_init.
126  logical, optional, intent(in) :: symmetrize !< If present and true, do extra calculations
127  !! of those values in visc that would be
128  !! calculated with symmetric memory.
129 
130  ! Local variables
131  real, dimension(SZIB_(G)) :: &
132  ustar, & ! The bottom friction velocity [Z T-1 ~> m s-1].
133  T_EOS, & ! The temperature used to calculate the partial derivatives
134  ! of density with T and S [degC].
135  s_eos, & ! The salinity used to calculate the partial derivatives
136  ! of density with T and S [ppt].
137  dr_dt, & ! Partial derivative of the density in the bottom boundary
138  ! layer with temperature [kg m-3 degC-1].
139  dr_ds, & ! Partial derivative of the density in the bottom boundary
140  ! layer with salinity [kg m-3 ppt-1].
141  press ! The pressure at which dR_dT and dR_dS are evaluated [Pa].
142  real :: htot ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
143  real :: htot_vel ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
144 
145  real :: Rhtot ! Running sum of thicknesses times the layer potential
146  ! densities [H kg m-3 ~> kg m-2 or kg2 m-5].
147  real, dimension(SZIB_(G),SZJ_(G)) :: &
148  D_u, & ! Bottom depth interpolated to u points [Z ~> m].
149  mask_u ! A mask that disables any contributions from u points that
150  ! are land or past open boundary conditions [nondim], 0 or 1.
151  real, dimension(SZI_(G),SZJB_(G)) :: &
152  D_v, & ! Bottom depth interpolated to v points [Z ~> m].
153  mask_v ! A mask that disables any contributions from v points that
154  ! are land or past open boundary conditions [nondim], 0 or 1.
155  real, dimension(SZIB_(G),SZK_(G)) :: &
156  h_at_vel, & ! Layer thickness at a velocity point, using an upwind-biased
157  ! second order accurate estimate based on the previous velocity
158  ! direction [H ~> m or kg m-2].
159  h_vel, & ! Arithmetic mean of the layer thicknesses adjacent to a
160  ! velocity point [H ~> m or kg m-2].
161  t_vel, & ! Arithmetic mean of the layer temperatures adjacent to a
162  ! velocity point [degC].
163  s_vel, & ! Arithmetic mean of the layer salinities adjacent to a
164  ! velocity point [ppt].
165  rml_vel ! Arithmetic mean of the layer coordinate densities adjacent
166  ! to a velocity point [kg m-3].
167 
168  real :: h_vel_pos ! The arithmetic mean thickness at a velocity point
169  ! plus H_neglect to avoid 0 values [H ~> m or kg m-2].
170  real :: ustarsq ! 400 times the square of ustar, times
171  ! Rho0 divided by G_Earth and the conversion
172  ! from m to thickness units [H kg m-3 ~> kg m-2 or kg2 m-5].
173  real :: cdrag_sqrt_Z ! Square root of the drag coefficient, times a unit conversion
174  ! factor from lateral lengths to vertical depths [Z m-1 ~> 1].
175  real :: cdrag_sqrt ! Square root of the drag coefficient [nondim].
176  real :: oldfn ! The integrated energy required to
177  ! entrain up to the bottom of the layer,
178  ! divided by G_Earth [H kg m-3 ~> kg m-2 or kg2 m-5].
179  real :: Dfn ! The increment in oldfn for entraining
180  ! the layer [H kg m-3 ~> kg m-2 or kg2 m-5].
181  real :: Dh ! The increment in layer thickness from
182  ! the present layer [H ~> m or kg m-2].
183  real :: bbl_thick ! The thickness of the bottom boundary layer [H ~> m or kg m-2].
184  real :: bbl_thick_Z ! The thickness of the bottom boundary layer [Z ~> m].
185  real :: C2f ! C2f = 2*f at velocity points [T-1 ~> s-1].
186 
187  real :: U_bg_sq ! The square of an assumed background
188  ! velocity, for calculating the mean
189  ! magnitude near the bottom for use in the
190  ! quadratic bottom drag [m2 s-2].
191  real :: hwtot ! Sum of the thicknesses used to calculate
192  ! the near-bottom velocity magnitude [H ~> m or kg m-2].
193  real :: hutot ! Running sum of thicknesses times the
194  ! velocity magnitudes [H m s-1 ~> m2 s-1 or kg m-1 s-1].
195  real :: Thtot ! Running sum of thickness times temperature [degC H ~> degC m or degC kg m-2].
196  real :: Shtot ! Running sum of thickness times salinity [ppt H ~> ppt m or ppt kg m-2].
197  real :: hweight ! The thickness of a layer that is within Hbbl
198  ! of the bottom [H ~> m or kg m-2].
199  real :: v_at_u, u_at_v ! v at a u point or vice versa [m s-1].
200  real :: Rho0x400_G ! 400*Rho0/G_Earth, times unit conversion factors
201  ! [kg T2 H m-3 Z-2 ~> kg s2 m-4 or kg2 s2 m-7].
202  ! The 400 is a constant proposed by Killworth and Edwards, 1999.
203  real, dimension(SZI_(G),SZJ_(G),max(GV%nk_rho_varies,1)) :: &
204  Rml ! The mixed layer coordinate density [kg m-3].
205  real :: p_ref(SZI_(G)) ! The pressure used to calculate the coordinate
206  ! density [Pa] (usually set to 2e7 Pa = 2000 dbar).
207 
208  real :: D_vel ! The bottom depth at a velocity point [H ~> m or kg m-2].
209  real :: Dp, Dm ! The depths at the edges of a velocity cell [H ~> m or kg m-2].
210  real :: a ! a is the curvature of the bottom depth across a
211  ! cell, times the cell width squared [H ~> m or kg m-2].
212  real :: a_3, a_12 ! a/3 and a/12 [H ~> m or kg m-2].
213  real :: C24_a ! 24/a [H-1 ~> m-1 or m2 kg-1].
214  real :: slope ! The absolute value of the bottom depth slope across
215  ! a cell times the cell width [H ~> m or kg m-2].
216  real :: apb_4a, ax2_3apb ! Various nondimensional ratios of a and slope.
217  real :: a2x48_apb3, Iapb, Ibma_2 ! Combinations of a and slope [H-1 ~> m-1 or m2 kg-1].
218  ! All of the following "volumes" have units of thickness because they are normalized
219  ! by the full horizontal area of a velocity cell.
220  real :: Vol_open ! The cell volume above which it is open [H ~> m or kg m-2].
221  real :: Vol_direct ! With less than Vol_direct [H ~> m or kg m-2], there is a direct
222  ! solution of a cubic equation for L.
223  real :: Vol_2_reg ! The cell volume above which there are two separate
224  ! open areas that must be integrated [H ~> m or kg m-2].
225  real :: vol ! The volume below the interface whose normalized
226  ! width is being sought [H ~> m or kg m-2].
227  real :: vol_below ! The volume below the interface below the one that
228  ! is currently under consideration [H ~> m or kg m-2].
229  real :: Vol_err ! The error in the volume with the latest estimate of
230  ! L, or the error for the interface below [H ~> m or kg m-2].
231  real :: Vol_quit ! The volume error below which to quit iterating [H ~> m or kg m-2].
232  real :: Vol_tol ! A volume error tolerance [H ~> m or kg m-2].
233  real :: L(SZK_(G)+1) ! The fraction of the full cell width that is open at
234  ! the depth of each interface, nondimensional.
235  real :: L_direct ! The value of L above volume Vol_direct [nondim].
236  real :: L_max, L_min ! Upper and lower bounds on the correct value for L.
237  real :: Vol_err_max ! The volume errors for the upper and lower bounds on
238  real :: Vol_err_min ! the correct value for L [H ~> m or kg m-2].
239  real :: Vol_0 ! A deeper volume with known width L0 [H ~> m or kg m-2].
240  real :: L0 ! The value of L above volume Vol_0 [nondim].
241  real :: dVol ! vol - Vol_0 [H ~> m or kg m-2].
242  real :: dV_dL2 ! The partial derivative of volume with L squared
243  ! evaluated at L=L0 [H ~> m or kg m-2].
244  real :: h_neglect ! A thickness that is so small it is usually lost
245  ! in roundoff and can be neglected [H ~> m or kg m-2].
246  real :: ustH ! ustar converted to units of H T-1 [H T-1 ~> m s-1 or kg m-2 s-1].
247  real :: root ! A temporary variable [H T-1 ~> m s-1 or kg m-2 s-1].
248 
249  real :: Cell_width ! The transverse width of the velocity cell [m].
250  real :: Rayleigh ! A nondimensional value that is multiplied by the layer's
251  ! velocity magnitude to give the Rayleigh drag velocity, times
252  ! a lateral to vertical distance conversion factor [Z L-1 ~> 1].
253  real :: gam ! The ratio of the change in the open interface width
254  ! to the open interface width atop a cell [nondim].
255  real :: BBL_frac ! The fraction of a layer's drag that goes into the
256  ! viscous bottom boundary layer [nondim].
257  real :: BBL_visc_frac ! The fraction of all the drag that is expressed as
258  ! a viscous bottom boundary layer [nondim].
259  real, parameter :: C1_3 = 1.0/3.0, c1_6 = 1.0/6.0, c1_12 = 1.0/12.0
260  real :: C2pi_3 ! An irrational constant, 2/3 pi.
261  real :: tmp ! A temporary variable.
262  real :: tmp_val_m1_to_p1
263  logical :: use_BBL_EOS, do_i(SZIB_(G))
264  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, m, n, K2, nkmb, nkml
265  integer :: itt, maxitt=20
266  type(ocean_OBC_type), pointer :: OBC => null()
267 
268  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
269  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
270  nkmb = gv%nk_rho_varies ; nkml = gv%nkml
271  h_neglect = gv%H_subroundoff
272  rho0x400_g = 400.0*(gv%Rho0 / (us%L_to_Z**2 * gv%g_Earth)) * gv%Z_to_H
273  vol_quit = 0.9*gv%Angstrom_H + h_neglect
274  c2pi_3 = 8.0*atan(1.0)/3.0
275 
276  if (.not.associated(cs)) call mom_error(fatal,"MOM_vert_friction(BBL): "//&
277  "Module must be initialized before it is used.")
278  if (.not.cs%bottomdraglaw) return
279 
280  if (present(symmetrize)) then ; if (symmetrize) then
281  jsq = js-1 ; isq = is-1
282  endif ; endif
283 
284  if (cs%debug) then
285  call uvchksum("Start set_viscous_BBL [uv]", u, v, g%HI, haloshift=1)
286  call hchksum(h,"Start set_viscous_BBL h", g%HI, haloshift=1, scale=gv%H_to_m)
287  if (associated(tv%T)) call hchksum(tv%T, "Start set_viscous_BBL T", g%HI, haloshift=1)
288  if (associated(tv%S)) call hchksum(tv%S, "Start set_viscous_BBL S", g%HI, haloshift=1)
289  endif
290 
291  use_bbl_eos = associated(tv%eqn_of_state) .and. cs%BBL_use_EOS
292  obc => cs%OBC
293 
294  u_bg_sq = cs%drag_bg_vel * cs%drag_bg_vel
295  cdrag_sqrt = sqrt(cs%cdrag)
296  cdrag_sqrt_z = us%m_to_Z * sqrt(cs%cdrag)
297  k2 = max(nkmb+1, 2)
298 
299 ! With a linear drag law, the friction velocity is already known.
300 ! if (CS%linear_drag) ustar(:) = cdrag_sqrt_Z*CS%drag_bg_vel
301 
302  if ((nkml>0) .and. .not.use_bbl_eos) then
303  do i=isq,ieq+1 ; p_ref(i) = tv%P_ref ; enddo
304  !$OMP parallel do default(shared)
305  do j=jsq,jeq+1 ; do k=1,nkmb
306  call calculate_density(tv%T(:,j,k), tv%S(:,j,k), p_ref, &
307  rml(:,j,k), isq, ieq-isq+2, tv%eqn_of_state)
308  enddo ; enddo
309  endif
310 
311  !$OMP parallel do default(shared)
312  do j=js-1,je ; do i=is-1,ie+1
313  d_v(i,j) = 0.5*(g%bathyT(i,j) + g%bathyT(i,j+1))
314  mask_v(i,j) = g%mask2dCv(i,j)
315  enddo ; enddo
316  !$OMP parallel do default(shared)
317  do j=js-1,je+1 ; do i=is-1,ie
318  d_u(i,j) = 0.5*(g%bathyT(i,j) + g%bathyT(i+1,j))
319  mask_u(i,j) = g%mask2dCu(i,j)
320  enddo ; enddo
321 
322  if (associated(obc)) then ; do n=1,obc%number_of_segments
323  if (.not. obc%segment(n)%on_pe) cycle
324  ! Use a one-sided projection of bottom depths at OBC points.
325  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
326  if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je)) then
327  do i = max(is-1,obc%segment(n)%HI%isd), min(ie+1,obc%segment(n)%HI%ied)
328  if (obc%segment(n)%direction == obc_direction_n) d_v(i,j) = g%bathyT(i,j)
329  if (obc%segment(n)%direction == obc_direction_s) d_v(i,j) = g%bathyT(i,j+1)
330  enddo
331  elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ie)) then
332  do j = max(js-1,obc%segment(n)%HI%jsd), min(je+1,obc%segment(n)%HI%jed)
333  if (obc%segment(n)%direction == obc_direction_e) d_u(i,j) = g%bathyT(i,j)
334  if (obc%segment(n)%direction == obc_direction_w) d_u(i,j) = g%bathyT(i+1,j)
335  enddo
336  endif
337  enddo ; endif
338  if (associated(obc)) then ; do n=1,obc%number_of_segments
339  ! Now project bottom depths across cell-corner points in the OBCs. The two
340  ! projections have to occur in sequence and can not be combined easily.
341  if (.not. obc%segment(n)%on_pe) cycle
342  ! Use a one-sided projection of bottom depths at OBC points.
343  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
344  if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je)) then
345  do i = max(is-1,obc%segment(n)%HI%IsdB), min(ie,obc%segment(n)%HI%IedB)
346  if (obc%segment(n)%direction == obc_direction_n) then
347  d_u(i,j+1) = d_u(i,j) ; mask_u(i,j+1) = 0.0
348  elseif (obc%segment(n)%direction == obc_direction_s) then
349  d_u(i,j) = d_u(i,j+1) ; mask_u(i,j) = 0.0
350  endif
351  enddo
352  elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= ie)) then
353  do j = max(js-1,obc%segment(n)%HI%JsdB), min(je,obc%segment(n)%HI%JedB)
354  if (obc%segment(n)%direction == obc_direction_e) then
355  d_v(i+1,j) = d_v(i,j) ; mask_v(i+1,j) = 0.0
356  elseif (obc%segment(n)%direction == obc_direction_w) then
357  d_v(i,j) = d_v(i+1,j) ; mask_v(i,j) = 0.0
358  endif
359  enddo
360  endif
361  enddo ; endif
362 
363  if (.not.use_bbl_eos) rml_vel(:,:) = 0.0
364 
365  !$OMP parallel do default(private) shared(u,v,h,tv,visc,G,GV,US,CS,Rml,is,ie,js,je,nz,nkmb, &
366  !$OMP nkml,Isq,Ieq,Jsq,Jeq,h_neglect,Rho0x400_G,C2pi_3, &
367  !$OMP U_bg_sq,cdrag_sqrt_Z,cdrag_sqrt,K2,use_BBL_EOS, &
368  !$OMP OBC,maxitt,Vol_quit,D_u,D_v,mask_u,mask_v)
369  do j=jsq,jeq ; do m=1,2
370 
371  if (m==1) then
372  if (j<g%Jsc) cycle
373  is = isq ; ie = ieq
374  do i=is,ie
375  do_i(i) = .false.
376  if (g%mask2dCu(i,j) > 0) do_i(i) = .true.
377  enddo
378  else
379  is = g%isc ; ie = g%iec
380  do i=is,ie
381  do_i(i) = .false.
382  if (g%mask2dCv(i,j) > 0) do_i(i) = .true.
383  enddo
384  endif
385 
386  if (m==1) then
387  do k=1,nz ; do i=is,ie
388  if (do_i(i)) then
389  if (u(i,j,k) *(h(i+1,j,k) - h(i,j,k)) >= 0) then
390  h_at_vel(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / &
391  (h(i,j,k) + h(i+1,j,k) + h_neglect)
392  else
393  h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
394  endif
395  endif
396  h_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
397  enddo ; enddo
398  if (use_bbl_eos) then ; do k=1,nz ; do i=is,ie
399  ! Perhaps these should be thickness weighted.
400  t_vel(i,k) = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
401  s_vel(i,k) = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
402  enddo ; enddo ; else ; do k=1,nkmb ; do i=is,ie
403  rml_vel(i,k) = 0.5 * (rml(i,j,k) + rml(i+1,j,k))
404  enddo ; enddo ; endif
405  else
406  do k=1,nz ; do i=is,ie
407  if (do_i(i)) then
408  if (v(i,j,k) * (h(i,j+1,k) - h(i,j,k)) >= 0) then
409  h_at_vel(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / &
410  (h(i,j,k) + h(i,j+1,k) + h_neglect)
411  else
412  h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
413  endif
414  endif
415  h_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
416  enddo ; enddo
417  if (use_bbl_eos) then ; do k=1,nz ; do i=is,ie
418  t_vel(i,k) = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
419  s_vel(i,k) = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
420  enddo ; enddo ; else ; do k=1,nkmb ; do i=is,ie
421  rml_vel(i,k) = 0.5 * (rml(i,j,k) + rml(i,j+1,k))
422  enddo ; enddo ; endif
423  endif
424 
425  if (associated(obc)) then ; if (obc%number_of_segments > 0) then
426  ! Apply a zero gradient projection of thickness across OBC points.
427  if (m==1) then
428  do i=is,ie ; if (do_i(i) .and. (obc%segnum_u(i,j) /= obc_none)) then
429  if (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_e) then
430  do k=1,nz
431  h_at_vel(i,k) = h(i,j,k) ; h_vel(i,k) = h(i,j,k)
432  enddo
433  if (use_bbl_eos) then
434  do k=1,nz
435  t_vel(i,k) = tv%T(i,j,k) ; s_vel(i,k) = tv%S(i,j,k)
436  enddo
437  else
438  do k=1,nkmb
439  rml_vel(i,k) = rml(i,j,k)
440  enddo
441  endif
442  elseif (obc%segment(obc%segnum_u(i,j))%direction == obc_direction_w) then
443  do k=1,nz
444  h_at_vel(i,k) = h(i+1,j,k) ; h_vel(i,k) = h(i+1,j,k)
445  enddo
446  if (use_bbl_eos) then
447  do k=1,nz
448  t_vel(i,k) = tv%T(i+1,j,k) ; s_vel(i,k) = tv%S(i+1,j,k)
449  enddo
450  else
451  do k=1,nkmb
452  rml_vel(i,k) = rml(i+1,j,k)
453  enddo
454  endif
455  endif
456  endif ; enddo
457  else
458  do i=is,ie ; if (do_i(i) .and. (obc%segnum_v(i,j) /= obc_none)) then
459  if (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_n) then
460  do k=1,nz
461  h_at_vel(i,k) = h(i,j,k) ; h_vel(i,k) = h(i,j,k)
462  enddo
463  if (use_bbl_eos) then
464  do k=1,nz
465  t_vel(i,k) = tv%T(i,j,k) ; s_vel(i,k) = tv%S(i,j,k)
466  enddo
467  else
468  do k=1,nkmb
469  rml_vel(i,k) = rml(i,j,k)
470  enddo
471  endif
472  elseif (obc%segment(obc%segnum_v(i,j))%direction == obc_direction_s) then
473  do k=1,nz
474  h_at_vel(i,k) = h(i,j+1,k) ; h_vel(i,k) = h(i,j+1,k)
475  enddo
476  if (use_bbl_eos) then
477  do k=1,nz
478  t_vel(i,k) = tv%T(i,j+1,k) ; s_vel(i,k) = tv%S(i,j+1,k)
479  enddo
480  else
481  do k=1,nkmb
482  rml_vel(i,k) = rml(i,j+1,k)
483  enddo
484  endif
485  endif
486  endif ; enddo
487  endif
488  endif ; endif
489 
490  if (use_bbl_eos .or. .not.cs%linear_drag) then
491  do i=is,ie ; if (do_i(i)) then
492 ! This block of code calculates the mean velocity magnitude over
493 ! the bottommost CS%Hbbl of the water column for determining
494 ! the quadratic bottom drag.
495  htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
496  thtot = 0.0 ; shtot = 0.0
497  do k=nz,1,-1
498 
499  if (htot_vel>=cs%Hbbl) exit ! terminate the k loop
500 
501  hweight = min(cs%Hbbl - htot_vel, h_at_vel(i,k))
502  if (hweight < 1.5*gv%Angstrom_H + h_neglect) cycle
503 
504  htot_vel = htot_vel + h_at_vel(i,k)
505  hwtot = hwtot + hweight
506 
507  if ((.not.cs%linear_drag) .and. (hweight >= 0.0)) then ; if (m==1) then
508  v_at_u = set_v_at_u(v, h, g, i, j, k, mask_v, obc)
509  hutot = hutot + hweight * sqrt(u(i,j,k)*u(i,j,k) + &
510  v_at_u*v_at_u + u_bg_sq)
511  else
512  u_at_v = set_u_at_v(u, h, g, i, j, k, mask_u, obc)
513  hutot = hutot + hweight * sqrt(v(i,j,k)*v(i,j,k) + &
514  u_at_v*u_at_v + u_bg_sq)
515  endif ; endif
516 
517  if (use_bbl_eos .and. (hweight >= 0.0)) then
518  thtot = thtot + hweight * t_vel(i,k)
519  shtot = shtot + hweight * s_vel(i,k)
520  endif
521  enddo ! end of k loop
522 
523  if (.not.cs%linear_drag .and. (hwtot > 0.0)) then
524  ustar(i) = cdrag_sqrt_z*us%T_to_s*hutot/hwtot
525  else
526  ustar(i) = cdrag_sqrt_z*us%T_to_s*cs%drag_bg_vel
527  endif
528 
529  if (use_bbl_eos) then ; if (hwtot > 0.0) then
530  t_eos(i) = thtot/hwtot ; s_eos(i) = shtot/hwtot
531  else
532  t_eos(i) = 0.0 ; s_eos(i) = 0.0
533  endif ; endif
534  endif ; enddo
535  else
536  do i=is,ie ; ustar(i) = cdrag_sqrt_z*us%T_to_s*cs%drag_bg_vel ; enddo
537  endif ! Not linear_drag
538 
539  if (use_bbl_eos) then
540  do i=is,ie
541  press(i) = 0.0 ! or = forces%p_surf(i,j)
542  if (.not.do_i(i)) then ; t_eos(i) = 0.0 ; s_eos(i) = 0.0 ; endif
543  enddo
544  do k=1,nz ; do i=is,ie
545  press(i) = press(i) + gv%H_to_Pa * h_vel(i,k)
546  enddo ; enddo
547  call calculate_density_derivs(t_eos, s_eos, press, dr_dt, dr_ds, &
548  is-g%IsdB+1, ie-is+1, tv%eqn_of_state)
549  endif
550 
551  do i=is,ie ; if (do_i(i)) then
552 ! The 400.0 in this expression is the square of a constant proposed
553 ! by Killworth and Edwards, 1999, in equation (2.20).
554  ustarsq = rho0x400_g * ustar(i)**2
555  htot = 0.0
556 
557 ! This block of code calculates the thickness of a stratification
558 ! limited bottom boundary layer, using the prescription from
559 ! Killworth and Edwards, 1999, as described in Stephens and Hallberg
560 ! 2000 (unpublished and lost manuscript).
561  if (use_bbl_eos) then
562  thtot = 0.0 ; shtot = 0.0 ; oldfn = 0.0
563  do k=nz,2,-1
564  if (h_at_vel(i,k) <= 0.0) cycle
565 
566  oldfn = dr_dt(i)*(thtot - t_vel(i,k)*htot) + &
567  dr_ds(i)*(shtot - s_vel(i,k)*htot)
568  if (oldfn >= ustarsq) exit
569 
570  dfn = (dr_dt(i)*(t_vel(i,k) - t_vel(i,k-1)) + &
571  dr_ds(i)*(s_vel(i,k) - s_vel(i,k-1))) * &
572  (h_at_vel(i,k) + htot)
573 
574  if ((oldfn + dfn) <= ustarsq) then
575  dh = h_at_vel(i,k)
576  else
577  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
578  endif
579 
580  htot = htot + dh
581  thtot = thtot + t_vel(i,k)*dh ; shtot = shtot + s_vel(i,k)*dh
582  enddo
583  if ((oldfn < ustarsq) .and. h_at_vel(i,1) > 0.0) then
584  ! Layer 1 might be part of the BBL.
585  if (dr_dt(i) * (thtot - t_vel(i,1)*htot) + &
586  dr_ds(i) * (shtot - s_vel(i,1)*htot) < ustarsq) &
587  htot = htot + h_at_vel(i,1)
588  endif ! Examination of layer 1.
589  else ! Use Rlay and/or the coordinate density as density variables.
590  rhtot = 0.0
591  do k=nz,k2,-1
592  oldfn = rhtot - gv%Rlay(k)*htot
593  dfn = (gv%Rlay(k) - gv%Rlay(k-1))*(h_at_vel(i,k)+htot)
594 
595  if (oldfn >= ustarsq) then
596  cycle
597  elseif ((oldfn + dfn) <= ustarsq) then
598  dh = h_at_vel(i,k)
599  else
600  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
601  endif
602 
603  htot = htot + dh
604  rhtot = rhtot + gv%Rlay(k)*dh
605  enddo
606  if (nkml>0) then
607  do k=nkmb,2,-1
608  oldfn = rhtot - rml_vel(i,k)*htot
609  dfn = (rml_vel(i,k) - rml_vel(i,k-1)) * (h_at_vel(i,k)+htot)
610 
611  if (oldfn >= ustarsq) then
612  cycle
613  elseif ((oldfn + dfn) <= ustarsq) then
614  dh = h_at_vel(i,k)
615  else
616  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
617  endif
618 
619  htot = htot + dh
620  rhtot = rhtot + rml_vel(i,k)*dh
621  enddo
622  if (rhtot - rml_vel(i,1)*htot < ustarsq) htot = htot + h_at_vel(i,1)
623  else
624  if (rhtot - gv%Rlay(1)*htot < ustarsq) htot = htot + h_at_vel(i,1)
625  endif
626  endif ! use_BBL_EOS
627 
628 ! The Coriolis limit is 0.5*ustar/f. The buoyancy limit here is htot.
629 ! The bottom boundary layer thickness is found by solving the same
630 ! equation as in Killworth and Edwards: (h/h_f)^2 + h/h_N = 1.
631 
632  if (m==1) then ; c2f = g%CoriolisBu(i,j-1) + g%CoriolisBu(i,j)
633  else ; c2f = g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j) ; endif
634 
635  if (cs%cdrag * u_bg_sq <= 0.0) then
636  ! This avoids NaNs and overflows, and could be used in all cases,
637  ! but is not bitwise identical to the current code.
638  usth = ustar(i)*gv%Z_to_H ; root = sqrt(0.25*usth**2 + (htot*c2f)**2)
639  if (htot*usth <= (cs%BBL_thick_min+h_neglect) * (0.5*usth + root)) then
640  bbl_thick = cs%BBL_thick_min
641  else
642  bbl_thick = (htot * usth) / (0.5*usth + root)
643  endif
644  else
645  bbl_thick = htot / (0.5 + sqrt(0.25 + htot*htot*c2f*c2f/ &
646  ((ustar(i)*ustar(i)) * (gv%Z_to_H**2)) ) )
647 
648  if (bbl_thick < cs%BBL_thick_min) bbl_thick = cs%BBL_thick_min
649  endif
650 ! If there is Richardson number dependent mixing, that determines
651 ! the vertical extent of the bottom boundary layer, and there is no
652 ! need to set that scale here. In fact, viscously reducing the
653 ! shears over an excessively large region reduces the efficacy of
654 ! the Richardson number dependent mixing.
655  if ((bbl_thick > 0.5*cs%Hbbl) .and. (cs%RiNo_mix)) bbl_thick = 0.5*cs%Hbbl
656 
657  if (cs%Channel_drag) then
658  ! The drag within the bottommost bbl_thick is applied as a part of
659  ! an enhanced bottom viscosity, while above this the drag is applied
660  ! directly to the layers in question as a Rayleigh drag term.
661  if (m==1) then
662  d_vel = d_u(i,j)
663  tmp = g%mask2dCu(i,j+1) * d_u(i,j+1)
664  dp = 2.0 * d_vel * tmp / (d_vel + tmp)
665  tmp = g%mask2dCu(i,j-1) * d_u(i,j-1)
666  dm = 2.0 * d_vel * tmp / (d_vel + tmp)
667  else
668  d_vel = d_v(i,j)
669  tmp = g%mask2dCv(i+1,j) * d_v(i+1,j)
670  dp = 2.0 * d_vel * tmp / (d_vel + tmp)
671  tmp = g%mask2dCv(i-1,j) * d_v(i-1,j)
672  dm = 2.0 * d_vel * tmp / (d_vel + tmp)
673  endif
674  if (dm > dp) then ; tmp = dp ; dp = dm ; dm = tmp ; endif
675 
676  ! Convert the D's to the units of thickness.
677  dp = gv%Z_to_H*dp ; dm = gv%Z_to_H*dm ; d_vel = gv%Z_to_H*d_vel
678 
679  a_3 = (dp + dm - 2.0*d_vel) ; a = 3.0*a_3 ; a_12 = 0.25*a_3
680  slope = dp - dm
681  ! If the curvature is small enough, there is no reason not to assume
682  ! a uniformly sloping or flat bottom.
683  if (abs(a) < 1e-2*(slope + cs%BBL_thick_min)) a = 0.0
684  ! Each cell extends from x=-1/2 to 1/2, and has a topography
685  ! given by D(x) = a*x^2 + b*x + D - a/12.
686 
687  ! Calculate the volume above which the entire cell is open and the
688  ! other volumes at which the equation that is solved for L changes.
689  if (a > 0.0) then
690  if (slope >= a) then
691  vol_open = d_vel - dm ; vol_2_reg = vol_open
692  else
693  tmp = slope/a
694  vol_open = 0.25*slope*tmp + c1_12*a
695  vol_2_reg = 0.5*tmp**2 * (a - c1_3*slope)
696  endif
697  ! Define some combinations of a & b for later use.
698  c24_a = 24.0/a ; iapb = 1.0/(a+slope)
699  apb_4a = (slope+a)/(4.0*a) ; a2x48_apb3 = (48.0*(a*a))*(iapb**3)
700  ax2_3apb = 2.0*c1_3*a*iapb
701  elseif (a == 0.0) then
702  vol_open = 0.5*slope
703  if (slope > 0) iapb = 1.0/slope
704  else ! a < 0.0
705  vol_open = d_vel - dm
706  if (slope >= -a) then
707  iapb = 1.0e30 ; if (slope+a /= 0.0) iapb = 1.0/(a+slope)
708  vol_direct = 0.0 ; l_direct = 0.0 ; c24_a = 0.0
709  else
710  c24_a = 24.0/a ; iapb = 1.0/(a+slope)
711  l_direct = 1.0 + slope/a ! L_direct < 1 because a < 0
712  vol_direct = -c1_6*a*l_direct**3
713  endif
714  ibma_2 = 2.0 / (slope - a)
715  endif
716 
717  l(nz+1) = 0.0 ; vol = 0.0 ; vol_err = 0.0 ; bbl_visc_frac = 0.0
718  ! Determine the normalized open length at each interface.
719  do k=nz,1,-1
720  vol_below = vol
721 
722  vol = vol + h_vel(i,k)
723  h_vel_pos = h_vel(i,k) + h_neglect
724 
725  if (vol >= vol_open) then ; l(k) = 1.0
726  elseif (a == 0) then ! The bottom has no curvature.
727  l(k) = sqrt(2.0*vol*iapb)
728  elseif (a > 0) then
729  ! There may be a minimum depth, and there are
730  ! analytic expressions for L for all cases.
731  if (vol < vol_2_reg) then
732  ! In this case, there is a contiguous open region and
733  ! vol = 0.5*L^2*(slope + a/3*(3-4L)).
734  if (a2x48_apb3*vol < 1e-8) then ! Could be 1e-7?
735  ! There is a very good approximation here for massless layers.
736  l0 = sqrt(2.0*vol*iapb) ; l(k) = l0*(1.0 + ax2_3apb*l0)
737  else
738  l(k) = apb_4a * (1.0 - &
739  2.0 * cos(c1_3*acos(a2x48_apb3*vol - 1.0) - c2pi_3))
740  endif
741  ! To check the answers.
742  ! Vol_err = 0.5*(L(K)*L(K))*(slope + a_3*(3.0-4.0*L(K))) - vol
743  else ! There are two separate open regions.
744  ! vol = slope^2/4a + a/12 - (a/12)*(1-L)^2*(1+2L)
745  ! At the deepest volume, L = slope/a, at the top L = 1.
746  !L(K) = 0.5 - cos(C1_3*acos(1.0 - C24_a*(Vol_open - vol)) - C2pi_3)
747  tmp_val_m1_to_p1 = 1.0 - c24_a*(vol_open - vol)
748  tmp_val_m1_to_p1 = max(-1., min(1., tmp_val_m1_to_p1))
749  l(k) = 0.5 - cos(c1_3*acos(tmp_val_m1_to_p1) - c2pi_3)
750  ! To check the answers.
751  ! Vol_err = Vol_open - a_12*(1.0+2.0*L(K)) * (1.0-L(K))**2 - vol
752  endif
753  else ! a < 0.
754  if (vol <= vol_direct) then
755  ! Both edges of the cell are bounded by walls.
756  l(k) = (-0.25*c24_a*vol)**c1_3
757  else
758  ! x_R is at 1/2 but x_L is in the interior & L is found by solving
759  ! vol = 0.5*L^2*(slope + a/3*(3-4L))
760 
761  ! Vol_err = 0.5*(L(K+1)*L(K+1))*(slope + a_3*(3.0-4.0*L(K+1))) - vol_below
762  ! Change to ...
763  ! if (min(Vol_below + Vol_err, vol) <= Vol_direct) then ?
764  if (vol_below + vol_err <= vol_direct) then
765  l0 = l_direct ; vol_0 = vol_direct
766  else
767  l0 = l(k+1) ; vol_0 = vol_below + vol_err
768  ! Change to Vol_0 = min(Vol_below + Vol_err, vol) ?
769  endif
770 
771  ! Try a relatively simple solution that usually works well
772  ! for massless layers.
773  dv_dl2 = 0.5*(slope+a) - a*l0 ; dvol = (vol-vol_0)
774  ! dV_dL2 = 0.5*(slope+a) - a*L0 ; dVol = max(vol-Vol_0, 0.0)
775 
776  ! The following code is more robust when GV%Angstrom_H=0, but it changes answers.
777  if (.not.cs%answers_2018) then
778  vol_tol = max(0.5*gv%Angstrom_H + gv%H_subroundoff, 1e-14*vol)
779  vol_quit = max(0.9*gv%Angstrom_H + gv%H_subroundoff, 1e-14*vol)
780  endif
781 
782  if ((.not.cs%answers_2018) .and. (dvol <= 0.0)) then
783  l(k) = l0
784  vol_err = 0.5*(l(k)*l(k))*(slope + a_3*(3.0-4.0*l(k))) - vol
785  elseif ( ((.not.cs%answers_2018) .and. &
786  (a*a*dvol**3 < vol_tol*dv_dl2**2 *(dv_dl2*vol_tol - 2.0*a*l0*dvol))) .or. &
787  (cs%answers_2018 .and. (a*a*dvol**3 < gv%Angstrom_H*dv_dl2**2 * &
788  (0.25*dv_dl2*gv%Angstrom_H - a*l0*dvol) )) ) then
789  ! One iteration of Newton's method should give an estimate
790  ! that is accurate to within Vol_tol.
791  l(k) = sqrt(l0*l0 + dvol / dv_dl2)
792  vol_err = 0.5*(l(k)*l(k))*(slope + a_3*(3.0-4.0*l(k))) - vol
793  else
794  if (dv_dl2*(1.0-l0*l0) < dvol + &
795  dv_dl2 * (vol_open - vol)*ibma_2) then
796  l_max = sqrt(1.0 - (vol_open - vol)*ibma_2)
797  else
798  l_max = sqrt(l0*l0 + dvol / dv_dl2)
799  endif
800  l_min = sqrt(l0*l0 + dvol / (0.5*(slope+a) - a*l_max))
801 
802  vol_err_min = 0.5*(l_min**2)*(slope + a_3*(3.0-4.0*l_min)) - vol
803  vol_err_max = 0.5*(l_max**2)*(slope + a_3*(3.0-4.0*l_max)) - vol
804  ! if ((abs(Vol_err_min) <= Vol_quit) .or. (Vol_err_min >= Vol_err_max)) then
805  if (abs(vol_err_min) <= vol_quit) then
806  l(k) = l_min ; vol_err = vol_err_min
807  else
808  l(k) = sqrt((l_min**2*vol_err_max - l_max**2*vol_err_min) / &
809  (vol_err_max - vol_err_min))
810  do itt=1,maxitt
811  vol_err = 0.5*(l(k)*l(k))*(slope + a_3*(3.0-4.0*l(k))) - vol
812  if (abs(vol_err) <= vol_quit) exit
813  ! Take a Newton's method iteration. This equation has proven
814  ! robust enough not to need bracketing.
815  l(k) = l(k) - vol_err / (l(k)* (slope + a - 2.0*a*l(k)))
816  ! This would be a Newton's method iteration for L^2:
817  ! L(K) = sqrt(L(K)*L(K) - Vol_err / (0.5*(slope+a) - a*L(K)))
818  enddo
819  endif ! end of iterative solver
820  endif ! end of 1-boundary alternatives.
821  endif ! end of a<0 cases.
822  endif
823 
824  ! Determine the drag contributing to the bottom boundary layer
825  ! and the Raleigh drag that acting on each layer.
826  if (l(k) > l(k+1)) then
827  if (vol_below < bbl_thick) then
828  bbl_frac = (1.0-vol_below/bbl_thick)**2
829  bbl_visc_frac = bbl_visc_frac + bbl_frac*(l(k) - l(k+1))
830  else
831  bbl_frac = 0.0
832  endif
833 
834  if (m==1) then ; cell_width = g%dy_Cu(i,j)
835  else ; cell_width = g%dx_Cv(i,j) ; endif
836  gam = 1.0 - l(k+1)/l(k)
837  rayleigh = us%m_to_Z * cs%cdrag * (l(k)-l(k+1)) * (1.0-bbl_frac) * &
838  (12.0*cs%c_Smag*h_vel_pos) / (12.0*cs%c_Smag*h_vel_pos + &
839  gv%m_to_H * cs%cdrag * gam*(1.0-gam)*(1.0-1.5*gam) * l(k)**2 * cell_width)
840  else ! This layer feels no drag.
841  rayleigh = 0.0
842  endif
843 
844  if (m==1) then
845  if (rayleigh > 0.0) then
846  v_at_u = set_v_at_u(v, h, g, i, j, k, mask_v, obc)
847  visc%Ray_u(i,j,k) = rayleigh*us%T_to_s*sqrt(u(i,j,k)*u(i,j,k) + &
848  v_at_u*v_at_u + u_bg_sq)
849  else ; visc%Ray_u(i,j,k) = 0.0 ; endif
850  else
851  if (rayleigh > 0.0) then
852  u_at_v = set_u_at_v(u, h, g, i, j, k, mask_u, obc)
853  visc%Ray_v(i,j,k) = rayleigh*us%T_to_s*sqrt(v(i,j,k)*v(i,j,k) + &
854  u_at_v*u_at_v + u_bg_sq)
855  else ; visc%Ray_v(i,j,k) = 0.0 ; endif
856  endif
857 
858  enddo ! k loop to determine L(K).
859 
860  bbl_thick_z = bbl_thick * gv%H_to_Z
861  if (m==1) then
862  visc%Kv_bbl_u(i,j) = max(cs%Kv_BBL_min, &
863  cdrag_sqrt*ustar(i)*bbl_thick_z*bbl_visc_frac)
864  visc%bbl_thick_u(i,j) = bbl_thick_z
865  else
866  visc%Kv_bbl_v(i,j) = max(cs%Kv_BBL_min, &
867  cdrag_sqrt*ustar(i)*bbl_thick_z*bbl_visc_frac)
868  visc%bbl_thick_v(i,j) = bbl_thick_z
869  endif
870 
871  else ! Not Channel_drag.
872 ! Here the near-bottom viscosity is set to a value which will give
873 ! the correct stress when the shear occurs over bbl_thick.
874  bbl_thick_z = bbl_thick * gv%H_to_Z
875  if (m==1) then
876  visc%Kv_bbl_u(i,j) = max(cs%Kv_BBL_min, cdrag_sqrt*ustar(i)*bbl_thick_z)
877  visc%bbl_thick_u(i,j) = bbl_thick_z
878  else
879  visc%Kv_bbl_v(i,j) = max(cs%Kv_BBL_min, cdrag_sqrt*ustar(i)*bbl_thick_z)
880  visc%bbl_thick_v(i,j) = bbl_thick_z
881  endif
882  endif
883  endif ; enddo ! end of i loop
884  enddo ; enddo ! end of m & j loops
885 
886 ! Offer diagnostics for averaging
887  if (cs%id_bbl_thick_u > 0) &
888  call post_data(cs%id_bbl_thick_u, visc%bbl_thick_u, cs%diag)
889  if (cs%id_kv_bbl_u > 0) &
890  call post_data(cs%id_kv_bbl_u, visc%kv_bbl_u, cs%diag)
891  if (cs%id_bbl_thick_v > 0) &
892  call post_data(cs%id_bbl_thick_v, visc%bbl_thick_v, cs%diag)
893  if (cs%id_kv_bbl_v > 0) &
894  call post_data(cs%id_kv_bbl_v, visc%kv_bbl_v, cs%diag)
895  if (cs%id_Ray_u > 0) &
896  call post_data(cs%id_Ray_u, visc%Ray_u, cs%diag)
897  if (cs%id_Ray_v > 0) &
898  call post_data(cs%id_Ray_v, visc%Ray_v, cs%diag)
899 
900  if (cs%debug) then
901  if (associated(visc%Ray_u) .and. associated(visc%Ray_v)) &
902  call uvchksum("Ray [uv]", visc%Ray_u, visc%Ray_v, g%HI, haloshift=0, scale=us%Z_to_m)
903  if (associated(visc%kv_bbl_u) .and. associated(visc%kv_bbl_v)) &
904  call uvchksum("kv_bbl_[uv]", visc%kv_bbl_u, visc%kv_bbl_v, g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
905  if (associated(visc%bbl_thick_u) .and. associated(visc%bbl_thick_v)) &
906  call uvchksum("bbl_thick_[uv]", visc%bbl_thick_u, &
907  visc%bbl_thick_v, g%HI, haloshift=0, scale=us%Z_to_m)
908  endif
909 

◆ set_viscous_ml()

subroutine, public mom_set_visc::set_viscous_ml ( 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(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 [m s-1].
[in]vThe meridional velocity [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 [s].
csThe control structure returned by a previous call to vertvisc_init.
[in]symmetrizeIf present and true, do extra calculations of those values in visc that would be calculated with symmetric memory.

Definition at line 1004 of file MOM_set_viscosity.F90.

1004  type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
1005  type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1006  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1007  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1008  intent(in) :: u !< The zonal velocity [m s-1].
1009  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1010  intent(in) :: v !< The meridional velocity [m s-1].
1011  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1012  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2].
1013  type(thermo_var_ptrs), intent(in) :: tv !< A structure containing pointers to any available
1014  !! thermodynamic fields. Absent fields have
1015  !! NULL ptrs.
1016  type(mech_forcing), intent(in) :: forces !< A structure with the driving mechanical forces
1017  type(vertvisc_type), intent(inout) :: visc !< A structure containing vertical viscosities and
1018  !! related fields.
1019  real, intent(in) :: dt !< Time increment [s].
1020  type(set_visc_CS), pointer :: CS !< The control structure returned by a previous
1021  !! call to vertvisc_init.
1022  logical, optional, intent(in) :: symmetrize !< If present and true, do extra calculations
1023  !! of those values in visc that would be
1024  !! calculated with symmetric memory.
1025  ! Local variables
1026  real, dimension(SZIB_(G)) :: &
1027  htot, & ! The total depth of the layers being that are within the
1028  ! surface mixed layer [H ~> m or kg m-2].
1029  thtot, & ! The integrated temperature of layers that are within the
1030  ! surface mixed layer [H degC ~> m degC or kg degC m-2].
1031  shtot, & ! The integrated salt of layers that are within the
1032  ! surface mixed layer [H ppt ~> m ppt or kg ppt m-2].
1033  rhtot, & ! The integrated density of layers that are within the surface mixed layer
1034  ! [H kg m-3 ~> kg m-2 or kg2 m-5]. Rhtot is only used if no
1035  ! equation of state is used.
1036  uhtot, & ! The depth integrated zonal and meridional velocities within
1037  vhtot, & ! the surface mixed layer [H m s-1 ~> m2 s-1 or kg m-1 s-1].
1038  idecay_len_tke, & ! The inverse of a turbulence decay length scale [H-1 ~> m-1 or m2 kg-1].
1039  dr_dt, & ! Partial derivative of the density at the base of layer nkml
1040  ! (roughly the base of the mixed layer) with temperature [kg m-3 degC-1].
1041  dr_ds, & ! Partial derivative of the density at the base of layer nkml
1042  ! (roughly the base of the mixed layer) with salinity [kg m-3 ppt-1].
1043  ustar, & ! The surface friction velocity under ice shelves [Z T-1 ~> m s-1].
1044  press, & ! The pressure at which dR_dT and dR_dS are evaluated [Pa].
1045  t_eos, & ! The potential temperature at which dR_dT and dR_dS are evaluated [degC]
1046  s_eos ! The salinity at which dR_dT and dR_dS are evaluated [ppt].
1047  real, dimension(SZIB_(G),SZJ_(G)) :: &
1048  mask_u ! A mask that disables any contributions from u points that
1049  ! are land or past open boundary conditions [nondim], 0 or 1.
1050  real, dimension(SZI_(G),SZJB_(G)) :: &
1051  mask_v ! A mask that disables any contributions from v points that
1052  ! are land or past open boundary conditions [nondim], 0 or 1.
1053  real :: h_at_vel(SZIB_(G),SZK_(G))! Layer thickness at velocity points,
1054  ! using an upwind-biased second order accurate estimate based
1055  ! on the previous velocity direction [H ~> m or kg m-2].
1056  integer :: k_massive(SZIB_(G)) ! The k-index of the deepest layer yet found
1057  ! that has more than h_tiny thickness and will be in the
1058  ! viscous mixed layer.
1059  real :: Uh2 ! The squared magnitude of the difference between the velocity
1060  ! integrated through the mixed layer and the velocity of the
1061  ! interior layer layer times the depth of the the mixed layer
1062  ! [H2 Z2 T-2 ~> m4 s-2 or kg2 m-2 s-2].
1063  real :: htot_vel ! Sum of the layer thicknesses up to some point [H ~> m or kg m-2].
1064  real :: hwtot ! Sum of the thicknesses used to calculate
1065  ! the near-bottom velocity magnitude [H ~> m or kg m-2].
1066  real :: hutot ! Running sum of thicknesses times the
1067  ! velocity magnitudes [H m s-1 ~> m2 s-1 or kg m-1 s-1].
1068  real :: hweight ! The thickness of a layer that is within Hbbl
1069  ! of the bottom [H ~> m or kg m-2].
1070  real :: tbl_thick_Z ! The thickness of the top boundary layer [Z ~> m].
1071 
1072  real :: hlay ! The layer thickness at velocity points [H ~> m or kg m-2].
1073  real :: I_2hlay ! 1 / 2*hlay [H-1 ~> m-1 or m2 kg-1].
1074  real :: T_lay ! The layer temperature at velocity points [degC].
1075  real :: S_lay ! The layer salinity at velocity points [ppt].
1076  real :: Rlay ! The layer potential density at velocity points [kg m-3].
1077  real :: Rlb ! The potential density of the layer below [kg m-3].
1078  real :: v_at_u ! The meridonal velocity at a zonal velocity point [m s-1].
1079  real :: u_at_v ! The zonal velocity at a meridonal velocity point [m s-1].
1080  real :: gHprime ! The mixed-layer internal gravity wave speed squared, based
1081  ! on the mixed layer thickness and density difference across
1082  ! the base of the mixed layer [L2 T-2 ~> m2 s-2].
1083  real :: RiBulk ! The bulk Richardson number below which water is in the
1084  ! viscous mixed layer, including reduction for turbulent
1085  ! decay. Nondimensional.
1086  real :: dt_Rho0 ! The time step divided by the conversion from the layer
1087  ! thickness to layer mass [s H m2 kg-1 ~> s m3 kg-1 or s].
1088  real :: g_H_Rho0 ! The gravitational acceleration times the conversion from H to m divided
1089  ! by the mean density [L2 m3 T-2 H-1 kg-1 ~> m4 s-2 kg-1 or m7 s-2 kg-2].
1090  real :: ustarsq ! 400 times the square of ustar, times
1091  ! Rho0 divided by G_Earth and the conversion
1092  ! from m to thickness units [H kg m-3 ~> kg m-2 or kg2 m-5].
1093  real :: cdrag_sqrt_Z ! Square root of the drag coefficient, times a unit conversion
1094  ! factor from lateral lengths to vertical depths [Z m-1 ~> 1].
1095  real :: cdrag_sqrt ! Square root of the drag coefficient [nondim].
1096  real :: oldfn ! The integrated energy required to
1097  ! entrain up to the bottom of the layer,
1098  ! divided by G_Earth [H kg m-3 ~> kg m-2 or kg2 m-5].
1099  real :: Dfn ! The increment in oldfn for entraining
1100  ! the layer [H kg m-3 ~> kg m-2 or kg2 m-5].
1101  real :: Dh ! The increment in layer thickness from
1102  ! the present layer [H ~> m or kg m-2].
1103  real :: U_bg_sq ! The square of an assumed background velocity, for
1104  ! calculating the mean magnitude near the top for use in
1105  ! the quadratic surface drag [m2 s-2].
1106  real :: h_tiny ! A very small thickness [H ~> m or kg m-2]. Layers that are less than
1107  ! h_tiny can not be the deepest in the viscous mixed layer.
1108  real :: absf ! The absolute value of f averaged to velocity points [T-1 ~> s-1].
1109  real :: U_star ! The friction velocity at velocity points [Z T-1 ~> m s-1].
1110  real :: h_neglect ! A thickness that is so small it is usually lost
1111  ! in roundoff and can be neglected [H ~> m or kg m-2].
1112  real :: Rho0x400_G ! 400*Rho0/G_Earth, times unit conversion factors
1113  ! [kg T2 H m-3 Z-2 ~> kg s2 m-4 or kg2 s2 m-7].
1114  ! The 400 is a constant proposed by Killworth and Edwards, 1999.
1115  real :: ustar1 ! ustar [H T-1 ~> m s-1 or kg m-2 s-1]
1116  real :: h2f2 ! (h*2*f)^2 [H2 T-2 ~> m2 s-2 or kg2 m-4 s-2]
1117  logical :: use_EOS, do_any, do_any_shelf, do_i(SZIB_(G))
1118  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, K2, nkmb, nkml, n
1119  type(ocean_OBC_type), pointer :: OBC => null()
1120 
1121  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1122  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1123  nkmb = gv%nk_rho_varies ; nkml = gv%nkml
1124 
1125  if (.not.associated(cs)) call mom_error(fatal,"MOM_vert_friction(visc_ML): "//&
1126  "Module must be initialized before it is used.")
1127  if (.not.(cs%dynamic_viscous_ML .or. associated(forces%frac_shelf_u) .or. &
1128  associated(forces%frac_shelf_v)) ) return
1129 
1130  if (present(symmetrize)) then ; if (symmetrize) then
1131  jsq = js-1 ; isq = is-1
1132  endif ; endif
1133 
1134  rho0x400_g = 400.0*(gv%Rho0/(us%L_to_Z**2 * gv%g_Earth)) * gv%Z_to_H
1135  u_bg_sq = cs%drag_bg_vel * cs%drag_bg_vel
1136  cdrag_sqrt = sqrt(cs%cdrag)
1137  cdrag_sqrt_z = us%m_to_Z * sqrt(cs%cdrag)
1138 
1139  obc => cs%OBC
1140  use_eos = associated(tv%eqn_of_state)
1141  dt_rho0 = dt/gv%H_to_kg_m2
1142  h_neglect = gv%H_subroundoff
1143  h_tiny = 2.0*gv%Angstrom_H + h_neglect
1144  g_h_rho0 = (gv%g_Earth*gv%H_to_Z) / gv%Rho0
1145 
1146  if (associated(forces%frac_shelf_u) .neqv. associated(forces%frac_shelf_v)) &
1147  call mom_error(fatal, "set_viscous_ML: one of forces%frac_shelf_u and "//&
1148  "forces%frac_shelf_v is associated, but the other is not.")
1149 
1150  if (associated(forces%frac_shelf_u)) then
1151  ! This configuration has ice shelves, and the appropriate variables need to
1152  ! be allocated.
1153  call safe_alloc_ptr(visc%tauy_shelf, g%isd, g%ied, g%JsdB, g%JedB)
1154  call safe_alloc_ptr(visc%tbl_thick_shelf_u, g%IsdB, g%IedB, g%jsd, g%jed)
1155  call safe_alloc_ptr(visc%tbl_thick_shelf_v, g%isd, g%ied, g%JsdB, g%JedB)
1156  call safe_alloc_ptr(visc%kv_tbl_shelf_u, g%IsdB, g%IedB, g%jsd, g%jed)
1157  call safe_alloc_ptr(visc%kv_tbl_shelf_v, g%isd, g%ied, g%JsdB, g%JedB)
1158 
1159  ! With a linear drag law, the friction velocity is already known.
1160 ! if (CS%linear_drag) ustar(:) = cdrag_sqrt_Z*CS%drag_bg_vel
1161  endif
1162 
1163  !$OMP parallel do default(shared)
1164  do j=js-1,je ; do i=is-1,ie+1
1165  mask_v(i,j) = g%mask2dCv(i,j)
1166  enddo ; enddo
1167  !$OMP parallel do default(shared)
1168  do j=js-1,je+1 ; do i=is-1,ie
1169  mask_u(i,j) = g%mask2dCu(i,j)
1170  enddo ; enddo
1171 
1172  if (associated(obc)) then ; do n=1,obc%number_of_segments
1173  ! Now project bottom depths across cell-corner points in the OBCs. The two
1174  ! projections have to occur in sequence and can not be combined easily.
1175  if (.not. obc%segment(n)%on_pe) cycle
1176  ! Use a one-sided projection of bottom depths at OBC points.
1177  i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
1178  if (obc%segment(n)%is_N_or_S .and. (j >= js-1) .and. (j <= je)) then
1179  do i = max(is-1,obc%segment(n)%HI%IsdB), min(ie,obc%segment(n)%HI%IedB)
1180  if (obc%segment(n)%direction == obc_direction_n) mask_u(i,j+1) = 0.0
1181  if (obc%segment(n)%direction == obc_direction_s) mask_u(i,j) = 0.0
1182  enddo
1183  elseif (obc%segment(n)%is_E_or_W .and. (i >= is-1) .and. (i <= je)) then
1184  do j = max(js-1,obc%segment(n)%HI%JsdB), min(je,obc%segment(n)%HI%JedB)
1185  if (obc%segment(n)%direction == obc_direction_e) mask_v(i+1,j) = 0.0
1186  if (obc%segment(n)%direction == obc_direction_w) mask_v(i,j) = 0.0
1187  enddo
1188  endif
1189  enddo ; endif
1190 
1191  !$OMP parallel do default(private) shared(u,v,h,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, &
1192  !$OMP h_neglect,h_tiny,g_H_Rho0,js,je,OBC,Isq,Ieq,nz, &
1193  !$OMP U_bg_sq,mask_v,cdrag_sqrt,cdrag_sqrt_Z,Rho0x400_G,nkml)
1194  do j=js,je ! u-point loop
1195  if (cs%dynamic_viscous_ML) then
1196  do_any = .false.
1197  do i=isq,ieq
1198  htot(i) = 0.0
1199  if (g%mask2dCu(i,j) < 0.5) then
1200  do_i(i) = .false. ; visc%nkml_visc_u(i,j) = nkml
1201  else
1202  do_i(i) = .true. ; do_any = .true.
1203  k_massive(i) = nkml
1204  thtot(i) = 0.0 ; shtot(i) = 0.0 ; rhtot(i) = 0.0
1205  uhtot(i) = dt_rho0 * forces%taux(i,j)
1206  vhtot(i) = 0.25 * dt_rho0 * ((forces%tauy(i,j) + forces%tauy(i+1,j-1)) + &
1207  (forces%tauy(i,j-1) + forces%tauy(i+1,j)))
1208 
1209  if (cs%omega_frac >= 1.0) then ; absf = 2.0*cs%omega ; else
1210  absf = 0.5*(abs(g%CoriolisBu(i,j)) + abs(g%CoriolisBu(i,j-1)))
1211  if (cs%omega_frac > 0.0) &
1212  absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
1213  endif
1214  u_star = max(cs%ustar_min, 0.5 * (forces%ustar(i,j) + forces%ustar(i+1,j)))
1215  idecay_len_tke(i) = ((absf / u_star) * cs%TKE_decay) * gv%H_to_Z
1216  endif
1217  enddo
1218 
1219  if (do_any) then ; do k=1,nz
1220  if (k > nkml) then
1221  do_any = .false.
1222  if (use_eos .and. (k==nkml+1)) then
1223  ! Find dRho/dT and dRho_dS.
1224  do i=isq,ieq
1225  press(i) = gv%H_to_Pa * htot(i)
1226  k2 = max(1,nkml)
1227  i_2hlay = 1.0 / (h(i,j,k2) + h(i+1,j,k2) + h_neglect)
1228  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
1229  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
1230  enddo
1231  call calculate_density_derivs(t_eos, s_eos, press, dr_dt, dr_ds, &
1232  isq-g%IsdB+1, ieq-isq+1, tv%eqn_of_state)
1233  endif
1234 
1235  do i=isq,ieq ; if (do_i(i)) then
1236 
1237  hlay = 0.5*(h(i,j,k) + h(i+1,j,k))
1238  if (hlay > h_tiny) then ! Only consider non-vanished layers.
1239  i_2hlay = 1.0 / (h(i,j,k) + h(i+1,j,k))
1240  v_at_u = 0.5 * (h(i,j,k) * (v(i,j,k) + v(i,j-1,k)) + &
1241  h(i+1,j,k) * (v(i+1,j,k) + v(i+1,j-1,k))) * i_2hlay
1242  uh2 = us%m_s_to_L_T**2*((uhtot(i) - htot(i)*u(i,j,k))**2 + (vhtot(i) - htot(i)*v_at_u)**2)
1243 
1244  if (use_eos) then
1245  t_lay = (h(i,j,k)*tv%T(i,j,k) + h(i+1,j,k)*tv%T(i+1,j,k)) * i_2hlay
1246  s_lay = (h(i,j,k)*tv%S(i,j,k) + h(i+1,j,k)*tv%S(i+1,j,k)) * i_2hlay
1247  ghprime = g_h_rho0 * (dr_dt(i) * (t_lay*htot(i) - thtot(i)) + &
1248  dr_ds(i) * (s_lay*htot(i) - shtot(i)))
1249  else
1250  ghprime = g_h_rho0 * (gv%Rlay(k)*htot(i) - rhtot(i))
1251  endif
1252 
1253  if (ghprime > 0.0) then
1254  ribulk = cs%bulk_Ri_ML * exp(-htot(i) * idecay_len_tke(i))
1255  if (ribulk * uh2 <= (htot(i)**2) * ghprime) then
1256  visc%nkml_visc_u(i,j) = real(k_massive(i))
1257  do_i(i) = .false.
1258  elseif (ribulk * uh2 <= (htot(i) + hlay)**2 * ghprime) then
1259  visc%nkml_visc_u(i,j) = real(k-1) + &
1260  ( sqrt(ribulk * uh2 / ghprime) - htot(i) ) / hlay
1261  do_i(i) = .false.
1262  endif
1263  endif
1264  k_massive(i) = k
1265  endif ! hlay > h_tiny
1266 
1267  if (do_i(i)) do_any = .true.
1268  endif ; enddo
1269 
1270  if (.not.do_any) exit ! All columns are done.
1271  endif
1272 
1273  do i=isq,ieq ; if (do_i(i)) then
1274  htot(i) = htot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k))
1275  uhtot(i) = uhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * u(i,j,k)
1276  vhtot(i) = vhtot(i) + 0.25 * (h(i,j,k) * (v(i,j,k) + v(i,j-1,k)) + &
1277  h(i+1,j,k) * (v(i+1,j,k) + v(i+1,j-1,k)))
1278  if (use_eos) then
1279  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))
1280  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))
1281  else
1282  rhtot(i) = rhtot(i) + 0.5 * (h(i,j,k) + h(i+1,j,k)) * gv%Rlay(k)
1283  endif
1284  endif ; enddo
1285  enddo ; endif
1286 
1287  if (do_any) then ; do i=isq,ieq ; if (do_i(i)) then
1288  visc%nkml_visc_u(i,j) = k_massive(i)
1289  endif ; enddo ; endif
1290  endif ! dynamic_viscous_ML
1291 
1292  do_any_shelf = .false.
1293  if (associated(forces%frac_shelf_u)) then
1294  do i=isq,ieq
1295  if (forces%frac_shelf_u(i,j)*g%mask2dCu(i,j) == 0.0) then
1296  do_i(i) = .false.
1297  visc%tbl_thick_shelf_u(i,j) = 0.0 ; visc%kv_tbl_shelf_u(i,j) = 0.0
1298  else
1299  do_i(i) = .true. ; do_any_shelf = .true.
1300  endif
1301  enddo
1302  endif
1303 
1304  if (do_any_shelf) then
1305  do k=1,nz ; do i=isq,ieq ; if (do_i(i)) then
1306  if (u(i,j,k) *(h(i+1,j,k) - h(i,j,k)) >= 0) then
1307  h_at_vel(i,k) = 2.0*h(i,j,k)*h(i+1,j,k) / &
1308  (h(i,j,k) + h(i+1,j,k) + h_neglect)
1309  else
1310  h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i+1,j,k))
1311  endif
1312  else
1313  h_at_vel(i,k) = 0.0 ; ustar(i) = 0.0
1314  endif ; enddo ; enddo
1315 
1316  do i=isq,ieq ; if (do_i(i)) then
1317  htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
1318  thtot(i) = 0.0 ; shtot(i) = 0.0
1319  if (use_eos .or. .not.cs%linear_drag) then ; do k=1,nz
1320  if (htot_vel>=cs%Htbl_shelf) exit ! terminate the k loop
1321  hweight = min(cs%Htbl_shelf - htot_vel, h_at_vel(i,k))
1322  if (hweight <= 1.5*gv%Angstrom_H + h_neglect) cycle
1323 
1324  htot_vel = htot_vel + h_at_vel(i,k)
1325  hwtot = hwtot + hweight
1326 
1327  if (.not.cs%linear_drag) then
1328  v_at_u = set_v_at_u(v, h, g, i, j, k, mask_v, obc)
1329  hutot = hutot + hweight * sqrt(u(i,j,k)**2 + &
1330  v_at_u**2 + u_bg_sq)
1331  endif
1332  if (use_eos) then
1333  thtot(i) = thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
1334  shtot(i) = shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
1335  endif
1336  enddo ; endif
1337 
1338  if ((.not.cs%linear_drag) .and. (hwtot > 0.0)) then
1339  ustar(i) = cdrag_sqrt_z*us%T_to_s*hutot/hwtot
1340  else
1341  ustar(i) = cdrag_sqrt_z*us%T_to_s*cs%drag_bg_vel
1342  endif
1343 
1344  if (use_eos) then ; if (hwtot > 0.0) then
1345  t_eos(i) = thtot(i)/hwtot ; s_eos(i) = shtot(i)/hwtot
1346  else
1347  t_eos(i) = 0.0 ; s_eos(i) = 0.0
1348  endif ; endif
1349  endif ; enddo ! I-loop
1350 
1351  if (use_eos) then
1352  call calculate_density_derivs(t_eos, s_eos, forces%p_surf(:,j), &
1353  dr_dt, dr_ds, isq-g%IsdB+1, ieq-isq+1, tv%eqn_of_state)
1354  endif
1355 
1356  do i=isq,ieq ; if (do_i(i)) then
1357  ! The 400.0 in this expression is the square of a constant proposed
1358  ! by Killworth and Edwards, 1999, in equation (2.20).
1359  ustarsq = rho0x400_g * ustar(i)**2
1360  htot(i) = 0.0
1361  if (use_eos) then
1362  thtot(i) = 0.0 ; shtot(i) = 0.0
1363  do k=1,nz-1
1364  if (h_at_vel(i,k) <= 0.0) cycle
1365  t_lay = 0.5 * (tv%T(i,j,k) + tv%T(i+1,j,k))
1366  s_lay = 0.5 * (tv%S(i,j,k) + tv%S(i+1,j,k))
1367  oldfn = dr_dt(i)*(t_lay*htot(i) - thtot(i)) + dr_ds(i)*(s_lay*htot(i) - shtot(i))
1368  if (oldfn >= ustarsq) exit
1369 
1370  dfn = (dr_dt(i)*(0.5*(tv%T(i,j,k+1)+tv%T(i+1,j,k+1)) - t_lay) + &
1371  dr_ds(i)*(0.5*(tv%S(i,j,k+1)+tv%S(i+1,j,k+1)) - s_lay)) * &
1372  (h_at_vel(i,k)+htot(i))
1373  if ((oldfn + dfn) <= ustarsq) then
1374  dh = h_at_vel(i,k)
1375  else
1376  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
1377  endif
1378 
1379  htot(i) = htot(i) + dh
1380  thtot(i) = thtot(i) + t_lay*dh ; shtot(i) = shtot(i) + s_lay*dh
1381  enddo
1382  if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0)) then
1383  t_lay = 0.5*(tv%T(i,j,nz) + tv%T(i+1,j,nz))
1384  s_lay = 0.5*(tv%S(i,j,nz) + tv%S(i+1,j,nz))
1385  if (dr_dt(i)*(t_lay*htot(i) - thtot(i)) + &
1386  dr_ds(i)*(s_lay*htot(i) - shtot(i)) < ustarsq) &
1387  htot(i) = htot(i) + h_at_vel(i,nz)
1388  endif ! Examination of layer nz.
1389  else ! Use Rlay as the density variable.
1390  rhtot = 0.0
1391  do k=1,nz-1
1392  rlay = gv%Rlay(k) ; rlb = gv%Rlay(k+1)
1393 
1394  oldfn = rlay*htot(i) - rhtot(i)
1395  if (oldfn >= ustarsq) exit
1396 
1397  dfn = (rlb - rlay)*(h_at_vel(i,k)+htot(i))
1398  if ((oldfn + dfn) <= ustarsq) then
1399  dh = h_at_vel(i,k)
1400  else
1401  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
1402  endif
1403 
1404  htot(i) = htot(i) + dh
1405  rhtot(i) = rhtot(i) + rlay*dh
1406  enddo
1407  if (gv%Rlay(nz)*htot(i) - rhtot(i) < ustarsq) &
1408  htot(i) = htot(i) + h_at_vel(i,nz)
1409  endif ! use_EOS
1410 
1411  !visc%tbl_thick_shelf_u(I,j) = GV%H_to_Z * max(CS%Htbl_shelf_min, &
1412  ! htot(I) / (0.5 + sqrt(0.25 + &
1413  ! (htot(i)*(G%CoriolisBu(I,J-1)+G%CoriolisBu(I,J)))**2 / &
1414  ! (ustar(i)*GV%Z_to_H)**2 )) )
1415  ustar1 = ustar(i)*gv%Z_to_H
1416  h2f2 = (htot(i)*(g%CoriolisBu(i,j-1)+g%CoriolisBu(i,j)) + h_neglect*cs%omega)**2
1417  tbl_thick_z = gv%H_to_Z * max(cs%Htbl_shelf_min, &
1418  ( htot(i)*ustar1 ) / ( 0.5*ustar1 + sqrt((0.5*ustar1)**2 + h2f2 ) ) )
1419  visc%tbl_thick_shelf_u(i,j) = tbl_thick_z
1420  visc%Kv_tbl_shelf_u(i,j) = max(cs%Kv_TBL_min, cdrag_sqrt*ustar(i)*tbl_thick_z)
1421  endif ; enddo ! I-loop
1422  endif ! do_any_shelf
1423 
1424  enddo ! j-loop at u-points
1425 
1426  !$OMP parallel do default(private) shared(u,v,h,tv,forces,visc,dt,G,GV,US,CS,use_EOS,dt_Rho0, &
1427  !$OMP h_neglect,h_tiny,g_H_Rho0,is,ie,OBC,Jsq,Jeq,nz, &
1428  !$OMP U_bg_sq,cdrag_sqrt,cdrag_sqrt_Z,Rho0x400_G,nkml,mask_u)
1429  do j=jsq,jeq ! v-point loop
1430  if (cs%dynamic_viscous_ML) then
1431  do_any = .false.
1432  do i=is,ie
1433  htot(i) = 0.0
1434  if (g%mask2dCv(i,j) < 0.5) then
1435  do_i(i) = .false. ; visc%nkml_visc_v(i,j) = nkml
1436  else
1437  do_i(i) = .true. ; do_any = .true.
1438  k_massive(i) = nkml
1439  thtot(i) = 0.0 ; shtot(i) = 0.0 ; rhtot(i) = 0.0
1440  vhtot(i) = dt_rho0 * forces%tauy(i,j)
1441  uhtot(i) = 0.25 * dt_rho0 * ((forces%taux(i,j) + forces%taux(i-1,j+1)) + &
1442  (forces%taux(i-1,j) + forces%taux(i,j+1)))
1443 
1444  if (cs%omega_frac >= 1.0) then ; absf = 2.0*cs%omega ; else
1445  absf = 0.5*(abs(g%CoriolisBu(i-1,j)) + abs(g%CoriolisBu(i,j)))
1446  if (cs%omega_frac > 0.0) &
1447  absf = sqrt(cs%omega_frac*4.0*cs%omega**2 + (1.0-cs%omega_frac)*absf**2)
1448  endif
1449 
1450  u_star = max(cs%ustar_min, 0.5 * (forces%ustar(i,j) + forces%ustar(i,j+1)))
1451  idecay_len_tke(i) = ((absf / u_star) * cs%TKE_decay) * gv%H_to_Z
1452 
1453  endif
1454  enddo
1455 
1456  if (do_any) then ; do k=1,nz
1457  if (k > nkml) then
1458  do_any = .false.
1459  if (use_eos .and. (k==nkml+1)) then
1460  ! Find dRho/dT and dRho_dS.
1461  do i=is,ie
1462  press(i) = gv%H_to_Pa * htot(i)
1463  k2 = max(1,nkml)
1464  i_2hlay = 1.0 / (h(i,j,k2) + h(i,j+1,k2) + h_neglect)
1465  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
1466  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
1467  enddo
1468  call calculate_density_derivs(t_eos, s_eos, press, dr_dt, dr_ds, &
1469  is-g%IsdB+1, ie-is+1, tv%eqn_of_state)
1470  endif
1471 
1472  do i=is,ie ; if (do_i(i)) then
1473 
1474  hlay = 0.5*(h(i,j,k) + h(i,j+1,k))
1475  if (hlay > h_tiny) then ! Only consider non-vanished layers.
1476  i_2hlay = 1.0 / (h(i,j,k) + h(i,j+1,k))
1477  u_at_v = 0.5 * (h(i,j,k) * (u(i-1,j,k) + u(i,j,k)) + &
1478  h(i,j+1,k) * (u(i-1,j+1,k) + u(i,j+1,k))) * i_2hlay
1479  uh2 = us%m_s_to_L_T**2*((uhtot(i) - htot(i)*u_at_v)**2 + (vhtot(i) - htot(i)*v(i,j,k))**2)
1480 
1481  if (use_eos) then
1482  t_lay = (h(i,j,k)*tv%T(i,j,k) + h(i,j+1,k)*tv%T(i,j+1,k)) * i_2hlay
1483  s_lay = (h(i,j,k)*tv%S(i,j,k) + h(i,j+1,k)*tv%S(i,j+1,k)) * i_2hlay
1484  ghprime = g_h_rho0 * (dr_dt(i) * (t_lay*htot(i) - thtot(i)) + &
1485  dr_ds(i) * (s_lay*htot(i) - shtot(i)))
1486  else
1487  ghprime = g_h_rho0 * (gv%Rlay(k)*htot(i) - rhtot(i))
1488  endif
1489 
1490  if (ghprime > 0.0) then
1491  ribulk = cs%bulk_Ri_ML * exp(-htot(i) * idecay_len_tke(i))
1492  if (ribulk * uh2 <= htot(i)**2 * ghprime) then
1493  visc%nkml_visc_v(i,j) = real(k_massive(i))
1494  do_i(i) = .false.
1495  elseif (ribulk * uh2 <= (htot(i) + hlay)**2 * ghprime) then
1496  visc%nkml_visc_v(i,j) = real(k-1) + &
1497  ( sqrt(ribulk * uh2 / ghprime) - htot(i) ) / hlay
1498  do_i(i) = .false.
1499  endif
1500  endif
1501  k_massive(i) = k
1502  endif ! hlay > h_tiny
1503 
1504  if (do_i(i)) do_any = .true.
1505  endif ; enddo
1506 
1507  if (.not.do_any) exit ! All columns are done.
1508  endif
1509 
1510  do i=is,ie ; if (do_i(i)) then
1511  htot(i) = htot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k))
1512  vhtot(i) = vhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * v(i,j,k)
1513  uhtot(i) = uhtot(i) + 0.25 * (h(i,j,k) * (u(i-1,j,k) + u(i,j,k)) + &
1514  h(i,j+1,k) * (u(i-1,j+1,k) + u(i,j+1,k)))
1515  if (use_eos) then
1516  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))
1517  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))
1518  else
1519  rhtot(i) = rhtot(i) + 0.5 * (h(i,j,k) + h(i,j+1,k)) * gv%Rlay(k)
1520  endif
1521  endif ; enddo
1522  enddo ; endif
1523 
1524  if (do_any) then ; do i=is,ie ; if (do_i(i)) then
1525  visc%nkml_visc_v(i,j) = k_massive(i)
1526  endif ; enddo ; endif
1527  endif ! dynamic_viscous_ML
1528 
1529  do_any_shelf = .false.
1530  if (associated(forces%frac_shelf_v)) then
1531  do i=is,ie
1532  if (forces%frac_shelf_v(i,j)*g%mask2dCv(i,j) == 0.0) then
1533  do_i(i) = .false.
1534  visc%tbl_thick_shelf_v(i,j) = 0.0 ; visc%kv_tbl_shelf_v(i,j) = 0.0
1535  else
1536  do_i(i) = .true. ; do_any_shelf = .true.
1537  endif
1538  enddo
1539  endif
1540 
1541  if (do_any_shelf) then
1542  do k=1,nz ; do i=is,ie ; if (do_i(i)) then
1543  if (v(i,j,k) * (h(i,j+1,k) - h(i,j,k)) >= 0) then
1544  h_at_vel(i,k) = 2.0*h(i,j,k)*h(i,j+1,k) / &
1545  (h(i,j,k) + h(i,j+1,k) + h_neglect)
1546  else
1547  h_at_vel(i,k) = 0.5 * (h(i,j,k) + h(i,j+1,k))
1548  endif
1549  else
1550  h_at_vel(i,k) = 0.0 ; ustar(i) = 0.0
1551  endif ; enddo ; enddo
1552 
1553  do i=is,ie ; if (do_i(i)) then
1554  htot_vel = 0.0 ; hwtot = 0.0 ; hutot = 0.0
1555  thtot(i) = 0.0 ; shtot(i) = 0.0
1556  if (use_eos .or. .not.cs%linear_drag) then ; do k=1,nz
1557  if (htot_vel>=cs%Htbl_shelf) exit ! terminate the k loop
1558  hweight = min(cs%Htbl_shelf - htot_vel, h_at_vel(i,k))
1559  if (hweight <= 1.5*gv%Angstrom_H + h_neglect) cycle
1560 
1561  htot_vel = htot_vel + h_at_vel(i,k)
1562  hwtot = hwtot + hweight
1563 
1564  if (.not.cs%linear_drag) then
1565  u_at_v = set_u_at_v(u, h, g, i, j, k, mask_u, obc)
1566  hutot = hutot + hweight * sqrt(v(i,j,k)**2 + &
1567  u_at_v**2 + u_bg_sq)
1568  endif
1569  if (use_eos) then
1570  thtot(i) = thtot(i) + hweight * 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
1571  shtot(i) = shtot(i) + hweight * 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
1572  endif
1573  enddo ; endif
1574 
1575  if (.not.cs%linear_drag) then ; if (hwtot > 0.0) then
1576  ustar(i) = cdrag_sqrt_z*us%T_to_s*hutot/hwtot
1577  else
1578  ustar(i) = cdrag_sqrt_z*us%T_to_s*cs%drag_bg_vel
1579  endif ; endif
1580 
1581  if (use_eos) then ; if (hwtot > 0.0) then
1582  t_eos(i) = thtot(i)/hwtot ; s_eos(i) = shtot(i)/hwtot
1583  else
1584  t_eos(i) = 0.0 ; s_eos(i) = 0.0
1585  endif ; endif
1586  endif ; enddo ! I-loop
1587 
1588  if (use_eos) then
1589  call calculate_density_derivs(t_eos, s_eos, forces%p_surf(:,j), &
1590  dr_dt, dr_ds, is-g%IsdB+1, ie-is+1, tv%eqn_of_state)
1591  endif
1592 
1593  do i=is,ie ; if (do_i(i)) then
1594  ! The 400.0 in this expression is the square of a constant proposed
1595  ! by Killworth and Edwards, 1999, in equation (2.20).
1596  ustarsq = rho0x400_g * ustar(i)**2
1597  htot(i) = 0.0
1598  if (use_eos) then
1599  thtot(i) = 0.0 ; shtot(i) = 0.0
1600  do k=1,nz-1
1601  if (h_at_vel(i,k) <= 0.0) cycle
1602  t_lay = 0.5 * (tv%T(i,j,k) + tv%T(i,j+1,k))
1603  s_lay = 0.5 * (tv%S(i,j,k) + tv%S(i,j+1,k))
1604  oldfn = dr_dt(i)*(t_lay*htot(i) - thtot(i)) + dr_ds(i)*(s_lay*htot(i) - shtot(i))
1605  if (oldfn >= ustarsq) exit
1606 
1607  dfn = (dr_dt(i)*(0.5*(tv%T(i,j,k+1)+tv%T(i,j+1,k+1)) - t_lay) + &
1608  dr_ds(i)*(0.5*(tv%S(i,j,k+1)+tv%S(i,j+1,k+1)) - s_lay)) * &
1609  (h_at_vel(i,k)+htot(i))
1610  if ((oldfn + dfn) <= ustarsq) then
1611  dh = h_at_vel(i,k)
1612  else
1613  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
1614  endif
1615 
1616  htot(i) = htot(i) + dh
1617  thtot(i) = thtot(i) + t_lay*dh ; shtot(i) = shtot(i) + s_lay*dh
1618  enddo
1619  if ((oldfn < ustarsq) .and. (h_at_vel(i,nz) > 0.0)) then
1620  t_lay = 0.5*(tv%T(i,j,nz) + tv%T(i,j+1,nz))
1621  s_lay = 0.5*(tv%S(i,j,nz) + tv%S(i,j+1,nz))
1622  if (dr_dt(i)*(t_lay*htot(i) - thtot(i)) + &
1623  dr_ds(i)*(s_lay*htot(i) - shtot(i)) < ustarsq) &
1624  htot(i) = htot(i) + h_at_vel(i,nz)
1625  endif ! Examination of layer nz.
1626  else ! Use Rlay as the density variable.
1627  rhtot = 0.0
1628  do k=1,nz-1
1629  rlay = gv%Rlay(k) ; rlb = gv%Rlay(k+1)
1630 
1631  oldfn = rlay*htot(i) - rhtot(i)
1632  if (oldfn >= ustarsq) exit
1633 
1634  dfn = (rlb - rlay)*(h_at_vel(i,k)+htot(i))
1635  if ((oldfn + dfn) <= ustarsq) then
1636  dh = h_at_vel(i,k)
1637  else
1638  dh = h_at_vel(i,k) * sqrt((ustarsq-oldfn)/dfn)
1639  endif
1640 
1641  htot(i) = htot(i) + dh
1642  rhtot = rhtot + rlay*dh
1643  enddo
1644  if (gv%Rlay(nz)*htot(i) - rhtot(i) < ustarsq) &
1645  htot(i) = htot(i) + h_at_vel(i,nz)
1646  endif ! use_EOS
1647 
1648  !visc%tbl_thick_shelf_v(i,J) = GV%H_to_Z * max(CS%Htbl_shelf_min, &
1649  ! htot(i) / (0.5 + sqrt(0.25 + &
1650  ! (htot(i)*(G%CoriolisBu(I-1,J)+G%CoriolisBu(I,J)))**2 / &
1651  ! (ustar(i)*GV%Z_to_H)**2 )) )
1652  ustar1 = ustar(i)*gv%Z_to_H
1653  h2f2 = (htot(i)*(g%CoriolisBu(i-1,j)+g%CoriolisBu(i,j)) + h_neglect*cs%omega)**2
1654  tbl_thick_z = gv%H_to_Z * max(cs%Htbl_shelf_min, &
1655  ( htot(i)*ustar1 ) / ( 0.5*ustar1 + sqrt((0.5*ustar1)**2 + h2f2 ) ) )
1656  visc%tbl_thick_shelf_v(i,j) = tbl_thick_z
1657  visc%Kv_tbl_shelf_v(i,j) = max(cs%Kv_TBL_min, cdrag_sqrt*ustar(i)*tbl_thick_z)
1658 
1659  endif ; enddo ! i-loop
1660  endif ! do_any_shelf
1661 
1662  enddo ! J-loop at v-points
1663 
1664  if (cs%debug) then
1665  if (associated(visc%nkml_visc_u) .and. associated(visc%nkml_visc_v)) &
1666  call uvchksum("nkml_visc_[uv]", visc%nkml_visc_u, &
1667  visc%nkml_visc_v, g%HI,haloshift=0)
1668  endif
1669  if (cs%id_nkml_visc_u > 0) &
1670  call post_data(cs%id_nkml_visc_u, visc%nkml_visc_u, cs%diag)
1671  if (cs%id_nkml_visc_v > 0) &
1672  call post_data(cs%id_nkml_visc_v, visc%nkml_visc_v, cs%diag)
1673