** PACKAGE P2MADJ  !"  物理過程(2) 湿潤対流調節
*
*"  [HIS]  90/05/19(numaguti)
*"         91/08/09(numaguti) バグ除去
*"         94/12/17(ishiwatari) gdq > c*gdqs ? の形に
**********************************************************************
      SUBROUTINE CUMADJ       !" 湿潤対流調節
     M         ( GDT   , GDQ   ,
     M           GQPR  , GQCNV , GSCNV ,
     I           GDP   , GDPM  , GDZ   , DELTA  )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#include        "zccom.F"                    !" 標準物理定数
#endif
*
*   [MODIFY] 
      REAL       GDT   ( IDIM*JDIM, KMAX )   !" 温度  Ｔ
      REAL       GDQ   ( IDIM*JDIM, KMAX )   !" 比湿  ｑ
      REAL       GQPR  ( IDIM*JDIM       )   !" 降水量
      REAL       GQCNV ( IDIM*JDIM, KMAX )   !" 凝結調節加熱率
      REAL       GSCNV ( IDIM*JDIM, KMAX )   !" 凝結調節比湿変化
*
*   [INPUT] 
      REAL       GDP    ( IDIM*JDIM, KMAX )  !" 気圧  Ｐ
      REAL       GDPM   ( IDIM*JDIM, KMAX+1 )!" 気圧（半整数）
      REAL       GDZ    ( IDIM*JDIM, KMAX )  !" 高度
      REAL       DELTA                       !" 2Δt
*
*   [INTERNAL PARM] 
      INTEGER    ITRMAX
      PARAMETER (ITRMAX=10)
      REAL       EPST ( ITRMAX )             !" 不安定の許容誤差
      REAL       RCRIT                       !" 臨界相対湿度
      DATA       EPST / 0.01  , 0.02  , 0.02  , 0.05  , 0.05  ,
     &                  0.10  , 0.10  , 0.20  , 0.20  , 0.40    /
      DATA       RCRIT/ 0.990 /
*
      NAMELIST  /NMMADJ/ EPST, RCRIT
*
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
      SAVE
*
*   [INTERNAL WORK] 
      COMMON    /COMWRK/
     &           ADJDA , ADJDB ,
     &           GDTB  , GDQB  ,
     &           DELP  , DPDP  , DPPMK ,
     &           GDQS  
      REAL       ADJDA ( IDIM*JDIM )         !" 今回調節されたか否か
      REAL       ADJDB ( IDIM*JDIM )         !" 前回調節されたか否か
      REAL       GDTB  ( IDIM*JDIM, KMAX )   !" 調節前の温度
      REAL       GDQB  ( IDIM*JDIM, KMAX )   !" 調節前の比湿
      REAL       DELP  ( IDIM*JDIM, KMAX )   !" ΔＰ
      REAL       DPDP  ( IDIM*JDIM, KMAX )   !" ファクター
      REAL       DPPMK ( IDIM*JDIM, KMAX )   !" ファクター
      REAL       GDQS  ( IDIM*JDIM, KMAX )   !" 飽和比湿
*
      INTEGER    IJ, K, ITR, I, J
      INTEGER    IFPAR, JFPAR
      REAL       ELF, AKAPPA
      REAL       ST, QEXE, DQSAT, DQSATB, DT, STEXE
      REAL       DHDT, DHDTB, DTB, ADJPTS
*      
*   [INTERNAL FUNC]
#ifdef SYS_IBMS
      INCLUDE   (ZQSAT)                      !" 飽和水蒸気圧
#else
#include        "zqsat.F"                    !" 飽和水蒸気圧
#endif
*
*   [ONCE] 
      IF ( OFIRST ) THEN
         WRITE ( 6,* ) ' MOIST CONVENCIVE ADJUSTMENT P2.1 DATE=91/01/12'
         OFIRST = .FALSE.
*
         CALL   REWNML ( IFPAR , JFPAR )
         READ   ( IFPAR, NMMADJ, END=190 )
  190    WRITE  ( JFPAR, NMMADJ )
*
      ENDIF
*
*"         < 1. 調節前温度, 比湿のセーブ >
*
      DO 1000 K = 1, KMAX
         DO 1000 IJ = 1, IDIM*JDIM
            GDTB ( IJ,K ) = GDT ( IJ,K )
            GDQB ( IJ,K ) = GDQ ( IJ,K )
 1000 CONTINUE
*
*"         < 2. ファクター計算 >
*
      ELF    = EL/CP
      AKAPPA = RAIR / CP
*
      DO 2100 K = 1, KMAX
         DO 2100 IJ = 1, IDIM*JDIM
            DELP  ( IJ,K ) = GDPM( IJ,K   ) - GDPM( IJ,K+1 )
 2100 CONTINUE
*
      DO 2200 K = 2, KMAX
         DO 2200 IJ = 1, IDIM*JDIM
            DPDP  ( IJ,K ) = DELP( IJ,K ) / DELP( IJ,K-1 )
            DPPMK ( IJ,K ) = AKAPPA
     &                       * ( GDP( IJ,K-1 ) - GDP( IJ,K   ) )
     &                       /  GDPM( IJ,K   ) / 2.  
 2200 CONTINUE
*
      DO 2300 K = 1, KMAX
         DO 2300 IJ = 1, IDIM*JDIM
*"
            IF ( GDT( IJ,K ) .LT. 0.0 ) THEN
               WRITE(6,*) 'P2CUMA(T<0),IJ,K,T',IJ,K,GDT(IJ,K)
            END IF
*"                       追加 95/05/23 (ishiwatari)

            GDQS  ( IJ,K ) = FQSAT( GDT( IJ,K ), GDP( IJ,K ) )
*"
            IF ( GDQS( IJ,K ) .EQ. 0.0 ) THEN
               WRITE(6,*) 'P2CUMA(q* <0),IJ,K,T',IJ,K,GDT(IJ,K)
            END IF
*"                       追加 94/12/17 (ishiwatari)
 2300 CONTINUE
*
*"         < 3. 調節 >
*
*
      CALL RESET ( ADJDB,  IDIM*JDIM )
      DO 3000 I = 1, IMAX
         DO 3000 J = 1, JMAX
            ADJDB( IDIM*(J-1)+I ) = 1.  
 3000 CONTINUE
*
*"         < 3.1 イテレーション >
*
      DO 3100 ITR = 1, ITRMAX
*
         CALL RESET ( ADJDA, IDIM*JDIM )
*
         DO 3110 K = 2, KMAX
*
            DO 3120 IJ = 1, IDIM*JDIM
*
               IF ( ADJDB( IJ ) .GT. 0.9   ) THEN
                  ST   =     GDT ( IJ,K-1 ) - GDT ( IJ,K )
     &                   + ( GDQS( IJ,K-1 ) - GDQS( IJ,K ) )*ELF
     &                   - DPPMK ( IJ,K )
     &                     * ( GDT( IJ,K-1 ) + GDT( IJ,K ) )
c                    write(6,*) 'P2CUMA:k,GDT(k-1),gdt(k)'
c                    write(6,*) k,gdt(ij,k-1),gdt(ij,k)
*
                  IF ( ST .GT. EPST( ITR ) ) THEN
*"                                   [ 不安定 ]
*
                     IF (     ( GDQ(IJ,K  ) .GE. RCRIT*GDQS(IJ,K) )
     &                   .AND.( GDQ(IJ,K-1).GE. RCRIT*GDQS(IJ,K-1) )
     &                  ) THEN
*"                                  変更 94/12/17 (ishiwatari)
*"                                   [ 飽和 ]
*
                        QEXE   = DELP( IJ,K-1 )
     &                           *( GDQ( IJ,K-1 )-GDQS( IJ,K-1 ) )
     &                         + DELP( IJ,K   )
     &                           *( GDQ( IJ,K   )-GDQS( IJ,K   ) )
*
                        DQSAT  = FDQSAT( GDT(IJ,K  ), GDQS(IJ,K  ) )
                        DQSATB = FDQSAT( GDT(IJ,K-1), GDQS(IJ,K-1) )
                        DHDT   = 1. + ELF * DQSAT
                        DHDTB  = 1. + ELF * DQSATB
*
                        STEXE = ST 
     &                        + ( 1. - DPPMK( IJ,K ) / DHDTB  )
     &                           * ELF * QEXE / DELP( IJ,K-1 )
*
*"                           < 温度の調節 >
*
                        DT  = STEXE
     &                      /(  ( 1. + DPDP( IJ,K ) )*DHDT
     &                        + DPPMK( IJ,K )
     &                          * ( 1. - DPDP( IJ,K )*DHDT/DHDTB ) )
*
                        DTB = - DHDT/DHDTB * DPDP( IJ,K ) * DT
     &                        + ELF * QEXE / ( DELP( IJ,K-1 ) * DHDTB )
*
                        GDT( IJ,K   ) = GDT( IJ,K   ) + DT
                        GDT( IJ,K-1 ) = GDT( IJ,K-1 ) + DTB
*
*"                           < 比湿の調節 >
*
                        GDQ ( IJ,K   ) = GDQS( IJ,K   ) + DQSAT  * DT
                        GDQ ( IJ,K-1 ) = GDQS( IJ,K-1 ) + DQSATB * DTB
*
                        GDQS( IJ,K   ) = GDQ( IJ,K   )
                        GDQS( IJ,K-1 ) = GDQ( IJ,K-1 )
*
                        ADJDA( IJ )     = 1.  
*
                     ENDIF
                  ENDIF
*
               ENDIF
 3120       CONTINUE
*
 3110    CONTINUE
*
         ADJPTS = 0.  
         DO 3130 IJ = 1, IDIM*JDIM
            ADJDB( IJ )     = ADJDA( IJ )
            ADJPTS = ADJPTS + ADJDA( IJ )
 3130    CONTINUE
*
         IF ( ADJPTS .LT. 1. ) GOTO 3200
*"                                  exit
*
 3100 CONTINUE
*
 3200 CONTINUE
*
*"         < 4. 降水量, 温度変化 >
*
      DO 4000 K = 1, KMAX
         DO 4000 IJ = 1, IDIM*JDIM
            GQCNV  ( IJ,K ) = GQCNV ( IJ,K )
     &                      + ( GDT ( IJ,K ) - GDTB( IJ,K ) ) / DELTA
*
            GSCNV  ( IJ,K ) = GSCNV ( IJ,K )
     &                      + ( GDQ ( IJ,K ) - GDQB( IJ,K ) ) / DELTA
*
            GQPR   ( IJ   ) = GQPR  ( IJ )
     &                      +( GDQB ( IJ,K ) - GDQ ( IJ,K ) )
     &                        * DELP( IJ,K ) / GRAV * EL / DELTA
 4000 CONTINUE
      RETURN
      END
