/*
      REAL FUNCTION SVEL (S,T,P0)                                              *
C-------------------------------------------------------------------------------
C
C     Function to calculate sound speed in sea water.                          *
C
C     Input parameters:                                                        *
C            S      -  Salinity PSS-78                                         *
C            T      -  Temperature in degrees C                                *
C            P0     -  Pressure in decibars                                    *
C     Output value:                                                            *
C            SVEL   -  Speed in meters/second                                  *
C
C     Reference:  Chen and Millero, 1977, JASA,62,1129-1135
C
C     Checkvalue: SVEL=1731.995 M/S, S=40 (PSS-78),T=40 DEG C, P=10000 Dbar
C
C-------------------------------------------------------------------------------
      EQUIVALENCE (A0,B0,C0),(A1,B1,C1),(A2,C2),(A3,C3)
*/

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

#define A0 a0
#define B0 a0
#define C0 a0
#define A1 a1
#define B1 a1
#define C1 a1
#define A2 a2
#define C2 a2
#define A3 a3
#define C3 a3

double SVEL (double S, double T, double P0)
{
     double a0, a1, a2, a3;
     double A, B, C, D, P, SR;
     double svel;
/*     =====>  Scale pressure to bars  */

      P=P0/10;
      SR = sqrt(fabs(S));
/*    =====>   S**2 term  */
      D = 1.727E-3 - 7.9836E-6*P;
/*    =====>  S**3/2 term; */
      B1 = 7.3637E-5 + 1.7945E-7*T;
      B0 = -1.922E-2 - 4.42E-5*T;
      B = B0 + B1*P;
/*    =====>  S**1 term */
      A3 = (-3.389E-13*T + 6.649E-12)*T + 1.100E-10;
      A2 = ((7.988E-12*T - 1.6002E-10)*T + 9.1041E-9)*T - 3.9064E-7;
      A1 = (((-2.0122E-10*T + 1.0507E-8)*T - 6.4885E-8)*T
           - 1.2580E-5)*T + 9.4742E-5;
      A0 = (((-3.21E-8*T + 2.006E-6)*T + 7.164E-5)*T - 1.262E-2)*T
           + 1.389;
      A = ((A3*P + A2)*P + A1)*P + A0;
/*    =====>  S**2 term */
      C3 = (-2.3643E-12*T + 3.8504E-10)*T - 9.7729E-9;
      C2 = (((1.0405E-12*T - 2.5335E-10)*T + 2.5974E-8)*T
           - 1.7107E-6)*T + 3.1260E-5;
      C1 = (((-6.1185E-10*T + 1.3621E-7)*T - 8.1788E-6)*T
            + 6.8982E-4)*T + 0.153563;
      C0 = ((((3.1464E-9*T - 1.47800E-6)*T + 3.3420E-4)*T
           - 5.80852E-2)*T + 5.03711)*T + 1402.388;
      C = ((C3*P + C2)*P + C1)*P + C0;
/*    =====>  Sound Speed....  */
      svel = C + (A + B*SR + D*S)*S;

     return(svel);
}
#undef A0
#undef B0
#undef C0
#undef A1
#undef B1
#undef C1
#undef A2
#undef C2
#undef A3
#undef C3

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

   if (argc != 4)
   {
      printf("\n USAGE: svel <salinity> <temperature> <pressure>");
      printf("\n EXAMPLE: svel 40.0 40.0 10000.0\n");
   }
   else
   {
      S = atof(argv[1]);
      T = atof(argv[2]);
      P = atof(argv[3]);
   }
   printf("\n sound speed = %f m/s", SVEL(S, T, P));
   printf("\n for S = %g (PSS-78), T = %g deg C, P = %g dbar\n\n", S, T, P);
   return 0;
}
#endif