**INTERFACE:**

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

(119) |

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

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 IMPLICIT NONE

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