/*
**  dynheight.c	1.03	Solaris 2.3 Unix	940906	SIO/ODF	fmd
**
**  Calculate dynamic height.
**
**  Reference:
**  =========
**
**    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.
*/
#include <math.h>
#include <odflib.h>

#ifndef	FORTRAN

/*
**  Calculate dynamic height (dynamic meters).
*/
void
DynamicHeight(n, p, pInc, spVolAn, spVolAnInc, dynHt, dynHtInc)
	int	n;		/* -> n intervals in p, spVolAn and dynHt */
	double	*p;		/* -> pressure series (decibars) */
	int	pInc;		/* -> byte increment for p */
	double	*spVolAn;	/* -> sp volume anomaly series (m**3/kg) */
	int	spVolAnInc;	/* -> byte increment for spVolAn */
	double	*dynHt;		/* <- dynamic height series (dynamic meters)
				      (or joules/(10*kg) or m**2/(10*sec**2)) */
	int	dynHtInc;	/* -> byte increment for dynHt */
{
	int	i;
	double	*pP, *nP, *pSp, *nSp, *pD, *nD;
    /*
    **  When calling DynamicHeight() repeatedly with a two-element
    **  series, the first call should be with a one-element series to
    **  initialize the starting value (see dynht_(), below).
    */
	if (n!=2) {			/* initialize the series */
	    /*
	    **  If the integration starts from > 15 db, calculate
	    **  dynamic height relative to starting place. Otherwise,
	    **  calculate from surface.
	    */
		if (p[0] > 15.0)
		    dynHt[0] = 0.0;
		else
		    dynHt[0] = spVolAn[0]*1000.0*p[0];
	}
    /*
    **  Calculate the rest of the series.
    **  Dynamic height is calculated by integration of specific
    **  volume anomaly with respect to pressure.
    */
	pP  = p; pSp = spVolAn; pD = dynHt;
	nP  = (double *)((char *)pP+pInc);
	nSp = (double *)((char *)pSp+spVolAnInc);
	nD  = (double *)((char *)pD+dynHtInc);
	for (i=1; i<n;
	    i++,
	    pP = nP, pSp = nSp, pD = nD,
	    nP  = (double *)((char *)pP+pInc),
	    nSp = (double *)((char *)pSp+spVolAnInc),
	    nD  = (double *)((char *)pD+dynHtInc))
	    nD[0] = pD[0]+(pSp[0]+nSp[0])*500.0*(nP[0]-pP[0]);
	     /* dynamic height in dynamic meters */

} /* DynamicHeight() */

#else	/* FORTRAN */

void
dynht_(prsdb, prvprs, delta, prvdel, prvdht, dynht)
	float	*prsdb, *prvprs, *delta, *prvdel, *prvdht, *dynht;
{
	double	p[2], spVolAn[2], dynHt[2];
	int	inc = sizeof (double);

    /*
    **  The FORTRAN interface is called once per dynamic height value.
    **  The first call should be made with *prvdht==quiet_nan() to initialize
    **  the surface (or first) value.
    */
	if (isnan(*prvdht)) {		/* initialize */
	    p[0]       = *prsdb;
	    spVolAn[0] = *delta;
	    DynamicHeight(1, p, inc, spVolAn, inc, dynHt, inc);
	    *dynht     = dynHt[0];
	} else {
	    p[0]       = *prvprs; p[1]       = *prsdb;
	    spVolAn[0] = *prvdel; spVolAn[1] = *delta;
	    dynHt[0]   = *prvdht;
	    DynamicHeight(2, p, inc, spVolAn, inc, dynHt, inc);
	    *dynht     = dynHt[1];
	}
}

void
ddynht_(prsdb, prvprs, delta, prvdel, prvdht, dynht)
	double	*prsdb, *prvprs, *delta, *prvdel, *prvdht, *dynht;
{
	double	p[2], spVolAn[2], dynHt[2];
	int	inc = sizeof (double);

    /*
    **  The FORTRAN interface is called once per dynamic height value.
    **  The first call should be made with *prvdht==quiet_nan() to initialize
    **  the surface (or first) value.
    */
	if (isnan(*prvdht)) {		/* initialize */
	    p[0]       = *prsdb;
	    spVolAn[0] = *delta;
	    DynamicHeight(1, p, inc, spVolAn, inc, dynHt, inc);
	    *dynht     = dynHt[0];
	} else {
	    p[0]       = *prvprs; p[1]       = *prsdb;
	    spVolAn[0] = *prvdel; spVolAn[1] = *delta;
	    dynHt[0]   = *prvdht;
	    DynamicHeight(2, p, inc, spVolAn, inc, dynHt, inc);
	    *dynht     = dynHt[1];
	}
}
#endif	/* FORTRAN */