c-----------------------------------------------------------------
      subroutine create_initial_grid        
c------------------------------------------------------------------
    
      include 'header.h'

      integer i

      nelt=1
      ntot=nelt*lx1*lx1*lx1 
      tree(1)=1
      mt_to_id(1)=1
      do i=1,7,2
        xc(i,1)=0.d0
        xc(i+1,1)=1.d0
      end do

      do i=1,2
        yc(i,1)=0.d0
        yc(2+i,1)=1.d0
        yc(4+i,1)=0.d0
        yc(6+i,1)=1.d0
      end do
     
      do i=1,4
        zc(i,1)=0.d0
        zc(4+i,1)=1.d0
      end do
  
      do i=1,6
        cbc(i,1)=0
      end do

      return

      end

c-----------------------------------------------------------------
      subroutine coef
c-----------------------------------------------------------------
c
c     generate 
c
c            - collocation points
c            - weights
c            - derivative matrices 
c            - projection matrices
c            - interpolation matrices 
c
c     associated with the 
c
c            - gauss-legendre lobatto mesh (suffix m1)
c
c----------------------------------------------------------------

      include 'header.h'

      integer i,j,k

c.....for gauss-legendre lobatto mesh (suffix m1)
c.....generate collocation points and weights 

      zgm1(1)=-1.d0
      zgm1(2)=-0.6546536707079771d0
      zgm1(3)=0.d0
      zgm1(4)= 0.6546536707079771d0
      zgm1(5)=1.d0
      wxm1(1)=0.1d0
      wxm1(2)=49.d0/90.d0
      wxm1(3)=32.d0/45.d0
      wxm1(4)=wxm1(2)
      wxm1(5)=0.1d0 

      do k=1,lx1
        do j=1,lx1
          do i=1,lx1
            w3m1(i,j,k)=wxm1(i)*wxm1(j)*wxm1(k)
          end do
        end do
      end do

c.....generate derivative matrices

      dxm1(1,1)=-5.0d0
      dxm1(2,1)=-1.240990253030982d0
      dxm1(3,1)= 0.375d0
      dxm1(4,1)=-0.2590097469690172d0
      dxm1(5,1)= 0.5d0
      dxm1(1,2)= 6.756502488724238d0
      dxm1(2,2)= 0.d0
      dxm1(3,2)=-1.336584577695453d0
      dxm1(4,2)= 0.7637626158259734d0
      dxm1(5,2)=-1.410164177942427d0
      dxm1(1,3)=-2.666666666666667d0
      dxm1(2,3)= 1.745743121887939d0
      dxm1(3,3)= 0.d0
      dxm1(4,3)=-dxm1(2,3)
      dxm1(5,3)=-dxm1(1,3)
      do j=4,lx1
        do i=1,lx1
          dxm1(i,j)=-dxm1(lx1+1-i,lx1+1-j)
        end do
      end do
      do j=1,lx1
        do i=1,lx1
          dxtm1(i,j)=dxm1(j,i)
        end do
      end do

c.....generate projection (mapping) matrices

      qbnew(1,1,1)=-0.1772843218615690d0
      qbnew(2,1,1)=9.375d-02
      qbnew(3,1,1)=-3.700139242414530d-02
      qbnew(1,2,1)= 0.7152146412463197d0
      qbnew(2,2,1)=-0.2285757930375471d0
      qbnew(3,2,1)= 8.333333333333333d-02
      qbnew(1,3,1)= 0.4398680650316104d0
      qbnew(2,3,1)= 0.2083333333333333d0
      qbnew(3,3,1)=-5.891568407922938d-02
      qbnew(1,4,1)= 8.333333333333333d-02
      qbnew(2,4,1)= 0.3561799597042137d0
      qbnew(3,4,1)=-4.854797457965334d-02
      qbnew(1,5,1)= 0.d0
      qbnew(2,5,1)=7.03125d-02
      qbnew(3,5,1)=0.d0
      
      do j=1,lx1
        do i=1,3
          qbnew(i,j,2)=qbnew(4-i,lx1+1-j,1)
        end do
      end do 

c.....generate interpolation matrices for mesh refinement

      ixtmc1(1,1)=1.d0
      ixtmc1(2,1)=0.d0
      ixtmc1(3,1)=0.d0
      ixtmc1(4,1)=0.d0
      ixtmc1(5,1)=0.d0 
      ixtmc1(1,2)= 0.3385078435248143d0
      ixtmc1(2,2)= 0.7898516348912331d0
      ixtmc1(3,2)=-0.1884018684471238d0
      ixtmc1(4,2)= 9.202967302175333d-02
      ixtmc1(5,2)=-3.198728299067715d-02
      ixtmc1(1,3)=-0.1171875d0
      ixtmc1(2,3)= 0.8840317166357952d0
      ixtmc1(3,3)= 0.3125d0    
      ixtmc1(4,3)=-0.118406716635795d0 
      ixtmc1(5,3)= 0.0390625d0   
      ixtmc1(1,4)=-7.065070066767144d-02
      ixtmc1(2,4)= 0.2829703269782467d0 
      ixtmc1(3,4)= 0.902687582732838d0
      ixtmc1(4,4)=-0.1648516348912333d0 
      ixtmc1(5,4)= 4.984442584781999d-02
      ixtmc1(1,5)=0.d0
      ixtmc1(2,5)=0.d0
      ixtmc1(3,5)=1.d0 
      ixtmc1(4,5)=0.d0
      ixtmc1(5,5)=0.d0  
      do j=1,lx1
        do i=1,lx1
          ixmc1(i,j)=ixtmc1(j,i)
        end do
      end do

      do j=1,lx1
        do i=1,lx1
          ixtmc2(i,j)=ixtmc1(lx1+1-i,lx1+1-j)
        end do
      end do

      do j=1,lx1
        do i=1,lx1
          ixmc2(i,j)=ixtmc2(j,i)
        end do
      end do

c.....solution interpolation matrix for mesh coarsening

      map2(1)=-0.1179652785083428d0
      map2(2)= 0.5505046330389332d0
      map2(3)= 0.7024534364259963d0
      map2(4)=-0.1972224518285866d0
      map2(5)= 6.222966087199998d-02

      do i=1,lx1
        map4(i)=map2(lx1+1-i)
      end do

      return
      end

c-------------------------------------------------------------------
      subroutine geom1
c-------------------------------------------------------------------
c
c     routine to generate elemental geometry information on mesh m1,
c     (gauss-legendre lobatto mesh).
c
c         xrm1_s   -   dx/dr, dy/dr, dz/dr
c         rxm1_s   -   dr/dx, dr/dy, dr/dz
c         g1m1_s  geometric factors used in preconditioner computation
c         g4m1_s  g5m1_s  g6m1_s :
c         geometric factors used in lapacian opertor
c         jacm1    -   jacobian
c         bm1      -   mass matrix
c         xfrac    -   will be used in prepwork for calculating collocation
c                      coordinates
c         idel     -   collocation points index on element boundaries 
c------------------------------------------------------------------

      include 'header.h'

      double precision temp,temp1,temp2,dtemp
      integer isize,i,j,k,ntemp,iel
 
      do i=1,lx1
        xfrac(i)=zgm1(i)*0.5d0 + 0.5d0
      end do

c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(ISIZE,TEMP,TEMP1,TEMP2,
c$OMP&  K,J,I,dtemp)
      do isize=1,refine_max
        temp=2.d0**(-isize-1)
        dtemp=1.d0/temp
        temp1=temp**3
        temp2=temp**2
        do k=1,lx1
          do j=1,lx1
            do i=1,lx1
              xrm1_s(i,j,k,isize)=dtemp
              jacm1_s(i,j,k,isize)=temp1
              rxm1_s(i,j,k,isize)=temp2
              g1m1_s(i,j,k,isize)=w3m1(i,j,k)*temp
              bm1_s(i,j,k,isize)=w3m1(i,j,k)*temp1
              g4m1_s(i,j,k,isize)=g1m1_s(i,j,k,isize)/wxm1(i)
              g5m1_s(i,j,k,isize)=g1m1_s(i,j,k,isize)/wxm1(j)
              g6m1_s(i,j,k,isize)=g1m1_s(i,j,k,isize)/wxm1(k)
            end do
          end do
        end do
      end do
c$OMP END PARALLEL DO

c$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(ntemp,i,j,iel)
      do iel = 1, lelt
        ntemp=lx1*lx1*lx1*(iel-1)
        do j = 1, lx1
          do i = 1, lx1
            idel(i,j,1,iel)=ntemp+(i-1)*lx1 + (j-1)*lx1*lx1+lx1
            idel(i,j,2,iel)=ntemp+(i-1)*lx1 + (j-1)*lx1*lx1+1
            idel(i,j,3,iel)=ntemp+(i-1)*1 + (j-1)*lx1*lx1+lx1*(lx1-1)+1
            idel(i,j,4,iel)=ntemp+(i-1)*1 + (j-1)*lx1*lx1+1
            idel(i,j,5,iel)=ntemp+(i-1)*1 + (j-1)*lx1+lx1*lx1*(lx1-1)+1
            idel(i,j,6,iel)=ntemp+(i-1)*1 + (j-1)*lx1+1
          end do
        end do
      end do
c$OMP END PARALLEL DO

      return
      end

c------------------------------------------------------------------
      subroutine setdef
c------------------------------------------------------------------
c     compute the discrete laplacian operators
c------------------------------------------------------------------

      include 'header.h'

      integer i,j,ip
 
      call r_init(wdtdr(1,1),lx1*lx1,0.d0)

      do i=1,lx1
        do j=1,lx1
          do ip=1,lx1
            wdtdr(i,j) = wdtdr(i,j) + wxm1(ip)*dxm1(ip,i)*dxm1(ip,j)
          end do
        end do
      end do

      return 
      end


c------------------------------------------------------------------
      subroutine prepwork
c------------------------------------------------------------------
c     mesh information preparations: calculate refinement levels of
c     each element, mask matrix for domain boundary and element 
c     boundaries
c------------------------------------------------------------------

      include 'header.h'

      integer i, j, iel, iface, cb
      double precision rdlog2

      ntot = nelt*nxyz
      rdlog2 = 1.d0/dlog(2.d0)

c$OMP PARALLEL DEFAULT(SHARED) PRIVATE(I,J,IEL,IFACE,CB)

c.....calculate the refinement levels of each element

c$OMP DO 
      do iel = 1, nelt
        size_e(iel)=-dlog(xc(2,iel)-xc(1,iel))*rdlog2+1.d-8
      end do
c$OMP END DO nowait

c.....mask matrix for element boundary

c$OMP DO
      do iel = 1, nelt
        call r_init(tmult(1,1,1,iel),nxyz,1.d0)   
        do iface=1,nsides
          call facev(tmult(1,1,1,iel),iface,0.0d0)
        end do
      end do
c$OMP END DO nowait

c.....masks for domain boundary at mortar 

c$OMP DO
      do iel=1,nmor
        tmmor(iel)=1.d0
      end do
c$OMP END DO

c$OMP DO
      do iel = 1, nelt
        do iface = 1,nsides
          cb=cbc(iface,iel)
          if(cb.eq.0) then
            do j=2,lx1-1
              do i=2,lx1-1
               tmmor(idmo(i,j,1,1,iface,iel))=0.d0
              end do
            end do

            j=1
            do i = 1, lx1-1
               tmmor(idmo(i,j,1,1,iface,iel))=0.d0
            end do

            if(idmo(lx1,1,1,1,iface,iel).eq.0)then
              tmmor(idmo(lx1,1,1,2,iface,iel))=0.d0
            else
              tmmor(idmo(lx1,1,1,1,iface,iel))=0.d0
              do i=1,lx1
                tmmor(idmo(i,j,1,2,iface,iel))=0.d0
              end do
            end if

            i=lx1
            if(idmo(lx1,2,1,2,iface,iel).eq.0)then
              do j=2,lx1-1
                tmmor(idmo(i,j,1,1,iface,iel))=0.d0
              end do
              tmmor(idmo(lx1,lx1,2,2,iface,iel))=0.d0
            else
              do j=2,lx1
                tmmor(idmo(i,j,1,2,iface,iel))=0.d0
              end do
              do j=1,lx1
                tmmor(idmo(i,j,2,2,iface,iel))=0.d0
              end do
            end if
            
            j=lx1
            tmmor(idmo(1,lx1,2,1,iface,iel))=0.d0
            if(idmo(2,lx1,2,1,iface,iel).eq.0)then
              do i=2,lx1-1
                tmmor(idmo(i,j,1,1,iface,iel))=0.d0
              end do
            else
              do i=2,lx1
                tmmor(idmo(i,j,2,1,iface,iel))=0.d0
              end do
              do i=1,lx1-1
                tmmor(idmo(i,j,2,2,iface,iel))=0.d0
              end do
            end if

            i=1
            do j=2,lx1-1
             tmmor(idmo(i,j,1,1,iface,iel))=0.d0
            end do
            if(idmo(1,lx1,1,1,iface,iel).ne.0)then
              tmmor(idmo(i,lx1,1,1,iface,iel))=0.d0
              do j=1,lx1-1
               tmmor(idmo(i,j,2,1,iface,iel))=0.d0
              end do
            end if

          endif
        end do
       end do
c$OMP END DO nowait
            
c$OMP END PARALLEL
      return
      end 
    

c------------------------------------------------------------------
      block data top_constants

c------------------------------------------------------------------
c.....We store some tables of useful topological constants
c------------------------------------------------------------------
      include 'header.h'

c     f_e_ef(e,f) returns the other face sharing the e'th local edge of face f.
      data f_e_ef/6,3,5,4, 6,3,5,4, 6,1,5,2, 6,1,5,2, 4,1,3,2, 4,1,3,2/

c.....e_c(n,j) returns n'th edge sharing the vertex j of an element
      data e_c /5,8,11, 1,4,11,  5,6,9, 1,2,9,
     &          7,8,12, 3,4,12, 6,7,10, 2,3,10/

c.....local_corner(n,i) returns the local corner index of vertex n on face i
      data local_corner /0,1,0,2,0,3,0,4, 1,0,2,0,3,0,4,0,
     &                   0,0,1,2,0,0,3,4, 1,2,0,0,3,4,0,0,
     &                   0,0,0,0,1,2,3,4, 1,2,3,4,0,0,0,0/

c.....cal_nnb(n,i) returns the neighbor elements neighbored by n'th edge
c     among the three edges sharing vertex i
c     the elements are the eight children elements ordered as 1 to 8.
      data cal_nnb/5,2,3, 6,1,4, 7,4,1, 8,3,2,
     &             1,6,7, 2,5,8, 3,8,5, 4,7,6/

c.....returns the opposite local corner index: 1-4,2-3
      data oplc /4,3,2,1/

c.....cal_iijj(i,n) returns the location of local corner number n on a face 
c     i =1  to get ii, i=2 to get jj
c     (ii,jj) is defined the same as in mortar location (ii,jj)
      data cal_iijj /1,1, 1,2, 2,1, 2,2/

c.....returns the adjacent(neighbored by a face) element's children,
c     assumming a vertex is shared by eight child elements 1-8. 
c     index n is local corner number on the face which is being 
c     assigned the mortar index number
      data cal_intempx /8,6,4,2, 7,5,3,1, 8,7,4,3, 
     $                  6,5,2,1, 8,7,6,5, 4,3,2,1/

c.....c_f(i,f) returns the vertex number of i'th local corner on face f
      data c_f /2,4,6,8, 1,3,5,7, 3,4,7,8, 1,2,5,6, 5,6,7,8, 1,2,3,4/

c.....on each face of the parent element, there are four children element.
c     le_arr(i,j,n) returns the i'th elements among the four children elements 
c     n refers to the direction: 1 for x, 2 for y and 3 for z direction. 
c     j refers to positive(0) or negative(1) direction on x, y or z direction.
c     n=1,j=0 refers to face 1 and n=1, j=1 refers to face 2, n=2,j=0 refers to
c     face 3.... 
c     The current eight children are ordered as 8,1,2,3,4,5,6,7 
      data    le_arr/8,2,4,6, 1,3,5,7, 
     $               8,1,4,5, 2,3,6,7, 
     $               8,1,2,3, 4,5,6,7/

c.....jjface(n) returns the face opposite to face n
      data jjface /2,1,4,3,6,5/

cc.....edgeface(n,f) returns OTHER face which shares local edge n on face f
c      integer edgeface(4,6)
c      data edgeface /6,3,5,4, 6,3,5,4, 6,1,5,2, 
c     $               6,1,5,2, 4,1,3,2, 4,1,3,2/

c.....e_face2(n,f) returns the local edge number of edge n on the
c     other face sharing local edge n on face f
      data e_face2 /2,2,2,2, 4,4,4,4, 3,2,3,2, 
     $              1,4,1,4, 3,3,3,3, 1,1,1,1/

c.....op(n) returns the local edge number of the edge which 
c     is opposite to local edge n on the same face
      data op /3,4,1,2/

c.....localedgenumber(f,e) returns the local edge number for edge e
c     on face f. A zero result value signifies illegal input
      data localedgenumber /1,0,0,0,0,2, 2,0,2,0,0,0, 3,0,0,0,2,0, 
     $                      4,0,0,2,0,0, 0,1,0,0,0,4, 0,2,4,0,0,0, 
     $                      0,3,0,0,4,0, 0,4,0,4,0,0, 0,0,1,0,0,3, 
     $                      0,0,3,0,3,0, 0,0,0,1,0,1, 0,0,0,3,1,0/

c.....edgenumber(e,f) returns the edge index of local edge e on face f
      data edgenumber / 1,2, 3,4,  5,6, 7,8,  9,2,10,6, 
     $                 11,4,12,8, 12,3,10,7, 11,1, 9,5/

c.....f_c(c,n) returns the face index of i'th face sharing vertex n 
      data f_c /2,4,6, 1,4,6, 2,3,6, 1,3,6,
     &          2,4,5, 1,4,5, 2,3,5, 1,3,5/

c.....if two elements are neighbor by one edge, 
c     e1v1(f1,f2) returns the smaller index of the two vertices on this 
c     edge on one element
c     e1v2 returns the larger index of the two vertices of this edge on 
c     on element. exfor a vertex on element 
c     e2v1 returns the smaller index of the two vertices on this edge on 
c     another element
c     e2v2 returns the larger index of the two vertiex on this edge on
c     another element
      data e1v1/0,0,4,2,6,2, 0,0,3,1,5,1, 4,3,0,0,7,3,
     &          2,1,0,0,5,1, 6,5,7,5,0,0, 2,1,3,1,0,0/
      data e2v1/0,0,1,3,1,5, 0,0,2,4,2,6, 1,2,0,0,1,5,
     &          3,4,0,0,3,7, 1,2,1,3,0,0, 5,6,5,7,0,0/
      data e1v2/0,0,8,6,8,4, 0,0,7,5,7,3, 8,7,0,0,8,4,
     &          6,5,0,0,6,2, 8,7,8,6,0,0, 4,3,4,2,0,0/
      data e2v2/0,0,5,7,3,7, 0,0,6,8,4,8, 5,6,0,0,2,6,
     &          7,8,0,0,4,8, 3,4,2,4,0,0, 7,8,6,8,0,0/

c.....children(n1,n)returns the four elements among the eight children 
c     elements to be merged on face n of the parent element
c     the IDs for the eight children are 1,2,3,4,5,6,7,8
      data children/2,4,6,8, 1,3,5,7, 3,4,7,8, 
     &              1,2,5,6, 5,6,7,8, 1,2,3,4/

c.....iijj(n1,n) returns the location of n's mortar on an element face
c     n1=1 refers to x direction location and n1=2 refers to y direction
      data iijj/1,1,1,2,2,1,2,2/

c.....v_end(n) returns the index of collocation points at two ends of each
c     direction
      data v_end /1,lx1/

c.....face_l1,face_l2,face_ld return for start,end,stride for a loop over faces 
c     used on subroutine  mortar_vertex
      data face_l1 /2,3,1/, face_l2 /3,1,2/, face_ld /1,-2,1/

      end