c    @(#)advectqg.f	1.1
      subroutine advectqg
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
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)
c
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
      fact=1./dxyml
      fact2=fact/2.
c
      r3(1)=fact*(r1(2)-r1(1)+r2(mp1)-r2(1))
     .       -wsurf(1)*r3(1)/hdpth(1,mlcur)
c
      do 10 ip=2,mm1
        r3(ip)=fact2*(r1(ip+1)-r1(ip-1))+fact*(r2(ip+m)-r2(ip))
     .           -wsurf(ip)*r3(ip)/hdpth(ip,mlcur)
   10 continue
c
      r3(m)=fact*(r1(m)-r1(mm1)+r2(m+m)-r2(m))
     .        -wsurf(m)*r3(m)/hdpth(m,mlcur)
c
      do 20 j=2,nm1
        ip=(j-1)*m +1
        r3(ip)=fact*(r1(ip+1)-r1(ip))+fact2*(r2(ip+m)-r2(ip-m))
     .           -wsurf(ip)*r3(ip)/hdpth(ip,mlcur)
        do 15 i=2,mm1
          ip=ip+1
          r3(ip)=fact2*(r1(ip+1)-r1(ip-1)+r2(ip+m)-r2(ip-m))
     .             -wsurf(ip)*r3(ip)/hdpth(ip,mlcur)
   15   continue
        ip=ip+1
        r3(ip)=fact*(r1(ip)-r1(ip-1))+fact2*(r2(ip+m)-r2(ip-m))
     .           -wsurf(ip)*r3(ip)/hdpth(ip,mlcur)
   20 continue
c
      ip=msq-m+1
      r3(ip)=fact*(r1(ip+1)-r1(ip)+r2(ip)-r2(ip-m))
     .         -wsurf(ip)*r3(ip)/hdpth(ip,mlcur)
      do 30 i=2,mm1
        ip=ip+1
        r3(ip)=fact2*(r1(ip+1)-r1(ip-1))+fact*(r2(ip)-r2(ip-m))
     .           -wsurf(ip)*r3(ip)/hdpth(ip,mlcur)
   30 continue
      ip=msq
      r3(ip)=fact*(r1(ip)-r1(ip-1)+r2(ip)-r2(ip-m))
     .         -wsurf(ip)*r3(ip)/hdpth(ip,mlcur)

      end
