!----------------------------------------------------------------------
!  This routine updates the following grid functions
!  u 
!----------------------------------------------------------------------
      subroutine update0(unp1,un,rhs,g1_Nx,g1_Ny,t,x,y,dt,dx,dy,mgeps,
     &  mgmaxcyc)
        implicit none

        include 'globals.inc'

        integer g1_Nx,g1_Ny
        real*8  unp1(1:g1_Nx,1:g1_Ny)
        real*8  un(1:g1_Nx,1:g1_Ny)
        real*8  rhs(1:g1_Nx,1:g1_Ny)
        real*8  t
        real*8  x(*)
        real*8  y(*)
        real*8  dt
        real*8  dx
        real*8  dy
        real*8  mgeps
        integer mgmaxcyc


      include    'mgcom.inc'
      real*8      dmnrm2

      integer     nx,   ny

      integer     sizeq
      parameter ( sizeq = 10 000 000 ) 
      real*8      q(sizeq)

      integer     nxmax
      parameter ( nxmax = 4097 ) 
      integer     nymax
      parameter ( nymax = 4097 ) 

      real*8      mg2v0
      real*8      mgresult,  tloc
      integer     ncycle,  icyc

      integer     i,      j

      logical     ltrace
      parameter ( ltrace = .true. )

      logical     vstrace
       parameter ( vstrace = .false. )

      integer     ushape(2)

      data        tloc / 0.0d0 /
      
      nx = g1_nx
      ny = g1_ny
      ushape(1) = nx
      ushape(2) = ny

      if( nx .gt. nxmax ) then
         write(0,*) 'update0: Insufficient internal storage'
         write(0,*) '   nx = ', nx, '  >  nxmax = ', nxmax
         stop
      end if
      if( ny .gt. nymax ) then
         write(0,*) 'update0: Insufficient internal storage'
         write(0,*) '   ny = ', ny, '  >  nymax = ', nymax
         stop
      end if

c-----------------------------------------------------------------------
c     Set common variables for communication with multigrid routines ...
c-----------------------------------------------------------------------
      cdt = dt
     
      if( ltrace ) then
         write(0,*) 'update0: nx =', nx
         write(0,*) 'update0: ny =', ny
         write(0,*) 'update0: dt =', dt
         write(0,*) 'update0: dx =', dx
         write(0,*) 'update0: |u| =', dmnrm2(un,nx,ny)
      end if
c-----------------------------------------------------------------------
c     Set up rhs vector and initialize solution ...
c-----------------------------------------------------------------------
      do j = 2 , ny - 1
         do i = 2 , nx - 1
              rhs(i,j) = un(i,j) + 0.5d0 * cdt * 
     &               ( un(i+1,j) - 2.0d0 * un(i,j) + un(i-1,j) +
     &                 un(i,j+1) - 2.0d0 * un(i,j) + un(i,j-1)
     &               ) / dx**2
            unp1(i,j) = un(i,j)
         end do
      end do

      call dmlbs(unp1,0.0d0,nx,ny)
      call dmlbs(rhs,0.0d0,nx,ny)


      if( ltrace ) then
         write(0,*) 'update0: |unp1| =', dmnrm2(unp1,nx,ny)
         write(0,*) 'update0: |urhs| =', dmnrm2(rhs,nx,ny)
      end if
c-----------------------------------------------------------------------
c     ... and invoke multigrid solver
c-----------------------------------------------------------------------
c     real*8 function mg2v0
c     &   (uf,rhsf,nxf,nyf,     hf, q, sizeq,
c     &    preswf,pstswf,presw,pstsw,ncycle)

      ncycle = 0 
      mgresult = 2.0d0 * mgeps
      call cres(q,unp1,rhs,2,ushape,dx)

      tloc = tloc + 1.0d0
       do while( mgresult .gt. mgeps )
c-----------------------------------------------------------------------
c        For this problem, and for lambda = 0.05, 2 pre-CGC and 2 
c        post-CGC sweeps seem to relatively optimal. 
c-----------------------------------------------------------------------
         mgresult = mg2v0(unp1,rhs,nx,ny,dx,q,sizeq,2,2,2,2,ncycle)
         call cres(q,unp1,rhs,1,nx,dx)
         if( ncycle .gt. mgmaxcyc ) then
            write(0,*) 'update0: ncycle ', ncycle, ' > mgmaxcyc ',
     &                 mgmaxcyc
            stop
         end if
         tloc = tloc + 1.0d0
         if( ltrace ) then 
            write(0,*) 'update0: ncycle =  ', ncycle
            write(0,*) 'update0: mgresult = ', mgresult
            write(0,*) 'update0: mgeps = ', mgeps
            write(0,*) 'update0: |res| = ', dmnrm2(q,nx,ny)
         end if
      end do

      if( ltrace ) then
         write(0,*) 'update0: mg2v0 returns ', mgresult
      end if
        return
      end

!----------------------------------------------------------------------
!  This routine updates the following grid functions
!  ures 
!----------------------------------------------------------------------
      subroutine update1(u_np1,u_n,ures_n,g1_Nx,g1_Ny,dt,dx,dy)
        implicit none

        include 'globals.inc'

        integer g1_Nx,g1_Ny
        real*8  u_np1(1:g1_Nx,1:g1_Ny)
        real*8  u_n(1:g1_Nx,1:g1_Ny)
        real*8  ures_n(1:g1_Nx,1:g1_Ny)
        real*8  dt
        real*8  dx
        real*8  dy
        integer i,j,k

          do j=2, g1_Ny-1, 1
          do i=2, g1_Nx-1, 1
          ures_n(i,j)=ures_n(i,j)-(ures_n(i,j)-((u_np1(i,j)-u_n(i,j))/dt
     &      -((u_np1(i+1,j)-2*u_np1(i,j)+u_np1(i-1,j))/dx**2+(u_n(i+1,j)
     &      -2*u_n(i,j)+u_n(i-1,j))/dx**2)/2-((u_np1(i,j+1)-2*u_np1(i,j)
     &      +u_np1(i,j-1))/dy**2+(u_n(i,j+1)-2*u_n(i,j)+u_n(i,j-1))/dy**
     &      2)/2))
          end do
          end do
          do j=1, g1_Ny, 1
          i=1
          ures_n(i,j)=ures_n(i,j)-ures_n(i,j)
          end do
          do j=1, g1_Ny, 1
          i=g1_Nx
          ures_n(i,j)=ures_n(i,j)-ures_n(i,j)
          end do
          j=1
          do i=1, g1_Nx, 1
          ures_n(i,j)=ures_n(i,j)-ures_n(i,j)
          end do
          j=g1_Ny
          do i=1, g1_Nx, 1
          ures_n(i,j)=ures_n(i,j)-ures_n(i,j)
          end do
        return
      end