       subroutine ctrest(lp,lev,p)
       use amr
       use u_interface
       type(LEVEL_INFO), dimension(maxlev) :: lp
       integer lev
       type(GF3TL), pointer :: p

        integer fstep

        type(GF3TL), pointer :: gfine,gcoars
        real :: dt,dr,dr2,dt2,idnrm
        integer :: nr,nr2,iflag


c       -----------------------------------------------------------
c       CREATE TWO COPIES OF THE GRID
c       -----------------------------------------------------------
!      ????????
       iflag=1

        allocate(gfine,stat=ierr)
        if(ierr>0) then
          write(6,*) 'CTREST(1): not enough memory'
          stop
        endif 

        allocate(gcoars,stat=ierr)
        if(ierr>0) then
          write(6,*) 'CTREST(2): not enough memory'
          stop
        endif 

        call copyGF3TL(p,gcoars)
        call copyGF3TL(p,gfine)

         write(6,*) gcoars%gfcn(w1)%tlev(1)%data(:,1,1)
         write(6,*) 'end.',nr

        dt=lp(lev)%ldt
        dr=lp(lev)%ldr(1)
        nx=p%gfcn(1)%tlev(1)%c(1)%n_coord
        ny=p%gfcn(1)%tlev(1)%c(2)%n_coord
        nz=p%gfcn(1)%tlev(1)%c(3)%n_coord
        nr=nx
        dt2=dt
        dr2=dr
        nr2=nr
c       -----------------------------------------------------------
c       COARSE STEP
c       -----------------------------------------------------------
        call elfppi(
     &  gcoars%gfcn(PPN)%tlev(1)%data(:,1,1),       !PPNP12
     &  p%gfcn(PPN)%tlev(1)%data(:,1,1),            !PPNP1
     &  p%gfcn(PPN)%tlev(2)%data(:,1,1),            !PPN
     &  p%gfcn(PIN)%tlev(1)%data(:,1,1),            !PINP1
     &  p%gfcn(1)%tlev(1)%c(1)%coord,        !R 
     &  p%gfcn(R2)%tlev(1)%data(:,1,1),             !R**2
     &  p%gfcn(R3)%tlev(1)%data(:,1,1),             !R**3
     &  nr,dr,-dt,2,nr-1,1)

        call elfpii(
     &  gcoars%gfcn(PIN)%tlev(1)%data(:,1,1),     !PINP12
     &  p%gfcn(PIN)%tlev(1)%data(:,1,1),     !PINP1
     &  p%gfcn(PIN)%tlev(2)%data(:,1,1),     !PIN
     &  p%gfcn(PPN)%tlev(1)%data(:,1,1),     !PPNP1
     &  p%gfcn(1)%tlev(1)%c(1)%coord,  !R 
     &  p%gfcn(R2)%tlev(1)%data(:,1,1),  !R**2
     &  p%gfcn(R3)%tlev(1)%data(:,1,1),  !R**3
     &  nr,dr,-dt,2,nr-1,1)


        jst=4
        jfin=nr-3-logic(mod(nr,2).eq.0)
        jinc=2

        call elfppi(
     &  gcoars%gfcn(PPN)%tlev(1)%data(:,1,1),  !PPNP1
     &  gcoars%gfcn(PPN)%tlev(2)%data(:,1,1),  !PPN
     &  gcoars%gfcn(PPN)%tlev(1)%data(:,1,1),  !PPNP1
     &  gcoars%gfcn(PIN)%tlev(2)%data(:,1,1),  !PIN
     &  gcoars%gfcn(1)%tlev(1)%c(1)%coord,  !R 
     &  gcoars%gfcn(R2)%tlev(1)%data(:,1,1),  !R**2
     &  gcoars%gfcn(R3)%tlev(1)%data(:,1,1),  !R**3
     &  nr,dr,2.d0*dt,jst,jfin,jinc)

        call elfpii(
     &  gcoars%gfcn(PIN)%tlev(1)%data(:,1,1),  !PINP1
     &  gcoars%gfcn(PIN)%tlev(2)%data(:,1,1),  !PIN
     &  gcoars%gfcn(PIN)%tlev(1)%data(:,1,1),  !PINP1
     &  gcoars%gfcn(PPN)%tlev(2)%data(:,1,1),  !PPN
     &  gcoars%gfcn(1)%tlev(1)%c(1)%coord,  !R 
     &  gcoars%gfcn(R2)%tlev(1)%data(:,1,1),  !R**2
     &  gcoars%gfcn(R3)%tlev(1)%data(:,1,1),  !R**3
     &  nr,dr,2.d0*dt,jst,jfin,jinc)

        call swptr1(gcoars)

c       -----------------------------------------------------------
c       TWO FINE STEPS
c       -----------------------------------------------------------

        do fstep = 1,2

         gfine%gfcn(W1)%tlev(1)%data(:,1,1)=
     &   gfine%gfcn(PPN)%tlev(1)%data(:,1,1)     ! W1=PPNP1

         call elfppi(
     &   gfine%gfcn(PPN)%tlev(1)%data(:,1,1),  !PPNP1
     &   gfine%gfcn(PPN)%tlev(2)%data(:,1,1),  !PPN
     &   gfine%gfcn(W1)%tlev(1)%data(:,1,1),   !W1
     &   gfine%gfcn(PIN)%tlev(2)%data(:,1,1),  !PIN
     &   gfine%gfcn(1)%tlev(1)%c(1)%coord,     !R 
     &   gfine%gfcn(R2)%tlev(1)%data(:,1,1),   !R**2
     &   gfine%gfcn(R3)%tlev(1)%data(:,1,1),   !R**3
     &   nr,dr,dt,fstep+1,nr-fstep,1)

         call adddis(
     &   gfine%gfcn(PPN)%tlev(1)%data(:,1,1),  !PPNP1
     &   gfine%gfcn(W1)%tlev(1)%data(:,1,1), !W1
     &   epsdis,fstep+4,nr-1)

         gfine%gfcn(W1)%tlev(1)%data(:,1,1)=
     &   gfine%gfcn(PIN)%tlev(1)%data(:,1,1)     ! W1=PINP1
         
         call elfpii(
     &   gfine%gfcn(PIN)%tlev(1)%data(:,1,1),  !PINP1
     &   gfine%gfcn(PIN)%tlev(2)%data(:,1,1),  !PIN
     &   gfine%gfcn(W1)%tlev(1)%data(:,1,1),   !W1
     &   gfine%gfcn(PPN)%tlev(2)%data(:,1,1),  !PPN
     &   gfine%gfcn(1)%tlev(1)%c(1)%coord,  !R 
     &   gfine%gfcn(R2)%tlev(1)%data(:,1,1),  !R**2
     &   gfine%gfcn(R3)%tlev(1)%data(:,1,1),  !R**3
     &   nr,dble(dr),dble(dt),fstep+1,nr-fstep,1)

         call adddis(
     &   gfine%gfcn(PIN)%tlev(1)%data(:,1,1),  !PINP1
     &   gfine%gfcn(W1)%tlev(1)%data(:,1,1), !W1
     &   epsdis,fstep+4,nr-1)

         call swptr1(gfine)

        enddo

c       -------------------------------------------------------
c       COMPUTE TRUNCATION ERROR ESTIMATE
c       -------------------------------------------------------
c        (now PPN, ...  correspond to coarse grid
c             PPN2, ... correspond to fine   grid)

        call hctr2(
     &  gcoars%gfcn(W1)%tlev(1)%data(:,1,1), !W1           (out)
     &  gcoars%gfcn(W2)%tlev(1)%data(:,1,1), !W2           (out)
     &  gcoars%gfcn(W3)%tlev(1)%data(:,1,1), !W3           (out)
     &  gcoars%gfcn(PPN)%tlev(2)%data(:,1,1),  !PPN           (in)
     &  gfine%gfcn(PPN)%tlev(2)%data(:,1,1),   !PPN2          (in)
     &  gcoars%gfcn(PIN)%tlev(2)%data(:,1,1),  !PIN           (in)
     &  gfine%gfcn(PIN)%tlev(2)%data(:,1,1),   !PIN2          (in)
     &  nr)

c       -------------------------------------------------------
c       RELOAD POINTERS
c       -------------------------------------------------------
c        (now PPN, ...  correspond to coarse grid
c             PPN2, ... correspond to original one)

         if(iflag.eq.1) then 
         gcoars%gfcn(W2)%tlev(1)%data(:,1,1)=     !W2
     &   dsqrt(
     &   (p%gfcn(PPN)%tlev(2)%data(:,1,1))**2 +   !PPN          
     &   (p%gfcn(PIN)%tlev(2)%data(:,1,1))**2     !PIN           
     &        )
         idnrm=dsqrt((1.0d0/nr)*DOT_PRODUCT(   
     &   gcoars%gfcn(W2)%tlev(1)%data(:,1,1),     !W2
     &   gcoars%gfcn(W2)%tlev(1)%data(:,1,1)      !W2
     &   ))      
         endif

c       ------------------------------------------------------- 
c       COMPUTE 'MEANED' RMS VALUE OF PP, PI FOR ESTIMATING ...
c       -------------------------------------------------------

         gcoars%gfcn(W2)%tlev(1)%data(:,1,1)=     !W2
     &   dsqrt(
     &   (p%gfcn(PPN)%tlev(2)%data(:,1,1))**2+    !PPN          
     &   (p%gfcn(PIN)%tlev(2)%data(:,1,1))**2     !PIN           
     &        )
         call dvglmf( 
     &   gcoars%gfcn(W2)%tlev(1)%data(:,1,1),     !W2
     &   gcoars%gfcn(W3)%tlev(1)%data(:,1,1),     !W3
     &   nr,2,idnrm)    
       
         write(6,*) gcoars%gfcn(W1)%tlev(1)%data(:,1,1)
         write(6,*) 'end.'

         p%gfcn(ESTRE)%tlev(1)%data(:,1,1)=  !ESTRE=W1/W3
     &   gcoars%gfcn(W1)%tlev(1)%data(:,1,1)/     !W1
     &   gcoars%gfcn(W3)%tlev(1)%data(:,1,1)      !W3

c       -------------------------------------------------------
c       FREE TEMPORARY STORAGE  
c       -------------------------------------------------------

        call deallocateGF3TL(gcoars)
        call deallocateGF3TL(gfine)
        deallocate(gcoars)
        deallocate(gfine)
        return

        CONTAINS 

        subroutine hctr2(trpp1,trpp,trpi,ppc,ppf,pic,pif,nr)
        integer nr
        real, dimension(:) :: trpp1,trpp,trpi,ppc,ppf,pic,pif

        integer logic, firstj,j,lastj

        trpp  = 0.d0
        trpi  = 0.d0
        trpp1 = 0.d0
        
cc        write(6,*) 'ppc',ppc(93:101)
cc        write(6,*) 'ppf',ppf(93:101)
cc        write(6,*) 'pic',pic(93:101)
cc        write(6,*) 'pif',pif(93:101)

        firstj=4
        lastj=nr-(3+logic(mod(nr,2).eq.0))

        trpp(firstj:lastj:2)=ppc(firstj:lastj:2)-ppf(firstj:lastj:2)
        trpi(firstj:lastj:2)=pic(firstj:lastj:2)-pif(firstj:lastj:2)

        firstj = firstj + 1  
        lastj  = lastj  - 1
     
        if(firstj.le.lastj) then 
        FORALL(j=firstj:lastj:2) trpp(j) = 0.5d0*(trpp(j-1)+trpp(j+1))
        FORALL(j=firstj:lastj:2) trpi(j) = 0.5d0*(trpi(j-1)+trpi(j+1))
        endif

        trpp1=sqrt(trpp**2+trpi**2)

cc        write(6,*) 'trpp1',trpp1(4:15) 

        return
        end subroutine hctr2

c-----------------------------------------------------------------------
c
c     Vector generalized mean with floor.
c
c     V2(j) := max(FLOOR,mean(V1(j-P),...V1(j),...V1(j+P)))
c
c-----------------------------------------------------------------------
c
      subroutine dvglmf(v1,v2,n,p,floor)
c
         integer        n
         real         v1(:),        v2(:)
         real         floor
         integer        p
c
         real         dvmax
         real         itwop1
         integer        i,            j
c
         itwop1 = 1.0d0 / (2 * p + 1)
         if( n .ge. (2 * p + 1) ) then
            v2=0.0d0
            do 20 i = p + 1 , n - p
               do 10 j = -p , p
                  v2(i) = v2(i) + v1(i+j)
 10            continue
 20         continue
            do 30 i = p + 1 , n - p
               v2(i) = max(floor,v2(i) * itwop1)
 30         continue
            do 40 i = 1 , p
               v2(i) = v2(p+1)
 40         continue
            do 50 i = n - p + 1 , n
               v2(i) = v2(n-p)
 50         continue
         else
            v2=max(floor,MAXVAL(v1))
         end if
c
         return
c
      end subroutine dvglmf




        subroutine swptr1(p)
        use adaptive
        type(GF3TL), pointer :: p 
        type(GF3), pointer :: p1

        allocate(p1,stat=ierr)
        if(ierr>0) then
          write(6,*) 'swptr1: not enough memory'
          stop
        endif 

!       the first 3 grid functions PPN, PIN and PHN have 2 time levels
        

c        if(nfcnpr.gt.0) then
         do i=1,p%n_gfcn 
          if(p%gfcn(i)%n_t.eq.2) then
           p1=p%gfcn(i)%tlev(1)
           p%gfcn(i)%tlev(1)=p%gfcn(i)%tlev(2)
           p%gfcn(i)%tlev(2)=p1
         endif
        enddo

        return 
        end subroutine swptr1
   
        end 










