c    @(#)getfld.f	1.1
c
c
c
      subroutine getfld(anapsi,anazeta,anatop,time)
c
c
      include'qg_par.inc'
      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
      common/field/psi(mn0,k0),zeta(mn0,k0),tlast,day1,tkeep,
     +  curpsi(mn0,k0),curzet(mn0,k0),psiprv(mn0,k0),zetprv(mn0,k0),
     +  topprv(mn0),topcur(mn0),topdens(mn0)
      common/unit/iin,iout,idff,idif,iorr,iowr,iplt,ifndd
      common/cnn/n,nm1,nm2,np1
      common/threed/hz(k0),scsig,sigz(k01),eigval(k0),ainv(k0,k0),
     +     amat(k0,k0)
      common/user/pram(16)
      dimension anapsi(m20,k0),anazeta(m20,k0),anatop(m20)
      dimension therm(m20,k0),dxxpsi(m20),dyypsi(m20)
      ifvort=pram(2)
      iftherm=pram(3)
c
        read(33)mbound,nbound,day1
         print 5,mbound,nbound,day1
    5    format( ' mbound,nbound,day1 = ',2i10,f10.1)
        day1=day1+0.03
        nsq=mbound*nbound
        do 10 k=1,kz
          read(33)(anapsi(i,k),i=1,nsq)
   10   continue
        if(ifvort.eq.0) then
          do 20 k=1,kz
            read(33)(anazeta(i,k),i=1,nsq)
   20     continue
        endif
       if(iftop.ne.0) read(33)(anatop(i),i=1,nsq)
       np=nbound
       mp=mbound 
      write(iout,111)mbound,nbound,time/day1
  111 format(2x,'infld called for m,n,time=',2i3,f5.2)
c
       if(iftherm.eq.1) call thermvt(anapsi,therm,kz,nsq)
       do 500 k=1,kz
       if(ifvort.eq.1) then
       call arelvor(anapsi(1,k),dxxpsi,dyypsi,mp,np,hy)
       call filter(dxxpsi,mp,np,4)
       call filter(dyypsi,mp,np,4)
       do 400 ip=1,nsq
 400   anazeta(ip,k)=dxxpsi(ip)+dyypsi(ip)
       endif  
       do 600 ip=1,nsq
       anazeta(ip,k)=anazeta(ip,k)+iftherm*therm(ip,k)
       if(k.eq.1.and.iftop.ne.0) anazeta(ip,k)=anazeta(ip,k)+
     1    iftherm*anatop(ip)/hz(1)
600    continue
500    continue
      if(day1.lt.0.05) rewind 33
      if(day1.lt.0.05) day1=9999.9
c
      return
      end
