      subroutine ctrest(lp,lev,origgrid)
      use amr
      type(GF3TL), pointer        :: origgrid,finegrid,coarsegrid
      integer lev
      type(LEVEL_INFO),dimension(maxlev) :: lp
      integer i,nx,ny,nz
      real idnorm,dr(3),dt,gflb(3),gfrb(3)
      integer, dimension(3) :: nr

      INTERFACE

      subroutine grid_step(grid,dt,dr,nr,gflb,gfrb)
      use user_param
      use u_interface
      use mesh 
      type(GF3TL), pointer :: grid
      real                 :: dt
      real, dimension(3)   :: dr,gflb,gfrb
      integer, dimension(3):: nr
      end subroutine grid_step

      subroutine expandcoarsegrid(coarsegrid,finegrid)
      use amr 
      type(GF3TL), pointer :: coarsegrid,finegrid
      integer i,ii,jj,kk
      end subroutine expandcoarsegrid  

      END INTERFACE
       write(6,*) 'start ctrest: +++++++++++++++++++++' 
      origgrid % gfcn(estre) % tlev(1) % data = 0.0

!     ASSUMES ENFORCED BOUNDARY CONDITIONS


      call copyGF3TL(origgrid,finegrid)
      call coarsecopyGF3TL(finegrid,coarsegrid)

      dr=lp(lev)%ldr
      dt=lp(lev)%ldt
! those two guys are still dummy
      gflb=0.0
      gfrb=0.0
! the following is not the ultimate solution either
! remember about refinement ratio of 2 (to be fixed)
      nr=finegrid%gfcn(1)%tlev(1)%ddim

      call grid_step(finegrid,dt,dr,nr,gflb,gfrb) 
      call grid_step(finegrid,dt,dr,nr,gflb,gfrb) 
      nr=finegrid%gfcn(1)%tlev(1)%ddim
      call grid_step(coarsegrid,2.0*dt,2.0*dr,nr/2,gflb,gfrb) 

      call expandcoarsegrid(coarsegrid,finegrid)
      do i=uxx,uzz
         nx = finegrid%gfcn(i)%tlev(1)%ddim(1) 
         nx = nx - (2 + mod(nx+1,2))
         ny = finegrid%gfcn(i)%tlev(1)%ddim(2) 
         ny = ny - (2 + mod(ny+1,2))
         nz = finegrid%gfcn(i)%tlev(1)%ddim(3) 
         nz = nz - (2 + mod(nz+1,2))
           
!        DIFFERENCE BTW FINE AND COARSE GOES IN TLEV 1

         finegrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,3:nz) =
     +      finegrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,3:nz) -
     +      finegrid%gfcn(i)%tlev(2)%data(3:nx,3:ny,3:nz) 
 
c        PUT NORMALIZATION IN TLEV 2

         finegrid%gfcn(i)%tlev(2)%data(3:nx,3:ny,3:nz) =
     +     (origgrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,3:nz) +
     +      origgrid%gfcn(i)%tlev(1)%data(2:nx-1,3:ny,3:nz) +
     +      origgrid%gfcn(i)%tlev(1)%data(4:nx+1,3:ny,3:nz) +
     +      origgrid%gfcn(i)%tlev(1)%data(3:nx,2:ny-1,3:nz) +
     +      origgrid%gfcn(i)%tlev(1)%data(3:nx,4:ny+1,3:nz) +
     +      origgrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,2:nz-1) +
     +      origgrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,4:nz+1) ) / 7.0

         idnorm = (SUM(ABS(origgrid%gfcn(i)%tlev(1)%data))/
     +   REAL(nx*ny*nz))/100.0

         WHERE (finegrid%gfcn(i)%tlev(2)%data < idnorm)
     +              finegrid%gfcn(i)%tlev(2)%data = idnorm

c        PUT ESTRE FOR EACH GFCN IN LEVEL 1

         finegrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,3:nz) =
     +      ABS(finegrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,3:nz) /
     +          finegrid%gfcn(i)%tlev(2)%data(3:nx,3:ny,3:nz) )

      enddo

      do i=qxx,qzz
         nx = finegrid%gfcn(i)%tlev(1)%ddim(1) 
         nx = nx - (2 + mod(nx+1,2))
         ny = finegrid%gfcn(i)%tlev(1)%ddim(2) 
         ny = ny - (2 + mod(ny+1,2))
         nz = finegrid%gfcn(i)%tlev(1)%ddim(3) 
         nz = nz - (2 + mod(nz+1,2))
           
c        DIFFERENCE BTW FINE AND COARSE GOES IN TLEV 1

         finegrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,3:nz) =
     +      finegrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,3:nz) -
     +      finegrid%gfcn(i)%tlev(2)%data(3:nx,3:ny,3:nz) 
 
c        PUT NORMALIZATION IN TLEV 2

         finegrid%gfcn(i)%tlev(2)%data(3:nx,3:ny,3:nz) =
     +     (origgrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,3:nz) +
     +      origgrid%gfcn(i)%tlev(1)%data(2:nx-1,3:ny,3:nz) +
     +      origgrid%gfcn(i)%tlev(1)%data(4:nx+1,3:ny,3:nz) +
     +      origgrid%gfcn(i)%tlev(1)%data(3:nx,2:ny-1,3:nz) +
     +      origgrid%gfcn(i)%tlev(1)%data(3:nx,4:ny+1,3:nz) +
     +      origgrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,2:nz-1) +
     +      origgrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,4:nz+1) ) / 7.0

         idnorm = (SUM(ABS(origgrid%gfcn(i)%tlev(1)%data))
     + /REAL(nx*ny*nz))/100.0

         WHERE (finegrid%gfcn(i)%tlev(2)%data < idnorm)
     +              finegrid%gfcn(i)%tlev(2)%data = idnorm

c        PUT ESTRE FOR EACH GFCN IN LEVEL 1

         finegrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,3:nz) =
     +      ABS(finegrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,3:nz) /
     +          finegrid%gfcn(i)%tlev(2)%data(3:nx,3:ny,3:nz) )

      enddo

c     ADD 'EM ALL UP, AND PUT 'EM IN ESTRE

      do i=uxx,uzz
         nx = finegrid%gfcn(i)%tlev(1)%ddim(1) 
         nx = nx - (2 + mod(nx+1,2))
         ny = finegrid%gfcn(i)%tlev(1)%ddim(2) 
         ny = ny - (2 + mod(ny+1,2))
         nz = finegrid%gfcn(i)%tlev(1)%ddim(3) 
         nz = nz - (2 + mod(nz+1,2))
         origgrid%gfcn(estre)%tlev(1)%data(3:nx,3:ny,3:nz) =
     +      origgrid%gfcn(estre)%tlev(1)%data(3:nx,3:ny,3:nz) +
     +      finegrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,3:nz)
      enddo

      do i=qxx,qzz
         nx = finegrid%gfcn(i)%tlev(1)%ddim(1) 
         nx = nx - (2 + mod(nx+1,2))
         ny = finegrid%gfcn(i)%tlev(1)%ddim(2) 
         ny = ny - (2 + mod(ny+1,2))
         nz = finegrid%gfcn(i)%tlev(1)%ddim(3) 
         nz = nz - (2 + mod(nz+1,2))
         origgrid%gfcn(estre)%tlev(1)%data(3:nx,3:ny,3:nz) =
     +      origgrid%gfcn(estre)%tlev(1)%data(3:nx,3:ny,3:nz) +
     +    finegrid%gfcn(i)%tlev(1)%data(3:nx,3:ny,3:nz)
      enddo

      origgrid%gfcn(estre)%tlev(1)%data = 
     &   origgrid%gfcn(estre)%tlev(1)%data / 12.0
 
      call deallocateGF3TL(finegrid)
      call deallocateGF3TL(coarsegrid)

      return
      end subroutine ctrest

c ***********************************************************************
c EXPANDCOARSEGRID 
c ***********************************************************************
c This routine expands (via extrapolation) 12 of the grid functions:
c    u's : 6 - 11
c    q's : 19 - 24
c The target grid is the lower timelevel (tlev = 2) of the fine grid
c ***********************************************************************

      subroutine expandcoarsegrid(coarsegrid,finegrid)
      use amr 
      type(GF3TL), pointer :: coarsegrid,finegrid
      integer i,ii,jj,kk

      do i = uxx,uzz
         do ii=1,coarsegrid % gfcn(i) % tlev(1) % ddim(1)
         do jj=1,coarsegrid % gfcn(i) % tlev(1) % ddim(2)
         do kk=1,coarsegrid % gfcn(i) % tlev(1) % ddim(3)
            finegrid%gfcn(i)%tlev(2)%data(2*ii-1,2*jj-1,2*kk-1) =
     +           coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk)
         enddo
         enddo
         enddo
      enddo
      do i = qxx,qzz
         do ii=1,coarsegrid % gfcn(i) % tlev(1) % ddim(1)
         do jj=1,coarsegrid % gfcn(i) % tlev(1) % ddim(2)
         do kk=1,coarsegrid % gfcn(i) % tlev(1) % ddim(3)
            finegrid%gfcn(i)%tlev(2)%data(2*ii-1,2*jj-1,2*kk-1) =
     +           coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk)
         enddo
         enddo
         enddo
      enddo

c     INTERPOLATION

      do i = uxx,uzz
         do ii=1,coarsegrid % gfcn(i) % tlev(1) % ddim(1) - 1 
         do jj=1,coarsegrid % gfcn(i) % tlev(1) % ddim(2) - 1
         do kk=1,coarsegrid % gfcn(i) % tlev(1) % ddim(3) - 1

c           X 
            finegrid%gfcn(i)%tlev(2)%data(2*ii,2*jj-1,2*kk-1) =
     +           (coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj,kk)) / 2.0

c           Y
            finegrid%gfcn(i)%tlev(2)%data(2*ii-1,2*jj,2*kk-1) =
     +           (coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj+1,kk)) / 2.0

c           Z
            finegrid%gfcn(i)%tlev(2)%data(2*ii-1,2*jj-1,2*kk) =
     +           (coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk+1)) / 2.0

c           X+Y
            finegrid%gfcn(i)%tlev(2)%data(2*ii,2*jj,2*kk-1) =
     +           (coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj+1,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj+1,kk) ) / 4.0
     
c           X+Z
            finegrid%gfcn(i)%tlev(2)%data(2*ii,2*jj-1,2*kk) =
     +           (coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk+1) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj,kk+1) ) / 4.0

c           Y+Z
            finegrid%gfcn(i)%tlev(2)%data(2*ii-1,2*jj,2*kk) =
     +           (coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj+1,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk+1) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj+1,kk+1) ) / 4.0

c           X+Y+Z
            finegrid%gfcn(i)%tlev(2)%data(2*ii,2*jj,2*kk) =
     +           (coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj+1,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk+1) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj+1,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj,kk+1) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj+1,kk+1) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj+1,kk+1) ) / 8.0

         enddo
         enddo
         enddo
      enddo
      
      do i = qxx,qzz
         do ii=1,coarsegrid % gfcn(i) % tlev(1) % ddim(1) - 1 
         do jj=1,coarsegrid % gfcn(i) % tlev(1) % ddim(2) - 1
         do kk=1,coarsegrid % gfcn(i) % tlev(1) % ddim(3) - 1

c           X 
            finegrid%gfcn(i)%tlev(2)%data(2*ii,2*jj-1,2*kk-1) =
     +           (coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj,kk)) / 2.0

c           Y
            finegrid%gfcn(i)%tlev(2)%data(2*ii-1,2*jj,2*kk-1) =
     +           (coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj+1,kk)) / 2.0

c           Z
            finegrid%gfcn(i)%tlev(2)%data(2*ii-1,2*jj-1,2*kk) =
     +           (coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk+1)) / 2.0

c           X+Y
            finegrid%gfcn(i)%tlev(2)%data(2*ii,2*jj,2*kk-1) =
     +           (coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj+1,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj+1,kk) ) / 4.0
     
c           X+Z
            finegrid%gfcn(i)%tlev(2)%data(2*ii,2*jj-1,2*kk) =
     +           (coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk+1) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj,kk+1) ) / 4.0

c           Y+Z
            finegrid%gfcn(i)%tlev(2)%data(2*ii-1,2*jj,2*kk) =
     +           (coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj+1,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk+1) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj+1,kk+1) ) / 4.0

c           X+Y+Z
            finegrid%gfcn(i)%tlev(2)%data(2*ii,2*jj,2*kk) =
     +           (coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj+1,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj,kk+1) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj+1,kk) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj,kk+1) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii,jj+1,kk+1) +
     +            coarsegrid%gfcn(i)%tlev(1)%data(ii+1,jj+1,kk+1) ) / 8.0

         enddo
         enddo
         enddo
      enddo

      return
      end subroutine expandcoarsegrid
      

















