      SUBROUTINE scbndy(tb,yb,rho,ijet)

      IMPLICIT NONE
      
      INCLUDE 'param.txt'
      INCLUDE 'param2.txt'

      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 alpha,beta,preexp,incrho,rhoi

      INTEGER jtmp
      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*************************************************************

      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)

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

c            rho(0,j,k) = 1.0/((incrho*alpha/(1.0-incrho*alpha))*1.0+rhoi)


c            rho(0,j,k) = 1.0/((incrho*alpha/(1.0-incrho*alpha))*tb(0,j,k
c     $           )+rhoi)
            yb(0,j,k)    = yb(1,j,k)

            tb(i1,j,k)   = tb(imax,j,k)
            rho(i1,j,k) = 1.0/((incrho*alpha/(1.0-incrho*alpha))*tb(i1,j
     $           ,k)+rhoi)
            yb(i1,j,k)   = yb(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,k1)  =  0.0 - tb(i,j,kmax) 
            yb(i,j,k1)  =  0.0 - yb(i,j,kmax)  ! y=zero

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

         END DO
      END DO


c.....conditions at lower jet exit.....

      DO i=0,ijet
         DO j=0,j1
               
c            tb(i,j,0)=0.0-tb(i,j,1) ! cold inlet

            tb(i,j,0)=2.0-tb(i,j,1) ! hot inlet

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

            yb(i,j,0)=2.0-yb(i,j,1) 

         END DO
      END DO

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

      DO i=ijet+1,i1-15
         DO j=0,j1
            
            tb(i,j,0)  = 0.0-tb(i,j,1) ! zero
            rho(i,j,0) = 1.0/((incrho*alpha/(1.0-incrho*alpha))
     $           *tb(i,j,0)+rhoi)

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

         END DO
      END DO

c....conditions at upper jet exit....

      DO i=0,ijet
         DO j=0,j1

            tb(i,j,k1) = 2.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)

         END DO
      END DO

c....conditions outside the upper jet exit...

      DO i=ijet+1,i1-15
         DO j=0,j1

            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)

         END DO
      END DO

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

      DO k=0,k1

         DO i=0,i1

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

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

         END DO

      END DO

      RETURN
      END
