C*********************************************************************
C*********************************************************************
C
      SUBROUTINE CLTURB
     O   ( SSTURB,
     I     U     , W     , TPOT  , QVAP  , QCLW  , QRAI  , CDTURB,
     I     IT1   , IT2   ,
     B     TPOT0 , QVAP0 , PPAI0 ,
     G     ZGRID , FDZ   , FDZM  , DX    ,
     C     GRAV  , CPDRY , ALATNT, RMVAP ,
     W     STABZ , EDFM1 , EDFM2 ,
     D     NX    , NZ    , NXB   , NZB   , NROT  )
C
C*******************************************************************
      IMPLICIT REAL*8 ( A-H, O-Z )
C*******************************************************************
C
      DIMENSION
     \  SSTURB(-2:NXB, -2:NZB)
C
      DIMENSION
     \  CDTURB(-2:NXB, -2:NZB, NROT),
     \    U   (-2:NXB, -2:NZB, NROT),
     \    W   (-2:NXB, -2:NZB, NROT),
     \    TPOT(-2:NXB, -2:NZB, NROT),
     \    QVAP(-2:NXB, -2:NZB, NROT),
     \    QCLW(-2:NXB, -2:NZB, NROT),
     \    QRAI(-2:NXB, -2:NZB, NROT)
C
      DIMENSION
     \    ZGRID(-2:NZB), QVAP0(-2:NZB), TPOT0(-2:NZB),
     \    PPAI0(-2:NZB),
     \    FDZM (-2:NZB), FDZ  (-2:NZB)
C
      DIMENSION
     \    EDFM1(-2:NXB,-2:NZB), EDFM2(-2:NXB,-2:NZB),
     \    STABZ(-2:NXB,-2:NZB)
C
C=================================================================
C     Velocity deformation.
C----------------------------------------------------------------------
C
C
      DO 11 IZ=0,NZ-1
      DO 11 IX=1,NX
          EDFM1(IX,IZ)
     \  = ( (  U(IX  ,IZ,IT2) - U(IX,IZ-1,IT2) ) / DX / FDZM(IZ)
     \     +(  W(IX  ,IZ,IT2) - W(IX-1,IZ,IT2) ) / DX     ) / 2.0D0
          EDFM2(IX,IZ)
     \  =    ( ( U(IX+1,IZ  ,IT2)-U(IX,IZ,IT2) ) / DX )**2
     \     + ( ( W(IX  ,IZ+1,IT2)-W(IX,IZ,IT2) ) / DX / FDZ(IZ) )**2
   11 CONTINUE
C
      DO 12 IX=1,NX
         EDFM1(IX,NZ) = 0.0
         EDFM2(IX,NZ) = 0.0
   12 CONTINUE
C
      CALL BOUND ( EDFM1, NX, NZ, NXB, NZB, 1,-1, 0, 1 )
      CALL BOUND ( EDFM2, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C-----------------------------------------------------------
C     Static stability.
C-----------------------------------------------------------
C
      DO 20 IX=1,NX
         STABZ(IX,NZ)=0.0
   20 CONTINUE
C
      DO 21 IZ = 0, NZ-1
      DO 21 IX = 1, NX
           STABZ(IX,IZ)
     \  =  GRAV / DX / FDZM(IZ) / ( (TPOT0(IZ)+TPOT0(IZ-1))/2.0D0 )
     \   *((  ( TPOT (IX,IZ,IT1) - TPOT (IX,IZ-1,IT1) )
     \      + ( TPOT0(   IZ   ) - TPOT0(   IZ-1   ) ) )
     \      + ALATNT/CPDRY / ((PPAI0(IZ)+PPAI0(IZ-1))/2.0D0)
     \        * ( 0.5D0 + SIGN ( 0.5D0,
     \          (QCLW(IX  ,IZ,IT1)+QCLW(IX  ,IZ-1,IT1) - 1.0D-4) ) )
     \       *( (QVAP (IX,IZ,IT1) - QVAP (IX,IZ-1,IT1) )
     \         +(QVAP0(   IZ   ) - QVAP0(   IZ-1   ) ) )          )
     \   + GRAV / DX / FDZM(IZ) * ( 1.0 / RMVAP - 1.0 )
     \       *( (QVAP (IX,IZ,IT1) - QVAP (IX,IZ-1,IT1) )
     \         +(QVAP0(   IZ   ) - QVAP0(   IZ-1   ) ) )
     \   + GRAV / DX / FDZM(IZ)
     \       * ( - QCLW(IX,IZ,IT1) + QCLW(IX,IZ-1,IT1)
     \           - QRAI(IX,IZ,IT1) + QRAI(IX,IZ-1,IT1)  )
   21 CONTINUE
C
      DO 22 IX = 1, NX
         STABZ(IX,NZ) = STABZ(IX,NZ-1)
   22 CONTINUE
C
      CALL BOUND(STABZ,NX, NZ, NXB,NZB,1,1,0,1)
C
C=====================================================================
C
      CMOM = 0.2D0
      CDIS = 0.2D0
      PRT  = 1.0D0
C
      DO 31 IZ = 0, NZ
C
      AMIX2 = ( DX * FDZ(IZ) ) ** 2
C
      DO 31 IX = 1, NX
         SSTURB(IX,IZ) =
     \      CMOM**2 * AMIX2 / 2.0D0 / PRT
     \    * ( - (  STABZ(IX  ,IZ)    + STABZ(IX  ,IZ+1) ) / 2.0D0
     \        + (  EDFM1(IX  ,IZ)**2 + EDFM1(IX  ,IZ+1)**2
     \           + EDFM1(IX+1,IZ)**2 + EDFM1(IX+1,IZ+1)**2 ) / 4.0D0
     \        +    EDFM2(IX  ,IZ)    )
     \    - CDIS * CDTURB(IX,IZ,IT2)**2 / 2.0D0 / CMOM / AMIX2
   31 CONTINUE
C
      CALL BOUND ( SSTURB, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C---------------------------------------------------------------
C
      RETURN
      E N D
