c***********************************************************************
c     This file contains all routines for the following B.C.'s:
c     rigid BC (u=v=w=0) in 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)
c***
      INCLUDE 'undefined.h'
      INCLUDE 'mpif.h'
      common /TIME/ t
      INTEGER nx, ny, nz, i, j, k, t
      INTEGER 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
      DOUBLE PRECISION fluxout,inlet, outlet, vout, outlet2, vout2
      h=1.d0/(nx-2)
c***
c     The 3 velocity components are zero on x=0, L_x
c***
      IF(NEIGHB(1).LT.0) THEN
      do k=2,nz-1
         do j=2,ny-1
          v(1,j,k)= -v(2,j,k) 
          w(1,j,k)= -w(2,j,k)
          u(1,j,k)= 0.d0
          u(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
          u(nx,j,k)= 0.d0
          v(nx,j,k)= v(nx-1,j,k)
          w(nx,j,k)= w(nx-1,j,k)
         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***
c$$$      outlet=0.
c$$$      outlet2=0.
c$$$      r0 = GP_DFETCH('radius_of_jet',1.d0,'SEVERE')
c$$$      u0 = GP_DFETCH('liquide_velocity',1.d0,'SEVERE')
c$$$      y0 = (ny-2)*h/2
c$$$      z0 = (nz-2)*h/2
c$$$      do k=2,nz-1
c$$$         zz = k*h
c$$$         do j=2,ny-1
c$$$            yy = j*h
c$$$               r = DSQRT((yy-y0)*(yy-y0) + (zz-z0)*(zz-z0))
c$$$               u(1,j,k)= inlet(r,nx)
c$$$               u(2,j,k)= inlet(r,nx)
c$$$c compute total flux
c$$$               outlet = outlet + inlet(r,nx)
c$$$               outlet2 = outlet2 + u(nx-1,j,k)
c$$$         enddo 
c$$$      enddo
c$$$c average entrance flux
c$$$      vout = outlet/ (ny-2)**2
c$$$c average exit flux
c$$$      vout2 = outlet2/ (ny-2)**2
c$$$      do k=2,nz-1
c$$$        do j=2,ny-1
c$$$            u(nx,j,k) = u(nx-1,j,k) + vout - vout2
c$$$        enddo
c$$$      enddo
      CALL BUPDAT3D(u,nx,ny,nz)
      CALL BUPDAT3D(v,nx,ny,nz)
      CALL BUPDAT3D(w,nx,ny,nz)
      return
      end
c***********************************************************************
      DOUBLE PRECISION FUNCTION inlet(r,nx)
      include 'undefined.h'
      include 'physics.h'
      integer nx
      double precision r, pi, r0, u0, scale
      DOUBLE PRECISION GP_DFETCH,tau_SI,timescale
      EXTERNAL GP_DFETCH
      h=1.d0/(nx-2)
      pi=4.d0*datan(1.d0)

      r0 = GP_DFETCH('radius',1.d0,'SEVERE')
      u0 = GP_DFETCH('liquid_velocity',1.d0,'SEVERE')
      
      scale = (1.0d0/vscale) * (tau/h)

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

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

      SUBROUTINE bc_tensor(S11,S22,S33,S12,S13,S23,nx,ny,nz)
c***
      INCLUDE 'undefined.h'
      INTEGER nx, ny, nz, i, j, k
      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***
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***
      CALL BUPDAT3D(S11,nx,ny,nz)
      CALL BUPDAT3D(S22,nx,ny,nz)
      CALL BUPDAT3D(S33,nx,ny,nz)
      CALL BUPDAT3D(S12,nx,ny,nz)
      CALL BUPDAT3D(S13,nx,ny,nz)
      CALL BUPDAT3D(S23,nx,ny,nz)
      return
      end
c***********************************************************************
      SUBROUTINE bc_scalar(scal,nx,ny,nz)
c***
      INCLUDE 'undefined.h'
      INTEGER nx,ny,nz, i, j, k
      DOUBLE PRECISION scal(nx,ny,nz)
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***
      CALL BUPDAT3D(scal,NX,NY,NZ)
      return
      end
c***********************************************************************
      SUBROUTINE bc_c(scal,nx,ny,nz)
c***
      INCLUDE 'undefined.h'
      INTEGER nx,ny,nz, i, j, k
      DOUBLE PRECISION scal(nx,ny,nz),y0,z0,zz,yy,r,h
      DOUBLE PRECISION GP_DFETCH
      EXTERNAL GP_DFETCH
      h=1.d0/(nx-2)
c***  
c     periodicity in the x direction
c***  
      IF(NEIGHB(1).LT.0) THEN
c      write(*,*)'BCC2 PPRANK',PPRANK,NEIGHB(1),NEIGHB(2)
      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***  
c      do k=2,nz-1
c         do i=1,nx
c            scal(i,1,k ) = scal(i,2   ,k)
c            scal(i,ny,k) = scal(i,ny-1,k)
c         enddo
c      enddo 
c      do j=1,ny
c         do i=1,nx
c            scal(i,j,1 ) = scal(i,j,2)
c            scal(i,j,nz) = scal(i,j,nz-1)
c         enddo
c      enddo 

      do j=2,ny-1
         do i=1,nx
            scal(i,j,1 ) = scal(i,j,2)
            scal(i,j,nz) = scal(i,j,nz-1)
         enddo
      enddo 
      
      do k=nz,1,-1
         do i=5,nx
            if (scal(i,2,k) .eq. 1) then
               scal(i-4,1,k ) = 1
               scal(i,ny,k) = scal(i,ny-1,k)
            else
               scal(i,1,k ) = scal(i,2   ,k)
               scal(i,ny,k) = scal(i,ny-1,k)
            endif               
         enddo
         do i=nx-3,nx
            scal(i,1,k ) = scal(i,2   ,k) 
            scal(i,ny,k) = scal(i,ny-1,k)
         enddo
      enddo 
c***
      CALL BUPDAT3D(scal,NX,NY,NZ)
      return
      end
c***********************************************************************
      SUBROUTINE bc_resid(scal,nx,ny,nz)
c***
      INCLUDE 'undefined.h'
      INTEGER nx,ny,nz, i, j, k
      DOUBLE PRECISION scal(nx,ny,nz)
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***  
      CALL BUPDAT3D(scal,NX,NY,NZ)
      return
      end
c***********************************************************************
      SUBROUTINE bc_coef(a,b,c,cc,r1,r2,nx,ny,nz)
c***
      INCLUDE 'undefined.h'
      INTEGER nx, ny, nz, i, j, k
      DOUBLE PRECISION a(nx,ny,nz), b(nx,ny,nz), c(nx,ny,nz)
      DOUBLE PRECISION cc(nx,ny,nz), r1, r2
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***
      CALL BUPDAT3D(a,nx,ny,nz)
      CALL BUPDAT3D(b,nx,ny,nz)
      CALL BUPDAT3D(c,nx,ny,nz)
      return 
      end
c***********************************************************************
      SUBROUTINE bc_flux(vof1,vof3,vel,nx,ny,nz,indx)
c***
      INCLUDE 'undefined.h'
      INTEGER nx, ny, nz, indx, i, j, k
      DOUBLE PRECISION vof1(nx,ny,nz), vof3(nx,ny,nz), vel(nx,ny,nz)
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***
      CALL BUPDAT3D(vof1,nx,ny,nz)
      CALL BUPDAT3D(vof3,nx,ny,nz)
      return
      end