!-----------------------------------------------------------------------
!
! (c) 2000-2001 Dale Choi
!
!      HISTORY: mg_main.F   (J. Gao's Unigrid MG solver)
!
!  Main Driver of MultiGrid Elliptic Solver
!  this code solves a poisson-type elliptic equation using LCS (Linear 
!  Correction Scheme), with octant symmentry and robin boundary condition
!
! (1) Tested in parADM Version 1
!
!   Input GF (gup,gam,Acoeff,\rho) and Output GF (\psi) defined at node-centers
!
! (2) Further modification to solve: 
!
!    gup11 * u,xx + gup22 * u,yy + gup33 * u,zz + acoeff * u^2 = rhs
!
! (3) Modification (being made) for AMR (non-uniform) grids
!
! Mar 14 2001  DONE: [1] reworked mg_gsrb_level.F / mg_gsrb_relax.F
!                    [2] OBC (mg_relaxb,mg_octant,mg_relaxb_robin)
!                    [3] check mg_lop.F (flux_matching)
!
! Thu Sep  6  2001   Stand Alone version (LCS) created
! Thu Sep 27  2001   2nd order convergence against analytic data verified!
!
! Sat Sep 22 16:08:38 EDT 2001
! FAS (Full Approximation Storage scheme) -- implemented & being tested now
!
!-----------------------------------------------------------------------
      program mg_main
  
      implicit none

#include "physicaldata.fh"

! Use this for PARAMESH V.1
!     include 'mpp/shmem.fh'
! Use this for PARAMESH V.2
#include "amr_shmem.fh"

      include 'tree.fh'
      include 'shmem_reduce.fh'
      include 'mg_data.fh'
      include 'mg_node.fh'
      include 'mg_oct.fh'

      integer mype, i, ir, iv, level,  ix,iy,iz,lb,  nprocs, maxprocs
      integer nrelax_up
      integer shmem_my_pe, shmem_n_pes
      save mype, nprocs

!-----------------------------------------------------------------------
!     initial paramesh
!-----------------------------------------------------------------------
!     call comm_start(maxprocs,nprocs,mype)
      mype  = shmem_my_pe()
      nprocs= shmem_n_pes()

      call amr_initialize
      call shmem_barrier_all()

!-----------------------------------------------------------------------
! read-in from parameter file and set up initial refinement mesh
!-----------------------------------------------------------------------
      call mg_init_param(mype)
      call shmem_barrier_all()
      call mg_init_grid(mype)
      call shmem_barrier_all()

!-----------------------------------------------------------------------
!     set up initial data
!          set RHS, RES here! (calls mg_lop_nf)
!-----------------------------------------------------------------------
      call mg_init_soln(lrefine_max) 
      call shmem_barrier_all()

!-----------------------------------------------------------------------
! Special treatment for LMAX = 1 (It *is* a simple relaxation!)
!-----------------------------------------------------------------------
      if( lrefine_max .eq. 1 ) then

! redistribute blocks among processors
       call amr_morton_order(nprocs,lnblocks)
       write(0,*)'mype,lnblocks = ',mype,lnblocks

       if( mg_type .eq. 0 .or. mg_type .eq. 2 ) then

        do i = 1, ns

          do ir = 0, 1
            call shmem_barrier_all()
            call mg_guardcell(iu,1)
            call shmem_barrier_all() !sep26 2001
#if OCTANT_TEST == 0
            call mg_octant(iu,1)
#else if OCTANT_TEST == 1
            call mg_octant2(iu,1)
#endif
            call shmem_barrier_all() !sep26 2001
            call mg_relax_swid(ir,1) !swid
          end do

          do ir = 0, 1
            call shmem_barrier_all()
            call mg_guardcell(iu,1)
            call shmem_barrier_all() !sep26 2001
#if OCTANT_TEST == 0
            call mg_octant(iu,1)
#else if OCTANT_TEST == 1
            call mg_octant2(iu,1)
#endif
           !call mg_relaxb(ir,1) ! original (DONOT USE!)
            call mg_fillc(itmp,1.0e0,1)        
            call shmem_barrier_all() !sep26 2001
            call mg_relaxb_robin(iu,itmp,ir,1) !to control 'rhs' for robin boundary eqn.
          end do

          call shmem_barrier_all()
         !call mg_comp_res(itmp,1,i,0)   ! use itmp as temp. storage
          call mg_comp_res(itmp,1,-100,i)

        end do

       else if( mg_type .eq. 1 ) then

        do i=1,ns

          do ir=0,1
            call shmem_barrier_all()
            call mg_guardcell1(iu,1) ! apr 3
            call mg_octant(iu,1)
            call mg_lapl_0(iu,irhs,1,ir) ! simple relaxation
          end do

          do ir=0,1
            call shmem_barrier_all()
            call mg_guardcell1(iu,1) ! apr 3
            call mg_octant(iu,1)
            call mg_relaxb_robin(iu,irhs,ir,1) 
          end do

          call shmem_barrier_all()     ! Mar 20, 2001
          call mg_comp_res(itmp,1,i,0)   ! use itmp as temp. storage

        enddo

       end if

       call mg_dump_1lev(iu,1,'u_')
       call shmem_barrier_all()

       stop

      endif

!-----------------------------------------------------------------------
!
!     V-Cycle                nvcycle := # of V-Cycle
!
!-----------------------------------------------------------------------
      do iv = 1, nvcycle

        if( mype .eq. 0 ) write(0,*)'V-cycle #= ',iv
        if( mype .eq. 0 ) write(59,*)'V-cycle #= ',iv
       !----------------------------------------------------------------
       ! Down leg of V-cycle
       !----------------------------------------------------------------
        do level = lrefine_max, 2, -1
          if( mype .eq. 0 ) write(0,*)' Down Leg: level = ',level
          if( mype .eq. 0 ) write(59,*)' Down Leg: level = ',level

         !--------------------------------------------------
         ! Relaxation on 'level'       nrelax := # of relaxation
         !--------------------------------------------------
         if( mg_type .eq. 0 ) then
          do i = 1, nrelax

            do ir=0,1
              call shmem_barrier_all()
              call mg_guardcell(iu,level)
              call shmem_barrier_all() !sep26 2001
#if OCTANT_TEST == 0
              call mg_octant(iu,level)
#else if OCTANT_TEST == 1
              call mg_octant2(iu,level)
#endif
              call shmem_barrier_all() !sep26 2001
              call mg_relax_swid(ir,level)
            end do

            do ir=0,1
              call shmem_barrier_all()
              call mg_guardcell(iu,level)
#if OCTANT_TEST == 0
              call mg_octant(iu,level)
#else if OCTANT_TEST == 1
              call mg_octant2(iu,level)
#endif
             !call mg_relaxb(ir,level) !original (DONOT USE!)
              call mg_fill0(itmp,level)                   !apr 21 2001
              if(level.eq.lrefine_max) 
     &            call mg_fillc(itmp,1.0e0,level)         !apr 21 2001
              call shmem_barrier_all() !sep26 2001
              call mg_relaxb_robin(iu,itmp,ir,level)      !apr 21 2001
            end do

            call shmem_barrier_all() !sep26 2001
            call mg_comp_res(itmp,level,0,iv)
            call shmem_barrier_all() !sep26 2001

          end do

         else if( mg_type .eq. 1 ) then

          call shmem_barrier_all()      ! Feb 27 2001
          call mg_copy(iusave,iu,level)

          call mg_guardcell1(iusave,level) ! Apr 10 2001
          call mg_octant(iusave,level)    ! Mar 14 2001
          call set_interface_from_c(iusave,level) !Apr 2 2001

          call mg_fill0(ierr,level-1)

          call mg_guardcell1(ires,level) ! Apr 10 2001
          call mg_octant(ires,level)    ! Mar 5 2001
          call set_interface_from_c(ires,level) !Apr 2 2001
          call shmem_barrier_all()      ! Mar 14 2001

          call mg_gsrb_level(ierr,level)
          call shmem_barrier_all()      ! Mar 14 2001

          call mg_guardcell1(ierr,level) ! Apr 10 2001
          call mg_octant(ierr,level)    ! Mar 5 2001
          call set_interface_from_c(ierr,level) !Apr 2 2001
          call shmem_barrier_all()      ! Mar 14 2001

          call mg_add(iu,iu,ierr,level) !ok
          call shmem_barrier_all()      ! Mar 14 2001

          call mg_guardcell1(iu,level)  ! Apr 10 2001
          call mg_octant(iu,level)     ! Mar 5 2001
          call set_interface_from_c(iu,level) !Apr 2 2001

         else if( mg_type .eq. 2 ) then

          do i = 1, nrelax

            do ir=0,1
              call shmem_barrier_all()
              call mg_guardcell(iu,level)
#if OCTANT_TEST == 0
              call mg_octant(iu,level)
#else if OCTANT_TEST == 1
              call mg_octant2(iu,level)
#endif
              call mg_relax_swid(ir,level)
            end do

            do ir=0,1
              call shmem_barrier_all()
              call mg_guardcell(iu,level)
#if OCTANT_TEST == 0
              call mg_octant(iu,level)
#else if OCTANT_TEST == 1
              call mg_octant2(iu,level)
#endif
              call mg_fillc(itmp,1.0e0,level)
              call mg_relaxb_robin(iu,itmp,ir,level)
            end do

            call mg_comp_res(itmp,level,0,iv)

          end do

         end if

         call shmem_barrier_all()

         !--------------------------------------------------
         ! Set up the coarse grid problem on 'level-1'
         !--------------------------------------------------
         if( mg_type .eq. 0 ) then

          call mg_guardcell(iu,level)
          call shmem_barrier_all()

#if OCTANT_TEST == 0
          call mg_octant(iu,level)
#else if OCTANT_TEST == 1
          call mg_octant2(iu,level)
#endif
          call shmem_barrier_all() !sep26 2001
          call mg_residue_swid(level)      ! Modified for swid(011901) & AMR(062101)
            !check mg_residue_swid again (might want to use ivar to compue residue?!?!
!         call set_interface_from_c(ires,level) !Jul 9, 2001

          call shmem_barrier_all()

          call mg_guardcell(ires,level)

#if OCTANT_TEST == 0
          call mg_octant(ires,level)
#else if OCTANT_TEST == 1
          call mg_octant2(ires,level)
#endif
          call shmem_barrier_all()

          call mg_restrict(ires,level-1)
          call shmem_barrier_all() ! sep26 2001
          call mg_guardcell(ires,level-1) !? need guardcell1 call for level-1 data? (Jun 21 2001)

!        !---------------------------------------------------
!        ! trying to restrict coefficients
!        !---------------------------------------------------
!         call mg_restrict(iacoeff,level-1) ! added on 020101
!         call mg_restrict(igup11c,level-1) ! added on 020101
!         call mg_restrict(igup22c,level-1) ! added on 020101
!         call mg_restrict(igup33c,level-1) ! added on 020101
! ok, but restriction is done only for non-guardcells (Jun 19, 2001)
! also over-write some data that causes problems!
!         ! Apr 19 2001   guardcell call really helps!
!         call mg_guardcell(iacoeff,level-1)  !no '1's here (jun 20, 2001)
!         call mg_guardcell(igup11c,level-1)
!         call mg_guardcell(igup22c,level-1)
!         call mg_guardcell(igup33c,level-1)
#if OCTANT_TEST == 0 
!         call mg_octant(iacoeff,level-1)
!         call mg_octant(igup11c,level-1)
!         call mg_octant(igup22c,level-1)
!         call mg_octant(igup33c,level-1)
#else if OCTANT_TEST == 1 
!         call mg_octant2(iacoeff,level-1)
!         call mg_octant2(igup11c,level-1)
!         call mg_octant2(igup22c,level-1)
!         call mg_octant2(igup33c,level-1)
#endif

          call mg_copy(irhs,ires,level-1)
          call mg_fill0(iu,level-1)

         else if( mg_type .eq. 1 ) then

          call mg_residual(level-1)       ! to test (mg_lop) only for AMR

          call shmem_barrier_all()        

          call mg_restrict(iacoeff,level-1) ! 022201
          call mg_restrict(igup11c,level-1) ! 022201
          call mg_restrict(igup22c,level-1) ! 022201
          call mg_restrict(igup33c,level-1) ! 022201
          call mg_restrict(irhs,level-1)    ! 022201
          call shmem_barrier_all()
          call mg_guardcell1(iacoeff,level-1) ! Mar 19 2001
          call mg_guardcell1(igup11c,level-1)
          call mg_guardcell1(igup22c,level-1)
          call mg_guardcell1(igup33c,level-1)
          call mg_guardcell1(irhs,level-1)
          call mg_octant(iacoeff,level-1) ! Mar 5 2001
          call mg_octant(igup11c,level-1) ! Mar 5 2001
          call mg_octant(igup22c,level-1) ! Mar 5 2001
          call mg_octant(igup33c,level-1) ! Mar 5 2001
          call mg_octant(irhs,level-1) ! Mar 5 2001

         else if( mg_type .eq. 2 ) then

!-----------------------------------------------------------
!           Compute truncation error estimate.
!
!               2h      2h  2h  h     2h  h  h
!           tau    :=  L   I   u  -  I   L  u
!               h           h         h
!
!           this implementation uses three temporary arrays for clarity (:-));
!           production implementation should economize.
!-----------------------------------------------------------
! restriction to compute I^{2h}_{h} u^{h} part
! u(l-1) contains I^{2h}_{h} u^{h}

          call shmem_barrier_all()
!         call mg_restrict(iu,level-1) ! dale's  this gives increase of residual vs. # of relax
          call inject(mype,iu,level-1) ! jimin's inject gives decrease of resi vs. # of relax
          call shmem_barrier_all()
          call mg_guardcell(iu,level-1)
          call mg_octant(iu,level-1)

! compute L^{2h} (I^{2h}_{h} u^{h}) from u(l-1) = I^{2h}_{h} u^{h}, 
! and store it in tmp(l-1)

          call shmem_barrier_all()
          call mg_fill0(itmp,level-1)
          call mg_lop2(iu,level-1)
          call shmem_barrier_all()
          call mg_guardcell(itmp,level-1)
          call mg_octant(itmp,level-1)

! compute L^{h} u^{h} part and store in tmp(l)

          call shmem_barrier_all()
          call mg_fill0(itmp,level)
          call mg_lop2(iu,level)
          call shmem_barrier_all()
          call mg_guardcell(itmp,level)
          call mg_octant(itmp,level)
!         call mg_fill_bnd(itmp,level,0.0e0)

! restriction of L^{h} u^{h} and store in tau(l-1)= I^{2h}_{h}(L^{h} u^{h})

          call shmem_barrier_all()
          call mg_fill0(itau,level-1)
          call mg_restrict2(itau,itmp,level-1)
          call shmem_barrier_all()
          call mg_guardcell(itau,level-1)
          call mg_octant(itau,level-1)

! compute truncation error and store in q(tau(l-1))

          call shmem_barrier_all()
          call mg_sub(itau,itmp,itau,level-1)
          call shmem_barrier_all()
          call mg_guardcell(itau,level-1)
          call mg_octant(itau,level-1)

! at this point tau contains truncation error and tmp L^{h} u^{h}
!         call mg_fill_bnd(itau,level-1,0.0e0)
          call mg_anorm(mype,itau,level-1)
          call mg_dump_1lev(itau,level-1,'lte_')

!-----------------------------------------------------------
!           Define right hand side and initialize unknown
!           for coarse grid problem.
!
!               2h       2h    h       2h
!           rhs     :=  I   rhs   + tau
!                        h              h
!
!            2h         2h  h
!           u       := I   u
!                       h
!-----------------------------------------------------------
! restrict rhs

          call shmem_barrier_all()
          call mg_restrict(irhs,level-1)
          call shmem_barrier_all()
          call mg_guardcell(irhs,level-1)
          call mg_octant(irhs,level-1)

! add truncation error to the rhs

          call shmem_barrier_all()
          call mg_add(irhs,irhs,itau,level-1)
          call shmem_barrier_all()        
          call mg_guardcell(irhs,level-1)
          call mg_octant(irhs,level-1)
!         call mg_fill_bnd(irhs,level-1,0.0e0)
          call mg_dump_1lev(irhs,level-1,'rhs_r')

! initial guess for the coarse grid problem
! restriction from the fine grid solution

          call shmem_barrier_all()        
          call mg_restrict(iu,level-1)
          call shmem_barrier_all()        
          call mg_guardcell(iu,level-1)
          call mg_octant(iu,level-1)

          call mg_copy(iuold,iu,level-1) ! jimin's

         endif

         call shmem_barrier_all()        

        end do
        
       !----------------------------------------------------------------
       ! solve on the coarest grid
       !----------------------------------------------------------------
       if( mg_type .eq. 0 ) then

        do i=1,ns

          do ir=0,1
            call shmem_barrier_all()   
            call mg_guardcell1(iu,1) !apr 19 2001
#if OCTANT_TEST == 0
            call mg_octant(iu,1)
#else if OCTANT_TEST == 1
            call mg_octant2(iu,1)
#endif
            call mg_relax_swid(ir,1) !swid
          end do

          do ir=0,1
            call shmem_barrier_all()   
            call mg_guardcell1(iu,1) !apr 19 2001
#if OCTANT_TEST == 0
            call mg_octant(iu,1)
#else if OCTANT_TEST == 1
            call mg_octant2(iu,1)
#endif
           !call mg_relaxb(ir,1) !original (DONOT USE!)
            call mg_fill0(itmp,1)               !apr 21 2001
            call mg_relaxb_robin(iu,itmp,ir,1)  !apr 21 2001
          end do

          call shmem_barrier_all()
          call mg_comp_res(itmp,1,i,iv)   ! use itmp as temp. storage

        end do

       else if( mg_type .eq. 1 ) then

        call mg_guardcell1(ires,1) ! apr 10
        call mg_octant(ires,1)    ! NEW Mar 14
        call set_interface_from_c(ires,level) !Apr 2 2001
        call shmem_barrier_all()  ! NEW Mar 14
       !call mg_dump_1lev(ierr,1,'error_')

        do i=1,ns

          do ir=0,1
            call shmem_barrier_all()   
            call mg_guardcell1(ierr,1) ! Apr 20 2001
            call mg_octant(ierr,1)
!           call set_interface_from_c(ierr,1) !Apr 10 2001
            call mg_lapl_0(ierr,ires,1,ir) ! simple relaxation
          end do

          do ir=0,1
            call shmem_barrier_all()
            call mg_guardcell1(ierr,1) ! Apr 10 2001
            call mg_octant(ierr,1)
!           call set_interface_from_c(ierr,1) !Apr 10 2001
            call mg_relaxb_robin(ierr,ires,ir,1) ! NEW Feb 27
          end do

          call shmem_barrier_all()
          call mg_add(iu,iu,ierr,1)   ! what iu??
          call shmem_barrier_all()
          call mg_comp_res(itmp,1,i,iv)  !can use itmp
!         call mg_anorm(mype,ierr,1)  ! apr 11 2001
          call shmem_barrier_all()
          call mg_sub(iu,iu,ierr,1)

        enddo

        call shmem_barrier_all()
        call mg_guardcell1(ierr,1) ! NEW Apr 10
        call mg_octant(ierr,1)    ! NEw Mar 14
        call shmem_barrier_all()

       else if( mg_type .eq. 2 ) then

        do i= 1, ns

          do ir=0,1
            call shmem_barrier_all()
            call mg_guardcell(iu,1) 
#if OCTANT_TEST == 0
            call mg_octant(iu,1)
#else if OCTANT_TEST == 1
            call mg_octant2(iu,1)
#endif
            call mg_relax_swid(ir,1)
          end do

          do ir=0,1
            call shmem_barrier_all()
            call mg_guardcell(iu,1)
#if OCTANT_TEST == 0
            call mg_octant(iu,1)
#else if OCTANT_TEST == 1
            call mg_octant2(iu,1)
#endif
            call mg_fillc(itmp,1.0e0,1)
            call mg_relaxb_robin(iu,itmp,ir,1)
          end do

          call shmem_barrier_all()
          call mg_guardcell(iu,1) ! NEW Oct 3 2001
          call mg_octant(iu,1)    ! NEw Oct 3 2001
          call mg_comp_res(itmp,1,i,iv)   ! use itmp as temp. storage

        end do

       end if

       call shmem_barrier_all()

       !----------------------------------------------------------------
       ! Up leg of V-cycle
       !----------------------------------------------------------------
        do level = 2, lrefine_max
         if( mype .eq. 0 ) write(0,*)' Up leg: level =', level
         if( mype .eq. 0 ) write(59,*)' Up leg: level =',level

         if( mg_type .eq. 0 ) then
          call shmem_barrier_all() !sep26 2001
          call mg_copy(itmp,iu,level-1)

          call shmem_barrier_all()

          call mg_prolong(mype,itmp,level)  !check this again!
          call shmem_barrier_all() !sep26 2001
          call mg_add(iu,iu,itmp,level)
!         call set_interface_from_c(iu,level) !Jul 9, 2001

          call shmem_barrier_all()        

          nrelax_up = nrelax
          if( i .eq. lrefine_max ) nrelax_up = nrelax_1
          do i = 1, nrelax_up

            do ir=0,1
              call shmem_barrier_all()
              call mg_guardcell(iu,level)
#if OCTANT_TEST == 0
              call mg_octant(iu,level)
#else if OCTANT_TEST == 1
              call mg_octant2(iu,level)
#endif
              call mg_relax_swid(ir,level) !swid
            end do

            do ir=0,1
              call shmem_barrier_all()
              call mg_guardcell(iu,level)
#if OCTANT_TEST == 0
              call mg_octant(iu,level)
#else if OCTANT_TEST == 1
              call mg_octant2(iu,level)
#endif
             !call mg_relaxb(ir,level) !original (DONOT USE!)
              call mg_fill0(itmp,level)              !apr 21 2001
              if(level.eq.lrefine_max) 
     &           call mg_fillc(itmp,1.0e0,level)     !apr 21 2001
              call mg_relaxb_robin(iu,itmp,ir,level) !apr 21 2001
            end do

            call mg_comp_res(itmp,level,0,iv)

          end do  

         else if( mg_type .eq. 1 ) then

          call shmem_barrier_all()

          call mg_guardcell1(ierr,level-1)        ! NEW Apr 10
          call mg_octant(ierr,level-1)            ! NEW Mar 14

          call shmem_barrier_all()

          call mg_prolong1(mype,itmp,ierr,level) !itmp=interp(err{l-1})

          call shmem_barrier_all()

          call mg_guardcell1(itmp,level)        ! NEW Apr 10
          call mg_octant(itmp,level)            ! NEW Mar 5

          call shmem_barrier_all()

          call mg_add(ierr,ierr,itmp,level)     !ierr=ierr+itmp

          call shmem_barrier_all()

          call mg_guardcell1(ierr,level)        ! NEW Apr 10
          call mg_octant(ierr,level)            ! NEW Mar 5

! should we do smoother here for ierr?? apr 11 2001
          call mg_guardcell1(ires,level)
          call mg_octant(ires,level)
          call mg_gsrb_level(ierr,level)        ! apr 11 2001

          call shmem_barrier_all()

          call mg_lop_nf(ierr,level)            !itmp=lop_nf(ierr)

          call shmem_barrier_all()

          call mg_guardcell1(itmp,level)        ! NEW Apr 10
          call mg_octant(itmp,level)            ! NEW Mar 14

          call shmem_barrier_all()

          call mg_sub(ires,ires,itmp,level)     !ires=ires-lop_nf

          call shmem_barrier_all()

          call mg_guardcell1(ires,level)        ! NEW Apr 10
          call mg_octant(ires,level)            ! NEW Mar 5

          call mg_fill0(iderr,level)
!         call mg_fill0(iderr,level-1)          ! NEW Apr 10

          call shmem_barrier_all()

          call mg_gsrb_level(iderr,level)

          call shmem_barrier_all()

          call mg_guardcell1(iderr,level)        ! apr 11
          call mg_octant(iderr,level)            ! Mar 14

          call shmem_barrier_all()

          call mg_add(ierr,ierr,iderr,level)

          call shmem_barrier_all()

          call mg_guardcell1(ierr,level)        ! apr 11
          call mg_octant(ierr,level)            ! Mar 14

          call shmem_barrier_all()

          call mg_guardcell1(iusave,level)        ! apr 11
          call mg_octant(iusave,level)            ! Mar 14

          call shmem_barrier_all()

          call mg_add(iu,iusave,ierr,level)

          call shmem_barrier_all()

          call mg_guardcell1(iu,level)        ! apr 11
          call mg_octant(iu,level)            ! NEW Mar 5

          call mg_comp_res(itmp,level,0,iv)

         else if( mg_type .eq. 2 ) then

c-----------------------------------------------------------
c           Linearly interpolate correction and update
c           level m unknown.
c
c            h       h    h    2h     2h  h
c           u    := u  + I  ( u   -  I   u  )
c                         2h          h
c-----------------------------------------------------------

! restrict previous fine grid solution and store in q(tau(l-1)) = I^{2h}_{h} u^{h}

           call shmem_barrier_all()
           call mg_fill0(itau,level-1)
           call mg_restrict2(itau,iu,level-1)
           call shmem_barrier_all()
           call mg_guardcell(itau,level-1)
           call mg_octant(itau,level-1)

! compute u^{2h} - I^{2h}_{h} u^{h} and store in q(tau(l-1))

           call shmem_barrier_all()
           call mg_sub(itau,iu,itau,level-1) !dale's
!          call mg_sub(itau,iu,iuold,level-1) !jimin's
           call shmem_barrier_all()
           call mg_guardcell(itau,level-1)
           call mg_octant(itau,level-1)

! interpolate to get correction to the previous solution; store in q(tmp(l))

           call shmem_barrier_all()
           call mg_fill0(itmp,level)
           call mg_prolong1(mype,itmp,itau,level)
           call shmem_barrier_all()
           call mg_guardcell(itmp,level)
           call mg_octant(itmp,level)

! add the correction to the previous solution

           call shmem_barrier_all()
           call mg_add(iu,iu,itmp,level)
           call mg_guardcell(iu,level)
           call mg_octant(iu,level)

! relax to remove high freq. components
           do i = 1, nrelax

            do ir=0,1
              call shmem_barrier_all()
              call mg_guardcell(iu,level)
#if OCTANT_TEST == 0
              call mg_octant(iu,level)
#else if OCTANT_TEST == 1
              call mg_octant2(iu,level)
#endif
              call mg_relax_swid(ir,level)
            end do

            do ir=0,1
              call shmem_barrier_all()
              call mg_guardcell(iu,level)
#if OCTANT_TEST == 0
              call mg_octant(iu,level)
#else if OCTANT_TEST == 1
              call mg_octant2(iu,level)
#endif
              call mg_guardcell(iu,level) !oct 3 2001
              call mg_fillc(itmp,1.0e0,level)
              call mg_relaxb_robin(iu,itmp,ir,level)
              call mg_guardcell(iu,level) !oct 3 2001
            end do

            call mg_comp_res(itmp,level,0,iv)

           end do

         end if ! mg_type

        end do !up leg of V-cycle

        call shmem_barrier_all() !sep26 2001
        call mg_comp_res(itmp,lrefine_max,-100,iv)

      end do !of V-Cycle

      call mg_dump_1lev(iu,lrefine_max,'u_')

      call mg_test(itmp,lrefine_max)
      call mg_dump_1lev(itmp,lrefine_max,'uan_')
      call mg_dump_1lev(irhs,lrefine_max,'rhs_')

      call mg_comp_admmass(iu,lrefine_max,0)
      call shmem_barrier_all()

      stop
      end