• VTF
  • FSI
  • AMROC
  • SFC
  • Motion
  • STLIB
  • Main Page
  • src/2d/equations/euler/rpznd/flgout2euznd.f

    c
    c     ==========================================================
          subroutine flgout2euznd(q,mx,my,lb,ub,qo,mxo,myo,lbo,ubo,
         &     lbr,ubr,shaper,meqn,nc,t)
    c     ==========================================================
    c
    c     # Computes primitives for ZND Euler equations for output 
    c     # and flagging.
    c
    c     # Copyright (C) 2002 Ralf Deiterding
    c     # Brandenburgische Universitaet Cottbus
    c
    c     # Copyright (C) 2003-2007 California Institute of Technology
    c     # Ralf Deiterding, ralf@cacr.caltech.edu
    c
          implicit double precision(a-h,o-z)
          common /param/  gamma,gamma1,q0
          common /PhysData/  Wk(2), RU, PA
    c
          integer meqn, mx, my, mxo, myo
          dimension q(meqn,mx,my), qo(mxo,myo)
    c
          integer  lb(2), ub(2), lbo(2), ubo(2), lbr(2), ubr(2), shaper(2), 
         &     mresult, stride, imin(2), imax(2), i, getindx, d
    c
          stride = (ub(1) - lb(1))/(mx-1)
          do 5 d = 1, 2
             imin(d) = max(lb(d), lbr(d))
             imax(d) = min(ub(d), ubr(d))
    
             if (mod(imin(d)-lb(d),stride) .ne. 0) then
                imin(d) = imin(d) + stride - mod(imin(d)-lb(d),stride) 
             endif
             imin(d) = getindx(imin(d), lb(d), stride)  
    
             if (mod(imax(d)-lb(d),stride) .ne. 0) then
                imax(d) = imax(d) - mod(imax(d)-lb(d),stride) 
             endif
             imax(d) = getindx(imax(d), lb(d), stride)  
     5    continue
    
          do 10 i = imin(1), imax(1)
             do 10 j = imin(2), imax(2)
    c           # Density
                if (nc.eq.1) qo(i,j) = q(1,i,j) + q(2,i,j) 
    c           # Velocity u
                if (nc.eq.2) qo(i,j) = q(3,i,j)/(q(1,i,j) + q(2,i,j))
    c           # Velocity v
                if (nc.eq.3) qo(i,j) = q(4,i,j)/(q(1,i,j) + q(2,i,j))
    c           # Total energy density
                if (nc.eq.4) qo(i,j) = q(5,i,j)
    c           # Temperature 
                if (nc.eq.5) then
                   rho = q(1,i,j) + q(2,i,j)
                   p = gamma1*(q(5,i,j) - q(2,i,j)*q0 -
         &              0.5d0*(q(3,i,j)**2+q(4,i,j)**2)/rho)
                   W = 1.d0/((q(1,i,j)/rho)/Wk(1) + (q(2,i,j)/rho)/Wk(2))
                   qo(i,j) = (p*W)/(rho*RU)
                endif
    c           # Pressure
                if (nc.eq.6) qo(i,j) = gamma1*(q(5,i,j) - q(2,i,j)*q0 -
         &           0.5d0*(q(3,i,j)**2+q(4,i,j)**2)/(q(1,i,j) + q(2,i,j)))
    c           # Gamma 
                if (nc.eq.7) qo(i,j) = gamma
    c           # Y1 
                if (nc.eq.8) qo(i,j) = q(1,i,j) / (q(1,i,j) + q(2,i,j))
    c           # Y2 
                if (nc.eq.9) qo(i,j) = q(2,i,j) / (q(1,i,j) + q(2,i,j))
     10   continue         
    
          return
          end
    

<