c    @(#)helm4.f	1.1
c
c
c
      subroutine helm4(klevel,xlam)
c
c --- subroutine to solve the poisson equation.
c
c --- this is a fourth-order poisson solver. to maintain
c --- accuracy for the vorticity equation, the solution of the
c --- poisson equation for the streamfunction must also be of
c --- fourth order. this is done using the method of deferred
c --- corrections. see haidvogel,schulman, and robinson p 12.
c
      include 'qg_par.inc'
      common/scrat/r(mn0),r1(mn0),r2(mn0),r3(mn0)
      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 /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/edges/icorn(4),ncorn(4),mast(4),nast(4)
      hx2=hy*hy
      cd=8.-xlam*hx2
      ce=xlam*xlam*hx2/12.
      one12=1./12
      cc2=1./(6.*hx2)
c
c --- call ncar poisson solver package
c
      dbasin=hy*nm1
      call pwscrt(0,0.,xbasin,mm1,1,bda,bdb,0.,dbasin,nm1,1,bdc,bdd
     +,xlam,stream(1,klevel),m,1,ierror,r1)
c
c --- set up the streamfunction array for the interior solution
c
      do 20 j=2,nm1
        ip=(j-1)*m+1
        do 10 i=2,mm1
          ip=ip+1
          a1=one12*(cd*r(ip)+r(ip-1)+r(ip+1)+r(ip+m)+r(ip-m))
          a2=cc2*(4.*stream(ip,klevel)+stream(ip-mp1,klevel)+
     1         stream(ip-mm1,klevel)+stream(ip+mm1,klevel)+
     2           stream(ip+mp1,klevel)
     3             -2.*(stream(ip-m,klevel)+stream(ip-1,klevel)+
     4               stream(ip+1,klevel)+stream(ip+m,klevel) ) )
          r1(ip)=a1-a2 +ce*stream(ip,klevel)
   10   continue
   20 continue
c
c --- place the interior of r1 into stream
c
      do 40 i=2,nm1
	ip=(i-1)*m+1
        do 30 j=2,mm1
	ip=ip+1
          stream(ip,klevel)=r1(ip)
   30   continue
   40 continue
c
c --- call ncar poisson solver package for second iteration
c
      call pwscrt(0,0.,xbasin,mm1,1,bda,bdb,0.,dbasin,nm1,1,bdc,bdd
     +,xlam,stream(1,klevel),m,1,ierror,r1)
      return
      end
