**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) |

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

use internal_pressure use variables_3d, only: 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)