* PACKAGE P2RADM   
*
*" 履歴 93/08/29 保坂征宏            長波放射変更
*********************************************************************
      SUBROUTINE RADFLX            
     O         ( RFLUX , RMTXS , DRFLUX,
     O           SFLUX ,
     I           GDT   , GDTG  , GDQ   , GDPM  ,
     I           IT    , IDATE , 
     C           ALON  , ALAT  , GRALB          )
*
*"       ｎ色大気
*
*   [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

*
*   [OUTPUT] 
      ReaL*8       RFLUX  ( IDIM*JDIM, KMAX+1 )     
      ReaL*8       RMTXS  ( IDIM*JDIM,  -1:1 )      
      ReaL*8       DRFLUX ( IDIM*JDIM, KMAX+1, 0:1 )
      ReaL*8       SFLUX  ( IDIM*JDIM, KMAX+1 )     
*
*   [INPUT] 
      ReaL*8       GDT    ( IDIM*JDIM, KMAX   ) 
*                                                
      ReaL*8       GDTG   ( IDIM*JDIM )         
      ReaL*8       GDQ    ( IDIM*JDIM, KMAX )   
      ReaL*8       GDPM   ( IDIM*JDIM, KMAX+1 ) 
      INTEGER    IT                           
      INTEGER    IDATE  ( 6 )                 
*
      ReaL*8       ALON   ( IDIM )              
      ReaL*8       ALAT   ( JDIM )              
      ReaL*8       GRALB  ( IDIM*JDIM )         
*
*   [INTERNAL PARM]
      CHARACTER  TUNIT *4                     
      INTEGER    INTLNG                       
      INTEGER    INTSHT                       
      DATA       TUNIT  / 'HOUR'/
      DATA       INTLNG / 3     /
      DATA       INTSHT / 1     /
*
      NAMELIST  /NMRADM/ TUNIT , INTLNG, INTSHT
*
*   [INTERNAL SAVE] 
      ReaL*8       GDTZ  ( IDIM*JDIM    )       
*                                                
*
      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*8       TAUQ   ( IDIM*JDIM, KMAX+1 ) 
      ReaL*8       TAUA   ( IDIM*JDIM, KMAX+1 ) 
      ReaL*8       TRNSM1 ( IDIM*JDIM, KMAX+1 ) 
      ReaL*8       TRNSM2 ( IDIM*JDIM, KMAX+1 ) 
*                                                 
      ReaL*8       TRANS  ( IDIM*JDIM, KMAX+1 ) 
      ReaL*8       PIB    ( IDIM*JDIM, KMAX+1 ) 
      ReaL*8       PIBGND ( IDIM*JDIM )         
      ReaL*8       SECZET ( IDIM*JDIM )         
*
      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. 長波フラックスの算出 >     
*
         CALL LNGFLX
     O         ( RFLUX , DRFLUX,
     I           GDT   , GDTG  , TAUQ  , TAUA  ,
     W           TRNSM1, TRNSM2, TRANS , PIB   , PIBGND  )
*
      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. 長波陰解用行列 >             
*
      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] 




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

*
*   [OUTPUT] 
      ReaL*8       RFLUXC ( IDIM*JDIM, KMAX+1 )    
*
*   [INPUT] 
      ReaL*8       RFLUX  ( IDIM*JDIM, KMAX+1 )    
      ReaL*8       DRFLUX ( IDIM*JDIM, KMAX+1, 0:1 )
      ReaL*8       GTT    ( IDIM*JDIM, KMAX )       
      ReaL*8       GTTG   ( IDIM*JDIM )             
      ReaL*8       DELT                             
*
*   [INTERNAL WORK] 
      INTEGER    IJ, K
*
*"  <  放射フラックス補正 >                    
*
      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] 





* 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

*
*   [OUTPUT] 
      ReaL*8       GTTRL  ( IDIM*JDIM, KMAX   ) 
      ReaL*8       GTTRS  ( IDIM*JDIM, KMAX   ) 
*
*   [INPUT] 
      ReaL*8       RFLUX  ( IDIM*JDIM, KMAX+1 ) 
      ReaL*8       SFLUX  ( IDIM*JDIM, KMAX+1 ) 
      ReaL*8       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
