c    @(#)loop3d.f	1.1
c
c
c
      subroutine loop3d
c
c --- subroutine loop3d is the main computation loop for this
c --- code.  it advances the vorticity and solves for the new
c --- streamfunction fields.
c
      include 'qg_par.inc'
      common/scrat/r(mn0),r1(mn0),r2(mn0),r3(mn0)
      common/cnn/n,nm1,nm2,np1
      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/fields/vort(mn0,k0),stream(mn0,k0),prvjac(mn0,k0),
     +  topden(mn0),botden(mn0),relief(mn0),prvtop(mn0),
     +  prvbot(mn0),topfct(k0),botfct(k0),pvort(mn0)
      common/edges/icorn(4),ncorn(4),mast(4),nast(4)
      common/threed/hz(k0),scsig,sigz(k01),eigval(k0),ainv(k0,k0),
     +      amat(k0,k0)
c
c ---advance the vorticity
c
      do 100 k=1,kz
        call vorteq(k)
  100 continue
c
c --- solve for three-dimensional streamfunction
c --- get boundary values for stream function.
c
      do 105 k=1,kz
        call infld(stream(1,k),k,tim,1,1)
  105 continue
c
c --- now transform them to ainv*stream function using r as temp. storage.
c
      do 150 j=1,2
	jm=0
	if(j.eq.2)jm=msq-m
	jmp1=jm+1
	jmpm=jm+m
        do 140 ip=jmp1,jmpm
          do 120 k=1,kz
            r(k)=0.
            do 110 k1=1,kz
              r(k)=r(k)+ainv(k,k1)*stream(ip,k1)
  110       continue
  120     continue
          do 130 k=1,kz
            stream(ip,k)=r(k)
  130     continue
  140   continue
  150 continue
c
      do 200 j=1,2
	jpm=j+m
	if(j.eq.2)jpm=m+m
	j1=1
	if(j.eq.2)j1=m
	jpmm=j1+nm2*m
        do 190 ip=jpm,jpmm,m
          do 170 k=1,kz
            r(k)=0.
            do 160 k1=1,kz
              r(k)=r(k)+ainv(k,k1)*stream(ip,k1)
  160       continue
  170     continue
          do 180 k=1,kz
            stream(ip,k)=r(k)
  180     continue
  190   continue
  200 continue
c
c --- solve helmholtz equation level by level.
c
      do 300 k=1,kz
c
c --- first use r to hold the transformed vorticity field for this level.
c
        do 250 ip=1,msq
          r(ip)=0.0
          if(iftop.ne.0)goto 211
          if(ifbot.ne.0)goto 221
          do 210 k1=1,kz
            r(ip)=r(ip)+ainv(k,k1)*vort(ip,k1)
  210     continue
          goto 250
  211     continue
          if(ifbot.ne.0)goto 231
          do 220 k1=1,kz
            r(ip)=r(ip)+ainv(k,k1)*(vort(ip,k1)-topfct(k1)*
     +                                  topden(ip))
  220     continue
          goto 250
  221     continue
          do 230 k1=1,kz
            r(ip)=r(ip)+ainv(k,k1)*(vort(ip,k1)+botfct(k1)*
     +                                  botden(ip))
  230     continue
          goto 250
  231     continue
          do 240 k1=1,kz
            r(ip)=r(ip)+ainv(k,k1)*(vort(ip,k1)-
     +            topfct(k1)*topden(ip)+botfct(k1)*botden(ip))
  240     continue
  250   continue
c
c --- now store the interior of the field in the interior of the stream fun
c
        mmm=nm2*m
        do 290 j=m,mmm,m
          jp2=j+2
          jpmm1=j+mm1
          do 280 ip=jp2,jpmm1
            stream(ip,k)=r(ip)
  280     continue
  290   continue
c
c --- call subroutine helm4 to solve the poisson equation
c
        call helm4(k,eigval(k))
c
  300 continue
c
c --- new streamfunction--point by point
c
      do 350 ip=1,msq
        do 340 k=1,kz
           r(k)=stream(ip,k)
           stream(ip,k)=0.
  340   continue
c
        do 345 k=1,kz
        do 345 j=1,kz
           stream(ip,k)=stream(ip,k)+amat(k,j)*r(j)
  345   continue
  350 continue
      return
      end
