*"   PACKAGE  TRPPAS       !"  kenkaimen
*"   
*"      Tue May 10 17:52:16 1994        Ishiwatari Masaki
***********************************************************************
      SUBROUTINE TRPPAS
     M        ( GDT, GDQ,
     O          TRPLEV,
     I          GDTG  ,
     I          GDPM, 
     I          RFLUX )
*   [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       GDTG   ( IDIM*JDIM )   !" 
*   [OUTPUT] 
      INTEGER    TRPLEV ( IDIM*JDIM )            !" 圏界面の高さ
*
*   [INPUT] 
      REAL       RFLUX  ( IDIM*JDIM, KMAX+1 )   !" 長波フラックス
      REAL       GDPM   ( IDIM*JDIM, KMAX+1 )   !" 
*
*   [INTERNAL WORK] 
c$$$      COMMON    /COMWRK/
c$$$     &           TAUQ  , TAUA   , TRNSM1, TRNSM2, TRANS ,
c$$$     &           PIB   , PIBGND , SECZET
      REAL       TAUQ   ( IDIM*JDIM, KMAX+1 ) !" 光学的厚さ：水
      REAL       TAUA   ( IDIM*JDIM, KMAX+1 ) !" 光学的厚さ：空気
      REAL       UQ   ( IDIM*JDIM, KMAX+1 ) !" 光路長：水
      REAL       UA   ( IDIM*JDIM, KMAX+1 ) !" 光路長：空気
      REAL       PIB    ( IDIM*JDIM, KMAX+1 ) !" πＢ＝σＴ**4
      REAL       PIBGND ( IDIM*JDIM )         !" 地表のπＢ
      INTEGER       LFNMAX ( IDIM*JDIM ) !" 
*
*   [WORK] 
      REAL       TRNSM1 ( IDIM*JDIM, KMAX+1 ) !" 1/2 レベルからの透過係数
      REAL       TRNSM2 ( IDIM*JDIM, KMAX+1 ) !" 3/2 レベルからの透過係数
*                                                !" 保坂 93/08/29 追加
      REAL       TRANS  ( IDIM*JDIM, KMAX+1 ) !" 透過係数
      REAL       WPIB   ( IDIM*JDIM, KMAX+1 ) !" πＢ＝σＴ**4
      REAL       WPIBG  ( IDIM*JDIM )         !" 地表のπＢ
      REAL       DRFLUX ( IDIM*JDIM, KMAX+1, 0:1 )!" 長波地表温度変化
*
*   [INTERNAL WORK] 
      REAL       ZZ
      REAL       ZZB
      INTEGER    IJ, K
*
*          < 1. パラメータ設定 >
*
      DO 1000 IJ = 1, IDIM*JDIM
         LFNMAX ( IJ ) =   1
 1000 CONTINUE
*
*          < 2. 圏界面の当たりをつける >
*
      write (6,*) 'in p2trps:table of rflux'
      do 111 ij=1,idim*jdim
         do 111 k=1,kmax+1
            write (6,*) ij,k,rflux(ij,k)
 111  continue 
*
      DO 2100 IJ=1, IDIM*JDIM
         DO 2000 K=KMAX,1,-1
            IF( RFLUX( IJ,K ) .LT. RFLUX( IJ,K+1 ) ) THEN
               LFNMAX(IJ) = K+1
               GOTO 8888
            END IF
 2000    CONTINUE 
 8888    CONTINUE 
 2100 CONTINUE 
      DO 2200 IJ=1, IDIM*JDIM
         TRPLEV ( IJ ) = LFNMAX ( IJ )
         write (6,*) 'in trppas:the level of max rflux'
         write (6,*) ij, trplev(ij)
 2200 CONTINUE 
*
*"         < 3. 微調整 >
*
      DO 3100 IJ = 1, IDIM*JDIM
         CALL STRTAU
     O        ( TAUQ, TAUA,
     I          GDQ, GDPM ,
     I          TRPLEV )
*
         PIB ( IJ,LFNMAX(IJ) ) = RFLUX ( IJ,LFNMAX(IJ) ) /2.0
     &                  * (1.5*TAUQ( IJ,LFNMAX(IJ) ) + 1.0 )
         ZZB = PIB ( IJ,LFNMAX(IJ) ) - STB * GDT( IJ,LFNMAX( IJ ) )**4
         write (6,*) 'in trppas:ij,trplev,rflux,tauq,pib,stb*T^4,zz'
         write (6,*) ij,TRPLEV(IJ), rflux ( IJ,TRPLEV(IJ) ),
     &               tauq( IJ,TRPLEV(IJ) )
         write (6,*) '  ', pib ( ij,TRPLEV(IJ) ), 
     &               stb * gdt( IJ,TRPLEV( IJ ) )**4,zzb
*
         IF( ZZB .LE. 0.0 ) THEN
            TRPLEV( IJ ) = LFNMAX( IJ )
         ELSE IF( ZZB .GT. 0.0 ) THEN
            DO 3200 K= LFNMAX(IJ)-1,1,-1
               TRPLEV( IJ ) = K
               CALL STRTAU
     O        ( TAUQ, TAUA,
     I          GDQ, GDPM ,
     I          TRPLEV )
*
               PIB ( IJ,K ) = RFLUX ( IJ,K ) /2.0
     &                     * (1.5*TAUQ( IJ,K ) + 1.0 )
               ZZ = PIB ( IJ,K ) - STB * GDT( IJ,K )**4
         write (6,*) 'in trppas:ij,trplev,rflux,tauq,pib,stb*T^4,zz'
         write (6,*) ij,TRPLEV(IJ), rflux ( IJ,TRPLEV(IJ) ),
     &               tauq( IJ,TRPLEV(IJ) )
         write (6,*) '  ', pib ( ij,TRPLEV(IJ) ), 
     &               stb * gdt( IJ,TRPLEV( IJ ) )**4,zz
               IF( ZZ .GT. 0.0 .AND. ZZ .LE. ZZB ) THEN
                    CONTINUE
               ELSE IF( ZZ .GT. 0.0 .AND. ZZ .GT. ZZB ) THEN
                    TRPLEV( IJ ) = K+1
                    GOTO 9999
               ELSE IF( ZZ .LE. 0.0 .AND. ABS(ZZ) .LE. ABS(ZZB) ) THEN
                    TRPLEV( IJ ) = K  
                    GOTO 9999
               ELSE IF( ZZ .LE. 0.0 .AND. ABS(ZZ) .GT. ABS(ZZB) ) THEN
                    TRPLEV( IJ ) = K+1
                    GOTO 9999
               END IF
               ZZB = ZZ
 3200       CONTINUE
 9999    CONTINUE       !" TRPLEV( IJ ) 決定した
         END IF
*"             一回目で決まらなかった場合終了.
*
 3100 CONTINUE 
*
      CALL STRTAU
     O        ( TAUQ, TAUA,
     I          GDQ, GDPM ,
     I          TRPLEV )
*
*"        < 4. 成層圏の上向きフラックス決定 >
*
      DO 4100 IJ = 1, IDIM*JDIM
         write (6,*) 'in trppas:the tropopause level'
         write (6,*) ij, trplev(ij)
         DO 4000 K = TRPLEV(IJ), KMAX+1
            RFLUX( IJ,K ) = 2.0 * STB 
     &                  * ( GDT( IJ,TRPLEV( IJ ) ) )**4
     &                  /(1.5 * TAUQ( IJ,TRPLEV( IJ ) ) + 1.0)
 4000    CONTINUE 
 4100 CONTINUE 
*
*"        < 5. 成層圏の温度構造 >
*
      DO 5100 IJ = 1, IDIM*JDIM
         DO 5000 K=TRPLEV(IJ),KMAX
            PIB ( IJ,K ) = RFLUX ( IJ,K ) /2.0
     &                       * (1.5*TAUQ( IJ,K ) + 1.0 )
            GDT ( IJ,K ) = ( PIB( IJ,K )/STB )**0.25
            GDQ ( IJ,K ) = GDQ ( IJ,TRPLEV ( IJ ) )
 5000    CONTINUE 
 5100 CONTINUE 
*

      RETURN
      END

