c    @(#)vorteq.f	1.1
c
c
c
      subroutine vorteq(klevel)
c
c --- subroutine vorteq performs the calulation of the vorticity field
c
c --- in the finite-element approximation, the order of the calculation
c --- of the vorticity is the following: interior points, regular
c --- boundary points, and corner points. afterwards, filter the field.
c
c --- array vort - vorticity
c --- arrays r,r1 - scratch arrays
c --- variable bet represents the beta effect
c
c --- fq is an arakawa jacobian which represents a finite-difference
c --- approximation. fq is anti-symmetric ith its arguments and
c --- conserves vorticity, energy, and enstrophy when integrated
c --- over a closed domain.
c
      include 'qg_par.inc'
      common/scrat/r(mn0),r1(mn0),r2(mn0),r3(mn0)
	common /cnn/n,nm1,nm2,np1
      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/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/philtr/iltord,iltfrq,iltcnt
      common/dynam/alpha,beta,rf
      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)
      data icnt/0/
c
c --- computation of rhs of finite-element vorticity eq.
c
c --- first the arakawa jacobian
c
      call arajac(stream(1,klevel),vort(1,klevel),r,m,n)
c
c --- normalize the jacobian by alpha
c
      do 20 ip=1,msq
        r(ip)=alpha*r(ip)
   20 continue
c
c --- include beta effect
c
      if(beta.eq.0)goto 40
      call arajac(stream(1,klevel),pvort,r2,m,n)
      do 30 ip=1,msq
        r(ip)=r(ip)+r2(ip)
   30 continue
   40 continue
c
c --- apply bottom friction. if no bottom friction wanted or if
c --- not on bottom level, skip this section.
c
      if(klevel.ne.kz)go to 65
      if(rf.eq.0.0)go to 65
c
      cf=(rf/36.0)*(hy**2)
      do 45 iside=1,4
        ip=icorn(iside)
        istp=mast(iside)
        il=ip-istp
        ir=ip+istp
        jstp=nast(iside)
        iul=il+jstp
        iu=ip+jstp
        iur=ir+jstp
c
c --- corner friction term.  mult by app. computat. molecule.
c
        r(ip)=r(ip)+cf*(vort(iur,klevel) + 2 * (vort(iu,klevel) +
     +        vort(ir,klevel)) + 4 * vort(ip,klevel))
c
c --- now regular boundary points.
c
	if(istp.eq.1.or.istp.eq.-1)mn1=mm1
	if(istp.eq.m.or.istp.eq.-m)mn1=nm1
        do 44 i=2,mn1
          ip=ip+istp
          il=il+istp
          ir=ir+istp
          iul=iul+istp
          iu=iu+istp
          iur=iur+istp
          r(ip)=r(ip)+cf*(vort(iul,klevel) + vort(iur,klevel) +
     +          2*(vort(il,klevel) + vort(ir,klevel)) +
     +          4*vort(iu,klevel) + 8*vort(ip,klevel))
   44   continue
   45 continue
c
c --- now interior points.
c
      do 50 irow=2,mn1
        ip1=(irow-1)*m+2
        ip2=(irow*m)-1
        do 49 ip=ip1,ip2
          r(ip)=r(ip) + cf*(vort(ip-m-1,klevel) + vort(ip-m+1,klevel) +
     +          vort(ip+m-1,klevel) + vort(ip+m+1,klevel) + 
     +          4*(vort(ip-m,klevel) + vort(ip+1,klevel) +
     +          vort(ip+m,klevel) + vort(ip-1,klevel)) + 
     +          16*vort(ip,klevel))
   49   continue
   50 continue
c
   65 continue
c
c --- now convert this term to the proper form for msolve
c --- and genrl to solve the adams-bashforth problem.
c
      call rhs(r,prvjac(1,klevel))
c
c --- if this is the initialization step, return.
c
      if(icalc.eq.0) go to 700
c
c --- call subroutines msolve, infld, genrl to
c --- do the inversion of the mass matrix
c
c --- msolve - matrix inversion for the interior field
c --- infld  - fill a scratch array with boundary values
c              used by subroutine genrl
c --- genrl  - invert the tri-diagonal system for boundary values
c
      call msolve
      call infld(r1,klevel,tim,2,1)
      call genrl(stream(1,klevel),vort(1,klevel),klevel)
c
c --- update the vorticity field
c
      do 70 ip=1,msq
        vort(ip,klevel)=vort(ip,klevel)+r(ip)
   70 continue
c
c --- definitions of the filter variables:
c
c --- icnt = counter (for frequency of application)
c --- iltcnt = frequency with which filter is applied
c --- iltfrq = no of times filter is applied per time step
c --- iltord = order of the filter
c
      icnt=icnt+1
      if(icnt.ne.iltcnt)go to 120
      if(iltord.eq.0) go to 120
      icnt=0
      do 110 i=1,iltfrq
        call filter(vort(1,klevel),m,n,iltord)
  110 continue
  120 continue
c
c --- come here (skipping the update of vorticity) if icalc .eq. 0
c
  700 continue
c
c --- if iftop=0 then do not want top density
c
      if(iftop.eq.0)goto 850
c
c --- section for prognostic computation of surface
c --- temperature anomaly.
c
c --- this calculation is done with the top layer only.
c
      if(klevel.ne.1)goto 850
c
c --- compute jacobian of top density
c
      call arajac(stream(1,1),topden,r,m,n)
      if(iftvv.ne.0)goto 720
      do 710 ip=1,msq
        r(ip)=alpha*r(ip)
  710 continue
      go to 740
c
c --- compute  -  top vertical velocity and add to alpha*r
c --- according to the mass matrix, storing the result in r
c
  720 continue
      call infld(r1,1,tim,5,0)
c
c --- now the top vertical velocity (w1) is in r1, so multiply by 
c --- the mass matrix and add to alpha*r(ip) which is alpha*
c --- jacobian (psi,top density anomaly)
c
      cv=(hy**2)/36
      do 725 iside=1,4
        ip=icorn(iside)
        istp=mast(iside)
        il=ip-istp
        ir=ip+istp
        jstp=nast(iside)
        iul=il+jstp
        iu=ip+jstp
        iur=ir+jstp
c
c --- corner velocity term.  mult by app. computat. molecule.
c
        r(ip)=alpha*r(ip)+cv*(r1(iur) + 2 * (r1(iu) +
     +        r1(ir)) + 4 * r1(ip))
c
c --- now regular boundary points.
c
	if(istp.eq.1.or.istp.eq.-1)mn1=mm1
	if(istp.eq.m.or.istp.eq.-m)mn1=nm1
        do 724 i=2,mn1
          ip=ip+istp
          il=il+istp
          ir=ir+istp
          iul=iul+istp
          iu=iu+istp
          iur=iur+istp
          r(ip)=alpha*r(ip)+cv*(r1(iul) + r1(iur) +
     +          2*(r1(il) + r1(ir)) +
     +          4*r1(iu) + 8*r1(ip))
  724   continue
  725 continue
c
c --- now interior points.
c
      do 730 irow=2,nm1
        ip1=(irow-1)*m+2
        ip2=(irow*m)-1
        do 729 ip=ip1,ip2
          r(ip)=alpha*r(ip) + cv*(r1(ip-m-1) + r1(ip-m+1) +
     +          r1(ip+m-1) + r1(ip+m+1) + 
     +          4*(r1(ip-m) + r1(ip+1) +
     +          r1(ip+m) + r1(ip-1)) + 
     +          16*r1(ip))
  729   continue
  730 continue
c
c
  740 continue
      call rhs(r,prvtop)
c
c --- if initialization, then finished with top density
c
      if(icalc.eq.0)goto 850
c
c --- now convert this term to the proper form for msolve
c --- and genrl to solve the adams-bashforth problem.
c
      call msolve
      call infld(r1,1,tim,7,1)
      call genrl(stream(1,1),topden,1)
c
c --- now advance the top density
c
      do 820 ip=1,msq
        topden(ip)=topden(ip)+r(ip)
  820 continue
      itcnt=itcnt+1
      if(itcnt.ne.iltcnt)goto 850
      if(iltord.eq.0)goto 850
      itcnt=0
      do 830 i=1,iltfrq
        call filter(topden,m,n,iltord)
  830 continue
c
c
c --- finished with top density calculations
c
  850 continue
c
c --- bottom density calculations are performed on the
c --- bottom layer only.
c
      if(ifbot.eq.0)goto 900
      if(klevel.ne.kz)goto 900
c
c --- form jacobian of alpha*botden + eta
c
      if(ifrel.eq.1)goto 865
      do 860 ip=1,msq
        r1(ip)=alpha*botden(ip)
  860 continue
      goto 875
  865 continue
      do 870 ip=1,msq
        r1(ip)=alpha*botden(ip)+relief(ip)
  870 continue
  875 continue
c
      fact=-hz(kz)/scsig/sigz(kz+1)/2.0
      do 877 ip=1,msq
        r2(ip)=stream(ip,kz)+fact*botden(ip)
  877 continue
      call arajac(r2,r1,r,m,n)
c
c --- now convert this term to the proper form for msolve
c --- and genrl to solve the adams-bashforth problem.
c
      call rhs(r,prvbot)
c
c --- if initialization, the done with bottom density
c
      if(icalc.eq.0)goto 900
c
c --- solve the equation
c
      call msolve
      call infld(r1,kz,tim,8,1)
      call genrl(r2,botden,0)
c
c --- advance the bottom density
c
      do 880 ip=1,msq
        botden(ip)=botden(ip)+r(ip)
  880 continue
c
c --- filter the bottom density
c
      ibcnt=ibcnt+1
      if(ibcnt.ne.iltcnt)goto 895
      if(iltord.eq.0)goto 895
      ibcnt=0
      do 890 i=1,iltfrq
        call filter(botden,m,n,iltord)
  890 continue
c
  895 continue
c
  900 continue
c
      return
      end
