PUBLIC INTERFACE ~ PUBLIC DATA ~ PUBLIC ROUTINES ~ NAMELIST ~ DIAGNOSTIC FIELDS ~ ERROR MESSAGES ~ REFERENCES ~ NOTES

Module ocean_lap_friction_mod

Contact:  S. M. Griffies
Reviewers: 
Change History: WebCVS Log


OVERVIEW

This module computes the thickness weighted time tendency for horizontal velocity arising from horizontal Laplacian friction.

This module computes the thickness weighted time tendency for horizontal velocity arising from horizontal Laplacian friction. The viscosity used to determine the strength of the tendency can be a general function of space and time as specified by the Smagorinsky approach as well as a grid-scale dependent background viscosity. The form of the friction operator can be isotropic or anisotropic in the horizontal plane.


OTHER MODULES USED

      constants_mod
diag_manager_mod
fms_mod
mpp_domains_mod
mpp_mod
ocean_domains_mod
ocean_operators_mod
ocean_types_mod

PUBLIC INTERFACE

ocean_lap_friction_init:
horz_lap_friction:
BDX_EU_smag:
BDY_NU_smag:
anisotropic_ncar:
horz_lap_viscosity_check:
horz_lap_reynolds_check:
compute_neptune_velocity:


PUBLIC DATA

None.


PUBLIC ROUTINES

  1. ocean_lap_friction_init

    DESCRIPTION
    Initialize the horizontal Laplacian friction module by registering fields for diagnostic output and performing some numerical checks to see that viscosity is set appropriately.


  2. horz_lap_friction

    DESCRIPTION
    This function computes thickness weighted time tendency for horizontal velocity (i.e., thickness weighted acceleration) from horizontal Laplacian friction.

    The algorithm is derived from a functional approach that ensures kinetic energy is consistenty dissipated for all flow configurations. The triad do-loops are expanded in order to enhance the ability of cache-based machines to keep most of the variables on-cache.

    Fundamental to the scheme are the rates of horizontal deformation
    horizontal tension = DT = (dy)(u/dy)_x - (dx)(v/dx)_y
    horizontal strain = DS = (dx)(u/dx)_y + (dy)(v/dy)_x
    Units of the tension and strain are sec^-1.

    Four tensions and four strains are computed for each velocity point,
    corresponding to the four triads surrounding the point.
    The following notation is used to distinguish the triads:
    (0,1)=northwest triad (1,1)=northeast triad,
    (0,0)=southwest triad, (1,0)=southeast triad

    A triad contributes when at least one of its velocities is not a land point. In order to obtain the correct tension and strain next to boundaries, tension and strain should not be masked with umask.



  3. BDX_EU_smag

    DESCRIPTION
    Compute backwards Derivative in X of a quantity defined on the east face of a U-cell. Slightly modified version of BDX_EU used in ocean_operators.F90. If input is a(i,j) then output is defined at (i-1/2,j) BDX_EU_smag changes dimensions by m^-3


    INPUT
    a    field defined on the east face of a U-cell
       [real, dimension(isd:ied,jsd:jed)]

  4. BDY_NU_smag

    DESCRIPTION
    Compute backwards Derivative in Y of a quantity defined on the north face of a U-cell. Slightly modified version of BDY_EU used in ocean_operators.F90. If input is a(i,j) then output is defined at (i,j-1/2) BDY_EU_smag changes dimensions by m^-3


    INPUT
    a    field defined on the north face of a U-cell
       [real, dimension(isd:ied,jsd:jed)]

  5. anisotropic_ncar

    DESCRIPTION
    Spatially-varying anisotropic viscosity initialization

    This routine defines NCOM-like spatial distributions of viscosity coefficients F_PARA and F_PERP. Uses NCAR CCSM2.0 algorithm with cm^2/sec --> m^2/sec.

    written by: Stephen Yeager 3/2000
    modified by: Gokhan Danabasoglu (08/2001)
    port to mom4: Stephen.Griffies@noaa.gov (9/2002)

    "A_viscosity" = F_PARA = Along = viscosity parallel to flow = max{0.5*visc_vel_scale(z)*A*max[dx,dy],vconst_6}

    where
    A = 0.425 * cos(pi*y*radian/30) + 0.575 for |y*radian| < 30
    A = 0.15 otherwise

    Here, A provides a horizontal variation for visc_vel_scale.

    "B_viscosity" = F_PERP = Across = viscosity perpendicular to flow = max( bu, bv)

    and
    F_PARA = min(F_PARA, AMAX_CFL),
    F_PERP = min(F_PERP, AMAX_CFL),
    F_PARA = max(F_PARA, F_PERP)
    are enforced

    In the above equations,

    bu = vconst_1 * ( 1 + vconst_2 * ( 1 + cos( 2*y + pi ) ) )
    bv = vconst_3 * beta_f * dx^3 * exp( - (vconst_4 * distance)^2 )
    with
    beta_f (x,y) = 2 * omega * cos(ULAT(i,j)) / radius
    distance (x,y,z) = actual distance to "vconst_5" points
    west of the nearest western boundary
    dx (x,y) = DXU(i,j)
    dy (x,y) = DYU(i,j)
    visc_vel_scale (z) = vconst_7 * exp(-zt(k)/visc_vel_scale_length)
    visc_vel_scale_length = e-folding scale ( = 1500.0e2 cm)
    y (x,y) = ULAT(i,j), latitude of "u/v" grid pts in radians
    In mom4, ULAT(radians) = xu*pi/180 with xu(i,j) the longitude of U grid points in degrees

    "vconst_#" are input parameters defined in namelist ocean_lap_friction_general_nml. "vconst_1", "vconst_6", and "vconst_4" have dimensions of cm^2/s, cm^2/s, and 1/cm, respectively. "vconst_5" is an INTEGER.

    NOTE: The nearest western boundary computations are done along the model i-grid lines. Therefore, viscosity based on these are only approximate in the high Northern Hemisphere when using generalized coordinates with coordinate pole(s) shifted onto land.



  6. horz_lap_viscosity_check

    DESCRIPTION
    Subroutine to perform linear stability check for the Laplacian operator given a value for the horizontal biharmonic viscosity.


  7. horz_lap_reynolds_check

    DESCRIPTION
    Subroutine to compute the LLaplacian grid Reynolds number. Large Reynolds numbers indicate regions where solution may experience some grid noise due to lack of enough horizontal friction.


    INPUT
    u    Horizontal velocity field at time tau
       [real, dimension(isd:ied,jsd:jed,nk,2)]

  8. compute_neptune_velocity

    DESCRIPTION
    Compute Neptune velocity.

    Method follows that used in MOM2/3 as implemented by Greg Holloway (zounds@ios.bc.ca) and Michael Eby (eby@uvic.ca)

    Neptune is calculated as an equilibrium streamfunction given by pnep = -f*snep*snep*ht and is applied through friction

    ht = depth of tracer cells snep = spnep + (senep-spnep)*(0.5 + 0.5*cos(2.0*latitude))

    Neptune length scale snep has a value of senep at the equator and smoothly changes to spnep at the poles

    Reference: Holloway, G., 1992: Representing topographic stress for large scale ocean models, J. Phys. Oceanogr., 22, 1033-1046




NAMELIST

&ocean_lap_friction_general_nml

lap_friction_on
Must be true to use this module.
[logical]
k_smag_iso
This is the dimensionless Smagorinsky coefficient used to set the scale of the Smagorinsky isotropic viscosity.
[real, units: dimensionless]
k_smag_aniso
This is the dimensionless Smagorinsky coefficient used to set the scale of the Smagorinsky anisotropic viscosity.
[real, units: dimensionless]
viscosity_ncar
Anisotropic background viscosities used by NCAR.
[logical]
vel_micom_iso
Velocity scale that is used for computing the MICOM isotropic viscosity.
[real, units: m/sec]
vel_micom_aniso
Velocity scale that is used for computing the MICOM anisotropic viscosity.
[real, units: m/sec]
equatorial_zonal
Orient the anisotropic friction within a latitudinal band according to zonal direction.
[logical]
equatorial_zonal_lat
Latitudinal band to use the zonal friction orientation.
[real]
ncar_isotropic_off_equator
Polewards of equatorial_zonal_lat, revert NCAR scheme to isotropic
[logical]
equatorial_no_smag
Turn smag off within equatorial_zonal_lat region.
[logical]
eq_vel_micom_iso
Velocity scale that is used for computing the MICOM isotropic viscosity within a user specified equatorial band.
[real]
eq_vel_micom_aniso
Velocity scale that is used for computing the MICOM anisotropic viscosity within a user specified equatorial band.
[real]
eq_lat_micom
Equatorial latitude band (degrees) within which the MICOM viscosity is set according to eq_vel_micom_iso and eq_vel_micom_aniso.
[real]
restrict_polar_visc
For restricting the background viscosity poleward of a latitude. This method may be useful for coupling to an ice model in which case the horizontal viscosity may need to be a bit smaller to maintain time step constraints. This is because the effective friction is larger than that just within the ocean.
[logical]
restrict_polar_visc_lat
Latitude poleward of which we restrict the viscosity.
[real]
restrict_polar_visc_ratio
Ratio of the normal critical value that we limit the viscosity to be no greater than. If restrict_polar_visc_ratio=1.0 then there is no special limitation of the viscosity beyond that of the one-dimensional stability constraint.
[real]
bottom_5point
To alleviate problems with small partial cells, it is often necessary to reduce the operator to the traditional 5-point Laplacian at the ocean bottom. This logical implements this mixing.
[logical]
neptune
Set to true for computing friction relative to Neptune barotropic velocity.
[logical]
neptune_length_eq
Length scale used to compute Neptune velocity at equator.
[real, units: m]
neptune_length_pole
Length scale used to compute Neptune velocity at pole.
[real, units: m]
vconst_1
Background viscosity for NCAR algorithm.
[real, units: cm^2/sec]
vconst_2
For NCAR viscosity algorithm.
[real]
vconst_3
For NCAR viscosity algorithm.
[real]
vconst_4
For NCAR viscosity algorithm.
[real, units: 1/cm]
vconst_5
For NCAR viscosity algorithm.
[integer]
vconst_6
For NCAR viscosity algorithm.
[real, units: cm^2/sec]
vconst_7
For NCAR viscosity algorithm.
[, units: cm/sec]
debug_ncar_A
Sets f_perp=f_para for debugging purposes with the NCAR scheme.
[logical]
debug_ncar_B
Sets f_para=f_perp for debugging purposes with the NCAR scheme.
[logical]


DATA SETS

None.


ERROR MESSAGES

None.


REFERENCES

  1. S.M. Griffies and R.W. Hallberg, Biharmonic friction with a Smagorinsky viscosity for use in large-scale eddy-permitting ocean models Monthly Weather Review vol 128 (2000) pages 2935--2946
  2. R. D. Smith and J. C. McWilliams, Anisotropic viscosity for ocean models Ocean Modelling submitted 2002


COMPILER SPECIFICS

None.


PRECOMPILER OPTIONS

None.


LOADER OPTIONS

None.


TEST PROGRAM

None.


KNOWN BUGS

None.


NOTES

The ocean model can generally run with both Laplacian and biharmonic friction enabled at the same time. Such has been found useful for some eddying ocean simulations.


FUTURE PLANS

None.


top