print*,'enter a temperature, salinity and depth'
      read*,temp,salt,dep

      eps = .001

      d0 = dens_unesco(temp,salt,dep)
      dt = dens_unesco(temp+eps,salt,dep)
      ds = dens_unesco(temp,salt+eps,dep)

      print*,d0,(dt-d0)/eps,(ds-d0)/eps

      stop
      end


      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, dens, h) 
c-----------------------------------------------------
csenq
      dimension t(npt,1), sal(npt,1), dens(npt,1), nzi(1), h(npt,1)
      include 'comm_new.h'

c
      do i = 1, npt
         do k = 1, nzi(i)
            POTND_BOT = amax1(POTND_BOT,pdens_pnt (t(i,k),sal(i,k)))
         enddo
      enddo
      POTND_BOT = amax1(POTND_BOT,pdens_pnt (TEMP_BOT, SALT_BOT))

      call situ_dens (npt, nz, nzi, t, sal, dens, h) 

      SITUD_BOT = POTND_BOT
      do i = 1, npt
         do k = 1, nzi(i)
            SITUD_BOT = amax1(SITUD_BOT,dens(i,k))
         enddo
      enddo
      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, 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,1), sal(npt,1), dens(npt,1), nzi(1)
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 dconv (npt,nz,nzi,u,v,uc,vc,h,t,sal,dens,tr,convn,zmld)
c------------------------------------------------------------------------------
c main aim: removing possible static instability in density.
c           No influence on layer depths in this version
c          (see isopycnal mixing model in Ragu's version)
c-----------------------------------------------------------------Senq Ltd. Co.
c
c      boundary condition at bottom changed - instead of t(i,nzi) = t_bot
c      i prefer t(i,nzi)_after = t(i,nzi)_before - NHN


      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      common /errors/ ioerr, nstep
c
      dimension u(npt,nz), v(npt,nz), uc(npt,nz), vc(npt,nz), h(npt,nz), 
     *        t(npt,nz), sal(npt,nz), dens(npt,nz), convn(npt,nz), nzi(npt),
     *        tr(npt,nz,1),zmld(npt)

      parameter (ERR_LEV = 1.e-4)

      do i=1,npt
         zbotint = 0.
         if_first = 1

         do k=1,nzi(i)-1
            kp = k+1
            zbotint = zbotint + h(i,k)

            rhok  = sdens_pnt(t(i,k) ,sal(i,k) ,zbotint)
            rhokp = sdens_pnt(t(i,kp),sal(i,kp),zbotint)

            if (rhok .gt. rhokp + ERR_LEV) then
               hinv = h(i,k) + h(i,kp)
               oldpotener = (dens(i,kp)*h(i,kp) + dens(i,k)*hinv)/2.

               hinv = 1.0 / hinv
               hki  = hinv * h(i,k)
               hkpi = hinv * h(i,kp)

               tmp = uc(i,k) + uc(i,kp)
               uc(i,k)  = tmp * hki
               uc(i,kp) = tmp * hkpi
               tmp = hinv * tmp
               u(i,k)  = tmp 
               u(i,kp) = tmp 

               tmp = vc(i,k) + vc(i,kp)
               vc(i,k)  = tmp * hki
               vc(i,kp) = tmp * hkpi
               tmp = hinv * tmp
               v(i,k)  = tmp 
               v(i,kp) = tmp                

               tmp1 = hki*t(i,k) + hkpi*t(i,kp)
               t(i,k)  = tmp1
               t(i,kp) = tmp1

               tmp = hki*sal(i,k) + hkpi*sal(i,kp)
               sal(i,k)  = tmp
               sal(i,kp) = tmp

               tmp = pdens_pnt (tmp1, tmp) 
               dens(i,k)  = tmp
               dens(i,kp) = tmp

               do m = 1, ntrac_tot
                  tmp = hki*tr(i,k,m) + hkpi*tr(i,kp,m)
                  tr(i,k,m)  = tmp
                  tr(i,kp,m) = tmp
               enddo

               znewpotener = dens(i,k)*(h(i,kp) + h(i,k)/2.)
               convn(i,k) =  oldpotener - znewpotener
            else
               if (if_first.eq.1) zmld(i) = max(zmld(i),zbotint)
               if_first = 0
            endif
         enddo
      enddo

      return
      end

c------------------------------------------------------------------------------
      subroutine dconv_cl (npt,nz,nzi,h,t,sal,dens)
c------------------------------------------------------------------------------
c main aim: removing possible static instability in density.
c           No influence on layer depths in this version
c          (see isopycnal mixing model in Ragu's version)
c-----------------------------------------------------------------Senq Ltd. Co.
c
c      boundary condition at bottom changed - instead of t(i,nzi) = t_bot
c      i prefer t(i,nzi)_after = t(i,nzi)_before - NHN

      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      common /errors/ ioerr, nstep
c
      dimension h(npt,nz), 
     *        t(npt,nz), sal(npt,nz), dens(npt,nz), nzi(npt)

      parameter (ERR_LEV = 1.e-4)
c

      do i = 1, npt
         do k = 1, nzi(i) - 1
            kp = k + 1

            if (dens(i,k)-dens(i,kp) .gt. ERR_LEV) then

               hinv = h(i,k) + h(i,kp)
               oldpotener = (dens(i,kp)*h(i,kp) + dens(i,k)*hinv)/2.

               hinv = 1.0 / hinv
               hki  = hinv * h(i,k)
               hkpi = hinv * h(i,kp)

               tmp1 = hki*t(i,k) + hkpi*t(i,kp)
               t(i,k)  = tmp1
               t(i,kp) = tmp1

               tmp = hki*sal(i,k) + hkpi*sal(i,kp)
               sal(i,k)  = tmp
               sal(i,kp) = tmp

               tmp = pdens_pnt (tmp1, tmp) 
               dens(i,k)  = tmp
               dens(i,kp) = tmp

               znewpotener = tmp*(h(i,kp) + h(i,k)/2.)
            endif
         enddo
      enddo

      return
      end

c---------------------------------------------------------------------
      subroutine drich_mix (npt, nz, nzi, h,u,v,uc,vc,tem,sal,tr, richn)
c---------------------------------------------------------------------
c Senya Basin, 1992

      implicit real(a-h,o-z),integer(i-n)
      include 'comm_para.h'
      include 'comm_new.h'

      dimension rnu(MAXNZ), rka(MAXNZ)
      dimension u(npt,nz), v(npt,nz), uc(npt,nz), vc(npt,nz), h(npt,nz), 
     *     tem(npt,nz), sal(npt,nz), nzi(npt), tr(npt,nz,1)
     *     ,richn(npt,nz)

      common /errors/ ioerr, nstep
      common /tria_loc/ ixy, ndim, alf(MAXNZ),bet(MAXNZ),gam(MAXNZ),aga(MAXNZ)

      parameter (R_COEF  = -0.5 * GRAVTY/1000.)
      parameter (R_CRIT  = 2.e5)
c      parameter (DUZ_0   = 1.e-5)
      parameter (DUZ_0   =  1.e-4)
      parameter (DRHO_0  = -1.e-6)
c
c  Ri = -(g/rho0) * d(rho)/dz / (du/dz**2 + du/dz**2)
c  
      do i = 1, npt
         ndim = nzi(i)
         rnu(ndim) = 0.
         rka(ndim) = 0.
         zbotint = 0.

         do k = 1, ndim-1

            kp = k + 1
            uu = u(i,k) - u(i,kp)
            vv = v(i,k) - v(i,kp)
            h12  = h(i,k) + h(i,kp)

            du2 = uu*uu + vv*vv
c            if (du2 .lt. DUZ_0) du2 = DUZ_0 
            du2 = du2 + DUZ_0

            zbotint = zbotint + h(i,k)
            rhom = sdens_pnt(tem(i,k ),sal(i,k ),zbotint)
            rhop = sdens_pnt(tem(i,kp),sal(i,kp),zbotint)

            drho = rhom - rhop 
            drho = min(drho, DRHO_0* h12)
            rich = R_COEF * h12 * drho / du2
            richn(i,k) = rich   ! save for output

            call visc_diff (rich, vnu, vka)
            tmp = DLT_MIX / h12
            rnu(k) = tmp * vnu
            rka(k) = tmp * vka
         enddo

         ixy = i
         call tria_init (npt, rnu, h)
         call tria_tem (npt, u, 0.)
         call tria_tem (npt, v, 0.)
            
         call tria_init (npt, rka, h)
         tb = tem(i,ndim)
         sb = sal(i,ndim)
         call tria_tem (npt, tem, tb)
         call tria_tem (npt, sal, sb)

         do m = 1, ntrac_tot
            trb = tr(i,ndim,m)
            call tria_tem (npt, tr(1,1,m), trb)
         enddo
         
      enddo
         
      do i = 1, npt 
         do k = 1, nzi(i)
            hi = h(i,k)
            uc(i,k) = u(i,k) * hi
            vc(i,k) = v(i,k) * hi
         enddo
      enddo
      
      return
      end

c---------------------------------------------
      subroutine visc_diff (Ri, rnu, rka)
c---------------------------------------------
c     eddy viscosity & diffusivity
c     a'la Pacanowski & Philander [1981]
c---------------------------------------------
c from PP-1981: 
c      parameter (GAMMA = 5., RNU_0 = 0.01, RNU_B = 5.e-5, RKA_B = 5.e-6)
c 
c suggested by Gilles to Senya (?):
      parameter (GAMMA = 5., RNU_0 = 0.05, RNU_B = 1.34e-5, RKA_B = 1.34e-6)

      parameter (RNU_NEG = RNU_B + RNU_0)
      parameter (RKA_NEG = RKA_B + RNU_NEG)

      if ( Ri .gt. 0. ) then
         tmp = 1. + GAMMA * Ri
         rnu = RNU_B + RNU_0 / (tmp * tmp)
         rka = RKA_B +   rnu / tmp
      else
         rnu = RNU_NEG
         rka = RKA_NEG
      endif

      return
      end

c------------------------------------------------------
      subroutine tria_init (npt, cappa, h)
c------------------------------------------------------
csenq
      include 'comm_para.h'
      real cappa(1), h(npt,1)
      common /tria_loc/ ixy, ndim, alf(MAXNZ),bet(MAXNZ),gam(MAXNZ),aga(MAXNZ)

      capk   = cappa(1)
      hi     = h(ixy,1)
      tmp    = capk + hi
      betk   = capk / tmp
      bet(1) = betk
      gam(1) = hi   / tmp

      do k = 2, ndim
         betkm1 = betk
         capkm1 = capk
         capk  = cappa(k) 
         hi    = h(ixy,k)

         tmp    = hi + capk + capkm1 - betkm1 * capkm1
         betk   = capk   / tmp
         bet(k) = betk
         aga(k) = capkm1 / tmp
         gam(k) = hi     / tmp
      enddo

      bet(ndim) = 0.

      return
      end 

      subroutine tria_run (npt, data)
c-------------------------------------
csenq
      include 'comm_para.h'

      real data(npt,1)
      common /tria_loc/ ixy, ndim, alf(MAXNZ),bet(MAXNZ),gam(MAXNZ),aga(MAXNZ)

      alfa = data(ixy,1)*gam(1)
      alf(1) = alfa

      do k = 2, ndim
         alfa   = data(ixy,k)*gam(k) + alfa*aga(k) 
         alf(k) = alfa
      enddo

      prev = 0.
      do k = ndim , 1, -1
         prev = alf(k) + bet(k) * prev
         data(ixy,k) = prev
      enddo

      return
      end

      subroutine tria_tem (npt, data, botval)
c--------------------------------------------
csenq
      include 'comm_para.h'

      real data(npt,1)
      common /tria_loc/ ixy, ndim, alf(MAXNZ),bet(MAXNZ),gam(MAXNZ),aga(MAXNZ)

      alfa = data(ixy,1)*gam(1)
      alf(1) = alfa

      do k = 2, ndim
         alfa   = data(ixy,k)*gam(k) + alfa*aga(k) 
         alf(k) = alfa
      enddo

      prev = botval
      do k = ndim, 1, -1
         prev = alf(k) + bet(k) * prev
         data(ixy,k) = prev
      enddo

      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 = 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


c     ------------------------------------------------------------------
      subroutine comp_bncy (npt,nzi,dens,bncy)
c     ------------------------------------------------------------------
c     Compute boyoncy as : b = -g(rho-rho_0)/rho_0
c     ------------------------------------------------------------------
      include 'comm_para.h'
      include 'comm_new.h'

      dimension dens(npt,1),bncy(npt,1),nzi(npt)
c
      coef = -GRAVTY/(1000. + POTND_BOT)
      do i = 1, npt
         do k = 1, nzi(i)
            bncy(i,k) = coef * (dens(i,k) - POTND_BOT)
         enddo
      enddo

      return
      end

c     ------------------------------------------------------------------
      subroutine cvmix(npt,nzi,h,t,s,b,u,v)
c     ------------------------------------------------------------------
c     Convective adjustment as of Dake Chen.
      include 'comm_new.h'

      dimension nzi(npt),h(npt,1),t(npt,1),s(npt,1),b(npt,1),u(npt,1),v(npt,1)
c
      do ns = 1, 2
         do ks = 1, 2
            do i = 1, npt
               do k = ks, nzi(i)-1, 2
                  k1 = k + 1
                  if (b(i,k) .lt. b(i,k1)) then
                     hik  = h(i,k)
                     hik1 = h(i,k1)
                     hsum1 = 1. / (hik + hik1)
                     t(i,k)  = (hik*t(i,k) + hik1*t(i,k1))*hsum1
                     t(i,k1) =  t(i,k)
                     s(i,k)  = (hik*s(i,k) + hik1*s(i,k1))*hsum1
                     s(i,k1) =  s(i,k)
                     b(i,k)  = (hik*b(i,k) + hik1*b(i,k1))*hsum1
                     b(i,k1) =  b(i,k)
                  endif
               enddo
            enddo
         enddo
      enddo
      
      return
      end

c     ----------------------------------------------------------------
      function tke0(wp,b0,br,h,hp)
c     ----------------------------------------------------------------
        hp2 = hp + hp
	hexp = h - hp2 + (h + hp2)*exp(-h/hp)
        bp   = b0*h + br*hexp
	tke0 = wp - bp
      return
      end

c     ---------------------------------------------------------------------
      subroutine ktmix(npt,nsig,ddt,h,t,s,b,u,v,q,qr,ep,taux,tauy,sigma,dh1)
c     ---------------------------------------------------------------------
c     Vertical Mixing Using Kraus-Turner Scheme.(Dake Chen, 1995)

      include 'comm_new.h'

      dimension u(npt,1),v(npt,1),h(npt,1),t(npt,1),s(npt,1),b(npt,1),
     +          dh1(1),q(1),qr(1),ep(1),taux(1),tauy(1),sigma(1)
      data alph/2.55e-4/, beta/7.6e-4/, gravty/9.8/, taumin/3.e-5/, hp/17./
c
      ga = alph*gravty
      gb = beta*gravty
c
      cm2 = 2.0 * cm_mix
      cn2 = cn_mix / 2.0

      do i = 1, npt
         h10 = h(i,1)
         tau = sqrt(taux(i)**2 + tauy(i)**2)
         tau = amax1(tau,taumin)
         ustar = sqrt(tau)
         wp = cm2 * ustar * tau
         br = ga*qr(i)
         b0  = ga*q(i) - gb*ep(i)
         dbh = (b(i,1) - b(i,2))*h10

         b0 = b0 - cn2 * (b0-abs(b0))
         dbh = amax1(dbh,1.e-5)

         tke = tke0(wp,b0,br,h10,hp)
         if( tke .lt. 0.) then
            h1 = h10
            h2 = 0.5*h1
            f1 = tke0(wp,b0,br,h1,hp)
            do iter = 1, 10
               f2   = tke0(wp,b0,br,h2,hp)
               hnew = h2 - f2*(h2-h1)/(f2-f1)
               err  = abs(hnew-h2)/h2
               if (err .lt. 1.e-4) goto 15
               h1 = h2
               h2 = hnew
               f1 = f2
            enddo
   15       continue
            hnew = amin1(h10, hnew)
         else
            h1 = h10
            h2 = h10 + h10
            hm = h10
            f1 = -ddt*tke0(wp,b0,br,hm,hp)
            do iter = 1, 10
               hm   = 0.5*(h2+h10)
               f2   = dbh*(h2-h10) - ddt*tke0(wp,b0,br,hm,hp)
               hnew = h2 - f2*(h2-h1)/(f2-f1)
               err  = abs(hnew-h2)/h2
               if (err .lt. 1.e-4) goto 25
               h1 = h2
               h2 = hnew
               f1 = f2
            enddo
   25       continue
            hnew = amax1(h10, hnew)
         endif

         hnew = amax1(hmin_mix, hnew)
         hnew = amin1(hmax_mix, hnew)
         dh     = hnew - h10
         adh    = amin1(abs(dh), h(i,2))
         dh1(i) = sign(adh, dh)
         
      enddo

      call impmix (npt,nsig,dh1,u,v,h,t,s,sigma)

      return
      end

c     ----------------------------------------------------------------
      subroutine impmix (npt,nz,dh1,u,v,h,t,s,sigma)
c     ----------------------------------------------------------------
c     Adjust vertical profiles according to the depth change of the 
c     mixed layer as of Dake Chen. /modified by Senya Basin, July 1995/

      include 'comm_para.h'
      include 'comm_new.h'

      dimension dh1(1),u(npt,1),v(npt,1),h(npt,1),
     *     t(npt,1),s(npt,1), sigma(1),
     *     u1(MAXNZ),u2(MAXNZ),v1(MAXNZ),v2(MAXNZ),t1(MAXNZ),t2(MAXNZ),
     *     aa(MAXNZ),bb(MAXNZ),cc(MAXNZ),dh(MAXNZ),we(MAXNZ),hr(MAXNZ),
     *     s1(MAXNZ),s2(MAXNZ)

      do i = 1, npt
         dh_mix = dh1(i)

         if (dh_mix .ne. 0.) then

            dh(1) = dh_mix
            do k = 2, nz
               dh(k) = -1.5*(sigma(k)+sigma(k+1))*dh_mix
            enddo
            
            we(1) = dh_mix
            do k = 2, nz-1
               we(k) = we(k-1) + dh(k)
            enddo

            do k = 1, nz
               hik = h(i,k)
               u1(k)  = u(i,k)*hik
               v1(k)  = v(i,k)*hik
               t1(k)  = t(i,k)*hik
               s1(k)  = s(i,k)*hik
               h(i,k) = hik + dh(k)
            enddo

           do k = 1, nz-1
              hr(k) = we(k)/(h(i,k)+h(i,k+1))
           enddo

           hr1 = 0.5*(dh_mix - abs(dh_mix)) / h(i,1)
           hr2 = 0.5*(dh_mix + abs(dh_mix)) / h(i,2)
           cc(1) = -hr2
           bb(1) = 1. - hr1
           aa(2) = hr1
           cc(2) = -hr(2)
           bb(2) = 1. - hr(2) + hr2
           aa(nz) = hr(nz-1)
           bb(nz) = 1. + aa(nz)

           hr_k = hr(2)
           do k = 3, nz-1
              hr_km1 = hr_k
              hr_k   = hr(k)
              aa(k) =  hr_km1
              cc(k) = -hr_k
              bb(k) = 1. + hr_km1 - hr_k
           enddo

           call tridiag(aa,bb,cc,u1,u2,nz)
           call tridiag(aa,bb,cc,v1,v2,nz)
           call tridiag(aa,bb,cc,t1,t2,nz)

           call tridiag(aa,bb,cc,s1,s2,nz)
           do k = 1, nz
              hik = 1./h(i,k)
              u(i,k) = hik*u2(k)
              v(i,k) = hik*v2(k)
              t(i,k) = hik*t2(k)
              s(i,k) = hik*s2(k)
           enddo
        endif
      enddo  

      return
      end

c     ------------------------------------------------------------------
      subroutine jpmix (npt,nz,nzi,h,t,s,b,u,v)
c     ------------------------------------------------------------------
c     Reduce gradient Richardson # instability using a Jim Price criterion.

      include 'comm_para.h'
      include 'comm_new.h'
      dimension u(npt,1),v(npt,1),h(npt,1),t(npt,1),s(npt,1),b(npt,1),
     *          gama(MAXNZ), ric(MAXNZ), nzi(npt)
      logical use_gamma

      save EPSILON, use_gamma, ifirst, ric, ric_inv

      data EPSLON/1.e-9/, ifirst/0/
c
      if ( ifirst .eq. 0 ) then
         ifirst = 1
         gama(1)    = gam1_mix
         ric(1)     = ric1_mix
         do k = 2, nz
            gama(k)    = gam2_mix
            ric(k)     = ric2_mix
         enddo
         use_gamma = (iuse_gam .eq. 1)
      endif

      do kn = 1, 2
         do ks = 1, 2
            do i = 1, npt
               do k = ks, nzi(i)-1, 2
                  
                  Rik      = ric(k)
                  Rik_inv  = 1./ric(k)
                  gamkt    = gama(k)
                  gamkv    = 0.5*(1.+gama(k))
                  
c     compute the gradient Richardson number
                  hm = h(i,k) + h(i,k+1)
                  bd = b(i,k) - b(i,k+1)
                  dd = 0.5 * hm * bd
                  
                  ud = u(i,k+1) - u(i,k)  
                  vd = v(i,k+1) - v(i,k) 
                  dv = ud*ud + vd*vd
                  Ri = amax1(dd,0.)/amax1(dv,EPSLON)
c     check to see if Ri is critical or not
                  
                  if (Ri .lt. Rik) then
c     partially mix layers k and k+1
                     if ( use_gamma ) then
                        ri2 = Ri*Rik_inv
                        ct = 1.- ri2**gamkt
                        cv = 1.- ri2**gamkv
                     else
                        ct = 1. - Ri*Rik_inv
                        cv = ct
                     endif
                     c1 = h(i,k)   / hm
                     c2 = h(i,k+1) / hm
                     td = (t(i,k+1)-t(i,k))*ct
                     sd = (s(i,k+1)-s(i,k))*ct
                     bd = -bd*ct
                     ud = ud*cv
                     vd = vd*cv
                     
                     t(i,k+1) = t(i,k+1) - td*c1
                     t(i,k)   = t(i,k)   + td*c2
                     
                     s(i,k+1) = s(i,k+1) - sd*c1
                     s(i,k)   = s(i,k)   + sd*c2
                     
                     b(i,k+1) = b(i,k+1) - bd*c1
                     b(i,k)   = b(i,k)   + bd*c2
                     
                     u(i,k+1) = u(i,k+1) - ud*c1
                     u(i,k)   = u(i,k)   + ud*c2
                     v(i,k+1) = v(i,k+1) - vd*c1
                     v(i,k)   = v(i,k)   + vd*c2
                  endif
               enddo
            enddo
         enddo
      enddo

      return
      end


c----------------------------------------------------------------------------
      subroutine tridiag(A,B,C,Y,X,N)
c----------------------------------------------------------------------------

c	Modified 2 Feb 1991

c	This routine solves a matrix equation of the form MX=Y where Y is
c	the know vector and M is an NxN tridiagonal matrix with a diagonal of
c	B, a lower diagonal of A and an upper diagonal of C.  It should be
c	noted that N will be equal to the number of layers (nz).  If there
c	are more than 100 the array size will exceed memory allocation.
c
c	This routine requires the matrix elements (A,B,C), the known vector
c	(Y) and the dimensions of both (N).
c
c	This routine supplies the new vector X.
c	The variables are defined as:
c
c	    Name	       Type		Description
c	    ====	       ====		===========
c
c	a(n)                  From      Lower diagonal elements
c			      dimpl 
c	b(n)		      From      Main diagonal elements
c			      dimpl
c	bit		      Internal  Storage variable
c	c(n)                  From      Upper diagonal elements
c			      dimpl   
c	gam(nmax)	      Internal  Storage array
c	nmax                  Internal  Maximum array size
c	x(n)                  Ret to    new velocity (temp) field
c			      dimpl
c	y(n)		      From      Original velocity (temp) field
c			      dimpl
c

c------------------ Establish variables and constants -----------------------
      parameter (nmax=100)
      dimension gam(nmax),a(n),b(n),c(n),y(n),x(n)

      bet  = b(1)
      x(1) = y(1)/bet
c-------------------------- Forward substitute ------------------------------
      do j = 2, n
         gam(j) = c(j-1)/bet
         bet    = b(j) - a(j)*gam(j)
         x(j)   = (y(j) - a(j)*x(j-1)) / bet
      enddo
c-------------------------- Back substitute ------------------------------
      do j = n-1, 1, -1
         x(j) = x(j) - gam(j+1)*x(j+1)
      enddo

      return
      end

c-----------------------------------------------------------
      subroutine windmix(npt,nz,dtmix,limp,delt,windm0,windz0,
     &              nzi,u,v,uc,vc,h,t,sal,dens,tr,convn,wnd,zmld)
c-----------------------------------------------------------
c  Kraus-Turner mixing AND convective adjustment (a la dconv)
c         for use with a level grid;  includes the effect of
c         wind stirring
c
c  Uses the Kraus-Turner vertical mixing algorithm described in:
c         Stefan Rahmstorf, 1991: A zonal-averaged model of the 
c         ocean's response to climate change, JGR Oceans, V.96,
c         p. 6951-6963;
c
      implicit real(a-h,o-z),integer(i-n)
      include 'comm_new.h'
      include 'comm_para.h'
c      include 'comm_tracer.h'
      common /errors/ ioerr, nstep
      real m,m0
      dimension u(npt,nz),v(npt,nz),uc(npt,nz),vc(npt,nz),
     &          h(npt,nz),t(npt,nz),sal(npt,nz),dens(npt,nz),
     &          convn(npt,nz),nzi(npt),tr(npt,nz,1),wnd(npt,3),zmld(npt)


      CD = 1.3e-3
      cd3 = CD**1.5
c                          windm0 from input parameter file;
      do i=1,npt
         zbotint = 0.
         tmix = t(i,1)
         smix = sal(i,1)

         wspeed = wnd(i,1)
         wspeed3 = wspeed*wspeed*wspeed
         ustar3  = wspeed3*cd3

         if_first = 1

         do k=1,nzi(i)-1
            kp = k+1
            zbotint = zbotint + h(i,k)
            rhomix = sdens_pnt(tmix,smix,zbotint)
            rhok =   sdens_pnt(t(i,kp),sal(i,kp),zbotint)
c  compute buoyancy [isn't this really potential energy since
c                    I'm multiplying buoyancy by depth???]
            buoy = GRAVTY/(1000.+SITUD_BOT)*(rhok-rhomix)*zbotint
c  compute wind energy

            m = windm0*exp(-zbotint/windz0)

            winden  = m*ustar3*(delt*limp*86400.)/zbotint

c  check for mixing

            if (winden .ge. buoy) then
               hinv = h(i,k) + h(i,kp)
               oldpotener = (dens(i,kp)*h(i,kp) + dens(i,k)*hinv)/2.

               hinv = 1.0 / hinv
               hki  = hinv * h(i,k)
               hkpi = hinv * h(i,kp)

               tmp = uc(i,k) + uc(i,kp)
               uc(i,k)  = tmp * hki
               uc(i,kp) = tmp * hkpi
               tmp = hinv * tmp
               u(i,k)  = tmp 
               u(i,kp) = tmp 

               tmp = vc(i,k) + vc(i,kp)
               vc(i,k)  = tmp * hki
               vc(i,kp) = tmp * hkpi
               tmp = hinv * tmp
               v(i,k)  = tmp 
               v(i,kp) = tmp                

               tmp1 = hki*t(i,k) + hkpi*t(i,kp)
               t(i,k)  = tmp1
               t(i,kp) = tmp1

               tmp = hki*sal(i,k) + hkpi*sal(i,kp)
               sal(i,k)  = tmp
               sal(i,kp) = tmp

               tmp = pdens_pnt (tmp1, tmp) 
               dens(i,k)  = tmp
               dens(i,kp) = tmp

               do m = 1, ntrac_tot
                  tmp = hki*tr(i,k,m) + hkpi*tr(i,kp,m)
                  tr(i,k,m)  = tmp
                  tr(i,kp,m) = tmp
               enddo

               znewpotener = dens(i,k)*(h(i,kp) + h(i,k)/2.)
               convn(i,k) =  oldpotener - znewpotener

               tmix = t(i,k)
               smix = sal(i,k)
            else
               if (if_first.eq.1) zmld(i) = max(zmld(i),zbotint)
               ustar3  = 0.
               if_first = 0
            endif
         enddo
      enddo

      return
      end