/*      REAL FUNCTION SVAN (S,T,P0,SIGMA)                                        *
C-------------------------------------------------------------------------------
C
C
C     Specific volume anomaly (steric anomaly) based on 1980 equation of       *
C     state for seawater AMD 1978 Practical Salinity Scale.                    *
C
C     Input parameters:                                                        *
C            S      -  Salinity    (PSS-78)                                    *
C            T      -  Temperature in degrees C (IPTS-68)                      *
C            P0     -  Pressure in decibars                                    *
      Output parameter:
C            SIGMA  -  Density anomaly in Kg/M**3                              *
C     Output value:                                                            *
C            SVAN   -  Specific volume anomaly in 1.0E-8 M**3/KG               *
C
C     References: Millero, et al 1980 Deep-Sea Research,27A,255-264
C                 Millero & Poisson, 1981, Deep-Sea Research,28A,625-695
C                 Both references are also found in UNESCO Report #38,1981
C
C     Checkvalue: SVAN = 981.30210 E-8 M**3/Kg for S = 40 (PSS-78) ,
C                 T = 40 deg C, P0=10000 dbars
C     Checkvalue: SIGMA = 59.82037 Kg/M**3 for S = 40 (PSS-78),
C                 T = 40 deg C, P0=10000 dbars
C
C-------------------------------------------------------------------------------
C     =====>  R4 is referred to as C in Millero and Poisson 1981
*/

#include <stdio.h>
#include <math.h>
#include "ocean.h"

#define R3500  1028.1063
#define R4     4.8314E-4
#define DR350  28.106331

double SVAN (double S, double T, double P0, double *SIGMA)
/****!! major difference between C and Fortran versions:
        SVAN must be called with a pointer to the variable
        SIGMA, since SIGMA is to be set by the function.
******/
{
     double P, SIG, SR, R1, R2, R3;
     double A,B,C,D,E,A1,B1,AW,BW,K0,KW,K35;
     double V350P, SVA, svan, DK, GAM, PK, DR35P, DVAN;
/*    EQUIVALENCE (E,D,B1),(BW,B,R3),(C,A1,R2)   (only to save space?)
      EQUIVALENCE (AW,A,R1),(KW,K0,K)
*/
/*  =====>  Convert pressure to bars and take square root salinity.... */
      P=P0/10.;
      SR= sqrt(fabs(S));
/*
C     =====>  Pure water density at atmospheric pressure....
C             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;
/*     =====>  Sea water density atm. pressure.
C             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;
/*     =====>  Specific volume at atmosperic pressure...*/
      V350P = 1.0/R3500;
      SVA = -SIG*V350P/(R3500+SIG);
      *SIGMA=SIG+DR350;
/*     =====>  Scale specific volume to normally reported units....*/
      svan=SVA*1.0E+8;
      if(P == 0.0) return (svan);
/*
C =====>  New high pressure equation of state for seawater <=====
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;

      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;

      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;
/*
C     =====>  Evaluate pressure polynomial....
C             K = The secant bulk modulus of seawater
C             DK=K(S,T,P)-K(35,0,P)
C             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 volume anomaly to normally reported units */
      svan=SVA*1.0E+8;
      V350P = V350P*PK;
/*
C     =====>  Compute density anomaly with respect to 1000.0 Kg/M**3
C             1) DR350: Density anomaly at 35 (PSS-78), 0 deg C & 0 dbars
C             2) DR35P: Density anomaly at 35 (PSS-78), 0 deg C,
C                       pressure variation
C             3) DVAN : Density anomaly variations involving specifc volume
C                       anomaly
C
C     Checkvalue: SIGMA = 59.82037 Kg/M**3 for S = 40 (PSS-78), T = 40 deg C,
C                 P0= 10000 dbars
*/
      DR35P=GAM/V350P;
      DVAN=SVA/(V350P*(V350P+SVA));
      *SIGMA=DR350+DR35P-DVAN;

      return(svan);
}
#undef R3500
#undef R4
#undef DR350

#ifdef EXE
int main(int argc, char *argv[])
{
   double sigma, T = 40.0, S = 40.0, P = 10000.0;

   if (argc != 4)
   {
      printf("\n USAGE: svan <salinity> <temperature> <pressure>");
      printf("\n EXAMPLE: svan 40.0 40.0 10000.0\n");
   }
   else
   {
      S = atof(argv[1]);
      T = atof(argv[2]);
      P = atof(argv[3]);
   }
   printf("\n specific volume anomaly = %f cl/t", SVAN(S, T, P, &sigma));
   printf("\n         density anomaly = %f kg/m**3", sigma);
   printf("\n for S = %g (PSS-78), T = %g deg C, P = %g dbar\n\n", S, T, P);
    /* Recall that the function returns specific volume anomaly
       in units of 1.e-8 m3/kg.
    */
   return 0;
}
#endif