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 i, j, k, m, ix, iy, iz
       double precision  xi, eta, zeta, Pface(5,3,2), Pxi, Peta, 
     >                   Pzeta, temp(5)

!DVM$ region
!DVM$ parallel (k,j,i) on u(*,i,j,k), private(zeta, eta, xi, ix, pxi, m,
!DVM$& pface, iy, peta, iz, pzeta, temp)
      do k = 0, problem_size-1
         do j = 0, problem_size-1
            do i = 0, problem_size-1
               u(1,i,j,k) = 1.0
               u(2,i,j,k) = 0.0
               u(3,i,j,k) = 0.0
               u(4,i,j,k) = 0.0
               u(5,i,j,k) = 1.0
			   
               zeta = dble(k) * dnzm1
               eta = dble(j) * dnym1
               xi = dble(i) * dnxm1		
			   
               do ix = 1, 2
                 Pxi = dble(ix-1)

                 do  m = 1, 5
                   Pface(m,1,ix) =  ce(m,1) +
     >             Pxi*(ce(m,2)+Pxi*(ce(m,5) +Pxi*(ce(m,8) +
     >             Pxi*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))))
                 end do
              end do

              do iy = 1, 2
                 Peta = 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)))) +
     >             Peta*(ce(m,3) +Peta*(ce(m,6) +Peta*(ce(m,9)+
     >             Peta*ce(m,12))))+
     >             zeta*(ce(m,4) + zeta*(ce(m,7) + zeta*(ce(m,10) + 
     >             zeta*ce(m,13))))
                end do
              end do

              do iz = 1, 2
                Pzeta = 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))))+
     >           Pzeta*(ce(m,4) + Pzeta*(ce(m,7) + Pzeta*(ce(m,10) + 
     >           Pzeta*ce(m,13))))
                end do
             end do

             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*Pzeta - Peta*Pzeta + 
     >             Pxi*Peta*Pzeta	
             end do	 
			 
             zeta = dble(k) * dnzm1
             eta = dble(j) * dnym1
             xi = 0.0d0			 
             if( i .eq. 0) then
!             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))))
              end do
              do  m = 1, 5
                u(m,i,j,k) = temp(m)
              end do
	         endif
			 
             xi = 1.0d0	
             if( i .eq. problem_size-1) then
!             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))))
              end do
              do  m = 1, 5
                u(m,i,j,k) = temp(m)
              end do
	         endif			 

             zeta = dble(k) * dnzm1
             eta = 0.0d0
             xi = dble(i) * dnxm1
	         if( j .eq. 0) then
!             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))))
               end do
               do  m = 1, 5
                 u(m,i,j,k) = temp(m)
               end do
             endif
			 
             eta = 1.0d0			 
	         if( j .eq. problem_size-1) then
!             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))))
               end do
               do  m = 1, 5
                 u(m,i,j,k) = temp(m)
               end do
             endif	
			 
             zeta = 0.0d0
             eta = dble(j) * dnym1
             xi = dble(i) *dnxm1			 
	         if( k .eq. 0) then
!             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))))
               end do
               do  m = 1, 5
                 u(m,i,j,k) = temp(m)
               end do
             endif	
			 
             zeta = 1.0d0	
	         if( k .eq. problem_size-1) then
!             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))))
               end do
               do  m = 1, 5
                 u(m,i,j,k) = temp(m)
               end do
             endif	
           end do
        end do
      end do
!DVM$ end region

       return
       end