c-----------------------------------------------------------
      subroutine adaptation (ifmortar,step)
c-----------------------------------------------------------
c     For 3-D mesh adaptation (refinement+ coarsening)
c-----------------------------------------------------------
      include 'header.h'
      
      logical if_coarsen,if_refine,ifmortar,ifrepeat
      integer iel,miel,irefine,icoarsen,neltold,step

      if (timeron) call timer_start(t_adaptation)
      ifmortar=.false.
c.....compute heat source center(x0,y0,z0)
      x0=x00+velx*time
      y0=y00+vely*time
      z0=z00+velz*time

c.....Search elements to be refined. Check with restrictions. Perform
c     refinement repeatedly until all desired refinements are done.

c.....ich(iel)=0 no grid change on element iel
c.....ich(iel)=2 iel is marked to be coarsened
c.....ich(iel)=4 iel is marked to be refined

c.....irefine records how many elements got refined
      irefine=0

c.....check whether elements need to be refined because they have overlap
c     with the  heat source
4     call find_refine(if_refine)

      if(if_refine) then
        ifrepeat=.true.
2       if(ifrepeat) then
c.........Check with restriction, unmark elements that cannot be refined.
c         Elements preventing desired refinement will be marked to be refined.
          call check_refine(ifrepeat) 
          go to 2
        end if
c.......perform refinement
        call do_refine(ifmortar,irefine)
        goto 4
      endif

c.....Search for elements to be coarsened. Check with restrictions,
c     Perform coarsening repeatedly until all possible coarsening
c     is done.

c.....icoarsen records how many elements got coarsened 
      icoarsen=0

c.....skip(iel)=.true. indicates an element no longer exists (because it
c     got merged)
c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(iel)
       do iel=1,nelt
        skip(iel)=.false.
      end do
c$OMP END PARALLEL DO

      neltold=nelt

c.....Check whether elements need to be coarsened because they don't have
c     overlap with the heat source. Only elements that don't have a larger 
c     size neighbor can be marked to be coarsened

5     call find_coarsen(if_coarsen,neltold)

      if(if_coarsen) then
c.......Perform coarsening, however subject to restriction. Only possible 
c       coarsening will be performed. if_coarsen=.true. indicates that
c       actual coarsening happened
        call do_coarsen(if_coarsen,icoarsen,neltold)
        if(if_coarsen) then
c.........ifmortar=.true. indicates the grid changed, i.e. the mortar points 
c         indices need to be regenerated on the new grid.
          ifmortar=.true.
          go to 5
        end if 
      end if

      write(*,1000) step, irefine, icoarsen, nelt
 1000 format('Step ',i4, ': elements refined, merged, total:',
     &       i6, 1X , i6, 1X, i6)

c.....mt_to_id(miel) takes as argument the morton index  and returns the actual 
c                    element index
c.....id_to_mt(iel)  takes as argument the actual element index and returns the 
c                    morton index
c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(miel,iel)
      do miel=1,nelt
        iel=mt_to_id(miel)
        id_to_mt(iel)=miel
      end do 
c$OMP END PARALLEL DO

c.....Reorder the elements in the order of the morton curve. After the move 
c     subroutine the element indices are  the same as the morton indices
      call move

c.....if the grid changed, regenerate mortar indices and update variables
c     associated to grid.
      if (ifmortar) then
        call mortar
        call prepwork
      endif
      if (timeron) call timer_stop(t_adaptation)

      return
      end 


c-----------------------------------------------------------
      subroutine do_coarsen(if_coarsen,icoarsen,neltold)
c---------------------------------------------------------------
c     Coarsening procedure: 
c     1) check with restrictions
c     2) perform coarsening
c---------------------------------------------------------------

      include 'header.h'

      logical if_coarsen, icheck,test,test1,test2,test3
      integer iel, ntp(8), ntempmin, ic, parent, mielnew, miel,
     &        icoarsen, i, index, num_coarsen, ntemp, ii, ntemp1, 
     &        neltold
      
      if_coarsen=.false.

c.....If an element has been merged, it will be skipped afterwards
c     skip(iel)=.true. for elements that will be skipped.
c     ifcoa_id(iel)=.true. indicates that element iel will be coarsened
c     ifcoa(miel)=.true. refers to element miel(mortar index) will be
c                        coarsened

c$OMP PARALLEL DEFAULT(SHARED) PRIVATE(iel)
c$OMP DO 
      do iel=1,nelt
        mt_to_id_old(iel)=mt_to_id(iel)
        mt_to_id(iel)=0
      end do
c$OMP END DO nowait
c$OMP DO 
      do iel=1,neltold 
        ifcoa_id(iel)=.false.
      end do
c$OMP END DO nowait
c$OMP END PARALLEL

c.....Check whether the potential coarsening will make neighbor, 
c     and neighbor's neighbor....break grid restriction

c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(miel,iel,ic,
c$OMP& ntp,parent,test,test1,i,test2,test3)
c$OMP& SHARED(if_coarsen)
      do miel=1,nelt
        ifcoa(miel)=.false.
        front(miel)=0
        iel=mt_to_id_old(miel)
c.......if an element is marked to be coarsened
        if(ich(iel).eq.2) then

c.........If the current  element is the "first" child (front-left-
c         bottom) of its parent (tree(iel) mod 8 equals 0), then 
c         find all its neighbors. Check whether they are from the same 
c         parent.

          ic=tree(iel)
          if(.not.btest(ic,0).and..not.btest(ic,1).and.
     &       .not.btest(ic,2)) then
            ntp(1)=iel
            ntp(2)=sje(1,1,1,iel)
            ntp(3)=sje(1,1,3,iel)
            ntp(4)=sje(1,1,1,ntp(3))
            ntp(5)=sje(1,1,5,iel)
            ntp(6)=sje(1,1,1,ntp(5))
            ntp(7)=sje(1,1,3,ntp(5))
            ntp(8)=sje(1,1,1,ntp(7))
 
            parent=ishft(tree(iel),-3)
            test=.false.

            test1=.true.
            do i=1,8
              if(ishft(tree(ntp(i)),-3).ne.parent)test1=.false.
            end do

c...........check whether all child elements are marked to be coarsened
            if(test1)then
              test2=.true.
              do i=1,8
                if(ich(ntp(i)).ne.2)test2=.false.
              end do

c.............check whether all child elements can be coarsened or not.
              if(test2)then
                test3=.true.
                do i=1,8
                  if(.not.icheck(ntp(i),i))test3=.false.
                end do
                if(test3)test=.true.
              end if
            end if
c...........if the eight child elements are eligible to be coarsened
c           mark the first children ifcoa(miel)=.true.
c           mark them all ifcoa_id()=.true.
c           front(miel) will be used to calculate (potentially in parallel) 
c                       how many elements with seuqnece numbers less than
c                       miel will be coarsened.
c           skip()      marks that an element will no longer exist after merge.

            if(test)then

              ifcoa(miel)=.true.
              do i=1,8
                ifcoa_id(ntp(i))=.true.
              end do
              front(miel)=1
              do i=1,7
                 skip(ntp(i+1))=.true.
              end do
              if(.not.if_coarsen) if_coarsen=.true.
            end if
          end if 
        end if 
      end do 
c$OMP END PARALLEL DO

c.....compute front(iel), how many elements will be coarsened before iel
c     (including iel)
      call parallel_add(front)

c.....num_coarsen is the total number of elements that will be coarsened
      num_coarsen=front(nelt)

c.....action(i) records the morton index of the i'th element (if it is an
c     element's front-left-bottom-child) to be coarsened.

c.....create array mt_to_id to convert actual element index to morton index
c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(miel,iel,mielnew)
      do miel=1,nelt
        iel=mt_to_id_old(miel)
        if(.not.skip(iel))then
          if(ifcoa(miel))then
            action(front(miel))=miel
            mielnew=miel-(front(miel)-1)*7
          else 
            mielnew=miel-front(miel)*7
          end if
          mt_to_id(mielnew)=iel
        end if
      end do
c$OMP END PARALLEL DO

c.....perform the coarsening procedure (potentially in parallel)
c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(index,miel,iel,ntp)
      do index=1,num_coarsen
        miel=action(index)
        iel=mt_to_id_old(miel)
c.......find eight child elements to be coarsened
        ntp(1)=iel
        ntp(2)=sje(1,1,1,iel)
        ntp(3)=sje(1,1,3,iel)
        ntp(4)=sje(1,1,1,ntp(3))
        ntp(5)=sje(1,1,5,iel)
        ntp(6)=sje(1,1,1,ntp(5))
        ntp(7)=sje(1,1,3,ntp(5))
        ntp(8)=sje(1,1,1,ntp(7))
c.......merge them to be the parent
        call merging(ntp)
      end do
c$OMP END PARALLEL DO
      nelt=nelt-num_coarsen*7
      icoarsen=icoarsen+num_coarsen*8

      return
      end

c-------------------------------------------------------
      subroutine do_refine(ifmortar,irefine)
c-------------------------------------------------------
c     Refinement procedure
c--------------------------------------------------------

      include 'header.h'

      logical ifmortar
      double precision xctemp(8), yctemp(8), zctemp(8), xleft, xright,
     &       yleft, yright, zleft, zright, ta1temp(lx1,lx1,lx1),
     &       xhalf, yhalf, zhalf
      integer iel, i, ii, jj, j, jface, 
     &        ntemp, ndir, facedir, k, le(4), ne(4), mielnew,
     &        miel, irefine,ntemp1, num_refine, index, treetemp,
     &        sjetemp(2,2,6), n1, n2, nelttemp,
     &        cb, cbctemp(6)

c.....initialize

c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(miel)
      do miel=1,nelt
        mt_to_id_old(miel)=mt_to_id(miel)
        mt_to_id(miel)=0
        action(miel)=0
        if(ich(mt_to_id_old(miel)).ne.4)then
          front(miel)=0
        else
          front(miel)=1
        end if
      end do
c$OMP END PARALLEL DO

c.....front(iel) records how many elements with sequence numbers less than
c     or equal to iel will be refined
      call parallel_add(front)

c.....num_refine is the total number of elements that will be refined
      num_refine=front(nelt)

c.....action(i) records the morton index of the  i'th element to be refined
c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(miel,iel)
      do miel=1,nelt
        iel=mt_to_id_old(miel)
        if(ich(iel).eq.4)then
          action(front(miel))=miel
        end if
      end do
c$OMP END PARALLEL DO

c.....Compute array mt_to_id to convert the element index to morton index.
c     ref_front_id(iel) records how many elements with index less than
c     iel (actual element index, not morton index), will be refined.
c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(miel,iel,ntemp,mielnew)
      do miel=1,nelt
        iel=mt_to_id_old(miel)
        if(ich(iel).eq.4)then
          ntemp=(front(miel)-1)*7
          mielnew=miel+ntemp
        else
          ntemp=front(miel)*7
          mielnew=miel+ntemp
        end if

        mt_to_id(mielnew)=iel
        ref_front_id(iel)=nelt+ntemp
      end do
c$OMP END PARALLEL DO


c.....Perform refinement (potentially in parallel): 
c       - Cut an element into eight children.
c       - Assign them element index  as iel, nelt+1,...., nelt+7.
c       - Update neighboring information.

      nelttemp=nelt

      if (num_refine .gt. 0) then
        ifmortar=.true.
      endif
c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(index,miel,mielnew,iel,nelt,
c$OMP& treetemp,xctemp,yctemp,zctemp,cbctemp,sjetemp,ta1temp,
c$OMP& ii,jj,ntemp,xleft,xright,xhalf,yleft,yright,yhalf,zleft,zright,
c$OMP& zhalf,ndir,facedir,jface,cb,le,ne,n1,n2,i,j,k)
      do index=1, num_refine  
c.......miel is old morton index and mielnew is new morton index after refinement.
        miel=action(index)
        mielnew=miel+(front(miel)-1)*7
        iel=mt_to_id_old(miel) 
        nelt=nelttemp+(front(miel)-1)*7 
c.......save iel's information in a temporary array
        treetemp=tree(iel)
        do i=1,8
          xctemp(i)=xc(i,iel)
          yctemp(i)=yc(i,iel)
          zctemp(i)=zc(i,iel)
        end do
        do i=1,6
          cbctemp(i)=cbc(i,iel)
          do jj=1,2
            do ii=1,2
              sjetemp(ii,jj,i)=sje(ii,jj,i,iel)
            end do
          end do
        end do
        call copy(ta1temp,ta1(1,1,1,iel),nxyz)


c.......zero out iel here
        tree(iel)=0
        call nr_init(cbc(1,iel),6,0)
        call nr_init(sje(1,1,1,iel),24,0)
        call nr_init(ijel(1,1,iel),12,0)
        call r_init(ta1(1,1,1,iel),nxyz,0.d0)


c.......initialize new child elements:iel and nelt+1~nelt+7
        do j=1,7 
          mt_to_id(mielnew+j)=nelt+j
          tree(nelt+j)=0
          call nr_init(cbc(1,nelt+j),6,0)
          call nr_init(sje(1,1,1,nelt+j),24,0)
          call nr_init(ijel(1,1,nelt+j),12,0)
          call r_init(ta1(1,1,1,nelt+j),nxyz,0.d0)
        end do
          
c.......update the tree()
        ntemp=ishft(treetemp,3)
        tree(iel)=ntemp
        do i=1,7
          tree(nelt+i)=ntemp+mod(i,8)
        end do   
c.......update the children's vertices' coordinates
        xhalf=xctemp(1)+(xctemp(2)-xctemp(1))/2.d0
        xleft=xctemp(1)
        xright=xctemp(2)
        yhalf=yctemp(1)+(yctemp(3)-yctemp(1))/2.d0
        yleft=yctemp(1)
        yright=yctemp(3)
        zhalf=zctemp(1)+(zctemp(5)-zctemp(1))/2.d0
        zleft=zctemp(1)
        zright=zctemp(5)
       
        do j=1,7,2
          do i=1,7,2
            xc(i,nelt+j)     = xhalf
            xc(i+1,nelt+j)   = xright 
          end do
        end do

        do j=2,6,2
          do i=1,7,2
            xc(i,nelt+j)   = xleft
            xc(i+1,nelt+j) = xhalf
          end do
        end do
         
        do i=1,7,2
          xc(i,iel)=xleft
          xc(i+1,iel)=xhalf
        end do

        do i=1,2
          yc(i,nelt+1)=yleft
          yc(i,nelt+4)=yleft
          yc(i,nelt+5)=yleft
          yc(i+4,nelt+1)=yleft
          yc(i+4,nelt+4)=yleft
          yc(i+4,nelt+5)=yleft
        enddo
        do i=3,4
          yc(i,nelt+1)=yhalf
          yc(i,nelt+4)=yhalf
          yc(i,nelt+5)=yhalf
          yc(i+4,nelt+1)=yhalf
          yc(i+4,nelt+4)=yhalf
          yc(i+4,nelt+5)=yhalf
        end do
        do j=2,3
          do i=1,2
            yc(i,nelt+j)=yhalf
            yc(i,nelt+j+4)=yhalf
            yc(i+4,nelt+j)=yhalf
            yc(i+4,nelt+j+4)=yhalf
          end do
          do i=3,4
            yc(i,nelt+j)=yright
            yc(i,nelt+j+4)=yright
            yc(i+4,nelt+j)=yright
            yc(i+4,nelt+j+4)=yright
          end do
        end do
          
        do i=1,2
          yc(i,iel)=yleft
          yc(i+4,iel)=yleft
        end do
        do i=3,4
          yc(i,iel)=yhalf
          yc(i+4,iel)=yhalf
        end do

        do j=1,3
          do i=1,4
            zc(i,nelt+j)=zleft
            zc(i+4,nelt+j)=zhalf
          end do
        end do
        do j=4,7
          do i=1,4
            zc(i,nelt+j)=zhalf
            zc(i+4,nelt+j)=zright
          end do
        end do
        do i=1,4
          zc(i,iel)=zleft
          zc(i+4,iel)=zhalf
        end do

c.......update the children's neighbor information

c.......ndir refers to the x,y,z directions, respectively.
c       facedir refers to the orientation of the face in each direction, 
c       e.g. ndir=1, facedir=0 refers to face 1,
c       and ndir =1, facedir=1 refers to face 2.

        do ndir = 1, 3
          do facedir = 0, 1
            i=2*ndir-1+facedir
            jface=jjface(i)
            cb=cbctemp(i)

c...........find the new element indices of the four children on each
c           face of the parent element
            do k = 1, 4
              le(k) = le_arr(k,facedir,ndir)+nelt
              ne(k) = le_arr(k,1-facedir,ndir)+nelt
            end do
            if(facedir.eq.0)then
              le(1)=iel
            else
              ne(1)=iel
            end if
c...........update neighbor information of the four child elements on each 
c           face of the parent element
            do k=1,4
              cbc(i,le(k))=2
              sje(1,1,i,le(k))=ne(k)
              ijel(1,i,le(k))=1
              ijel(2,i,le(k))=1
            end do

c...........if the face type of the parent element is type 2
            if(cb.eq.2) then
              ntemp=sjetemp(1,1,i)

c.............if the neighbor ntemp is not marked to be refined
              if(ich(ntemp).ne.4)then
                cbc(jface,ntemp)=3
                ijel(1,jface,ntemp)=1
                ijel(2,jface,ntemp)=1
  
                do k=1,4
                  cbc(i,ne(k))=1
                  sje(1,1,i,ne(k))=ntemp
                  if(k.eq.1) then
                    ijel(1,i,ne(k))=1
                    ijel(2,i,ne(k))=1
                    sje(1,1,jface,ntemp)=ne(k)
                  elseif(k.eq.2) then
                    ijel(1,i,ne(k))=1
                    ijel(2,i,ne(k))=2
                    sje(1,2,jface,ntemp)=ne(k)
                  elseif(k.eq.3) then
                    ijel(1,i,ne(k))=2
                    ijel(2,i,ne(k))=1
                    sje(2,1,jface,ntemp)=ne(k)
                  elseif(k.eq.4) then
                    ijel(1,i,ne(k))=2
                    ijel(2,i,ne(k))=2
                    sje(2,2,jface,ntemp)=ne(k)
                  end if
                end do

c.............if the neighbor ntemp is also marked to be refined
              else
                n1=ref_front_id(ntemp)
                 
                do k=1,4
                  cbc(i,ne(k))=2
                  n2=n1+le_arr(k,facedir,ndir)
                  if(n2.eq.n1+8)n2=ntemp
                  sje(1,1,i,ne(k))=n2
                  ijel(1,i,ne(k))=1
                end do

              endif
c...........if the face type of the parent element is type 3
            elseif(cb.eq.3) then
              do k=1,4
                cbc(i,ne(k))=2
                if(k.eq.1) then
                  ntemp=sjetemp(1,1,i)
                elseif (k.eq.2) then
                  ntemp=sjetemp(1,2,i)
                elseif(k.eq.3) then
                  ntemp=sjetemp(2,1,i)
                elseif(k.eq.4) then
                  ntemp=sjetemp(2,2,i)
                end if
                ijel(1,i,ne(k))=1
                ijel(2,i,ne(k))=1
                sje(1,1,i,ne(k))=ntemp
                cbc(jface,ntemp)=2
                sje(1,1,jface,ntemp)=ne(k)
                ijel(1,jface,ntemp)=1
                ijel(2,jface,ntemp)=1
              end do

c...........if the face type of the parent element is type 0
            elseif(cb.eq.0) then
              do k=1,4
                cbc(i,ne(k))=cb
              end do
            end if

          end do 
        end do 

c.......map solution from parent element to children
        call remap(ta1(1,1,1,iel),ta1(1,1,1,ref_front_id(iel)+1),
     &             ta1temp(1,1,1))
      end do
c$OMP ENDPARALLEL DO

      nelt=nelttemp+num_refine*7
      irefine=irefine+num_refine
      ntot=nelt*lx1*lx1*lx1
      return
      end

c-----------------------------------------------------------
       logical function ifcor(n1,n2,i,iface)
c-----------------------------------------------------------
c      returns whether element n1's face i and element n2's 
c      jjface(iface) have intersections, i.e. whether n1 and 
c      n2 are neighbored by an edge.
c-----------------------------------------------------------

       include 'header.h'

       integer n1,n2,i,iface
       logical ifsame

       ifcor=.false.

       if(ifsame(n1,e1v1(iface,i),n2,e2v1(iface,i)).or.
     &    ifsame(n1,e1v2(iface,i),n2,e2v2(iface,i))) then
          ifcor=.true.
       end if

       return
       end

c-----------------------------------------------------------
      logical function icheck(ie,n)
c-----------------------------------------------------------
c     Check whether element ie's three faces (sharing vertex n)
c     are nonconforming. This will prevent it from being coarsened.
c     Also check ie's neighbors on those three faces, whether ie's
c     neighbors by only an edge have a size smaller than ie's,
c     which also prevents ie from being coarsened.
c-----------------------------------------------------------
      include 'header.h'

      integer ie, n, iside, ntemp1, ntemp2, ntemp3, n1, n2, n3,
     &cb2_1,cb3_1,cb1_2,cb3_2,cb1_3,cb2_3

      icheck=.true.
      cb2_1=0
      cb3_1=0
      cb1_2=0
      cb3_2=0
      cb1_3=0
      cb2_3=0

      n1=f_c(1,n)
      n2=f_c(2,n)
      n3=f_c(3,n)
      if((cbc(n1,ie).eq.3) .or. (cbc(n2,ie).eq.3) .or.
     &   (cbc(n3,ie).eq.3)) then
         icheck=.false.
      else
        ntemp1=sje(1,1,n1,ie)
        ntemp2=sje(1,1,n2,ie)
        ntemp3=sje(1,1,n3,ie)
        if(ntemp1.ne.0)then
           cb2_1=cbc(n2,ntemp1)
           cb3_1=cbc(n3,ntemp1)
        end if
        if(ntemp2.ne.0)then
           cb3_2=cbc(n3,ntemp2)
           cb1_2=cbc(n1,ntemp2)
        end if
        if(ntemp3.ne.0)then
           cb1_3=cbc(n1,ntemp3)
           cb2_3=cbc(n2,ntemp3)
        end if
        if((cbc(n1,ie).eq.2.and.(cb2_1.eq.3.or.
     &                               cb3_1.eq.3)).or.
     &     (cbc(n2,ie).eq.2.and.(cb3_2.eq.3.or.
     &                               cb1_2.eq.3)).or.
     &     (cbc(n3,ie).eq.2.and.(cb1_3.eq.3.or.
     &                              cb2_3.eq.3)))then
          icheck=.false.
        end if
      end if

      return
      end 

c-----------------------------------------------------------
      subroutine find_coarsen(if_coarsen,neltold)
c-----------------------------------------------------------
c     Search elements to be coarsened. Check with restrictions.
c     This subroutine only checks the element itself, not its
c     neighbors.
c-----------------------------------------------------------
      
      include 'header.h'

      logical if_coarsen, iftemp, iftouch
      integer iel,i,neltold

      if_coarsen=.false.

c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(iel,i,iftemp)
c$OMP& SHARED(if_coarsen)
      do iel=1,neltold
        if(.not.skip(iel))then
          ich(iel)=0
          if(.not.iftouch(iel)) then
            iftemp=.false.
            do i=1,nsides
c.............if iel has a larger size than its face neighbors, it
c             can not be coarsened
              if(cbc(i,iel).eq.3) then
                iftemp=.true.
              endif
            enddo
            if(.not.iftemp) then
              if(.not.if_coarsen) if_coarsen=.true.
              ich(iel)=2
            end if
          end if
        endif
      enddo
c$OMP END PARALLEL DO

      return
      end

c-----------------------------------------------------------
      subroutine find_refine(if_refine)
c-----------------------------------------------------------
c     search elements to be refined based on whether they
c     have overlap with the heat source
c-----------------------------------------------------------

      include 'header.h'

      logical if_refine, iftouch
      integer iel

      if_refine=.false.

c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(iel)
c$OMP& SHARED(if_refine)
      do iel=1,nelt
        ich(iel)=0
        if(iftouch(iel)) then
          if((xc(2,iel)-xc(1,iel)).gt.dlmin) then
            if(.not.if_refine) if_refine=.true.
            ich(iel)=4
          end if
        end if
      enddo
c$OMP END PARALLEL DO

      return
      end

c-----------------------------------------------------------------
      subroutine check_refine(ifrepeat)
c-----------------------------------------------------------------
c     Check whether the potential refinement will violate the
c     restriction. If so, mark the neighbor and unmark the
c     original element, and set ifrepeat true. i.e. this procedure
c     needs to be repeated until no further check is needed
c-----------------------------------------------------------------

      include 'header.h'
 
      logical ifrepeat,ifcor
      integer iel,iface,ntemp,nntemp,i,jface

      ifrepeat=.false.

c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(iel,i,jface,ntemp,
c$OMP& iface,nntemp) SHARED(ifrepeat)
      do iel=1,nelt
c.......if iel is marked to be refined
        if(ich(iel).eq.4) then
c.........check its six faces
          do i=1,nsides
            jface=jjface(i)
            ntemp=sje(1,1,i,iel)
c...........if one face neighbor is larger in size than iel
            if(cbc(i,iel).eq.1) then
c.............unmark iel
              ich(iel)=0
c.............the large size neighbor ntemp is marked to be refined
              if(ich(ntemp).ne.4) then
                if(.not.ifrepeat) ifrepeat=.true.
                ich(ntemp)=4
              end if
c.............check iel's neighbor, neighbored by an edge on face i, which
c             must be a face neighbor of ntemp
              do iface=1,nsides
                if(iface.ne.i.and.iface.ne.jface) then
c................if edge neighbors are larger than iel, mark them to be refined
                  if(cbc(iface,ntemp).eq.2) then
                    nntemp=sje(1,1,iface,ntemp)
c..................ifcor is to make sure the edge neighbor exist
                    if(ich(nntemp).ne.4.and.
     &                 ifcor(iel,nntemp,i,iface))then
                      ich(nntemp)=4
                    end if
                  end if
                end if
              end do
c...........if face neighbor are of the same size of iel, check edge neighbors
            elseif(cbc(i,iel).eq.2)then
              do iface=1,nsides
                if(iface.ne.i.and.iface.ne.jface) then
                  if(cbc(iface,ntemp).eq.1)then
                    nntemp=sje(1,1,iface,ntemp)
                    ich(nntemp)=4
                    ich(iel)=0
                    if(.not.ifrepeat) ifrepeat=.true.
                  end if
                end if
              end do
            end if
          enddo
        end if
      end do
c$OMP END PARALLEL DO

      return
      end

c-----------------------------------------------------------------
      logical function iftouch(iel)
c-----------------------------------------------------------------
c     check whether element iel has overlap with the heat source
c-----------------------------------------------------------------

      include 'header.h'

      double precision dis, dis1, dis2, dis3, alpha2
      integer iel

      alpha2 = alpha*alpha

      if     (x0 .lt. xc(1,iel)) then
        dis1 = xc(1,iel) - x0
      elseif (x0 .gt. xc(2,iel)) then
        dis1 = x0 - xc(2,iel)
      else
        dis1 = 0.d0
      endif

      if     (y0 .lt. yc(1,iel)) then
        dis2 = yc(1,iel) - y0
      elseif (y0 .gt. yc(3,iel)) then
        dis2 = y0 - yc(3,iel)
      else
        dis2 = 0.d0
      endif

      if     (z0 .lt. zc(1,iel)) then
        dis3 = zc(1,iel) - z0
      elseif (z0 .gt. zc(5,iel)) then
        dis3 = z0 - zc(5,iel)
      else
       dis3 = 0.d0
      endif

      dis = dis1**2+dis2**2+dis3**2

      if (dis .lt. alpha2) then
       iftouch=.true.
      else
       iftouch=.false.
      end if

      return
      end


c-----------------------------------------------------------------
      subroutine remap (y,y1,x) 
c-----------------------------------------------------------------
c     After a refinement, map the solution  from the parent (x) to
c     the eight children. y is the solution on the first child
c     (front-bottom-left) and y1 is the solution on the next 7 
c     children.
c-----------------------------------------------------------------

      include 'header.h'

      double precision x(lx1,lx1,lx1),y(lx1,lx1,lx1),y1(lx1,lx1,lx1,7),
     &       yone(lx1,lx1,lx1,2), ytwo(lx1,lx1,lx1,4)
      integer i, iz, ii, jj, kk

      call r_init(y,lx1*lx1*lx1,0.d0)
      call r_init(y1,lx1*lx1*lx1*7,0.d0)
      call r_init(yone,lx1*lx1*lx1*2,0.d0)
      call r_init(ytwo,lx1*lx1*lx1*4,0.d0)

      do  i=1,lx1
        do kk = 1, lx1
          do jj = 1, lx1
            do ii = 1, lx1
              yone(ii,jj,i,1) = yone(ii,jj,i,1) +ixmc1(ii,kk)*x(kk,jj,i)
              yone(ii,jj,i,2) = yone(ii,jj,i,2) +ixmc2(ii,kk)*x(kk,jj,i)
            end do
          end do
        end do

        do kk = 1, lx1
          do jj = 1, lx1
            do ii = 1, lx1
              ytwo(ii,i,jj,1) = ytwo(ii,i,jj,1) + 
     &                          yone(ii,kk,i,1)*ixtmc1(kk,jj)
              ytwo(ii,i,jj,2) = ytwo(ii,i,jj,2) + 
     &                          yone(ii,kk,i,1)*ixtmc2(kk,jj)
              ytwo(ii,i,jj,3) = ytwo(ii,i,jj,3) + 
     &                          yone(ii,kk,i,2)*ixtmc1(kk,jj)
              ytwo(ii,i,jj,4) = ytwo(ii,i,jj,4) + 
     &                          yone(ii,kk,i,2)*ixtmc2(kk,jj)
            end do
          end do
        end do
      end do

      do  iz=1,lx1
        do kk = 1, lx1
          do jj = 1, lx1
            do ii = 1, lx1
              y(ii,iz,jj) = y(ii,iz,jj) +
     &                        ytwo(ii,kk,iz,1)*ixtmc1(kk,jj)
              y1(ii,iz,jj,1) = y1(ii,iz,jj,1) +
     &                        ytwo(ii,kk,iz,3)*ixtmc1(kk,jj)
              y1(ii,iz,jj,2) = y1(ii,iz,jj,2) +
     &                        ytwo(ii,kk,iz,2)*ixtmc1(kk,jj)
              y1(ii,iz,jj,3) = y1(ii,iz,jj,3) +
     &                        ytwo(ii,kk,iz,4)*ixtmc1(kk,jj)
              y1(ii,iz,jj,4) = y1(ii,iz,jj,4) +
     &                        ytwo(ii,kk,iz,1)*ixtmc2(kk,jj)
              y1(ii,iz,jj,5) = y1(ii,iz,jj,5) +
     &                        ytwo(ii,kk,iz,3)*ixtmc2(kk,jj)
              y1(ii,iz,jj,6) = y1(ii,iz,jj,6) +
     &                        ytwo(ii,kk,iz,2)*ixtmc2(kk,jj)
              y1(ii,iz,jj,7) = y1(ii,iz,jj,7) +
     &                        ytwo(ii,kk,iz,4)*ixtmc2(kk,jj)            
            end do
          end do
        end do
      end do

      return
      end


c=======================================================================
      subroutine merging(iela)
c-----------------------------------------------------------------------
c     This subroutine is to merge the eight child elements and map 
c     the solution from eight children to the  merged element. 
c     iela array records the eight elements to be merged.
c-----------------------------------------------------------------------

      include 'header.h'

      double precision x1,x2,y1,y2,z1,z2
      integer ielnew,i,ntemp,jface,ii,cb,ntempa(4),iela(8),ielold,
     &        ntema(4)

      ielnew=iela(1)

      tree(ielnew)=ishft(tree(ielnew),-3)   

c.....element vertices 
      x1=xc(1,iela(1))
      x2=xc(2,iela(2))
      y1=yc(1,iela(1))
      y2=yc(3,iela(3))
      z1=zc(1,iela(1))
      z2=zc(5,iela(5))

      do i=1,7,2
        xc(i,ielnew)=x1
      end do
      do i=2,8,2
        xc(i,ielnew)=x2
      end do
      do i=1,2
        yc(i,ielnew)=y1
        yc(i+4,ielnew)=y1
      end do
      do i=3,4
        yc(i,ielnew)=y2
        yc(i+4,ielnew)=y2
      end do
      do i=1,4
        zc(i,ielnew)=z1
      end do
      do i=5,8
        zc(i,ielnew)=z2
      end do

c.....update neighboring information
      do i=1,nsides
        jface=jjface(i)
        ielold=iela(children(1,i))
        do ii=1,4
          ntempa(ii)=iela(children(ii,i))
        end do

        cb=cbc(i,ielold)
       
        if (cb.eq.2) then
c.........if the neighbor elements also will be coarsened
          if(ifcoa_id(sje(1,1,i,ielold)))then
            if (i.eq.2 .or. i.eq. 4 .or. i.eq.6) then
              ntemp=sje(1,1,i,sje(1,1,i,ntempa(1)))
            else
              ntemp=sje(1,1,i,ntempa(1))
            end if 
            sje(1,1,i,ielnew)=ntemp
            ijel(1,i,ielnew)=1
            ijel(2,i,ielnew)=1
            cbc(i,ielnew)=2

c.........if the neighbor elements will not be coarsened
          else
            do ii=1,4
              ntema(ii)=sje(1,1,i,ntempa(ii)) 
              cbc(jface,ntema(ii))=1
              sje(1,1,jface,ntema(ii))=ielnew
              ijel(1,jface,ntema(ii))=iijj(1,ii)
              ijel(2,jface,ntema(ii))=iijj(2,ii)
              sje(iijj(1,ii),iijj(2,ii),i,ielnew)=ntema(ii)
              ijel(1,i,ielnew)=1
              ijel(2,i,ielnew)=1
            end do
            cbc(i,ielnew)=3
          end if       

        else if(cb.eq.1)then

          ntemp=sje(1,1,i,ielold)
          cbc(jface,ntemp)=2
          ijel(1,jface,ntemp)=1
          ijel(2,jface,ntemp)=1
          sje(1,1,jface,ntemp)=ielnew
          sje(1,2,jface,ntemp)=0
          sje(2,1,jface,ntemp)=0
          sje(2,2,jface,ntemp)=0
           
          cbc(i,ielnew)=2
          ijel(1,i,ielnew)=1
          ijel(2,i,ielnew)=1
          sje(1,1,i,ielnew)=ntemp
         
        else if(cb.eq.0)then
          cbc(i,ielnew)=0
          sje(1,1,i,ielnew)=0
          sje(1,2,i,ielnew)=0
          sje(2,1,i,ielnew)=0
          sje(2,2,i,ielnew)=0
        endif

      end do

c.....map solution from children to the merged element
      call remap2(iela, ielnew)
      
      return
      end

c-----------------------------------------------------------------
      subroutine remap2(iela, ielnew)
c-----------------------------------------------------------------
c     Map the solution from the children to the parent.
c     iela array records the eight elements to be merged.
c     ielnew is the element index of the merged element.
c-----------------------------------------------------------------

      include 'header.h'
      integer iela(8), ielnew

      double precision temp1(lx1,lx1,lx1),
     &       temp2(lx1,lx1,lx1),temp3(lx1,lx1,lx1),temp4(lx1,lx1,lx1),
     &       temp5(lx1,lx1,lx1),temp6(lx1,lx1,lx1)

      call remapx(ta1(1,1,1,iela(1)),ta1(1,1,1,iela(2)),temp1)
      call remapx(ta1(1,1,1,iela(3)),ta1(1,1,1,iela(4)),temp2)
      call remapx(ta1(1,1,1,iela(5)),ta1(1,1,1,iela(6)),temp3)
      call remapx(ta1(1,1,1,iela(7)),ta1(1,1,1,iela(8)),temp4)
      call remapy(temp1,temp2,temp5)
      call remapy(temp3,temp4,temp6)
      call remapz(temp5,temp6,ta1(1,1,1,ielnew))

      return
      end       

c-----------------------------------------------------------------
      subroutine remapz(x1,x2,y)
c-----------------------------------------------------------------
c     z direction mapping after the merge.
c     Map solution from x1 & x2 to y.
c-----------------------------------------------------------------

      include 'header.h'

      double precision x1(lx1,lx1,lx1),x2(lx1,lx1,lx1),y(lx1,lx1,lx1)
      integer ix, iy, ip

      do iy=1,lx1
        do ix=1,lx1
          y(ix,iy,1)=x1(ix,iy,1)

          y(ix,iy,2)=0.d0
          do ip=1,lx1
            y(ix,iy,2)=y(ix,iy,2)+map2(ip)*x1(ix,iy,ip)
          end do

          y(ix,iy,3)=x1(ix,iy,lx1)

          y(ix,iy,4)=0.d0
          do ip=1,lx1
            y(ix,iy,4)=y(ix,iy,4)+map4(ip)*x2(ix,iy,ip)
          end do

          y(ix,iy,lx1)=x2(ix,iy,lx1)
        end do
      end do

      return
      end      

c-----------------------------------------------------------------
      subroutine remapy(x1,x2,y)
c-----------------------------------------------------------------
c     y direction mapping after the merge.
c     Map solution from x1 & x2 to y.
c-----------------------------------------------------------------

      include 'header.h'

      double precision x1(lx1,lx1,lx1),x2(lx1,lx1,lx1),y(lx1,lx1,lx1)
      integer ix, iz, ip

      do iz=1,lx1
        do ix=1,lx1
          y(ix,1,iz)=x1(ix,1,iz)

          y(ix,2,iz)=0.d0
          do ip=1,lx1
            y(ix,2,iz)=y(ix,2,iz)+map2(ip)*x1(ix,ip,iz)
          end do

          y(ix,3,iz)=x1(ix,lx1,iz)

          y(ix,4,iz)=0.d0
          do ip=1,lx1
            y(ix,4,iz)=y(ix,4,iz)+map4(ip)*x2(ix,ip,iz)
          end do

          y(ix,lx1,iz)=x2(ix,lx1,iz)
        end do
      end do

      return
      end      

c-----------------------------------------------------------------
      subroutine remapx(x1,x2,y)
c-----------------------------------------------------------------
c     x direction mapping after the merge.
c     Map solution from x1 & x2 to y.
c-----------------------------------------------------------------

      include 'header.h'

      double precision x1(lx1,lx1,lx1),x2(lx1,lx1,lx1),y(lx1,lx1,lx1)
      integer iy, iz, ip

      do iz=1,lx1
        do iy=1,lx1
          y(1,iy,iz)=x1(1,iy,iz)

          y(2,iy,iz)=0.d0
          do ip=1,lx1
            y(2,iy,iz)=y(2,iy,iz)+map2(ip)*x1(ip,iy,iz)
          end do

          y(3,iy,iz)=x1(lx1,iy,iz)

          y(4,iy,iz)=0.d0
          do ip=1,lx1
            y(4,iy,iz)=y(4,iy,iz)+map4(ip)*x2(ip,iy,iz)
          end do

          y(lx1,iy,iz)=x2(lx1,iy,iz)
        end do
      end do

      return
      end