      SUBROUTINE bound(ub,vb,wb,tb,yb,rho,time,ijet,dt,zp,zw,advout,ru,rp)

      IMPLICIT NONE

      INCLUDE 'param.txt'
      INCLUDE 'param2.txt'

      REAL ub(0:i1,0:j1,0:k1)
      REAL vb(0:i1,0:j1,0:k1)
      REAL wb(0:i1,0:j1,0:k1)
      REAL tb(0:i1,0:j1,0:k1)
      REAL yb(0:i1,0:j1,0:k1)
      REAL rho(0:i1,0:j1,0:k1)

      REAL time,dt
      REAL zp(0:k1),zw(0:k1)
      REAL ru(0:i1),rp(0:i1)

      REAL alpha,beta,preexp,incrho,rhoi,wmean(imax)

      INTEGER advout

      INTEGER jtmp
      INTEGER iseed
      REAL*4 ran,ranf,tmp1,tmp2

      INTEGER ijet

      COMMON /heatrl/ alpha,beta,preexp,incrho,rhoi
c---------------------------------------------------------------

c*************************************************************
c
c     Set boundary conditions for i=0 and i=i1. For the boundary
c     conditions at i=0 we assume Neumann (no-flux) BC's. Mostly,
c     this B.C. is not used and hence, treated as dummies.

c*************************************************************
      wmean = 0
      do k=1,kmax
	 do j=1,jmax
	   do i=1,imax	
 	      wmean(i)=wmean(i)+wb(i,j,kmax)/(jmax*kmax)
           enddo
         enddo
      enddo

      DO k=1,kmax

         DO j=0,j1

            if ( j .le. jmax/2 ) then
               jtmp = j + jmax/2
            else
               jtmp = j - jmax/2
            endif

c *** Impermeability condition for the radial velocity component for
c *** i = imax:

            tb(0,j,k)    = tb(1,j,k)
            yb(0,j,k)    = yb(1,j,k)

            rho(0,j,k)   = rho(1,j,k)
            

            ub(0,j,k)    =   ub(1,j,k)
            vb(0,j,k)    =   vb(1,j,k)
            wb(0,j,k)    =   wb(1,j,k)



            tb(i1,j,k)   = tb(imax,j,k)
            yb(i1,j,k)   = yb(imax,j,k)

            rho(i1,j,k) = 1.0/((incrho*alpha/(1.0-incrho*alpha))*tb(i1,j
     $           ,k)+rhoi)

            ub(i1,j,k)   =  ub(imax,j,k)
            vb(i1,j,k)   =  vb(imax,j,k)
            wb(i1,j,k)   =  wb(imax,j,k)

         END DO

      END DO

c*************************************************************
c
c     Set boundary conditions for k=0 and k=k1.
c
c*************************************************************

      DO j=0,j1

         DO i=0,i1

            tb(i,j,0)  = 0.0 - tb(i,j,1)     ! t=zero
            yb(i,j,0)  = 0.0 - yb(i,j,1)     ! y=zero

            rho(i,j,0) = 1.0/((incrho*alpha/(1.0-incrho*alpha))*tb(i,j,0
     $           )+rhoi)

            ub(i,j,0)     = -ub(i,j,1)     ! zero velocity
            vb(i,j,0)     = -vb(i,j,1)     ! zero velocity
            wb(i,j,0)     = 0.0            ! zero velocity


            tb(i,j,k1)  =  0.0 - tb(i,j,kmax) 
            yb(i,j,k1)  = 0.0 - yb(i,j,kmax)  

            rho(i,j,k1) = 1.0/((incrho*alpha/(1.0-incrho*alpha))*tb(i,j
     $           ,k1)+rhoi)

c            ub(i,j,k1)    =  -ub(i,j,kmax) 
c            vb(i,j,k1)    =  -vb(i,j,kmax) 

            ub(i,j,k1)    =  ub(i,j,kmax) 
            vb(i,j,k1)    =  vb(i,j,kmax) 
            wb(i,j,kmax)  = wb(i,j,kmax)-dt*wmean(i)*(wb(i,j
     $              ,kmax)-wb(i,j,kmax-1))/(zw(kmax)-zw(kmax-1))


         END DO

      END DO

c.....conditions at lower and upper jet exits.....


c      CALL inlet(ub,vb,wb,yb,tb,rho,time,ru,rp,ijet)
cc      CALL jets(umb,vmb,wmb,ub,vb,wb,yb,tb,rho,time)

c.....conditions outside the lower jet exit....     

      DO i=i,i1
         DO j=0,j1
                  
            tb(i,j,0)  = 0.0 - tb(i,j,1) ! zero
            yb(i,j,0)  = 0.0 - yb(i,j,1) ! zero

            rho(i,j,0) = 1.0/((incrho*alpha/(1.0-incrho*alpha))
     $           *tb(i,j,0)+rhoi)

            wb(i,j,0)  = 0.

         END DO
      END DO

c.....conditions at the jet exit....     
      DO i=0,ijet
         DO j=0,j1

            tb(i,j,k1) =  tb(i,j,kmax)
            yb(i,j,k1) =  yb(i,j,kmax) 
            wb(i,j,0)  =  1.
            tb(i,j,0)  =  2 - tb(i,j,1)
            yb(i,j,0)  =  2 - yb(i,j,1)
            rho(i,j,0) = 1.0/((incrho*alpha/(1.0-incrho*alpha))
     $           *tb(i,j,0)+rhoi)
         ENDDO
      ENDDO
            

c*************************************************************************
c
c   boundary conditions on j=0 and j=j1 - periodic.
c
c
c*************************************************************************

      DO k=0,k1

         DO i=0,i1

            tb(i,0,k)  = tb(i,jmax,k)
            yb(i,0,k)  = yb(i,jmax,k)

            rho(i,0,k) = rho(i,jmax,k)

            ub(i,0,k) = ub(i,jmax,k)
            vb(i,0,k) = vb(i,jmax,k)
            wb(i,0,k) = wb(i,jmax,k)



            tb(i,j1,k) = tb(i,1,k)
            yb(i,j1,k) = yb(i,1,k)

            rho(i,j1,k) = rho(i,1,k)

            ub(i,j1,k) = ub(i,1,k)
            vb(i,j1,k) = vb(i,1,k)
            wb(i,j1,k) = wb(i,1,k)

         END DO

      END DO

      RETURN
      END
