vtf-logo

clawpack/applications/euler/2d/VortexRotGFM/src/ip2navrfl.f

c     -----------------------------------------------------
c     Internal reflecting physical boundary conditions
c     for Navier-Stokes equations
c     -----------------------------------------------------
c
c     Copyright (C) 2003-2007 California Institute of Technology
c     Ralf Deiterding, ralf@amroc.net
c
c     Transformation of vector of conserved quantities
c     into primitives (rho,u,v,p)
c
c     =====================================================
      subroutine it2eurfl(mx,my,meqn,q,qt)
c     =====================================================
      implicit none
c     
      common /param/  gamma,gamma1
      double precision    gamma,gamma1
      integer   i, j, mx, my, meqn
      double precision    q(meqn,mx,my), qt(meqn,mx,my)
c
      do 10 j = 1, my
         do 10 i = 1, mx 
            qt(1,i,j) = q(1,i,j)
            qt(2,i,j) = q(2,i,j)/q(1,i,j)
            qt(3,i,j) = q(3,i,j)/q(1,i,j)
            qt(4,i,j) = gamma1*(q(4,i,j) - 0.5d0*(q(2,i,j)**2 + 
     &           q(3,i,j)**2)/q(1,i,j))
 10   continue
c         
      return
      end
c
c     -----------------------------------------------------
c
c     Construction of reflective boundary conditions from
c     mirrored primitive values and application in
c     conservative form in local patch
c
c     =====================================================
      subroutine ip2eurfl(q,mx,my,lb,ub,meqn,nc,idx,
     &     qex,xc,phi,vn,maux,auex,dx,time)
c     =====================================================

      implicit none

      common /param/  gamma,gamma1
      double precision    gamma,gamma1
      integer   mx, my, meqn, maux, nc, idx(2,nc), lb(2), 
     &     ub(2)
      double precision    q(meqn, mx, my), qex(meqn,nc), xc(2,nc), 
     &     phi(nc), vn(2,nc), auex(maux,nc), dx(2), time
c
c     Local variables
c
      integer   i, j, n, stride, getindx
      double precision    rho, u, v, p, vl
c
      stride = (ub(1) - lb(1))/(mx-1)
c
      do 100 n = 1, nc

         i = getindx(idx(1,n), lb(1), stride)
         j = getindx(idx(2,n), lb(2), stride)
c
         rho =  qex(1,n)
         u   = -qex(2,n)       
         v   = -qex(3,n)
         p   =  qex(4,n)
c
c        # Add boundary velocities if available
         if (maux.ge.2) then
            u = u + auex(1,n)
            v = v + auex(2,n)
         endif
c
         u = qex(2,n) + 2.d0*u 
         v = qex(3,n) + 2.d0*v
c
         q(1,i,j) = rho
         q(2,i,j) = u*rho
         q(3,i,j) = v*rho
         q(4,i,j) = p/gamma1 + 0.5d0*rho*(u**2 + v**2)
c
 100  continue
c
      return
      end
c