* PACKAGE P2CCLL    
*
***********************************************************************
      SUBROUTINE CCLLEV     
     O         ( GDTC  , GDQC  , KBTM  , KTOP  ,
     I           GDT   , GDQ   , GDP   , GDZ   ,
     I           IREUSE,
     W           ADIAF                           )
*
*"       雲底下端 = LCL
*"       エントレインメント
*
*   [PARAM] 





* PACKAGE  ZXT10L32 !" 格子点数，波数 X-Z
*
      INTEGER     IMAX           !" 経度格子数
      INTEGER     JMAX           !" 緯度格子数
      INTEGER     IDIM           !" 経度配列サイズ
      INTEGER     JDIM           !" 緯度配列サイズ
      INTEGER     MMAX           !" 最大東西波数
      INTEGER     LMAX           !" 最大南北波数
      INTEGER     NMAX           !" 最大全波数
      INTEGER     MINT           !" 東西波数間隔
      INTEGER     KMAX           !" 鉛直層数
      INTEGER     KDIM           !" 鉛直配列サイズ
      INTEGER     NMDIM          !" 水平波数成分数
      INTEGER     JMXHF          !" JMAX/2+1
      INTEGER     NWORK          !" ワーク領域大きさ
*
      PARAMETER ( IMAX  = 32,     JMAX =  1,
     &            IDIM  = IMAX+1, JDIM = JMAX,
     &            MMAX  = 10,     LMAX =  0,    NMAX = MMAX, MINT = 1,
     &            KMAX  = 32,     KDIM = KMAX,
     &            NMDIM = (MMAX/MINT+1)*(2*(NMAX+1)-MMAX)
     &                  - (NMAX-LMAX)/MINT*(NMAX-LMAX+1),
     &            JMXHF = JMAX/2+1,
     &            NWORK = IDIM*JDIM*KMAX*12                       )

* PACKAGE ZHDIM  !" 標準物理定数(common)  
*
*"  [HIS]  90/08/31(numaguti)
*
      ReaL*8       ER             !" 地球半径
      ReaL*8       GRAV                        !" 重力加速度
      ReaL*8       CP                          !" 大気定圧比熱
      ReaL*8       RAIR                        !" 大気気体定数
      ReaL*8       EL                          !" 水の凝結の潜熱
      ReaL*8       CPVAP                       !" 水蒸気定圧比熱
      ReaL*8       RVAP                        !" 水蒸気気体定数
      ReaL*8       DH2O                        !" 水の密度
      ReaL*8       EPSV                        !" 水と大気の分子量比ε
      ReaL*8       ES0                         !" ０℃飽和蒸気圧：Ｐａ
      ReaL*8       STB                         !" Stefan-Boltzman
      ReaL*8       FKARM                       !" Karman 定数
      ReaL*8       EPSVT                       !" 1/ε-1
*
      COMMON   /COMCON/
     &          ER    , GRAV  ,
     &          CP    , RAIR  ,
     &          EL    , CPVAP , RVAP  , DH2O  , EPSV  , ES0   ,
     &          STB   , FKARM , EPSVT

*
*   [OUTPUT] 
      ReaL*8       GDTC  ( IDIM*JDIM, KMAX )   
      ReaL*8       GDQC  ( IDIM*JDIM, KMAX )   
      INTEGER    KBTM  ( IDIM*JDIM )         
      INTEGER    KTOP  ( IDIM*JDIM )         
*
*   [INPUT] 
      ReaL*8       GDT   ( IDIM*JDIM, KMAX )   
      ReaL*8       GDQ   ( IDIM*JDIM, KMAX )   
      ReaL*8       GDP   ( IDIM*JDIM, KMAX )   
      ReaL*8       GDZ   ( IDIM*JDIM, KMAX )   
*
      INTEGER    IREUSE                      
*
*   [WORK] 
      ReaL*8       ADIAF ( IDIM*JDIM, KMAX )   
*
*   [INTERNAL PARM] 
      INTEGER    KBMIN
      INTEGER    KBMAX
      DATA       KBMIN  / 1    /
      DATA       KBMAX  / KMAX /
      ReaL*8       EPST                        
      DATA       EPST  / 0.0   /
      ReaL*8       ENTR                        
      DATA       ENTR  / 0.0   /
      NAMELIST  /NMCCLL/  KBMIN, KBMAX, EPST, ENTR
*
      LOGICAL    OFIRST, OSET
      DATA       OFIRST, OSET  / .TRUE., .FALSE. /
      SAVE
*
*   [INTERNAL WORK] 
      INTEGER    KBTMZ ( IDIM*JDIM )         
      INTEGER    KTOPZ ( IDIM*JDIM )         
*     
      INTEGER    IJ, K, KUNDEF
      INTEGER    IFPAR, JFPAR
      ReaL*8       AKAPPA, ELF
      ReaL*8       RENTR, QSAT, DQSAT, QEXE, DT, DQ
*




* PACKAGE ZQSAT     !" Clausis-Clapeyron の式(文関数)
*
*"    QSAT:  saturation water vapour mixing ratio
*"    DQSAT: d(QSAT)/d(T)
*
      ReaL*8       FQSAT
      ReaL*8       FDQSAT
      ReaL*8       T, P, QS
*
c$$$      FQSAT ( T,P )    = EPSV * ES0
c$$$     &                   * EXP( EL/RVAP *( 1./273. - 1./T ) )
c$$$     &                   / P
      FQSAT ( T,P )    = EPSV
     &                   * 1.4E+11* EXP( - 5250.87/T )
     &                   / P
*"         <---------  Nakajima version
      FDQSAT( T,QS )   = EL * QS / ( RVAP * T*T )

*
      IF ( OFIRST ) THEN
         WRITE ( 6,* ) ' CONVECTIVE CLOUD LAYER P2.3 DATE=92/02/07'
         OFIRST = .FALSE.
*
         CALL   REWNML ( IFPAR , JFPAR )
         READ   ( IFPAR, NMCCLL, END=90 )
   90    WRITE  ( JFPAR, NMCCLL )
*
         KUNDEF = -1
      ENDIF
*
*"         < 0. IREUSE = 1: 前値再利用 (shallow convection) >
*
      IF ( ( IREUSE .EQ. 1 ) .AND. OSET ) THEN
         DO 100 IJ = 1, IDIM*JDIM
            KBTM ( IJ ) = KBTMZ ( IJ )
            KTOP ( IJ ) = KTOPZ ( IJ )
  100    CONTINUE
         OSET = .FALSE.
      ENDIF
*
      IF ( IREUSE .NE. 1 ) OSET = .TRUE.
*
*
*"         < 1. ファクター計算 >
*
      AKAPPA = RAIR / CP
      DO 1100 K = 1, KMAX-1
         DO 1100 IJ = 1, IDIM*JDIM
            ADIAF( IJ,K ) = ( GDP( IJ,K+1 )/GDP( IJ,K ) ) **AKAPPA
 1100 CONTINUE
*
      ELF    = EL / CP
*
*"         < 2. 初期設定 >
*
      DO 2100 IJ = 1, IDIM*JDIM
         KBTM  ( IJ ) = KUNDEF
         KTOP  ( IJ ) = KUNDEF
 2100 CONTINUE
*
      DO 2200 IJ = 1, IDIM*JDIM
         GDTC( IJ,1 ) = GDT( IJ,1 )
         GDQC( IJ,1 ) = GDQ( IJ,1 )
 2200 CONTINUE
*
*"         < 3. 雲底・雲頂 >
*
      DO 3100 K = 2, KMAX
         DO 3110 IJ = 1, IDIM*JDIM
*
            IF ( KBTM( IJ ) .EQ. KUNDEF ) THEN
*
*"               < 3.1 持ち上げ凝結高度 >
*
               GDTC( IJ,K ) = GDTC( IJ,K-1 ) * ADIAF( IJ,K-1 )
               GDQC( IJ,K ) = GDQC( IJ,K-1 )
               QSAT         = FQSAT( GDTC( IJ,K ),GDP( IJ,K ) )
*
               IF (  ( GDQC( IJ,K ) .GE. QSAT ).AND.
     &               ( K .GE. KBMIN )          .AND.
     &               ( K .LE. KBMAX )                ) THEN
                  KBTM( IJ )   = K
                  GDTC( IJ,K ) = GDT( IJ,K )
                  GDQC( IJ,K ) = FQSAT( GDTC( IJ,K ), GDP( IJ,K ) )
               ENDIF
*
            ELSE IF ( KTOP( IJ ) .EQ. KUNDEF )  THEN
*
*"               < 3.2 湿潤上昇 >
*
               RENTR        = ENTR * ( GDZ( IJ,K ) - GDZ( IJ,K-1 ) )
               RENTR        = MIN( RENTR, 1. )
*
               GDTC( IJ,K ) = ( 1.-RENTR ) 
     &                         * GDTC( IJ,K-1 ) * ADIAF( IJ,K-1 )
     &                       + RENTR 
     &                         * ( GDT( IJ,K-1 ) + GDT( IJ,K ) )/2.  
               GDQC( IJ,K ) = ( 1.  -RENTR ) * GDQC( IJ,K-1 )
     &                       + RENTR 
     &                         * ( GDQ( IJ,K-1 ) + GDQ( IJ,K ) )/2.  
*
               QSAT         = FQSAT( GDTC( IJ,K ),GDP( IJ,K ) )
*
               IF ( GDQC( IJ,K ) .GT. QSAT ) THEN
                  DQSAT        = FDQSAT( GDTC( IJ,K ), QSAT  )
                  QEXE         = GDQC ( IJ,K ) - QSAT
                  DT           = ELF * QEXE / ( 1. + ELF * DQSAT )
                  DQ           = DQSAT * DT
*
                  GDTC( IJ,K ) = GDTC( IJ,K ) + DT
                  GDQC( IJ,K ) = QSAT         + DQ
               ENDIF
*
*"               < 3.3 安定度チェック >
*
               IF ( GDTC( IJ,K ) .LT. GDT( IJ,K )-EPST ) THEN
*
                  IF ( K .EQ. KBTM( IJ )+1  ) THEN
                     KBTM( IJ   ) = K
                     GDTC( IJ,K ) = GDT( IJ,K )
                     GDQC( IJ,K ) = FQSAT( GDTC( IJ,K ), GDP( IJ,K ) )
                  ELSE
                     KTOP( IJ ) = K - 1
                  ENDIF
*
               ELSE IF ( K .EQ. KMAX ) THEN
                  KTOP( IJ ) = K
               ENDIF
*
            ENDIF
*
 3110    CONTINUE
 3100 CONTINUE
*
      DO 3300 IJ = 1, IDIM*JDIM
         IF ( KTOP( IJ ) .EQ. KUNDEF ) THEN
            KBTM( IJ ) = KUNDEF
         ENDIF
*
         KBTMZ ( IJ ) = KBTM ( IJ )
         KTOPZ ( IJ ) = KTOP ( IJ )
 3300 CONTINUE
*
      RETURN
      END
