!=================================================================
!  一次元火星大気ダストの放射計算プログラム;
!  ダストの可視光放射計算
!
!  2000/10/15  小高正嗣
!  2000/11/09  小高正嗣; 2D モデルへ
!=================================================================
      
      subroutine CLDIRAD
     I   ( TPOT  , PSFC  , ICRAD ,
     I     QTAU  , QDTAU , 
     I     DINU1 , DINU2 , DELNU ,
     O     QRAD  , 
     B     PRES0 , DENS0 , TEMP0 , TPOT0 , 
     C     QALB  , QASM  , 
     C     RDRY  , CPDRY , TSFC  , GRAV  , PI   , 
     W     TEMP  , BBRAD , BBRADS, 
     W     FXRADU, FXRADD, FXRADN, 
     W     FXRADUW,FXRADDW, 
     W     ASR   , BSR   , CSR   , DSR  ,
     W     RSR   , SSR   ,
     G     FDZ   , FDZM  , DX    ,
     D     NX    , NXB   , NZ    , NZB  , NROT   )

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

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

      dimension 
     \   QTAU (-2:NXB,-2:NZB), QDTAU(-2:NXB,-2:NZB) 

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

      dimension
     \   TEMP(-2:NXB,-2:NZB)

      dimension
     \   FXRADU(-2:NXB,-2:NZB), FXRADD(-2:NXB,-2:NZB), 
     \   FXRADN(-2:NXB,-2:NZB), QRAD  (-2:NXB,-2:NZB)

      dimension
     \   TSFC( -2:NXB )

      dimension 
     \   ASR(0:NZ,0:NZ), BSR(0:NZ,0:NZ), CSR(0:NZ,0:NZ), 
     \   DSR(0:NZ,0:NZ), RSR(0:NZ), SSR(0:NZ)

      dimension
     \   FDZ(-2:NZB), FDZM(-2:NZB)

      dimension
     \   FXRADUW(-2:NXB,-2:NZB), FXRADDW(-2:NXB,-2:NZB), 
     \   BBRAD  (-2:NXB,-2:NZB)

      dimension
     \   BBRADS (-2:NXB)

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

      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 ix = -2, NXB
         call PLANK 
     I      ( DINU1, DINU2, TSFC(IX), DELNU,
     C        PI , 
     O        BBRADS(IX) )
      end do

      do iz = -2, NZB
      do ix = -2, NXB
         call PLANK 
     I      ( DINU1, DINU2, TEMP(IX,IZ), DELNU,
     C        PI , 
     O        BBRAD(ix,iz) )
      end do
      end do

!-----------------------------------------------------------------      
! Eddington 近似のパラメータ
      
      call EDDIN
     I   ( QASM  , QALB  , 1.0D0,  
     O     gamma1, gamma2, gamma3 )

!-----------------------------------------------------------------      
      DO IX = -2, NXB
!-----------------------------------------------------------------      
! tau = 0 の場合

      if ( QTAU(IX,0) .eq. 0.0 ) then

         do iz = 0, NZ
            FXRADU(IX,IZ) = BBRADS(IX)
            FXRADD(IX,IZ) = 0.0D0
         end do

!-----------------------------------------------------------------      
      else
!-----------------------------------------------------------------      
! 散乱行列の計算

      call CLSCMTRX
     I   ( QTAU  , QDTAU  , gamma1, gamma2, gamma3,  
     I     BBRAD , 1.0    , IX    , 
     C     ALBEDO, QALB   ,
     O     ASR   , BSR    , CSR   , DSR   ,
     O     RSR   , SSR    ,
     D     NX    , NXB    , NZ  , NZB              )

! 境界条件の変更

      BSR(0,0) = 0.0D0
      RSR(0) = BBRADS(IX)

      do iz = 1, NZ
         RSR(iz) = 2.0D0 * QDTAU(IX,iz-1) * 
     \        ( 1.0D0 - QALB ) * BBRAD(IX,iz-1)
      end do

      do iz = 0, NZ-1
         SSR(iz) = RSR(iz+1)
      end do
      SSR(NZ) = 0.0D0

!-----------------------------------------------------------------      
! 散乱光の計算

      call CLSCRAD
     I   ( ASR, BSR, CSR, DSR, 
     I     RSR, SSR, 
     I     IX ,
     O     FXRADU , FXRADD ,
     W     FXRADUW, FXRADDW,
     D     NX     , NXB    , NZ   , NZB      )

      end if
!-----------------------------------------------------------------      

      FXRADU(IX,NZ+1) = FXRADU(IX,NZ)
      FXRADD(IX,NZ)   = 0.0D0 

!-----------------------------------------------------------------      
      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     QRAD  ,
     C     CPDRY , DX   , GRAV ,
     D     NX    , NXB  , NZ   , NZB          )

!-----------------------------------------------------------------      
      return
      end

