C*********************************************************************
C*********************************************************************
C
      SUBROUTINE CLVISC
     I   (   U,   V,   W,
     I     CD1, CD2, EDFM1,
     I     ICL,
     W     DIV,
     O     DUVIS, DVVIS, DWVIS,
     B     ZGRID, DENS0, FDZ, FDZM, CD00,
     G     DX,
     D     NX, NZ, NXB, NZB, NROT  )
C
C*********************************************************************
      IMPLICIT REAL*8 ( A-H, O-Z )
C*********************************************************************
C
      DIMENSION
     \  U   ( -2:NXB, -2:NZB, NROT),
     \  V   ( -2:NXB, -2:NZB, NROT),
     \  W   ( -2:NXB, -2:NZB, NROT)
C
      DIMENSION
     \  CD1 ( -2:NXB, -2:NZB ), CD2 ( -2:NXB, -2:NZB ),
     \  EDFM1( -2:NXB, -2:NZB )
C
      DIMENSION
     \   DUVIS( -2:NXB, -2:NZB ),
     \   DVVIS( -2:NXB, -2:NZB ),
     \   DWVIS( -2:NXB, -2:NZB ),
     \   DIV   ( -2:NXB, -2:NZB )
C
      DIMENSION
     \   DENS0( -2:NZB  ),
     \   ZGRID( -2:NZB  ),
     \   FDZ  ( -2:NZB  ),
     \   FDZM ( -2:NZB  ),
     \   CD00 ( -2:NZB  )
C
C====================================================================
C
      CD0MAX = 1.0D+3
      ZC1 = 15000.
      ZC2 = 20000.
C
      DO 10 IZ = 0, NZ-1
         CD00 (IZ) = CD0MAX
     \     * MAX( 0.0D0 , MIN( 1.0D0 ,
     \         ( ( ZGRID(IZ) - ZC1 ) / ( ZC2 - ZC1 ) )   ))
   10 CONTINUE
C
      DO 11 IZ = 0, NZ-1
      DO 11 IX = 1, NX
         DIV(IX,IZ) =
     \      ( U(IX+1,IZ,ICL) - U(IX,IZ,ICL) ) / DX
     \    + ( W(IX,IZ+1,ICL) - W(IX,IZ,ICL) ) / DX / FDZ(IZ)
   11 CONTINUE
C
      CALL BOUND ( DIV, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
      DO 21 IZ = 0, NZ-1
      DO 21 IX = 1, NX
         DUVIS(IX,IZ) =
     \  ( (  ( CD00(IZ) + CD2(IX,IZ) )
     \          * ( U(IX+1,IZ,ICL) - U(IX,IZ,ICL) )
     \     - ( CD00(IZ) + CD2(IX-1,IZ) )
     \          * ( U(IX,IZ,ICL) - U(IX-1,IZ,ICL) )  )
     \  + (  ( DENS0(IZ+1) + DENS0(IZ) ) / 2.0D0
     \     * ( ( CD00(IZ+1) + CD00(IZ) ) / 2.0D0 + CD1(IX,IZ+1) )
     \     * ( U(IX,IZ+1,ICL) - U(IX,IZ,ICL) ) / FDZM(IZ+1)
     \    -  ( DENS0(IZ) + DENS0(IZ-1) ) / 2.0D0
     \     * ( ( CD00(IZ) + CD00(IZ-1) ) / 2.0D0 + CD1(IX,IZ) )
     \     * ( U(IX,IZ,ICL) - U(IX,IZ-1,ICL) ) / FDZM(IZ  ) )
     \                     / FDZ(IZ) / DENS0(IZ)  )  / DX**2
C        DUVIS(IX,IZ) =
C    \  ( (  ( CD00(IZ) + CD2(IX,IZ) )
C    \          * (   2.0D0 * ( U(IX+1,IZ,ICL) - U(IX,IZ,ICL) ) / DX
C    \              - 2.0/3.0 * DIV(IX,IZ)                      )
C    \     - ( CD00(IZ) + CD2(IX-1,IZ) )
C    \          * (   2.0D0 * ( U(IX,IZ,ICL) - U(IX-1,IZ,ICL) ) / DX
C    \              - 2.0/3.0 * DIV(IX-1,IZ)                    ) ) /
C    \  + (  ( DENS0(IZ+1) + DENS0(IZ) ) / 2.0D0
C    \     * ( ( CD00(IZ+1) + CD00(IZ) ) / 2.0D0 + CD1(IX,IZ+1) )
C    \                    * EDFM1(IX,IZ+1)
C    \    -  ( DENS0(IZ) + DENS0(IZ-1) ) / 2.0D0
C    \     * ( ( CD00(IZ) + CD00(IZ-1) ) / 2.0D0 + CD1(IX,IZ) )
C    \                    * EDFM1(IX,IZ)                  ) / DX
C    \                     / FDZ(IZ) / DENS0(IZ)               )
   21 CONTINUE
C
      CALL BOUND ( DUVIS, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C
      DO 22 IZ = 0, NZ-1
      DO 22 IX = 1, NX
         DVVIS(IX,IZ) =  1.0D0 / DX**2
     \ * (  ( (CD1(IX+1,IZ) + CD1(IX+1,IZ+1)) / 2.0D0 + CD00(IZ) )
     \    * ( V(IX+1,IZ,ICL) - V(IX,IZ,ICL) )
     \  -   ( (CD1(IX,IZ)   + CD1(IX,IZ+1 ) ) / 2.0D0 + CD00(IZ) )
     \    * ( V(IX,IZ,ICL) - V(IX-1,IZ,ICL) )
     \  + ( ( DENS0(IZ) + DENS0(IZ+1) ) / 2.0D0
     \     *(  ( CD1(IX+1,IZ+1) + CD1(IX,IZ+1))
     \        +( CD00(IZ+1) + CD00(IZ)  )         ) / 2.0D0
     \     * ( V(IX,IZ+1,ICL) - V(IX,IZ,ICL) ) / FDZM(IZ+1)
     \    - ( DENS0(IZ) + DENS0(IZ-1) ) / 2.0D0
     \     *(  ( CD1(IX+1,IZ  ) + CD1(IX,IZ  ))
     \        +( CD00(IZ) + CD00(IZ-1) )          ) / 2.0D0
     \     * ( V(IX,IZ,ICL) - V(IX,IZ-1,ICL) ) / FDZM(IZ)
     \            ) / FDZ(IZ) / DENS0(IZ)                    )
   22 CONTINUE
C
      CALL BOUND( DVVIS, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C
      DO 23 IZ = 0, NZ
      DO 23 IX = 1, NX
          DWVIS(IX,IZ) =
     \     (    CD1(IX+1,IZ) * ( W(IX+1,IZ,ICL) - W(IX,IZ,ICL) )
     \       -  CD1(IX  ,IZ) * ( W(IX,IZ,ICL) - W(IX-1,IZ,ICL) )
     \   +   ( CD00(IZ) + CD00(IZ-1) ) / 2.0D0
     \     * ( W(IX+1,IZ,ICL) - 2.0D0*W(IX,IZ,ICL) + W(IX-1,IZ,ICL) )
     \   + (   DENS0(IZ) * ( CD2(IX,IZ) + CD00(IZ) )
     \       * ( W(IX,IZ+1,ICL) - W(IX,IZ,ICL) ) / FDZ(IZ)
     \       - DENS0(IZ-1) * ( CD2(IX,IZ-1) + CD00(IZ-1) )
     \       * ( W(IX,IZ,ICL) - W(IX,IZ-1,ICL) ) / FDZ(IZ-1) )
     \    / ( (DENS0(IZ)+DENS0(IZ-1))/2.0D0 ) / FDZM(IZ)     ) / DX**2
C         DWVIS(IX,IZ) =
C    \     (    CD1(IX+1,IZ) * EDFM1(IX+1,IZ)
C    \       -  CD1(IX  ,IZ) * EDFM1(IX  ,IZ) ) / DX
C    \   +   ( CD00(IZ) + CD00(IZ-1) ) / 2.0D0
C    \     * ( EDFM1(IX+1,IZ) - EDFM1(IX,IZ) ) / DX
C    \   + (   DENS0(IZ) * ( CD2(IX,IZ) + CD00(IZ) )
C    \       * (   2.0D0 * ( W(IX,IZ+1,ICL) - W(IX,IZ,ICL) ) / DX
C    \           - 2.0/3.0 * DIV(IX,IZ)  )
C    \       - DENS0(IZ-1) * ( CD2(IX,IZ-1) + CD00(IZ-1) )
C    \       * (   2.0D0 * ( W(IX,IZ,ICL) - W(IX,IZ-1,ICL) ) / DX
C    \           - 2.0/3.0 * DIV(IX,IZ-1)  )                 )
C    \    / ( (DENS0(IZ)+DENS0(IZ-1))/2.0D0 )  / FDZM(IZ)
   23 CONTINUE
C
C     ALPHA = 1.0D-2
C
C     DO 24 IX = 1, NX
C        DWVIS(IX,NZ-1) = DWVIS(IX,NZ-1) - ALPHA * W(IX,NZ-1,ICL)
C  24 CONTINUE
C
      CALL BOUND ( DWVIS, NX, NZ, NXB, NZB, 1, -1, 0, 1 )
C
C====================================================================
C
      RETURN
      E N D
