!=================================================================
!  CO2 大気放射過程計算プログラム
!
!  2000/05/12 小高正嗣; 射出率スキーム用
!  2000/05/12 小高正嗣; バンドモデル
!  2000/05/23 小高正嗣; 高速化
!=================================================================

      subroutine CLRAD
     I   ( TPOT  , TAUN  , PSFC ,
     I     ICRAD ,
     O     DTPRAD,
     B     TPOT0 , TEMP0 , PRES0 , DENS0 , FDZ  ,
     C     WNU   ,
     C     TSFC  , CPDRY , RDRY  , DX    , DELNU,  GRAV,
     W     TEMP  , BBRAD , BBRAD_S, 
     W     FXRADU, FXRADD, FXRADN, FXRADU_NU, FXRADD_NU,
     D     NX    , NXB   , NZ    , NZB   , NROT  , NBND  )

!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      

      dimension
     \   TPOT  ( -2:NXB, -2:NZB, NROT ),
     \   DTPRAD( -2:NXB, -2:NZB )

      dimension
     \   TAUN( 0:NZ+1, 0:NZ+1, NBND ),
     \   WNU(NBND)

      dimension
     \   TSFC  ( -2:NXB )

      dimension
     \   TPOT0 ( -2:NZB ), TEMP0( -2:NZB ), 
     \   PRES0 ( -2:NZB ), DENS0( -2:NZB ), FDZ  ( -2:NZB ) 

      dimension
     \   TEMP  ( -2:NXB, -2:NZB ),
     \   BBRAD ( -2:NXB, -2:NZB ),
     \   BBRAD_S( -2:NXB ),
     \   FXRADU( -2:NXB, -2:NZB ),
     \   FXRADD( -2:NXB, -2:NZB ),
     \   FXRADN( -2:NXB, -2:NZB )

      DIMENSION
     \   FXRADU_NU(-2:NXB,-2:NZB,NBND), 
     \   FXRADD_NU(-2:NXB,-2:NZB,NBND)

      parameter ( SIGMA = 5.67D-8 )

!-----------------------------------------------------------------      
! 初期化

      do iz = -2, NZB
      do ix = -2, NXB
         FXRADU(ix,iz) = 0.0D0
         FXRADD(ix,iz) = 0.0D0
         FXRADN(ix,iz) = 0.0D0
         BBRAD (ix,iz) = 0.0D0
      end do
      end do

      do inb = 1, NBND
      do iz = -2, NZB
      do ix = -2, NXB 
         FXRADU_NU(ix,iz,inb) = 0.0D0
         FXRADD_NU(ix,iz,inb) = 0.0D0
      end do
      end do
      end do

      do ix = -2, NXB
         BBRAD_S(ix) = 0.0D0
      end do

!-----------------------------------------------------------------      
! 温度場の計算

      do iz = 0, NZ
      do ix = -2, NXB
         TEMP( ix, iz ) = ( TPOT( ix, iz, ICRAD ) + TPOT0( iz ) )  
     \        * ( PRES0( iz ) / PSFC ) ** ( RDRY / CPDRY )
      end do
      end do

      do ix = -2, NXB
         TEMP( ix, -1 )   = TEMP( ix, 0 )
         TEMP( ix, -2 )   = TEMP( ix, 0 )
         TEMP( ix, NZ+1 ) = TEMP( ix, NZ )
         TEMP( ix, NZ+2 ) = TEMP( ix, NZ )
      end do

!-----------------------------------------------------------------      
! バンドごとの寄与を計算

      do inb = 1, NBND

!-----------------------------------------------------------------      
! 放射場の計算; 黒体放射

      call CLBBRAD 
     I   ( TEMP , WNU(inb),
     O     BBRAD, BBRAD_S,
     C     TSFC , 
     D     NX   , NXB,  
     D     NZ   , NZB  )

!-----------------------------------------------------------------      
! 放射場の計算; 上向き放射束密度

      call CLFXRDU
     I   ( BBRAD , BBRAD_S, TAUN, inb , 
     O     FXRADU,
     D     NZ    , NZB    , NX  , NXB , NBND )

!-----------------------------------------------------------------      
! 放射場の計算; 下向き放射束密度

      call CLFXRDD
     I   ( BBRAD , TAUN , inb ,  
     O     FXRADD,
     D     NZ    , NZB  , NX  , NXB , NBND )

         do iz =  0, NZB
         do ix = -2, NXB
            FXRADU_NU(ix,iz,inb) = FXRADU(ix,iz)
            FXRADD_NU(ix,iz,inb) = FXRADD(ix,iz)
         end do
         end do

         do iz = -2, NZB
         do ix = -2, NXB
            FXRADU(ix,iz) = 0.0D0
            FXRADD(ix,iz) = 0.0D0
         end do
         end do

      end do

!-----------------------------------------------------------------      
! 波長積分

      do iz = -2, NZB
      do ix = -2, NXB
      do inb = 1, NBND/2
         FXRADU(ix,iz) = FXRADU(ix,iz) + DELNU * 
     \        ( FXRADU_NU(ix,iz,inb) + FXRADU_NU(ix,iz,NBND+1-inb) )
         FXRADD(ix,iz) = FXRADD(ix,iz) + DELNU * 
     \        ( FXRADD_NU(ix,iz,inb) + FXRADD_NU(ix,iz,NBND+1-inb) )
      end do
      end do
      end do


! 地表面からの寄与

      do iz = -2, NZB
      do ix = -2, NXB
         FXRADU(ix,iz) = sigma * TSFC(ix)**4 + FXRADU(ix,iz)
      end do
      end do

!-----------------------------------------------------------------      
! 放射場の計算; 正味の放射束密度

      call CLFXRDN
     I   ( FXRADU, FXRADD,
     O     FXRADN,
     D     NX    , NXB   , 
     D     NZ    , NZB     )

!-----------------------------------------------------------------      
! 加熱率の計算

      call CLQRAD
     I   ( FXRADN, DENS0, TEMP0, TPOT0, PRES0, 
     I     FDZ   , 
     O     DTPRAD,
     C     CPDRY , DX   , GRAV ,
     D     NX    , NXB  , NZ   , NZB          )

      return

      end
