* PACKAGE P2RADL  !" 物理過程(2) 長波放射(多色)
*
*" 履歴 93/08/29 保坂征宏            長波放射変更
***********************************************************************
      SUBROUTINE LNGFLX           !" 長波放射フラックス
     O         ( RFLUX , DRFLUX,
     I           GDT   , GDTG  , TAUQ  , TAUA  ,
     W           TRNSM1, TRNSM2, TRANS , PIB   , PIBGND  )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                       !" 格子点数，波数
      INCLUDE   (ZCCOM)                       !" 標準物理定数
#else
#include        "zcdim.F"                     !" 格子点数，波数
#include        "zccom.F"                     !" 標準物理定数
#endif
*
*   [OUTPUT] 
      REAL       RFLUX  ( IDIM*JDIM, KMAX+1 )     !" 長波フラックス
      REAL       DRFLUX ( IDIM*JDIM, KMAX+1, 0:1 )!" 長波地表温度変化
*
*   [INPUT] 
      REAL       GDT    ( IDIM*JDIM, KMAX   ) !" 温度(整数)
      REAL       GDTG   ( IDIM*JDIM )         !" 地表面温度
      REAL       TAUQ   ( IDIM*JDIM, KMAX+1 ) !" 光学的厚さ：水
      REAL       TAUA   ( IDIM*JDIM, KMAX+1 ) !" 光学的厚さ：大気
*
*   [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       PIB    ( IDIM*JDIM, KMAX+1 ) !" πＢ＝σＴ**4
      REAL       PIBGND ( IDIM*JDIM )         !" 地表のπＢ
*
*   [INTERNAL PARM] 
      INTEGER    NCLNGD
      PARAMETER (NCLNGD = 10)
      INTEGER    NCLNG                        !" 長波バンド数
      REAL       ABSQL  ( NCLNGD )            !" 水の吸収係数
      REAL       ABSAL  ( NCLNGD )            !" 空気の吸収係数
      REAL       RWGTL  ( NCLNGD )            !" バンドウェイト
      REAL       PATHC                        !" 光路長のファクタ
      DATA       NCLNG  / 4 /
      DATA       ABSQL  / 8.0  , 1.0  , 0.1  , 0.0  , 6*0.   /
      DATA       ABSAL  / 0.0  , 0.0  , 0.0  , 5.E-5, 6*0.   /
      DATA       RWGTL  / 0.2  , 0.1  , 0.1  , 0.6  , 6*0.   /
      DATA       PATHC  / 1.5   /
*
      NAMELIST  /NMRADL/ NCLNG , ABSQL , ABSAL , RWGTL , PATHC
      SAVE
*
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
      SAVE       OFIRST
*
*   [INTERNAL WORK] 
      INTEGER    IJ, K, KK, IC
      INTEGER    IFPAR, JFPAR
      REAL       SUMWGL
*
*   [ONCE] 
      IF ( OFIRST ) THEN
         WRITE ( 6,* ) ' LONGWAVE RADIATION N-BAND DATE=93/08/29'
         OFIRST = .FALSE.
*
         CALL   REWNML ( IFPAR , JFPAR )
         READ   ( IFPAR, NMRADL, END=190 )
  190    WRITE  ( JFPAR, NMRADL )
*
         SUMWGL = 0.  
         DO 210 IC = 1, NCLNG
            SUMWGL = SUMWGL + RWGTL( IC )
  210    CONTINUE
         DO 220 IC = 1, NCLNG
            RWGTL( IC ) =  RWGTL( IC ) / SUMWGL
  220    CONTINUE
*
      ENDIF
*
*"         < 1. πＢの計算 >
*
      DO 1100 K = 1, KMAX
         DO 1110 IJ = 1, IDIM*JDIM
            PIB  ( IJ,K ) = STB * ( GDT( IJ,K )**4 )
 1110    CONTINUE
 1100 CONTINUE
*
      DO 1200 IJ = 1, IDIM*JDIM
         PIBGND( IJ ) = STB * ( GDTG( IJ )**4 )
 1200 CONTINUE
*
*"         < 2. 透過関数計算 >
*
      DO 2000 K = 1, KMAX+1
*
         CALL RESET ( TRANS, IDIM*JDIM*(KMAX+1) )
*
         DO 2100 IC = 1, NCLNG
            DO 2110 KK = 1, KMAX+1
               DO 2120 IJ = 1, IDIM*JDIM
                  TRANS( IJ,KK )
     &            = TRANS( IJ,KK )
     &              +  RWGTL ( IC )
     &                  * EXP( - PATHC
     &                         * (  ABSQL(IC)*ABS(  TAUQ( IJ,KK )
     &                                            - TAUQ( IJ,K  )  )
     &                            + ABSAL(IC)*ABS(  TAUA( IJ,KK )
     &                                            - TAUA( IJ,K  )  )  )
     &                       )
 2120          CONTINUE
 2110       CONTINUE
 2100    CONTINUE
*
*"      < 3.   放射フラックス計算  >       !"  保坂 93/08/29 変更
*
         DO 3100 IJ = 1, IDIM*JDIM
            RFLUX ( IJ,K ) = PIBGND ( IJ ) * TRANS ( IJ, 1 )
 3100    CONTINUE

         DO 3200 KK = 1, KMAX
            DO 3200 IJ = 1, IDIM*JDIM
                RFLUX( IJ,K ) 
     &             = RFLUX ( IJ, K ) 
     &               - PIB ( IJ, KK ) 
     &                 * ( TRANS ( IJ, KK ) - TRANS ( IJ , KK+1 ) ) 
 3200   CONTINUE
*
*"       < 4.  補正項計算用透過関数 >       !"  保坂 93/08/29 変更
*
        DO 4000 IJ= 1 ,IDIM*JDIM
            TRNSM1 ( IJ,K ) = TRANS( IJ,1 )
            TRNSM2 ( IJ,K ) = TRANS( IJ,2 )
 4000   CONTINUE
*
 2000 CONTINUE
*
*"         < 5. 地表面補正用 >               !"  保坂 93/08/29 変更
*
      DO 5100 K = 1, KMAX+1
         DO 5110 IJ = 1, IDIM*JDIM
            DRFLUX( IJ,K,0 ) = 4. * PIBGND( IJ   ) / GDTG  ( IJ )
     &                            * TRNSM1( IJ,K )
            DRFLUX( IJ,K,1 ) = 4. * PIB   ( IJ,1 ) / GDT  ( IJ,1 )
     &                         * ( TRNSM2 ( IJ,K ) - TRNSM1 ( IJ,K ) )      
 5110    CONTINUE
 5100 CONTINUE
*
      RETURN
      END
