applications/euler/2d/Shockbubble/src/init2.f

c
c
c     =====================================================
       subroutine ic(maxmx,maxmy,meqn,mbc,mx,my,x,y,dx,dy,q)
c     =====================================================
       implicit double precision (a-h,o-z)
c
       include  "cuser.i"
c
       dimension q(meqn,1-mbc:maxmx+mbc, 1-mbc:maxmy+mbc)
       dimension x(1-mbc:maxmx+mbc),y(1-mbc:maxmy+mbc)
c
c
       pi   = 4.d0*datan(1.d0)
       pi2  = 2.d0*pi
       disp = 0.0d0
c
c     # circle of concentration 1 against background concentration 0:
       do 60 i = 1, mx

          if(x(i).ge.disp) then
             
             idisc = 2
             r0    = 0.2d0
             x0    = .4d0
             y0    = .0d0
             
             rhol  = 0.1d0
             rhor  = 1.0d0
             pl    = 1.d0
             pr    = 1.d0
             ur    = 0.d0
             ul    = 0.d0
             vr    = 0.d0
             vl    = 0.d0
             
          else
             idisc = 1
             x0    = 0.0d0
             y0    = 0.0d0
             alf   = 1.0d0
             beta  = 0.0d0   
          
             rhor  = 1.0d0
             rhol  = 3.81d0
             pl    = 10.d0            
             pr    = 1.d0
             ul    = 2.58d0
             ur    = 0.d0
             vr    = 0.d0
             vl    = 0.d0
          end if

          do 60 j = 1, my
             
             call cellave(x(i)-dx/2.d0,y(j)-dy/2.d0,dx,dy,wl)
             q(1,i,j) = wl*rhol + (1.d0-wl)*rhor
             q(2,i,j) = q(1,i,j) * (wl*ul + (1.d0-wl)*ur)
             q(3,i,j) = q(1,i,j) * (wl*vl + (1.d0-wl)*vr)
             q(4,i,j) = wl*(pl/gamma1) + (1.d0-wl)*(pr/gamma1) +
     $            .5d0*(q(2,i,j)*q(2,i,j) + 
     $            q(3,i,j)*q(3,i,j)) / q(1,i,j)  
 60       continue
         
       return
       end