*********************************************************************   
      subroutine shtube  
*   
*      Theoretical shock tube results   
*   
*      called by model (shock)  
*   
      include 'cparam.inc'
      include 'cfiles.inc'
      include 'cdata.inc'
      include 'cwork.inc'
*   
      f(p,p0) = (p-1.)**2/(4.*p+1) - 5.*(1.-(p/p0)**0.2)**2 
      df(p,p0) = 2.*(p-1.)*(2*p+3.)/(4.*p+1)**2 
     1    + (2./p0)*(1.0-(p/p0)**0.2)*(p0/p)**0.8   
*
	print *,'shtube0, 1.1, 10-may-90'
*   
*  Physical constants
*
      boltz=1.38062e-16
      amass=1.660531e-24
      uue= (uul/uut)**2
      rgas=boltz/amass/uue
      xmu=1.2
      gamma=5./3.
*
      p0 = pp(1,nz)/pp(1,1)
      p = 1.1   
      do 10 i=1,20  
      fp = f(p,p0)  
      dfp = df(p,p0)
      dfp1 = fp/dfp 
      if(abs(dfp1).lt.1.e-5) go to 2
      p = p - dfp1  
   10 continue  
  701 format(i10,1p5e12.2)  
*      write (*,711) dfp1,p,p0 
  711 format(1x,35hno convergence for shock tube: dp =,1pe11.2,2e13.2)  
*   
    2 s0 = sqrt(gamma*pp(1,1)/rho(1,1)) 
      rp5 = (p/p0)**0.2 
      vcd = 3.0*s0*(1.0-rp5)
      vsh = vcd*((8./3.)*p+(2./3.))/(2.*(p-1.)) 
      dnp = rho(1,nz)
      dnm = rho(1,1)  
      dnl = dnp*(p*p0**(2./3.))**(3./5.)
      dnr = dnp*((5./3.)*(p+1.)+(p-1.))/((5./3.)*(p+1.)-(p-1.)) 
      tgp = tt(1,nz)  
      tgm = tt(1,1)  
      tgr = tgp*p*(dnp/dnr) 
      tgl = tgp*p*(dnp/dnl) 
	eep = ee(1,nz)
	eem1 = ee(1,1)
      eer = eep*p*(dnp/dnr) 
      eel = eep*p*(dnp/dnl) 
      vrp = vcd - s0*rp5
      vrm = - s0
      ppm1=pp(1,1)
      ppp=pp(1,nz)
      dnr=dnm*((gamma+1.)*p/(gamma-1.)+1.)/((gamma+1.)/(gamma-1.)+p)
	dnl=p*ppm1/eel/(gamma-1.)
c   
      write (*,702) vsh,vcd,vrp,vrm,dnm,dnr,dnl,dnp,eem1,eer,eel,eep 
  702 format (' shock tube results: vsh =',1pe11.3,'  vcd =',   
     1    e11.3,'  vrp =',e11.3,'  vrm =',e11.3/5x,'dnp =',e11.3, 
     2    '  dnr =',e11.3,'  dnl =',e11.3,'  dnm =',e11.3
     3    /5x,'eep =',e11.3,'  eer =',e11.3,'  eel =',e11.3,   
     4    '  eem =',e11.3)
c   
*  Store in spares in cdata
*
      scra(5)=vsh
      scra(6)=vcd
      scra(7)=vrp
      scra(8)=vrm
	scra(1)=dnm
	scra(2)=dnr
	scra(3)=dnl
	scra(4)=dnp
	scrb(1)=eem1
	scrb(2)=eer
	scrb(3)=eel
	scrb(4)=eep
*
*     zcd0 = z(n0)  
c      if(itime.ne.0) go to 9   
c      dtn = dtn/p0 
c      dtnm = dtn   
c   
    9 return
      end   
