vtf-logo

weno/applications/euler/3d/OneSphere/src/init3.f

c
c     =====================================================
      subroutine ic(maxmx,maxmy,maxmz,meqn,mbc,mx,my,mz,
     &              x,y,z,dx,dy,dz,q)
c     =====================================================
c     
c     Copyright (C) 2003-2007 California Institute of Technology
c     Ralf Deiterding, ralf@amroc.net
c
      implicit double precision (a-h,o-z)
c     
      include 'air.i'
      include 'cuser.i'
c     
      dimension q(meqn, 1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc, 
     &            1-mbc:maxmz+mbc)
      dimension x(1-mbc:maxmx+mbc),y(1-mbc:maxmy+mbc),
     &          z(1-mbc:maxmz+mbc)
c
      do 60 k = 1, mz
         do 60 j = 1, my
            do 60 i = 1, mx
               call cellave(x(i)-dx/2.d0,y(j)-dy/2.d0,dx,dy,wl)
               rho =    (1.d0-wl)*rhoamb+ wl*rhoshk
               u =      (1.d0-wl)*uamb  + wl*ushk
               v =      (1.d0-wl)*vamb  + wl*vshk
               w =      (1.d0-wl)*wamb  + wl*wshk
               p =      (1.d0-wl)*pamb  + wl*pshk
               q(1,i,j,k) = rho
               q(2,i,j,k) = rho * uamb
               q(3,i,j,k) = rho * v
               q(4,i,j,k) = rho * w
               q(5,i,j,k) = p/gamma1 + .5d0*rho*(u*u+v*v+w*w)
               q(6,i,j,k) = (1.d0-wl)*s1amb+ wl*s1shk
               q(7,i,j,k) = (1.d0-wl)*s2amb+ wl*s2shk
               q(8,i,j,k) = (1.d0-wl)*dcamb+ wl*dcshk
 60   continue
         
      return
      end