/*
**  depth.c	1.1	Solaris 2.3 Unix	940906	SIO/ODF	fmd
**
**  Functions for calculating depth:
**
**  void Depth(int n, double g, double *p, int pInc, double *rho,
**		int rInc, double *depth, int dInc) 
**	calculate depth by integration of insitu density.
**  double PtoZ(double grav, double dynz, double prsdb)
**	calculate depth from pressure using Saunders-Mantyla method.
**  double ZtoP(double grav, double dynz, double depth)
**	calculate pressure from depth using Saunders-Mantyla method.
**
**  Plus, soon-to-be obsolete Fortran interface functions.
**
**  References:
**  =========
**
**	Sverdrup, H. U.,Johnson, M. W., and Fleming, R. H., 1942.
**	The Oceans, Their Physics, Chemistry and General Biology.
**	Prentice-Hall, Inc., Englewood Cliff, N.J.
**
**	Saunders, P. M., 1981. Practical Conversion of Pressure to Depth.
**	Journal of Physical Oceanography 11, 573-574.
**	Mantyla, A. W., 1982-1983. Private correspondence.
*/
#include <math.h>
#include <odflib.h>

#ifndef	FORTRAN

/*
**  Calculate depth by integration of insitu density.
*/
void
Depth(n, g, p, pInc, rho, rInc, depth, dInc)
	int	n;	/* -> number of intervals in p, rho and depth */
	double	g;	/* -> local gravity (m/sec**2) @ 0.0 db. */
	double	*p;	/* -> pressure series (decibars) */
	int	pInc;	/* -> byte increment for pressure series */
	double	*rho;	/* -> insitu density series (kg/m**3) */
	int	rInc;	/* -> byte increment for insitu density series */
	double	*depth;	/* <- depth series (meters) */
	int	dInc;	/* -> byte increment for depth series */
{
	int	i;
	double	*pP, *nP, *pR, *nR, *pD, *nD;
    /*
    **  When calling Depth() repeatedly with a two-element
    **  series, the first call should be with a one-element series to
    **  initialize the starting value (see depth_(), below).
    */
	if (n!=2) {			/* initialize the series */
	/*
	**  If the integration starts from > 15 db, calculate
	**  depth relative to starting place. Otherwise,
	**  calculate from surface.
	*/
	    if (p[0] > 15.0)
		depth[0] = 0.0;
	    else
		depth[0] = p[0]/(rho[0]*10000.0*(g+2.184e-6*p[0]));
	}
    /*
    **  Calculate the rest of the series.
    */
	pP = p; pR = rho; pD = depth;
	nP = (double *)((char *)pP + pInc);
	nR = (double *)((char *)pR + rInc);
	nD = (double *)((char *)pD + dInc);
	for (i=1; i<n;
	    i++,
	    pP = nP, pR = nR, pD = nD,
	    nP = (double *)((char *)pP + pInc),
	    nR = (double *)((char *)pR + rInc),
	    nD = (double *)((char *)pD + dInc))
	    nD[0] = pD[0]+(nP[0]-pP[0])/
		((nR[0]+pR[0])*5000.0*(g+2.184e-6*nP[0]));
	     /* depth in meters */

} /* Depth() */

/*
**  Calculate depth from pressure using Saunders-Mantyla method.
*/
double				/* <- depth (meters) */
PtoZ(grav, dynz, prsdb)
	double grav;		/* -> local gravity at 0db pressure */
	double dynz;		/* -> dynamic height (dynamic meters) */
	double prsdb;		/* -> pressure (decibars) */
{
/*
**  Method:
**  ======
**
**       a standard ocean (salinity=35.0,t=0.0) is assumed. A dynamic
**       height correction to account for deviations from the standard
**       ocean is then applied. depth is calculated from pressure by:
**
**       z = p/(b*c) + dynz/b
**
**       where:
**
**          p    = insitu pressure (decibars)
**          b    = insitu gravity as a function of latitude and pressure
**                 (decameters/sec/sec)
**          c    = insitu mean density rho(35.0,0.0,p)
**                 (grams/centimeter**3)
**          dynz = dynamic height in dynamic meters
**          z    = depth in meters
*/
	double	a, b, c, depth;

	a     = prsdb * 2.184e-6;	/* pressure correction to gravity */
	b     = (grav + a * 0.5) * 0.1;	/* insitu gravity decameters/sec**2 */
	c     = a + 1.0285;		/* insitu density rho(35.0,0.0,p) */
	if (prsdb != 0.0)
	    depth = prsdb/(b*c);
	else
	    depth = 0.0;
	if (isnan(dynz) == 0)
	    depth += dynz/b;
	return (depth);
} /* PtoZ() */

/*
**  Calculate pressure from depth using Saunders-Mantyla method.
*/
double				/* <- pressure (decibars) */
ZtoP(grav, dynz, depth)
	double grav;		/* -> local gravity at 0db pressure */
	double dynz;		/* -> dynamic height (dynamic meters) */
	double depth;		/* -> depth (meters) */
{
	double	a, b, c, prsdb;

	grav *= 0.1;
	a     = depth * 2.42e-13;
	b     = depth * (grav * 2.184e-6 + 1.13135e-7) - 1.0;
	c     = grav * depth;
	if (isnan(dynz) == 0) {
	    b -= dynz * 2.2e-6;
	    c -= dynz;
	}
	c    *= 1.0285;
	if (a != 0.0)
	    prsdb = (-b - sqrt(b*b - 4.0*a*c))/(a+a);
	else
	    prsdb = 0.0;
	return (prsdb);
} /* ZtoP() */

#else	/* FORTRAN */

void
depth_(grav, prsdb, prvprs, rho, prvrho, prvdpt, depthm)
	float	*grav, *prsdb, *prvprs, *rho, *prvrho, *prvdpt, *depthm;
{
	double	g, p[2], drho[2], depth[2];
	int	inc = sizeof (double);
    /*
    **  depth_() is called once per depth value.
    **  The first call should be made with *prvdpt==quiet_nan() to initialize
    **  the surface (or first) value.
    */
	g = *grav - 2.184e-6*(*prsdb);  /* convert in-situ gravity to 0.0 db */
	if (isnan(*prvdpt)) {		/* initialize */
	    p[0]       = *prsdb;
	    drho[0]    = *rho;
	    Depth(1, g, p, inc, drho, inc, depth, inc);
	    *depthm    = depth[0];
	} else {
	    p[0]       = *prvprs; p[1]    = *prsdb;
	    drho[0]    = *prvrho; drho[1] = *rho;
	    depth[0]   = *prvdpt;
	    Depth(2, g, p, inc, drho, inc, depth, inc);
	    *depthm    = depth[1];
	}
}
void
pzcon_(cnvrsn, grav, dynz, prsdb, depth)
	int	*cnvrsn;
	float	*grav, *dynz, *prsdb, *depth;
{
	if (*cnvrsn==0)
	    *depth = PtoZ(*grav, *dynz, *prsdb);
	else
	    *prsdb = ZtoP(*grav, *dynz, *depth);
}

void
dpzcon_(cnvrsn, grav, dynz, prsdb, depth)
	int	*cnvrsn;
	double	*grav, *dynz, *prsdb, *depth;
{
	if (*cnvrsn==0)
	    *depth = PtoZ(*grav, *dynz, *prsdb);
	else
	    *prsdb = ZtoP(*grav, *dynz, *depth);
}
#endif	/* FORTRAN */