C***********************************************************
C condensation scheme follows Klemp & Wilhelmson (1978)
C***********************************************************
C
      SUBROUTINE CLPHYS
     I   (  TPOT1, QVAP1, QCLW1, QRAI1,
     I      DTI  , ITMIC, IOUT ,
     M      TPOT , QVAP , QCLW , QRAI ,
     O      PCCON, PVEVC, PVEVR,
     O      PRMIC, PCONV, PCOAG,
     B      DENS0, PPAI0, QVAP0, pres0, tpot0,
     C      CPDRY, ALATNT, rmvap,
     D      NX, NZ, NXB, NZB, NROT  )
C
C***********************************************************
      IMPLICIT REAL*8 ( A-H, O-Z )
C***********************************************************
C
      DIMENSION
     \   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
     \   TPOT1 ( -2:NXB, -2:NZB ),
     \   QVAP1 ( -2:NXB, -2:NZB ),
     \   QCLW1 ( -2:NXB, -2:NZB ),
     \   QRAI1 ( -2:NXB, -2:NZB )
C
      DIMENSION
     \   PCCON ( -2:NXB, -2:NZB ),
     \   PVEVC ( -2:NXB, -2:NZB ),
     \   PVEVR ( -2:NXB, -2:NZB ),
     \   PRMIC ( -2:NXB, -2:NZB ),
     \   PCONV ( -2:NXB, -2:NZB ),
     \   PCOAG ( -2:NXB, -2:NZB )
C
      DIMENSION
     \   DENS0( -2:NZB ),
     \   PPAI0( -2:NZB ),
     \   PRES0( -2:NZB ),
     \   TPOT0( -2:NZB ),
     \   QVAP0( -2:NZB )
C
C
C===========================================================
C
C --- Parameters for microphysical calculations.
C
      RCONV = 1.0 D-2
      DCRIT = 0.0 D-3
C
C
      RCOAG = 2.2 D+0 
      REVAP = 4.85D-2
C
C --- Conversion of cloud into rain. ---
C
      DO 11 IZ = 0, NZ
      DO 11 IX = 1, NX
C
         PCONV(IX,IZ) =
     \      RCONV * DTI
     \    * MAX( 0.0D0, ( QCLW(IX,IZ,ITMIC) - DCRIT/DENS0(IZ) ) )
C
C
         PCOAG(IX,IZ) =
     \      RCOAG * DTI
     \    *   MAX( 0.0D0, QCLW(IX,IZ,ITMIC) )
     \    * ( MAX( 0.0D0, QRAI(IX,IZ,ITMIC) ) * DENS0(IZ) ) ** 0.875
C
         PRMIC(IX,IZ) =
     \      MIN (   MAX( 0.0D0, QCLW1(IX,IZ) ),
     \             ( PCONV(IX,IZ) + PCOAG(IX,IZ) )  )
C
         QCLW1(IX,IZ) = QCLW1(IX,IZ) - PRMIC(IX,IZ)
         QRAI1(IX,IZ) = QRAI1(IX,IZ) + PRMIC(IX,IZ)
C
   11 CONTINUE
C
C --- Saturation adjustment including the evaporation of rain.
C
      dummy = evsat_p ( e0, c1, t0, t1 )

c     write (0,*) 'clphys_x e0 c1 t0 t1', e0, c1, t0, t1

      DO 21 IZ = 0, NZ-1
      DO 21 IX = 1, NX
C
C *** new scheme follows Qsat linearization in Klemp & WIlmhelmson
C
         tdegc1 = ppai0(iz) * ( tpot0(iz) + tpot1(ix,iz) ) - t0

         qvsat1 =   e0 / pres0(iz) * rmvap 
     &       * EXP( c1 * tdegc1 / (tdegc1 + (t0 - t1 ) ) ) 

         gamlam =   alatnt * qvsat1 * c1 * (t0 - t1 )
     &            / cpdry / ( tdegc1 + (t0 - t1 ) )**2

         qchng =   
     &        ( (qvap1(ix,iz) + qvap0(iz) ) - qvsat1 )
     &        / (1.0d0 + gamlam )

*
* end of new scheme 
*

         PCCON(IX,IZ) =
     \      MAX( 0.0D0 ,
     \           MIN(  ( QVAP1(IX,IZ) + QVAP0(IZ) ),
     \                   QCHNG                       ) )
         PVEVC(IX,IZ) =
     \      MAX( 0.0D0 , MIN( QCLW1(IX,IZ),  - QCHNG )   )
C
C
         QUSAT = MAX( 0.0D0, ( - QCHNG - PVEVC(IX,IZ) ) )
C
         EVRAI =   REVAP * QUSAT * DTI
     \     * ( MAX( 0.0D0, QRAI(IX,IZ,ITMIC) ) * DENS0(IZ) ) ** 0.65
         PVEVR(IX,IZ) = MIN(             EVRAI, 
     &                       MAX( 0.0D0, QRAI1(IX,IZ) )  )
C
C
         QVAP(IX,IZ,IOUT) = QVAP1(IX,IZ)
     \           - PCCON(IX,IZ) + PVEVC(IX,IZ) + PVEVR(IX,IZ)
         QCLW(IX,IZ,IOUT) = QCLW1(IX,IZ)
     \           + PCCON(IX,IZ) - PVEVC(IX,IZ)
         QRAI(IX,IZ,IOUT) = QRAI1(IX,IZ)       - PVEVR(IX,IZ)
         TPOT(IX,IZ,IOUT) = TPOT1(IX,IZ)
     \        +    ALATNT / CPDRY / PPAI0(IZ)
     \         * ( PCCON(IX,IZ) - PVEVC(IX,IZ) - PVEVR(IX,IZ) )
C
   21 CONTINUE
C
C-----------------------------------------------------------
C
      CALL BOUND ( TPOT, NX, NZ, NXB, NZB, NROT, 1, 1, IOUT )
      CALL BOUND ( QVAP, NX, NZ, NXB, NZB, NROT, 1, 1, IOUT )
      CALL BOUND ( QCLW, NX, NZ, NXB, NZB, NROT, 1, 1, IOUT )
      CALL BOUND ( QRAI, NX, NZ, NXB, NZB, NROT, 1, 1, IOUT )
C
      DO 31 IX = -2, NXB
c        TPOT(IX,NZ-1,IOUT) = 0.0D0
c        TPOT(IX,NZ  ,IOUT) = 0.0D0
c        TPOT(IX,NZ-2,IOUT) = 0.0D0
c        TPOT(IX,NZ+1,IOUT) = 0.0D0
c        QVAP(IX,NZ-1,IOUT) = 0.0D0
c        QVAP(IX,NZ  ,IOUT) = 0.0D0
c        QVAP(IX,NZ-2,IOUT) = 0.0D0
c        QVAP(IX,NZ+1,IOUT) = 0.0D0
c        QCLW(IX,NZ-1,IOUT) = 0.0D0
c        QCLW(IX,NZ  ,IOUT) = 0.0D0
c        QCLW(IX,NZ-2,IOUT) = 0.0D0
c        QCLW(IX,NZ+1,IOUT) = 0.0D0
c        QRAI(IX,NZ-1,IOUT) = 0.0D0
c        QRAI(IX,NZ  ,IOUT) = 0.0D0
c        QRAI(IX,NZ-2,IOUT) = 0.0D0
c        QRAI(IX,NZ+1,IOUT) = 0.0D0
   31 CONTINUE
C
      CALL BOUND ( PRMIC, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
      CALL BOUND ( PCONV, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
      CALL BOUND ( PCOAG, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
      CALL BOUND ( PCCON, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
      CALL BOUND ( PVEVC, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
      CALL BOUND ( PVEVR, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C===========================================================
C
      RETURN
      E N D
