/*
**  transport.c	1.01	BSD 4.3	920110	SIO/ODF	fmd
**
**  Calculate vertical transport function.
**
**  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.
*/
#include <math.h>
#include <odflib.h>

#ifndef	FORTRAN

/*
**  Calculate vertical transport.
*/
void
Transport(n, z, zInc, dynHt, dInc, transport, tInc)
	int	n;		/* -> n intervals in z, dynHt and transport */
	double	*z;		/* -> depth series (meters) */
	int	zInc;		/* -> byte increment of depth series */
	double	*dynHt;		/* -> dynamic height series (dynamic meters) */
	int	dInc;		/* -> byte increment of dynamic height series */
	double	*transport;	/* <- transport series (m**3/(10*sec**2)) */
	int	tInc;		/* -> byte increment of transport series */
{
	int	i;
	double	*pZ, *nZ, *pD, *nD, *pT, *nT;
    /*
    **  When calling Transport() repeatedly with a two-element
    **  series, the first call should be with a one-element series to
    **  initialize the starting value (see xprtf_(), below).
    */
	if (n!=2) {			/* initialize the series */
	    /*
	    **  If the integration starts from > 15 db, calculate
	    **  transport relative to starting place. Otherwise,
	    **  calculate from surface.
	    */
		if (z[0] > 15.0)
		    transport[0] = 0.0;
		else
		    transport[0] = dynHt[0]*z[0];
	}
    /*
    **  Calculate the rest of the series by integrating dynamic
    **  height with respect to depth. 
    */
	pZ = z; pD = dynHt; pT = transport;
	nZ = (double *)((char *)pZ+zInc);
	nD = (double *)((char *)pD+dInc);
	nT = (double *)((char *)pT+tInc);
	for (i=1; i<n;
	    i++,
	    pZ = nZ, pD = nD, pT = nT,
	    nZ = (double *)((char *)pZ+zInc),
	    nD = (double *)((char *)pD+dInc),
	    nT = (double *)((char *)pT+tInc))
	    nT[0] = pT[0]+(pD[0]+nD[0])*0.5*(nZ[0]-pZ[0]);
	     /* transport function m**3/(10*sec**2) */

} /* Transport() */

#else	/* FORTRAN */

void
xprtf_(depth, prvdep, dynht, prvdht, prvtrn, transp)
	float	*depth, *prvdep, *dynht, *prvdht, *prvtrn, *transp;
{
	double	z[2], dynHt[2], transport[2];
	int	inc = sizeof (double);

    /*
    **  The FORTRAN interface is called once per transport value.
    **  The first call should be made with *prvtrn==quiet_nan() to initialize
    **  the surface (or first) value.
    */
	if (isnan(*prvtrn)) {		/* initialize */
	    z[0]     = *depth;
	    dynHt[0] = *dynht;
	    Transport(1, z, inc, dynHt, inc, transport, inc);
	    *transp  = transport[0];
	} else {
	    z[0]         = *prvdep; z[1]     = *depth;
	    dynHt[0]     = *prvdht; dynHt[1] = *dynht;
	    transport[0] = *prvtrn;
	    Transport(2, z, inc, dynHt, inc, transport, inc);
	    *transp      = transport[1];
	}
}
#endif	/* FORTRAN */