      function dens_unesco (temp, sal, pres)
c--------------------------------------------
      real*8 rh, rh0

ccccccccccccccccccccccccccccccccccccccccccccccc
c Gilles's variant:
c     situ = theta (pres, temp, sal, 0.)
c then, for insitu dens (for poten. dens pres = 0.):
c     call dens_eos (pres, situ, sal, rh0, rh)
c     dens_unesco = 1.e3 * real(rh - 1.d0) 
c
c "sigth" variant:       
c      call dens_eos (pres, thet, sal, rh0, rh)
c      dens_unesco = 1.e3 * real(rh0 - 1.d0) 
ccccccccccccccccccccccccccccccccccccccccccccccc

      call dens_eos (pres, temp, sal, rh0, rh)
      dens_unesco = 1.e3 * real(rh - 1.d0) 
      
      return
      end

      function sdens_pnt (temp, sal, pres)
c-------------------------------------------------
      include 'comm_para.h'
      include 'comm_new.h'

      if     (isalt .eq. 1) then
         sdens_pnt = SIGMA0 - TCOEF * (temp - TEMP_BOT) 
      elseif (isalt .eq. 2) then
         sdens_pnt = pdens1 (temp)
      elseif (isalt .eq. 3) then
         sdens_pnt = pdens4 (temp, sal, pres)
      elseif (isalt .eq. 4) then
         sdens_pnt = pdens12 (temp, sal, pres)
      elseif (isalt .eq. 5) then
         sdens_pnt = pdens17 (temp, sal, pres)
      elseif (isalt .eq. 6) then
         situ = theta_eos (pres, temp, sal, 0.) 
         sdens_pnt = dens_unesco (situ, sal, pres)
      elseif (isalt .eq. 7) then
         sdens_pnt = dens_marina (temp)
      endif

      return
      end

      function pdens_pnt (thet, sal)
c-------------------------------------------------
      include 'comm_para.h'
      include 'comm_new.h'

      if     (isalt .eq. 1) then
         pdens_pnt = SIGMA0 - TCOEF * (thet - TEMP_BOT) 
      elseif (isalt .eq. 2) then
         pdens_pnt = pdens1 (thet)
      elseif (isalt .eq. 3) then
         pdens_pnt = pdens4 (thet, sal, 0.)
      elseif (isalt .eq. 4) then
         pdens_pnt = pdens012 (thet, sal)
      elseif (isalt .eq. 5) then
         pdens_pnt = pdens017 (thet, sal)
      elseif (isalt .eq. 6) then
         pdens_pnt = dens_unesco (thet, sal, 0.)
      elseif (isalt .eq. 7) then
         pdens_pnt = dens_marina (thet)
      endif

      return
      end

c-----------------------------------------------------
      subroutine dens_init (npt, nz, nzi, t, sal, pdens, dens, h) 
c-----------------------------------------------------
csenq
      dimension t(npt,1), sal(npt,1), pdens(npt,1),dens(npt,1),nzi(1),h(npt,1)
      include 'comm_new.h'

      call situ_dens (npt, nz, nzi, t, sal, dens, h) 
      call potn_dens (npt, nz, nzi, t, sal, pdens) 
c
c      do i = 1, npt
c         do k = 1, nzi(i)
c            POTND_BOT = amax1(POTND_BOT,pdens_pnt (t(i,k),sal(i,k)))
c         enddo
c      enddo
c      POTND_BOT = amax1(POTND_BOT,pdens_pnt (TEMP_BOT, SALT_BOT))


c      SITUD_BOT = POTND_BOT
c      do i = 1, npt
c         do k = 1, nzi(i)
c            SITUD_BOT = amax1(SITUD_BOT,dens(i,k))
c         enddo
c      enddo
c      SITUD_BOT = amax1(SITUD_BOT,sdens_pnt (TEMP_BOT, SALT_BOT, dep_max))


      write (iout, *) 'TEMP_BOT   = ', TEMP_BOT
      write (iout, *) 'SALT_BOT   = ', SALT_BOT
      write (iout, *) 'MAX_DEPTH  = ', dep_max
      write (iout, *) 'DENS_BOT(in situ)    = ', SITUD_BOT
      write (iout, *) 'DENS_BOT(potential)  = ', POTND_BOT

      return
      end

c----------------------------------------------------------------
      subroutine situ_dens (npt, nz, nzi, tem, sal, dens, h)
c----------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_para.h'
      include 'comm_new.h'
      dimension tem(npt,nz), sal(npt,nz), dens(npt,nz), h(npt,nz), nzi(npt)
c in situ SIGMA density as a function of potential temperature, 
c salinity & pressure. sigma == 1000. * (rho - 1.); 
c Senya Basin 1992
c
      if     (isalt .eq. 1) then
         do i = 1, npt
            do k = 1, nzi(i)
               dens(i,k) = SIGMA0 - TCOEF * (tem(i,k) - TEMP_BOT) 
            enddo
         enddo
      elseif (isalt .eq. 2) then
         do i = 1, npt
            do k = 1, nzi(i)
               dens(i,k) = pdens1 (tem(i,k))
            enddo
         enddo
      elseif (isalt .eq. 3) then
         do i = 1, npt
            pp  = h(i,1)/2.
            pp1 =  pp
            pp0 = -pp
            do k = 1, nzi(i)
               dens(i,k) = pdens4 (tem(i,k), sal(i,k), pp) 
               pp = pp0 + 2.*h(i,k)
               pp0 = pp1
               pp1 = pp
            enddo
         enddo
      elseif (isalt .eq. 4) then
         do i = 1, npt
            pp  = h(i,1)/2.
            pp1 =  pp
            pp0 = -pp
            do k = 1, nzi(i)
               dens(i,k) = pdens12 (tem(i,k), sal(i,k), pp) 
               pp = pp0 + 2.*h(i,k)
               pp0 = pp1
               pp1 = pp
            enddo
         enddo
      elseif (isalt .eq. 5) then
         do i = 1, npt
            pp  = h(i,1)/2.
            pp1 =  pp
            pp0 = -pp
            do k = 1, nzi(i)
               dens(i,k) = pdens17 (tem(i,k), sal(i,k), pp) 
               pp = pp0 + 2.*h(i,k)
               pp0 = pp1
               pp1 = pp
            enddo
         enddo
      elseif (isalt .eq. 6) then
         do i = 1, npt
            pp  = h(i,1)/2.
            pp1 =  pp
            pp0 = -pp
            do k = 1, nzi(i)
               salt = sal(i,k)
c Gilles variant:
               situ = theta_eos (pp, tem(i,k), salt, 0.)
               dens(i,k) = dens_unesco (situ, salt, pp)
               pp = pp0 + 2.*h(i,k)
               pp0 = pp1
               pp1 = pp
            enddo
         enddo
      elseif (isalt .eq. 7) then
         do i = 1, npt
            do k = 1, nzi(i)
               dens(i,k) = dens_marina (tem(i,k))
            enddo
         enddo
      endif

      return
      end

c----------------------------------------------------------------
      subroutine potn_dens (npt, nz, nzi, tem, sal, dens)
c----------------------------------------------------------------
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_para.h'
      include 'comm_new.h'
      dimension tem(npt,nz), sal(npt,nz), dens(npt,nz), nzi(npt)
c Potential SIGMA density as a function of potential temperature & salinity
c sigma == 1000. * (rho - 1.); 
c Senya Basin 1992
c
      if     (isalt .eq. 1) then
         do i = 1, npt
            do k = 1, nzi(i)
               dens(i,k) = SIGMA0 - TCOEF * (tem(i,k) - TEMP_BOT) 
            enddo
         enddo
      elseif (isalt .eq. 2) then
         do i = 1, npt
            do k = 1, nzi(i)
               dens(i,k) = pdens1 (tem(i,k))
            enddo
         enddo
      elseif (isalt .eq. 3) then
         do i = 1, npt
            do k = 1, nzi(i)
               dens(i,k) = pdens4 (tem(i,k), sal(i,k), 0.) 
            enddo
         enddo
      elseif (isalt .eq. 4) then
         do i = 1, npt
            do k = 1, nzi(i)
               dens(i,k) = pdens012 (tem(i,k), sal(i,k)) 
            enddo
         enddo
      elseif (isalt .eq. 5) then
         do i = 1, npt
            do k = 1, nzi(i)
               dens(i,k) = pdens017 (tem(i,k), sal(i,k)) 
            enddo
         enddo
      elseif (isalt .eq. 6) then
         do i = 1, npt
            do k = 1, nzi(i)
               dens(i,k) = dens_unesco (tem(i,k), sal(i,k), 0.)
            enddo
         enddo
      elseif (isalt .eq. 7) then
         do i = 1, npt
            do k = 1, nzi(i)
               dens(i,k) = dens_marina (tem(i,k))
            enddo
         enddo
      endif

      return
      end

c****************************************************************************
      subroutine dens_eos(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****************************************************************************
      function theta_eos(p0, t0, s, pf)
c****************************************************************************
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_eos = t + (xk - 2.0 * q) / 6.0

      return
      end
c****************************************************************************


c****************************************************************************
      subroutine sbulk(pr, t, s, kk)
c****************************************************************************
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****************************************************************************
      function atg(p, t, s)
c****************************************************************************
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

