vtf-logo

src/1d/equations/euler/rp/bcouteu1.f

!-----------------------------------------------------------------------
! Physical boundary conditions for 1d Euler equations.
! Zero order outflow at both sides.
! Interface:
!   mx,my    := shape of grid function
!
!   u(,) := grid function
!
!   lb(1) := lower bound for grid
!   ub(1) := upper bound for grid
!   lbbnd(1) := lower bound for boundary region
!   ubnd(1) := upper bound for boundary region
!   shapebnd(1) := shape of boundary region 
!   xc(1) := lower left corner of grid
!   dx(1) := grid spacing
!   dir := at which side of the grid is the boundary?
!   bnd(,2,1) := lower left and upper right corner of global grid and 
!      of mb-1 internal boundary regions 
!
! Copyright (C) 2002 Ralf Deiterding
! Brandenburgische Universitaet Cottbus
!
! Copyright (C) 2003-2007 California Institute of Technology
! Ralf Deiterding, ralf@amroc.net
!
!-----------------------------------------------------------------------

      subroutine physbd(u,mx,lb,ub,lbbnd,ubbnd,shapebnd,
     &     xc,dx,dir,bnd,mb,time,meqn)

      implicit none

      integer   mb, dir, meqn, mx 
      integer   lb(1), ub(1), lbbnd(1), ubbnd(1), shapebnd(1)
      double precision u(meqn, mx), xc(1), dx(1), bnd(mb,2,1), time

!      Local variables
      integer   i, m, imin, imax, isym
      integer   stride, stridebnd, getindx

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!      See definition of member-function extents() in BBox.h 
!      for calculation of stride
         
      stride = (ub(1) - lb(1))/(mx-1)

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!     Find working domain
         
      imin = getindx(max(lbbnd(1), lb(1)), lb(1), stride)
      imax = getindx(min(ubbnd(1), ub(1)), lb(1), stride)
      
      if(imax .gt. mx .or. imin .lt. 1) then
         write(0,*)'INDEX ERROR in physbd'
      end if

!        Left Side --- Outflow
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      if (dir.eq.0) then
         do 200 i = imax, imin, -1
            do 200 m = 1, meqn
               u(m,i) = u(m,i+1)
 200     continue
         return
      end if

!        Right Side --- Outflow
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      if (dir.eq.1) then
         do 210 i = imin, imax
            do 210 m = 1, meqn
               u(m,i) = u(m,i-1)
 210     continue
         return
      end if      
c     
      return
      end

<