/*
**  potentemp.c	1.02	Solaris 2.3 Unix	940906	SIO/ODF	fmd
**
**  Calculate local potential temperature of sea water for an
**  arbitrary reference pressure.
**
**  References:
**  ==========
**
**	Fofonoff, N. P., 1977. Computation of Potential Temperature of
**	Seawater for an Arbitrary Reference Pressure. Deep-Sea Research
**	24, 489-491.
**	Bryden, H. L., 1973. New Polynomials for Thermal Expansion,
**	Adiabatic Temperature Gradient, Deep-Sea Research 20, 401-408.
*/

#include <math.h>
#include <odflib.h>

#ifndef FORTRAN

#ifdef __GNUC__
#define INLINESTATIC	__inline static
#else  __GNUC__
#define INLINESTATIC	static
#endif __GNUC__

/*
**  Calculate adiabatic temperature gradient
*/
INLINESTATIC
double				/* <- adiabatic temp gradient deg C/decibar */
atg(s, t, p)
	double s;		/* -> salinity (psu) */
	double t;		/* -> insitu temperature (degrees C) */
	double p;		/* -> pressure (decibars) */
{
	return (((((-2.1687e-16)*t +
		     1.8676e-14)*t
		    -4.6206e-13)*p +
		  (( 2.7759e-12 *t
		    -1.1351e-10)*s +
		 (((-5.4481e-14)*t +
		     8.7330e-12)*t
		    -6.7795e-10)*t +
		     1.8741e-8))*p +
		   (-4.2393e-8  *t +
		     1.8932e-6) *s +
		  (((6.6228e-10)*t
		    -6.8360e-8) *t +
		     8.5258e-6) *t +
		     3.5803e-5);
}

/*
**  Calculate potential temperature.
*/
double				/* <- potential temperature (degrees C) */
PotentialTemp(s, t, p, rp)
	double s;		/* -> salinity (psu) */
	double t;		/* -> temperature (degrees C) */
	double p;		/* -> pressure (decibars) */
	double rp;		/* -> reference pressure (decibars) */
{
	double	q, x, s1, dp;

	s1 =  s - 35.0;
	dp = rp - p;
    /*
    **  Runge-Kutta 4th-order integration.
    */
	x  = dp*atg(s1,t,p);

	t += 0.5*x;
	q  = x;
	p += 0.5*dp;
	x  = dp*atg(s1,t,p);

	t += 0.292893220*(x-q);
	q  = 0.585786440*x + 0.121320344*q;
	x  = dp*atg(s1,t,p);

	t += 1.707106781*(x-q);
	q  = 3.414213562*x - 4.121320344*q;
	p += 0.5*dp;
	x  = dp*atg(s1,t,p);

	t += (x-q-q)/6.0;

	return (t);

} /* PotentialTemp() */

#else  /* FORTRAN */
void
potmp_(press, temp, salt, rp, t)
	float	*press, *temp, *salt, *rp, *t;
{
	*t = PotentialTemp(*salt, *temp, *press, *rp);
}
#endif /* FORTRAN */