!=================================================================
!  一次元火星大気ダストの放射計算プログラム;
!  ダストの光学的厚さの計算
!
!  2000/10/15  小高正嗣
!  2000/11/09  小高正嗣; 2D モデルへ
!=================================================================

      subroutine SETQTAU
     I   ( DENS0, QVAP , 
     I     ITQ  , 
     O     QTAU , QDTAU,
     O     QTAUIR1, QDTAUIR1, QTAUIR2, QDTAUIR2,
     C     QREFF  , QDENS   , QALB   , QASM    , 
     C     QEXT_SR, QEXT_IR1, QEXT_IR2, PI     ,
     G     FDZ  , DX,
     D     NX   , NXB  , NZ   , NZB   , NROT    )

!-----------------------------------------------------------------      
      implicit real*8 ( a-h,o-z )
!-----------------------------------------------------------------      
      dimension
     \   DENS0(-2:NZB), QVAP(-2:NXB,-2:NZB,NROT)

      dimension 
     \   QTAU   (-2:NXB,-2:NZB), QDTAU   (-2:NXB,-2:NZB),
     \   QTAUIR1(-2:NXB,-2:NZB), QDTAUIR1(-2:NXB,-2:NZB),
     \   QTAUIR2(-2:NXB,-2:NZB), QDTAUIR2(-2:NXB,-2:NZB)

      dimension
     \   FDZ(-2:NZB)

!-----------------------------------------------------------------      
! 可視光に対する光学的厚さ

      BETA = ( 3.0D0 * QEXT_SR ) / ( 4.0D0 * PI * QREFF * QDENS ) 

      do iz = -2, NZB
         do ix = -2, NXB
            QDTAU(ix,iz) = BETA * ( 1.0D0 - QALB * QASM**2.0D0 ) * 
     \           DENS0(iz) * QVAP(ix,iz,ITQ) * FDZ(iz) * DX
         end do
      end do

      do ix = -2, NXB
         QTAU(ix,NZ) = 0.0D0 
         do iz = NZ-1, 0, -1
            QTAU(ix,iz) = QTAU(ix,iz+1) + QDTAU(ix,iz)
         end do
      end do

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

!-----------------------------------------------------------------      
! 赤外に対する光学的厚さ

      do iz = -2, NZB
         do ix = -2, NXB
            QTAUIR1 (ix,iz) = QTAU (ix,iz) * 0.5
            QDTAUIR1(ix,iz) = QDTAU(ix,iz) * 0.5
         end do
      end do

      do iz = -2, NZB
         do ix = -2, NXB
            QTAUIR2 (ix,iz) = QTAUIR1 (ix,iz) * QEXT_IR2 / QEXT_IR1
            QDTAUIR2(ix,iz) = QDTAUIR1(ix,iz) * QEXT_IR2 / QEXT_IR1
         end do
      end do


      return
      end
