C a copy of CAL81 routines from B Huber
C	fixed range check bugs 5 May 1987
C	routines kept in /usr/lib/libeos.a
c  hydrographic subroutines to compute various oceanographic
c  parameters from measured values of pressure, temp, conductivity,
c  and salinity.   routines re-written & reorganized october, 1982
c  largely from &calcs, &hysub, & &rhosb.
c
c  to conform as closely as possible to unesco '81 standards
c
c  programmers: s rennie, b huber, p mele, c greengrove 
c                 & some whoi types referenced below
c
c  group i -- mostly functions returning 1 value
c             no double precision except zeta & zeta2
c
c   call sal81( p, t, c ,s) - convert conductivity => salin (practical sal'78)
c   potmp( p, t, s )     - potential temperature (bryden)
c   sigmt( t, s )        - new sigma-t function (calls dens)
c   sigth( p, t, s )     - new sigma-theta function (call dens)
c   sigp(p, t, s, pref)  - function for sigma-p (calls atg & dens )
c   spvan(p, t, s, spvl) - new function for spec vol anomaly (calls spvol)
c   alfbt(p, t, s, alpha, beta, dwrtt, dwrts) - subr for dens partial deriv
c                                               (calls dens)
c   zeta(p, zlat, zdep)  - function to convert pressure to depth (approx)
c   zeta2(p, zlat, zdep) - sub to convert pressure to depth (approx)
c   bvun(p, t, s, sn2)  - brunt-vaisalla freq in cycles/hour (changed from
c                                                             bvfof 7/19/84)
c
c  group ii -- used by routines in group i, may have double precision args
c
c   atg(p, t, s)                - adiabatic temp grad (bryden 1973)
c   dens(p, t, s, r0ts, rpts )  - unecso'81 density rho(sigm-t) & rho(insitu)
c   spvol(p, t, s, spv0, spv )  - unesco'81 specific volume
c   sbulk(p, t, s, kk)          - unesco'81 secant bulk modulas
c   theta( p, t, s, pref )      - local potential temp (fofonoff) - uses atg
c
c****************************************************************************
c

      function potmp (prs, temp, sal)

c
c potential temperature according to bryden
c convert to 1948 temperature scale
c

      tmp = temp
      tmp = temp + 4.4e-6 * temp * (100.0 - temp)

      s0 = sal - 35.0

      a3 = 0.50484e-14 * tmp - 0.16056e-12
      a2 = (0.21987e-11 * tmp - 0.31628e-9) * tmp + 0.89309e-8 -
     *      0.41057e-10 * s0
      a11 = -0.29778e-7 * tmp + 0.17439e-5
      a10 = ((0.40274e-9 * tmp - 0.54065e-7) * tmp +
     *        0.83198e-5) * tmp + 0.36504e-4
      a1 = a11 * s0 + a10

      potmp = temp - ((a3 * prs + a2) * prs + a1) * prs


      return
      end


c****************************************************************************

      function zeta(p, zlat, zdep)

c
c  depth from pressure -- ignores dynamic height anamoly
c  zlat (input) is lat in radians
c
c  from saunders & fofonoff  -  dsr 23, aug 1978
c

      real*8 zgrav, zdep, zlat

      zgrav = 978.0318 * (1.0 + 5.3024d-3 * dsin(zlat)**2 - 5.9d-6 *
     *       dsin(2.0 * zlat)**2)  
c  cm/s**2

c      zdep = 0.712953 * p + 1.113d-7 * p**2 - 3.434d-12 * p**3 +
c     *      14190.7 * dlog(1.0 + 1.83d-5 * p)
      zdep = (( -3.434d-12*p + 1.113d-7)*p + 0.712953)*p 
     *       + 14190.7 * dlog(1.0 + 1.83d-5 * p)

      zdep = zdep / (zgrav + 1.113d-4 * p) * 1000.0


      return
      end




c****************************************************************************
c
      subroutine zeta2(p, zlat, zdep)
c
c  compute depth from pressure & lat (in radians)
c  ignores dynamic height anomoly
c  from saunders  -  jpo 11/1, apr 81
c

      real*8 zlat, zdep, zpres



      zpres = p

      zdep = ((1.0d0 - (5.92d0 + 5.25d0 * dsin(zlat)**2) * 1d-3) 
     *        - 2.21d-6 * zpres)*zpres


      return
      end


c****************************************************************************

      subroutine sal81(pres, temp, cond, sal)
c
c sal81 salinity subroutine derived from
c sal78 subr ********** oct 24 1979 *************
c
c subroutine to convert conductivity to salinity
c
c algorithms recommended by jpots using the 1978 practical
c salinity scale and ipts-68 for temperature
c
c n fofonoff
c
c            code basically from rte aqui
c            it was found that subr. salin formerly used in
c            in plt78 produced values that were too low by
c            approx. .006 o/oo, possibly due to use of dauphinee
c            correction. this version incorporates the newest
c            recommendations on how to derive salinity from
c            conductivity (unesco 78), but still uses old
c            pressure term. update when new pressure data
c            available in the literature. for now, this
c            routine should produce results consistent with
c            the aquisition program.
c            bah nov 8 '81
c            mikhail somov
c

      sfn(xr, xt) = ((((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)
c
c rt35
c
      rt35(xt) = (((1.0031e-9 * xt - 6.9698e-7) * xt + 1.104259e-4) *
     *           xt + 2.00564e-2) * xt + 0.6766097
c
c cba
c
      c(xp) = ((3.989e-15 * xp - 6.370e-10) * xp + 2.070e-5) * xp
      b(xt) = (4.464e-4 * xt + 3.426e-2) * xt + 1.0
      a(xt) = -3.107e-3 * xt + 0.4215
c
c prog
c
      dt = temp - 15.0
      r = cond / 42.909
      rt = r / (rt35(temp) * (1.0 + c(pres) / (b(temp) + a(temp) * r)))
c
c  avoid neg arg to sqrt
c
      if (rt .le. 0.0) rt = 0.0
      rt = sqrt(rt)

      sal = sfn(rt, dt)


      return
      end


c****************************************************************************

      function sigmt ( t, s )    
c sigma-t from r0ts

      double precision r0ts, rho

      p = 0.0
      call dens( p, t, s, r0ts, rho )
      sigmt  = (r0ts - 1.d0) * 1.d3


      return
      end


c****************************************************************************

      function sigth( p, t, s )    
c sigma-theta from r0(potmp)s

      double precision r0ts, rho

      theta = potmp( p, t, s)
      call dens( p, theta, s, r0ts, rho )
      sigth = (r0ts - 1.d0) * 1.d3


      return
      end


c****************************************************************************

      function sigp(p, t, s, pref)

c
c * sigp ******  potential density fcn (was sigz2)*** bah 8/82
c
c          compute density of parcel moved adiabatically from
c          pressure p to pref.   uses bryden (73) polynomial (in atg )
c

      double precision r0ts, rho

      pincr = 100.0
      nmax  = 10000.0 / pincr  
c  max pr allowed is 10,000 db
      n = 0
      pi = p
      tp = t


      if(pref .lt. p) pincr = -pincr  
c  moving up or down?
c
c    compute local atg and pot. temp at p + pincr.  repeat in
c    pincr dbar increments until we reach pref.
c
      do 10 n = 1, nmax
        if ( abs(pi - pref) .lt. abs(pincr) ) then
          tp = tp + atg( pi, tp, s ) * (pref - pi)
          go to 20  
c  reached pref, we're done
        end if

        tp = tp + atg( pi,tp,s) * pincr
        pi = pi + pincr
 10   continue

      sigp = -99.999  
c  if we get here then something is wrong
      return

c    now compute density at pref, tp, s

 20   call dens( pref, tp, s, r0ts, rho )

      sigp = (rho - 1.d0) * 1.0d3
  

      return
      end


c****************************************************************************

      subroutine sbulk(pr, t, s, kk)

c
c  subroutines to calculate density, spec vol, secant bulk
c  modulas and alpha & beta
c  based on unesco 1981 report
c  equation of state for seawater - millero
c  programmer - c. greengrove, jan 1982
c  modified for hp - p mele, sep '82
c
c  range:
c    s =  0 to 42 (practical salinity)
c    t = -4 to 40 (c)
c    pr =  0 to 10000 (decibars)
c
c  other units:
c    density = kg/m3 **3
c    bulk deni mod.(k) = bars
c
c
c  kk is secant bulk modulas - returned
c

      implicit double precision (a-z)

      real*4  t, s, pr, s12  
c  single precision


       parameter
     1 (e0=19652.21d+00,e1=148.4206d+00,e2=-2.327105d+00,
     2  e3=1.360477d-02,e4=-5.155288d-05,

     3  f0=54.6746d+00,f1=-.603459d+00,f2=1.09987d-02,f3=-6.167d-05,

     4  g0=7.944d-02,g1=1.6483d-02,g2=-5.3009d-04,

     5  h0=3.239908d+00,h1=1.43713d-03,h2=1.16092d-04,h3=-5.77905d-07,

     6  i0=2.2838d-03,i1=-1.0981d-05,i2=-1.6078d-06,

     7  j=1.91075d-04,

     8  k0=8.50935d-05,k1=-6.12293d-06,k2=5.2787d-08,

     9  m0=-9.9348d-07,m1=2.0816d-08,m2=9.1697d-10)



      if (t.lt.-4.0 .or. t.gt.40.0) then  
c  range specifications
        kk = -99.9
        return
      else if (s.lt.0.0 .or. s.gt.42.0) then
        kk = -99.9
        return
      else if (pr.lt.0.0 .or. pr.gt.10000.0) then
        kk = -99.9
        return
      end if

      p = pr / 10.0  
c  convert to bars
 
c  define sqrt(s)
      s12=sqrt(s)
c
c  secant bulk modulas (k) of seawater
c
c  pure water terms of sbm are w series
c

      kw = (((e4*t + e3)*t + e2)*t + e1)*t + e0
      aw = ((h3*t + h2)*t + h1)*t + h0
      bw = (k2*t + k1)*t + k0
c
c  coeff for final equation
c
      aa = aw + s*((i2*t + i1)*t + i0 + j*s12)
      bb = bw + s*((m2*t + m1)*t + m0)
c
c  sbm at p = 0 first term in the final eq
c
      ko = kw + s*(((f3*t + f2)*t + f1)*t + f0)
     *     + s*s12*((g2*t + g1)*t + g0)
c
c  final eq sbm
c
      kk = (bb*p + aa)*p + ko


      return
      end


c****************************************************************************

      subroutine dens(pr, t, s, r0, rr)
c
c  sub to compute density
c  calls sub 'sbulk', for secant bulk modulas
c
c  r0 is density at p = 0  -  returned in gr cm**3
c  rr is in situ density   -  returned
c

      implicit double precision (a-z)

      real*4  t, s, pr  

c      dimension a(0:5), b(0:4), c(0:2)


       parameter
     1 (a0=999.842594d+00,a1=6.793952d-02,a2=-9.095290d-03,
     2   a3=1.001685d-04,a4=-1.120083d-06,a5=6.536332d-09,

     3  b0=8.24493d-01,b1=-4.0899d-03,b2=7.6438d-05,
     4  b3=-8.2467d-07,b4=5.3875d-09,

     5  c0=-5.72466d-03,c1=1.0227d-04,c2=-1.6546d-06,

     6  d=4.8314d-04)



      if (t.lt.-4.0 .or. t.gt.40.0) then
        r0 = -99.9
        rr = -99.9
        return
      else if (s.lt.0.0 .or. s.gt.42.0) then
        r0 = -99.9
        rr = -99.9
        return
      else if (pr.lt.0.0 .or. pr.gt.10000.0) then
        r0 = -99.9
        rr = -99.9
        return
      end if


      call sbulk(pr, t, s, kk)  
c  secant bulk modulas (k) of seawater
c
c  density of smow
c
      rw = ((((a5*t + a4)*t + a3)*t + a2)*t + a1)*t +a0
c
c  density at p = 0
c
      r0 = rw + s*((((b4*t + b3)*t + b2)*t + b1)*t + b0)
     *     + s*sqrt(s)*((c2*t + c1)*t + c0) + s*s*d 
c
c  in situ density
c
      p = pr / 10.0  
c  p is in bars
      rr = r0 / (1.d0 - p / kk)

      rr = rr / 1.d3  
c  densities are returned in
      r0 = r0 / 1.d3  
c  grams / cubic centimeter


      return
      end


c****************************************************************************

      subroutine spvol(pr, t, s, spv0, spv)

c
c  sub to compute specific volume
c  calls sub dens, for density
c
c  spv0 is specific volume at p = 0  -  returned
c  spv is in situ specific volume    -  returned
c

      implicit double precision (a-z)

      real*4  t, s, pr  
c  single precision



      if (t.lt.-4.0 .or. t.gt.40.0) then
        spv0 = -99.9
        spv = -99.9
        return
      else if (s.lt.0.0 .or. s.gt.42.0) then
        spv0 = -99.9
        spv = -99.9
        return
      else if (pr.lt.0.0 .or. pr.gt.10000.0) then
        spv0 = -99.9
        spv0 = -99.9
        return
      end if

      call dens(pr, t, s, r0, rr)  
c  this sub calls 'sbulk'


      spv0 = 1.0 / r0  
c  specific volume at p = 0

      spv = 1.0 / rr  
c  in situ specific volume


      return
      end


c****************************************************************************

      subroutine alfbt(p, t, s, alph, beta, dwrtt, dwrts)

c
c  sub to compute alpha & beta
c  calls subs:  sbulk   -  secant bulk modulas
c               dens    -  density
c               spvol   -  specific volume
c

      implicit double precision (a-z)

      real*4  p, t, s, alph, beta, dwrtt, dwrts, s12
c  arguments single precision

      dimension a(0:5), b(0:4), c(0:2), e(0:4), f(0:3), g(0:2),
     *          h(0:3), i(0:2), k(0:2), m(0:2)

       data
     1 a / 999.842594d+00, 6.793952d-02, -9.095290d-03,
     2     1.001685d-04, -1.120083d-06, 6.536332d-09 /,

     3 b / 8.24493d-01, -4.0899d-03, 7.6438d-05,
     4    -8.2467d-07, 5.3875d-09 /,

     5 c / -5.72466d-03, 1.0227d-04, -1.6546d-06 /,

     6 d / 4.8314d-04 /,

     7 e / 19652.21d+00, 148.4206d+00, -2.327105d+00,
     8     1.360477d-02, -5.155288d-05 /,

     9 f / 54.6746d+00, -.603459d+00, 1.09987d-02, -6.167d-05 /


       data
     1 g / 7.944d-02, 1.6483d-02, -5.3009d-04 /,

     2 h / 3.239908d+00, 1.43713d-03, 1.16092d-04, -5.77905d-07 /,

     3 i / 2.2838d-03, -1.0981d-05, -1.6078d-06 /,

     4 j / 1.91075d-04 /,

     5 k / 8.50935d-05, -6.12293d-06, 5.2787d-08 /,

     6 m / -9.9348d-07, 2.0816d-08, 9.1697d-10 /



      if (t.lt.-4.0 .or. t.gt.40.0) then
        alph = -99.9
        beta = -99.9
        return
      else if (s.lt.0.0 .or. s.gt.42.0) then
        alph = -99.9
        beta = -99.9
        return
      else if (p.lt.0.0 .or. p.gt.10000.0) then
        alph = -99.9
        beta = -99.9
        return
      end if

      call sbulk(p, t, s, kk)                    
c  need kk

      call dens(p, t, s, r0, rr)                 
c  need r0

      call spvol(p, t, s, spv0, spv)             
c  need spv0 & spv

c
c  compute sqrt(s)
       s12=sqrt(s)
c  derivatives working toward alpha and beta
c

      dbt = k(1) + 2 * k(2) * t + (m(1) + 2 * m(2) * t) * s   
c  derv b wrt t

      dbs = m(0) + t*(m(1) + m(2)*t)                     
c  derv b wrt s

      dat = h(1) + t*(2*h(2) +3*h(3)*t) + (i(1) +   
c  derv a wrt t
     *      2 * i(2) * t) * s

      das = i(0) + t*(i(1)*t + i(2)*t) + 1.5 * j * s12     
c  derv a wrt s

c      dkot = e(1) + 2 * e(2) * t + 3 * e(3) * t**2 + 4 *      
c     *       e(4) * t**3 + (f(1) + 2 * f(2) * t + 3 * f(3) *
c     *       t**2) * s + (g(1) + 2 * g(2) * t) * s**1.5
      dkot = ((4*e(4)*t + 3*e(3))*t + 2*e(2))*t + e(1)  
c  derv ko wrt t
     *        + s*(f(1) + (3*f(3)*t + 2*f(2))*t + s12*
     *        (2*g(2)*t + g(1)))

      dkos = f(0) + ((f(3)*t + f(2))*t + f(1))*t 
c  derv ko wrt s
     *       + 1.5 * (g(0) + t*(g(1) + g(2)*t)) * s12

      drt = a(1)+(((5*a(5)*t + 4*a(4))*t + 3*a(3))*t + 2*a(2))*t
c  derv dens
     *     + s*(b(1) + ((4*b(4)*t + 3*b(3))*t + 2* b(2))*t 
c  (p = 0) wrt t
     *     + s12*(c(1) + 2*c(2)*t))

      dwrtt = drt           
c 
 argument returned

      drs = b(0)+ (((b(4)*t + b(3))*t + b(2))*t + b(1))*t
c  derv dens
     *      + 1.5*s12*(c(0) + (c(2)*t + c(1))*t) + 2*d*s
c  (p = 0) wrt s
	dwrts = drs          
c 
 argument returned
      r0sq=r0*r0

      dvot = (-1.0 / r0sq) * drt              
c  derv spec vol (p = 0) wrt t

      dvos = (-1.0 / r0sq) * drs              
c  derv spec vol (p = 0) wrt s

      pbar = p/10.

      dkt = dkot+pbar*(dat+dbt*pbar)                  
c  derv k wrt t

      dks = dkos+pbar*(das+dbs*pbar)                  
c  derv k wrt s
      kksq=kk*kk
      fact1=(1.-pbar/kk)
      fact2=spv0*pbar/(kk*kk)

      dspvt = dvot*fact1 + dkt*fact2
c  derv spec vol wrtt


      dspvs = dvos*fact1 + dks*fact2
c  derv spec vol wt
c
c  alpha & beta
c
      alph =  dspvt/spv
      beta = -dspvs/spv


      return
      end


c****************************************************************************

      function atg(p, t, s)

c
c adiabatic temperature gradient (bryden 1973)
c

      ds = s - 35.0
      atg = (((-2.1687e-16 * t + 1.8676e-14) * t - 4.6206e-13) * p +
     *      ((2.7759e-12 * t - 1.1351e-10) * ds + ((-5.4481e-14 * t +
     *      8.733e-12) * t - 6.7795e-10) * t + 1.8741e-8)) * p +
     *      (-4.2393e-8 * t + 1.8932e-6) * ds + ((6.6228e-10 * t -
     *      6.836e-8) * t + 8.5258e-6) * t + 3.5803e-5


      return
      end


c****************************************************************************

      function theta(p0, t0, s, pf)

c
c to compute local potential temperature at pf
c
c oct 12 1975 n. fofonoff
c

      p = p0
      t = t0
      h = pf - p
      xk = h * atg(p, t, s)

      t = t + 0.5 * xk
      q = xk
      p = p + 0.5 * h
      xk = h*atg(p,t,s)

      t = t + 0.29298322*(xk-q)
      q = 0.58578644*xk + 0.121320344*q
      xk = h*atg(p,t,s)

      t = t + 1.707106781*(xk-q)
      q = 3.414213562*xk - 4.121320344*q
      p = p + 0.5*h
      xk = h*atg(p,t,s)

      theta = t + (xk - 2.0 * q) / 6.0

      return
      end


c****************************************************************************

      function spvan(p,t,s,spv)

c
c specific volume anomaly*1e5 unesco routines, simple-minded approach
c oct ,1982, s.rennie
c

      double precision xspv0,xspv, xstnd

      call spvol( p, t, s, xspv0, xspv)
      call spvol( p, 0.0, 35., xspv0, xstnd )
      spvan = 1d5*( xspv - xstnd )

      spv = xspv       
c return single prec. arg for spec.vol.


      return
      end


c****************************************************************************

      function svel(pr,t,sal)

c
c svel wilson
c wilson oct sound speed (m/sec) jasa,1960,32,(10),1357
c

      p = 0.1019716*(pr+10.1325)
      sd = sal - 35.
   10 a = (((7.9851e-6*t-2.6045e-4)*t-4.4532e-2)*t+4.5721)*t
     x+1449.14
      svel = (7.7711e-7*t-1.1244e-2)*t+1.39799

      v0 = (1.69202e-3*sd+svel)*sd+a
      a = ((4.5283e-8*t+7.4812e-6)*t-1.8607e-4)*t+.160272
      svel = (1.579e-9*t+3.158e-8)*t+7.7016e-5

      v1 = svel*sd+a
      a = (1.8563e-9*t-2.5294e-7)*t+1.0268e-5
      svel = -1.2943e-7*sd+a

      a = -1.9646e-10*t+3.5216e-9
      svel = (((-3.3603e-12*p+a)*p+svel)*p+v1)*p+v0


      return
      end


c****************************************************************************

      function oxsat(t,s)

c
c oxygen saturation (ml/l) weiss,1970 dsr 17,(4);721
c

      x = (t+273.16)/100.0
      oxsat = exp(((-21.8492*x-173.4292)*x+249.6339)/x
     *+s*((-0.0017*x+0.014259)*x-0.033096)+143.3483*alog(x))


      return
      end


c****************************************************************************

      subroutine ctdo2 ( p,t,s,oxc,oxt,ox,pcor,tcor,c2)

c
c    ctdo2       oxygen sensor algorithm
c     uses weiss ( dsr,17,(4); 721,1970 )
c     formula for saturation.  units (ml/l)
c

      ox = oxc*exp(tcor*(t+c2*(oxt-t))+pcor*p)
      ox = ox*oxsat(t,s)


      return
      end


c****************************************************************************

      function bvun(p,t,s,  sn2)

c
c bvun -- from bvfof ***** brunt-vaisala freq *****
c ************************************
c sept 25 1976 n fofonoff
c
c computes n in cycles per hour,n**2 in rad/sec**2
c

       double precision  rlast, rho, r0 , e
	SAVE
      if ( p.eq. 0.0 ) then
	sn2 = 0.0
	bvun = 0.0
	go to 90
      end if

      pav = 0.5*(p + plast)
      call dens( pav,theta(plast,tlast,slast,pav),slast, r0,rlast)
      call dens(pav,theta(p,t,s,pav),s, r0, rho)

cc bvfof      e   = 38.467369d+0 *( rho - rlast) /
cc bvfof     @                ((p-plast)* (2.0*rlast*rho + rlast + rho)**2 )
cc      changed 7/19/84 to correct approx. factor of 2 difference from version
cc      on hp & previous calculations

	e = 9.8/rho * (rho - rlast )/ (p - plast)

      sn2 = e

      bvun = 572.9578 * dsign( dsqrt(dabs(e)),e)         
c rem: e is double prc

 90   plast = p       
c warning: routine needs to 'remember' the last
      tlast = t       
c values. if bvfof is in a segment, the
      slast = s       
c last values will be lost when swapped.


      return
      end