!=================================================================
!  一次元火星大気ダストの放射計算プログラム;
!  散乱行列の計算
!
!  2000/10/15  小高正嗣
!  2000/11/09  小高正嗣; 2D モデルへ
!=================================================================
      
      subroutine CLSCMTRX
     I   ( QTAU  , QDTAU  , gamma1, gamma2, gamma3,  
     I     SXSOL , COSZ   , IXRAD ,
     C     ALBEDO, QALB_SR,
     O     ASR   , BSR    , CSR   , DSR   ,
     O     RSR   , SSR    ,
     D     NX    , NXB    , NZ  , NZB              )

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

      dimension 
     \   QTAU (-2:NXB,-2:NZB), QDTAU(-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)

!-----------------------------------------------------------------      
      do jz = 0, NZ
         do iz = 0, NZ
            ASR(iz,iz) = 0.0D0
            BSR(iz,iz) = 0.0D0
            CSR(iz,iz) = 0.0D0
            DSR(iz,iz) = 0.0D0
         end do
      end do
!-----------------------------------------------------------------      

      do iz = 1, NZ
         ASR(iz,iz) = - 0.5D0 * gamma1 * QDTAU(IXRAD,iz-1)
      end do

      do iz = 1, NZ
         ASR(iz,iz-1) = 1.0D0 + ASR(iz,iz)
      end do

!-----------------------------------------------------------------      

      BSR(0,0) = ALBEDO 
      do iz = 1, NZ
         BSR(iz,iz) = 0.5D0 * gamma2 * QDTAU(IXRAD,iz-1) 
      end do

      do iz = 1, NZ
         BSR(iz,iz-1) = BSR(iz,iz)
      end do

!-----------------------------------------------------------------      

      do iz = 0, NZ-1
         CSR(iz,iz) = 0.5D0 * gamma2 * QDTAU(IXRAD,iz)
      end do

      do iz = 0, NZ-2
         CSR(iz,iz+1) = CSR(iz,iz)
      end do

!-----------------------------------------------------------------      

      do iz = 0, NZ-1
         DSR(iz,iz) = - 0.5D0 * gamma1 * QDTAU(IXRAD,iz) 
      end do

      do iz = 0, NZ-2
         DSR(iz,iz+1) = 1.0D0 + DSR(iz,iz)
      end do

!-----------------------------------------------------------------      

      RSR(0) = ALBEDO * SXSOL(IXRAD,0)
      do iz = 1, NZ
         RSR(iz) = gamma3 * QALB_SR 
     \        * QDTAU(IXRAD,iz-1) * SXSOL(IXRAD,NZ+1) 
     \        * exp ( - 0.5D0 * 
     \               ( QTAU(IXRAD,iz) + QTAU(IXRAD,iz-1) ) / COSZ )
      end do

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

!-----------------------------------------------------------------      

      return
      end
