; $Id: boundary.pro,v 1.1.1.1 2000/06/05 18:00:57 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.
;
;-----------------------------------------------------------------------
PRO extrap_full_log, f
@cparam
@cwork
      scr = f
      scr(lb:ub,*,*) = amax1(f(lb:ub,*,*),1e-20)
      scr(lb:ub,*,*) = alog(scr(lb:ub,*,*))
      extrap_full, scr
      scr = amax1(scr,-40.)
      f = exp(scr)

END
;-----------------------------------------------------------------------
PRO extrap_full, f
@cparam

      for i=1,mb do begin
        f(lb-i,*,*) = 2.*f(lb  ,*,*)-f(lb+i,*,*)
        f(ub+i,*,*) = 2.*f(ub  ,*,*)-f(ub-i,*,*)
      endfor

END
;-----------------------------------------------------------------------
PRO extrap_half, f
@cparam
@cwork

      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,*,*))
      for i=1,mb do begin
        f(lb-i+1,*,*) = flb - f(lb+i  ,*,*)
        f(ub+i  ,*,*) = fub - f(ub-i+1,*,*)
      endfor
      i=mb+1
      f(lb-i+1,*,*) = flb - f(lb+i  ,*,*)

END
;-----------------------------------------------------------------------
PRO extrap_half_scalar, f
@cparam
@cwork

      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))
      for i=1,mb do begin
        f(lb-i+1) = flb - f(lb+i  )
        f(ub+i  ) = fub - f(ub-i+1)
      endfor
      i=mb+1
      f(lb-i+1) = flb - f(lb+i)

END
;-----------------------------------------------------------------------
PRO extrap_lower_half, f
@cparam
@cwork

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

END
;-----------------------------------------------------------------------
PRO antisymm_full, f
@cparam

      for i=1,mb do begin
        f(lb-i,*,*) = -f(lb+i,*,*)
        f(ub+i,*,*) = -f(ub-i,*,*)
      endfor
      f(lb  ,*,*) = 0.0
      f(ub  ,*,*) = 0.0

END
;---------------------------------------------------------------------
PRO antisymm_half_scalar, f
@cparam

      for i=1,mb do begin
        f(lb-i+1) = - f(lb+i  )
        f(ub+i  ) = - f(ub-i+1)
      endfor
      i=mb+1
      f(lb-i+1) = - f(lb+i)

END
;---------------------------------------------------------------------
PRO antisymm_half, f
@cparam

      for i=1,mb do begin
        f(lb-i+1,*,*) = - f(lb+i  ,*,*)
        f(ub+i  ,*,*) = - f(ub-i+1,*,*)
      endfor
      i=mb+1
      f(lb-i+1,*,*) = - f(lb+i,*,*)

END
;---------------------------------------------------------------------
PRO symmetric_full, f
@cparam

      for i=1,mb do begin
        f(lb-i,*,*) = f(lb+i,*,*)
        f(ub+i,*,*) = f(ub-i,*,*)
      endfor

END
;---------------------------------------------------------------------
PRO symmetric_half, f
@cparam

      for i=1,mb do begin
        f(lb-i+1,*,*) = f(lb+i  ,*,*)
        f(ub+i  ,*,*) = f(ub-i+1,*,*)
      endfor
      i=mb+1
      f(lb-i+1,*,*) = f(lb+i,*,*)

END
;---------------------------------------------------------------------
PRO boundary_force, f, e
@cparam
;
; correct pressure for external force at the upper and lower boundary
;
      for i=1,mb do begin
        f(lb-i,*,*) = f(lb+i,*,*) + i*2.*e(lb  ,*,*)
        f(ub+i,*,*) = f(ub-i,*,*) - i*2.*e(ub  ,*,*)
      endfor

END
;---------------------------------------------------------------------
PRO boundary_grav, f, rho, gg
@cparam
@cdata
;
; correct pressure for external force at the upper and lower boundary
;
      for i=1,mb do begin
        f(lb-i,*,*) = f(lb+i,*,*) $
           + (i*2.*dx*grav)*rho(lb  ,*,*)
        f(ub+i,*,*) = f(ub-i,*,*) $
           - (i*2.*dx*grav)*rho(ub  ,*,*)
      endfor

END
;--------------------------------------------------------------------
PRO antisymm_half_b, f, flux_lower, flux_upper
@cparam
;
; add a mean flux throught the lower and upper bdry
;
      flb = 2*flux_lower
      fub = 2*flux_upper
      for i=1,mb do begin
        f(lb-i+1,*,*) = flb - f(lb+i  ,*,*)
        f(ub+i  ,*,*) = fub - f(ub-i+1,*,*)
      endfor
      i=mb+1
      f(lb-i+1,*,*) = flb - f(lb+i,*,*)

END
;---------------------------------------------------------------------
PRO  open_upper_full, f
@cparam
@cdata
@cwork
;
; symmetric around the position mx-4
;
;      cb=0.1
      cb = 0.0
      for i=1,mb do begin
        f(ub+i,*,*) = (1.-cb*2^(i-1))*f(ub-i,*,*)
      endfor

END
;---------------------------------------------------------------------
PRO  open_upper_half, f
@cparam
@cdata
@cwork
;
; symmetric around the position mx-4+1/2, with down-trend
;
;      cb=0.1
      cb = 0.0
      for i=1,mb do begin
        f(ub+i,*,*) = (1.-cb*2^(i-1))*f(ub-i+1,*,*)
      endfor

END
;---------------------------------------------------------------------
PRO antisymm_upper_half, f
@cparam
@cdata
@cwork
;
; antisymmetric around the position mx-4+1/2
;
      for i=1,mb do begin
        f(lb-i+1,*,*) = - f(lb+i,*,*)
      endfor
      f(lb,*,*) = 0.
      i=mb+1
      f(lb-i+1,*,*) = - f(lb+i,*,*)

END
;---------------------------------------------------------------------
PRO extrap_half_log, f
@cparam
@cwork
;
; 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,*,*))
      extrap_half, scr
      scr = amax1(scr,-40.)
      f = exp(scr)

END
;---------------------------------------------------------------------
PRO extrap_half_b, f, tmp1, tmp2
@cparam
@cdata
@cwork

      for i=1,mb do begin
        f(lb-i+1,*,*) = 2.*tmp1(lb,*,*) - f(lb+i  ,*,*)
        f(ub+i  ,*,*) = 2.*tmp2(ub,*,*) - f(ub-i+1,*,*)
      endfor
      i=mb+1
      f(lb-i+1,*,*) = 2.*tmp1(lb,*,*) - f(lb+i,*,*)

END