c
c    @(#)advect3d.f	1.6
c
      subroutine advbg(f,base)

      include'mladv.inc'

      dimension f(mnmax,lmax)
      dimension adv3d(lmax), vadj(lmax)
      external base

      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
                

      do 80 ip=1,msq

      ifvadj = 0

      do 70 k=1,l

        curz=zml(k)
        acurz=abs(curz)

c
c --- vertical derivatives
c
      if(k.eq.1) then
        k00=1
        k1 = 1
        k2= 2
c -- assume surface vertical velocity is zero
        c1 = 0
        c2 = zml(1)/(zml(2) - zml(1))/zml(2)
      elseif(k.eq.l) then
        k00=k
        k1=k-1
        k2=k
c -- one sided difference
        c1 = 1/(zml(l-1)-zml(l))
        c2 = 0
      else

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

        if (ifvadv.eq.1) then

           c1=(zml(k2)-zml(k00))/(zml(k1)-zml(k00))/(zml(k2)-zml(k1))
           c2=-(zml(k1)-zml(k00))/(zml(k2)-zml(k00))/(zml(k2)-zml(k1))

         else if (ifvadv.ge.2) then

c    upwind differencing
           if (ifvadv .eq. 2) then
              kl = -1
             else
              kl = kmlofq(hdpth(ip,mlcur),zml,l)
            endif

           if (wbg(ip,k) .gt. 0.0) then

c    upwelling
              c1 = 0.0
              if (k00.eq.kl-1) then
                 c2 = 0.0
                else if (k00.eq.kl) then
                 zl = zml(k00) - 0.5*dzml(k00)
                 zm = zl + 0.5*(-zl - hdpth(ip,mlcur))
                 c2 = 1.0 / (zml(k2) - zm)
                 c2 = 1.0 / (zml(k2) - zml(k00))
                else
                 c2 = 1.0 / (zml(k2) - zml(k00))
               endif

            else
c    downwelling
              c2 = 0.0
              if (k00.eq.kl) then
                 c1 = 0.0
                else if (k1.eq.kl) then
                 zl = zml(k1) - 0.5*dzml(k1)
                 zm = zl + 0.5*(-zl - hdpth(ip,mlcur))
                 c1 = 1.0 / (zm - zml(k00))
                 c1 = 1.0 / (zml(k1) - zml(k00))
                else
                 c1 = 1.0 / (zml(k1) - zml(k00))
               endif

            endif

          else
              write(*,*)'Invalid ifadv'
         endif

      endif
c

      c0=-(c2+c1)

78    continue

        if  ( (ifvadv.ne.0) .and. (zcross(ip,1).ne.0.0) ) then

             ifvadj = 1
             kold = int(zcross(ip,2))

             if ( (zcross(ip,1).eq.-1.0) .and. (k.eq.kold+1) ) then

                vadj(k) = f(ip,k) - f(ip,k-1)

              else if ((zcross(ip,1).eq.-1.0) .and. (k.eq.kold+2)) then
                dz1 = 0.5 * dzml(kold+1) - (hdpth(ip,mlcur) 
     .                                   + zml(kold+1))
                dz2 = dzml(kold+2)
                dz3 = dz1 + dz2
                vadj1 = (f(ip,k-1) * dz1 + f(ip,k) * dz2 ) / dz3
                vadj(k) = f(ip,k) - vadj1

              else if ( (zcross(ip,1).eq.1.0) .and. (k.eq.kold-1) ) then
                dz1 = 0.5 * dzml(kold) + (hdpth(ip,mlcur) + zml(kold))
                dz2 = dzml(kold - 1)
                dz3 = dz1 + dz2
                vadj1 = (f(ip,k+1) * dz1 + f(ip,k) * dz2 ) / dz3
                vadj(k) = f(ip,k) - vadj1

              else if ( (zcross(ip,1).eq.1.0) .and. (k.eq.kold) ) then
                vadj(k) = f(ip,k) - f(ip,k+1)

              else
                vadj(k) = 0.0

             endif

          endif

      adv3d(k) = -qgts * wbg(ip,k) * 
     .                     (c0*f(ip,k00)+c1*f(ip,k1)+c2*f(ip,k2))

   70 continue

c
c  add advection terms and adjust profile
c

      do 120 k=1,l
         f(ip,k) = f(ip,k) + adv3d(k)
 120   continue

      if (ifvadj.eq.1) then

         do 130 k=1,l
            if (vadj(k).ne.0.0) then
               f(ip,k) = f(ip,k) - vadj(k) - adv3d(k)
             endif
 130      continue

       endif

      if (wbg(ip,k) .gt. 0.0) f(ip,l) = base(ip)

 80   continue

      return
      end
