!---------------------------------------------------------------------
!---------------------------------------------------------------------
!---------------------------------------------------------------------
!     This subroutine initializes the field variable u using 
!     tri-linear transfinite interpolation of the boundary values     
!---------------------------------------------------------------------
      subroutine initialize ()

      include 'header3d.h'
      integer  i,j,k,m,ix,iy,iz
      double precision  xi,eta,zeta,pface(5,3,2),pxi,peta,pzeta,temp(5),
     &xi1,yi1,zi1
      xi = 0.0
      eta = 0.0
      zeta = 0.0

!---------------------------------------------------------------------
!  Later (in compute_rhs) we compute 1/u for every element. A few of 
!  the corner elements are not used, but it convenient (and faster) 
!  to compute the whole thing with a simple loop. Make sure those 
!  values are nonzero by initializing the whole thing here. 
!---------------------------------------------------------------------
!DVM$ region out(u)
!DVM$ PARALLEL (k,j,i)  ON  u(*,i,j,k), SHADOW_COMPUTE, private(m)
      do  k = 0,imax - 1
         do  j = 0,imax - 1
            do  i = 0,imax - 1
               do  m = 1,5
                  u(m,i,j,k) = 1.0
               enddo  
            enddo  
         enddo  
      enddo  

!DVM$ PARALLEL (k,j,i)  ON  u(*,i,j,k), private(m,zeta,eta, xi,ix,
!DVM$& iy,iz,Pxi,Peta,Pzeta,Pface,xi1,yi1,zi1,temp),SHADOW_COMPUTE
      do  k = 0,grid_points(3) - 1
         do  j = 0,grid_points(2) - 1
            do  i = 0,grid_points(1) - 1
               zeta = dble (k) * dnzm1
               eta = dble (j) * dnym1
               xi = dble (i) * dnxm1
               do  ix = 1,2

!                  call exact_solution(dble(ix-1), eta, zeta, Pface(1,1,ix))
                  xi1 = dble (ix - 1)
                  do  m = 1,5
                     pface(m,1,ix) = ce(m,1) + xi1 * (ce(m,2) + xi1 * (c
     &e(m,5) + xi1 * (ce(m,8) + xi1 * ce(m,11)))) + eta * (ce(m,3) + eta
     & * (ce(m,6) + eta * (ce(m,9) + eta * ce(m,12)))) + zeta * (ce(m,4)
     & + zeta * (ce(m,7) + zeta * (ce(m,10) + zeta * ce(m,13))))
                  enddo  
               enddo  
               do  iy = 1,2

!                  call exact_solution(xi, dble(iy-1) , zeta, Pface(1,2,iy))
                  yi1 = dble (iy - 1)
                  do  m = 1,5
                     pface(m,2,iy) = ce(m,1) + xi * (ce(m,2) + xi * (ce(
     &m,5) + xi * (ce(m,8) + xi * ce(m,11)))) + yi1 * (ce(m,3) + yi1 * (
     &ce(m,6) + yi1 * (ce(m,9) + yi1 * ce(m,12)))) + zeta * (ce(m,4) + z
     &eta * (ce(m,7) + zeta * (ce(m,10) + zeta * ce(m,13))))
                  enddo  
               enddo  
               do  iz = 1,2

!                  call exact_solution(xi, eta, dble(iz-1), Pface(1,3,iz))
                  zi1 = dble (iz - 1)
                  do  m = 1,5
                     pface(m,3,iz) = ce(m,1) + xi * (ce(m,2) + xi * (ce(
     &m,5) + xi * (ce(m,8) + xi * ce(m,11)))) + eta * (ce(m,3) + eta * (
     &ce(m,6) + eta * (ce(m,9) + eta * ce(m,12)))) + zi1 * (ce(m,4) + zi
     &1 * (ce(m,7) + zi1 * (ce(m,10) + zi1 * ce(m,13))))
                  enddo  
               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,i,j,k) = pxi + peta + pzeta - pxi * peta - pxi * p
     &zeta - peta * pzeta + pxi * peta * pzeta
               enddo  
			   
			   if(i .eq. 0) then  
	               do  m = 1,5
                  temp(m) = ce(m,1) + xi * (ce(m,2) + xi * (ce(m,5) + xi
     & * (ce(m,8) + xi * ce(m,11)))) + eta * (ce(m,3) + eta * (ce(m,6) +
     & eta * (ce(m,9) + eta * ce(m,12)))) + zeta * (ce(m,4) + zeta * (ce
     &(m,7) + zeta * (ce(m,10) + zeta * ce(m,13))))
               enddo  
               do  m = 1,5
                  u(m,i,j,k) = temp(m)
               enddo 
               endif
			   if(i .eq. grid_points(1) - 1) then
			   xi = 1.0d0
               do  m = 1,5
                  temp(m) = ce(m,1) + xi * (ce(m,2) + xi * (ce(m,5) + xi
     & * (ce(m,8) + xi * ce(m,11)))) + eta * (ce(m,3) + eta * (ce(m,6) +
     & eta * (ce(m,9) + eta * ce(m,12)))) + zeta * (ce(m,4) + zeta * (ce
     &(m,7) + zeta * (ce(m,10) + zeta * ce(m,13))))
               enddo  
               do  m = 1,5
                  u(m,i,j,k) = temp(m)
               enddo  			   
			   endif
			   
			   if(j .eq. 0) then
			   zeta = dble (k) * dnzm1
               xi = dble (i) * dnxm1
			   eta = 0.0d0
			   
			   do  m = 1,5
                  temp(m) = ce(m,1) + xi * (ce(m,2) + xi * (ce(m,5) + xi
     & * (ce(m,8) + xi * ce(m,11)))) + eta * (ce(m,3) + eta * (ce(m,6) +
     & eta * (ce(m,9) + eta * ce(m,12)))) + zeta * (ce(m,4) + zeta * (ce
     &(m,7) + zeta * (ce(m,10) + zeta * ce(m,13))))
               enddo  
               do  m = 1,5
                  u(m,i,j,k) = temp(m)
               enddo  
			   endif
			   
               if(j .eq. grid_points(2) - 1) then
               zeta = dble (k) * dnzm1
               xi = dble (i) * dnxm1
               eta = 1.0d0
!            call exact_solution(xi, eta, zeta, temp)
               do  m = 1,5
                  temp(m) = ce(m,1) + xi * (ce(m,2) + xi * (ce(m,5) + xi
     & * (ce(m,8) + xi * ce(m,11)))) + eta * (ce(m,3) + eta * (ce(m,6) +
     & eta * (ce(m,9) + eta * ce(m,12)))) + zeta * (ce(m,4) + zeta * (ce
     &(m,7) + zeta * (ce(m,10) + zeta * ce(m,13))))
               enddo  
               do  m = 1,5
                  u(m,i,j,k) = temp(m)
               enddo 			   
			   endif
			   
			   if(k .eq. 0) then
               zeta = 0.0d0
               xi = dble (i) * dnxm1
               eta = dble (j) * dnym1

!            call exact_solution(xi, eta, zeta, temp)
               do  m = 1,5
                  temp(m) = ce(m,1) + xi * (ce(m,2) + xi * (ce(m,5) + xi
     & * (ce(m,8) + xi * ce(m,11)))) + eta * (ce(m,3) + eta * (ce(m,6) +
     & eta * (ce(m,9) + eta * ce(m,12)))) + zeta * (ce(m,4) + zeta * (ce
     &(m,7) + zeta * (ce(m,10) + zeta * ce(m,13))))
               enddo  
               do  m = 1,5
                  u(m,i,j,k) = temp(m)
               enddo                			   
			   endif
			   
			   if(k .eq. grid_points(3) - 1) then
               zeta = 1.0d0
               xi = dble (i) * dnxm1
               eta = dble (j) * dnym1

!            call exact_solution(xi, eta, zeta, temp)
               do  m = 1,5
                  temp(m) = ce(m,1) + xi * (ce(m,2) + xi * (ce(m,5) + xi
     & * (ce(m,8) + xi * ce(m,11)))) + eta * (ce(m,3) + eta * (ce(m,6) +
     & eta * (ce(m,9) + eta * ce(m,12)))) + zeta * (ce(m,4) + zeta * (ce
     &(m,7) + zeta * (ce(m,10) + zeta * ce(m,13))))
               enddo  
               do  m = 1,5
                  u(m,i,j,k) = temp(m)
               enddo 			   
			   endif
            enddo  
         enddo  
      enddo 
!DVM$ end region

      return
      end
