subroutine ip_shchepetkin_mcwilliams()

Here, the pressure gradient is calculated according to the method and the algorithm suggested by Shchepetkin and McWilliams, 2003. This method uses a nonconservative Density-Jacobian scheme, based on cubic polynomial fits for the bouyancy "buoy" and "zz", the vertical position of rho-points, as functions of its respective array indices. The cubic polynomials are monotonized by using harmonic mean instead of linear averages to interpolate slopes. Exact anti-symmetry of the density Jacobian

$\displaystyle J(rho,zz)=-J(zz,rho)$ (119)

is retained for the density/bouyancy Jacobian in the pressure gradient formulation in x-direction for a non aligned vertical coordinate $ \sigma $, the atmospheric pressure $ p_0$ and the sea surface elevation $ \eta$:

$\displaystyle - {1 \over \rho_0} \partial_x p = \underbrace{-{1 \over \rho_0} \...
...brace{buoy(\eta) \partial_x\eta + \int_z^\eta J(buoy,zz)\mbox{d}\sigma}_{idpdx}$ (120)

Details about the calculation of the integral over the Jacobian in (120) can be found in Shchepetkin and McWilliams, 2003.

If parameter OneFifth (below) is set to zero, the scheme should become identical to standard Jacobian. USES:

    use internal_pressure
    use variables_3d, only: hn,buoy,sseo
    use domain, only: H,az,au,av
  $ use omp_lib
    Original author(s): Richard Hofmeister
    integer                   :: i,j,k
    REALTYPE                  :: dR(I3DFIELD)
    REALTYPE                  :: dZ(I3DFIELD)
    REALTYPE                  :: P(I3DFIELD)
    REALTYPE                  :: dxm1,dym1,cff,cff1,cff2
    REALTYPE                  :: AJ
    REALTYPE                  :: eps=1.e-10
    REALTYPE                  :: OneFifth = 0.2
    REALTYPE                  :: FC(I2DFIELD)
    REALTYPE                  :: dZx(I2DFIELD)
    REALTYPE                  :: dRx(I2DFIELD)

kklingbe 2017-10-02