*"   PACKAGE  
* PACKAGE P2RDLN  !" 物理過程(2) 長波放射(多色) 
*
*" 履歴 93/08/29 保坂征宏            長波放射変更
*" 履歴 94/07/16 石渡正樹            上向き・下向き放射計算
***********************************************************************
      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       ODEPTH ( IDIM*JDIM, KMAX+1 ) !" 光学的深さ
      REAL       GDTM   ( 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       RUFLUX  ( IDIM*JDIM, KMAX+1 )     !" 上向きフラックス
      REAL       RDFLUX  ( IDIM*JDIM, KMAX+1 )     !" 下向きフラックス
      REAL       UFLUX  ( IDIM*JDIM, KMAX+1 )     !" 上向きフラックス
      REAL       DFLUX  ( IDIM*JDIM, KMAX+1 )     !" 下向きフラックス
      REAL       SUMWGL
      SAVE       RUFLUX,RDFLUX
*
*   [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
*
      CALL RESET ( RUFLUX, IDIM*JDIM*(KMAX+1) )
      CALL RESET ( RDFLUX, IDIM*JDIM*(KMAX+1) )
*
*"         < 1.5 光学的深さの計算 >
*
      CALL RESET ( ODEPTH, IDIM*JDIM*(KMAX+1) )
*
c$$$      write (6,*) 'in p2rdld:ij,k,odepth'
      DO 1500 IC = 1, NCLNG
         DO 1510 IJ = 1, IDIM*JDIM
            DO 1520 K = 1, KMAX+1
                   ODEPTH( IJ,K )
     &         =   PATHC * (   ABSQL(IC) * TAUQ( IJ,K )
     &                       + ABSAL(IC) * TAUA( IJ,K ) )
c$$$                   write (6,*) ij,k,odepth( ij,k )
 1520       CONTINUE
 1510    CONTINUE
 1500 CONTINUE
*
*"        < 温度の補間 : 正しくすべし!!!!>
      DO 1610 IJ = 1,IDIM*JDIM
         GDTM(IJ,1) =  GDTG( IJ )
         DO 1600 K = 2,KMAX
            GDTM(IJ,K) = ( GDT(IJ,K) + GDT(IJ,K-1))/2.0
 1600    CONTINUE
         GDTM(IJ,KMAX+1) = GDT(IJ,KMAX)
 1610 CONTINUE 
*
*"         < 2. 透過関数計算 >
*"                 光学的距離の計算追加・変更 94/08/23 石渡
*
c$$$      open(unit=71,file='transU',form='unformatted')
c$$$      open(unit=72,file='transD',form='unformatted')
      open(unit=71,file='transU')
      open(unit=72,file='transD')
*
      DO 2000 K = 1, KMAX+1
*
c         write (6,*) 'in p2rdln:k',k
*
         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( - ABS ( ODEPTH( IJ,KK )- ODEPTH( IJ,K ) ) )
*"                      weight のかけ方はここで良い???
*"                      良いわけがない. なんとかするべし.
 2120          CONTINUE
 2110       CONTINUE
 2100    CONTINUE
*
*"      < 3.   放射フラックス計算  >       !"  保坂 93/08/29 変更
*
         DO 3100 IJ = 1, IDIM*JDIM
            RUFLUX ( IJ,K ) = PIBGND ( IJ ) * TRANS ( IJ, 1 )
 3100    CONTINUE

c$$$         write (6,*) 'in p2rdld:before flux loop'
c         DO 3200 KK = 1,K-1
            DO 3200 IJ = 1, IDIM*JDIM
c$$$                RUFLUX( IJ,K ) 
c$$$     &             = RUFLUX ( IJ, K ) 
c$$$     &               - PIB ( IJ, KK ) 
c$$$     &                 * ( TRANS ( IJ, KK ) - TRANS ( IJ , KK+1 ) ) 
               IF( K .EQ. 1) THEN
                   RUFLUX( IJ,1 ) = RUFLUX ( IJ, 1 ) 
c$$$                write (6,*) 'upword:ij,k',ij,k
c$$$                write (6,*)
c$$$     &             - PIBGND ( IJ ) * EXP ( -  
c$$$     &                ( ODEPTH( IJ,1 ) - ODEPTH( IJ,1) )   ),
c$$$     &             PIBGND ( IJ ),
c$$$     &             stb*( GDTM(IJ,2)**4 - GDTM( IJ,1 )**4 )
c$$$     &               / ( ODEPTH(IJ,2) - ODEPTH( IJ,1 ) )
c$$$     &               *( - ( ODEPTH(IJ,1) - ODEPTH( IJ,1)  + 1.0 ) )
c$$$     &               * EXP ( -( ODEPTH(IJ,1) - ODEPTH(IJ,1) ) ),
c$$$     &             stb*( GDTM(IJ,2)**4 - GDTM( IJ,1 )**4 )
c$$$     &               / ( ODEPTH(IJ,2) - ODEPTH( IJ,1 ) )
                ELSE
                   RUFLUX( IJ,K ) 
     &             = RUFLUX ( IJ, K ) 
     &               + PIB ( IJ, K-1 )
     &                * ( - EXP ( - ( ODEPTH( IJ,1 ) - ODEPTH( IJ,K) )
     &                           ) + 1.0 )
     &               + STB * ( GDTM(IJ,K)**4 - GDTM( IJ,K-1 )**4 )
     &                     / ( ODEPTH(IJ,K) - ODEPTH( IJ,K-1) )
     &                 * ( - ( ODEPTH(IJ,1) - ODEPTH( IJ,K) + 1.0 )
     &                       * EXP ( - ( ODEPTH( IJ,1 ) - ODEPTH(IJ,K) )
     &                             )         + 1.0 )
c$$$                write (6,*) 'upword:ij,k',ij,k
c$$$                write (6,*)
c$$$     &             PIB ( IJ, K-1 )
c$$$     &                * EXP ( - ( ODEPTH( IJ,1 ) - ODEPTH( IJ,K) ) ),
c$$$     &             PIB ( IJ, K-1 ),
c$$$     &             stb*( GDTM(IJ,K)**4 - GDTM( IJ,K-1 )**4 )
c$$$     &               / ( ODEPTH(IJ,K) - ODEPTH( IJ,K-1 ) )
c$$$     &               *(  - ( ODEPTH(IJ,1) - ODEPTH( IJ,K)
c$$$     &                   + 1.0 ) )
c$$$     &               * EXP ( - ( ODEPTH(IJ,1) - ODEPTH(IJ,K) ) ),
c$$$     &             stb*( GDTM(IJ,K)**4 - GDTM( IJ,K-1 )**4 )
c$$$     &               / ( ODEPTH(IJ,K) - ODEPTH( IJ,K-1 ) )
             END IF
c$$$                if (ij .eq. 32 ) then
c$$$                   write (71,*) ij,k,kk,-(trans(ij,kk)-trans(ij,kk+1))
c$$$                endif
 3200   CONTINUE
*
c         DO 3300 KK = K, KMAX
            DO 3300 IJ = 1, IDIM*JDIM
c$$$                RDFLUX( IJ,K ) 
c$$$     &             = RDFLUX ( IJ, K ) 
c$$$     &               + PIB ( IJ, KK ) 
c$$$     &                 * ( TRANS ( IJ, KK ) - TRANS ( IJ , KK+1 ) ) 
               IF( K .EQ. 1) THEN
                   RDFLUX( IJ,K ) 
     &          = - PIBGND ( IJ )
     &               * ( EXP ( - ODEPTH( IJ,1 ) ) - 1.0 )
     &            - STB *   ( GDTM(IJ,2)**4 - GDTM( IJ,1 )**4 )
     &                    / ( ODEPTH(IJ,2)  - ODEPTH( IJ,1 )  )
     &               * ( - ( ODEPTH(IJ,1) + 1.0 )
     &                     * EXP ( - ODEPTH( IJ,1 ) )
     &                   + 1.0 )
c$$$                write (6,*) 'downword:ij,k',ij,k
c$$$                write (6,*)
c$$$     &           - pibgnd(ij) * EXP ( - ODEPTH( IJ,1 ) ),
c$$$     &             pibgnd(ij),
c$$$     &           - stb*( GDTM(IJ,2)**4 - GDTM( IJ,1 )**4 )
c$$$     &                /( ODEPTH(IJ,2) - ODEPTH( IJ,1 ) )
c$$$     &               *(- ( ODEPTH(IJ,1) + 1.0)
c$$$     &                   * EXP ( - ODEPTH( IJ,1 ) ) ),
c$$$     &           - stb*( GDTM(IJ,2)**4 - GDTM( IJ,1 )**4 )
c$$$     &                /( ODEPTH(IJ,2) - ODEPTH( IJ,1 ) )
               ELSE
                   RDFLUX( IJ,K ) 
     &          = - PIB ( IJ, K-1 )
     &               * ( EXP ( - ODEPTH( IJ,K ) ) - 1.0 )
     &            - STB *   ( GDTM(IJ,K)**4 - GDTM( IJ,K-1 )**4 )
     &                    / ( ODEPTH(IJ,K)  - ODEPTH( IJ,K-1 )  )
     &               * ( - ( ODEPTH(IJ,K) + 1.0 )
     &                     * EXP ( - PATHC * ODEPTH( IJ,K ) )
     &                   + 1.0 )
c$$$                write (6,*) 'downword:ij,k',ij,k
c$$$                write (6,*)
c$$$     &           - pib(ij,k-1)
c$$$     &               * EXP ( - ODEPTH( IJ,K ) ),
c$$$     &             pib(ij,k-1),
c$$$     &           - stb*( GDTM(IJ,K)**4 - GDTM( IJ,K-1 )**4 )
c$$$     &                /( ODEPTH(IJ,K) - ODEPTH( IJ,K-1 ) )
c$$$     &               *(- ( ODEPTH(IJ,K) + 1.0)
c$$$     &                   * EXP ( - ODEPTH( IJ,K ) ) ),
c$$$     &           - stb*( GDTM(IJ,K)**4 - GDTM( IJ,K-1 )**4 )
c$$$     &                /( ODEPTH(IJ,K) - ODEPTH( IJ,K-1 ) )
               END IF
c$$$                if (ij .eq. 32 ) then
c$$$                   write (72,*) ij,k,kk,trans(ij,kk)-trans(ij,kk+1)
c$$$                endif
 3300   CONTINUE
*
         DO 3400 KK = K, KMAX+1
            DO 3400 IJ = 1, IDIM*JDIM
                RFLUX( IJ,K ) 
     &             = RUFLUX( IJ, K ) - RDFLUX ( IJ, K ) 
c$$$                write (6,*) 'ij,k,ru,rd',ij,k,ruflux(ij,k),rdflux(ij,k)
 3400   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
*
      close(71)
      close(72)
*
      RETURN
***************************************************************
      ENTRY UDFLUX
     O      (UFLUX,DFLUX)
*
*"                            石渡 94/08/05 追加 
*
      CALL COPY (UFLUX,RUFLUX,IDIM*JDIM*(KMAX+1))
      CALL COPY (DFLUX,RDFLUX,IDIM*JDIM*(KMAX+1))
*
      END
