vtf-logo

fsi/sfc-amroc/TubeCJBurnFlaps/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)
      include  "cuser.i"
c
      parameter( lin = 5, lmechout = 11, lrin = 13  )
      character *16 cwork
      dimension cwork(2)
      character *256 filename
c
      cwork(1)= 'Product'
      cwork(2)= 'Reactant'
c      
      open(unit=lin,status='old',form='formatted',file='init.dat')
      read (lin, *) gamma, qr
      read (lin, *) ar, pr, ao, po
      read (lin, *) sloc, Nden, NCJ
      read (lin, *) Wk(1), Wk(2), RU, PA
      read (lin, *) Nread, filename, sshift
      read (lin, *) rf, zt0, zt1
      read (lin, *) rfback, Nback
      close (lin)
c
      open(unit=lmechout, status='unknown', form='formatted', 
     &     file='chem.dat')
      write (lmechout,400) RU
      write (lmechout,401) PA
      write (lmechout,402) (k,cwork(k),k=1,2)        
      write (lmechout,403) (k,Wk(k),k=1,2)        
      close (lmechout)
c
 400  format('RU       ',e16.8)
 401  format('PA       ',e16.8)
 402  format('Sp(',i2.2,')     ',a16)
 403  format('W(',i2.2,')    ',e16.8)
c
      gamma1 = gamma-1.d0
c
      Yo = 0.d0
      uo = 0.d0
      if (Nden.eq.1) then
         rhoo = ao
      else
         To   = ao
         Wo = 1.d0/((1.d0-Yo)/Wk(1) + Yo/Wk(2))
         rhoo = (po*Wo)/(RU*To)
      end if
c
      Yr = 1.d0
      ur = 0.d0
      if (Nden.eq.1) then
         rhor = ar
      else
         Tr   = ar
         Wr = 1.d0/((1.d0-Yr)/Wk(1) + Yr/Wk(2))
         rhor = (pr*Wr)/(RU*Tr)
      end if
      V0 = 1.d0/rhor
      P0 = pr
c
      if (NCJ.eq.1) then
         Dj = qr
         Dj = Dj/dsqrt(P0*V0)
         qn = 0.5d0*((Dj+gamma/Dj)**2-4.d0*gamma)/
     &        (gamma**2-1.d0)
         q0 = qn*P0*V0
      else
         q0 = qr
         qn=q0/(P0*V0)
         Dj=dsqrt((gamma**2-1.d0)*qn/2.d0)
     &        +dsqrt((gamma**2-1.d0)*qn/2.d0+gamma)
      end if
      Vj=gamma*(1.d0+Dj**2)/((gamma+1.d0)*Dj**2)
      Pj=(1.d0+Dj**2)/(gamma+1.d0)
      Uj=-Dj*Vj+Dj
c
      write (6,*) 'Normalized CJ state:',
     &     Dj,Pj,1.d0/Vj,Uj,qn
c
      Dj = Dj*dsqrt(P0*V0)
      Pj = Pj*P0
      Vj = Vj*V0
      Uj = Uj*dsqrt(P0*V0)
c
      write (6,*) 'Actual CJ state:',
     &     Dj,Pj,1.d0/Vj,Uj,q0
c
      rhol = 1.d0/Vj
      ul = Uj
      pl = Pj
      Yl = 0.d0
c
      Yact = 0.01d0
c
      if (Nread.gt.0) then
         do k=1,256
            if (filename(k:k).eq.' ') goto 100
         enddo
         Nread = Nread + 2
 100     open(unit=lrin, status='old', form='formatted', 
     &        file=filename(1:k-1))
         do i=2,Nread-1
            read (lrin,*) xtab(i),(qtab(k,i),k=1,4)
            xtab(i) = xtab(i)-sshift
         enddo
         close (lrin)
         xtab(1) = xtab(2) - 0.5d0*(xtab(3)-xtab(2))
         alpha = (xtab(1)-xtab(2))/(xtab(3)-xtab(2))
         do k=1,4
            qtab(k,1) = (1.d0-alpha)*qtab(k,2) + alpha*qtab(k,3)
         enddo
c
         xtab(Nread) = xtab(Nread-1)+0.5d0*(xtab(Nread-1)-xtab(Nread-2))
         alpha = (xtab(Nread)-xtab(Nread-1))/
     &        (xtab(Nread-2)-xtab(Nread-1))
         do k=1,4
            qtab(k,Nread) = (1.d0-alpha)*qtab(k,Nread-1) + 
     &           alpha*qtab(k,Nread-2)
         enddo
      endif
c
      return
      end