/*
**  cfcsat.c 1.03	Solaris 2.3 Unix	940906	SIO/ODF	fmd
**
**  Calculate CFC-12, CFC-13 saturation values in sea water:
**
**  void CFCSetAirVals(double f11AirVal, double f12AirVal)
**	specify cfc-11 and cfc-12 air values.
**  double CFC11Saturation(double theta, double salt)
**	calculate cfc-11 saturation (solubility) in seawater.
**  double CFC12Saturation(double theta, double salt)
**	calculate cfc-12 saturation (solubility) in seawater.
**
**  Plus soon-to-be obsolete Fortran interface routines.
**
**  References:
**  ==========
**
**       Warner, M., Weiss, R. F., 1988
*/
#include <math.h>
#include <odflib.h>

#ifndef FORTRAN

#define F11DefAirVal	230.0
#define F12DefAirVal	410.0

static	double	F11AirVal=F11DefAirVal, F12AirVal=F12DefAirVal;

void
CFCSetAirVals(f11AirVal, f12AirVal)
	double f11AirVal;		/* cfc-11 air value pM/kg */
	double f12AirVal;		/* cfc-12 air value pM/kg */
{
	F11AirVal = f11AirVal;
	F12AirVal = f12AirVal;
}

#define	Kelvin(celsiusT)	((celsiusT)+273.15)

/*
**  Calculate cfc-11 saturation (solubility) in seawater.
*/
double				/* <- cfc-11 saturation (pM/kg) */
CFC11Saturation(theta, salt)
	double theta;		/* -> potential temperature (degrees C) */
	double salt;		/* salinity (psu) */
{
	double	a, b, kelvin, k100;
#define t11c1	-232.0411
#define t11c2	 322.5546
#define t11c3	 120.4956
#define t11c4	  -1.39165
#define s11c1	-0.146531
#define s11c2	 0.093621
#define s11c3	-0.0160693

	kelvin = Kelvin(theta);
	k100   = kelvin * 0.01;
	a      = t11c1 + t11c2/k100 + t11c3*log(k100) + t11c4*k100*k100;
	b      = salt*(s11c1 + k100*(s11c2 + s11c3*k100));
	return (exp(a+b)*F11AirVal);
} /* CFC11Saturation() */

/*
**  Calculate cfc-12 saturation (solubility) in seawater.
*/
double				/* <- cfc-12 saturation (pM/kg) */
CFC12Saturation(theta, salt)
	double theta;		/* -> potential temperature (degrees C) */
	double salt;		/* salinity (psu) */
{
	double	a, b, kelvin, k100;
#define t12c1	-220.2120
#define t12c2	 301.8695
#define t12c3	 114.8533
#define t12c4	  -1.39165
#define s12c1	-0.147718
#define s12c2	 0.093175
#define s12c3	-0.0157340

	kelvin = Kelvin(theta);
	k100   = kelvin * 0.01;
	a      = t12c1 + t12c2/k100 + t12c3*log(k100) + t12c4*k100*k100;
	b      = salt*(s12c1 + k100*(s12c2 + s12c3*k100));
	return (exp(a+b)*F12AirVal);
} /* CFC12Saturation() */

#else /* FORTRAN */

void
frsetav_(f12av, f11av)
	float	*f12av, *f11av;
{
	CFCSetAirVals(*f11av, *f12av);
}

void
frsat_(fr12, fr11, theta, salt, f12sat, f11sat, f12spc, f11spc)
	float	*fr12,   *fr11,   *theta,  *salt,
		*f12sat, *f11sat, *f12spc, *f11spc;
{
	double	df12sat, df11sat;

	*f11sat =
	df11sat = CFC11Saturation( *theta, *salt);
	*f12sat =
	df12sat = CFC12Saturation( *theta, *salt);

	if (isnan(*fr12))
	    *f12spc = quiet_nan();
	else
	    *f12spc = *fr12/df12sat*100.0;

	if (isnan(*fr11))
	    *f11spc = quiet_nan();
	else
	    *f11spc = *fr11/df11sat*100.0;
}
#endif /* FORTRAN */