vtf-logo

clawpack/applications/euler_chem/1d/RealDetonation/src/ownout1eurhok.f

c
c     ==========================================================
      subroutine out1eurhok(q,mx,lb,ub,qo,mxo,lbo,ubo,
     &     lbr,ubr,shaper,meqn,nc,time)
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)
      include  "ck.i"    
c
      integer meqn, mx, mxo
      dimension q(meqn,mx), qo(mxo)
c
      integer  lb(1), ub(1), lbo(1), ubo(1), lbr(1), ubr(1), shaper(1), 
     &     mresult, stride, imin(1), imax(1), i, getindx
c
      stride = (ub(1) - lb(1))/(mx-1)
      imin(1) = max(lb(1), lbo(1), lbr(1))
      imax(1) = min(ub(1), ubo(1), ubr(1))
c
      if (mod(imin(1)-lb(1),stride) .ne. 0) then
         imin(1) = imin(1) + stride - mod(imin(1)-lb(1),stride) 
      endif
      imin(1) = getindx(imin(1), lb(1), stride)  

      if (mod(imax(1)-lb(1),stride) .ne. 0) then
         imax(1) = imax(1) - mod(imax(1)-lb(1),stride) 
      endif
      imax(1) = getindx(imax(1), lb(1), stride)  
c
      do 10 i = imin(1), imax(1)
         rho  = 0.d0
         rhoW = 0.d0
         do k = 1, Nsp
            rho  = rho  + q(k,i)
            rhoW = rhoW + q(k,i)/Wk(k)
         enddo
c
c        # Total density 
         if (nc.eq.1) then 
            qo(i) = rho
c           # Convert g/cm**3 into kg/m**3
            if (ckunits) qo(i) = qo(i)*1.d3
         endif
c
c        # Velocity
         if (nc.eq.2) then
            qo(i) = q(Nsp+1,i)/rho
c           # Convert cm/sec into m/sec
            if (ckunits) qo(i) = qo(i)*1.d-2
         endif
c
c        # Total energy density
         if (nc.eq.3) then
            qo(i) = q(Nsp+2,i)
c           # Convert ergs/cm**3 into J/m**3
            if (ckunits) qo(i) = qo(i)*1.d-1
         endif
c
c        # Temperature - No unit conversion
         if (nc.eq.4) qo(i) = q(Nsp+3,i)
c
c        # Pressure
         if (nc.eq.5) then
            qo(i) = rhoW*RU*q(Nsp+3,i)
c           # Convert dynes/cm**2 into Pa = J/m**2
            if (ckunits) qo(i) = qo(i)*1.d-1
         endif
c
c        # Gamma - No unit conversion
         if (nc.eq.6) then
            rhoCp = avgtabip(q(Nsp+3,i),q(1,i),cpk,Nsp)
            qo(i) = RU/(rhoCp/rhoW-RU)+1.d0       
         endif
c
c        # Mass fractions - No unit conversion
         if (nc.ge.7.and.nc.le.Nsp+7) qo(i) = q(nc-6,i)/rho
c
c        # Number of succesful time steps in ODE-solver in last step
         if (nc.eq.Nsp+8) qo(i) = q(Nsp+5,i)
c
c        # Number of unsuccesful time steps in ODE-solver in last step
         if (nc.eq.Nsp+9) qo(i) = q(Nsp+6,i)
c         
 10   continue         

      return
      end