      SUBROUTINE viscos(ib,ie,jb,je,kb,ke)

      IMPLICIT NONE

      INCLUDE 'param.txt'
      INCLUDE 'common.txt'

      INTEGER ib,ie,jb,je,kb,ke

      REAL tmp
c--------------------------------------------------------------------

      DO k=kb,ke
         DO j=jb,je
            DO i=ib,ie

               tmp=((tnew(i,j,k)-tnew(i,j,k-1))/(zp(k)-zp(k-1)))**2
     $              +((tnew(i,j,k)-tnew(i-1,j,k))/(rp(i)-rp(i-1)))**2
     $              +((tnew(i,j,k)-tnew(i,j-1,k))/(rp(i)*dphi))**2

               tmp=SQRT(tmp)

               munew(i,j,k)  = 1.0/(Re) + rhonew(i,j,k)*mut(i,j,k)
               lmbnew(i,j,k) = 1.0/(Re*Pr) + rhonew(i,j,k)*mut(i,j,k)
     $              *Cl*tmp
               difnew(i,j,k) = 1.0/(Re*Pr*Le) + rhonew(i,j,k)*mut(i,j,k)
     $              *Cl*tmp


               if(munew(i,j,k) .lt. 0.0)then
                  munew(i,j,k) = 1.0/Re
               end if


            END DO
         END DO
      END DO

c*************************************************************
c
c     Set boundary boundary condition for K-closure factors.
c
c*************************************************************

      DO k=1,kmax

c *** Theta-direction *** Periodic direction ***

         DO i=1,imax

            munew(i,0,k)  = munew(i,jmax,k)
            munew(i,j1,k) = munew(i,1,k)

            lmbnew(i,0,k)  = lmbnew(i,jmax,k)
            lmbnew(i,j1,k) = lmbnew(i,1,k)

            difnew(i,0,k)  = difnew(i,jmax,k)
            difnew(i,j1,k) = difnew(i,1,k)

         END DO

c *** R-direction ***

         DO j=0,j1

            munew(0,j,k)  =   munew(1,j,k)
c            munew(i1,j,k) = - munew(imax,j,k)
            munew(i1,j,k) = munew(imax,j,k)


            lmbnew(0,j,k)  =   lmbnew(1,j,k)
            lmbnew(i1,j,k) = - lmbnew(imax,j,k)

            difnew(0,j,k)  =   difnew(1,j,k)
            difnew(i1,j,k) = - difnew(imax,j,k)

         END DO

      END DO

c *** Z-direction ***

      DO j=0,j1
        DO i=0,i1

           munew(i,j,0)  = munew(i,j,kmax)
           munew(i,j,k1) = munew(i,j,1)

           lmbnew(i,j,0)  = lmbnew(i,j,kmax)
           lmbnew(i,j,k1) = lmbnew(i,j,1)

           difnew(i,j,0)  = difnew(i,j,kmax)
           difnew(i,j,k1) = difnew(i,j,1)

        END DO
      END DO

      RETURN
      END
