       subroutine regrid(lp,curlev,finlev,first)
       use amr
       type(LEVEL_INFO), dimension(maxlev) :: lp
       integer curlev,finlev
       logical first
       
       interface

       subroutine flag_candidates(lp,curlev,finlev)
       use amr
       type(LEVEL_INFO), dimension(maxlev) :: lp
       integer curlev,finlev   
       end subroutine flag_candidates

       subroutine rgrid2(lp,curlev,finlev)
       use amr
       type(LEVEL_INFO), dimension(maxlev) :: lp
       integer curlev,finlev
       end subroutine rgrid2

       end interface

       if(curlev.eq.maxlev) then
          write(6,*) 'REGRID: should not be calling routine',
     &    ' at finest allowable level'
          write(6,*) curlev,maxlev
          stop
          return
       endif

       write(6,*) ' -------- START REGRIDDING ------------' 

       call clear_cluster_structure(lp)
       call flag_candidates(lp,curlev,finlev)

       call rgrid2(lp,curlev,finlev)

       return
       end


       subroutine flag_candidates(lp,curlev,finlev)
       use amr
       type(LEVEL_INFO), dimension(maxlev) :: lp
       integer curlev,finlev

       type(ADAPTIVE_MESH), pointer :: p
       type(CLUSTER),pointer        :: cp
       real*8 dx

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

       subroutine merge_clusters(lp,lev)
       use amr
       type(LEVEL_INFO),dimension(maxlev) :: lp
       integer lev
       end subroutine merge_clusters

       END INTERFACE

       nullify(p)

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

       do   lev=min(finlev,maxlev-1),curlev, -1
            write(6,*) 'flag candidates on level=',lev
            p=>lp(lev)%ml
            DO WHILE(associated(p))
             call ctrest(lp,lev,p%grid)

             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

c          ----------------------------------------------------------
c          ZERO ELEMENTS OF ESTRE CORRESPONDING TO SPATIAL LOCATIONS
c          WHICH ARE SPECIFFICALLY EXCLUDED FROM REFINEMENT AT THIS
c          LEVEL
c          ----------------------------------------------------------
c            rcp/rcpi structures; ignored at this time    
c          ----------------------------------------------------------
c          ZERO ELEMENTS OF ESTRE IN BUFFER REGIONS
c          ----------------------------------------------------------

c          this works for 1-d system only

           if(lev.ge.crslev+1) then
             o=lp(lev)%lbwx
             if(p%gflb(1).eq.flbin) then
              p%grid%gfcn(ESTRE)%tlev(1)%data(1:o+1,1,1)=0.0
             endif
             if(p%gfrb(1).eq.flbin) then
              p%grid%gfcn(ESTRE)%tlev(1)%data(nr-o-1:nr,1,1)=0.0
             endif
           endif  
c          ----------------------------------------------------------         
c          FORCE PROPER NESTING
c          ---------------------------------------------------------- 

c           write(6,*) 'force nesting on level',lev
           IF(lev.le. maxlev-2) THEN 
             if(associated(lp(lev+2)%clhead)) then
c               write(6,*) lev+2,' exists'
               cp=>lp(lev+2)%clhead
               DO WHILE(associated(cp))
                  dx=lp(lev)%ldr(1)
                  nx=p%gnr(1)

c             write(6,*) 'level',lev+2, cp%clrmin(1),cp%clrmax(1)
c             ixx=1
c             do while(p%grid%gfcn(ESTRE)%tlev(1)%data(ixx,1,1)<0.5d0)
c             ixx=ixx+1
c             enddo
c             write(6,*) 'ixx=',ixx
c         write(6,*) 'rmin=',p%grid%gfcn(ESTRE)%tlev(1)%c(1)%coord(ixx)
c         write(6,*) p%grid%gfcn(ESTRE)%tlev(1)%data(:,1,1)
c             iyy=ixx+2 
c             do while(p%grid%gfcn(ESTRE)%tlev(1)%data(iyy,1,1)>0.5d0)
c             iyy=iyy+1
c             enddo
c             write(6,*) 'iyy=',iyy
              
                  if(cp%clrmax(1).gt.p%grmin(1).and.
     &               cp%clrmin(1).lt.p%grmax(1)) then 
           i1=min(nx,int(max((cp%clrmin(1)-p%grmin(1))/dx,0d0))+1)
           i2=min(nx,int(max((cp%clrmax(1)-p%grmin(1))/dx,0d0))+2)

           p%grid%gfcn(ESTRE)%tlev(1)%data(i1:i2+1,1,1)=1.d0

                  endif 
        
                 cp=>cp%clnext
               END DO
             endif
           ENDIF

           p=>p%neighbor
           ENDDO

c       write(6,*) '---isolate---'
       call isolate_clusters(lp,lev)
       call merge_clusters(lp,lev)
       write(6,*) '-done with level',lev,'-----------------'
       enddo

       return
       end

       subroutine isolate_clusters(lp,lev)
       use amr
       type(LEVEL_INFO),dimension(maxlev) :: lp
       integer lev

       type(CLUSTER),pointer              :: clp,clpnew
       
       type(ADAPTIVE_MESH),pointer        :: p
       type(CLUSTER),pointer              :: cp
       real*8 rmin,rmax
       integer npt,start 

       INTERFACE
       subroutine rgexcln(p,lp,rmin,rmax,start,npt)
       use amr
       type(ADAPTIVE_MESH), pointer           :: p
       type(LEVEL_INFO), dimension(maxlev)    :: lp 
       real*8                                 :: rmin,rmax
       integer                                :: start,npt
       end subroutine rgexcln

       END INTERFACE

       if(.not.associated(lp(lev)%ml)) return
       p=>lp(lev)%ml
       nullify(clp) 

       DO WHILE(associated(p))
         start=1

 100     continue
         call rgexcln(p,lp,rmin,rmax,start,npt)
         if(start.gt.0) then
C???           if(.not.rcebl(lp,p,rmin,rmax,1)) goto 100
           IF(lev.lt.maxlev) then
                 write(6,*) 'clusters on level',lev+1
                 write(6,*) 'new cluster:',rmin,rmax
                 allocate(clpnew)
                 nullify(clpnew%clnext)
                 nullify(clpnew%clprnt)
                 clpnew%clrmin=0
                 clpnew%clrmax=0
                 clpnew%clrmin(1)=rmin
                 clpnew%clrmax(1)=rmax
                 clpnew%clprnt=>p
                 if(associated(clp)) then
                     clp%clnext=>clpnew
                 else
                     lp(lev+1)%clhead=>clpnew
                 endif
                 clp=>clpnew
           ELSE
              write(6,*) 'isolate clusters: maximum levels exceeded'
           ENDIF   
         endif
         if(start.gt.0) goto 100 
         p=>p%neighbor
        
       END DO             
  
       return
       end

       subroutine merge_clusters(lp,lev)
c      --------------------------------------------------------
c      ADD BUFFER REGIONS AND MERGE
c      --------------------------------------------------------
       use amr
       type(LEVEL_INFO),dimension(maxlev) :: lp
       integer lev
  
       type(CLUSTER),pointer :: clp,clpnext

       if(associated(lp(lev+1)%clhead)) then
         clp=>lp(lev+1)%clhead

       DO WHILE(associated(clp))

c         write(6,*) 'before:',clp%clrmin(1)

         clp%clrmin(1)=max(clp%clrmin(1) - lp(lev)%ldr(1)*lbwx,
     &    clp%clprnt%grmin(1), gxmin)

c         write(6,*) 'after:',clp%clrmin(1)

         if(clp%clrmin(1)-glrmin.lt.rminxtnd*lp(lev)%ldr(1))
     &    clp%clrmin(1)=gxmin

         clp%clrmax(1)=min(
     &    clp%clrmax(1) + lp(lev)%ldr(1)*lbwx,
     &    clp%clprnt%grmax(1),
     &    gxmax)

c         write(6,*) "AA",gxmax
c         write(6,*) "AB",clp%clrmax(1)
c         write(6,*) "AC",rmaxxtnd
c         write(6,*) "AD",lp(lev)%ldr(1)

         if(gxmax-clp%clrmax(1).lt.rmaxxtnd*lp(lev)%ldr(1))
     &    clp%clrmax(1)=gxmax

         clp=>clp%clnext 
       ENDDO

         if(.not.associated(lp(lev+1)%clhead)) goto 345
         clp=>lp(lev+1)%clhead
         if(.not.associated(clp%clnext)) goto 345

 340     continue       

         clpnext=>clp%clnext
 
         IF(clpnext%clrmin(1)-clp%clrmax(1).le.
     &     lp(lev)%ldr(1)*mrgwid) THEN
           if(clpnext%clprnt%g.eq.clp%clprnt%g) then
              clp%clrmax(1)=clpnext%clrmax(1)
              clp%clnext => clpnext%clnext
              deallocate(clpnext)
              if(associated(clp%clnext)) goto 340
           endif
         ENDIF 
         if(associated(clpnext%clnext)) then
           clp=>clpnext   
           goto 340
         endif
 345     continue
       endif  

       return
       end

       subroutine rgrid2(lp,curlev,finlev)
       use amr
       type(LEVEL_INFO), dimension(maxlev) :: lp
       integer curlev,finlev

       type(ADAPTIVE_MESH),pointer :: parent,child,head,current,old
       type(CLUSTER),pointer       :: clp
       character(200)              :: chname = 'child'  
       integer g 

c      ------------------------------------------------------- 
c      DETERMINE NEW FINEST LEVEL OF REFINEMENT 
c      -------------------------------------------------------

       finlev=curlev
       do lev = curlev + 1, maxlev
          if(associated(lp(lev)%clhead)) finlev = lev
       enddo

c       zrob kopie biezacej struktury 

       IF(curlev+1.le.finlev) THEN       
       do lev=curlev+1,finlev 

       nullify(child)
       nullify(head)
       nullify(parent)
       nullify(current)

       clp=>lp(lev)%clhead
 
       DO WHILE(associated(clp))

          call find_parent(lp(lev-1)%ml,parent,clp%clrmin,clp%clrmax)
          if(.not.associated(parent)) then
             write(6,*) 'regrid: parent not found'
             stop
          endif

           call initialize_grid(child,lev,chname,clp%clrmin(1),
     &     clp%clrmin(2),clp%clrmin(3),clp%clrmax(1),clp%clrmax(2),
     &     clp%clrmax(3),lp(lev)%ldr)

         
          child%parent=>parent 
c          nullify(child%neighbor)
c          nullify(child%sibiling)
c          nullify(child%child)

          is=int(lp(lev)%rcxshp)
          call initialize_data(child,parent,lp(lev)%ml,lp(lev)%ldr,
     &    lp(lev-1)%ldr,lp(lev)%ldt,is)

          child%parent=>parent
          if(associated(current)) then
             if(parent%g == current%parent%g) then
                current%sibiling=>child
             endif
          endif

          if(associated(head)) then             
             current%neighbor=>child
             current=>current%neighbor           
             nullify(current%neighbor)
          else
             head=>child
             current=>child
          endif 

          clp=>clp%clnext

         ENDDO

        current=>lp(lev)%ml
        DO WHILE(associated(current))
         old=>current%neighbor
         call deallocateGF3TL(current%grid)
         g=get_g(current%g)
         deallocate(current)
         current=>old
        ENDDO

        lp(lev)%ml=>head
        ENDDO 

       ENDIF

c       write(6,*) ' existing levels after regridding'
c       do i=1,maxlev
c       if(associated(lp(i)%ml)) write(6,*) 'ml level',i
c       enddo

       return 
       end subroutine rgrid2



       subroutine rgexcln(p,lp,rmin,rmax,start,npt)
       use amr
       type(ADAPTIVE_MESH), pointer           :: p
       type(LEVEL_INFO), dimension(maxlev)    :: lp 
       real*8                                 :: rmin,rmax
       integer                                :: start,npt

C      -------------------------------------------------------------
C      BEGINNING AT START-TH POINT ON GRID, SCANS FOR CLUSTERS OF
C      POINTS. IF ONE IS FOUND, RMIN AND RMAX ARE RETURNED ALONG 
C      WITH THE UPDATED VALUE OF START (-1 IF NO CLUSTER FOUND),
C      AND NUMBER OF POINTS, NPT, IN CLUSTER
C      -------------------------------------------------------------       
       nr=p%gnr(1)
       npt=0
       level=p%level
 100   continue
c       write(6,*) 'rgexcln:',nr,start
         if(start.gt.nr) then      
            start=-1
            return
         endif

         if(p%grid%gfcn(ESTRE)%tlev(1)%data(start,1,1).eq.1.0d0) then
            rmin=p%grmin(1)+lp(level)%ldr(1)*(start-1)
 200        continue
c       write(6,*) 'rgexcln: 200'
            
            start=start+1
            if(start.gt.nr) then
               rmax=p%grmax(1)
               npt=int(rmax-rmin)/lp(level)%ldr(1) + 1
               return
             endif
         if(p%grid%gfcn(ESTRE)%tlev(1)%data(start,1,1).eq.0.0d0) then
            rmax=p%grmin(1)+lp(level)%ldr(1)*(start-1)
            npt=int((rmax-rmin)/lp(level)%ldr(1)) + 1
            return             
         endif
         goto 200
         endif
         start=start+1
         goto 100

      end subroutine rgexcln








































