! $Id: boundary.f,v 1.9 2001/04/17 06:19:25 aake Exp $
!***********************************************************************
!
!  Boundary conditions, implemented by setting values in "ghost zones".
!  The range of indices for zone-centered variables is lb:ub, while the
!  range for face-centered variables is lb+1:ub.
!
!  In order to allow two succesive operators (each involving three 
!  nearest neighbors), one needs to allow 5+5 ghost zones, with the
!  following arrangement, with lb=6 and ub=mx-5.
!
!                        L
!    1...2...3...4...5...6... ... 
!  1...2...3...4...5...6... ... 
!
!                     ... n-7...n-6...n-5...n-4...n-3...n-2...n-1...n-0
!                  ... n-7...n-6...n-5...n-4...n-3...n-2...n-1...n-0
!                                      U
!
!  There are 10 ghost zones for zone centered variable and 11 for face
!  centered variables.  Thus, there are 5 boundary zone values to be
!  set for both zone- and face-centered variables.  Their indices are
!  1:lb-1 and ub+1:n for zone centered vars.  For face centered vars they
!  are 1:lb and ub+1:n.  
!
!  Previously, the choice was lb=3 and ub=mx-4, which required restoring
!  ghost zone values after every application of an operator:
!
!                L
!    1...2...3...4...5...6... ... 
!  1...2...3...4...5...6... ... 
!
!                     ... n-7...n-6...n-5...n-4...n-3...n-2...n-1...n-0
!                  ... n-7...n-6...n-5...n-4...n-3...n-2...n-1...n-0
!                                            U
!
!  The current version, with ifs on each line, covers both cases.  For
!  the mb=5 case one may remove the ifs and change the loop lenths to
!  1,mb.
!
!-----------------------------------------------------------------------
      subroutine extrap_full_log (f)
      include 'cparam.inc'
      include 'cwork.inc'
      real, dimension(mx,my,mz) :: f, scr
!hpf$ align with TheCube:: f, scr

      scr(lb:ub,:,:) = amax1(f(lb:ub,:,:),1e-20)
      scr(lb:ub,:,:) = alog(scr(lb:ub,:,:))
      call extrap_full(scr)
      scr = amax1(scr,-40.)
      if (lb .gt.  1) f( 1:lb-1,:,:) = exp(scr( 1:lb-1,:,:))
      if (ub .lt. mx) f(ub+1:mx,:,:) = exp(scr(ub+1:mx,:,:))

      end
!-----------------------------------------------------------------------
      subroutine extrap_full(f)
      include 'cparam.inc'
      real, dimension(mx,my,mz) :: f
!hpf$ align with TheCube:: f

      do i=1,mb
        f(lb-i,:,:) = 2.*f(lb  ,:,:)-f(lb+i,:,:)
        f(ub+i,:,:) = 2.*f(ub  ,:,:)-f(ub-i,:,:)
      end do
      
      end
!-----------------------------------------------------------------------
      subroutine extrap_half (f)
      include 'cparam.inc'
      include 'cwork.inc'
      real, dimension(mx,my,mz) :: f
!hpf$ align with TheCube:: f
      real, dimension(my,mz) :: flb, fub

      flb = 2./8.*(15.*f(lb+1,:,:)-10.*f(lb+2,:,:)+3.*f(lb+3,:,:))
      fub = 2./8.*(15.*f(ub  ,:,:)-10.*f(ub-1,:,:)+3.*f(ub-2,:,:))
      do i=1,mb
        f(lb-i+1,:,:) = flb - f(lb+i  ,:,:)
        f(ub+i  ,:,:) = fub - f(ub-i+1,:,:)
      end do
      i=mb+1
      f(lb-i+1,:,:) = flb - f(lb+i  ,:,:)

      end
*-----------------------------------------------------------------------
      subroutine extrap_half_scalar (f)
      include 'cparam.inc'
      include 'cwork.inc'
      real, dimension(mx) :: f

      flb = 2./8.*(15.*f(lb+1)-10.*f(lb+2)+3.*f(lb+3))
      fub = 2./8.*(15.*f(ub  )-10.*f(ub-1)+3.*f(ub-2))
      do i=1,mb
        f(lb-i+1) = flb - f(lb+i  )
        f(ub+i  ) = fub - f(ub-i+1)
      end do
      i=mb+1
      f(lb-i+1) = flb - f(lb+i)

      end
*-----------------------------------------------------------------------
      subroutine extrap_lower_half (f)
      include 'cparam.inc'
      include 'cwork.inc'
      real, dimension(mx,my,mz) :: f
!hpf$ align with TheCube:: f
      real, dimension(my,mz) :: fub

      fub = 2./8.*(15.*f(ub+1,:,:)-10.*f(ub+2,:,:)+3.*f(ub+3,:,:))
      do i=1,mb
        f(ub+i,:,:) = fub - f(lb-1+1,:,:)
      end do

      end
*-----------------------------------------------------------------------
      subroutine antisymm_full (f)
      include 'cparam.inc'
      real, dimension(mx,my,mz) :: f
!hpf$ align with TheCube:: f

      do i=1,mb
        f(lb-i,:,:) = -f(lb+i,:,:)
        f(ub+i,:,:) = -f(ub-i,:,:)
      end do
      f(lb  ,:,:) = 0.0
      f(ub  ,:,:) = 0.0

      end
*---------------------------------------------------------------------
      subroutine antisymm_half_scalar (f)
      include 'cparam.inc'
      real, dimension(ub+3) :: f

      do i=1,mb
        f(lb-i+1) = - f(lb+i  )
        f(ub+i  ) = - f(ub-i+1)
      end do
      i=mb+1
      f(lb-i+1) = - f(lb+i)

      end
*---------------------------------------------------------------------
      subroutine antisymm_half (f)
      include 'cparam.inc'
      real, dimension(mx,my,mz) :: f
!hpf$ align with TheCube:: f

      do i=1,mb
        f(lb-i+1,:,:) = - f(lb+i  ,:,:)
        f(ub+i  ,:,:) = - f(ub-i+1,:,:)
      end do
      i=mb+1
      f(lb-i+1,:,:) = - f(lb+i,:,:)

      end
*---------------------------------------------------------------------
      subroutine symmetric_full(f)
      include 'cparam.inc'
      real, dimension(mx,my,mz) :: f
!hpf$ align with TheCube:: f

      do i=1,mb
        f(lb-i,:,:) = f(lb+i,:,:)
        f(ub+i,:,:) = f(ub-i,:,:)
      end do

      end
*---------------------------------------------------------------------
      subroutine symmetric_half(f)
      include 'cparam.inc'
      real, dimension(mx,my,mz) :: f
!hpf$ align with TheCube:: f

      do i=1,mb
        f(lb-i+1,:,:) = f(lb+i  ,:,:)
        f(ub+i  ,:,:) = f(ub-i+1,:,:)
      end do
      i=mb+1
      f(lb-i+1,:,:) = f(lb+i,:,:)

      end
*---------------------------------------------------------------------
      subroutine boundary_force(f, e)
      include 'cparam.inc'
      real, dimension(mx,my,mz) :: f, e
!hpf$ align with TheCube:: f, e
!
! correct pressure for external force at the upper and lower boundary
!
      do i=1,mb
        f(lb-i,:,:) = f(lb+i,:,:) + i*2.*e(lb  ,:,:)
        f(ub+i,:,:) = f(ub-i,:,:) - i*2.*e(ub  ,:,:)
      end do

      end
*---------------------------------------------------------------------
      subroutine boundary_grav (f, rho, gg)
      include 'cparam.inc'
      include 'cdata.inc'
      real, dimension(mx,my,mz) :: f, rho
!hpf$ align with TheCube:: f, rho
!
! correct pressure for external force at the upper and lower boundary
!
      do i=1,mb
        f(lb-i,:,:) = f(lb+i,:,:) 
     &      + (i*2.*dx*grav)*rho(lb  ,:,:)
        f(ub+i,:,:) = f(ub-i,:,:) 
     &      - (i*2.*dx*grav)*rho(ub  ,:,:)
      end do

      end
*--------------------------------------------------------------------
      subroutine antisymm_half_b (f, flux_lower, flux_upper)
      include 'cparam.inc'
      real, dimension(mx,my,mz) :: f
!hpf$ align with TheCube:: f
!
! add a mean flux throught the lower and upper bdry
!
      flb = 2*flux_lower
      fub = 2*flux_upper
      do i=1,mb
        f(lb-i+1,:,:) = flb - f(lb+i  ,:,:)
        f(ub+i  ,:,:) = fub - f(ub-i+1,:,:)
      end do
      i=mb+1
      f(lb-i+1,:,:) = flb - f(lb+i,:,:)

      end
*---------------------------------------------------------------------
      subroutine  open_upper_full (f)
      include 'cparam.inc'
      include 'cdata.inc'
      include 'cwork.inc'
      real, dimension(mx,my,mz) :: f
!hpf$ align with TheCube:: f
!
! symmetric around the position mx-4
!
!      cb=0.1
      cb = 0.0
      do i=1,mb
        f(ub+i,:,:) = (1.-cb*2**(i-1))*f(ub-i,:,:)
      end do

      END
*---------------------------------------------------------------------
      subroutine  open_upper_half (f)
      include 'cparam.inc'
      include 'cdata.inc'
      include 'cwork.inc'
      real, dimension(mx,my,mz) :: f
!hpf$ align with TheCube:: f
*
* symmetric around the position mx-4+1/2, with down-trend
*
!      cb=0.1
      cb = 0.0
      do i=1,mb
        f(ub+i,:,:) = (1.-cb*2**(i-1))*f(ub-i+1,:,:)
      end do

      END
*---------------------------------------------------------------------
      subroutine antisymm_upper_half (f)
      include 'cparam.inc'
      include 'cdata.inc'
      include 'cwork.inc'
      real, dimension(mx,my,mz) :: f
!hpf$ align with TheCube:: f
!
! antisymmetric around the position mx-4+1/2
!
      do i=1,mb
        f(lb-i+1,:,:) = - f(lb+i,:,:)
      end do
      i=mb+1
      f(lb-i+1,:,:) = - f(lb+i,:,:)

      END
*---------------------------------------------------------------------
      subroutine extrap_half_log (f)
      include 'cparam.inc'
      include 'cwork.inc'
      real, dimension(mx,my,mz) :: f, scr
!hpf$ align with TheCube:: f, scr
*
* extrapolate the values around the position 4 and mx-4
*
      scr(lb+1:ub,:,:) = amax1( f(lb+1:ub,:,:),1e-20)
      scr(lb+1:ub,:,:) = alog(scr(lb+1:ub,:,:))
      call extrap_half (scr)
      scr = amax1(scr,-40.)
      f = exp(scr)

      end
*---------------------------------------------------------------------
      subroutine extrap_half_b (f, tmp1, tmp2)
      include 'cparam.inc'
      include 'cdata.inc'
      include 'cwork.inc'
      include 'stagger.ifc'
      include 'df.ifc'
      real, dimension(mx,my,mz) :: f, tmp1, tmp2
!hpf$ align with TheCube:: f, tmp1, tmp2

      do i=1,mb
        f(lb-i+1,:,:) = 2.*tmp1(lb,:,:) - f(lb+i  ,:,:)
        f(ub+i  ,:,:) = 2.*tmp2(ub,:,:) - f(ub-i+1,:,:)
      end do
      i=mb+1
      f(lb-i+1,:,:) = 2.*tmp1(lb,:,:) - f(lb+i,:,:)

      end
*---------------------------------------------------------------------
      subroutine set_full(f,a,b)
      include 'cparam.inc'
      real, dimension(mx,my,mz) :: f
!hpf$ align with TheCube:: f

      do i=0,mb
        f(lb-i,:,:) = a
        f(ub+i,:,:) = b
      end do

      end
*---------------------------------------------------------------------
      subroutine set_half(f,a,b)
      include 'cparam.inc'
      real, dimension(mx,my,mz) :: f
!hpf$ align with TheCube:: f

      do i=1,mb
        f(lb-i+1,:,:) = a
        f(ub+i  ,:,:) = b
      end do
      i=mb+1
      f(lb-i+1,:,:) = a

      end
*---------------------------------------------------------------------
      subroutine mixed_full(f,a)
      include 'cparam.inc'
      real, dimension(mx,my,mz) :: f, scr
!hpf$ align with TheCube:: f

! Outflow
      scr = amax1(f,1e-30)
      scr = alog(scr)
      do i=1,mb
        scr(lb-i,:,:) = 2.*scr(lb  ,:,:)-scr(lb+i,:,:)
!       f(lb-i,:,:) = f(lb,:,:)
      enddo
      scr = amax1(scr,-70.)
      f(1:lb-1,:,:)=exp(scr(1:lb-1,:,:))
! Inflow
      do i=0,mb
        f(ub+i,:,:) = a
      end do

	end
*---------------------------------------------------------------------
      subroutine mixed_half(f,a)
      include 'cparam.inc'
      real, dimension(mx,my,mz) :: f, scr
!hpf$ align with TheCube:: f

! Outflow
      do i=1,mb
!       scr(lb-i,:,:) = 2.*scr(lb  ,:,:)-scr(lb+i,:,:)
!       f(lb-i,:,:) = f(lb+i,:,:)
        f(lb-i,:,:) = f(lb+i-1,:,:)
      enddo
! Inflow
      do i=1,mb
        f(ub+i  ,:,:) = f(ub-i+1,:,:)
      end do

      end
*---------------------------------------------------------------------
      subroutine periodic(f)
      include 'cparam.inc'
      real, dimension(mx,my,mz) :: f
!hpf$ align with TheCube:: f

      end
*---------------------------------------------------------------------