subroutine SetMarsDustInit

    ! ファイル入出力補助
    ! File I/O support
    !
    use dc_iounit, only: FileOpen

    ! NAMELIST ファイル入力に関するユーティリティ
    ! Utilities for NAMELIST file input
    !
    use namelist_util, only: namelist_filename, NmlutilMsg, NmlutilAryValid

    ! メッセージ出力
    ! Message output
    !
    use dc_message, only: MessageNotify


    ! 宣言文 ; Declaration statements
    !

    integer:: unit_nml        ! NAMELIST ファイルオープン用装置番号.
                              ! Unit number for NAMELIST file open
    integer:: iostat_nml      ! NAMELIST 読み込み時の IOSTAT.
                              ! IOSTAT of NAMELIST read

    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /set_Mars_dust_nml/ DustSenario, DOD067, DustVerDistCoef
          !
          ! デフォルト値については初期化手続 "rad_Mars_V1#RadMarsV1Init"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "rad_Mars_V1#RadMarsV1Init" for the default values.
          !


    ! デフォルト値の設定
    ! Default values settings
    !
    DustSenario     = 'Const'
    DOD067          = 0.2_DP
    DustVerDistCoef = 0.01_DP

    ! NAMELIST の読み込み
    ! NAMELIST is input
    !
    if ( trim(namelist_filename) /= '' ) then
      call FileOpen( unit_nml, namelist_filename, mode = 'r' ) ! (in)

      rewind( unit_nml )
      read( unit_nml, nml = set_Mars_dust_nml, iostat = iostat_nml )             ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
    end if


    if ( DustSenario == 'Const' ) then
      IDDustSenario = IDDustSenarioConst
    else if ( DustSenario == 'VikingNoDS' ) then
      IDDustSenario = IDDustSenarioVikingNoDS
    else if ( DustSenario == 'Viking' ) then
      IDDustSenario = IDDustSenarioViking
    else if ( DustSenario == 'MGS' ) then
      IDDustSenario = IDDustSenarioMGS
    else
      call MessageNotify( 'E', module_name, 'DustSenario of %c is not supported.', c1 = trim( DustSenario ) )
    end if


    ! Initialization of modules used in this module
    !


    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'DustSenario     = %c', c1 = trim( DustSenario ) )
    call MessageNotify( 'M', module_name, 'DOD067          = %f', d  = (/ DOD067      /) )
    call MessageNotify( 'M', module_name, 'DustVerDistCoef = %f', d  = (/ DustVerDistCoef /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )


    set_Mars_dust_inited = .true.

  end subroutine SetMarsDustInit

  subroutine SetMarsDustSetDOD067( Ls, xyr_Press, xyz_Press, xyr_DOD067 )
    !
    ! 
    !
    ! Set dust optical depth at 0.67 micron
    !

    ! モジュール引用 ; USE statements
    !

    ! メッセージ出力
    ! Message output
    !
    use dc_message, only: MessageNotify

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: Grav


    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in ):: Ls
                              ! Ls
    real(DP), intent(in ):: xyr_Press    (0:imax-1, 1:jmax, 0:kmax)
                              ! Pressure
    real(DP), intent(in ):: xyz_Press    (0:imax-1, 1:jmax, 1:kmax)
                              ! Pressure
    real(DP), intent(out):: xyr_DOD067   (0:imax-1, 1:jmax, 0:kmax)
                              ! Optical depth

    ! 作業変数
    ! Work variables
    !
    real(DP)            :: DOD
    real(DP)            :: xy_DOD067    (0:imax-1, 1:jmax)
                              ! Dust optical depth at 0.67 micron
    real(DP)            :: xyz_MixRtDust(0:imax-1, 1:jmax, 1:kmax)
    real(DP)            :: xy_DODFac    (0:imax-1, 1:jmax)

    real(DP), parameter :: DustOptDepRefPress  = 610.0_DP
    real(DP), parameter :: DustVerDistRefPress = 610.0_DP

    real(DP)            :: MixRtDust0

    integer:: k               ! 鉛直方向に回る DO ループ用作業変数
                              ! Work variables for DO loop in vertical direction

    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. set_Mars_dust_inited ) call SetMarsDustInit


    select case ( IDDustSenario )
    case ( IDDustSenarioConst )
      xy_DOD067 = DOD067
    case ( IDDustSenarioVikingNoDS )

      call SetMarsDustDODVikingNoDS( Ls, DOD )
      xy_DOD067 = DOD

    case ( IDDustSenarioViking )

      call SetMarsDustDODViking( Ls, DOD )
      xy_DOD067 = DOD

    case ( IDDustSenarioMGS )
      call MessageNotify( 'E', module_name, 'DustSenario of %c is not supported.', c1 = trim( DustSenario ) )
    case default
      call MessageNotify( 'E', module_name, 'DustSenario of %c is not supported.', c1 = trim( DustSenario ) )
    end select


    MixRtDust0      =   1.0_DP

    xyz_MixRtDust = MixRtDust0 * exp( DustVerDistCoef * ( 1.0_DP - ( DustVerDistRefPress / xyz_Press ) ) )

    k = kmax
    xyr_DOD067(:,:,k) = 0.0_DP
    do k = kmax-1, 0, -1
      xyr_DOD067(:,:,k) = xyr_DOD067(:,:,k+1) + xyz_MixRtDust(:,:,k+1) * ( xyr_Press(:,:,k) - xyr_Press(:,:,k+1) ) / Grav
    end do

    xy_DODFac = xy_DOD067 * xyr_Press(:,:,0) / DustOptDepRefPress / xyr_DOD067(:,:,0)
    do k = 0, kmax
      xyr_DOD067(:,:,k) = xyr_DOD067(:,:,k) * xy_DODFac
    end do


  end subroutine SetMarsDustSetDOD067

  subroutine SetMarsDustDODViking( Ls, DOD )

    real(DP), intent(in ) :: Ls
    real(DP), intent(out) :: DOD


    !
    ! Local variables
    !
    real(DP) :: DODDS1
    real(DP) :: DODDS2
    real(DP) :: DSLs
    real(DP) :: MaxDOD
    real(DP) :: DSDTC


    call SetMarsDustDODVikingNoDS( Ls, DOD )

    ! Add two dust storms
    !
    DSLs   = 210.0_DP
    MaxDOD = 2.7_DP
    DSDTC  = 50.0_DP
    call SetMarsDustDSExp( Ls, DSLs, MaxDOD, DSDTC, DOD )

    DSLs   = 280.0_DP
    MaxDOD = 4.0_DP
    DSDTC  = 50.0_DP
    call SetMarsDustDSExp( Ls, DSLs, MaxDOD, DSDTC, DOD )

    DOD = max( DOD, DODDS1, DODDS2 )


  end subroutine SetMarsDustDODViking

  subroutine SetMarsDustDODVikingNoDS( Ls, DOD )

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants0, only: PI

    real(DP), intent(in ) :: Ls
    real(DP), intent(out) :: DOD


    ! This expression is obtained from Lewis et al. [1999].
    !
    DOD = 0.7_DP + 0.3_DP * cos( ( Ls + 80.0_DP ) * PI / 180.0_DP )


  end subroutine SetMarsDustDODVikingNoDS

  subroutine SetMarsDustDSExp( Ls, DSLs, MaxDOD, DSDTC, DOD )

    real(DP), intent(in ) :: Ls
    real(DP), intent(in ) :: DSLs
    real(DP), intent(in ) :: MaxDOD
    real(DP), intent(in ) :: DSDTC
    real(DP), intent(out) :: DOD


    ! Local variables
    !
    real(DP) :: TMPLs

    if( Ls < DSLs ) then
      TMPLs = Ls + 360.0_DP
    else
      TMPLs = Ls
    endif

    DOD = MaxDod * exp( -( TMPLs - DSLs ) / DSDTC )


  end subroutine SetMarsDustDSExp