vtf-logo

clawpack/applications/euler_chem/3d/ConvReflDet/src/combl.f

c     =====================================================
      subroutine combl()
c     =====================================================
c
c     Create and initialize application specific common-blocks.
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  "ck.i"
      include  "cuser.i" 
c
      dimension Xk(leNsp), Yk(leNsp)
c
      call cheminit()      
c     
c     SI-units
c     
      rhol = 0.072d0 
      ul   = 0.d0
      vl   = 0.d0
      wl   = 0.d0
      pl   = 7123.d0
c     
      rhor = 0.18075d0 
      ur   = 0.d0
      vr   = 0.d0
      wr   = 0.d0
      pr   = 35594.d0
c     
c     Convert kg/m**3 into g/cm**3
      rhol = rhol * 1.d-3
      rhor = rhor * 1.d-3
c     Convert Pa = J/m**2 into dynes/cm**2
      pl = pl * 1.d1
      pr = pr * 1.d1
c     Convert m/sec into cm/sec
      ul = ul * 1.d2
      ur = ur * 1.d2
c     Convert m/sec into cm/sec
      vl = vl * 1.d2
      vr = vr * 1.d2
      u  = u  * 1.d2
c
      do k=1,Nsp
         Xk(k) = 0.d0
      enddo
      Xk(1) = 1.d0
      Xk(6) = 2.d0
      Xk(9) = 7.d0
c
      sum = 0.d0
      do k = 1, Nsp
         sum = sum + Xk(k)
      enddo
      do k=1,Nsp
         Xk(k) = Xk(k)/sum
      enddo
c 
      call ckxty (Xk,iwork,rwork,Yk)
      call ckmmwy (Yk,iwork,rwork,W)
c
c     # data in left state:
      Tl = (pl*W)/(rhol*RU)
      rhoul = rhol * ul
      rhovl = rhol * vl
      rhowl = rhol * wl
      el = rhol * (0.5d0*(ul**2+vl**2+wl**2) + 
     &     avgtabip(Tl,Yk,hms,Nsp)) - pl
c
      do k=1, Nsp
         qin(k) = rhol*Yk(k)
      enddo
      qin(Nsp+1) = rhoul 
      qin(Nsp+2) = rhovl 
      qin(Nsp+3) = rhowl
      qin(Nsp+4) = el
      qin(Nsp+5) = Tl
c
      do k=1,Nsp 
         Xk(k) = 0.d0
      enddo
      Xk(2) = 1.d0
      Xk(3) = 2.d0
      Xk(9) = 7.d0
c
      sum = 0.d0
      do k = 1, Nsp
         sum = sum + Xk(k)
      enddo
      do k=1,Nsp
         Xk(k) = Xk(k)/sum
      enddo
c 
      call ckxty (Xk,iwork,rwork,Yk)
      call ckmmwy (Yk,iwork,rwork,W)
c
c     # data in right state:
      Tr = (pr*W)/(rhor*RU)
      rhour = rhor * ur
      rhovr = rhor * vr
      rhowr = rhor * wr
      er = rhor * (0.5d0*(ur**2+vr**2+wr**2) + 
     &     avgtabip(Tr,Yk,hms,Nsp)) - pr
c
      write(6,*) 'left state'
      write(6,*) rhol, ul, vl, wl, pl, Tl
      write(6,*) 'right state'
      write(6,*) rhor, ur, vr, wr, pr, Tr
c     
      x0 = 0.d0
      y0 = 0.d0
      z0 = 0.d0 
      ra = 6.0d0
c
      do k=1, Nsp
         qout(k) = rhor*Yk(k)
      enddo
      qout(Nsp+1) = rhour
      qout(Nsp+2) = rhovr
      qout(Nsp+3) = rhowr
      qout(Nsp+4) = er
      qout(Nsp+5) = Tr
c
      pi = 4.d0*atan(1.d0)
      idisc = 2
c
c     Geometry for reclecting boundary conditions
c
      angpl = pi/4.d0
      rdi = 0.5d0
c
c     Outer radius - extrapolation
c
      rd = 11.5d0
c
      return
      end