/*
      REAL FUNCTION SAL78 (CND,T,P,M)                                          *
C-------------------------------------------------------------------------------
C
C                            Function = SAL78
C
C     Function to convert conductivity ratio to salinity (M=0) or              *
C     salinity to conductivity ratio (M=1).                                    *
C
C     Input parameters:                                                        *
C            CND    -  Conductivity Ratio (M=0)                                *
C            T      -  Temperature in degrees C                                *
C            P      -  Pressure in decibars                                    *
C            M      -  Flag indicating either ratio or salinity                *
C
C
C     The conductivity ratio (CND = 1.0000000 for salinity = 35, PSS-78
C     temperature = 15.0 deg C, and atmospheric pressure.
C
C     References:  Also located in UNESCO Report No. 37, 1981
C                  Practical Salinity Scale 1978: E.L. Lewis,IEEE Ocean
C                  Engineering, Jan., 1980
C
C     Checkvalues:
C            SAL78=1.888091 :CND= 40.00,T=40 deg C,P=10000 dbars: M=1
C            SAL78=40.00000 :CND=1.888091,T=40 deg C,P=10000 dbars: M=0
C
C     SAL78 Ratio: Returns zero for conductivity ratio:  < 0.0005
C     SAL78: Returns zero for salinity:  < 0.02
C
C---------------------------------------------------------------------
C     =====>  Internal functions...
C             Practical salinity scale 1978 definition with temperature
C             correction
C             XT=T-15.0 : XR=SQRT(RT)
*/
#include <stdio.h>
#include <math.h>
#include "ocean.h"

double SAL(double XR, double XT)
{
     return( ((((2.7081*XR-7.0261)*XR+14.0941)*XR+25.3851)*XR
                  -0.1692)* XR+0.0080
                  +(XT/(1.0+0.0162*XT))*(((((-0.0144*XR+
                 0.0636)*XR-0.0375)*XR-0.0066)*XR-0.0056)*XR+0.0005) );
}
/*  =====>  DSAL(XR,XT) function for derivative of SAL(XR,XT) with XR.*/
double DSAL(double XR, double XT)
{
     return( ((((13.5405*XR-28.1044)*XR+42.2823)*XR+50.7702)*XR
        -0.1692)+(XT/(1.0+0.0162*XT))*((((-0.0720*XR+0.2544)*XR
        -0.1125)*XR-0.0132)*XR-0.0056) );
}

/* ====> Function RT35 : C(35,T,0)/C(35,15,0) variation with temperature */
double RT35(double XT)
{
     return( (((1.0031E-9*XT-6.9698E-7)*XT+1.104259E-4)*XT
         + 2.00564E-2)*XT  + 0.6766097 );
}
/*  =====>  Polynomials of RP: C(S,T,P)/C(S,T,0) variation with pressure
            C(XP) polynomials corresponds to A1-A3 constants:  Lewis 1980
*/
double C(double XP)
{
     return( ((3.989E-15*XP-6.370E-10)*XP+2.070E-5)*XP );
}
double B(double XT)
{
     return( (4.464E-4*XT+3.426E-2)*XT + 1.0 );
}
/* =====>  A(XT) Polynomial corresponds to B3 & B4 constants: Lewis 1980
*/
double A(double XT)
{
     return( -3.107E-3*XT + 0.4215 );
}
/*********** Main subroutine, SAL78, starts here. *********************/
double SAL78(double CND, double T, double P, double M)
{
   /* declarations of internal functions */
   double A(), B(), C(), RT35(), SAL(), DSAL();
     double sal78 = 0.0;
     double DT;
     double R, RT, SI, dels;
     double RTT, AT, BT, CP;
     int n;
/*
C     =====>  Zero salinity/conductivity trap
*/
     if ( (M == 0) && (CND <= 5E-4) ) return(sal78);
     if ( (M == 1) && (CND <= 0.02) ) return(sal78);
     DT = T - 15.0;
/*
C     =====>  Select branch for salinity (M=0) or conductivity (M=1)
*/
     if ( M == CR_TO_SAL)     /* Convert  conductivity to salinity */
     {
          R = CND;
          RT = R/(RT35(T)*(1.0 + C(P)/(B(T) + A(T)*R)));
          RT = sqrt(fabs(RT));
          sal78 = SAL(RT,DT);
          return(sal78);
     }
/*
C     =====>  Invert salinity to conductivity by the Newton-Raphson
C             iterative method.
C
C     =====>  First approximation
*/
      RT = sqrt(CND/35.0);
      SI = SAL(RT,DT);
/*
C     =====>  Iteration loop begins here with a maximum of 10 cycles
*/
     for( n=0, dels=1.0; (n<10) && (dels >= 1.E-4); n++)
     {
          RT = RT + (CND - SI)/DSAL(RT,DT);
          SI = SAL(RT,DT);
          dels = fabs(SI - CND);
     }
/*
C     =====>  Compute conductivity ratio
*/
     RTT = RT35(T)*RT*RT;
     AT = A(T);
     BT = B(T);
     CP = C(P);
     CP = RTT*(CP + BT);
     BT = BT - RTT*AT;
/*
C     =====>  Solve quadratic equation for R: R = RT35 * RT * (1+C/AR+B)
*/
     R = sqrt( fabs(BT*BT + 4.0*AT*CP) ) - BT;
/*    ======>  Conductivity return      */
     sal78 = 0.5*R/AT;
     return(sal78);
}
/************* end of function SAL78 ****************************/
#ifdef EXE
void main(void)
{
   printf("SAL78=1.888091 :CND= 40.00,T=40 deg C,P=10000 dbars: M=SAL_TO_CR\n");
   printf( "compare conductivity ratio: %f\n",
                        SAL78(40.0,40.0,10000,SAL_TO_CR) );
   printf("SAL78=40.00000 :CND=1.888091,T=40 deg C,P=10000 dbars: M=CR_TO_SAL\n");
   printf( "compare salinity: %f\n", SAL78(1.888091,40.0,10000,CR_TO_SAL) );
}
#endif