* PACKAGE P2RADM   !" 物理過程(2) 放射(多色)*
*
*" 履歴 93/08/29 保坂征宏            長波放射変更
*" 履歴 94/09/10 石渡正樹            温度半整数に変更
*********************************************************************
      SUBROUTINE RADFLX            !" 放射フラックス
     O         ( RFLUX , RMTXS , DRFLUX,
     O           SFLUX ,
     I           GDTM  , GDTG  , GDQ   , GDPM  ,
     I           IT    , IDATE , 
     C           ALON  , ALAT  , GRALB          )
*
*"       ｎ色大気
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                      !" 格子点数，波数
      INCLUDE   (ZCCOM)                      !" 標準物理定数
#else
#include        "zcdim.F"                    !" 格子点数，波数
#include        "zccom.F"                    !" 標準物理定数
#endif
*
*   [OUTPUT] 
      REAL       RFLUX  ( IDIM*JDIM, KMAX+1 )     !" 長波フラックス
      REAL       RMTXS  ( IDIM*JDIM,  -1:1 )      !" Ｔ陰解行列：地表
      REAL       DRFLUX ( IDIM*JDIM, KMAX+1, 0:1 )!" 長波地表温度変化
      REAL       SFLUX  ( IDIM*JDIM, KMAX+1 )     !" 日射フラックス
*
*   [INPUT] 
      REAL       GDTM   ( IDIM*JDIM, KMAX+1 ) !" 温度(半整数)  
*                                                !" 石渡 93/09/10 GDT と入替 
      REAL       GDTG   ( IDIM*JDIM )         !" 地表面温度
      REAL       GDQ    ( IDIM*JDIM, KMAX )   !" 比湿  ｑ
      REAL       GDPM   ( IDIM*JDIM, KMAX+1 ) !" 気圧(半整数)
      INTEGER    IT                           !" 時刻
      INTEGER    IDATE  ( 6 )                 !" 日付
*
      REAL       ALON   ( IDIM )              !" 経度
      REAL       ALAT   ( JDIM )              !" 緯度
      REAL       GRALB  ( IDIM*JDIM )         !" 地表アルベド
*
*   [INTERNAL PARM]
      CHARACTER  TUNIT *4                     !" INTLNG等の単位
      INTEGER    INTLNG                       !" 長波計算間隔
      INTEGER    INTSHT                       !" 短波計算間隔
      DATA       TUNIT  / 'HOUR'/
      DATA       INTLNG / 3     /
      DATA       INTSHT / 1     /
*
      NAMELIST  /NMRADM/ TUNIT , INTLNG, INTSHT
*
*   [INTERNAL SAVE] 
      REAL       GDTZ  ( IDIM*JDIM    )       !" 前のステップでの温度 
*                                                !" 保坂 93/08/29 追加
*
      INTEGER    IORADL                       !" 長波計算間隔：秒
      INTEGER    IORADS                       !" 短波計算間隔：秒
      INTEGER    ITLSTL                       !" 以前長波計算した時刻
      INTEGER    ITLSTS                       !" 以前短波計算した時刻
*
      LOGICAL    OFIRST
      DATA       OFIRST / .TRUE. /
      SAVE
*
*   [INTERNAL WORK] 
      COMMON    /COMWRK/
     &           TAUQ  , TAUA   , TRNSM1, TRNSM2, TRANS ,
     &           PIB   , PIBGND , SECZET
      REAL       TAUQ   ( IDIM*JDIM, KMAX+1 ) !" 光学的厚さ：水
      REAL       TAUA   ( IDIM*JDIM, KMAX+1 ) !" 光学的厚さ：空気
      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 )         !" 地表のπＢ
      REAL       SECZET ( IDIM*JDIM )         !" sec(入射角)
*
      INTEGER    IJ, K
      INTEGER    IFPAR, JFPAR
*
*   [ONCE]
      IF ( OFIRST ) THEN
         WRITE ( 6,* ) ' RADIATION PACKAGE - DATE=93/08/29 pochaka '
         OFIRST = .FALSE.
*
         CALL   REWNML ( IFPAR , JFPAR )
         READ   ( IFPAR, NMRADM, END=190 )
  190    WRITE  ( JFPAR, NMRADM )
*
         CALL ACTIME
     I      ( INTLNG, TUNIT  , '#'    ,
     O        IORADL                    )
         CALL ACTIME
     I      ( INTSHT, TUNIT  , '#'    ,
     O        IORADS                   )
         ITLSTL = IT - IORADL
         ITLSTS = IT - IORADS
      ENDIF
*
*"         < 1. τ の計算 >
*
      CALL RESET ( TAUQ ( 1,KMAX+1 ), IDIM*JDIM )
      CALL RESET ( TAUA ( 1,KMAX+1 ), IDIM*JDIM )
      DO 1110  K = KMAX, 1, -1
         DO 1120 IJ = 1, IDIM*JDIM
            TAUQ ( IJ,K ) = TAUQ ( IJ,K+1 )
     &                       + GDQ( IJ,K )
     &                         * ( GDPM( IJ,K )-GDPM( IJ,K+1 ) )/ GRAV
            TAUA ( IJ,K ) = TAUA ( IJ,K+1 )
     &                       +   ( GDPM( IJ,K )-GDPM( IJ,K+1 ) )/ GRAV
 1120       CONTINUE
 1110    CONTINUE
 1100 CONTINUE
*
*
      IF (  ( IT - ITLSTL ) .GE. IORADL ) THEN
*
         ITLSTL = IT
*
*"         < 2. 長波フラックスの算出 >     !"  保坂 93/08/29 変更
*
         CALL LNGFLX
     O         ( RFLUX , DRFLUX,
     I           GDTM  , GDTG  , TAUQ  , TAUA  ,
     W           TRNSM1, TRNSM2, TRANS , PIB   , PIBGND  )
*
         write (6,*) 'in p2radm;ij,k,rflux'
         do 1 ij=1,idim*jdim
            do 1 k=1,kmax+1
               write (6,*) ij,k,rflux(ij,k)
 1       continue 
*
      ELSE 
*
         DO 2000 K = 1, KMAX+1             
            DO 2000 IJ = 1, IDIM*JDIM
               RFLUX ( IJ,K ) = RFLUX ( IJ,K ) 
     &                        + DRFLUX ( IJ,K,1 ) 
     &                           * ( GDT ( IJ,1 ) - GDTZ ( IJ ) )
*                       "! しかるべく半整数に変更せよ
               DRFLUX ( IJ,K,1 ) = DRFLUX ( IJ,K,1 ) 
     &                             / GDTZ ( IJ   )**3
     &                             * GDT  ( IJ,1 )**3
 2000    CONTINUE

      ENDIF
*
      DO 2100 IJ = 1, IDIM*JDIM
         GDTZ ( IJ ) = GDT ( IJ , 1 )
 2100 CONTINUE
*
*"         < 3. 長波陰解用行列 >             !" 保坂 93/08/29 変更
*
      DO 3100 IJ = 1, IDIM*JDIM
         RMTXS ( IJ,0  ) =   DRFLUX( IJ,1,0 )
         RMTXS ( IJ,1  ) =   DRFLUX( IJ,1,1 )
         RMTXS ( IJ,-1 ) = 0.  
 3100 CONTINUE
*
*
      IF (  ( IT - ITLSTS ) .GE. IORADS ) THEN
*
         ITLSTS = IT
*
*"         < 4. 短波入射 >
*
         CALL SHTINS
     O         ( SFLUX(1,KMAX+1), SECZET ,
     I           IDATE ,
     C           ALON  , ALAT             )
*
*"         < 5. 短波フラックス >
*
         CALL SHTFLX
     M         ( SFLUX ,
     I           TAUQ  , TAUA  , SECZET,
     C           GRALB                  )
*
      ENDIF
*
      RETURN
      END
************************************************************************
      SUBROUTINE RADCOR                      !" 放射フラックス補正
     O         ( RFLUXC ,
     I           RFLUX  , DRFLUX, GTT   , GTTG  , DELT  )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                       !" 格子点数，波数
#else
#include        "zcdim.F"                     !" 格子点数，波数
#endif
*
*   [OUTPUT] 
      REAL       RFLUXC ( IDIM*JDIM, KMAX+1 )    !" 補正長波フラックス 
*
*   [INPUT] 
      REAL       RFLUX  ( IDIM*JDIM, KMAX+1 )    !" 長波フラックス
      REAL       DRFLUX ( IDIM*JDIM, KMAX+1, 0:1 )!" 長波地表温度変化
      REAL       GTT    ( IDIM*JDIM, KMAX )       !" 温度時間変化項
      REAL       GTTG   ( IDIM*JDIM )             !" 地表温度変化率
      REAL       DELT                             !" 時間刻みΔt
*
*   [INTERNAL WORK] 
      INTEGER    IJ, K
*
*"  <  放射フラックス補正 >                    !" 93/08/29 変更
*
      DO 1100 K = 1, KMAX+1
         DO 1100 IJ = 1, IDIM*JDIM
            RFLUXC ( IJ,K ) = RFLUX ( IJ,K )
     &                      + (  GTTG ( IJ )   * DRFLUX ( IJ,K,0 ) 
     &                         + GTT  ( IJ,1 ) * DRFLUX ( IJ,K,1 ) )
     &                        * DELT 
 1100 CONTINUE
*
      RETURN
      END
*********************************************************************
      SUBROUTINE RADTND               !" 放射加熱率
     O         ( GTTRL , GTTRS ,
     I           RFLUX , SFLUX , GDPM    )
*
*   [PARAM] 
#ifdef SYS_IBMS
      INCLUDE   (ZCDIM)                       !" 格子点数，波数
      INCLUDE   (ZCCOM)                       !" 標準物理定数
#else
#include        "zcdim.F"                     !" 格子点数，波数
#include        "zccom.F"                     !" 標準物理定数
#endif
*
*   [OUTPUT] 
      REAL       GTTRL  ( IDIM*JDIM, KMAX   ) !" 長波加熱率
      REAL       GTTRS  ( IDIM*JDIM, KMAX   ) !" 短波加熱率
*
*   [INPUT] 
      REAL       RFLUX  ( IDIM*JDIM, KMAX+1 ) !" 長波フラックス
      REAL       SFLUX  ( IDIM*JDIM, KMAX+1 ) !" 短波フラックス
      REAL       GDPM   ( IDIM*JDIM, KMAX+1 ) !" 気圧
*
*   [INTERNAL WORK] 
      INTEGER    IJ, K
*
*
      DO 1100 K = 1, KMAX
         DO 1110 IJ = 1, IDIM*JDIM
            GTTRL ( IJ,K )
     &         =  (  RFLUX( IJ,K ) - RFLUX( IJ,K+1 )   )
     &           /(  GDPM ( IJ,K ) - GDPM ( IJ,K+1 )   )
     &           / CP * GRAV
 1110    CONTINUE
 1100 CONTINUE
*
      DO 2100 K = 1, KMAX
         DO 2110 IJ = 1, IDIM*JDIM
            GTTRS ( IJ,K )
     &         =  (  SFLUX( IJ,K ) - SFLUX( IJ,K+1 )  )
     &           /(  GDPM ( IJ,K ) - GDPM ( IJ,K+1 )   )
     &           / CP * GRAV
 2110    CONTINUE
 2100 CONTINUE
*
      RETURN
      END
