!= 
!
!= Dust distribution is set
!
! Authors::   Yoshiyuki O. Takahashi
! Version::   $Id: set_Mars_dust.f90,v 1.13 2013/09/21 14:40:52 yot Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!
module set_Mars_dust
  !
  != 
  !
  != Dust distribution is set
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 
  !
  ! 
  !
  !
  !== References
  !
  !  Lewis, S. R., Collins, M., Forget, F., 
  !    Mars climate database v3.0 detailed design document, 
  !    Technical Note. Contract 11369/95/NL/JG. Work Package 7, ESA, 2001.
  !
  !== Procedures List
  !
!!$  ! RadiationFluxDennouAGCM :: 放射フラックスの計算
!!$  ! RadiationDTempDt        :: 放射フラックスによる温度変化の計算
!!$  ! RadiationFluxOutput     :: 放射フラックスの出力
!!$  ! RadiationFinalize       :: 終了処理 (モジュール内部の変数の割り付け解除)
!!$  ! ------------            :: ------------
!!$  ! RadiationFluxDennouAGCM :: Calculate radiation flux
!!$  ! RadiationDTempDt        :: Calculate temperature tendency with radiation flux
!!$  ! RadiationFluxOutput     :: Output radiation fluxes
!!$  ! RadiationFinalize       :: Termination (deallocate variables in this module)
  !
  !== NAMELIST
  !
!!$  ! NAMELIST#radiation_DennouAGCM_nml
  !

  ! USE statements
  !

  !
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! Double precision.
    &                 STRING, &  ! Strings.
    &                 TOKEN      ! Keywords.

  ! 格子点設定
  ! Grid points settings
  !
  use gridset, only: imax, & ! 経度格子点数.
                             ! Number of grid points in longitude
    &                jmax, & ! 緯度格子点数.
                             ! Number of grid points in latitude
    &                kmax    ! 鉛直層数.
                             ! Number of vertical level

  implicit none

  private


  ! 公開変数
  ! Public variables
  !
  logical, save, public:: set_Mars_dust_inited = .false.
                              ! 初期設定フラグ.
                              ! Initialization flag

  ! Private variables
  !
  real(DP), save :: DustExtEff
  real(DP), save :: REff
  real(DP), save :: RhoDust


  character(STRING), save :: DustScenario
  real(DP)         , save :: DOD067       ! Dust optical depth at 0.67 micron.
  real(DP)         , save :: DustVerDistCoef

  real(DP)         , save :: DustOptDepRefPress
  real(DP)         , save :: DustVerDistRefPress

  integer          , save      :: IDDustScenario
  integer          , parameter :: IDDustScenarioConst          = 1
  integer          , parameter :: IDDustScenarioVikingNoDS     = 2
  integer          , parameter :: IDDustScenarioViking         = 3
  integer          , parameter :: IDDustScenarioMGS            = 4
  integer          , parameter :: IDDustScenarioMGSDODFromFile = 5

  character(STRING), save      :: DODFileName
  character(STRING), save      :: DODVarName

  public :: SetMarsDustCalcDOD067
  public :: SetMarsDustSetDOD067
  public :: SetMarsDustInit

  character(*), parameter:: module_name = 'set_Mars_dust'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: set_Mars_dust.f90,v 1.13 2013/09/21 14:40:52 yot Exp $'
                              ! モジュールのバージョン
                              ! Module version

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

contains

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

  subroutine SetMarsDustCalcDOD067( &
    & xyr_Press, xyz_QDust,         & ! (in)
    & xyr_DOD067                    & ! (out)
    & )
    !
    ! 
    !
    ! Calculate dust optical depth at 0.67 micron
    !

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

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

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

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


    ! 宣言文 ; Declaration statements
    !

    real(DP), intent(in ):: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
                              ! Pressure
    real(DP), intent(in ):: xyz_QDust (0:imax-1, 1:jmax, 1:kmax)
                              ! Dust mixing ratio
    real(DP), intent(out):: xyr_DOD067(0:imax-1, 1:jmax, 0:kmax)
                              ! Optical depth

    ! 作業変数
    ! Work variables
    !
    real(DP)            :: xyz_DelDOD(0:imax-1, 1:jmax, 1:kmax)

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

    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. set_Mars_dust_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    do k = 1, kmax
      xyz_DelDOD(:,:,k) =                                            &
        &   3.0_DP / 4.0_DP * DustExtEff / ( REff * RhoDust * Grav ) &
        & * xyz_QDust(:,:,k)                                         &
        & * ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) )
    end do

    k = kmax
    xyr_DOD067(:,:,k) = 0.0_DP
    do k = kmax-1, 0, -1
      xyr_DOD067(:,:,k) = xyr_DOD067(:,:,k+1) + xyz_DelDOD(:,:,k+1)
    end do


    ! ヒストリデータ出力
    ! History data output
    !


  end subroutine SetMarsDustCalcDOD067

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

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

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

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

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoPut

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN

    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, only: &
      & PI                    ! $ \pi $.
                              ! 円周率. Circular constant

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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & y_Lat    ! $ \varphi $ [rad.] . 緯度. Latitude

    ! 時系列データの読み込み
    ! Reading time series
    !
    use read_time_series, only: SetValuesFromTimeSeriesWrapper


    ! 宣言文 ; 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)            :: xy_MaxHeightDust(0:imax-1, 1:jmax)

    real(DP)            :: MixRtDust0

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

    ! 実行文 ; Executable statement
    !

    ! 初期化
    ! Initialization
    !
    if ( .not. set_Mars_dust_inited ) then
      call MessageNotify( 'E', module_name, 'This module has not been initialized.' )
    end if


    select case ( IDDustScenario )
    case ( IDDustScenarioConst )
      xy_DOD067 = DOD067

      ! Height of dust top
      xy_MaxHeightDust = 70.0d3

    case ( IDDustScenarioVikingNoDS )

      call SetMarsDustDODVikingNoDS( &
        & Ls, & ! (in)
        & DOD & ! (out)
        & )
      xy_DOD067 = DOD

      ! Height of dust top
!!$      xy_MaxHeightDust = 70.0d3
      !
      do j = 1, jmax
        xy_MaxHeightDust(:,j) =                                 &
          &   60.0d3                                            &
          & + 18.0d3 * sin( ( Ls - 158.0_DP ) * PI / 180.0_DP ) &
          & - 22.0d3 * sin( y_Lat(j) )**2
      end do

    case ( IDDustScenarioViking )

      call SetMarsDustDODViking( &
        & Ls, & ! (in)
        & DOD & ! (out)
        & )
      xy_DOD067 = DOD

      ! Height of dust top
!!$      xy_MaxHeightDust = 70.0d3
      !
      do j = 1, jmax
        xy_MaxHeightDust(:,j) =                                 &
          &   60.0d3                                            &
          & + 18.0d3 * sin( ( Ls - 158.0_DP ) * PI / 180.0_DP ) &
          & - 22.0d3 * sin( y_Lat(j) )**2
      end do

    case ( IDDustScenarioMGS )

      call SetMarsDustDODMGS(          &
        & Ls,                          &
        & xy_DOD067, xy_MaxHeightDust  &
        & )

    case ( IDDustScenarioMGSDODFromFile )

      call SetMarsDustDODMGS(          &
        & Ls,                          &
        & xy_DOD067, xy_MaxHeightDust  &
        & )

      call SetValuesFromTimeSeriesWrapper( &
        & 'DOD',                   &
        & DODFileName, DODVarName, &
        & xy_DOD067                &               ! (inout)
        & )

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


    MixRtDust0      =   1.0_DP

    do k = 1, kmax
      xyz_MixRtDust(:,:,k) = MixRtDust0 &
        & * exp( DustVerDistCoef        &
        &        * ( 1.0_DP - ( DustVerDistRefPress / xyz_Press(:,:,k) )**(70.0d3/xy_MaxHeightDust) ) &
        &      )
    end do
    xyz_MixRtDust = min( xyz_MixRtDust, MixRtDust0 )


    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


    ! ヒストリデータ出力
    ! History data output
    !
    call HistoryAutoPut( TimeN, 'DustPresc'    , xyz_MixRtDust    )
    call HistoryAutoPut( TimeN, 'DustMaxHeight', xy_MaxHeightDust )


  end subroutine SetMarsDustSetDOD067

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

  subroutine SetMarsDustDODViking( &
    & Ls, & ! (in)
    & DOD & ! (out)
    & )

    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, & ! (in)
      & DOD & ! (out)
      & )

    ! Add two dust storms
    !
    DSLs   = 210.0_DP
    MaxDOD = 2.7_DP
    DSDTC  = 50.0_DP
    call SetMarsDustDSExp( &
      & Ls, DSLs, MaxDOD, DSDTC, & ! (in)
      & DODDS1                   & ! (out)
      & )

    DSLs   = 280.0_DP
    MaxDOD = 4.0_DP
    DSDTC  = 50.0_DP
    call SetMarsDustDSExp( &
      & Ls, DSLs, MaxDOD, DSDTC, & ! (in)
      & DODDS2                   & ! (out)
      & )

    DOD = max( DOD, DODDS1, DODDS2 )


  end subroutine SetMarsDustDODViking

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

  subroutine SetMarsDustDODVikingNoDS( &
    & Ls, & ! (in)
    & DOD & ! (out)
    & )

    ! 物理定数設定
    ! 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 SetMarsDustDODMGS( &
    & Ls,                       &
    & xy_DOD, xy_MaxHeight      &
    & )

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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & y_Lat    ! $ \varphi $ [rad.] . 緯度. Latitude

    real(DP), intent(in ) :: Ls
    real(DP), intent(out) :: xy_DOD      (0:imax-1, 1:jmax)
    real(DP), intent(out) :: xy_MaxHeight(0:imax-1, 1:jmax)


    ! Local variables

    real(DP) :: DODEq
    real(DP) :: DODSouth
    real(DP) :: DODNorth
    real(DP) :: LsFactor

    integer  :: j


    DODEq    =                &
      &   0.2_DP              &
      & + ( 0.5_DP - 0.2_DP ) &
      &   * cos( ( Ls - 250.0_DP ) / 2.0_DP * PI / 180.0_DP )**14
    DODSouth =                &
      &   0.1_DP              &
      & + ( 0.5_DP - 0.1_DP ) &
      &   * cos( ( Ls - 250.0_DP ) / 2.0_DP * PI / 180.0_DP )**14
    DODNorth = 0.1_DP


    LsFactor = sin( ( Ls - 160.0_DP ) * PI / 180.0_DP )

    do j = 1, jmax

      if( y_Lat(j) > 0.0_DP ) then
        ! wrong
!!$        xy_DOD(:,j) = DODNorth              &
!!$          & + 0.5_DP * ( DODEq - DODNorth ) &
!!$          &   * ( 1.0_DP + tanh( ( 45.0_DP * PI / 180.0_DP - y_Lat(j) ) / 10.0_DP ) )
        xy_DOD(:,j) = DODNorth              &
          & + 0.5_DP * ( DODEq - DODNorth ) &
          &   * ( 1.0_DP + tanh( ( 45.0_DP * PI / 180.0_DP - y_Lat(j) ) * 10.0_DP ) )
      else
        ! wrong
!!$        xy_DOD(:,j) = DODSouth              &
!!$          & + 0.5_DP * ( DODEq - DODSouth ) &
!!$          &   * ( 1.0_DP + tanh( ( 45.0_DP * PI / 180.0_DP + y_Lat(j) ) / 10.0_DP ) )
        xy_DOD(:,j) = DODSouth              &
          & + 0.5_DP * ( DODEq - DODSouth ) &
          &   * ( 1.0_DP + tanh( ( 45.0_DP * PI / 180.0_DP + y_Lat(j) ) * 10.0_DP ) )
      end if

      xy_MaxHeight(:,j) =                                         &
        &   60.0_DP                                               &
        & + 18.0_DP * LsFactor                                    &
        & - ( 32.0_DP + 18.0_DP * LsFactor ) * sin( y_Lat(j) )**4 &
        & - 8.0_DP * LsFactor * sin( y_Lat(j) )**5
      xy_MaxHeight(:,j) = xy_MaxHeight(:,j) * 1.0d3

    end do


  end subroutine SetMarsDustDODMGS

  !----------------------------------------------------------------------------
!!$
!!$  subroutine SetMarsDustDODMGS( &
!!$    & Ls, &
!!$    & Lat, &
!!$    & DOD  &
!!$    & )
!!$
!!$    real(DP), intent(in ) :: Ls
!!$    real(DP), intent(in ) :: Lat
!!$    real(DP), intent(out) :: DOD
!!$
!!$
!!$    ! Local variables
!!$    ! lat2     : Temporary variable for latitude
!!$    ! interc   : Intercept
!!$    ! amp      :
!!$    ! phase    : Phase (degree)
!!$
!!$    real(DP) :: lat2
!!$    real(DP) :: interc
!!$    real(DP) :: amp, phase
!!$
!!$
!!$    if( lat .lt. -40.0d0 ) then
!!$      lat2 = -40.0d0
!!$    else if( lat .gt. 40.0d0 ) then
!!$      lat2 = 40.0d0
!!$    else
!!$      lat2 = lat
!!$    endif
!!$
!!$    if( lat2 .gt. 0.0d0 ) then
!!$      interc = 0.16d0 - 1.8d-3 * lat2
!!$    else
!!$      interc = 0.16d0 + 1.3d-3 * lat2
!!$    endif
!!$    amp   = 0.0623d0 - 0.015d0 * atan( ( lat2 + 2.0d0 ) / 5.0d0 )
!!$    phase = 258.0d0 + 1.8d-1 * lat2
!!$    dod   = interc + amp * cos( ( ls - phase ) * d2r )
!!$
!!$
!!$    ! Dust optical depth at 9 micron is converted to that at 0.67 micron.
!!$
!!$    dod = dod * 2.0d0
!!$
!!$
!!$  end subroutine SetMarsDustDODMGS
!!$
!!$    !**************************************************************************
!!$
!!$    subroutine dodMGS_1yr( ls, lat, dod )
!!$
!!$      use mars_const, only : d2r
!!$
!!$      real(dp), intent(in ) :: ls
!!$      real(dp), intent(in ) :: lat
!!$      real(dp), intent(out) :: dod
!!$
!!$
!!$      ! Local variables
!!$
!!$      real(dp) :: dodds
!!$      real(dp) :: dsls, maxdod, dsdtc
!!$
!!$
!!$      call dodMGS( ls, lat, dod )
!!$
!!$
!!$      !*****Add dust storms
!!$      !-----First year
!!$      dsls   = 227.0d0
!!$      maxdod = 0.475d0
!!$      dsdtc  = 35.0d0
!!$      call duststormexp( ls, dodds, dsls, maxdod, dsdtc )
!!$      dodds  = dodds * exp( -( lat - ( -20.0d0 ) )**2 / ( 60.0d0 )**2 )
!!$      dodds  = dodds * 2.0d0
!!$      dod    = max( dod, dodds )
!!$
!!$      dsls   = 235.0d0
!!$      maxdod = 0.5d0
!!$      dsdtc  = 50.0d0
!!$      call duststormexp( ls, dodds, dsls, maxdod, dsdtc )
!!$      dodds  = dodds * exp( -( lat - ( -60.0d0 ) )**2 / ( 60.0d0 )**2 )
!!$      dodds  = dodds * 2.0d0
!!$      dod    = max( dod, dodds )
!!$
!!$      dsls   = 259.0d0
!!$      maxdod = 0.4d0
!!$      dsdtc  = 70.0d0
!!$      call duststormexp( ls, dodds, dsls, maxdod, dsdtc )
!!$      dodds  = dodds * exp( -( lat - ( -80.0d0 ) )**2 / ( 30.0d0 )**2 )
!!$      dodds  = dodds * 2.0d0
!!$      dod    = max( dod, dodds )
!!$
!!$      !-----Second year
!!$!      dsls   = 360.0d0 + 190.0d0
!!$!      maxdod = 1.7d0
!!$!      dsdtc  = 40.0d0
!!$!      call duststormexp( ls, dodds, dsls, maxdod, dsdtc )
!!$!      dodds  = dodds * exp( -( lat - ( -20.0d0 ) )**2 / ( 60.0d0 )**2 )
!!$!      if( abs( lat ) .gt. 60.0d0 ) dodds = 0.0d0
!!$!      dod    = max( dod, dodds )
!!$
!!$
!!$    end subroutine dodMGS_1yr
!!$
!!$    !**************************************************************************
!!$
!!$    subroutine duststormlin( ls, dod, ls0, x1, y1, x2, y2 )
!!$
!!$      real(dp), intent(in ) :: ls
!!$      real(dp), intent(out) :: dod
!!$      real(dp), intent(in ) :: ls0
!!$      real(dp), intent(in ) :: x1, y1, x2, y2
!!$
!!$
!!$      ! Local variables
!!$      !
!!$      real(dp) :: a, b
!!$      real(dp) :: tmpls
!!$
!!$
!!$      a = ( y2 - y1 ) / ( x2 - x1 )
!!$      b = y1 - ( y2 - y1 ) / ( x2 - x1 ) * x1
!!$
!!$      if( ls .lt. ls0 ) then
!!$         tmpls = ls + 360.0d0
!!$      else
!!$         tmpls = ls
!!$      endif
!!$
!!$      dod = a * tmpls + b
!!$
!!$      dod = max( dod, 0.0d0 )
!!$
!!$
!!$    end subroutine duststormlin

  !**************************************************************************
  ! dustsstormexp
  !**************************************************************************
  ! ls      : Areocentric solar longitude (degree)
  ! dod     : Derived dust optical depth at 0.67 micron
  ! dsls    : Areocentric solar longitude at the initiation of dust storm
  !         : (degree)
  ! maxdod  : Maximum value of dust optical depth at 0.67 micron
  ! dsdtc   : Decay time constant of dust storm in unit of areocentric
  !         : solar longitude (degree)
  !**************************************************************************

  subroutine SetMarsDustDSExp( &
    & Ls, DSLs, MaxDOD, DSDTC, & ! (in)
    & DOD                      & ! (out)
    & )

    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

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

  subroutine SetMarsDustRegDSExp( &
    & Ls, DSLs, MaxDOD, DSDTC, & ! (in)
    & xy_DOD                   & ! (out)
    & )

    ! 物理・数学定数設定
    ! Physical and mathematical constants settings
    !
    use constants0, only: &
      & PI                    ! $ \pi $.
                              ! 円周率. Circular constant

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only: &
      & x_Lon, & !
      & y_Lat    ! $ \varphi $ [rad.] . 緯度. Latitude

    real(DP), intent(in ) :: Ls
    real(DP), intent(in ) :: DSLs
    real(DP), intent(in ) :: MaxDOD
    real(DP), intent(in ) :: DSDTC
    real(DP), intent(out) :: xy_DOD(0:imax-1, 1:jmax)


    ! Local variables
    !
    real(DP) :: TMPLs
    integer  :: i
    integer  :: j


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

    do j = 1, jmax
      if ( ( -75.0_DP * PI / 180.0_DP <= y_Lat(j) ) .and. &
        &  ( y_Lat(j) <= -15.0_DP * PI / 180.0_DP ) ) then
        do i = 0, imax-1
          if ( ( 300.0_DP * PI / 180.0_DP <= x_Lon(i) ) .or. &
            &  ( x_Lon(i) <= 60.0_DP * PI / 180.0_DP  ) ) then
            xy_DOD(i,j) = 1.0_DP
          else
            xy_DOD(i,j) = 0.0_DP
          end if
        end do
      else
        xy_DOD(:,j) = 0.0_DP
      end if
    end do

    xy_DOD = xy_DOD * MaxDOD * exp( -( TMPLs - DSLs ) / DSDTC )


  end subroutine SetMarsDustRegDSExp


    !**************************************************************************
!!$
!!$    subroutine mars_setdust_vdist( gph, gp, grho, dod610, qdust, ijs, ije )
!!$
!!$      use maparam   , only : im => imax, jm => jmax, km => kmax
!!$      use maconst   , only : grav
!!$      use mars_const, only : pi
!!$
!!$      real(dp)    , intent(in ) :: gph( im, jm, km+1 ), gp( im, jm, km ), &
!!$           grho( im, jm, km )
!!$      real(dp)    , intent(in ) :: dod610( im, jm )
!!$      real(dp)    , intent(out) :: qdust ( im, jm, km )
!!$      integer(i4b), intent(in ) :: ijs, ije
!!$
!!$
!!$      !
!!$      ! local variables
!!$      !
!!$      ! dod067   : Dust optical depth at 0.67 micron meter
!!$      !          : This is a local variable.
!!$      ! qdust0   : Constant for Use of Dust Mixing Ratio profile
!!$      !          : profile is obtained from Conrath [1975]
!!$      ! qdust1   : Constant for Use of Dust Mixing Ratio profile
!!$
!!$      real(dp) :: dod067( im, jm, km+1 )
!!$      real(dp) :: qdust0, qdust1
!!$
!!$      ! refp       : reference pressure (refp is set to 610 Pa)
!!$      ! p0         : reference pressure (p0 is set to 610 Pa)
!!$      !
!!$      real(dp)     :: refp = 610.0d0, p0 = 610.0d0
!!$
!!$      real(dp)     :: dodtmp
!!$
!!$      integer(i4b) :: ij, k
!!$
!!$
!!$      qdust0=1.0d0
!!$!      qdust1=0.007d0
!!$!      qdust1=0.03d0
!!$      qdust1 = dust_nu_coef
!!$
!!$      do k = 1, km
!!$         do ij = ijs, ije
!!$            qdust( ij, 1, k ) = qdust0 &
!!$                 * exp( qdust1 * ( 1.0d0 - ( p0 / gp( ij, 1, k ) ) ) )
!!$         end do
!!$      end do
!!$
!!$      call calcdod067( gph, grho, qdust, dod067, ijs, ije )
!!$
!!$      do ij = ijs, ije
!!$         dodtmp = dod067( ij, 1, km+1 ) * refp / gph( ij, 1, km+1 )
!!$         qdust0 = 1.0d0
!!$         qdust0 = qdust0 * dod610( ij, 1 ) / dodtmp
!!$         do k = 1, km
!!$            qdust( ij, 1, k ) = qdust( ij, 1, k ) * qdust0
!!$         end do
!!$         do k = 1, km+1
!!$            dod067( ij, 1, k ) = dod067( ij, 1, k ) * qdust0
!!$         end do
!!$      end do
!!$
!!$
!!$    end subroutine mars_setdust_vdist
!!$
  !--------------------------------------------------------------------------------------

  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

    ! ヒストリデータ出力
    ! History data output
    !
    use gtool_historyauto, only: HistoryAutoAddVariable


    ! 宣言文 ; 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/ &
      & DustExtEff,         &
      & REff,               &
      & RhoDust,            &
      & DustScenario,       &
      & DODFileName,        &
      & DODVarName,         &
      & DOD067,             &
      & DustVerDistCoef,    &
      & DustOptDepRefPress, &
      & DustVerDistRefPress
          !
          ! デフォルト値については初期化手続 "rad_Mars_V1#RadMarsV1Init"
          ! のソースコードを参照のこと.
          !
          ! Refer to source codes in the initialization procedure
          ! "rad_Mars_V1#RadMarsV1Init" for the default values.
          !


    ! デフォルト値の設定
    ! Default values settings
    !

    DustExtEff = 3.04_DP   ! Ockert-Bell et al. (1997)
    REff       = 1.85d-6   ! Ockert-Bell et al. (1997)
    RhoDust    = 2500.0_DP ! Pettengill and Ford (2000)


    DustScenario    = 'Const'

    DODFileName     = ''
    DODVarName      = ''

    DOD067          = 0.2_DP
!!$    DustVerDistCoef = 0.01_DP
    DustVerDistCoef = 0.007_DP

!!$    DustOptDepRefPress  = 610.0_DP
!!$    DustVerDistRefPress = 610.0_DP
    DustOptDepRefPress  = 700.0_DP
    DustVerDistRefPress = 700.0_DP

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

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

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


    if ( DustScenario == 'Const' ) then
      IDDustScenario = IDDustScenarioConst
    else if ( DustScenario == 'VikingNoDS' ) then
      IDDustScenario = IDDustScenarioVikingNoDS
    else if ( DustScenario == 'Viking' ) then
      IDDustScenario = IDDustScenarioViking
    else if ( DustScenario == 'MGS' ) then
      IDDustScenario = IDDustScenarioMGS
    else if ( DustScenario == 'MGSDODFromFile' ) then
      IDDustScenario = IDDustScenarioMGSDODFromFile
    else
      call MessageNotify( 'E', module_name, 'DustScenario of %c is not supported.', c1 = trim( DustScenario ) )
    end if



    ! Initialization of modules used in this module
    !


    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'DustPresc' , &
      & (/ 'lon ', 'lat ', 'sig ', 'time'/),   &
      & 'DustPresc', '1' )

    call HistoryAutoAddVariable( 'DustMaxHeight' , &
      & (/ 'lon ', 'lat ', 'time'/),               &
      & 'DustMaxHeight', 'm' )



    ! 印字 ; Print
    !
    call MessageNotify( 'M', module_name, '----- Initialization Messages -----' )
    call MessageNotify( 'M', module_name, 'DustExtEff          = %f', d  = (/ DustExtEff /) )
    call MessageNotify( 'M', module_name, 'REff                = %f', d  = (/ REff /) )
    call MessageNotify( 'M', module_name, 'RhoDust             = %f', d  = (/ RhoDust /) )
    call MessageNotify( 'M', module_name, 'DustScenario        = %c', c1 = trim( DustScenario ) )
    call MessageNotify( 'M', module_name, 'DODFileName         = %c', c1 = trim( DODFileName ) )
    call MessageNotify( 'M', module_name, 'DODVarName          = %c', c1 = trim( DODVarName ) )
    call MessageNotify( 'M', module_name, 'DOD067              = %f', d  = (/ DOD067      /) )
    call MessageNotify( 'M', module_name, 'DustVerDistCoef     = %f', d  = (/ DustVerDistCoef /) )
    call MessageNotify( 'M', module_name, 'DustOptDepRefPress  = %f', d  = (/ DustOptDepRefPress /) )
    call MessageNotify( 'M', module_name, 'DustVerDistRefPress = %f', d  = (/ DustVerDistRefPress /) )
    call MessageNotify( 'M', module_name, '-- version = %c', c1 = trim(version) )

    set_Mars_dust_inited = .true.

  end subroutine SetMarsDustInit

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

end module set_Mars_dust