c    
c @(#)sink.f	1.2
c
c
      subroutine sink(f,fadd,fnew)
c
c
c
      include'mladv.inc'
      include'trc.inc'
c
c
c
      dimension f(mnmax,lmax+1),fadd(mnmax,lmax+1),fnew(mnmax,lmax+1)
c
      common/parm/kz,m,mm1,mm2,mp1,msq,tim,tmax,tstart,hy,rea,dt,
     +   xbasin,yy(mx0),icalc,kzp1,kzm1,iftop,ifbot,iftvv,ifrel,theta,
     +   ifsbl,ifeva
c
c
c
      isave = icurtrace

      do 80 ip=1,msq

      if (wsinkfac(isave) .ne. 0.0) then

         fave = 0.0

         do 85 k=1,l

            icurtrace = 1
            s  = trace(ip,k,mlcur,icurtrace) + trcqgofk(ip,k)
            icurtrace = 4
            s2 = trace(ip,k,mlcur,icurtrace) + trcqgofk(ip,k)

            q1 = (s * exp(-kn * s2)) / (ks + s)
            q2 = s2 / (ks2 + s2)
            fave = fave + (q1 / (q1 + q2))

  85      continue

         fave = fave / float(l)

        endif

      sinkvel = wsink(isave) + fave * wsinkfac(isave)

      call findkml(hdpth(ip,mlcur),zml,lmax,kml)

      do 70 k=kml+1,l+1

c
c --- vertical derivatives
c

      if (k .eq. l) then

         k00 = k
         k1  = k-1
         k2  = k

         if (mlstep * sinkvel .gt. 0.0) then
            c0 = 1.0/dzml(l)
            c1 = 0.0
            c2 = 0.0
           else
            c2 = 0.0
            c1 = 1.0 / (zml(k1) - zml(k00))
            c0 = -(c2+c1)
          endif

       else if (k .eq. l+1) then

         k00 = l+1
         k1  = l
         k2  = l+1

         if (mlstep * sinkvel .gt. 0.0) then
            c0 = 0.0
            c1 = 0.0
            c2 = 0.0
           else
            c0 = 0.0
            c1 = 1.0/dzml(l)
            c2 = 0.0
          endif

       else

         k00 = k
         k1  = k-1
         k2  = k+1

         if (mlstep * sinkvel .gt. 0.0) then
             c1 = 0.0
             c2 = 1.0 / (zml(k2) - zml(k00))
           else
             c2 = 0.0
             c1 = 1.0 / (zml(k1) - zml(k00))
          endif
                 
        endif

      if (k .lt. l) then
         c0=-(c2+c1)
       endif

c
c --- now advect
c

        fnew(ip,k) = fadd(ip,k) - mlstep * sinkvel *
     .                 (c0*f(ip,k00)+c1*f(ip,k1)+c2*f(ip,k2))

   70  continue

         k00 = kml
         k1  = kml+1
         k2  = kml+1
 
         if (mlstep * sinkvel .gt. 0.0) then
             c0 = 0.0
             c1 = -1.0 / dzml(k1)
             c2 = 0.0
           else
             c0 = -1.0 / dzml(k00)
             c1 = 0.0
             c2 = 0.0
          endif

        fnew(ip,kml) = fadd(ip,kml) - mlstep * sinkvel *
     .                 (c0*f(ip,k00)+c1*f(ip,k1)+c2*f(ip,k2))


   80  continue
c
      icurtrace = isave
c
      return
      end
