c  Begin MG test 1

      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               S11(i,j,k) =  w3(i,j,k)
               S12(i,j,k) =  u(i,j,k)
               S13(i,j,k) =  v(i,j,k)
               S22(i,j,k) =  w(i,j,k)
               S23(i,j,k) = w1(i,j,k)
            enddo
         enddo
      enddo

c  End MG test 1


c Begin Test MG 2

      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               w3(i,j,k) =  S11(i,j,k)
               u(i,j,k) = S12(i,j,k)
               v(i,j,k) = S13(i,j,k) 
               w(i,j,k) = S22(i,j,k)
               p(i,j,k) = S23(i,j,k)
               relax_times = relax_times_old
            enddo
         enddo
      enddo
c***
c     check calls
c***
      if (t .le. 2) then
         write (6,*) "IN TEST MODE t = ",t," divergence sum = ", SUMFIELD(w3,nx,ny,nz)
      endif

      if (MOD(t,tcheck) .eq. 0) then
         write (*,*) "IN TEST MODE Norm2 div before = ",h*FINDN2S(w3,nx,ny,nz)
      endif 
c***
c     call the pressure solver
c***
      call mglin()
c***  
c     make the velocity field about divergence free
c***  
      do k=2,nz-1
         do j=2,ny-1
            do i=2,nx-1
               u(i,j,k) = u(i,j,k) - (p(i,j,k)-p(i-1,j,k))*a(i-1,j,k)
               v(i,j,k) = v(i,j,k) - (p(i,j,k)-p(i,j-1,k))*b(i,j-1,k)
               w(i,j,k) = w(i,j,k) - (p(i,j,k)-p(i,j,k-1))*c(i,j,k-1)
            enddo
         enddo
      enddo
c***
      call bc_vector(u,v,w,nx,ny,nz)
c***
c     check calls
c***
      if (MOD(t,tcheck) .eq. 0) then
         do k=2,nz-1
            do j=2,ny-1
               do i=2,nx-1
                  w3(i,j,k) = (u(i+1,j,k) - u(i,j,k) + v(i,j+1,k) -
     %                         v(i,j,k) + w(i,j,k+1) - w(i,j,k))*hi2
               enddo
            enddo
         enddo
         write (*,*) "IN TEST MODE Norm2 div after = ",
     %   h*FINDN2S(w3,nx,ny,nz)
      endif

c End MG Tetsting 2