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 [L T-1 ~> m s-1] or other units.
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]iThe i-index of the u-location to work on.
[in]jThe j-index of the u-location to work on.
[in]kThe k-index of the u-location to work on.
[in]mask2dcuA multiplicative mask of the u-points
obcA pointer to an open boundary condition structure
Returns
The return value of u at v points in the same units as u, i.e. [L T-1 ~> m s-1] or other units.

Definition at line 971 of file MOM_set_viscosity.F90.

971  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
972  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
973  intent(in) :: u !< The zonal velocity [L T-1 ~> m s-1] or other units.
974  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
975  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
976  integer, intent(in) :: i !< The i-index of the u-location to work on.
977  integer, intent(in) :: j !< The j-index of the u-location to work on.
978  integer, intent(in) :: k !< The k-index of the u-location to work on.
979  real, dimension(SZIB_(G),SZJ_(G)), &
980  intent(in) :: mask2dCu !< A multiplicative mask of the u-points
981  type(ocean_OBC_type), pointer :: OBC !< A pointer to an open boundary condition structure
982  real :: set_u_at_v !< The return value of u at v points in the
983  !! same units as u, i.e. [L T-1 ~> m s-1] or other units.
984 
985  ! This subroutine finds a thickness-weighted value of u at the v-points.
986  real :: hwt(-1:0,0:1) ! Masked weights used to average u onto v [H ~> m or kg m-2].
987  real :: hwt_tot ! The sum of the masked thicknesses [H ~> m or kg m-2].
988  integer :: i0, j0, i1, j1
989 
990  do j0 = 0,1 ; do i0 = -1,0 ; i1 = i+i0 ; j1 = j+j0
991  hwt(i0,j0) = (h(i1,j1,k) + h(i1+1,j1,k)) * mask2dcu(i1,j1)
992  enddo ; enddo
993 
994  if (associated(obc)) then ; if (obc%number_of_segments > 0) then
995  do j0 = 0,1 ; do i0 = -1,0 ; if ((obc%segnum_u(i+i0,j+j0) /= obc_none)) then
996  i1 = i+i0 ; j1 = j+j0
997  if (obc%segment(obc%segnum_u(i1,j1))%direction == obc_direction_e) then
998  hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dcu(i1,j1)
999  elseif (obc%segment(obc%segnum_u(i1,j1))%direction == obc_direction_w) then
1000  hwt(i0,j0) = 2.0 * h(i1+1,j1,k) * mask2dcu(i1,j1)
1001  endif
1002  endif ; enddo ; enddo
1003  endif ; endif
1004 
1005  hwt_tot = (hwt(-1,0) + hwt(0,1)) + (hwt(0,0) + hwt(-1,1))
1006  set_u_at_v = 0.0
1007  if (hwt_tot > 0.0) set_u_at_v = &
1008  ((hwt(0,0) * u(i,j,k) + hwt(-1,1) * u(i-1,j+1,k)) + &
1009  (hwt(-1,0) * u(i-1,j,k) + hwt(0,1) * u(i,j+1,k))) / hwt_tot
1010 

Referenced by set_viscous_bbl(), and set_viscous_ml().

Here is the caller graph for this function:

◆ 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 [L T-1 ~> m s-1]
[in]hLayer thicknesses [H ~> m or kg m-2]
[in]iThe i-index of the u-location to work on.
[in]jThe j-index of the u-location to work on.
[in]kThe k-index of the u-location to work on.
[in]mask2dcvA multiplicative mask of the v-points
obcA pointer to an open boundary condition structure
Returns
The return value of v at u points points in the same units as u, i.e. [L T-1 ~> m s-1] or other units.

Definition at line 927 of file MOM_set_viscosity.F90.

927  type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
928  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
929  intent(in) :: v !< The meridional velocity [L T-1 ~> m s-1]
930  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
931  intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
932  integer, intent(in) :: i !< The i-index of the u-location to work on.
933  integer, intent(in) :: j !< The j-index of the u-location to work on.
934  integer, intent(in) :: k !< The k-index of the u-location to work on.
935  real, dimension(SZI_(G),SZJB_(G)),&
936  intent(in) :: mask2dCv !< A multiplicative mask of the v-points
937  type(ocean_OBC_type), pointer :: OBC !< A pointer to an open boundary condition structure
938  real :: set_v_at_u !< The return value of v at u points points in the
939  !! same units as u, i.e. [L T-1 ~> m s-1] or other units.
940 
941  ! This subroutine finds a thickness-weighted value of v at the u-points.
942  real :: hwt(0:1,-1:0) ! Masked weights used to average u onto v [H ~> m or kg m-2].
943  real :: hwt_tot ! The sum of the masked thicknesses [H ~> m or kg m-2].
944  integer :: i0, j0, i1, j1
945 
946  do j0 = -1,0 ; do i0 = 0,1 ; i1 = i+i0 ; j1 = j+j0
947  hwt(i0,j0) = (h(i1,j1,k) + h(i1,j1+1,k)) * mask2dcv(i1,j1)
948  enddo ; enddo
949 
950  if (associated(obc)) then ; if (obc%number_of_segments > 0) then
951  do j0 = -1,0 ; do i0 = 0,1 ; if ((obc%segnum_v(i+i0,j+j0) /= obc_none)) then
952  i1 = i+i0 ; j1 = j+j0
953  if (obc%segment(obc%segnum_v(i1,j1))%direction == obc_direction_n) then
954  hwt(i0,j0) = 2.0 * h(i1,j1,k) * mask2dcv(i1,j1)
955  elseif (obc%segment(obc%segnum_v(i1,j1))%direction == obc_direction_s) then
956  hwt(i0,j0) = 2.0 * h(i1,j1+1,k) * mask2dcv(i1,j1)
957  endif
958  endif ; enddo ; enddo
959  endif ; endif
960 
961  hwt_tot = (hwt(0,-1) + hwt(1,0)) + (hwt(1,-1) + hwt(0,0))
962  set_v_at_u = 0.0
963  if (hwt_tot > 0.0) set_v_at_u = &
964  ((hwt(0,0) * v(i,j,k) + hwt(1,-1) * v(i+1,j-1,k)) + &
965  (hwt(1,0) * v(i+1,j,k) + hwt(0,-1) * v(i,j-1,k))) / hwt_tot
966 

Referenced by set_viscous_bbl(), and set_viscous_ml().

Here is the caller graph for this function:

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

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

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

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

References mom_cvmix_ddiff::cvmix_ddiff_is_used(), mom_kappa_shear::kappa_shear_at_vertex(), mom_kappa_shear::kappa_shear_is_used(), mom_error_handler::mom_error(), and mom_restart::register_restart_field_as_obsolete().

Here is the call graph for this function:

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

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

References mom_cvmix_conv::cvmix_conv_is_used(), mom_cvmix_shear::cvmix_shear_is_used(), mom_kappa_shear::kappa_shear_at_vertex(), and mom_kappa_shear::kappa_shear_is_used().

Referenced by mom::initialize_mom().

Here is the call graph for this function:
Here is the caller graph for this function:

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

References mom_error_handler::mom_error(), set_u_at_v(), and set_v_at_u().

Here is the call graph for this function:

◆ 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 [L T-1 ~> m s-1].
[in]vThe meridional velocity [L T-1 ~> m s-1].
[in]hLayer thicknesses [H ~> m or kg m-2].
[in]tvA structure containing pointers to any available thermodynamic fields. Absent fields have NULL ptrs.
[in]forcesA structure with the driving mechanical forces
[in,out]viscA structure containing vertical viscosities and related fields.
[in]dtTime increment [T ~> s].
csThe control structure returned by a previous call to 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 1019 of file MOM_set_viscosity.F90.

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

References mom_error_handler::mom_error(), set_u_at_v(), and set_v_at_u().

Referenced by mom_dynamics_split_rk2::step_mom_dyn_split_rk2(), mom_dynamics_unsplit::step_mom_dyn_unsplit(), and mom_dynamics_unsplit_rk2::step_mom_dyn_unsplit_rk2().

Here is the call graph for this function:
Here is the caller graph for this function: