c***********************************************************************
c     This file contains all routines for the following B.C.'s:
c     inlet and outlet conditions in the x direction;
c     in the y and z directions: free-slip conditions on the wall (at
c     j=1,k=1), symmetry at the center (at j=ny, k=nz).
c     To include directly the four vertices of each horizontal plane,
c     the B.C.'s are imposed in the following way:
c     y-z limiting planes: k: 2->nz-1; j: 2->ny-1;
c     x-z limiting planes: k: 2->nz-1; i: 1->nx;
c     x-y limiting planes: j: 1->ny;   i: 1->nx;
c***********************************************************************
      SUBROUTINE bc_vector (u,v,w,nx,ny,nz,PARINT,PARDPR,COMM2D)
C
C     Name: bc_vector
C     Author: ?? - Anthony Leboissetier (for the MPI spatial version) (10/01)
C     Objective: BC on the veolcity field
C     Called by: time_step & momentum
C     It calls: inlet
C     Modifications :
C
      INCLUDE 'undefined.h'
      INCLUDE 'mpif.h'
      common /TIME/ t
      INTEGER nx, ny, nz, i, j, k, t, PARINT(100),NEIGHB(2),NPROCS
      INTEGER COMM2D,IERR00,STAT00(MPI_STATUS_SIZE)
      DOUBLE PRECISION u(nx,ny,nz), v(nx,ny,nz), w(nx,ny,nz)
      DOUBLE PRECISION h,x,r0,u0,r,zz,yy,y0,z0,PARDPR(100)
      DOUBLE PRECISION fluxout,inlet, outlet, vout, outlet2, vout2
      h=1.d0/(nx-2)
      NEIGHB(1) = PARINT(61)
      NEIGHB(2) = PARINT(62)
      NPROCS = PARINT(56)
      r0 = PARDPR(10)
      u0 = PARDPR(11)
c***
c     periodicity of the 3 velocity components in the x direction
c***
      IF(NEIGHB(1).LT.0) THEN
      do k=2,nz-1
         do j=2,ny-1
          v(1,j,k)= 0.d0
          v(2,j,k)= 0.d0 
          w(1,j,k)= 0.d0
          w(2,j,k)= 0.d0
         enddo
      enddo
      ENDIF
c***
      IF(NEIGHB(2).LT.0) THEN
      do k=2,nz-1
         do j=2,ny-1
          v(nx,j,k)= 0.d0
          w(nx,j,k)= 0.d0
         enddo
      enddo
      ENDIF
c***
c     y direction: free-slip on the wall for the parallel components,
c:    to zero the perpendicular one: 
c     grad u = grad w = 0;      v(.,1,.) = v(.,2,.) = 0;
c     symmetry at the center (w becomes the radial component!): 
c     grad u = grad w = 0;      v(.,ny,.) = 0 
c***
      do k=2,nz-1
         do i=1,nx
            u(i,1,k) = u(i,2,k)
            u(i,ny,k) = u(i,ny-1,k)
            v(i,1 ,k) = 0.d0
            v(i,2 ,k) = 0.d0
            v(i,ny,k) = 0.d0
            w(i,1 ,k) = w(i,2,k)
            w(i,ny,k) = w(i,ny-1,k)
         enddo 
      enddo 
c***
c     z direction: same as before with v <-> w 
c     grad u = grad v = 0;      w(.,.,1) = w(.,.,2) = 0;
c     grad u = grad v = 0;      w(.,.,nz) = 0 
c***
      do j=1,ny
         do i=1,nx
            u(i,j,1 ) = u(i,j,2)
            u(i,j,nz) = u(i,j,nz-1)
            v(i,j,1 ) = v(i,j,2)
            v(i,j,nz) = v(i,j,nz-1)
            w(i,j,1 ) = 0.d0
            w(i,j,2 ) = 0.d0
            w(i,j,nz) = 0.d0
         enddo 
      enddo 
c***
      outlet=0.
      outlet2=0.
      y0 = (ny-2)*h/2
      z0 = (nz-2)*h/2
      IF(NEIGHB(1).LT.0) THEN
        do k=2,nz-1
          zz = k*h
          do j=2,ny-1
            yy = j*h
            r = DSQRT((yy-y0)*(yy-y0) + (zz-z0)*(zz-z0))
            u(1,j,k)= inlet(r,nx,PARDPR)
            u(2,j,k)= inlet(r,nx,PARDPR)
c compute total entrance flux
            outlet = outlet + inlet(r,nx)
          enddo 
        enddo
c average entrance flux
        vout = outlet/ (ny-2)**2
        CALL MPI_SEND(vout,1,MPI_DOUBLE_PRECISION,NPROCS-1,100,COMM2D,
     &                IERR00)
c      WRITE(*,*)'BCV 0 envoie',vout,' à',NPROCS-1
      ENDIF
c
      IF(NEIGHB(2).LT.0) THEN
        do k=2,nz-1
          do j=2,ny-1
c compute total exit flux
            outlet2 = outlet2 + u(nx-1,j,k)
          enddo 
        enddo
c average exit flux
        vout2 = outlet2/ (ny-2)**2
c
        CALL MPI_RECV(vout,1,MPI_DOUBLE_PRECISION,0,100,COMM2D,STAT00,
     &                IERR00)
c      WRITE(*,*)'BCV',NPROCS-1,' reçoit',vout,' de 0'
        do k=2,nz-1
          do j=2,ny-1
            u(nx,j,k) = u(nx-1,j,k) + vout - vout2
c	    write(*,*)PARINT(55),nx,j,k,u(nx,j,k)
          enddo
        enddo
      ENDIF
c
      return
      end
c***********************************************************************
      DOUBLE PRECISION FUNCTION inlet(r,nx,PARDPR)
      include 'undefined.h'
      integer nx
c      common /INOUT/ tau
      double precision r, tau, h, pi, r0, u0, vscale
      DOUBLE PRECISION PARDPR(100)
      h=1.d0/(nx-2)
      pi=4.d0*datan(1.d0)

      tau = PARDPR(1)
      r0 = PARDPR(10)
      u0 = PARDPR(11)
      vscale = PARDPR(2)

      vscale = (1.0d0/vscale) * (tau/h)

      if (r.lt.r0) then
         inlet=u0*vscale
      else 
         inlet=0.
      endif  

c      print*,u0,tau,h,r,r0

      return
      end
c***********************************************************************

      SUBROUTINE bc_tensor(S11,S22,S33,S12,S13,S23,nx,ny,nz,
     &                      PARINT)
C
C     Name: bc_tensor
C     Author: ?? - Anthony Leboissetier (for the MPI spatial version) (10/01)
C     Objective: BC on the tensor field S
C     Called by: momentum
C     It calls: 
C     Modifications :
C
      INCLUDE 'undefined.h'
      INTEGER nx, ny, nz, i, j, k, PARINT(100),NEIGHB(2)
      DOUBLE PRECISION S11(nx,ny,nz), S22(nx,ny,nz), S33(nx,ny,nz), 
     %         S12(nx,ny,nz), S13(nx,ny,nz), S23(nx,ny,nz)
c***
      NEIGHB(1) = PARINT(61)
      NEIGHB(2) = PARINT(62)
c***
c     periodicity in the x direction of all components of the stress 
c     tensor [ Note: there are a few commented out, this is because 
c                    they are not used in the code ]
c***
      IF(NEIGHB(1).LT.0) THEN
      do k=2,nz-1
         do j=2,ny-1
            S11(1 ,j,k) = 0.d0
         enddo 
      enddo
      ENDIF
c***
      IF(NEIGHB(2).LT.0) THEN
      do k=2,nz-1
         do j=2,ny-1
            S12(nx,j,k) = 0.d0
            S13(nx,j,k) = 0.d0
         enddo 
      enddo
      ENDIF
c***
c     y direction: free slip at the wall and simmetry at the center 
c     for indices 11, 33, 13 (this means zero gradient for both 
c     conditions); to zero all the rest
c***
      do k=2,nz-1
         do i=1,nx
c---            S11(i,1 ,k) = S11(i,2   ,k)
c---            S11(i,ny,k) = S11(i,ny-1,k)
            S22(i,1 ,k) = 0.d0
c---            S22(i,ny,k) = 0.d0
c---            S33(i,1 ,k) = S33(i,2,k)
c---            S33(i,ny,k) = S33(i,ny-1,k)
c---            S12(i,1 ,k) = 0.d0
            S12(i,2 ,k) = 0.d0
            S12(i,ny,k) = 0.d0
c---            S13(i,1 ,k) = S13(i,2,k)
c---            S13(i,ny,k) = S13(i,ny-1,k)
c---            S23(i,1 ,k) = 0.d0
            S23(i,2 ,k) = 0.d0
            S23(i,ny,k) = 0.d0
         enddo 
      enddo
c***
c     z direction: free slip at the wall and simmetry at the center 
c     for indices 11, 22, 12 (this means zero gradient for both 
c     conditions); to zero all the rest
c***
      do j=1,ny
         do i=1,nx
c---            S11(i,j,1 ) = S11(i,j,2)
c---            S11(i,j,nz) = S11(i,j,nz-1)
c---            S22(i,j,1 ) = S22(i,j,2)
c---            S22(i,j,nz) = S22(i,j,nz-1)
            S33(i,j,1 ) = 0.d0
c---            S33(i,j,nz) = 0.d0
c---            S12(i,j,1 ) = S12(i,j,2)
c---            S12(i,j,nz) = S12(i,j,nz-1)
c---            S13(i,j,1 ) = 0.d0
            S13(i,j,2 ) = 0.d0
            S13(i,j,nz) = 0.d0
c---            S23(i,j,1 ) = 0.d0
            S23(i,j,2 ) = 0.d0
            S23(i,j,nz) = 0.d0
         enddo 
      enddo
c***
      return
      end
c***********************************************************************
      SUBROUTINE bc_scalar(scal,nx,ny,nz,PARINT)
C
C     Name: bc_scalar
C     Author: ?? - Anthony Leboissetier (for the MPI spatial version) (10/01)
C     Objective: BC on the pressure field
C     Called by: rstrct, interp, addint & relax
C     It calls: 
C     Modifications :
C
      INCLUDE 'undefined.h'
      INTEGER nx,ny,nz, i, j, k, PARINT(100),NEIGHB(2)
      DOUBLE PRECISION scal(nx,ny,nz)
c***
      NEIGHB(1) = PARINT(61)
      NEIGHB(2) = PARINT(62)
c***
c     periodicity in the x direction
c***  
      IF(NEIGHB(1).LT.0) THEN
      do k=2,nz-1
         do j=2,ny-1
            scal(1 ,j,k) = scal(2,j,k)
         enddo
      enddo
      ENDIF
c***  
      IF(NEIGHB(2).LT.0) THEN
      do k=2,nz-1
         do j=2,ny-1
            scal(nx,j,k) = scal(nx-1,j,k)
         enddo
      enddo
      ENDIF
c***  
c     gradient equal to zero in the y and z directions
c***  
      do k=2,nz-1
         do i=1,nx
            scal(i,1,k ) = scal(i,2   ,k)
            scal(i,ny,k) = scal(i,ny-1,k)
         enddo
      enddo 
      do j=1,ny
         do i=1,nx
            scal(i,j,1 ) = scal(i,j,2)
            scal(i,j,nz) = scal(i,j,nz-1)
         enddo
      enddo 
c***
      return
      end
c***********************************************************************
      SUBROUTINE bc_c(scal,nx,ny,nz,PARINT,PARDPR)
C
C     Name: bc_c
C     Author: ?? - Anthony Leboissetier (for the MPI spatial version) (10/01)
C     Objective: BC on the volume fraction field
C     Called by: timestep & swp
C     It calls: 
C     Modifications :
C
      INCLUDE 'undefined.h'
      INTEGER nx,ny,nz, i, j, k, PARINT(100),NEIGHB(2),PPRANK
      DOUBLE PRECISION scal(nx,ny,nz),r0,y0,z0,zz,yy,r,h, PARDPR(100)
      h=1.d0/(nx-2)
      PPRANK = PARINT(55)
      NEIGHB(1) = PARINT(61)
      NEIGHB(2) = PARINT(62)
c      write(*,*)'BCC PPRANK',PPRANK,NEIGHB(1),NEIGHB(2)
      r0 = PARDPR(10)
c***  
c     periodicity in the x direction
c***  
      IF(NEIGHB(1).LT.0) THEN
c      write(*,*)'BCC2 PPRANK',PPRANK,NEIGHB(1),NEIGHB(2)
      y0 = (ny-1)*h/2
      z0 = (nz-1)*h/2
      do k=2,nz-1
         zz = k*h
         do j=2,ny-1
            yy = j*h
            r = DSQRT((yy-y0)*(yy-y0) + (zz-z0)*(zz-z0))
            if (dabs(r) .le. r0) then
              scal(2,j,k) = 1.d0
c	      write(*,*)'BC_C PPRANK',PPRANK,2,j,k,scal(2,j,k)
            else
              scal(2,j,k)= 0.d0
            endif
            scal(1 ,j,k) = scal(2,j,k)
         enddo 
      enddo
      ENDIF
c***  
      IF(NEIGHB(2).LT.0) THEN
      do k=2,nz-1
         do j=2,ny-1
            scal(nx,j,k) = scal(nx-1,j,k)
         enddo 
      enddo
      ENDIF
c***  
c     gradient equal to zero in the y and z directions
c***  
      do k=2,nz-1
         do i=1,nx
            scal(i,1,k ) = scal(i,2   ,k)
            scal(i,ny,k) = scal(i,ny-1,k)
         enddo
      enddo 
      do j=1,ny
         do i=1,nx
            scal(i,j,1 ) = scal(i,j,2)
            scal(i,j,nz) = scal(i,j,nz-1)
         enddo
      enddo 
c***
      return
      end
c***********************************************************************
      SUBROUTINE bc_resid(scal,nx,ny,nz,PARINT)
C
C     Name: bc_resid
C     Author: ?? - Anthony Leboissetier (for the MPI spatial version) (10/01)
C     Objective: BC on the the pressure equation
C     Called by: resid
C     It calls: 
C     Modifications :
C
      INCLUDE 'undefined.h'
      INTEGER nx,ny,nz, i, j, k,PARINT(100),NEIGHB(2)
      DOUBLE PRECISION scal(nx,ny,nz)
c***
      NEIGHB(1) = PARINT(61)
      NEIGHB(2) = PARINT(62)
c*** 
c     periodicity in the x direction
c***
      IF(NEIGHB(1).LT.0) THEN
      do k=2,nz-1
         do j=2,ny-1
            scal(1 ,j,k) = scal(2,j,k)
         enddo
      enddo
      ENDIF
c***
      IF(NEIGHB(2).LT.0) THEN
      do k=2,nz-1
         do j=2,ny-1
            scal(nx,j,k) = scal(nx-1,j,k)
         enddo
      enddo
      ENDIF
c***  
c     residue equal to zero in the y and z directions
c***  
      do k=2,nz-1
         do i=1,nx
            scal(i,1 ,k) = 0.d0
            scal(i,ny,k) = 0.d0
         enddo
      enddo
      do j=1,ny
         do i=1,nx
            scal(i,j,1 ) = 0.d0
            scal(i,j,nz) = 0.d0
         enddo
      enddo
c***  
      return
      end
c***********************************************************************
      SUBROUTINE bc_coef(a,b,c,cc,r1,r2,nx,ny,nz,PARINT)
C
C     Name: bc_coef
C     Author: ?? - Anthony Leboissetier (for the MPI spatial version) (10/01)
C     Objective: BC on the a, b & c fields
C     Called by: momentum
C     It calls: 
C     Modifications :
C
      INCLUDE 'undefined.h'
      INTEGER nx, ny, nz, i, j, k, PARINT(100), NEIGHB(2)
      DOUBLE PRECISION a(nx,ny,nz), b(nx,ny,nz), c(nx,ny,nz)
      DOUBLE PRECISION cc(nx,ny,nz), r1, r2
c***
      NEIGHB(1) = PARINT(61)
      NEIGHB(2) = PARINT(62)
c***
c     need to define indices i=1 for a, j=1 for b and k=1 for c: use
c     periodic boundary conditions in x direction, simple average
c     in y and z directions
c***
      IF(NEIGHB(1).LT.0) THEN
      do k=2,nz-1
         do j=2,ny-1
            a(1 ,j,k) = a(2,j,k)
         enddo
      enddo
      ENDIF
      IF(NEIGHB(2).LT.0) THEN
      do k=2,nz-1
         do j=2,ny-1
            a(nx,j,k) = a(nx-1,j,k)
         enddo
      enddo
      ENDIF
      do k=2,nz-1
         do i=2,nx-1
            b(i,1,k) = 1.d0 / (r1*(cc(i,1,k) + cc(i,2,k)) + r2)
         enddo
      enddo
      do j=2,ny-1
         do i=2,nx-1
            c(i,j,1) = 1.d0 / (r1*(cc(i,j,1) + cc(i,j,2)) + r2)
         enddo
      enddo
c***
      return 
      end
c***********************************************************************
      SUBROUTINE bc_flux(vof1,vof3,vel,nx,ny,nz,indx,PARINT)
C
C     Name: ?? - bc_flux
C     Author: Anthony Leboissetier (for the MPI spatial version) (10/01)
C     Objective: BC on the the fluxes
C     Called by: swp
C     It calls: 
C     Modifications :
C
      INCLUDE 'undefined.h'
      INTEGER nx, ny, nz, indx, i, j, k, PARINT(100), NEIGHB(2)
      DOUBLE PRECISION vof1(nx,ny,nz), vof3(nx,ny,nz), vel(nx,ny,nz)
c***
      NEIGHB(1) = PARINT(61)
      NEIGHB(2) = PARINT(62)
c***
c     along x direction
c***
      if (indx .eq. 1) then
        IF(NEIGHB(1).LT.0) THEN
          do k=2,nz-1
            do j=2,ny-1
               vof3(1 ,j,k) = 0.d0
            enddo
          enddo
        ENDIF
        IF(NEIGHB(2).LT.0) THEN
          do k=2,nz-1
            do j=2,ny-1
               vof1(nx,j,k) = 0.d0
            enddo
          enddo
        ENDIF
c***
c     along y direction
c***
      elseif (indx .eq. 2) then
         do k=2,nz-1
            do i=2,nx-1
               vof3(i,1 ,k) = DMAX1(0.0d0,vel(i,2,k))
               vof1(i,ny,k) = 0.d0
            enddo
         enddo
c***
c     along z direction
c***
      elseif (indx .eq. 3) then
         do j=2,ny-1
            do i=2,nx-1
               vof3(i,j,1) = DMAX1(0.0d0,vel(i,j,2))
               vof1(i,j,nz) = 0.0d0
            enddo
         enddo
c***  
c     wrong value
c***  
      else 
         stop 'wrong value for indx in one of the swaps'
      endif 
c***
      return
      end