function [SVAN,SIGMA]=swstate(S,T,P0,ftn);
% SWSTATE State equation for seawater
%
% [SVAN,SIGMA]=SWSTATE(S,T,P) returns the specific volume
% anomaly SVAN (m^3/kg*1e-8) and the density anomaly SIGMA (kg/m^3)
% given the salinity S (ppt), temperature T (deg C) and pressure
% P (dbars).
%
% [dVdT,dRdT]=SWSTATE(S,T,P,'dT') returns derivatives w.r.t.
% temperature of the volume and density.
%
% [dVdS,dRdS]=SWSTATE(S,T,P,'dS') returns derivatives w.r.t.
% salinity.
%
% [dVdP,dRdP]=SWSTATE(S,T,P,'dP') returns derivatives w.r.t.
% pressure.
%
% All elements can be scalars, vectors, or matrices but should be
% the same size.
%Notes: RP (WHOI 2/dec/91)
%
% This stuff is directly copied from the UNESCO algorithms, with some
% minor changes to make it Matlab compatible (like adding ";" and changing
% "*" to ".*" when necessary.
%
% RP (WHOI 3/dec/91)
%
% Added first derivative calculations.
derivT=0;
derivS=0;
derivP=0;
if (nargin==4),
if (ftn=='dT'), derivT=1;
elseif (ftn=='dS'), derivS=1;
elseif (ftn=='dP'), derivP=1;
else error('swstate: Unrecognized option!');
end;
end;
% ******************************************************
% SPECIFIC VOLUME ANOMALY (STERIC ANOMALY) BASED ON 1980 EQUATION
% OF STATE FOR SEAWATER AND 1978 PRACTICAL SALINITY SCALE.
% REFERENCES
% MILLERO, ET AL (1980) DEEP-SEA RES.,27A,255-264
% MILLERO AND POISSON 1981,DEEP-SEA RES.,28A PP 625-629.
% BOTH ABOVE REFERENCES ARE ALSO FOUND IN UNESCO REPORT 38 (1981)
% UNITS:
% PRESSURE P0 DECIBARS
% TEMPERATURE T DEG CELSIUS (IPTS-68)
% SALINITY S (IPSS-78)
% SPEC. VOL. ANA. SVAN M**3/KG *1.0E-8
% DENSITY ANA. SIGMA KG/M**3
% ******************************************************************
% CHECK VALUE: SVAN=981.3021 E-8 M**3/KG. FOR S = 40 (IPSS-78) ,
% T = 40 DEG C, P0= 10000 DECIBARS.
% CHECK VALUE: SIGMA = 59.82037 KG/M**3 FOR S = 40 (IPSS-78) ,
% T = 40 DEG C, P0= 10000 DECIBARS.
%HECK VALUE: FOR S = 40 (IPSS-78) , T = 40 DEG C, P0= 10000 DECIBARS.
% DR/DP DR/DT DR/DS
% DRV(1,7) DRV(2,3) DRV(1,8)
%
% FINITE DIFFERENCE WITH 3RD ORDER CORRECTION DONE IN DOUBLE PRECSION
%
% 3.46969238E-3 -.43311722 .705110777
%
% EXPLICIT DIFFERENTIATION SINGLE PRECISION FORMULATION EOS80
%
% 3.4696929E-3 -.4331173 .7051107
%
% (RP...I think this ---------^^^^^^ should be -.4431173!);
% *******************************************************
% DATA
R3500=1028.1063;
R4=4.8314E-4;
DR350=28.106331;
% CONVERT PRESSURE TO BARS AND TAKE SQUARE ROOT SALINITY.
P=P0/10.;
SAL=S;
SR = sqrt(abs(S));
% *********************************************************
% PURE WATER DENSITY AT ATMOSPHERIC PRESSURE
% BIGG P.H.,(1967) BR. J. APPLIED PHYSICS 8 PP 521-537.
%
R1 = ((((6.536332E-9*T-1.120083E-6).*T+1.001685E-4).*T ...
-9.095290E-3).*T+6.793952E-2).*T-28.263737;
% SEAWATER DENSITY ATM PRESS.
% COEFFICIENTS INVOLVING SALINITY
R2 = (((5.3875E-9*T-8.2467E-7).*T+7.6438E-5).*T-4.0899E-3).*T+8.24493E-1;
R3 = (-1.6546E-6*T+1.0227E-4).*T-5.72466E-3;
% INTERNATIONAL ONE-ATMOSPHERE EQUATION OF STATE OF SEAWATER
SIG = (R4*S + R3.*SR + R2).*S + R1;
% SPECIFIC VOLUME AT ATMOSPHERIC PRESSURE
V350P = 1.0/R3500;
SVA = -SIG*V350P./(R3500+SIG);
SIGMA=SIG+DR350;
V0 = 1.0./(1000.0 + SIGMA);
% SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS
SVAN=SVA*1.0E+8;
if (derivS), % These are derivatives for (S,T,0).
R4S=9.6628E-4;
RHO1 = 1000.0 + SIGMA;
RHOS=R4S*SAL+1.5.*R3.*SR+R2;
V0S=-RHOS./(RHO1.*RHO1);
elseif (derivT),
R1 =(((3.268166E-8*T-4.480332E-6).*T+3.005055E-4).*T...
-1.819058E-2).*T+6.793952E-2;
R2 = ((2.155E-8*T-2.47401E-6).*T+1.52876E-4).*T-4.0899E-3;
R3 = -3.3092E-6*T+1.0227E-4;
RHO1 = 1000.0 + SIGMA;
RHOT = (R3.*SR + R2).*SAL + R1;
V0T = -RHOT./(RHO1.*RHO1);
end;
% ******************************************************************
% ****** NEW HIGH PRESSURE EQUATION OF STATE FOR SEAWATER ********
% ******************************************************************
% MILLERO, ET AL , 1980 DSR 27A, PP 255-264
% CONSTANT NOTATION FOLLOWS ARTICLE
%********************************************************
% COMPUTE COMPRESSION TERMS
E = (9.1697E-10*T+2.0816E-8).*T-9.9348E-7;
BW = (5.2787E-8*T-6.12293E-6).*T+3.47718E-5;
B = BW + E.*S; % Bulk Modulus (almost)
% CORRECT B FOR ANAMOLY BIAS CHANGE
Bout = B + 5.03217E-5;
if (derivS),
DBDS=E;
elseif (derivT),
BW = 1.05574E-7*T-6.12293E-6;
E = 1.83394E-9*T +2.0816E-8;
BT = BW + E.*SAL;
end;
%
D = 1.91075E-4;
C = (-1.6078E-6*T-1.0981E-5).*T+2.2838E-3;
AW = ((-5.77905E-7*T+1.16092E-4).*T+1.43713E-3).*T-0.1194975;
A = (D*SR + C).*S + AW;
% CORRECT A FOR ANAMOLY BIAS CHANGE
Aout = A + 3.3594055;
if (derivS),
DADS=2.866125E-4*SR+C;
elseif (derivT),
C = -3.2156E-6*T -1.0981E-5;
AW = (-1.733715E-6*T+2.32184E-4).*T+1.43713E-3;
AT = C.*SAL + AW;
end;
B1 = (-5.3009E-4*T+1.6483E-2).*T+7.944E-2;
A1 = ((-6.1670E-5*T+1.09987E-2).*T-0.603459).*T+54.6746;
KW = (((-5.155288E-5*T+1.360477E-2).*T-2.327105).*T+148.4206).*T-1930.06;
K0 = (B1.*SR + A1).*S + KW;
if (derivS),
K0S=1.5*B1.*SR+A1;
KS=(DBDS.*P+DADS).*P+K0S;
elseif (derivT),
B1 = -1.06018E-3*T+1.6483E-2;
% APRIL 9 1984 CORRECT A1 BIAS FROM -.603457 !!!
A1 = (-1.8501E-4*T+2.19974E-2).*T-0.603459;
KW = ((-2.0621152E-4*T+4.081431E-2).*T-4.65421).*T+148.4206;
K0T = (B1.*SR+A1).*SAL + KW;
KT = (BT.*P + AT).*P + K0T;
end;
% EVALUATE PRESSURE POLYNOMIAL
% ***********************************************
% K EQUALS THE SECANT BULK MODULUS OF SEAWATER
% DK=K(S,T,P)-K(35,0,P)
% K35=K(35,0,P)
% ***********************************************
DK = (B.*P + A).*P + K0;
K35 = (5.03217E-5*P+3.359406).*P+21582.27;
GAM=P./K35;
PK = 1.0 - GAM;
SVA = SVA.*PK + (V350P+SVA).*P.*DK./(K35.*(K35+DK));
% SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS
SVAN=SVA*1.0E+8; % Volume anomaly
V350P = V350P.*PK;
% ****************************************************
% COMPUTE DENSITY ANAMOLY WITH RESPECT TO 1000.0 KG/M**3
% 1) DR350: DENSITY ANAMOLY AT 35 (IPSS-78), 0 DEG. C AND 0 DECIBARS
% 2) DR35P: DENSITY ANAMOLY 35 (IPSS-78), 0 DEG. C , PRES. VARIATION
% 3) DVAN : DENSITY ANAMOLY VARIATIONS INVOLVING SPECFIC VOL. ANAMOLY
% ********************************************************************
% CHECK VALUE: SIGMA = 59.82037 KG/M**3 FOR S = 40 (IPSS-78),
% T = 40 DEG C, P0= 10000 DECIBARS.
% *******************************************************
DR35P=GAM./V350P;
DVAN=SVA./(V350P.*(V350P+SVA));
SIGMA=DR350+DR35P-DVAN; % Density anomaly
K=K35+DK;
VP=1.0-P./K;
V = (1.) ./(SIGMA+1000.0);
if (derivS),
VS=V0S.*VP+V0.*P.*KS./(K.*K);
SVAN=VS; % dVdS
SIGMA=-VS./(V.*V); % dRdS
elseif (derivT),
VT = V0T.*VP + V0.*P.*KT./(K.*K);
SVAN=VT; % dVdT
SIGMA=-VT./(V.*V); % dRdT
elseif (derivP),
DKDP = 2.0*Bout.*P + Aout;
% CORRECT DVDP TO PER DECIBAR BY MULTIPLE *.1
DVDP = -.1*V0.*(1.0 - P.*DKDP./K)./K;
SVAN=DVDP; % dVdP
SIGMA=-DVDP./(V.*V); % dRdP
end;