subroutine flag_candidates(lp,curlev,finlev)
!------------------------------------------------------------------------
! For all meshes, in a loop over leves from finlev down to curlev
! computes truncation error estimate (subroutine ctrest)
! and then generates a new cluster structure which will serve as a
! template to build a new mesh hierarchy (rgrid2.f90).
!
! The truncation error estimate is stored as a grid function ESTRE. The
! output from subroutine ctrest is then converted to a "binary form",
! i.e., ESTRE is set to 1.0 for all points with truncation error
! estimate larger than a theshold (tremax defined in module rc_param)
! otherwise is set to 0. In this form it it passed to subroutine cluster
! that identifies "rectangular" (1d, 2d or 3d) regions (clusters).
!
! NOTE: there is a predefined maximum number of clusters that can be
! returned by subroutine cluster (parameter maxclust in module rc_param)
!
! User/or application dependent code is buried inside subroutine ctrest.
! clustering routine is taken from PAMRSL, developed by N. Chrisochoides
!
! AUTHOR: Tom Haupt, NPAC at Syracuse University, haupt@npac.syr.edu
!------------------------------------------------------------------------
use types
use amr
type(LEVEL_INFO), dimension(maxlev) :: lp
integer(I4B) curlev,finlev

!local variables
type(ADAPTIVE_MESH), pointer     :: p,parent
type(CLUSTERS),pointer           :: cp,clp,clpnew
real(SP),dimension(3)            :: cmin,cmax

!variables needed for clustering routine        
integer(I4B), dimension(1) :: npatches
integer(I4B) patchno,startpt
integer(I4B), dimension(:), pointer :: patches
integer(I4B), dimension(1,1,1)      :: n1,n2,n3,dim,level
 
INTERFACE
       subroutine ctrest(lp,lev,p)
       use types
       use amr
       type(LEVEL_INFO), dimension(maxlev) :: lp
       integer(I4B) lev
       type(GF3TL), pointer :: p
       end subroutine ctrest 
        
END INTERFACE

write(64,*) '%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%START FLAGGING on level',curlev
  nullify(p)

! estimate truncation error on each grid from the finest level 
! to the current one

LEVELS: do lev=min(finlev,maxlev-1),curlev, -1
write(64,*) 'FLAG lev',lev

    nullify(clp)
    p=>lp(lev)%ml

    DO WHILE(associated(p))
       call ctrest(lp,lev,p%grid)  !truncation error estimate, returns ESTRE

write(71,*) 'ctrest'
write(71,*) p%grid%gfcn(ESTRE)%tlev(1)%data(5,:,1)


       !convert it to a binary form
       WHERE(abs(p%grid%gfcn(ESTRE)%tlev(1)%data).ge.tremax)
        p%grid%gfcn(ESTRE)%tlev(1)%data=1.0d0
       ELSEWHERE
        p%grid%gfcn(ESTRE)%tlev(1)%data=0.0
       ENDWHERE


       ! anything to regrid?
       if(ALL(p%grid%gfcn(ESTRE)%tlev(1)%data<0.5)) then
         write(6,*) 'nothing to regrid!!!!'
         return
       endif 
icont=COUNT(p%grid%gfcn(ESTRE)%tlev(1)%data<0.5)
write(6,*) 'ICONT=',icont
if(icont>389) then      
  write(6,*) 'EMERGENCY RETURN',icont
  return
endif

write(64,*) 'ISOLATE'
!      ------------------------------------------------------------------   
!      ISOLATE CLUSTERS
!      ------------------------------------------------------------------

       ! interface to subroutine cluster written in C
       npatches=0
       n1=p%gnr(1)
       n2=p%gnr(2)
       n3=p%gnr(3)
       if(maxdim<3) n3=1
       level=lev
       allocate(patches(maxclust*6))
       patches=0
       ndx=3 ! for 2d case we use (:,:,1) and 1d (:,1,1) so always ndx=3

write(36,*) n1,n2,n3,level
write(36,*) int(p%grid%gfcn(ESTRE)%tlev(1)%data) 
write(6,*) n1,n2,n3,level

       call cluster(ndx,n1,n2,n3,int(p%grid%gfcn(ESTRE)%tlev(1)%data),&
                    npatches,patches,level)

       if(npatches(1)>maxclust) then
         write(6,*) 'subroutine CLUSTER returned',npatches,' clusters'
         write(6,*) 'i.e., more than allowed maximum of',maxclust
         write(6,*) 'adjust parameter maxclust in module rc_param and', &
                    'run the program again'
         stop
       endif 


!     ---------------------------------------------------------------- 
!     CREATE CLUSTER STRUCTURE AT A FINER LEVEL
!     ----------------------------------------------------------------

      ! unpack data returned by the C routine cluster         
      if(npatches(1)>0) then
      do patchno = 1, npatches(1)
      if(lev.lt.maxlev) then

!        ----------------------------------------------------------------
!        TRANSLATE ARRAY INDICES INTO COORDINATES AND ADD BUFFERS
!        ----------------------------------------------------------------
         startpt=1+(patchno-1)*6
         write(6,*) lev,lp(lev)%lbwx,lp(lev)%lbwy,lp(lev)%lbwz

         i1=patches(startpt)
         i2=i1+patches(startpt+3)+1     ! a bug in clustering routine?
         if(i1==0) i1=1                 ! a bug in clustering routine?
         cmin(1)=p%grid%c(1)%coord(i1) - lp(lev+1)%lbwx
         cmax(1)=p%grid%c(1)%coord(i2) + lp(lev+1)%lbwx
         if(cmin(1) < gxmin) cmin(1)=gxmin
         if(cmax(1) > gxmax) cmax(1)=gxmax

         i1=patches(startpt+1)
         i2=i1+patches(startpt+4)
         cmin(2)=p%grid%c(2)%coord(i1) - lp(lev+1)%lbwy
         cmax(2)=p%grid%c(2)%coord(i2) + lp(lev+1)%lbwy
         if(cmin(2) < gymin) cmin(2)=gymin
         if(cmax(2) > gymax) cmax(2)=gymax

write(71,*) 'cluster'
write(71,*) i1,i2,cmin(2),cmax(2)

         i1=patches(startpt+2)
         i2=i1+patches(startpt+5)-1      ! a bug?
         cmin(3)=p%grid%c(3)%coord(i1) - lp(lev+1)%lbwz
         cmax(3)=p%grid%c(3)%coord(i2) + lp(lev+1)%lbwz
         if(cmin(3) < gzmin) cmin(3)=gzmin
         if(cmax(3) > gzmax) cmax(3)=gzmax

        ! eliminate small clusters
   if((cmax(1)-cmin(1)+lp(lev)%ldr(1))/lp(lev)%ldr(1)<mcsize) cycle     
   if(maxdim>1 .and. (cmax(2)-cmin(2)+lp(lev)%ldr(2))/lp(lev)%ldr(2)<mcsize) cycle 
   if(maxdim>2 .and. (cmax(3)-cmin(3)+lp(lev)%ldr(3))/lp(lev)%ldr(3)<mcsize) cycle

!     ----------------------------------------------------------------- 
!     ADD BUFFERS  
!     -----------------------------------------------------------------

        !cluster accepted: create corresponding structure 
        allocate(clpnew)
        nullify(clpnew%clnext)
        nullify(clpnew%clprnt)
        clpnew%clrmin=cmin   
        clpnew%clrmax=cmax   
        clpnew%clprnt=>p

        if(associated(clp)) then
         clp%clnext=>clpnew       !next cluster at this level
        else
         lp(lev+1)%clhead=>clpnew !first cluster at this level

        endif

        if(maxdim<3) then
         clpnew%clrmin(3)=1
         clpnew%clrmax(3)=1
        endif

       if(maxdim<2) then
        clpnew%clrmin(2)=1
        clpnew%clrmax(2)=1
       endif

       clp=>clpnew
      else
        write(6,*) 'isolate clusters: maximum levels exceeded'
        write(6,*) 'this should never happen!'
        stop
      endif
      enddo !patchno          
      endif !npatches>0

      deallocate(patches) 
      p=>p%neighbor
      ENDDO !process next grid at this level (while associated(p))
write(64,*) '======= end level',lev
enddo  LEVELS

return
end




