/*
**  o2.c 1.3	Solaris 2.4 Unix	950830	SIO/ODF	fmd
**
**  Various dissolved O2 conversion functions for seawater:
**
**  double BunsenO2(double s, double t)
**	calculate the Bunsen absorption coefficient for O2
**  double Poynting(double s, double t, double p)
**	calculate the Poynting correction for pressure
**  double O2PerLiterToPerKg(double o2mlpl, double rho)
**	convert dissolved O2 concentration (ml/l) to micro-moles/Kg.
**  double O2PerKgToPerLiter(double o2umpkg, double rho)
**	convert dissolved O2 concentration (micro-moles/Kg) to ml/l.
**  double O2MllToPartPress(double o2mlpl, double s, double t, double p)
**	convert an O2 concentration in ml/l to a partial-pressure in atm.
**  double O2PartPressToMll(double o2pp, double s, double t, double p)
**	convert an O2 partial-pressure in atm to a concentration in ml/l.
**  double O2Saturation(double t, double s)
**	calculate the O2 saturation value in ml/l.
**  double O2SaturationPartPress(double s, double t, double p)
**	calculate the O2 saturation value in atm.
**
**  Plus, soon-to-be-obsolete Fortran interface routines:
**
**  double o22pp_(float *o2mlpl, float *s, float *t, float *p)
**	convert an O2 concentration in ml/l to a partial-pressure in atm.
**  double o22ml_(float *o2pp, float *s, float *t, float *p)
**	convert an O2 partial-pressure in atm to a concentration in ml/l.
**  double o2sat_(float *t, float *s)
**	calculate the O2 saturation value in ml/l.
**  double o2satpp_(float *s, float *t, float *p)
**	calculate the O2 saturation value in atm.
**  double bunseno2_(float *s, float *t)
**	calculate the Bunsen absorption coefficient for O2
**
**  References:
**  ==========
**
**	Weast, R.C., Editor, 1971, "CRC Handbook of Chemistry and Physics",
**	The Chemical Rubber Co., 1971.
**
**	Weiss, R. F., 1970, The solubility of nitrogen, oxygen and argon
**	in water and seawater. Deep-Sea Research 17, 721-735.
**
**	Dymond, J.H. and Smith, E.B., The Virial Coefficients of Pure Gases
**	and Mixtures, Oxford, 518 pp., 1980.
**
**  Poynting pressure correction from:
**
**	Hitchman, M.L., "Measurement of Dissolved Oxygen",
**	pg. 30, eqn 2.28,  John Wiley & Sons, Inc., 1978.
**
*/

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

#ifndef FORTRAN

#define	Kelvin(celsius)	((celsius)+273.15)
#define	M	  31.9988	/* Molecular weight of O2 */
#define	R	 831.432	/* Gas constant, X 10 J Kmole-1 K-1 */
#define	D	   1.42905481	/* density O2 g/l @ 0 C */

/*
**  Calculate the Bunsen absorption coefficient of dissolved O2 
**  in seawater @ s,t
*/
double					/* <- Bunsen coeff for O2 */
BunsenO2(s, t)
	double		s;		/* -> salinity (PSS 78) */
	double		t;		/* -> potential temperature (deg C) */
{
	double		kelvin, k100, bunsen;
	double		a1 =	 -58.3877;
	double		a2 =	  85.8079;
	double		a3 =	  23.8439;
	double		b1 =	  -0.034892;
	double		b2 =	   0.015568;
	double		b3 =	  -0.0019387;

    kelvin = Kelvin(t);			/* absolute potential temperature */
    k100   = kelvin       * 0.01;
    bunsen = exp(a1 + a2/k100 + a3*log(k100) +
                 s*(b1 + k100*(b2 + b3*k100)))*1000.0;
    return (bunsen);
} /* BunsenO2() */

/*
**  Calculate the Poynting correction for pressure to account
**  for isothermal variation of fugacity in a potential field.
*/
double					/* <- Poynting correction */
Poynting(s, t, p)
	double		s;		/* -> salinity (PSS 78) */
	double		t;		/* -> potential temperature (deg C) */
	double		p;		/* -> pressure (decibars) */
{
	double		IESRho();	/* insitu density @ S,T,P */
	double		kelvin = Kelvin(t),
			rhostp = IESRho(s, t, p)*0.001;	/* insitu density
							   g/cm**3 */
	return (exp(M*(p/rhostp)/(R*kelvin)));
} /* Poynting() */

/*
**  Convert dissolved O2 concentration (ml/l) to micro-moles/Kg.
*/
double					/* <- O2 uM/Kg */
O2PerLiterToPerKg(o2mlpl, rho)
	double		o2mlpl;		/* -> O2 concentration in ml/l */
	double		rho;		/* -> density (S,T,0) (Kg/M**3) */
{
	double		C =	M/D*0.001;	/* ml/M */
	return (o2mlpl/(C*rho*0.001));
}

/*
**  Convert dissolved O2 concentration (micro-moles/Kg) to ml/l.
*/
double					/* <- O2 ml/l */
O2PerKgToPerLiter(o2umpkg, rho)
	double		o2umpkg;	/* -> O2 concentration in uM/Kg */
	double		rho;		/* -> density (S,T,0) (Kg/M**3) */
{
	double		C =	M/D*0.001;	/* ml/M */
	return (o2umpkg*(C*rho*0.001));
}

/*
**  Convert dissolved O2 concentration (ml/l) to partial-pressure.
*/
double					/* <- O2 part-press (atm) */
O2MllToPartPress(o2mlpl, s, t, p)
	double		o2mlpl;		/* -> oxygen concentration in ml/l */
	double		s;		/* -> salinity (PSS 78) */
	double		t;		/* -> potential temperature (deg C) */
	double		p;		/* -> pressure (decibars) */
{
	return (o2mlpl/(BunsenO2(s, t)*Poynting(s,t,p)));
} /* O2MllToPartPress() */

/*
**  Convert dissolved O2 partial-pressure to concentration (ml/l).
*/
double					/* <- O2 ml/l */
O2PartPressToMll(o2pp, s, t, p)
	double		o2pp;		/* -> oxygen part-press (atm) */
	double		s;		/* -> salinity (PSS 78) */
	double		t;		/* -> potential temperature (deg C) */
	double		p;		/* -> pressure (decibars) */
{
	return (o2pp*BunsenO2(s, t)*Poynting(s,t,p));
} /* O2PartPressToMll() */

/*
**  Calculate O2 saturation value (solubility) in seawater.
*/
double					/* <-  O2 saturation ml/l */
O2Saturation(t, s)
	double		t;		/* -> potential temperature (deg C) */
	double		s;		/* -> salinity (PSS 78) */
{
	double		kelvin, k100;
	double		a1 =	-173.4292;
	double		a2 =	 249.6339;
	double		a3 =	 143.3483;
	double		a4 =	 -21.8492;
	double		b1 =	 - 0.033096;
	double		b2 =	   0.014259;
	double		b3 =	  -0.0017000;

	kelvin = Kelvin(t);
	k100   = kelvin *   0.01;
	return (exp(a1 + a2/k100 + a3*log(k100) + a4*k100 + 
		s*(b1 + k100*(b2 + b3*k100))));
} /* O2Saturation() */

/*
**  Calculate O2 saturation value partial-pressure in seawater.
*/
double					/* <- O2 sat part-press */
O2SaturationPartPress(s, t, p)
	double		s;		/* -> salinity (PSS 78) */
	double		t;		/* -> potential temperature (deg C) */
	double		p;		/* -> pressure (decibars) */
{
	return (O2MllToPartPress(O2Saturation(t, s), s, t, p));
} /* O2SaturationPartPress() */

#else	/* FORTRAN */

double
o22pp_(o2mlpl, s, t, p)
	float	*o2mlpl, *s, *t, *p;
{
	return (O2MllToPartPress(*o2mlpl, *s, *t, *p));
}

double
o22ml_(o2pp, s, t, p)
	float	*o2pp, *s, *t, *p;
{
	return (O2PartPressToMll(*o2pp, *s, *t, *p));
}
double
o2sat_(t, s)
	float	*t, *s;
{
	return (O2Saturation(*t, *s));
}
double
o2satpp_(s, t, p)
	float	*s, *t, *p;
{
	return (O2SaturationPartPress(*s, *t, *p));
}
double
bunseno2_(s, t)
	float *s, *t;
{
	return(BunsenO2(*s,*t));
}
#endif	/* FORTRAN */