*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
*
* [FILE]       boundaries.cpp
*
* [VERSION]    H3expresso  (c) 1995 Joan Masso, NCSA & UIB
*
* [PURPOSE]    Routines for the boundary conditions.
*
* [ROUTINES]   Boundaries
*              Onebound
*
* [COMMENTS]   H3expresso has VERY limited boundary control.
*
*$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

*==============================================================================
*
*  [ROUTINE NAME] Boundaries
*  [AUTHOR] Joan Masso, NCSA & UIB
*
*  [PURPOSE] Provide boundary conditions to all the metric variables that
*            have been evolved in the interior. H3expresso simplifies the
*            proccess by enforcing the same boundary condition on all
*            quantities thru a call to Onebound for each variable.
*
*  [ARGUMENTS] 
*     [INPUT]
*        nx,ny,nz   : grid sizes of the 3d cube.
*     [INPUT/OUTPUT] All the metric structure. Only the boundaries 
*                       of the evolved quantities are modified.
*  [INCLUDES]  metric.h  declares all the grid and metric arrays.
*
*  [VARIABLES]  ibound: hardwired choice for periodic or flat 
*                            Boundary conditions.
* 
*  [CALLED BY]  H3
*
*  [CALLS TO]   Onebound
*               invert 
*
*  [WARNINGS] The boundaries are NOT applied to the covariant 
*             down metric since it is NOT one of the evolved
*             quantities. Therefore, a final call to invert recomputes
*             all the down metric.
*
*             Note also that the grid vectors and time information
*             are passed but not used. As more sophisticated boundaries
*             are implemented, these quantities are necessary... voila.
*
*<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<


      subroutine Boundaries(
     &     nx,ny,nz,time,dt,dx,dy,dz,
     &     x,y,z,r,psi,
     &     alp,cux,cuy,cuz,rg,
     &     uxx,uxy,uxz,uyy,uyz,uzz,
     &     gxx,gxy,gxz,gyy,gyz,gzz,
     &     qxx,qxy,qxz,qyy,qyz,qzz,
     &     dxuxx,dxuxy,dxuxz,dxuyy,dxuyz,dxuzz,
     &     dyuxx,dyuxy,dyuxz,dyuyy,dyuyz,dyuzz,
     &     dzuxx,dzuxy,dzuxz,dzuyy,dzuyz,dzuzz
     &     )
      implicit none
         integer nx,ny,nz
         real time,dt,dx,dy,dz
         real, dimension(:,:,:) ::x,y,z,r,psi,
     &     alp,cux,cuy,cuz,rg,
     &     uxx,uxy,uxz,uyy,uyz,uzz,
     &     gxx,gxy,gxz,gyy,gyz,gzz,
     &     qxx,qxy,qxz,qyy,qyz,qzz,
     &     dxuxx,dxuxy,dxuxz,dxuyy,dxuyz,dxuzz,
     &     dyuxx,dyuxy,dyuxz,dyuyy,dyuyz,dyuzz,
     &     dzuxx,dzuxy,dzuxz,dzuyy,dzuyz,dzuzz


       INTERFACE 

       subroutine Onebound(ibound,nx,ny,nz,var)

       integer ibound
       integer nx,ny,nz
       real var(:,:,:)
         
        end subroutine Onebound

       subroutine invert(flag, 
     &     in_xx,in_xy,in_xz,in_yy,in_yz,in_zz,
     &     rg, out_xx,out_xy,out_xz,out_yy,out_yz,out_zz)
c      integer nx,ny,nz
      integer flag

      real rg(:,:,:)
      real in_xx(:,:,:),in_xy(:,:,:),in_xz(:,:,:),
     &     in_yy(:,:,:),in_yz(:,:,:),in_zz(:,:,:) 

      real out_xx(:,:,:),out_xy(:,:,:),out_xz(:,:,:),
     &     out_yy(:,:,:),out_yz(:,:,:),out_zz(:,:,:)   

       end subroutine invert
    
        END INTERFACE

*******************************************************************************
*
* [INCLUDE FILE] metric.h
* [VERSION]      H3expresso  (c) 1995 Joan Masso, NCSA & UIB
*
* [PURPOSE]      Declare the "Data Structures" of H3expresso, that is,
*                all the grid and metric variables
*
* [INCLUDED IN]  H3
*                Initial
*                Method
*                Sources
*                Fluxes
*                Boundaries 
*
******************************* 3-D vars **************************************

c $x,y,z,r$ : Spatial vectors and radius
c      real x(mx,my,mz),y(mx,my,mz),z(mx,my,mz),r(mx,my,mz)

c $\Psi$ : conformal factor Psi 
c      real psi(mx,my,mz)

c $\alpha$ : Lapse 
c      real alp(mx,my,mz)

c $\Gamma^i$ : Momentum constraint related var  
c      real cux(mx,my,mz),cuy(mx,my,mz),cuz(mx,my,mz)

c $g^{ij}$ : Contravariant "upper" metric 
c      real uxx(mx,my,mz),uxy(mx,my,mz),uxz(mx,my,mz)
c     &    ,uyy(mx,my,mz),uyz(mx,my,mz),uzz(mx,my,mz)

c $g_{ij}$ : Covariant "down" metric
c           ==> It is NOT evolved but computed from the contravariant metric.
c      real gxx(mx,my,mz),gxy(mx,my,mz),gxz(mx,my,mz)
c     &    ,gyy(mx,my,mz),gyz(mx,my,mz),gzz(mx,my,mz)
c $\sqrt{g}$ : root of the covariant metric determinant 
c      real rg(mx,my,mz)

c $q^{ij}  = 2 \sqrt{g} k^{ij} $ : Extrimsic curvature related var 
c      real qxx(mx,my,mz),qxy(mx,my,mz),qxz(mx,my,mz)
c     &    ,qyy(mx,my,mz),qyz(mx,my,mz),qzz(mx,my,mz)

c $\partial_k g^{ij}$ : Derivatives of the contravariant upper metric 
c     x derivs
c      real dxuxx(mx,my,mz),dxuxy(mx,my,mz),dxuxz(mx,my,mz)
c     &    ,dxuyy(mx,my,mz),dxuyz(mx,my,mz),dxuzz(mx,my,mz)
c     y derivs
c      real dyuxx(mx,my,mz),dyuxy(mx,my,mz),dyuxz(mx,my,mz)
c     &    ,dyuyy(mx,my,mz),dyuyz(mx,my,mz),dyuzz(mx,my,mz)
c     z derivs
c      real dzuxx(mx,my,mz),dzuxy(mx,my,mz),dzuxz(mx,my,mz)
c     &    ,dzuyy(mx,my,mz),dzuyz(mx,my,mz),dzuzz(mx,my,mz)

*******************************************************************************




      integer ibound

c     H3expresso: Hard-wired boundary selection: see Onebound

c      write (*,*) 'BOundaries : entered '            
      ibound = 1

      call Onebound(ibound,nx,ny,nz, alp)
      
      call Onebound(ibound,nx,ny,nz, uxx)
      call Onebound(ibound,nx,ny,nz, uxy)
      call Onebound(ibound,nx,ny,nz, uxz)
      call Onebound(ibound,nx,ny,nz, uyy)
      call Onebound(ibound,nx,ny,nz, uyz)
      call Onebound(ibound,nx,ny,nz, uzz)

      call Onebound(ibound,nx,ny,nz, cux)
      call Onebound(ibound,nx,ny,nz, cuy)
      call Onebound(ibound,nx,ny,nz, cuz)
      
      call Onebound(ibound,nx,ny,nz, qxx)
      call Onebound(ibound,nx,ny,nz, qxy)
      call Onebound(ibound,nx,ny,nz, qxz)
      call Onebound(ibound,nx,ny,nz, qyy)
      call Onebound(ibound,nx,ny,nz, qyz)
      call Onebound(ibound,nx,ny,nz, qzz)
      
      call Onebound(ibound,nx,ny,nz, dxuxx)
      call Onebound(ibound,nx,ny,nz, dxuxy)
      call Onebound(ibound,nx,ny,nz, dxuxz)
      call Onebound(ibound,nx,ny,nz, dxuyy)
      call Onebound(ibound,nx,ny,nz, dxuyz)
      call Onebound(ibound,nx,ny,nz, dxuzz)
      
      call Onebound(ibound,nx,ny,nz, dyuxx)
      call Onebound(ibound,nx,ny,nz, dyuxy)
      call Onebound(ibound,nx,ny,nz, dyuxz)
      call Onebound(ibound,nx,ny,nz, dyuyy)
      call Onebound(ibound,nx,ny,nz, dyuyz)
      call Onebound(ibound,nx,ny,nz, dyuzz)
      
      call Onebound(ibound,nx,ny,nz, dzuxx)
      call Onebound(ibound,nx,ny,nz, dzuxy)
      call Onebound(ibound,nx,ny,nz, dzuxz)
      call Onebound(ibound,nx,ny,nz, dzuyy)
      call Onebound(ibound,nx,ny,nz, dzuyz)
      call Onebound(ibound,nx,ny,nz, dzuzz)

c     Compute covariant metric
      call invert(1,
     &     uxx,uxy,uxz,uyy,uyz,uzz,
     &     rg,gxx,gxy,gxz,gyy,gyz,gzz)
c       write(6,*) maxval(gxx),
c     &           minval(gxx)
      
      return
      end
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
      


*==============================================================================
*
*  [ROUTINE NAME] Onebound
*  [AUTHOR] Joan Masso, NCSA & UIB
*
*  [PURPOSE] Simple boundary conditions on one variable
*
*  [ARGUMENTS] 
*     [INPUT]
*        ibound     : FLag controlling the boundary applied
*                        1 = Periodic 
*                        2 = flat: copying the neighbouring plane.
*        nx,ny,nz   : grid sizes of the 3d cube.
*     [INPUT/OUTPUT] 
*        var        :  The 3d variable whose interior is used to
*                      compute the exterior boundaries.
* 
*  [CALLED BY]  Boundaries
*
*<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      
      subroutine Onebound(ibound,nx,ny,nz,var)
      integer ibound
      integer nx,ny,nz
      real var(:,:,:)
      

      
c     Periodic Boundaries.......................................
c      write(6,*) 'onebound entered'
      if (ibound.eq.1) then
         var(1,:,:) = var(nx-1,:,:)
         var(:,1,:) = var(:,ny-1,:)
         var(:,:,1) = var(:,:,nz-1)
         
         var(nx,:,:) = var(2,:,:)
         var(:,ny,:) = var(:,2,:)
         var(:,:,nz) = var(:,:,2)
      endif
c     Flat Boundaries (zero order extrapolation)................
      if (ibound.eq.2) then
         var(1,:,:) = var(2,:,:)
         var(:,1,:) = var(:,2,:)
         var(:,:,1) = var(:,:,2)
         
         var(nx,:,:) = var(nx-1,:,:)
         var(:,ny,:) = var(:,ny-1,:)
         var(:,:,nz) = var(:,:,nz-1)
      endif
           
      return
      end
*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>















