c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine  initialize

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     This subroutine initializes the field variable u using 
c     tri-linear transfinite interpolation of the boundary values     
c---------------------------------------------------------------------

      include 'header.h'
      
      integer c, i, j, k, m, ii, jj, kk, ix, iy, iz
      double precision  xi, eta, zeta, Pface(5,3,2), Pxi, Peta, 
     >     Pzeta, temp(5)

c---------------------------------------------------------------------
c  Later (in compute_rhs) we compute 1/u for every element. A few of 
c  the corner elements are not used, but it convenient (and faster) 
c  to compute the whole thing with a simple loop. Make sure those 
c  values are nonzero by initializing the whole thing here. 
c---------------------------------------------------------------------
      do c = 1, ncells
         do kk = -1, KMAX
            do jj = -1, JMAX
               do ii = -1, IMAX
                  do m = 1, 5
                     u(m, ii, jj, kk, c) = 1.0
                  end do
               end do
            end do
         end do
      end do
c---------------------------------------------------------------------


c---------------------------------------------------------------------
c     first store the "interpolated" values everywhere on the grid    
c---------------------------------------------------------------------
      do c=1, ncells
         kk = 0
         do k = cell_low(3,c), cell_high(3,c)
            zeta = dble(k) * dnzm1
            jj = 0
            do j = cell_low(2,c), cell_high(2,c)
               eta = dble(j) * dnym1
               ii = 0
               do i = cell_low(1,c), cell_high(1,c)
                  xi = dble(i) * dnxm1
                  
                  do ix = 1, 2
                     call exact_solution(dble(ix-1), eta, zeta, 
     >                    Pface(1,1,ix))
                  enddo

                  do iy = 1, 2
                     call exact_solution(xi, dble(iy-1) , zeta, 
     >                    Pface(1,2,iy))
                  enddo

                  do iz = 1, 2
                     call exact_solution(xi, eta, dble(iz-1),   
     >                    Pface(1,3,iz))
                  enddo

                  do m = 1, 5
                     Pxi   = xi   * Pface(m,1,2) + 
     >                    (1.0d0-xi)   * Pface(m,1,1)
                     Peta  = eta  * Pface(m,2,2) + 
     >                    (1.0d0-eta)  * Pface(m,2,1)
                     Pzeta = zeta * Pface(m,3,2) + 
     >                    (1.0d0-zeta) * Pface(m,3,1)
                     
                     u(m,ii,jj,kk,c) = Pxi + Peta + Pzeta - 
     >                    Pxi*Peta - Pxi*Pzeta - Peta*Pzeta + 
     >                    Pxi*Peta*Pzeta

                  enddo
                  ii = ii + 1
               enddo
               jj = jj + 1
            enddo
            kk = kk+1
         enddo
      enddo

c---------------------------------------------------------------------
c     now store the exact values on the boundaries        
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c     west face                                                  
c---------------------------------------------------------------------
      c = slice(1,1)
      ii = 0
      xi = 0.0d0
      kk = 0
      do k = cell_low(3,c), cell_high(3,c)
         zeta = dble(k) * dnzm1
         jj = 0
         do j = cell_low(2,c), cell_high(2,c)
            eta = dble(j) * dnym1
            call exact_solution(xi, eta, zeta, temp)
            do m = 1, 5
               u(m,ii,jj,kk,c) = temp(m)
            enddo
            jj = jj + 1
         enddo
         kk = kk + 1
      enddo

c---------------------------------------------------------------------
c     east face                                                      
c---------------------------------------------------------------------
      c  = slice(1,ncells)
      ii = cell_size(1,c)-1
      xi = 1.0d0
      kk = 0
      do k = cell_low(3,c), cell_high(3,c)
         zeta = dble(k) * dnzm1
         jj = 0
         do j = cell_low(2,c), cell_high(2,c)
            eta = dble(j) * dnym1
            call exact_solution(xi, eta, zeta, temp)
            do m = 1, 5
               u(m,ii,jj,kk,c) = temp(m)
            enddo
            jj = jj + 1
         enddo
         kk = kk + 1
      enddo

c---------------------------------------------------------------------
c     south face                                                 
c---------------------------------------------------------------------
      c = slice(2,1)
      jj = 0
      eta = 0.0d0
      kk = 0
      do k = cell_low(3,c), cell_high(3,c)
         zeta = dble(k) * dnzm1
         ii = 0
         do i = cell_low(1,c), cell_high(1,c)
            xi = dble(i) * dnxm1
            call exact_solution(xi, eta, zeta, temp)
            do m = 1, 5
               u(m,ii,jj,kk,c) = temp(m)
            enddo
            ii = ii + 1
         enddo
         kk = kk + 1
      enddo


c---------------------------------------------------------------------
c     north face                                    
c---------------------------------------------------------------------
      c = slice(2,ncells)
      jj = cell_size(2,c)-1
      eta = 1.0d0
      kk = 0
      do k = cell_low(3,c), cell_high(3,c)
         zeta = dble(k) * dnzm1
         ii = 0
         do i = cell_low(1,c), cell_high(1,c)
            xi = dble(i) * dnxm1
            call exact_solution(xi, eta, zeta, temp)
            do m = 1, 5
               u(m,ii,jj,kk,c) = temp(m)
            enddo
            ii = ii + 1
         enddo
         kk = kk + 1
      enddo

c---------------------------------------------------------------------
c     bottom face                                       
c---------------------------------------------------------------------
      c = slice(3,1)
      kk = 0
      zeta = 0.0d0
      jj = 0
      do j = cell_low(2,c), cell_high(2,c)
         eta = dble(j) * dnym1
         ii = 0
         do i =cell_low(1,c), cell_high(1,c)
            xi = dble(i) *dnxm1
            call exact_solution(xi, eta, zeta, temp)
            do m = 1, 5
               u(m,ii,jj,kk,c) = temp(m)
            enddo
            ii = ii + 1
         enddo
         jj = jj + 1
      enddo

c---------------------------------------------------------------------
c     top face     
c---------------------------------------------------------------------
      c = slice(3,ncells)
      kk = cell_size(3,c)-1
      zeta = 1.0d0
      jj = 0
      do j = cell_low(2,c), cell_high(2,c)
         eta = dble(j) * dnym1
         ii = 0
         do i =cell_low(1,c), cell_high(1,c)
            xi = dble(i) * dnxm1
            call exact_solution(xi, eta, zeta, temp)
            do m = 1, 5
               u(m,ii,jj,kk,c) = temp(m)
            enddo
            ii = ii + 1
         enddo
         jj = jj + 1
      enddo

      return
      end


c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine lhsinit

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      include 'header.h'
      
      integer i, j, k, d, c, m, n

c---------------------------------------------------------------------
c     loop over all cells                                       
c---------------------------------------------------------------------
      do c = 1, ncells

c---------------------------------------------------------------------
c     first, initialize the start and end arrays
c---------------------------------------------------------------------
         do d = 1, 3
            if (cell_coord(d,c) .eq. 1) then
               start(d,c) = 1
            else 
               start(d,c) = 0
            endif
            if (cell_coord(d,c) .eq. ncells) then
               end(d,c) = 1
            else
               end(d,c) = 0
            endif
         enddo

c---------------------------------------------------------------------
c     zero the whole left hand side for starters
c---------------------------------------------------------------------
         do k = 0, cell_size(3,c)-1
            do j = 0, cell_size(2,c)-1
               do i = 0, cell_size(1,c)-1
                  do m = 1,5
                     do n = 1, 5
                        lhsc(m,n,i,j,k,c) = 0.0d0
                     enddo
                  enddo
               enddo
            enddo
         enddo

      enddo

      return
      end


c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine lhsabinit(lhsa, lhsb, size)
      implicit none

      integer size
      double precision lhsa(5, 5, -1:size), lhsb(5, 5, -1:size)

      integer i, m, n

c---------------------------------------------------------------------
c     next, set all diagonal values to 1. This is overkill, but convenient
c---------------------------------------------------------------------
      do i = 0, size
         do m = 1, 5
            do n = 1, 5
               lhsa(m,n,i) = 0.0d0
               lhsb(m,n,i) = 0.0d0
            enddo
            lhsb(m,m,i) = 1.0d0
         enddo
      enddo

      return
      end