#======================================================================

#          L I B E O S 8 3 . P L 

#          doc: Mon Mar  8 08:22:05 1999

#          dlm: Tue Feb  5 16:42:28 2013

#          (c) 1999 A.M. Thurnherr

#          uE-Info: 20 59 NIL 0 0 72 2 2 4 NIL ofnI

#======================================================================

# Perl Implementation of UNESCO Eqn of State 1983

# Notes:

# 	- copied from eos83calc.c which was in turn...

#	- copied from pexec_v4 (incl comments)

#	- only subset of functions implemented

#	- some attempt at cleaning up the code has been made

#	- pressures in dbar throughout

#	- no temperature scale assumed; set PARAM ITS=90|68

#	- T90*1.00024=T68

#	- check values calculated with T68

#	- no conductivity unit assumed; set PARAM cond.unit=S/m|mS/cm

# HISTORY:

#	Mar 08, 1999: - translated by hand from krc

#				  - added &alpha(), &beta(), &Rrho(), &TurnerAngle()

#	Mar 13, 1999: - cosmetic changes

#	Mar 14, 1999: - BUG NAN instead of NaN

#	Mar 21, 1999: - make $sigmaR optional in &sVolAnom()

#	Mar 31, 1999: - alias &potemp() = &theta(); &podens() = &sigma()

#				  - added &temp() to calc in-situ from potemp

#	Sep 18, 1999: - parameter typechecking

#	Aug 28, 2000: - added PARAM T68 check

#	Aug 29, 2000: - forced temp_scale param check --- set to T68 or T90

#	Sep 25, 2000: - changed temp_scale to 68 or 90 (easier check)

#				  - check for temp_scale during loading

#	Nov 07, 2000: - added &dynHt()

#				  - strictified

#	Nov 13, 2000: - removed temp_scale check during loading (STUPID,

#					because PARAMs are only available after header

#					is read)

#	Feb 28, 2001: - added &rho(S,T,P)

#	Mar 28, 2001: - changed Rrho() to use podens

#				  - optimized &theta() for P == Pref

#	Apr  2, 2001: - &TurnerAngle() disabled pending adaptation to new Rrho()

#	Apr  3, 2001: - updated &TurnerAngle()

#	Jul 17, 2001: - cosmetics

#				  - added &grav(), &potErgAnom()

#	Jul 18, 2001: - cosmetics

#	Jul 20, 2001: - changed temp_scale to ITS

#	Nov 26, 2001: - added &g(), &f()

#	Dec 26, 2005: - update notes

#   Jul  1, 2006: - Version 3.3 [HISTORY]

#				  - nan must be quoted for use strict

#	Nov  5, 2006: - added K15toSalin

#	Oct 18, 2007: - adapted to new &antsFunUsage()

#				  - changed order of &salin() params

#	Oct 19, 2007: - continued

#				  - removed &grav(), podens(), potemp()

#	Dec  1, 2007: - made theta(), sigma(), rho() return nan on nan input

#	Dec 21, 2007: - BUG: grav() was still used

#	Jan 20, 2008: - BUG: theta(), sigma(), rho() still generated error on nan in

#				  - made depth(), dynht(), potErgAnom() return nan on nan in

#	Jan  4, 2011: - maded salin(), sVel() return nan on nan in

#	Oct  6, 2011: - added %cond.unit (analogous to %ITS)

require "$ANTS/libvec.pl";

use strict;

#======================================================================

# PART I: stuff taken from PEXEC

#======================================================================

{ # BEGIN STATIC SCOPE

	my($TCONV);

	sub TCONV()

	{

		unless (defined($TCONV)) {

			my($ITS) = &antsRequireParam('ITS');

			if ($ITS == 68) {

				$TCONV = 1;

			} elsif ($ITS == 90) {

				$TCONV = 1.00024;

			} else {

				croak("$0: illegal PARAM-value ITS=$ITS\n");

	        }

	    }

	    return $TCONV;

	}

} # END STATIC SCOPE

#---------------------------------------------------------------------- 

#	ADIABATIC TEMPERATURE GRADIENT DEG C/BAR							 

#	REF: BRYDEN,H.,1973,DEEP-SEA RES.,20,401-408						 

#	CHECK VALUE: ATGR80=3.255976E-3 FOR S=40 NSU,T=40 DEG C,			 

#	PIN=10000 DECIBARS

#	NB: real check value appears to be 3.255976E-4; this makes sense

#	    as check value of potential temperature below is correct. Also

#       note the 0.1 factor on return.

#---------------------------------------------------------------------- 

{ my(@fc);

	sub adiaTempGrad(@)

	{

		my($S,$T,$P) = &antsFunUsage(3,'fff','[salin, temp, press(db)]',

									 \@fc,'salin','temp','press',@_);

		my($DS);

		my($TCONV) = &TCONV();

		$T *= $TCONV;		# use T68 

		$P *= 0.1;

		$DS = $S - 35.0;

		return 0.1 * ((((-2.1687E-13*$T+1.8676E-11)*$T-4.6206E-10)*$P

					+((2.7759E-10*$T-1.1351E-8)*$DS+((-5.4481E-12*$T

					+8.733E-10)*$T-6.7795E-8)*$T+1.8741E-6))*$P

					+(-4.2393E-7*$T+1.8932E-5)*$DS

					+((6.6228E-9*$T-6.836E-7)*$T+8.5258E-5)*$T+3.5803E-4) / $TCONV;

	}

}

#---------------------------------------------------------------------- 

# TO COMPUTE LOCAL POTENTIAL TEMPERATURE AT PR							 

#	USING BRYDEN 1973 POLYNOMIAL FOR ADIABATIC LAPSE RATE				 

#	AND RUNGE-KUTTA 4-TH ORDER INTEGRATION ALGORITHM.					 

#	REF: BRYDEN,H.,1973,DEEP-SEA RES.,20,401-408						 

#	FOFONOFF,N.,1977,DEEP-SEA RES.,24,489-491							 

#	CHECK VALUE: PTMP83 =36.89072 FOR S=40 NSU,T=40 DEG C,				 

#	P(MEASURED PRESSURE)=10000 DECIBARS,Pref=0 DECIBARS

#---------------------------------------------------------------------- 

{ my(@fc);

	sub theta(@)

	{

		my($S,$T,$P,$Pref) =

			&antsFunUsage(4,'....','[salin, temp, press(db),] refpress(db)',

						  \@fc,'salin','temp','press',undef,@_);

		return 'nan' unless numberp($S) && numberp($T) && numberp($P);

		my($H,$XK,$Q);

		my($TCONV) = &TCONV();

		return $T if ($P == $Pref);

		$T *= $TCONV;	# use T68 

		$H = $Pref - $P;

		$XK = $H * &adiaTempGrad($S,$T/$TCONV,$P)*$TCONV;

		$T += 0.5 * $XK;

		$Q = $XK;

		$P += 0.5 * $H;

		$XK = $H * &adiaTempGrad($S,$T/$TCONV,$P)*$TCONV;

		$T += 0.29289322 * ($XK-$Q);

		$Q = 0.58578644*$XK + 0.121320344*$Q;

		$XK = $H * &adiaTempGrad($S,$T/$TCONV,$P)*$TCONV;

		$T += 1.707106781 * ($XK-$Q);

		$Q = 3.414213562*$XK - 4.121320344*$Q;

		$P += 0.5*$H;

		$XK = $H * &adiaTempGrad($S,$T/$TCONV,$P)*$TCONV;

		return ($T + ($XK-2.0*$Q)/6.0)/$TCONV;

	}

}

# can't easily do default fields because theta field name is not unique

sub temp(@)

{

	my($S,$T,$P,$Pref) =

		&antsFunUsage(4,"ffff","salin, potemp, press(db), refpress(db)",@_);

	return &theta($S,$T,$Pref,$P);

}

#---------------------------------------------------------------------- 

#	Calculation of specific volume anomaly

#		Gill (1982), p.215: - specific volume = inverse of density

#							- SVA referenced to same press, T=0, S=35

#																		 

#	 *************************************************					 

#	 *** USES EQUATION OF STATE FOR SEA WATER 1980 ***					 

#	 *************************************************					 

#																		 

#	Check value. SVAN83=981.30210E-8, SIGMA=59.82037					 

#	for S=40, T=40, P=10000												 

#																		 

# 	NB: - this one was a bitch to translate, used f2c and it shows :-)	 

#		- the 1e-8 factor in the check value is wrong; yup: see gill p.215

#		  (typical units are 1e-8m^3kg^-1)

#		- check val appears to be slightly off, I get 981.30187319561

#	  

#---------------------------------------------------------------------- 

sub sVolAnom(@)

{

	my($S,$T,$P,$sigmaR) =

		&antsFunUsage(-3,'fff','salin, temp, press(db)[, ref to sigma var]',@_);

	my($r3500) = 1028.1063;

	my($r4)	   = 4.8314e-4;

	my($dr350) = 28.106331;

	my($temp0,$temp1,$temp2,$temp3,$temp4);

	my($dvan,$dr35p,$p,$t,$dk,$pk,$sr,$gam,$sig,$rk35,$sva,$v350p);

	my($TCONV) = &TCONV();

	$t = $T*$TCONV;								# use T68 

	$p = $P / 10.;

#	$sr = sqrt(abs($S));

	$sr = sqrt($S);

	$temp3 = (((($t * 6.536332e-9 - 1.120083e-6) * $t + 

	  1.001685e-4) * $t - .00909529) * $t + .06793952) * $t 

	  - 28.263737;

	$temp2 = ((($t * 5.3875e-9 - 8.2467e-7) * $t + 7.6438e-5) 

	  * $t - .0040899) * $t + .824493;

	$temp1 = ($t * -1.6546e-6 + 1.0227e-4) * $t - .00572466;

	$sig = ($r4 * $S + $temp1 * $sr + $temp2) * $S + $temp3;

	$v350p = 1. / $r3500;

	$sva = -$sig * $v350p / ($r3500 + $sig);

	$$sigmaR = $sig + $dr350 if (defined($sigmaR));

	return $sva * 1e8 if ($p == 0.0);

	$temp0 = ($t * 9.1697e-10 + 2.0816e-8) * $t - 9.9348e-7;

	$temp1 = ($t * 5.2787e-8 - 6.12293e-6) * $t + 3.47718e-5;

	$temp1 = $temp1 + $temp0 * $S;

	$temp0 = 1.91075e-4;

	$temp2 = ($t * -1.6078e-6 - 1.0981e-5) * $t + .0022838;

	$temp3 = (($t * -5.77905e-7 + 1.16092e-4) * $t + .00143713) * $t - .1194975;

	$temp3 = ($temp0 * $sr + $temp2) * $S + $temp3;

	$temp0 = ($t * -5.3009e-4 + .016483) * $t + .07944;

	$temp2 = (($t * -6.167e-5 + .0109987) * $t - .603459) * $t + 54.6746;

	$temp4 = ((($t * -5.155288e-5 + .01360477) * $t - 2.327105) * $t + 148.4206) * $t - 1930.06;

	$temp4 = ($temp0 * $sr + $temp2) * $S + $temp4;

	$dk = ($temp1 * $p + $temp3) * $p + $temp4;

	$rk35 = ($p * 5.03217e-5 + 3.359406) * $p + 21582.27;

	$gam = $p / $rk35;

	$pk = 1. - $gam;

	$sva = $sva * $pk + ($v350p + $sva) * $p * $dk / ($rk35 * ($rk35 + $dk));

	$v350p *= $pk;

	$dr35p = $gam / $v350p;

	$dvan = $sva / ($v350p * ($v350p + $sva));

	$$sigmaR = $dr350 + $dr35p - $dvan if (defined($sigmaR));

	return $sva * 1e8;

}

#----------------------------------------------------------------------

# Dynamic Height (from dyht83.F)

#	Usage:

#		- dynHt(salin,temp,press[,idx])

#	Check Values from PEXEC:

#		- dynHt(40,18,20)  -> -.01879

#		- dynHt(??,16,50)  -> -.03194

#		- dynHt(??,14,100) -> -.00255

#	NB: - check values 2 & 3 do not appear to check out

#		- checked against WHOI-supplied BBTRE data indicates max

#		  diff of 1e-3 dyn.m at 5141m

#		- use different idx to use multiple times in single program

#		  (c.f. [gshear])

#		- behaves the same as the pexec version with PFIRST=0

#		- dynamic height seems to be defined (1e-5 factor) to allow

#		  station distance to be in km and velocities in cm/s

#----------------------------------------------------------------------

{ # BEGIN static scope

	my(@lastP,@lastAnom,@ht);

	sub dynHt(@)

	{

		my($S,$T,$P,$idx) =

			&antsFunUsage(-3,"...","salin, temp, press(db)[, <idx>]",@_);

		return 'nan' unless numbersp($S,$T,$P);

		my($anom) = sVolAnom($S,$T,$P) * 1e-5;

		if (!defined($ht[$idx])) {					# first call

			$ht[$idx] = $P * $anom;

		} else {									# successive calls

		    croak("$0: pressure not increasing monotonically ($lastP[$idx] -> $P)\n")

		    	if ($P < $lastP[$idx]);

			$ht[$idx] += ($P-$lastP[$idx]) * ($anom+$lastAnom[$idx])/2;

		}

		$lastP[$idx] = $P; $lastAnom[$idx] = $anom;

		return $ht[$idx];

	}

} # END static scope

#----------------------------------------------------------------------

# Local Gravity (from grav83.F)

#	Usage:

#		- g(press,lat)

#	Check Values:

#		- g(10000,30) = 9.804160

#----------------------------------------------------------------------

{ my(@fc);

	sub g(@)

	{

		my($P,$lat) = &antsFunUsage(2,'ff','[press(db), lat]',\@fc,'press','%lat',@_);

		my($x) = sin($lat/57.29578) ** 2;				    

		my($g) = 9.780318*(1.0+(5.2788e-3+2.36e-5*$x)*$x);	# global variation

		return $g + 1.092e-6 * $P;							# pressure correction

	}

}

#----------------------------------------------------------------------

# Coriolis frequency

#	Usage:

#		- f(%lat)

#----------------------------------------------------------------------

{ my(@fc);

	sub f(@)

	{

		my($lat) = &antsFunUsage(1,'f','[lat]',\@fc,'%lat',@_);

		my($Omega) = 7.292e-5;								# Gill (1982)

		return 2 * $Omega * sin(rad($lat));

	}

}

#----------------------------------------------------------------------

# Potential Energy Anomaly (from pean83.F)

#	Usage:

#		- potErgAnom(salin,temp,press,refPress,lat)

#	Check Values from PEXEC: (NB: sequence of calls)

#		- potErgAnom(40,18, 20,0,49.2235) = -1.91462

#		- potErgAnom(38,16, 50,0,49.2235) = -4.31092

#		- potErgAnom(36,14,100,0,49.2235) = 24.82922

#	Units:

#		- Gill (1982) p.45: m^2/s^2 = J/kg

#	NB:

#		- geopotential Phi is PEA wrt zero pressure

#		- PE of volume is volume integral of rho*PEA (c.f. Gill, p.80)

#		- check value lat calculated so that g(0,lat) ~ 9.81

#		- IF Pref<P ON 1ST CALL, HEIGHT WILL BE CALCULATED FROM LEVEL

#		  Pref.  HENCE IF Pref = 0.0, FROM SEA SURFACE DOWN.

#		- different to pexec, this version sets values for P<Pref to nan

#----------------------------------------------------------------------

{ # BEGIN STATIC SCOPE

	my($gzero,$H,$lastZ,$lastP);

	sub potErgAnom(@)

	{

		my($S,$T,$P,$Pref,$lat) =

	        &antsFunUsage(5,".....","salin, temp, press(db), refpress(db), lat",@_);

	    return 'nan' unless numbersp($S,$T,$P,$Pref,$lat);

		return 'nan' if ($P < $Pref);

		$gzero = g(0,$lat) unless (defined($gzero));		# 1st time

		my($g) = $gzero + 1.113e-4*$P;

		my($anom) = sVolAnom($S,$T,$P) * 1e-3;

		my($Z) = $anom * $P/$g;

		if (defined($H)) {									# not 1st time

		    croak("$0: pressure not increasing monotonically\n")

		    	if ($P < $lastP);

			$H += ($Z+$lastZ) * ($P-$lastP)*0.5;

		} else {

			$H = (($anom*$Pref/$g)+$Z)*($P-$Pref)*0.5;

        }

		$lastP = $P; $lastZ = $Z;

		return $H;

    }

} # END static scope

#---------------------------------------------------------------------- 

# Density at pressure P (P - Measured pressure)						 

# 						(PREF - Reference pressure)					 

# Equation of state for seawater proposed by JPOTS 1980				 

# References															 

#		Millero et al 1980, Deep Sea Res.,27A,255-264					 

#		Jpots Ninth Report 1978, Tenth Report 1980						 

# Units:																 

#	   Pressure   P  Decibars									 

#	   Temperature T  Deg Celcius (IPTS-68)						 

#	   Salinity   S  NSU (IPSS-78)								 

#	   Density   RHO KG/M**3										 

#	   Spec. Vol.  EOS80 M**3/KG									 

#	check value. 43.331642 for P=10000,PREF=5000,T=40,S=40				 

# NB: check value appears to be 42.33164!!!							 

#---------------------------------------------------------------------- 

{ my(@fc); 

	sub sigma(@)

	{

		my($S,$T,$P,$Pref) =

			&antsFunUsage(4,'....','[salin, temp, press(db),] refpress(db)',

						  \@fc,'salin','temp','press',undef,@_);

		return 'nan' unless numberp($S) && numberp($T) && numberp($P);

		my($sig);

		&sVolAnom($S,&theta($S,$T,$P,$Pref),$Pref,\$sig);

		return $sig;

	}

}

{ my(@fc);

	sub rho(@)

	{

		my($S,$T,$P) = &antsFunUsage(3,'...','[salin, temp, press(db)]',

									 \@fc,'salin','temp','press',@_);

		return 'nan' unless numberp($S) && numberp($T) && numberp($P);

		return 1000 + &sigma($S,$T,$P,$P);

	}

}

#---------------------------------------------------------------------- 

#	FUNCTION TO CONVERT CONDUCTIVITY TO SALINITY ACCORDING TO THE		 

#	ALGORITHMS RECOMMMENDED BY JPOTS USING THE 1978 PRACTICAL 			 

#	SALINITY SCALE (IPSS-78) AND IPTS-68 FOR TEMPERATURE. 				 

#	C1535=CONDUCTIVITY (FROM CULKIN&SMITH,1980,J.OCEAN.ENG.VOL.5 		 

#    PP22-23) = 42.9140 AT    PRES=0.0, 							 

#	SALINITY 35 NSU AND TEMPERATURE 15 DEG CELSIUS (IPTS-68) 			 

#	PRESSURE=DECIBARS 													 

#	RETURNS ZERO FOR CND<.0005 		 									 

#	CHECKVALUES.... 													 

#	SAL83 =40.000000 FOR CND=81.025545,T=40 DEG C,PIN=10000 			 

#   DECIBARS														 

#	WRITTEN BY N FOFONOFF; REVISED OCT 6 1980 							 

#	FCN SAL83, XR=SQRT(RT) 												 

#	DERIVATIVE WRT XR ; DSAL/DXR 										 

#	RT35 																 

#	C,B,A, POLYNOMIALS 													 

#	NB: - done using f2c (couldn't be bothered)							 

#		- removed SAL->COND												 

#---------------------------------------------------------------------- 

{ my(@fc);

  my($cond_scale);	# 1 for mS/cm, 10 for S/m

	sub salin(@)

	{

		my($C,$T,$P) = &antsFunUsage(3,'...','[cond, temp, press(db)]',

									 \@fc,'cond','temp','press',@_);

		return 'nan' unless numberp($C) && numberp($T) && numberp($P);

		my($r__,$dt,$rt,$c1535);

		my($TCONV) = &TCONV();

		unless (defined($cond_scale)) {				# deal with different conductivity units

			my($cu) = &antsRequireParam('cond.unit');

			if    ($cu eq 'S/m')   { $cond_scale = 10; }

			elsif ($cu eq 'mS/cm') { $cond_scale = 1;  }

			else { croak("$0: illegal PARAM-value cond.unit=$cu\n"); }

		}

		$C *= $cond_scale;

		return 0.0 if ($C <= 5e-4); 				# zero salinity trap 

		$T *= $TCONV;								# use T68 scale 

		$c1535 = 42.914;

		$dt = $T - 15.0;

		$P *= .1;									# convert pressure to bars 

		$r__ = $C / $c1535; 						# convert cond to salin

		$rt = $r__ / ((((($T * 1.0031e-9 - 6.9698e-7) * $T + 

		  1.104259e-4) * $T + .0200564) * $T + .6766097) * ((($P 

		  * 3.989e-12 - 6.37e-8) * $P + 2.07e-4) * $P / (

		  ($T * 4.464e-4 + .03426) * $T + 1. + ($T *

		  -.003107 + .4215) * $r__) + 1.));

		$rt = sqrt(abs($rt));

		return (((($rt * 2.7081 - 7.0261) * $rt + 14.0941) *

		   $rt + 25.3851) * $rt - .1692) * $rt + .008 + 

		  $dt / ($dt * .0162 + 1.) * ((((($rt * -.0144 + 

		  .0636) * $rt - .0375) * $rt - .0066) * $rt - 

		  .0056) * $rt + 5e-4);

	}

}

#----------------------------------------------------------------------

# K15toSalin; see http://ioc.unesco.org/oceanteacher/OceanTeacher2/01_GlobOcToday/02_CollDta/02_OcDtaFunda/03_T&SScales/TemperatureAndSalinityScales.htm

#----------------------------------------------------------------------

sub K15toSalin(@)

{

	my($K15) = &antsFunUsage(1,'f','K15',@_);

	return  0.0080 -  0.1692 * $K15**0.5

				   + 25.3851 * $K15**1.0

				   + 14.0941 * $K15**1.5

				   -  7.0261 * $K15**2.0

				   +  2.7081 * $K15**2.5;

}

#---------------------------------------------------------------------- 

#	DEPTH IN METERS FROM PRESSURE IN DECIBARS USING						 

#	SAUNDERS AND FOFONOFF'S METHOD.										 

#	DEEP SEA RES., 1976,23,109-111.										 

#	FORMULA REFITTED FOR EOS80											 

#	CHECK VALUE. 9712.654 M FOR PIN=10000 DECIBARS,LATITUDE=30 DEG		 

#	..CONVERT PRESSURE TO BARS

#	NB: results are closer to Saunders, JPO 11, 1981 than ref given

#	    above. They are not exact, either, however.

#---------------------------------------------------------------------- 

{ my(@fc);

	sub depth(@)

	{

		my($P,$lat) = &antsFunUsage(2,'..','[press(db), lat]',\@fc,'press','%lat',@_);

		return 'nan' unless numbersp($P,$lat);

		my($x) = sin($lat/57.29578);

		$P *= 0.1;

		$x = $x*$x;

		return ((((-1.82E-11*$P+2.279E-7)*$P-2.2512E-3)*$P+97.2659)*$P) /

			  (9.780318*(1.0+(5.2788E-3+2.36E-5*$x)*$x) + 1.092E-5*$P);

	}

}

#----------------------------------------------------------------------

#	Convert depth to pressure using 2nd order polynomial

#	From file pdepth.F (pexec)

#	This is based on Saunders, JPO 11, 1981

#----------------------------------------------------------------------

{ my(@fc);

	sub press(@)

	{

		my($d,$lat) = &antsFunUsage(2,'ff','[depth, lat]',\@fc,'depth','%lat',@_);

		my($c1) = 5.92E-3 + 5.25E-3 * sin($lat/57.29578) ** 2;

		my($c2) = 2.21E-6;

		my($press) = (1 - $c1 - sqrt((1 - $c1)**2 - 4*$d*$c2)) / (2*$c2);

		my($derr) = abs(&depth($press,$lat) - $d);

		&antsInfo("WARNING (libEOS83.pl): %.1gm depth error due to pressure " .

				  "approximation", $derr) if ($derr >= 1);

		return $press;

	}

}

#---------------------------------------------------------------------- 

#	SOUND SPEED SEAWATER (CHEN & MILLERO 1977, JASA,62,1129-1135) 		 

#	SPEED IN M/S, P IN DECIBARS, T DEG C (IPTS-68), S NSU(IPSS-78)