function [SVAN,SIGMA]=svan(S,T,P0)
% function [SVAN,SIGMA]=svan(S,T,P0)
%  MODIFIED RCM
% ******************************************************
% 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.
% ***
%

      R3500 = 1028.1063;
      R4 = 4.8314e-4;
      DR350 = 28.106331;
%   R4 IS REFERED TO AS  C  IN MILLERO AND POISSON 1981
% CONVERT PRESSURE TO BARS AND TAKE SQUARE ROOT SALINITY.
      P=P0/10.;					
                clear P0
      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 = A   IN NOTATION OF MILLERO AND POISSON 1981
      R2 = (((5.3875E-9*T-8.2467E-7).*T+7.6438E-5).*T-4.0899E-3).*T...
      +8.24493E-1; 
%  R3 = B  IN NOTATION OF MILLERO AND POISSON 1981
      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;
                clear R1 R2 R3
% SPECIFIC VOLUME AT ATMOSPHERIC PRESSURE
      V350P = 1.0/R3500;
      SVA = -SIG.*V350P./(R3500+SIG);
      SIGMA=SIG+DR350;
%  SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS
      SVAN=SVA*1.0E+8;

% added second any so that program would function properly with 
% two dimensional data.  dail - 7/11/97
if any(any(P ~= 0.0)),

% ***
% ******  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;
                clear E BW
%             
      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;
                clear D C AW
%             
      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;
                clear B1 A1 KW
% 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));
                clear K35 
%  SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS
      SVAN=SVA*1.0E+8;
      V350P = V350P.*PK;
                clear 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.
% ***si
      DR35P=GAM./V350P;
      DVAN=SVA./(V350P.*(V350P+SVA));
      SIGMA=DR350+DR35P-DVAN;
end

