c
c    @(#)advect3d.f	1.6
c
c
c
c
      subroutine advect3d(f,fadd,fnew,base)
c
c
c leonard walstad 5/84
c
c
c this routine will advect level k of field f using u,v,wz
c as defined the the corresponding arrays.  the result is stored 
c in af in units of f/s.
c
c----------------------------------------------------------------------
c
c
      include'mladv.inc'
c
c
c
      dimension f(mnmax,lmax), fadd(mnmax,lmax), fnew(mnmax,lmax)
      dimension adv3d(lmax), vadj(lmax)
      external base
c
c
c  qg model quantities needed  ( the arrays in common /scrat/ are
c  used as scratch arrays.)
c
      common /scrat/ r(mn0), r1(mn0), r2(mn0), r3(mn0)
      dimension r4(mn0,lmax)
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
      common /cnn/n,nm1,nm2,np1
c
      common /extpar/ rlat0, rlng0, t0, v0, dhor, ht, time0,
     +                r0,f0star

c
c
c----------------------------------------------------------------------
c
c   first horizontally advect by the quasigeostrophic fieldss.......
c   this is accomplished by using the arakawa jacobian .....
c
c   darajac is the usual arakawa jacobian using sf to scale and 
c   not integrated, therefore the weights are changed on the 
c   sides and corners.
c
c
      sf=v0/(hy*dxyml)

      do 1 k=1,l

         if ((k.gt.kmlmin).and.(k.le.kmlmax)) then
            call advfix(k,1)
           else
            call darajac(r,f(1,k),r1,m,n,sf)
          endif

         if (ifadvqg .eq. 1) then
            do 2 ip=1,msq
               r4(ip,k) = r1(ip)
 2           continue
          else
            do 3 ip=1,msq
               r4(ip,k) = 0.0
 3           continue
         endif

 1     continue

c
c --- now advect
c

      do 80 ip=1,msq

      ifvadj = 0

      do 70 k=1,l

        wtot = wsurf(ip) + wzvel(ip) * zml(k) + wbg(ip,k)

        if (k.ge.kmlmin) then
          tmp=0.0
        else
          tmp=1.0
        endif

        if (ifadvml .eq. 0) then
           tmp = 0.0
         endif

c
c --- vertical derivatives
c
      if (ifvadv .eq. 0) then
         c0 = 0.0
         c1 = 0.0
         c2 = 0.0
         goto 78
       endif
 
      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 (mlstep * wtot .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) = fadd(ip,k) - fadd(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 = (fadd(ip,k-1) * dz1 + fadd(ip,k) * dz2 ) / dz3
                vadj(k) = fadd(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 = (fadd(ip,k+1) * dz1 + fadd(ip,k) * dz2 ) / dz3
                vadj(k) = fadd(ip,k) - vadj1

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

              else
                vadj(k) = 0.0

             endif

          endif

c
c  compute vertical advections
c

      if (ifvadv.eq.0) then

         vadv = 0.0

        else

         vadv = wtot * (c0*f(ip,k00)+c1*f(ip,k1)+c2*f(ip,k2))

       endif

      adv3d(k) = - mlstep*(vadv + r4(ip,k) + tmp * r3(ip))
c
c  end loop over all levels k
c
   70 continue

c
c  add advection terms and adjust profile
c

         call findkml(hdpth(ip,mlcur),zml,l,k1)
         call findkml(hdpth(ip,mladd),zml,l,k2)

         if (k1 .ne. k2) then

c  mixed layer
 
             do 115 k=1,k2
                if (k .gt. k1) then
                   kadd = k1
                  else
                   kadd = k
                 endif
                fnew(ip,k) = fadd(ip,k) + adv3d(kadd)
 115          continue

c  remnant layer

             do 116 k=k2+1,l
                if (k .gt. k1) then
                   kadd = k
                  else
                   kadd = k1 + 1
                 endif
                fnew(ip,k) = fadd(ip,k) + adv3d(kadd)
 116          continue

           else

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

          endif

      if (ifvadj.eq.1) then

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

       endif

      if (mlstep * wtot .gt. 0.0) fnew(ip,l)=base(ip)

 80   continue

      return
      end
