** PACKAGE P2ADBT  
*
*"  [HIS]  94/05/02(ishiwatari)
*"  [HIS]  95/12/03(ishiwatari)
**********************************************************************
      SUBROUTINE ADIABT       
     M         ( GDT   , GDQ   ,
     I           GDP   , GDPM  , GDZ   ,
     I           GDTG  , GDQG  )
*
*" 湿潤断熱構造を引く.
*" 一般場合で 1 次近似(線形項だけを取り出したもの)を伸ばしていく
*
*   [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

*
*   [MODIFY] 
      ReaL*8       GDT   ( IDIM*JDIM, KMAX )   
      ReaL*8       GDQ   ( IDIM*JDIM, KMAX )   
      ReaL*8       GQPR  ( IDIM*JDIM       )   
      ReaL*8       GQCNV ( IDIM*JDIM, KMAX )   
      ReaL*8       GSCNV ( IDIM*JDIM, KMAX )   
*
*   [INPUT] 
      ReaL*8       GDP    ( IDIM*JDIM, KMAX )  
      ReaL*8       GDPM   ( IDIM*JDIM, KMAX+1 )
      ReaL*8       GDZ    ( IDIM*JDIM, KMAX )  
      ReaL*8       GDTG   ( IDIM*JDIM )  
      ReaL*8       GDQG   ( IDIM*JDIM )  
C      ReaL*8       DELTA                       
*
*   [INTERNAL PARM] 
      INTEGER    ITRMAX
      PARAMETER (ITRMAX=10)
      ReaL*8       RH                       
      DATA       RH / 1.0 /
*
      NAMELIST  /NMMADJ/ RH
*
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
      SAVE
*
*   [INTERNAL WORK] 
      COMMON    /COMWRK/
     &           ADJDA , ADJDB ,
     &           GDTB  , GDQB  ,
     &           DELP  , DPDP  , DPPMK ,
     &           GDQS  
      ReaL*8       ADJDA ( IDIM*JDIM )         
      ReaL*8       ADJDB ( IDIM*JDIM )         
      ReaL*8       GDTB  ( IDIM*JDIM, KMAX )   
      ReaL*8       GDQB  ( IDIM*JDIM, KMAX )   
      ReaL*8       DELP  ( IDIM*JDIM, KMAX )   
      ReaL*8       DPDP  ( IDIM*JDIM, KMAX )   
      ReaL*8       DPPMK ( IDIM*JDIM, KMAX )   
      ReaL*8       GDQS  ( IDIM*JDIM, KMAX )   
      ReaL*8       GDPM34( IDIM*JDIM       )   
*
      INTEGER    IJ, K, ITR, I, J
      INTEGER    IFPAR, JFPAR
      ReaL*8       ELF, AKAPPA
      ReaL*8       ST, QEXE, QSAT, DQSAT, DQSATB, DT, STEXE
      ReaL*8       DHDT, DHDTB, DTB, ADJPTS
*
*" adibat 計算用
      ReaL*8       QSAT1, QSAT2
      ReaL*8       MM
      ReaL*8       D1, D2
      ReaL*8       TP, TM
      ReaL*8       M1, M2
      ReaL*8       F
      ReaL*8       QP, QM
      ReaL*8       A,C
      ReaL*8       PP1, PP2
*      
*   [INTERNAL FUNC]




* 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 )

*
*   [ONCE] 
      IF ( OFIRST ) THEN
         WRITE ( 6,* ) ' MOIST ADIABAT DATE=96/06/07 ishiwatari'
         OFIRST = .FALSE.
*
         CALL   REWNML ( IFPAR , JFPAR )
         READ   ( IFPAR, NMMADJ, END=190 )
  190    WRITE  ( JFPAR, NMMADJ )
*
      ENDIF
*
*"         < 2. ファクター計算 >
*
      ELF    = EL/CP
      AKAPPA = RAIR / CP
      DO 2000 IJ=1,IDIM*JDIM
         GDPM34( IJ ) = (   (   GDPM( IJ,1 ) ** (1. + AKAPPA )
     &                        - GDP ( IJ,1 ) ** (1. + AKAPPA ) )
     &                    / ( (   GDPM( IJ,1 ) - GDP( IJ,1 ) )
     &                      * (1. + AKAPPA) )
     &                  )
     &                   **( 1./AKAPPA )
 2000 CONTINUE 
*"             これインチキ. 本当は補間公式にすべし!?
*
c$$$      write(6,*) 'in p2adbt:k,gdp,gdpm'
c$$$      do 110 ij=1,idim*jdim
c$$$         do 100 k=1,5
c$$$            write(6,*) ij,k,gdp(ij,k),gdpm(ij,k)
c$$$ 100     continue 
c$$$            write(6,*) gdpm34( ij )
c$$$ 110  continue 
*
      DO 2100 K = 1, KMAX
         DO 2100 IJ = 1, IDIM*JDIM
            DELP  ( IJ,K ) = GDPM( IJ,K   ) - GDPM( IJ,K+1 )
 2100 CONTINUE
*
      DO 2300 IJ = 1, IDIM*JDIM
         DO 2200 K = 2, KMAX
            DPDP  ( IJ,K ) = DELP( IJ,K ) / DELP( IJ,K-1 )
*"              DPDP は後使ってない 
            DPPMK ( IJ,K ) = AKAPPA
     &                       * ( GDP( IJ,K-1 ) - GDP( IJ,K   ) )
     &                       /  GDPM( IJ,K   ) / 2.  
 2200    CONTINUE
         DPPMK ( IJ,1 ) = AKAPPA
     &                       * ( GDPM( IJ,1 ) - GDP( IJ,1   ) )
     &                       /  GDPM34( IJ ) / 2.  
 2300 CONTINUE
*
*"         < 3. T, q 構造決定>
*
*
*"         < 3.1 イテレーション(Newton 法) >
*
      DO 3100 IJ = 1, IDIM*JDIM
         DO 3110 K = 1, KMAX
c         DO 3110 K = 2, KMAX
            DO 3120 ITR = 1, ITRMAX
*
               IF(ITR .EQ. 1) THEN
                  IF( K .EQ. 1) THEN
                     GDT( IJ,1 ) = GDTG( IJ )
                  ELSE
                     GDT(IJ,K) = GDT(IJ,K-1)
c                     write(6,*) '*220*', ij,k,gdt(ij,k),gdt(ij,k-1)
                  END IF
               END IF
*
               IF (K .EQ. 1) THEN
                  QSAT1 = FQSAT ( GDTG( IJ ), GDPM( IJ,1 ) )
                  QSAT2 = FQSAT ( GDT( IJ,K   ), GDP( IJ,K   ) )
                  D1 = FDQSAT( GDTG( IJ   ), QSAT1  )
                  D2 = FDQSAT( GDT ( IJ,K ), QSAT2  )
                  TP = GDTG( IJ ) + GDT( IJ,K )
                  TM = GDTG( IJ ) - GDT( IJ,K )
                  write(6,*) 'in adbt', GDPM(IJ,1),GDP(IJ,1)
                  PP1 = GDPM(IJ,1)/((GDPM(IJ,1)+GDP(IJ,1))/2.0)
                  PP2 = GDP (IJ,1)/((GDPM(IJ,1)+GDP(IJ,1))/2.0)
               ELSE
                  QSAT1 = FQSAT ( GDT( IJ,K-1 ), GDP( IJ,K-1 ) )
                  QSAT2 = FQSAT ( GDT( IJ,K   ), GDP( IJ,K   ) )
                  D1 = FDQSAT( GDT( IJ,K-1   ), QSAT1  )
                  D2 = FDQSAT( GDT( IJ,K     ), QSAT2  )
                  TP = GDT( IJ,K-1 ) + GDT( IJ,K )
                  TM = GDT( IJ,K-1 ) - GDT( IJ,K )
                  PP1 = GDP(IJ,K-1)/GDPM(IJ,K)
                  PP2 = GDP(IJ,K  )/GDPM(IJ,K)
               END IF
               M1 = 1.0 - QSAT1
               M2 = 1.0 - QSAT2
               MM = 1.0 - 0.5 * QSAT1 - 0.5 * QSAT2
               QP = QSAT1 + QSAT2
               QM = QSAT1 - QSAT2
*
               ST =   AKAPPA/4.0 * TP**2 * MM
     &               * (PP1 *M1 - PP2*M2)
     &             - 0.5 * MM**2 * TP * TM
     &             + ELF/2.0 * QP * MM * TM
     &             - ELF/2.0 * TP * QM
               C = - AKAPPA/4.0
     &               * (   TP**2 * MM * PP2 * D2
     &                   + (PP1 * M1 - PP2 * M2)
     &                     * ( - 0.5 * TP**2 * D2
     &                         + 2.0 * TP * MM)     )
     &             + 0.5 * (   MM**2 *(-TP + TM)
     &                       - MM * TP * TM * D2 )
     &             - ELF/2.0 * ( - QP * MM - 0.5 * QP * TM * D2
     &                   + MM * TM * D2 )
     &             + ELF/2.0 * ( - TP * D2 + QM )
               DT = ST/C
*
               GDT  ( IJ,K ) = GDT( IJ,K ) + DT
               GDQ  ( IJ,K ) = RH * FQSAT ( GDT( IJ,K ), GDP( IJ,K ) )
*
c$$$               write(6,*) 'at the end of p2adbt:ij,k,dt,gdt,gdq'
c$$$               write(6,*) ij,k,dt,gdt( ij,k ), gdq( ij,k )
*
 3120       CONTINUE
 3110    CONTINUE
 3100 CONTINUE
*
      RETURN
      END
