      SUBROUTINE makebcf(u,v,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION u(nx,ny),v(nx,ny),h,x,r1,r2
      DOUBLE PRECISION fluxout,inlet, outlet, vout
c     , sum
      INTEGER i,j
      r1=0.25d0
      r2 =0.35d0
      h=1.d0/(nx-2)

      do i=1,nx
         u(i,1)= 0.d0
         u(i,2)= 0.d0
         u(i,ny)= 0.d0
      enddo

      outlet=0.
c      sum = 0.

      do i=2,nx-1
         x=(dfloat(i-1)/(nx-2))
         v(i,1)= inlet(x,nx)
         v(i,2)= inlet(x,nx)
c compute total flux
         outlet = outlet + inlet(x,nx)*(i-1.5d0)
c         sum = sum + (i-1.5d0)
      enddo
c average flux
      vout = outlet/ ( nx*(nx-1)*0.5d0 - 1.d0 - 1.5d0*(nx-2))
c      fluxout = 0. 
c impose constant outlet velocity with same flux
      do i=2,nx-1
         v(i,ny)= vout
c         fluxout = fluxout + v(i,ny)*(i-1.5d0)
      enddo
c      write(6,*) outlet, fluxout
c      write(6,*) "sum", sum,
c      stop
      do i=1,ny
         u(1,i) = 0.d0
         u(2,i) = 0.d0
         v(1,i) = v(2,i)
         u(nx,i) = 0.d0
         v(nx,i) = v(nx-1,i)
      enddo

      return
      end

      DOUBLE PRECISION FUNCTION inlet(x,nx)
      include 'undefined.h'
      integer nx
      double precision x, vel1, vel2, tau, h
      common /INOUT/ vel1, vel2, tau
      double precision r1,r2,r3,d1,d2
      INTEGER i, ifetch, dfetch
      h=1.d0/(nx-2)

      i=dfetch("inner_nozzle_radius",r1)
      i=dfetch("separator_width",d1)
      i=dfetch("outer_nozzle_width",d2) 
      if(r1.le.0.) stop'r1<=0'
      if(d1.lt.0.) stop'd1<0'
      if(d2.lt.0.) stop'd2<0'

      r2=r1+d1
      r3=r2+d2
      
      if(x.lt.r1) then
         inlet=vel1*tau/h
      else if(x.lt.r2) then
         inlet=0.
      else if(x.lt.r3) then
         inlet=vel2*tau/h
      else 
         inlet=0.
      endif
      return
      end

      SUBROUTINE makebcsf(S11,S22,S12,nx,ny)
      include 'undefined.h'
      integer nx,ny
      DOUBLE PRECISION s11(nx,ny),s22(nx,ny),s12(nx,ny)
      INTEGER i,j

      do i=1,nx
         S12(i,1)= 0.d0
         S12(i,2) = 0.d0
         S12(i,ny) = 0.d0
      enddo
      do j=1,ny
         S12(nx,j) = 0.d0
         S12(2,j) = 0.d0
         S12(1,j) = 0.d0
      enddo

      return
      end

      SUBROUTINE makebccf(cc,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION cc(nx,ny),x,r1
      INTEGER dfetch
      INTEGER i,j
      i=dfetch("inner_nozzle_radius",r1)
      do i=1,ny
         cc(1,i) = cc(2,i)
         cc(nx,i) = cc(nx-1,i)
      enddo
      do i=1,nx
      x =(dfloat(i-1)/(nx-2))
       if (dabs(x) .le. r1) then
         cc(i,2) = 1.d0
       else
         cc(i,2)= 0.d0
       endif
         cc(i,1)= cc(i,2)
         cc(i,ny) = cc(i,ny-1)
      enddo

      return
      end

      SUBROUTINE makebcpf(cc,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION cc(nx,ny)
      INTEGER i,j

      do i=1,nx
         cc(i,1) = cc(i,2)
         cc(i,ny) = cc(i,ny-1)
      enddo
      do i=1,ny
         cc(1,i) = cc(2,i)
         cc(nx,i) = cc(nx-1,i)
      enddo

      return
      end

      SUBROUTINE makebcresf(res,nx,ny)
      include 'undefined.h'
      INTEGER nx,ny
      DOUBLE PRECISION res(nx,ny)
      INTEGER i,j

      do i=1,nx
         res(i,1) = 0.d0
         res(i,ny) = 0.d0
      enddo
      do i=1,ny
         res(1,i) = 0.d0
         res(nx,i) = 0.d0
      enddo
 
      return
      END

      SUBROUTINE makebcxvof(vof1,vof3,nx,ny)
      include 'undefined.h'
      INTEGER i,nx,ny
      DOUBLE PRECISION vof1(nx,ny),vof3(nx,ny)

      do i=2,ny-1
         vof3(1,i)=0.d0
         vof1(nx,i)=0.d0
      enddo

      return
      END

      SUBROUTINE makebczvof(vof1,vof3,nx,ny)
      include 'undefined.h'
      INTEGER i,nx,ny
      DOUBLE PRECISION vof1(nx,ny),vof3(nx,ny)
 
      do i=2,nx-1
         vof3(i,1)=0.0d0
         vof1(i,ny)=0.0d0
      enddo
 
      return
      END

c***********************************************************************
      SUBROUTINE makebcflux(flxm,flxp,flzm,flzp,nx,nz)
c***
      INCLUDE 'undefined.h'
      INTEGER i,k,nx,nz
      DOUBLE PRECISION flxm(nx,nz),flxp(nx,nz),flzm(nx,nz),flzp(nx,nz)
c***
      do k=2,nz-1
         flxp(1,k) = 0.d0
         flxm(nx,k) = 0.d0
c         flxp(1,k) = flxp(nx-1,k)
c         flxm(nx,k) = flxm(2,k)
      enddo 
      do i=2,nx-1
         flzp(i,1) = 0.d0
         flzm(i,nz) = 0.d0
      enddo 
c***
      return
      end

