c
c
c =====================================================
subroutine ic(maxmx,maxmy,maxmz,meqn,mbc,mx,my,mz,
& x,y,z,dx,dy,dz,q)
c =====================================================
implicit none
c
include "cuser.i"
c
integer maxmx, maxmy, maxmz, meqn, mbc, mx, my, mz
double precision q(meqn, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc,
& 1-mbc:maxmz+mbc)
double precision x(1-mbc:maxmx+mbc),y(1-mbc:maxmy+mbc),
& z(1-mbc:maxmz+mbc), dx, dy, dz
double precision uq(32,32,32,5)
integer i, j, k, l
integer il, jl, kl, ip, jp, kp, im, jm, km
double precision dxsrc, eta1, eta2, eta3
character*(128) name
c
c Homogeneous Turbulence test case
do 50 k = 1-mbc,maxmz+mbc
do 50 j = 1-mbc,maxmy+mbc
do 50 i = 1-mbc,maxmx+mbc
! --- this insures the any extended bits
! --- in the vector of state will start as zero
do l = 1,meqn
q(l,i,j,k) = 0.0d0
end do
50 continue
name = 'INIT_DATA/M0.488_Re375_t0.83_32.bin'//CHAR(0)
call load_binary_file(name, 32, 32, 32, uq)
dxsrc = 6.283185307d0/32
do 60 k = 1, mz
kl=z(k)/dxsrc+1
kp = kl+1
km = kl-1
if ( kp .gt. 32 ) kp = 1
if ( km .lt. 1 ) km = 32
eta3 = 0.5d0*(z(k)-dxsrc*(kl-0.5d0))
do 60 j = 1, my
jl=y(j)/dxsrc+1
jp = jl+1
jm = jl-1
if ( jp .gt. 32 ) jp = 1
if ( jm .lt. 1 ) jm = 32
eta2 = 0.5d0*(y(j)-dxsrc*(jl-0.5d0))
do 60 i = 1, mx
il=x(i)/dxsrc+1
ip = il+1
im = il-1
if ( ip .gt. 32 ) ip = 1
if ( im .lt. 1 ) im = 32
eta1 = 0.5d0*(x(i)-dxsrc*(il-0.5d0))
do l=1,5
q(l,i,j,k) = uq(il,jl,kl,l)
$ +eta1*(uq(ip,jl,kl,l)-uq(im,jl,kl,l))
$ +eta2*(uq(il,jp,kl,l)-uq(il,jm,kl,l))
$ +eta3*(uq(il,jl,kp,l)-uq(il,jl,km,l))
enddo
60 continue
return
end