cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine adjmd2cr(nx,ny,phi,cf)
c     ****     
c     ****     adjust righthand side in cf(i,j,10) for boundary 
c     ****     conditions
c     ****     
      implicit real*8 (a-h,o-z)
      integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,
     +             maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,
     +             kcycle,iprer,ipost,intpol,kps
      integer nx,ny,i,j
      dimension cf(nx,ny,10),phi(0:nx+1,0:ny+1)
      common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,
     +               iguess, maxcy,method,nwork,lwork,itero,ngrid,
     +               klevel,kcur,kcycle,iprer,ipost,intpol,kps
      common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
c     ****     
c     ****     set specified boundaries in rhs from phi
c     ****     
      if (nxa.eq.1) then
	i = 1
	do j=1,ny
	  cf(i,j,10) = phi(i,j)
	end do
      end if
      if (nxb.eq.1) then
	i = nx
	do j=1,ny
	  cf(i,j,10) = phi(i,j)
	end do
      end if
      if (nyc.eq.1) then
	j = 1
	do i=1,nx
	  cf(i,j,10) = phi(i,j)
	end do
      end if
      if (nyd.eq.1) then
	j = ny
	do i=1,nx
	  cf(i,j,10) = phi(i,j)
	end do
      end if
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine cor2(nx,ny,phif,ncx,ncy,phic,nxa,nxb,nyc,nyd,intpol,
     +                phcor)
c     ****     
c     ****     add coarse grid correction in phic to fine grid 
c     ****     approximation in phif using linear or cubic 
c     ****     interpolation
c     ****     
      implicit real*8 (a-h,o-z)
      integer i,j,nx,ny,ncx,ncy,nxa,nxb,nyc,nyd,intpol,ist,ifn,jst,jfn
      dimension phif(0:nx+1,0:ny+1),phic(0:ncx+1,0:ncy+1)
      dimension phcor(0:nx+1,0:ny+1)
      do j=0,ny+1
	do i=0,nx+1
	  phcor(i,j) = 0.0
	end do
      end do
c     ****     
c     ****     lift correction in phic to fine grid in phcor
c     ****     
      call prolon2(ncx,ncy,phic,nx,ny,phcor,nxa,nxb,nyc,nyd,intpol)
c     ****     
c     ****     add correction in phcor to phif on nonspecified 
c     ****     boundaries
c     ****     
      ist = 1
      ifn = nx
      jst = 1
      jfn = ny
      if (nxa.eq.1) ist = 2
      if (nxb.eq.1) ifn = nx-1
      if (nyc.eq.1) jst = 2
      if (nyd.eq.1) jfn = ny-1
      do j=jst,jfn
	do i=ist,ifn
	  phif(i,j) = phif(i,j) + phcor(i,j)
	end do
      end do
c
c     add periodic points if necessary
c
      if (nyc.eq.0) then
	do i=ist,ifn
	  phif(i,0) = phif(i,ny-1)
	  phif(i,ny+1) = phif(i,2)
	end do
      end if
      if (nxa.eq.0) then
	do j=jst,jfn
	  phif(0,j) = phif(nx-1,j)
	  phif(nx+1,j) = phif(2,j)
	end do
      end if
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine dismd2cr(nx,ny,cf,tx,ty,wk,ier)
      INCLUDE "param.h" 
c     ****     
c     ****     discretize elliptic pde for mud2cr, set nonfatal errors
c     ****     
      integer nx,ny,i,j,l,im1,jm1,ier,nnx,nny
      dimension cf(nx,ny,10),tx(nx,ny,*),ty(ny,nx,*)
      dimension wk(*)
      integer        intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,
     +               iguess, maxcy,method,nwork,lwork,itero,ngrid,
     +               klevel,kcur,kcycle,iprer,ipost,intpol,kps
      common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,
     +               iguess, maxcy,method,nwork,lwork,itero,ngrid,
     +               klevel,kcur,kcycle,iprer,ipost,intpol,kps
      common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
      include "eceee.h"
c     ****     
c     ****     CHECK FOR CONSISTENCYT WRT KLEVEL
c     ****     
      NNX = ixp*2**(KLEVEL-1)+1
      NNY = jyq*2**(KLEVEL-1)+1
      IF(NNX.NE.NX.OR.NNY.NE.NY)THEN
	WRITE(6,100)NX,NY,NNX,NNY,ixp,jyq,KLEVEL
  100   FORMAT(' INCONSISTENCY WRT LEVEL. NX,NY,NNX,NNY,ixp,jyq,
     +  klevel = ',7I5)
	STOP
      ENDIF
      CALL SEEE(CEE(NC(6-KLEVEL)),NX,NY,CF)
c     ****     
c     ****     set coefficient for specified boundaries
c     ****     
      if (nxa.eq.1) then
	i = 1
	do j=1,ny
	  do l=1,9
	    cf(i,j,l) = 0.0
	  end do
	  cf(i,j,9) = 1.0
	end do
      end if
      if (nxb.eq.1) then
	i = nx
	do j=1,ny
	  do l=1,9
	    cf(i,j,l) = 0.0
	  end do
	  cf(i,j,9) = 1.0
	end do
      end if
      if (nyc.eq.1) then
	j = 1
	do i=1,nx
	  do l=1,9
	    cf(i,j,l) = 0.0
	  end do
	  cf(i,j,9) = 1.0
	end do
      end if
      if (nyd.eq.1) then
	j = ny
	do i=1,nx
	  do l=1,9
	    cf(i,j,l) = 0.0
	  end do
	  cf(i,j,9) = 1.0
	end do
      end if
c     ****     
c     ****     set and factor tridiagonal matrices for line 
c     ****     relaxation(s) if flagged
c     ****     
      if (method.eq.1.or.method.eq.3) then
	if (nxa.ne.0) then
c         ****     
c         ****     nonperiodic x line relaxation
c         ****     
	  do i=1,nx
	    im1 = max0(i-1,1)
	    do j=1,ny
	      tx(im1,j,1) = cf(i,j,5)
	      tx(i,j,2) = cf(i,j,9)
	      tx(i,j,3) = cf(i,j,1)
	    end do
	  end do
	  call factri(ny,nx,tx(1,1,1),tx(1,1,2),tx(1,1,3))
	else
c         ****     
c         ****    periodic x line relaxation
c         ****     
	  if (nx .gt. 3) then
c           ****     
c           ****     set and factor iff nx > 3
c           ****     
	    do i=1,nx-1
	      do j=1,ny
		tx(i,j,1) = cf(i,j,5)
		tx(i,j,2) = cf(i,j,9)
		tx(i,j,3) = cf(i,j,1)
	      end do
	    end do
	    call factrp(ny,nx,tx,tx(1,1,2),tx(1,1,3),tx(1,1,4),
     +                  tx(1,1,5),wk(kps))
	  end if
	end if
      end if
      if (method.eq.2.or.method.eq.3) then
	if (nyc.ne.0) then
c         ****     
c         ****     nonperiodic y line relaxation
c         ****     
	  do j=1,ny
	    jm1 = max0(j-1,1)
	    do i=1,nx
	      ty(jm1,i,1) = cf(i,j,7)
	      ty(j,i,2) = cf(i,j,9)
	      ty(j,i,3) = cf(i,j,3)
	    end do
	  end do
	  call factri(nx,ny,ty(1,1,1),ty(1,1,2),ty(1,1,3))
	else
c         ****     
c         ****     periodic y line relaxation
c         ****     
	  if (ny .gt. 3) then
c           ****     
c           ****     set and factor iff ny > 3
c           ****     
	    do j=1,ny-1
	      do i=1,nx
		ty(j,i,1) = cf(i,j,7)
		ty(j,i,2) = cf(i,j,9)
		ty(j,i,3) = cf(i,j,3)
	      end do
	    end do
	    call factrp(nx,ny,ty,ty(1,1,2),ty(1,1,3),ty(1,1,4),
     +                  ty(1,1,5),wk(kps))
	  end if
	end if
      end if
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine factri(m,n,a,b,c)
c     ****     
c     ****     factor the m simultaneous tridiagonal systems of order n
c     ****     
      implicit real*8 (a-h,o-z)
      integer m,n,i,j
      dimension a(n,m),b(n,m),c(n,m)
      do i=2,n
	do j=1,m
	  a(i-1,j) = a(i-1,j)/b(i-1,j)
	  b(i,j) = b(i,j)-a(i-1,j)*c(i-1,j)
       end do
      end do
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine factrp(m,n,a,b,c,d,e,sum)
c     ****     
c     ****     factor the m simultaneous "tridiagonal" systems of order 
c     ****     n from discretized periodic system (leave out periodic n 
c     ****     point) (so sweeps below only go from i=1,2,...,n-1) 
c     ****     n > 3 is necessary
c     ****     
      implicit real*8 (a-h,o-z)
      integer m,n,i,j
      dimension a(n,m),b(n,m),c(n,m),d(n,m),e(n,m),sum(m)
      do j=1,m
	d(1,j) = a(1,j)
      end do
      do i=2,n-2
	do j=1,m
	  a(i,j) = a(i,j)/b(i-1,j)
	  b(i,j) = b(i,j)-a(i,j)*c(i-1,j)
	  d(i,j) = -a(i,j)*d(i-1,j)
       end do
      end do
c     ****     
c     ****     correct computation of last d element
c     ****     
      do j=1,m
	d(n-2,j) = c(n-2,j)+d(n-2,j)
      end do
      do j=1,m
	e(1,j) = c(n-1,j)/b(1,j)
      end do
      do i=2,n-3
	do j=1,m
	  e(i,j) = -e(i-1,j)*c(i-1,j)/b(i,j)
	end do
      end do
      do j=1,m
	e(n-2,j) = (a(n-1,j)-e(n-3,j)*c(n-3,j))/b(n-2,j)
      end do
c     ****     
c     ****     compute  inner product (e,d) for each j in sum(j)
c     ****     
      do j=1,m
	sum(j) = 0.
      end do
      do i=1,n-2
	do j=1,m
	  sum(j) = sum(j)+e(i,j)*d(i,j)
	end do
      end do
c     ****     
c     ****     set last diagonal element
c     ****     
      do j=1,m
	b(n-1,j) = b(n-1,j)-sum(j)
      end do
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine kcymd2cr(wk)
c     ****     
c     ****     execute multigrid k cycle from kcur grid level
c     ****     kcycle=1 for v cycles, kcycle=2 for w cycles
c     ****     
      implicit real*8 (a-h,o-z)
      dimension wk(*)
      integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,
     +             maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,
     +             kcycle,iprer,ipost,intpol,kps
      integer nx,ny,ip,ic,ipc,irc,itx,ity,ncx,ncy,l,nrel
      integer kpbgn,kcbgn,ktxbgn,ktybgn,nxk,nyk,isx,jsy
      common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,
     +               iguess, maxcy,method,nwork,lwork,itero,ngrid,
     +               klevel,kcur,kcycle,iprer,ipost,intpol,kps
      common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
      common/mud2crc/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50),
     +nxk(50),nyk(50),isx,jsy
      integer kount(50)
      klevel = kcur
      nx = nxk(klevel)
      ny = nyk(klevel)
      ip = kpbgn(klevel)
      ic = kcbgn(klevel)
      itx = ktxbgn(klevel)
      ity = ktybgn(klevel)
c     ****     
c     ****     prerelax at current finest grid level
c     ****     
      do l=1,iprer
        call relmd2cr(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps))
      end do
      if (kcur .eq. 1) go to 5
c     ****     
c     ****     restrict residual to kcur-1 level
c     ****     
      ipc = kpbgn(klevel-1)
      ncx = nxk(klevel-1)
      ncy = nyk(klevel-1)
      irc = kcbgn(klevel-1)+9*ncx*ncy
      call resmd2cr(nx,ny,wk(ip),ncx,ncy,wk(ipc),wk(irc),wk(ic),wk(kps))
c     ****     
c     ****     set counter for grid levels to zero
c     ****     
      do l = 1,kcur
	kount(l) = 0
      end do
c     ****     
c     ****     set new grid level and continue k-cycling
c     ****     
      klevel = kcur-1
      nrel = iprer
c     ****     
c     ****     kcycle control point
c     ****     
   10 continue
c     ****     
c     ****     post relax when kcur revisited
c     ****     
      if (klevel .eq. kcur) go to 5
c     ****     
c     ****     count hit at current level
c     ****     
      kount(klevel) = kount(klevel)+1
c     ****     
c     ****     relax at current level
c     ****     
      nx = nxk(klevel)
      ny = nyk(klevel)
      ip = kpbgn(klevel)
      ic = kcbgn(klevel)
      itx = ktxbgn(klevel)
      ity = ktybgn(klevel)
      do l=1,nrel
        call relmd2cr(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps))
      end do
      if (kount(klevel) .eq. kcycle+1) then
c       ****     
c       ****     kcycle complete at klevel
c       ****     
	ipc = ip
	ip = kpbgn(klevel+1)
	ncx = nxk(klevel)
	ncy = nyk(klevel)
	nx = nxk(klevel+1)
	ny = nyk(klevel+1)
c       ****     
c       ****     inject correction to finer grid
c       ****     
	call cor2(nx,ny,wk(ip),ncx,ncy,wk(ipc),nxa,nxb,nyc,nyd,
     +            intpol,wk(kps))
c       ****     
c       ****     reset counter to zero
c       ****     
	kount(klevel) = 0
c       ****     
c       ****     ascend to next higher level and set to postrelax there
c       ****     
	klevel = klevel+1
	nrel = ipost
	go to 10
      else
	if (klevel .gt. 1) then
c         ****     
c         ****     kcycle not complete so descend unless at coarsest 
c         ****     grid
c         ****     
	  ipc = kpbgn(klevel-1)
	  ncx = nxk(klevel-1)
	  ncy = nyk(klevel-1)
	  irc = kcbgn(klevel-1)+9*ncx*ncy
           call resmd2cr(nx,ny,wk(ip),ncx,ncy,wk(ipc),wk(irc),wk(ic),
     +                   wk(kps))
c         ****     
c         ****     prerelax at next coarser level
c         ****     
	  klevel = klevel-1
	  nrel = iprer
	  go to 10
	else
c         ****     
c         ****     postrelax at coarsest level
c         ****     
	  do l=1,ipost
    	    call relmd2cr(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps))
	  end do
	  ipc = ip
	  ip = kpbgn(2)
	  ncx = nxk(1)
	  ncy = nyk(1)
	  nx = nxk(2)
	  ny = nyk(2)
c         ****     
c         ****     inject correction to level 2
c         ****     
          call cor2(nx,ny,wk(ip),ncx,ncy,wk(ipc),nxa,nxb,nyc,nyd,
     +              intpol,wk(kps))
c         ****     
c         ****     set to postrelax at level 2
c         ****     
	  nrel = ipost
	  klevel = 2
	  go to 10
	end if
      end if
    5 continue
c     ****     
c     ****     post relax at current finest grid level
c     ****     
      nx = nxk(kcur)
      ny = nyk(kcur)
      ip = kpbgn(kcur)
      ic = kcbgn(kcur)
      itx = ktxbgn(kcur)
      ity = ktybgn(kcur)
      do l=1,ipost
        call relmd2cr(nx,ny,wk(ip),wk(ic),wk(itx),wk(ity),wk(kps))
      end do
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine mud(pe,jntl)
      implicit real*8 (a-h,o-z)
      integer jntl
c     ****     
c     ****     set grid size params
c     ****     
      integer iixp,jjyq,iiex,jjey,nnx,nny,llwork
      parameter (iixp = 5 , jjyq = 3, iiex = 9, jjey = 9 )
      parameter (nnx=iixp*2**(iiex-1)+1, nny=jjyq*2**(jjey-1)+1)
      integer imx0, jmx0, imx1, jmx1,imx2,jmx2,imx3,jmx3,imx4,jmx4
      PARAMETER(IMX0=NNX,JMX0=NNY,IMX1=(IMX0+1)/2,JMX1=(JMX0+1)/2,
     1  IMX2=(IMX1+1)/2,JMX2=(JMX1+1)/2,IMX3=(IMX2+1)/2,
     2  JMX3=(JMX2+1)/2,IMX4=(IMX3+1)/2,JMX4=(JMX3+1)/2)
      integer nc0,nc1,nc2,nc3,nc4,ncee,nny2
      PARAMETER (NC0=1,NC1=NC0+10*(IMX0*JMX0),NC2=NC1+9*(IMX1*JMX1),
     1  NC3=NC2+9*(IMX2*JMX2),NC4=NC3+9*(IMX3*JMX3),
     2  NCEE=NC4+9*(IMX4*JMX4)-1)
      PARAMETER (NNY2=2*NNY+1)
      include "eceee.h"
c     ****     
c     ****     estimate work space for point relaxation (see mud2cr.d)
c     ****     
      parameter (llwork=(7*(nnx+2)*(nny+2)+76*nnx*nny)/3 )
      dimension phi(nnx,nny),rhs(nnx,nny),work(llwork)
c     ****     
c     ****     put integer and floating point argument names in 
c     ****     contiguous storage for labelling in vectors iprm,fprm
c     ****     
      integer iprm(17),mgopt(4)
      dimension fprm(6)
      integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nx,ny,
     +              iguess,maxcy,method,nwork,lwrkqd,itero
      common/itmud2cr/
     + intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nx,ny,
     +              iguess,maxcy,method,nwork,lwrkqd,itero
      common/ftmud2cr/
     + xa,xb,yc,yd,tolmax,relmax
      equivalence(intl,iprm)
      equivalence(xa,fprm)
      integer i,j,ierror
      dimension PE(NNX,1)
      integer maxcya
      DATA MAXCYA/150/
      integer mm,nn,jj,jjj
c     ****     
c     ****     set input integer arguments
c     ****     
      MM = NNX
      NN = NNY
      PI = 4*DATAN(1.D0)
C     ****     
C     ****     SET INPUT INTEGER PARAMETERS
C     ****     
      INTL = JNTL
c     ****     
c     ****     set boundary condition flags
c     ****     
      nxa = 0
      nxb = 0
      nyc = 2
      nyd = 1
c     ****     
c     ****     set grid sizes from parameter statements
c     ****     
      ixp = iixp
      jyq = jjyq
      iex = iiex
      jey = jjey
      nx = nnx
      ny = nny
c     ****     
c     ****     set multigrid arguments (w(2,1) cycling with fully 
c     ****     weighted residual restriction and cubic prolongation)
c     ****     
      mgopt(1) = 2
      mgopt(2) = 2
      mgopt(3) = 1
      mgopt(4) = 3
c     ****     
c     ****     set for one cycle
c     ****     
      maxcy = maxcya
c     ****     
c     ****     set no initial guess forcing full multigrid cycling
c     ****     
      iguess = 0
c     ****     
c     ****     set work space length approximation from parameter 
c     ****     statement
c     ****     
      nwork = llwork
c     ****     
c     ****     set line z relaxation
c     ****     
      method = 3
c     ****     
c     ****     set end points of solution rectangle in (x,y) space
c     ****     
      xa = -pi
      xb =  pi
      yc = 0.0
      yd = 0.5*pi
c     ****     
c     ****     set error control flag
c     ****     
      tolmax = 0.01
c     ****     
c     ****     set right hand side in rhs
c     ****     initialize phi to zero
c     ****     
      do i=1,nx
	do j=1,ny
	  RHS(I,J) = CEE(I+(J-1)*NX+9*NX*NY)
	  phi(i,j) = 0.0
	end do
      end do
c     ****     
c     ****     set specified boundaries in phi
c     ****     
      DO I=1,NX
        PHI(I,NY) = RHS(I,NY)/CEE(I+(NY-1)*NX+8*NX*NY)
      END DO
c     ****     
c     ****     intialization call
c     ****     
       call mud2cr(iprm,fprm,work,rhs,phi,mgopt,ierror)
c     ****     
c     ****     attempt solution
c     ****     
      intl = 1
      call mud2cr(iprm,fprm,work,rhs,phi,mgopt,ierror)
C     ****     
C     ****     COPY PHI TO PE
C     ****     
      DO J = 1,NY
	JJ = NY+J-1
	JJJ = NY+1-J
	DO I = 1,NX
	  PE(I,JJ) = PHI(I,J)
	  PE(I,JJJ) = PHI(I,J)
	END DO
      END DO
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine mud2cr(iparm,fparm,work,rhs,phi,mgopt,
     +                  ierror)
      implicit real*8 (a-h,o-z)
      integer iparm,mgopt,ierror
      integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,
     +             maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,
     +             kcycle,iprer,ipost,intpol,kps
      integer kpbgn,kcbgn,ktxbgn,ktybgn,nxk,nyk,isx,jsy
      integer int,iw,k,kb,nx,ny,ic,itx,ity
      dimension iparm(17),fparm(6),mgopt(4)
      dimension phi(*),rhs(*),
     +           work(*)
      common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,
     +               iguess, maxcy,method,nwork,lwork,itero,ngrid,
     +               klevel,kcur,kcycle,iprer,ipost,intpol,kps
      common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
      common/mud2crc/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50),
     +nxk(50),nyk(50),isx,jsy
      data int / 0 /
      save int
      ierror = 1
      intl = iparm(1)    ! set and check intl on all calls
      if (intl*(intl-1).ne.0) return
      if (int.eq.0) then
	int = 1
	if (intl.ne.0) return  ! very first call is not intl=0
      end if
      ierror = 0
c     ****     
c     ****     set  arguments internally
c     ****     these will not be rechecked if intl=1!
c     ****     
      nxa = iparm(2)
      nxb = iparm(3)
      nyc = iparm(4)
      nyd = iparm(5)
      ixp = iparm(6)
      jyq = iparm(7)
      iex = iparm(8)
      jey = iparm(9)
      ngrid = max0(iex,jey)
      nfx = iparm(10)
      nfy = iparm(11)
      iguess = iparm(12)
      maxcy = iparm(13)
      method = iparm(14)
      nwork = iparm(15)
      kcycle = mgopt(1)
      if (kcycle .eq. 0) then
c       ****     
c       ****     set defaults
c       ****     
	kcycle = 2
	iprer = 2
	ipost = 1
	intpol = 3
      else
	iprer = mgopt(2)
	ipost = mgopt(3)
	intpol = mgopt(4)
      end if
      xa = fparm(1)
      xb = fparm(2)
      yc = fparm(3)
      yd = fparm(4)
      tolmax = fparm(5)
      if (intl .eq. 0) then  ! intialization call
c       ****     
c       ****     check input arguments
c       ****     
	ierror = 2   ! check boundary condition flags
	if (max0(nxa,nxb,nyc,nyd).gt.2) return
	if (min0(nxa,nxb,nyc,nyd).lt.0) return
	if (nxa.eq.0.and.nxb.ne.0) return
	if (nxa.ne.0.and.nxb.eq.0) return
	if (nyc.eq.0.and.nyd.ne.0) return
	if (nyc.ne.0.and.nyd.eq.0) return
	ierror = 3   ! check grid sizes
	if (ixp.lt.2) return
	if (jyq.lt.2) return
	ierror = 4
	ngrid = max0(iex,jey)
	if (iex.lt.1) return
	if (jey.lt.1) return
	if (ngrid.gt.50) return
	ierror = 5
	if (nfx.ne.ixp*2**(iex-1)+1) return
	if (nfy.ne.jyq*2**(jey-1)+1) return
	ierror = 6
	if (iguess*(iguess-1).ne.0) return
	ierror = 7
	if (maxcy.lt.1) return
	ierror = 8
	if (method.lt.0 .or. method.gt.3) return
	ierror = 9
c       ****     	
c       ****     compute and test minimum work space
c       ****     
	isx = 0
	if (method.eq.1 .or. method.eq.3) then
	  if (nxa.ne.0) isx = 3
	  if (nxa.eq.0) isx = 5
	end if
	jsy = 0
	if (method.eq.2 .or. method.eq.3) then
	  if (nyc.ne.0) jsy = 3
	  if (nyc.eq.0) jsy = 5
	end if
	kps = 1
	do k=1,ngrid
c         ****     	
c         ****     set subgrid sizes
c         ****     
	  nxk(k) = ixp*2**(max0(k+iex-ngrid,1)-1)+1
	  nyk(k) = jyq*2**(max0(k+jey-ngrid,1)-1)+1
	  nx = nxk(k)
	  ny = nyk(k)
	  kps = kps+(nx+2)*(ny+2)+nx*ny*(10+isx+jsy)
	end do
	iparm(16) = kps+(nfx+2)*(nfy+2)   ! exact minimum work space
	lwork = iparm(16)
	if (lwork .gt. nwork) return
	ierror = 10   ! check solution region
	if (xb.le.xa .or. yd.le.yc) return
	ierror = 11
	if (tolmax .lt. 0.0) return
	ierror = 12   ! multigrid parameters
	if (kcycle.lt.0) return
	if (min0(iprer,ipost).lt.1) return
	if ((intpol-1)*(intpol-3).ne.0) return
	if (max0(kcycle,iprer,ipost).gt.2) then
	  ierror = -5   ! inefficient multigrid cycling
	end if
	if (ierror .gt. 0) ierror = 0   ! no fatal errors
c       ****     
c       ****     set work space pointers and discretize pde at each 
c       ****     grid level
c       ****     
	iw = 1
	do kb=1,ngrid
	  k = ngrid-kb+1
	  nx = nxk(k)
	  ny = nyk(k)
	  kpbgn(k) = iw
	  kcbgn(k) = kpbgn(k)+(nx+2)*(ny+2)
	  ktxbgn(k) = kcbgn(k)+10*nx*ny
	  ktybgn(k) = ktxbgn(k)+isx*nx*ny
	  iw = ktybgn(k)+jsy*nx*ny
	  ic = kcbgn(k)
	  itx = ktxbgn(k)
	  ity = ktybgn(k)
	  klevel = k
	 call dismd2cr(nx,ny,work(ic),work(itx),work(ity),
     +                  work,ierror)
	end do
	return
      end if   ! end of intl=0 initialization call block
      nx = nfx
      ny = nfy
      call mud2cr1(nx,ny,rhs,phi,work)
      iparm(17) = itero
      if (tolmax.gt.0.0) then   ! check for convergence
	fparm(6) = relmax
	if (relmax.gt.tolmax) ierror = -1   ! flag convergenc failure
      end if
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine mud2cr1(nx,ny,rhsf,phif,wk)
      implicit real*8 (a-h,o-z)
      integer nx,ny
      dimension phif(nx,ny),rhsf(nx,ny),wk(*)
      integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,
     +             maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,
     +             kcycle,iprer,ipost,intpol,kps
      integer kpbgn,kcbgn,ktxbgn,ktybgn,nxk,nyk,isx,jsy
      integer k,kb,ip,ic,ir,ipc,irc,icc
      integer ncx,ncy,jj,ij,i,j,iter
      common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,
     +               iguess, maxcy,method,nwork,lwork,itero,ngrid,
     +               klevel,kcur,kcycle,iprer,ipost,intpol,kps
      common/fmud2cr/xa,xb,yc,yd,tolmax,relmax
      common/mud2crc/kpbgn(50),kcbgn(50),ktxbgn(50),ktybgn(50),
     +nxk(50),nyk(50),isx,jsy
      nx = nxk(ngrid)
      ny = nyk(ngrid)
      ip = kpbgn(ngrid)
      ic = kcbgn(ngrid)
      ir = ic+9*nx*ny
c     ****     
c     ****     set phif,rhsf in wk and adjust right hand side
c     ****     
      call swk2(nx,ny,phif,rhsf,wk(ip),wk(ir))
      if (iguess.eq.0) then
c       ****     
c       ****     no initial guess at finest grid level!
c       ****     
	do kb=2,ngrid
	  k = ngrid-kb+1
	  nx = nxk(k+1)
	  ny = nyk(k+1)
	  ip = kpbgn(k+1)
	  ir = kcbgn(k+1)+9*nx*ny
	  ncx = nxk(k)
	  ncy = nyk(k)
	  ipc = kpbgn(k)
	  icc = kcbgn(k)
	  irc = icc+9*ncx*ncy
c         ****     
c         ****     transfer down to all grid levels
c         ****     
	  call trsfc2(nx,ny,wk(ip),wk(ir),ncx,ncy,
     +                wk(ipc),wk(irc))
	end do
c       ****     
c       ****     adjust right hand side at all grid levels in case
c       ****     rhs or specified b.c. in phi or gbdy changed
c       ****     
	do k=1,ngrid
	  nx = nxk(k)
	  ny = nyk(k)
	  ip = kpbgn(k)
	  ic = kcbgn(k)
  	  call adjmd2cr(nx,ny,wk(ip),wk(ic))
	end do
c       ****     
c       ****     execute one full multigrid cycle
c       ****     
	do k=1,ngrid-1
	  kcur = k
	  call kcymd2cr(wk)
	  nx = nxk(k+1)
	  ny = nyk(k+1)
	  ip = kpbgn(k+1)
	  ipc = kpbgn(k)
	  ncx = nxk(k)
	  ncy = nyk(k)
c         ****     
c         ****     lift or prolong approximation from k to k+1
c         ****     
	  call prolon2(ncx,ncy,wk(ipc),nx,ny,wk(ip),nxa,nxb,
     +                 nyc,nyd,intpol)
        end do
      else
c       ****     
c       ****     adjust rhs at finest grid level only
c       ****     
	nx = nxk(ngrid)
	ny = nyk(ngrid)
	ip = kpbgn(ngrid)
	ic = kcbgn(ngrid)
	call adjmd2cr(nx,ny,wk(ip),wk(ic))
      end if
c     ****     
c     ****     execute maxcy more multigrid k cycles from finest level
c     ****     
      kcur = ngrid
      do iter=1,maxcy
	itero = iter
	call kcymd2cr(wk)
	if (tolmax.gt.0.0) then
c         ****     
c         ****     error control
c         ****     
	  relmax = 0.0
	  phmax = 0.0
	  do j=1,nfy
	    jj = j*(nfx+2)
	    do i=1,nfx
	      ij = jj+i+1
	      phmax = DMAX1(phmax,DABS(wk(ij)))
	      relmax = DMAX1(relmax,DABS(wk(ij)-phif(i,j)))
	      phif(i,j) = wk(ij)
	    end do
	  end do
c         ****     
c         ****     set maximum relative difference and check for 
c         ****     convergence
c         ****     
	  if (phmax.gt.0.0) relmax = relmax/phmax
	  if (relmax.le.tolmax) return
	end if
      end do
c     ****     
c     ****     set final interate after maxcy cycles in phif
c     ****     
      do j=1,nfy
	jj = j*(nfx+2)
	do i=1,nfx
	  ij = jj+i+1
	  phif(i,j) = wk(ij)
	end do
      end do
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine prolon1(ncx,p,nx,q,nxa,nxb,intpol)
      implicit real*8 (a-h,o-z)
      integer intpol,nxa,nxb,ncx,nx,i,ic,ist,ifn,ioddst,ioddfn
      dimension p(0:ncx+1),q(0:nx+1)
      ist = 1
      ioddst = 1
      ifn = nx
      ioddfn = nx
      if (nxa.eq.1) then
	ist = 2
	ioddst = 3
      end if
      if (nxb.eq.1) then
	ifn = nx-1
	ioddfn = nx-2
      end if
      if (intpol.eq.1 .or. ncx.lt.4) then
c       ****     
c       ****     linear interpolation in x
c       ****     
	if (ncx .lt. nx) then
c         ****     
c         ****     every other point of nx grid is ncx grid
c         ****     
	  do i=ioddst,ioddfn,2
	    ic = (i+1)/2
	    q(i) = p(ic)
	  end do
	  do i=2,ifn,2
	    q(i) = 0.5*(q(i-1)+q(i+1))
	  end do
	else
c         ****     
c         ****     nx grid equals ncx grid
c         ****     
	  do i=ist,ifn
	    q(i) = p(i)
	  end do
	end if
c       ****     
c       ****     set virtual end points if periodic
c       ****     
	if (nxa.eq.0) then
	  q(0) = q(nx-1)
	  q(nx+1) = q(2)
	end if
	return
      else
c       ****     
c       ****     cubic interpolation in x
c       ****     
	if (ncx .lt. nx) then
	  do i=ioddst,ioddfn,2
	    ic = (i+1)/2
	    q(i) = p(ic)
	  end do
c         ****     
c         ****       set deep interior with symmetric formula
c         ****     
	  do i=4,nx-3,2
	    q(i)=(-q(i-3)+9.*(q(i-1)+q(i+1))-q(i+3))*.0625
	  end do
c         ****     
c         ****     interpolate from q at i=2 and i=nx-1
c         ****     
	  if (nxa.ne.0) then
c           ****     
c           ****     asymmetric formula near nonperiodic bndys
c           ****     
  	    q(2)=(5.*q(1)+15.*q(3)-5.*q(5)+q(7))*.0625
	    q(nx-1)=(5.*q(nx)+15.*q(nx-2)-5.*q(nx-4)+q(nx-6))*.0625
	  else
c           ****     
c           ****     periodicity in x alows symmetric formula near bndys
c           ****     
	    q(2) = (-q(nx-2)+9.*(q(1)+q(3))-q(5))*.0625
	    q(nx-1) = (-q(nx-4)+9.*(q(nx-2)+q(nx))-q(3))*.0625
	    q(nx+1) = q(2)
	    q(0) = q(nx-1)
	  end if
	  return
	else
c         ****     
c         ****     ncx grid equals nx grid
c         ****     
	  do i=ist,ifn
	    q(i) = p(i)
	  end do
	  if (nxa.eq.0) then
	    q(0) = q(nx-1)
	    q(nx+1) = q(2)
	  end if
	  return
	end if
      end if
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine prolon2(ncx,ncy,p,nx,ny,q,nxa,nxb,nyc,nyd,intpol)
      implicit real*8 (a-h,o-z)
      integer ncx,ncy,nx,ny,intpol,nxa,nxb,nyc,nyd
      dimension p(0:ncx+1,0:ncy+1),q(0:nx+1,0:ny+1)
      integer i,j,jc,ist,ifn,jst,jfn,joddst,joddfn
      ist = 1
      ifn = nx
      jst = 1
      jfn = ny
      joddst = 1
      joddfn = ny
      if (nxa.eq.1) then
	ist = 2
      end if
      if (nxb.eq.1) then
	ifn = nx-1
      end if
      if (nyc.eq.1) then
	jst = 2
	joddst = 3
      end if
      if (nyd.eq.1) then
	jfn = ny-1
	joddfn = ny-2
      end if
      if (intpol.eq.1 .or. ncy.lt.4) then
c       ****     
c       ****     linearly interpolate in y
c       ****     
	if (ncy .lt. ny) then
c         ****     
c         ****     ncy grid is an every other point subset of ny grid
c         ****     set odd j lines interpolating in x and then set even
c         ****     j lines by averaging odd j lines
c         ****     
	  do j=joddst,joddfn,2
	    jc = j/2+1
	    call prolon1(ncx,p(0,jc),nx,q(0,j),nxa,nxb,intpol)
	  end do
	  do j=2,jfn,2
	    do i=ist,ifn
	      q(i,j) = 0.5*(q(i,j-1)+q(i,j+1))
	    end do
	  end do
c         ****     
c         ****     set periodic virtual boundaries if necessary
c         ****     
	  if (nyc.eq.0) then
	    do i=ist,ifn
	      q(i,0) = q(i,ny-1)
	      q(i,ny+1) = q(i,2)
	    end do
	  end if
	  return
	else
c         ****     
c         ****     ncy grid is equals ny grid so interpolate in x only
c         ****     
	  do j=jst,jfn
	    jc = j
	    call prolon1(ncx,p(0,jc),nx,q(0,j),nxa,nxb,intpol)
	  end do
c         ****     
c         ****     set periodic virtual boundaries if necessary
c         ****     
	  if (nyc.eq.0) then
	    do i=ist,ifn
	      q(i,0) = q(i,ny-1)
	      q(i,ny+1) = q(i,2)
	    end do
	  end if
	  return
	end if
      else
c       ****     
c       ****     cubically interpolate in y
c       ****     
	if (ncy .lt. ny) then
c         ****     
c         ****     set every other point of ny grid by interpolating 
c         ****     in x
c         ****     
	  do j=joddst,joddfn,2
	    jc = j/2+1
	    call prolon1(ncx,p(0,jc),nx,q(0,j),nxa,nxb,intpol)
	  end do
c         ****     
c         ****     set deep interior of ny grid using values just
c         ****     generated and symmetric cubic interpolation in y
c         ****     
	  do j=4,ny-3,2
	    do i=ist,ifn
	    q(i,j)=(-q(i,j-3)+9.*(q(i,j-1)+q(i,j+1))-q(i,j+3))*.0625
	    end do
	  end do
c         ****     
c         ****     interpolate from q at j=2 and j=ny-1
c         ****     
	  if (nyc.ne.0) then
c           ****     
c           ****     asymmetric formula near nonperiodic y boundaries
c           ****     
	    do i=ist,ifn
	      q(i,2)=(5.*q(i,1)+15.*q(i,3)-5.*q(i,5)+q(i,7))*.0625
	      q(i,ny-1)=(5.*q(i,ny)+15.*q(i,ny-2)-5.*q(i,ny-4)+
     +                  q(i,ny-6))*.0625
	    end do
	  else
c           ****     
c           ****     periodicity in y alows symmetric formula near bndys
c           ****     
	    do i=ist,ifn
	      q(i,2) = (-q(i,ny-2)+9.*(q(i,1)+q(i,3))-q(i,5))*.0625
	      q(i,ny-1)=(-q(i,ny-4)+9.*(q(i,ny-2)+q(i,ny))-q(i,3))*.0625
	      q(i,ny+1) = q(i,2)
	      q(i,0) = q(i,ny-1)
	    end do
	  end if
	  return
	else
c        ****     
c        ****     ncy grid is equals ny grid so interpolate in x only
c        ****     
	  do j=jst,jfn
	    jc = j
	    call prolon1(ncx,p(0,jc),nx,q(0,j),nxa,nxb,intpol)
	  end do
c         ****     
c         ****     set periodic virtual boundaries if necessary
c         ****     
	  if (nyc.eq.0) then
	    do i=ist,ifn
	      q(i,0) = q(i,ny-1)
	      q(i,ny+1) = q(i,2)
	    end do
	  end if
	  return
	end if
      end if
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine relmd2cr(nx,ny,phi,cof,tx,ty,sum)
c     ****     
c     ****     relaxation for mud2
c     ****     
      implicit real*8 (a-h,o-z)
      integer nx,ny
      dimension phi(*),cof(*),tx(*),ty(*),sum(*)
      integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,
     +             maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,
     +             kcycle,iprer,ipost,intpol,kps
      common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,
     +               iguess, maxcy,method,nwork,lwork,itero,ngrid,
     +               klevel,kcur,kcycle,iprer,ipost,intpol,kps
      if (method.eq.0) then                ! point relaxation
	call relmd2crp(nx,ny,phi,cof)
      else if (method.eq.1) then           ! line x relaxation
	call slxmd2cr(nx,ny,phi,cof,tx,sum)
      else if (method.eq.2) then           ! line y relaxation
	call slymd2cr(nx,ny,phi,cof,ty,sum)
      else if (method.eq.3) then           ! line x&y relaxation
	call slxmd2cr(nx,ny,phi,cof,tx,sum)
	call slymd2cr(nx,ny,phi,cof,ty,sum)
      end if
      return
      end      
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine relmd2crp(nx,ny,phi,cof)
c     ****     
c     ****     gauss-seidel four color point relaxation
c     ****     
      implicit real*8 (a-h,o-z)
      integer nx,ny,i,j,lcolor,i1,i2,i3,i4,it
      integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,
     +             maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,
     +             kcycle,iprer,ipost,intpol,kps
      common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,
     +               iguess, maxcy,method,nwork,lwork,itero,ngrid,
     +               klevel,kcur,kcycle,iprer,ipost,intpol,kps
      dimension phi(0:nx+1,0:ny+1),cof(nx,ny,10)
      i1 = 1
      i2 = 4
      i3 = 3
      i4 = 2
c     ****     
c     ****     sweep four colored grid points
c     ****     
      do lcolor=1,4
	do j=1,ny,4
	  do i=i1,nx,4
	      phi(i,j) = (cof(i,j,10) - (
     +                    cof(i,j,1)*phi(i+1,j)   +
     +                    cof(i,j,2)*phi(i+1,j+1) +
     +                    cof(i,j,3)*phi(i,j+1)   +
     +                    cof(i,j,4)*phi(i-1,j+1) +
     +                    cof(i,j,5)*phi(i-1,j)   +
     +                    cof(i,j,6)*phi(i-1,j-1) +
     +                    cof(i,j,7)*phi(i,j-1)   +
     +                    cof(i,j,8)*phi(i+1,j-1)))/cof(i,j,9)
	  end do
	end do
	do j=2,ny,4
	  do i=i2,nx,4
	      phi(i,j) = (cof(i,j,10) - (
     +                    cof(i,j,1)*phi(i+1,j)   +
     +                    cof(i,j,2)*phi(i+1,j+1) +
     +                    cof(i,j,3)*phi(i,j+1)   +
     +                    cof(i,j,4)*phi(i-1,j+1) +
     +                    cof(i,j,5)*phi(i-1,j)   +
     +                    cof(i,j,6)*phi(i-1,j-1) +
     +                    cof(i,j,7)*phi(i,j-1)   +
     +                    cof(i,j,8)*phi(i+1,j-1)))/cof(i,j,9)
	  end do
	end do
	do j=3,ny,4
	  do i=i3,nx,4
	      phi(i,j) = (cof(i,j,10) - (
     +                    cof(i,j,1)*phi(i+1,j)   +
     +                    cof(i,j,2)*phi(i+1,j+1) +
     +                    cof(i,j,3)*phi(i,j+1)   +
     +                    cof(i,j,4)*phi(i-1,j+1) +
     +                    cof(i,j,5)*phi(i-1,j)   +
     +                    cof(i,j,6)*phi(i-1,j-1) +
     +                    cof(i,j,7)*phi(i,j-1)   +
     +                    cof(i,j,8)*phi(i+1,j-1)))/cof(i,j,9)
	  end do
	end do
	do j=4,ny,4
	  do i=i4,nx,4
	      phi(i,j) = (cof(i,j,10) - (
     +                    cof(i,j,1)*phi(i+1,j)   +
     +                    cof(i,j,2)*phi(i+1,j+1) +
     +                    cof(i,j,3)*phi(i,j+1)   +
     +                    cof(i,j,4)*phi(i-1,j+1) +
     +                    cof(i,j,5)*phi(i-1,j)   +
     +                    cof(i,j,6)*phi(i-1,j-1) +
     +                    cof(i,j,7)*phi(i,j-1)   +
     +                    cof(i,j,8)*phi(i+1,j-1)))/cof(i,j,9)
	  end do
	end do
c       ****     
c       ****     set periodic virtual boundaries as necessary
c       ****     
	if (nxa.eq.0) then
	  do j=1,ny
	    phi(0,j) = phi(nx-1,j)
	    phi(nx+1,j) = phi(2,j)
	  end do
	end if
	if (nyc.eq.0) then
	  do i=1,nx
	    phi(i,0) = phi(i,ny-1)
	    phi(i,ny+1) = phi(i,2)
	  end do
	end if
c       ****     
c       ****     permute (i1,i2,i3,i4) for next color
c       ****     
	it = i4
	i4 = i3
	i3 = i2
	i2 = i1
	i1 = it
      end do
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine res2(nx,ny,resf,ncx,ncy,rhsc,nxa,nxb,nyc,nyd)
      implicit real*8 (a-h,o-z)
      integer nx,ny,ncx,ncy,nxa,nxb,nyc,nyd
      integer i,j,ic,jc,im1,ip1,jm1,jp1,ix,jy
c     ****     
c     ****     restrict fine grid residual in resf to coarse grid in 
c     ****     rhsc using full weighting for all real 2d codes
c     ****     
      dimension resf(nx,ny),rhsc(ncx,ncy)
c     ****     
c     ****     set x,y coarsening integer subscript scales
c     ****     
      ix = 1
      if (ncx.eq.nx) ix = 0
      jy = 1
      if (ncy.eq.ny) jy = 0
c     ****     
c     ****     restrict on interior
c     ****     
      if (ncy.lt.ny .and. ncx.lt.nx) then
c       ****     
c       ****     coarsening in both directions
c       ****     
	do jc=2,ncy-1
	  j = jc+jc-1
	  do ic=2,ncx-1
	    i = ic+ic-1
	    rhsc(ic,jc) = (resf(i-1,j-1)+resf(i+1,j-1)+resf(i-1,j+1)+
     +                     resf(i+1,j+1)+2.*(resf(i-1,j)+resf(i+1,j)+
     +                     resf(i,j-1)+resf(i,j+1))+4.*resf(i,j))*.0625
	  end do
	end do
      else if (ncy.eq.ny) then
c       ****     
c       ****     no coarsening in y but coarsening in x
c       ****     
	do jc=2,ncy-1
	  j = jc
	  do ic=2,ncx-1
	    i = ic+ic-1
	    rhsc(ic,jc) = (resf(i-1,j-1)+resf(i+1,j-1)+resf(i-1,j+1)+
     +                     resf(i+1,j+1)+2.*(resf(i-1,j)+resf(i+1,j)+
     +                     resf(i,j-1)+resf(i,j+1))+4.*resf(i,j))*.0625
	  end do
	end do
      else
c       ****     
c       ****     no coarsening in x but coarsening in y
c       ****     
	do jc=2,ncy-1
	  j = jc+jc-1
	  do ic=2,ncx-1
	    i = ic
	    rhsc(ic,jc) = (resf(i-1,j-1)+resf(i+1,j-1)+resf(i-1,j+1)+
     +                     resf(i+1,j+1)+2.*(resf(i-1,j)+resf(i+1,j)+
     +                     resf(i,j-1)+resf(i,j+1))+4.*resf(i,j))*.0625
	  end do
	end do
      end if
c     ****     
c     ****     set residual on boundaries
c     ****     
      do jc=1,ncy,ncy-1
c       ****     
c       ****     y=yc,yd boundaries
c       ****     
	j = jc+jy*(jc-1)
	jm1 = max0(j-1,2)
	jp1 = min0(j+1,ny-1)
	if (j.eq.1 .and. nyc.eq.0) jm1 = ny-1
	if (j.eq.ny .and. nyc.eq.0) jp1 = 2
c       ****     
c       ****     y=yc,yd and x=xa,xb cornors
c       ****     
	do ic=1,ncx,ncx-1
	  i = ic+ix*(ic-1)
	  im1 = max0(i-1,2)
	  ip1 = min0(i+1,nx-1)
	  if (i.eq.1 .and. nxa.eq.0) im1 = nx-1
	  if (i.eq.nx .and. nxa.eq.0) ip1 = 2
	  rhsc(ic,jc) = (resf(im1,jm1)+resf(ip1,jm1)+resf(im1,jp1)+
     +               resf(ip1,jp1)+2.*(resf(im1,j)+resf(ip1,j)+
     +               resf(i,jm1)+resf(i,jp1))+4.*resf(i,j))*.0625
	end do
c       ****     
c       ****     set y=yc,yd interior edges
c       ****     
	do ic=2,ncx-1
	  i = ic+ix*(ic-1)
	  rhsc(ic,jc) = (resf(i-1,jm1)+resf(i+1,jm1)+resf(i-1,jp1)+
     +               resf(i+1,jp1)+2.*(resf(i-1,j)+resf(i+1,j)+
     +               resf(i,jm1)+resf(i,jp1))+4.*resf(i,j))*.0625
        enddo
      end do
c     ****     
c     ****     set x=xa,xb interior edges
c     ****     
      do ic=1,ncx,ncx-1
	i = ic+ix*(ic-1)
	im1 = max0(i-1,2)
	ip1 = min0(i+1,nx-1)
	if (i.eq.1 .and. nxa.eq.0) im1 = nx-1
	if (i.eq.nx .and. nxa.eq.0) ip1 = 2
	do jc=2,ncy-1
	  j = jc+jy*(jc-1)
	  rhsc(ic,jc) = (resf(im1,j-1)+resf(ip1,j-1)+resf(im1,j+1)+
     +               resf(ip1,j+1)+2.*(resf(im1,j)+resf(ip1,j)+
     +               resf(i,j-1)+resf(i,j+1))+4.*resf(i,j))*.0625
	end do
      end do
c     ****     
c     ****     set coarse grid residual zero on specified boundaries
c     ****     
      if (nxa.eq.1) then
	do jc=1,ncy
	  rhsc(1,jc) = 0.0
	end do
      end if
      if (nxb.eq.1) then
	do jc=1,ncy
	  rhsc(ncx,jc) = 0.0
	end do
      end if
      if (nyc.eq.1) then
	do ic=1,ncx
	  rhsc(ic,1) = 0.0
	end do
      end if
      if (nyd.eq.1) then
	do ic=1,ncx
	  rhsc(ic,ncy) = 0.0
	end do
      end if
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc      
      subroutine resmd2cr(nx,ny,phi,ncx,ncy,phic,rhsc,cof,resf)
c     ****     
c     ****     restrict residual from fine to coarse mesh using fully 
c     ****     weighted residual restriction
c     ****     
      implicit real*8 (a-h,o-z)
      integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,
     +             maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,
     +             kcycle,iprer,ipost,intpol,kps
      integer nx,ny,ncx,ncy,i,j,ic,jc
      common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,
     +               iguess, maxcy,method,nwork,lwork,itero,ngrid,
     +               klevel,kcur,kcycle,iprer,ipost,intpol,kps
      dimension rhsc(ncx,ncy),resf(nx,ny)
      dimension phi(0:nx+1,0:ny+1),phic(0:ncx+1,0:ncy+1)
      dimension cof(nx,ny,10)
c     ****     
c     ****     set phic zero
c     ****     
      do jc=0,ncy+1
	do ic=0,ncx+1
	  phic(ic,jc) = 0.0
	end do
      end do
c     ****     
c     ****     compute residual on fine mesh in resf
c     ****     
      do j=1,ny
	do i=1,nx
	  resf(i,j) = cof(i,j,10)-(
     +                cof(i,j,1)*phi(i+1,j)+
     +                cof(i,j,2)*phi(i+1,j+1)+
     +                cof(i,j,3)*phi(i,j+1)+
     +                cof(i,j,4)*phi(i-1,j+1)+
     +                cof(i,j,5)*phi(i-1,j)+
     +                cof(i,j,6)*phi(i-1,j-1)+
     +                cof(i,j,7)*phi(i,j-1)+
     +                cof(i,j,8)*phi(i+1,j-1)+
     +                cof(i,j,9)*phi(i,j))
	end do
      end do
c     ****     
c     ****     restrict resf to coarse mesh in rhsc
c     ****     
      call res2(nx,ny,resf,ncx,ncy,rhsc,nxa,nxb,nyc,nyd)
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      SUBROUTINE SEEE(CEE,NX,NY,CF)
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION CEE(NX,NY,1),CF(NX,NY,1)
      DO 1 N = 1,9
        DO J = 1,NY
          DO I = 1,NX
            CF(I,J,N) = CEE(I,J,N)
          enddo
       enddo
    1 CONTINUE
      RETURN
      END
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine slxmd2cr(nx,ny,phi,cof,tx,sum)
c     ****     
c     ****     line relaxation in the x direction (periodic or 
c     ****     nonperiodic)
c     ****     
      implicit real*8 (a-h,o-z)
      integer nx,ny,i,ib,j
      integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,
     +             maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,
     +             kcycle,iprer,ipost,intpol,kps
      common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,
     +               iguess, maxcy,method,nwork,lwork,itero,ngrid,
     +               klevel,kcur,kcycle,iprer,ipost,intpol,kps
      dimension phi(0:nx+1,0:ny+1),cof(nx,ny,10),tx(nx,ny,*),sum(ny)
c     ****     
c     ****     replace line x with point gauss-seidel if
c     ****     x direction is periodic and nx = 3 (coarsest)
c     ****     
      if (nxa .eq. 0 .and. nx .eq. 3) then
	call relmd2crp(nx,ny,phi,cof)
	return
      end if
c     ****     
c     ****     set periodic y virtual boundary if necessary
c     ****     
      if (nyc.eq.0) then
	do i=1,nx
	  phi(i,0) = phi(i,ny-1)
	  phi(i,ny+1) = phi(i,2)
	end do
      end if
      if (nxa.ne.0) then
c       ****     
c       ****     x direction not periodic, sweep odd j lines
c       ****     
	do j=1,ny,2
	  do i=1,nx
	    phi(i,j) = cof(i,j,10)-(cof(i,j,2)*phi(i+1,j+1)+
     +                              cof(i,j,3)*phi(i,j+1)+
     +                              cof(i,j,4)*phi(i-1,j+1)+
     +                              cof(i,j,6)*phi(i-1,j-1)+
     +                              cof(i,j,7)*phi(i,j-1)+
     +                              cof(i,j,8)*phi(i+1,j-1))
	  end do
c         ****     
c         ****     forward sweep
c         ****     
	  do i=2,nx
	    phi(i,j) = phi(i,j)-tx(i-1,j,1)*phi(i-1,j)
	  end do
c         ****     
c         ****     backward sweep
c         ****     
	  phi(nx,j) = phi(nx,j)/tx(nx,j,2)
	  do ib=2,nx
	    i = nx-ib+1
	    phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j))/tx(i,j,2)
	  end do
	end do
c       ****     
c       ****     sweep even j lines forward and back
c       ****     
	do j=2,ny,2
	  do i=1,nx
	    phi(i,j) = cof(i,j,10)-(cof(i,j,2)*phi(i+1,j+1)+
     +                              cof(i,j,3)*phi(i,j+1)+
     +                              cof(i,j,4)*phi(i-1,j+1)+
     +                              cof(i,j,6)*phi(i-1,j-1)+
     +                              cof(i,j,7)*phi(i,j-1)+
     +                              cof(i,j,8)*phi(i+1,j-1))
	  end do
	  do i=2,nx
	    phi(i,j) = phi(i,j)-tx(i-1,j,1)*phi(i-1,j)
	  end do
	  phi(nx,j) = phi(nx,j)/tx(nx,j,2)
	  do ib=2,nx
	    i = nx-ib+1
	    phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j))/tx(i,j,2)
	  end do
	end do
      else
c       ****     
c       ****     x direction periodic
c       ****     
	do j=1,ny
	  sum(j) = 0.0
	  phi(0,j) = phi(nx-1,j)
	  phi(nx+1,j) = phi(2,j)
	end do
c       ****     
c       ****     sweep odd lines forward and back
c       ****     
	do j=1,ny,2
	  do i=1,nx-1
	    phi(i,j) = cof(i,j,10)-(cof(i,j,2)*phi(i+1,j+1)+
     +                              cof(i,j,3)*phi(i,j+1)+
     +                              cof(i,j,4)*phi(i-1,j+1)+
     +                              cof(i,j,6)*phi(i-1,j-1)+
     +                              cof(i,j,7)*phi(i,j-1)+
     +                              cof(i,j,8)*phi(i+1,j-1))
	  end do
c         ****     
c         ****     forward sweep
c         ****     
	  do i=2,nx-2
	    phi(i,j) = phi(i,j)-tx(i,j,1)*phi(i-1,j)
	  end do
	  do i=1,nx-2
	    sum(j) = sum(j)+tx(i,j,5)*phi(i,j)
	  end do
	  phi(nx-1,j) = phi(nx-1,j)-sum(j)
c         ****     
c         ****     backward sweep
c         ****     
	  phi(nx-1,j) = phi(nx-1,j)/tx(nx-1,j,2)
	  phi(nx-2,j) = (phi(nx-2,j)-tx(nx-2,j,4)*phi(nx-1,j))/
     +                   tx(nx-2,j,2)
	  do ib=4,nx
	    i = nx-ib+1
	    phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j)-tx(i,j,4)*
     +                 phi(nx-1,j))/tx(i,j,2)
	  end do
	end do
c       ****     
c       ****     set periodic and virtual points for j odd
c       ****     
	do j=1,ny,2
	  phi(nx,j) = phi(1,j)
	  phi(0,j) = phi(nx-1,j)
	  phi(nx+1,j) = phi(2,j)
	end do
c       ****     
c       ****     sweep even j lines
c       ****     
	do j=2,ny,2
	  do i=1,nx-1
	    phi(i,j) = cof(i,j,10)-(cof(i,j,2)*phi(i+1,j+1)+
     +                              cof(i,j,3)*phi(i,j+1)+
     +                              cof(i,j,4)*phi(i-1,j+1)+
     +                              cof(i,j,6)*phi(i-1,j-1)+
     +                              cof(i,j,7)*phi(i,j-1)+
     +                              cof(i,j,8)*phi(i+1,j-1))
	  end do
c         ****     
c         ****     forward sweep
c         ****     
	  do i=2,nx-2
	    phi(i,j) = phi(i,j)-tx(i,j,1)*phi(i-1,j)
	  end do
	  do i=1,nx-2
	    sum(j) = sum(j)+tx(i,j,5)*phi(i,j)
	  end do
	  phi(nx-1,j) = phi(nx-1,j)-sum(j)
c         ****     
c         ****     backward sweep
c         ****     
	  phi(nx-1,j) = phi(nx-1,j)/tx(nx-1,j,2)
	  phi(nx-2,j) = (phi(nx-2,j)-tx(nx-2,j,4)*phi(nx-1,j))/
     +                   tx(nx-2,j,2)
	  do ib=4,nx
	    i = nx-ib+1
	    phi(i,j) = (phi(i,j)-tx(i,j,3)*phi(i+1,j)-tx(i,j,4)*
     +                 phi(nx-1,j))/tx(i,j,2)
	  end do
	end do
c       ****     
c       ****     set periodic and virtual points for j even
c       ****     
	do j=2,ny,2
	  phi(nx,j) = phi(1,j)
	  phi(0,j) = phi(nx-1,j)
	  phi(nx+1,j) = phi(2,j)
	end do
       do j=1,ny
        do i=1,nx
        enddo
      enddo
      end if
c     ****     
c     ****     set periodic y virtual boundaries if necessary
c     ****     
      if (nyc.eq.0) then
	do i=1,nx
	  phi(i,0) = phi(i,ny-1)
	  phi(i,ny+1) = phi(i,2)
	end do
      end if
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine slymd2cr(nx,ny,phi,cof,ty,sum)
      implicit real*8 (a-h,o-z)
      integer nx,ny,i,j,jb
      integer intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,iguess,
     +             maxcy,method,nwork,lwork,itero,ngrid,klevel,kcur,
     +             kcycle,iprer,ipost,intpol,kps
      common/imud2cr/intl,nxa,nxb,nyc,nyd,ixp,jyq,iex,jey,nfx,nfy,
     +               iguess, maxcy,method,nwork,lwork,itero,ngrid,
     +               klevel,kcur,kcycle,iprer,ipost,intpol,kps
      dimension phi(0:nx+1,0:ny+1),cof(nx,ny,10),ty(ny,nx,*),sum(nx)
c     ****     
c     ****     replace line y with point gauss-seidel if
c     ****     y direction is periodic and ny = 3
c     ****     
      if (nyc .eq. 0 .and. ny .eq. 3) then
	call relmd2crp(nx,ny,phi,cof)
	return
      end if
c     ****     
c     ****     set periodic and virtual x boundaries if necessary
c     ****     
      if (nxa.eq.0) then
	do j=1,ny
	  phi(0,j) = phi(nx-1,j)
	  phi(nx,j) = phi(1,j)
	  phi(nx+1,j) = phi(2,j)
	end do
      end if
      if (nyc.ne.0) then
c       ****     
c       ****     y direction not periodic
c       ****     
	do i=1,nx,2
	  do j=1,ny
	    phi(i,j) = cof(i,j,10)-(cof(i,j,1)*phi(i+1,j)+
     +                              cof(i,j,2)*phi(i+1,j+1)+
     +                              cof(i,j,4)*phi(i-1,j+1)+
     +                              cof(i,j,5)*phi(i-1,j)+
     +                              cof(i,j,6)*phi(i-1,j-1)+
     +                              cof(i,j,8)*phi(i+1,j-1))
	  end do
c         ****     
c         ****     forward sweep thru odd x lines
c         ****     
	  do j=2,ny
	    phi(i,j) = phi(i,j)-ty(j-1,i,1)*phi(i,j-1)
	  end do
c         ****     
c         ****     backward sweep
c         ****     
	  phi(i,ny) = phi(i,ny)/ty(ny,i,2)
	  do jb=2,ny
	    j = ny-jb+1
	    phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1))/ty(j,i,2)
	  end do
	end do
c       ****     
c       ****     forward sweep even x lines
c       ****     
	do i=2,nx,2
	  do j=1,ny
	    phi(i,j) = cof(i,j,10)-(cof(i,j,1)*phi(i+1,j)+
     +                              cof(i,j,2)*phi(i+1,j+1)+
     +                              cof(i,j,4)*phi(i-1,j+1)+
     +                              cof(i,j,5)*phi(i-1,j)+
     +                              cof(i,j,6)*phi(i-1,j-1)+
     +                              cof(i,j,8)*phi(i+1,j-1))
	  end do
	  do j=2,ny
	    phi(i,j) = phi(i,j)-ty(j-1,i,1)*phi(i,j-1)
	  end do
c         ****     
c         ****     backward sweep
c         ****     
	  phi(i,ny) = phi(i,ny)/ty(ny,i,2)
	  do jb=2,ny
	    j = ny-jb+1
	    phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1))/ty(j,i,2)
	  end do
	end do
      else
c       ****     
c       ****     y direction periodic
c       ****     
	do i=1,nx
	  sum(i) = 0.0
	  phi(i,0) = phi(i,ny-1)
	  phi(i,ny) = phi(i,1)
	  phi(i,ny+1) = phi(i,2)
	end do
c       ****     
c       ****     forward sweep odd x lines
c       ****     
	do i=1,nx,2
	  do j=1,ny-1
	    phi(i,j) = cof(i,j,10)-(cof(i,j,1)*phi(i+1,j)+
     +                              cof(i,j,2)*phi(i+1,j+1)+
     +                              cof(i,j,4)*phi(i-1,j+1)+
     +                              cof(i,j,5)*phi(i-1,j)+
     +                              cof(i,j,6)*phi(i-1,j-1)+
     +                              cof(i,j,8)*phi(i+1,j-1))
	  end do
	  do j=2,ny-2
	    phi(i,j) = phi(i,j)-ty(j,i,1)*phi(i,j-1)
	  end do
	  do j=1,ny-2
	    sum(i) = sum(i)+ty(j,i,5)*phi(i,j)
	  end do
	  phi(i,ny-1) = phi(i,ny-1)-sum(i)
c         ****     
c         ****     backward sweep
c         ****     
	  phi(i,ny-1) = phi(i,ny-1)/ty(ny-1,i,2)
	  phi(i,ny-2) = (phi(i,ny-2)-ty(ny-2,i,4)*phi(i,ny-1))/
     +                   ty(ny-2,i,2)
	  do jb=4,ny
	    j = ny-jb+1
	    phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1)-ty(j,i,4)*
     +                  phi(i,ny-1))/ty(j,i,2)
	  end do
	end do
c       ****     
c       ****     set odd periodic and virtual y boundaries
c       ****     
	do i=1,nx,2
	  phi(i,0) = phi(i,ny-1)
	  phi(i,ny) = phi(i,1)
	  phi(i,ny+1) = phi(i,2)
	end do
c       ****     
c       ****     forward sweep even x lines
c       ****     
	do i=2,nx,2
	  do j=1,ny-1
	    phi(i,j) = cof(i,j,10)-(cof(i,j,1)*phi(i+1,j)+
     +                              cof(i,j,2)*phi(i+1,j+1)+
     +                              cof(i,j,4)*phi(i-1,j+1)+
     +                              cof(i,j,5)*phi(i-1,j)+
     +                              cof(i,j,6)*phi(i-1,j-1)+
     +                              cof(i,j,8)*phi(i+1,j-1))
	  end do
	  do j=2,ny-2
	    phi(i,j) = phi(i,j)-ty(j,i,1)*phi(i,j-1)
	  end do
	  do j=1,ny-2
	    sum(i) = sum(i)+ty(j,i,5)*phi(i,j)
	  end do
	  phi(i,ny-1) = phi(i,ny-1)-sum(i)
c         ****     
c         ****     backward sweep
c         ****     
	  phi(i,ny-1) = phi(i,ny-1)/ty(ny-1,i,2)
	  phi(i,ny-2) = (phi(i,ny-2)-ty(ny-2,i,4)*phi(i,ny-1))/
     +                   ty(ny-2,i,2)
	  do jb=4,ny
	    j = ny-jb+1
	    phi(i,j) = (phi(i,j)-ty(j,i,3)*phi(i,j+1)-ty(j,i,4)*
     +                  phi(i,ny-1))/ty(j,i,2)
	  end do
	end do
c       ****     
c       ****     set even periodic and virtual y boundaries
c       ****     
	do i=2,nx,2
	  phi(i,0) = phi(i,ny-1)
	  phi(i,ny) = phi(i,1)
	  phi(i,ny+1) = phi(i,2)
	end do
      end if
c     ****     
c     ****     set periodic and virtual x boundaries if necessary
c     ****     
      if (nxa.eq.0) then
	do j=1,ny
	  phi(0,j) = phi(nx-1,j)
	  phi(nx+1,j) = phi(2,j)
	end do
      end if
      return
      end

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine swk2(nfx,nfy,phif,rhsf,phi,rhs)
c     ****     
c     ****     set phif,rhsf input in arrays which include
c     ****     virtual boundaries for phi (for all 2-d real codes)
c     ****     
      implicit real*8 (a-h,o-z)
      integer nfx,nfy,i,j
      dimension phif(nfx,nfy),rhsf(nfx,nfy)
      dimension phi(0:nfx+1,0:nfy+1),rhs(nfx,nfy)
      do j=1,nfy
	do i=1,nfx
	  phi(i,j) = phif(i,j)
	  rhs(i,j) = rhsf(i,j)
	end do
      end do
c     ****     
c     ****     set virtual boundaries in phi to zero
c     ****     
      do j=0,nfy+1
	phi(0,j) = 0.0
	phi(nfx+1,j) = 0.0
      end do
      do i=0,nfx+1
	phi(i,0) = 0.0
	phi(i,nfy+1) = 0.0
      end do
      return
      end
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine trsfc2(nx,ny,phi,rhs,ncx,ncy,phic,rhsc)
c     ****     
c     ****     transfer fine grid to coarse grid
c     ****     
      implicit real*8 (a-h,o-z)
      integer nx,ny,ncx,ncy,i,j,ic,jc
      dimension phi(0:nx+1,0:ny+1),rhs(nx,ny)
      dimension phic(0:ncx+1,0:ncy+1),rhsc(ncx,ncy)
c     ****     
c     ****     set virtual boundaries in phic to zero
c     ****     
      do jc=0,ncy+1
	phic(0,jc) = 0.0
	phic(ncx+1,jc) = 0.0
      end do
      do ic=0,ncx+1
	phic(ic,0) = 0.0
	phic(ic,ncy+1) = 0.0
      end do
      if (ncx.lt.nx .and. ncy.lt.ny) then
c       ****     
c       ****     coarsening in both x and y
c       ****     
	do jc=1,ncy
	  j = jc+jc-1
	  do ic=1,ncx
	    i = ic+ic-1
	    phic(ic,jc) = phi(i,j)
	    rhsc(ic,jc) = rhs(i,j)
	  end do
	end do
      else if (ncx.lt.nx .and. ncy.eq.ny) then
c       ****     
c       ****     coarsening in x only
c       ****     
	do jc=1,ncy
	  j = jc
	  do ic=1,ncx
	    i = ic+ic-1
	    phic(ic,jc) = phi(i,j)
	    rhsc(ic,jc) = rhs(i,j)
	  end do
	end do
      else
c       ****     
c       ****     coarsening in y only
c       ****     
	do jc=1,ncy
	  j = jc+jc-1
	  do ic=1,ncx
	    i = ic
	    phic(ic,jc) = phi(i,j)
	    rhsc(ic,jc) = rhs(i,j)
	  end do
	end do
      end if
      return
      end
                 
