c
c    @(#)svan.f	1.1
c
C  JUNE 24 1982
      real function svan(s,t,p0,sigma)
C  MODIFIED RCM
C ******************************************************
C SPECIFIC VOLUME ANOMALY (STERIC ANOMALY) BASED ON 1980 EQUATION
C OF STATE FOR SEAWATER AND 1978 PRACTICAL SALINITY SCALE.
C REFERENCES
C MILLERO, ET AL (1980) DEEP-SEA RES.,27A,255-264
C MILLERO & POISSON 1981,DEEP-SEA RES.,28A PP 625-629.
C BOTH ABOVE REFERENCES ARE ALSO FOUND IN UNESCO REPORT 38 (1981)
C UNITS:      
C       PRESSURE        P0       DECIBARS
C       TEMPERATURE     T        DEG CELSIUS (IPTS-68)
C       SALINITY        S        PSU (IPSS-78)
C       SPEC. VOL. ANA. SVAN     M**3/KG *1.0E-8
C       DENSITY ANA.    SIGMA    KG/M**3
C CHECK VALUE: SVAN = 981.30190 M**3/KG FOR S = 40 PSU,
C T = 40 DEG C, P0= 10000 DECIBARS.
C ******************************************************************
C *****************************            
C CHECK VALUE: SIGMA = 59.82037  KG/M**3 FOR S = 40 PSU,
C T = 40 DEG C, P0= 10000 DECIBARS.
C *******************************************************
      REAL P,T,S,SIG,SR,R1,R2,R3,R4
      REAL A,B,C,D,E,A1,B1,AW,BW,K,K0,KW,K35
C EQUIV       
      EQUIVALENCE (E,D,B1),(BW,B,R3),(C,A1,R2) 
      EQUIVALENCE (AW,A,R1),(KW,K0,K)
C ********************
C DATA
      DATA R3500,R4/1028.1063,4.8314E-4/
      DATA DR350/28.106331/
C   R4 IS REFERED TO AS  C  IN MILLERO & POISSON 1981
C CONVERT PRESSURE TO BARS AND TAKE SQUARE ROOT SALINITY.
      P=P0/10.
      SR = SQRT(ABS(S)) 
C *********************************************************
C PURE WATER DENSITY AT ATMOSPHERIC PRESSURE
C   BIGG P.H.,(1967) BR. J. APPLIED PHYSICS 8 PP 521-537.
C
      R1 = ((((6.536332E-9*T-1.120083E-6)*T+1.001685E-4)*T 
     X-9.095290E-3)*T+6.793952E-2)*T-28.263737
C SEAWATER DENSITY ATM PRESS. 
C  COEFFICIENTS INVOLVING SALINITY
C  R2 = A   IN NOTATION OF MILLERO & POISSON 1981
      R2 = (((5.3875E-9*T-8.2467E-7)*T+7.6438E-5)*T-4.0899E-3)*T
     X+8.24493E-1 
C  R3 = B  IN NOTATION OF MILLERO & POISSON 1981
      R3 = (-1.6546E-6*T+1.0227E-4)*T-5.72466E-3
C  INTERNATIONAL ONE-ATMOSPHERE EQUATION OF STATE OF SEAWATER
      SIG = (R4*S + R3*SR + R2)*S + R1 
C SPECIFIC VOLUME AT ATMOSPHERIC PRESSURE
      V350P = 1.0/R3500
      SVA = -SIG*V350P/(R3500+SIG)
      SIGMA=DR350-(SVA/(V350P*(V350P+SVA)))
C  SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS
      SVAN=SVA*1.0E+8
      IF(P.EQ.0.0) RETURN
C ******************************************************************
C ******  NEW HIGH PRESSURE EQUATION OF STATE FOR SEAWATER ********
C ******************************************************************
C        MILLERO, ET AL , 1980 DSR 27A, PP 255-264
C               CONSTANT NOTATION FOLLOWS ARTICLE
C********************************************************
C 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
C             
      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 
     X-0.1194975
      A = (D*SR + C)*S + AW 
C             
      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 
     X+148.4206)*T-1930.06
      K0 = (B1*SR + A1)*S + KW
C EVALUATE PRESSURE POLYNOMIAL 
C ***********************************************
C   K EQUALS THE SECANT BULK MODULUS OF SEAWATER
C   DK=K(S,T,P)-K(35,0,P)
C  K35=K(35,0,P)
C ***********************************************
      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))
C  SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS
      SVAN=SVA*1.0E+8
      V350P = V350P*PK
C  ****************************************************
C COMPUTE DENSITY ANAMOLY WITH RESPECT TO 1000.0 KG/M**3
C  1) DR350: DENSITY ANAMOLY AT 35 PSU 0 DEG. C & ATMOSPHERIC PRES.
C  2) DR35P: DENSITY ANAMOLY 35 PSU 0 DEG. C , WITH PRES. VARIATION
C  3) DVAN : DENSITY ANAMOLY VARIATIONS INVOLVING SPECFIC VOL. ANAMOLY
C ********************************************************************
C CHECK VALUE: SIGMA = 59.82037  KG/M**3 FOR S = 40 PSU,
C T = 40 DEG C, P0= 10000 DECIBARS.
C *******************************************************
      D350=GAM/PK
      DR35P=R3500*D350
      DVAN=SVA/(V350P*(V350P+SVA))
      SIGMA=DR350+DR35P-DVAN
      RETURN  
      END     
