subroutine conjug(rem,x,x0,p,y,ab,nxy,nx,ny)
c*********************************************************************
c
c     solve the linear system A x = rem
c     by the unconditioned conjugate gradient method
c
c     p,y     :      work space
c     nxy     ->     dimension of rem and x
c     ab      ->     lambda in CN scheme matrix
c
c     matrix multiplication by routine atimex
c     accuracy controlled by:
c            || e || 2 <= eps.
c
c*********************************************************************
      integer nxy,nx,ny,i,its,itmax
      double precision rem(nxy),x(nxy),x0(nxy),p(nxy),y(nxy),
     +                 ab,eps,rnosq,rnnsq,papt,ak,bk
      parameter (eps=1.D-10)
c---  initialize
      do 5 i=1,nxy
         x(i) = x0(i)
 5    continue
      call atimex(x,y,ab,nxy,nx,ny)
      itmax = nxy
      rnosq = 0.D+00
      do 10 i=1,nxy
         rem(i) = rem(i)-y(i)
         p(i) = rem(i)
         rnosq = rnosq + rem(i)*rem(i)
   10 continue
c---  starting iterations
      do 50 its=1,itmax
         call atimex(p,y,ab,nxy,nx,ny)
         papt = 0.D+00
         do 20 i=1,nxy
            papt = papt + p(i)*y(i)
   20    continue
         if( papt .eq. 0.D+00 ) return
         ak    = rnosq/papt
         rnnsq = 0.D+00
         do 30 i=1,nxy
            x(i)   = x(i)+ak*p(i)
            rem(i) = rem(i)-ak*y(i)
            rnnsq  = rnnsq + rem(i)*rem(i)
   30    continue
         if( dsqrt(rnnsq/dble(nxy)) .le. eps ) return
         bk = rnnsq/rnosq
         do 40 i=1,nxy
            p(i) = rem(i)+bk*p(i)
   40    continue
   50 rnosq=rnnsq
         pause 'maximum iterations exceeded'
c      write(6,*) 'CNCG =', its, '  rnnsq=',rnnsq
      end

      subroutine atimex(x,y,ab,nxy,nx,ny)
c**********************************************************************
c
c     matrix multiplication y = A x                                 
c     
c     x         ->   vector of length nxy
c     ab        ->   0.5*dt/h/h/rey 
c     y         <-   vector of length nxy
c
c     A = discretization of the operator (I-dt/2/rey*grad^2)
c
c**********************************************************************
      integer nxy,nx,ny,nxm,nym,mxym
      include 'gridsz'
      double precision x(nxy),y(nxy),ab,roa
      common/density/roa(0:nxm+1,0:nym+1)
      integer i,j,k,k1,k2
      double precision tn,sv,c2,c3
      parameter (sv=7.D0, tn=10.D0, c2=2.D0, c3=-0.2D0)
c**********************************************************
c     interior points.
      k=nx+2
      do 10 j=2,ny-1
      do 15 i=2,nx-1
         y(k) = x(k) + ab*(4.D0*x(k) - x(k-1) - x(k+1) -
     +                     x(k-nx) - x(k+nx) )/roa(i,j)
         k = k+1
 15   continue
         k = k+2
 10   continue
c**********************************************************
c     boundary points.
      k1=2
      k2=nx*(ny-1)+2
      do 20 i=2,nx-1
         y(k1) = x(k1) + ab*(sv*x(k1) - x(k1-1) - x(k1+1) -
     +                       c2*x(k1+nx) - c3*x(k1+2*nx) )/
     +                       roa(i,1)
         y(k2) = x(k2) + ab*(sv*x(k2) - x(k2-1) - x(k2+1) -
     +                       c2*x(k2-nx) - c3*x(k2-2*nx) )/
     +                       roa(i,ny)
         k1=k1+1
         k2=k2+1
 20   continue
      k1=nx+1
      k2=2*nx
      do 30 j=2,ny-1
         y(k1) = x(k1) + ab*(sv*x(k1) - x(k1-nx) - x(k1+nx) -
     +                       c2*x(k1+1) - c3*x(k1+2))/
     +                       roa(1,j)
         y(k2) = x(k2) + ab*(sv*x(k2) - x(k2-nx) - x(k2+nx) -
     +                       c2*x(k2-1) - c3*x(k2-2))/
     +                       roa(nx,j)
         k1=k1+nx
         k2=k2+nx
 30   continue
c***********************************************************
c     corners
      k=1
      y(k) = x(k) + ab*(tn*x(k) - c2*(x(k+1)+x(k+nx)) -
     +                  c3*(x(k+2)+x(k+2*nx)) )/roa(1,1)
      k=nx
      y(k) = x(k) + ab*(tn*x(k) - c2*(x(k-1)+x(k+nx)) -
     +                  c3*(x(k-2)+x(k+2*nx)))/roa(nx,1)
      k=nx*(ny-1)+1
      y(k) = x(k) + ab*(tn*x(k) - c2*(x(k-nx)+x(k+1)) -
     +                  c3*(x(k-2*nx)+x(k+2)))/roa(1,ny)
      k=nx*ny
      y(k) = x(k) + ab*(tn*x(k) - c2*(x(k-1)+x(k-nx)) -
     +                  c3*(x(k-2)+x(k-2*nx)))/roa(nx,ny)
      end