!=================================================================
!  一次元火星大気ダストの放射計算プログラム;
!  ダストの可視光放射計算
!
!  2000/10/15  小高正嗣
!  2000/11/09  小高正嗣; 2D モデルへ
!=================================================================
      
      subroutine CLDSRAD
     I   ( QTAU  , QDTAU  , S0    , COSZ  , 
     O     QRAD  , SXSOL  , 
     B     TPOT0 , TEMP0  , PRES0 , DENS0 ,
     C     ALBEDO, QALB_SR, QASM_SR,
     C     DELNU , CPDRY  , TSOL   , SIGMA, 
     C     GRAV  , PI     , 
     W     SXRADU, SXRADD , SXRADN , 
     W     SXRADUW, SXRADDW,
     W     ASR   , BSR    , CSR    , DSR  ,
     W     RSR   , SSR    ,
     G     FDZ   , FDZM   , DX     ,
     D     NX    , NXB    , NZ     , NZB           )

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

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

      dimension
     \   SXRADU(-2:NXB,-2:NZB), SXRADD(-2:NXB,-2:NZB), 
     \   SXRADN(-2:NXB,-2:NZB), QRAD  (-2:NXB,-2:NZB),
     \   SXSOL (-2:NXB,-2:NZB)
      
      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
     \   SXRADUW(-2:NXB,-2:NZB), SXRADDW(-2:NXB,-2:NZB)

!-----------------------------------------------------------------      
! 可視波長域における波長積分されたプランク関数の計算
      
      DSNU1 = 2000.0D0
      DSNUM = 1.0D+4
      DSNU2 = 1.0D+5 

      DELNU1 = DELNU * 40.0 
      DELNU2 = DELNU * 400.0

      call PLANK
     I   ( DSNU1, DSNUM, TSOL, DELNU1,
     C     PI , 
     O     SXSOL1 )

      call PLANK
     I   ( DSNUM, DSNU2, TSOL, DELNU2,
     C     PI , 
     O     SXSOL2 )

      do IX = -2, NXB

         SXSOL(IX,NZ+1) = 0.0D0 ! 初期化
         SXSOL(IX,NZ+1) = ( SXSOL1 + SXSOL2 )*
     \        ( S0 / (SIGMA * TSOL ** 4.0 ) ) 

         do iz = 0, NZ
            SXSOL(IX,iz) = SXSOL(IX,NZ+1) * 
     \           exp ( - QTAU(IX,iz) / COSZ )
         end do

      end do
!-----------------------------------------------------------------      
! Eddington 近似のパラメータ
      
      call EDDIN
     I   ( QASM_SR, QALB_SR, COSZ, 
     O     gamma1, gamma2, gamma3 )

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

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

         do iz = 0, NZ
            SXRADU(IX,IZ) = ALBEDO * SXSOL(IX,0) * COSZ
            SXRADD(IX,IZ) = 0.0D0
         end do

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

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

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

      call CLSCRAD
     I   ( ASR, BSR, CSR, DSR, 
     I     RSR, SSR, 
     I     IX ,
     O     SXRADU, SXRADD,
     W     SXRADUW, SXRADDW,
     D     NX     , NXB    , NZ   , NZB      )

      end if

      SXRADU(IX,NZ+1) = SXRADU(IX,NZ)
      SXRADD(IX,NZ+1) = 0.0D0

!-----------------------------------------------------------------      
! 直達光の寄与を加算

      do iz = 0, NZ+1
         SXRADD(IX,iz) = SXRADD(IX,iz) + SXSOL(IX,iz) * COSZ
      end do

!-----------------------------------------------------------------      
      END DO
!-----------------------------------------------------------------      
! 正味のフラックスを計算

      call CLFXRDN
     I   ( SXRADU, SXRADD,
     O     SXRADN,
     D     NX    , NXB   , 
     D     NZ    , NZB     )

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

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

