!= 重力沈降過程
!
!= Gravitational sedimentation process
!
! Authors::   Yoshiyuki O. TAKAHASHI, Shin-ichi Takehiro
! Version::   $Id: grav_sed.f90,v 1.3 2025/09/17 22:00:00 takepiro Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008-2025. All rights reserved.
! License::   See COPYRIGHT[link:../../../COPYRIGHT]
!
module grav_sed
  !
  != 重力沈降過程
  !
  != Gravitational sedimentation process
  !
  ! <b>Note that Japanese and English are described in parallel.</b>
  !
  ! 重力沈降過程を計算するモジュールです.
  !
  ! This module calculates gravitational sedimentation. 
  !
  !== Procedures List
  !
  ! GravSed         :: 計算
  ! GravSedInit     :: 初期化
  ! --------------- :: ------------
  ! GravSed         :: Calculation
  ! GravSedInit     :: Initialization
  !
  !== NAMELIST
  !
  ! NAMELIST#grav_sed_nml
  !
  !== References
  !
  ! * Conrath, B. J., 1975:
  !   Thermal structure of the Martian atmosphere during the dissipation of 
  !   the dust storm of 1971, 
  !   <i>Icarus</i>, <b>24</b>, 36--46.
  !
  ! * Lin, S.-J., and R. B. Rood, 1996:
  !   Multidimensional flux-form semi-Lagrangian transport scheme, 
  !   <i>Mon. Wea. Rev.</i>, <b>124</b>, 2046--2070.
  !

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

  ! 種別型パラメタ
  ! Kind type parameter
  !
  use dc_types, only: DP, &      ! 倍精度実数型. Double precision.
    &                 STRING     ! 文字列.       Strings.

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

  ! 組成に関わる配列の設定
  ! Settings of array for atmospheric composition
  !
  use composition, only: ncmax

  ! 格子点設定
  ! 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


  ! 宣言文 ; Declaration statements
  !
  implicit none
  private

  ! 公開手続き
  ! Public procedure
  !
  public :: GravSed
  public :: GravSedInit


  ! 公開変数
  ! Public variables
  !

  ! 非公開変数
  ! Private variables
  !
  real(DP), save :: RadiusMarsH2OCloud
  real(DP), save :: IceNumRatio
  real(DP), save :: DOD067ForMarsH2OCloud

  logical, save :: grav_sed_inited = .false.
                              ! 初期設定フラグ.
                              ! Initialization flag


  character(*), parameter:: module_name = 'grav_sed'
                              ! モジュールの名称.
                              ! Module name
  character(*), parameter:: version = &
    & '$Name:  $' // &
    & '$Id: grav_sed.f90,v 1.3 2025/09/17 22:00:00 takepiro Exp $'
                              ! モジュールのバージョン
                              ! Module version


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

contains

  !------------------------------------------------------------------------------------
  ! Gravitational sedimentation is calculated.
  ! This routine works as adjustment routine, i.e., mass mixing ratio is 
  ! updated in this routine. 
  ! Current version calculates sedimentation of dust, only. 
  ! In this routine, gravitational sedimentation is calculated by the use of 
  ! a method following flux-form semi-Lagrangian transport scheme (Lin and 
  ! Rood, 1996). 
  !

  subroutine GravSed(                    &
    & SpcName,                           & ! (in )
    & xyr_Press, xyr_Height,             & ! (in )
    & xyz_QMix,                          & ! (out)
    & xy_SurfGravSedFlux                 & ! (out) optional
    & )

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

    ! 時刻管理
    ! Time control
    !
    use timeset, only: &
      & TimeN, &              ! ステップ $ t $ の時刻. Time of step $ t $.
      & TimesetClockStart, TimesetClockStop, &
      & DelTime                 ! $ \Delta t $ [s]

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

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & Grav
                              ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration


    character(*), intent(in   ) :: SpcName
    real(DP)    , intent(in   ) :: xyr_Press  (0:imax-1, 1:jmax, 0:kmax)
    real(DP)    , intent(in   ) :: xyr_Height (0:imax-1, 1:jmax, 0:kmax)
    real(DP)    , intent(inout) :: xyz_QMix   (0:imax-1, 1:jmax, 1:kmax)
    real(DP)    , intent(out  ), optional :: xy_SurfGravSedFlux(0:imax-1, 1:jmax)


    !
    ! local variables
    !
    ! rhod : dust density
    ! mfp  : mean free path
    ! mvc  : modecular viscosity coefficient
    ! rdia : particle diameter
    !
    real(DP) :: PartDen
    real(DP) :: xyr_PartDia(0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: MolVisCoef

    real(DP) :: MeanFreePathRef
    real(DP) :: PressLambdaRef


    real(DP) :: xyz_DelAtmMass  (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_DelCompMass (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_DelZ        (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyr_SedVel      (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_FracSedDist (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_Dist        (0:imax-1, 1:jmax, 0:kmax)
    integer  :: xyr_KIndex      (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_QMixFlux    (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_IntQMixFlux (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_FracQMixFlux(0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyz_DQMixDt     (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: xyz_QMixA       (0:imax-1, 1:jmax, 1:kmax)
    real(DP) :: LogPress
    real(DP) :: Press

    real(DP), parameter :: AMU = 1.6605655e-27_DP

    real(DP) :: DustExtEff
    real(DP) :: DustREff
    real(DP) :: DOD067Ref
    real(DP) :: MeanMolMass
    real(DP) :: PressRef
    real(DP) :: DustNumRatio
!!$    real(DP) :: xyr_NumDensDust(0:imax-1, 1:jmax, 0:kmax)
!!$    real(DP) :: xyr_NumDensIce (0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyz_PartDia    (0:imax-1, 1:jmax, 1:kmax)

    integer  :: i
    integer  :: j
    integer  :: k
    integer  :: kk


    ! 実行文 ; Executable statement
    !

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



    ! Calculation of mass in each layer and layer thickness in unit of meter
    !   Layer thickness is calculated by using mass of a layer.
!!$    xyz_Rho = xyz_Press / ( GasRDry * xyz_VirTemp )
    do k = 1, kmax
      xyz_DelAtmMass(:,:,k) = ( xyr_Press(:,:,k-1) - xyr_Press(:,:,k) ) / Grav
    end do
!!$    xyz_DelZ = xyz_DelAtmMass / xyz_Rho
    do k = 1, kmax
      xyz_DelZ(:,:,k) = xyr_Height(:,:,k) - xyr_Height(:,:,k-1)
    end do


    ! Calculation of mass of constituents in a layer
    xyz_DelCompMass = xyz_QMix * xyz_DelAtmMass


    !
    ! calculation of sedimentation terminal velocity
    !
    if ( SpcName == 'MarsDust' ) then

      !
      ! The values below are obtained from Conrath (1975). 
      ! Particle radius of 1e-6 m is assumed. 
      !
      MolVisCoef      = 1.5e-4_DP * 1.0e-3_DP * 1.0e2_DP
      MeanFreePathRef = 2.2e-4_DP * 1.0e-2_DP
      PressLambdaRef  = 25.0e2_DP
      PartDen         = 3.0e3_DP
      xyr_PartDia     = 2.0_DP * 1.0e-6_DP

      xyr_SedVel = &
        & aaa_SedVel( &
        &   MolVisCoef, MeanFreePathRef, PressLambdaRef, PartDen, xyr_PartDia, &
        &   max( xyr_Press, 1.0e-20_DP )                                       &
        & )
      k = kmax
      xyr_SedVel(:,:,k) = 0.0_DP

    else if ( SpcName == 'MarsH2OCloud' ) then

      !
      ! The values below are obtained from Conrath (1975). 
      !
      MolVisCoef      = 1.5e-4_DP * 1.0e-3_DP * 1.0e2_DP
      MeanFreePathRef = 2.2e-4_DP * 1.0e-2_DP
      PressLambdaRef  = 25.0e2_DP
      !
      ! Particle radius of 1e-6 m is assumed. 
      !
      PartDen         = 1.0e3_DP
!!$      PartDia         = 2.0_DP * 1.0e-6_DP
      if ( RadiusMarsH2OCloud > 0.0_DP ) then
        xyr_PartDia         = 2.0_DP * RadiusMarsH2OCloud
      else
        MeanMolMass  = 44.0_DP * AMU
        ! DOD067Ref : Dust optical depth at 0.67 um
        if ( IceNumRatio < 0.0_DP ) then
          DustExtEff   = 3.04_DP    ! Ockert-Bell et al. (1997)
          DustREff     = 1.85e-6_DP ! Ockert-Bell et al. (1997)
          DOD067Ref    = DOD067ForMarsH2OCloud
          PressRef     = 700.0_DP
          ! DustNumRatio : numb. dens. of dust / num. dens. of atm. molecules.
          DustNumRatio = DOD067Ref * MeanMolMass &
            & / ( DustExtEff * PI * DustREff**2 * PressRef / Grav )
          IceNumRatio = DustNumRatio
!!$        xyz_NumDensDust = DustNumRatio * xyz_Rho / MeanMolMass
!!$        xyz_NumDensIce  = xyz_NumDensDust
        end if
        ! calculate radius, first
!!$        xyz_PartDia = &
!!$          & (                                                       &
!!$          &     xyz_QMix * xyz_Rho                                  &
!!$          &   / ( xyz_NumDensIce * PartDen * 4.0_DP / 3.0_DP * PI ) &
!!$          & )**(1.0_DP/3.0_DP)
        xyz_PartDia =                                               &
          & (                                                       &
          &     xyz_QMix * MeanMolMass                              &
          &   / ( IceNumRatio * PartDen * 4.0_DP / 3.0_DP * PI )    &
          & )**(1.0_DP/3.0_DP)
        ! calculate diameter
        xyz_PartDia = 2.0_DP * xyz_PartDia
        !
        xyr_PartDia(:,:,0:kmax-1) = xyz_PartDia(:,:,1:kmax)
        xyr_PartDia(:,:,kmax) = 0.0_DP
      end if

      xyr_SedVel = &
        & aaa_SedVel(                                                          &
        &   MolVisCoef, MeanFreePathRef, PressLambdaRef, PartDen, xyr_PartDia, &
        &   max( xyr_Press, 1.0e-20_DP )                                       &
        & )
      k = kmax
      xyr_SedVel(:,:,k) = 0.0_DP

    else
      call MessageNotify( 'E', module_name, 'Specie %c is inappropriate', c1 = trim( SpcName ) )
    end if


    ! Calculation of sedimentation distance during a time step of 2 * DelTime
    xyr_Dist = abs( xyr_SedVel ) * 2.0_DP * DelTime
    do k = 0, kmax-1
      do j = 1, jmax
        do i = 0, imax-1

          ! A k index in which all mass of the layer does not fall is 
          ! searched. In addition, fractional sedimentation velocity is 
          ! calculated. 
          xyr_KIndex(i,j,k) = -1
          do kk = k+1, kmax-1
            ! If sedimentation velocity (distance) is positive, and all of 
            ! mass in kk layer does not fall, KIndex is kk.
            if ( ( xyr_Dist(i,j,k) >= 0.0_DP ) .and. &
              &  ( xyr_Dist(i,j,k) <= xyz_DelZ(i,j,kk) ) ) then
              xyr_KIndex     (i,j,k) = kk
              xyr_FracSedDist(i,j,k) = xyr_Dist(i,j,k)
            end if
            ! Sedimentation distance is decreased for preparation for next 
            ! layer.
            ! If xyz_Dist become negative, any mass of the upper layer does 
            ! not fall.
            xyr_Dist(i,j,k) = xyr_Dist(i,j,k) - xyz_DelZ(i,j,kk)
          end do
          ! Calculation for upper most layer.
          kk = kmax
          if ( xyr_Dist(i,j,k) >= 0.0_DP ) then
            xyr_KIndex     (i,j,k) = kk
            xyr_FracSedDist(i,j,k) = min( xyr_Dist(i,j,k), xyz_DelZ(i,j,kk) )
          end if

        end do
      end do
    end do
    ! K index and fractional sedimentation velocity at model top.
    ! No flux is assumed at the model top. 
    k = kmax
    xyr_KIndex     (:,:,k) = -1
    xyr_FracSedDist(:,:,k) = 0.0_DP


    ! Calculation of integer mass flux.
    xyr_IntQMixFlux = 0.0_DP
    do k = 0, kmax-1
      do j = 1, jmax
        do i = 0, imax-1

          do kk = k+1, xyr_KIndex(i,j,k)-1
            xyr_IntQMixFlux(i,j,k) = xyr_IntQMixFlux(i,j,k) &
              & + xyz_DelCompMass(i,j,kk)
          end do
          xyr_IntQMixFlux(i,j,k) = xyr_IntQMixFlux(i,j,k) / ( 2.0_DP * DelTime )

        end do
      end do
    end do

    ! Add sign of sedimentation velocity.
    ! This is equivalent to mulplying -1.
    xyr_IntQMixFlux = sign( 1.0_DP, xyr_SedVel ) * xyr_IntQMixFlux


    ! Calculation of fractional mass flux
    k = kmax
    xyr_FracQMixFlux(:,:,k) = 0.0_DP
    do k = kmax-1, 0, -1
      do j = 1, jmax
        do i = 0, imax-1
          kk = xyr_KIndex(i,j,k)
          !-----
          ! Simple method
!!$            xyrf_FracQMixFlux(i,j,k,n) =                       &
!!$              &   xyrf_FracSedDist(i,j,k,n) / xyz_DelZ(i,j,kk) &
!!$              & * xyzf_DelCompMass(i,j,kk,n)
          !-----
          ! Method considering exponential distribution of mass with height
          if ( xyr_Press(i,j,kk) == 0.0_DP ) then
            LogPress =                                                      &
              &   log( xyr_Press(i,j,kk-1) * 1.0e-1_DP / xyr_Press(i,j,kk-1) ) &
              & / xyz_DelZ(i,j,kk) * xyr_FracSedDist(i,j,k)                 &
              & + log( xyr_Press(i,j,kk-1) )
            Press = exp( LogPress )
            xyr_FracQMixFlux(i,j,k) =                                     &
              &   ( xyr_Press(i,j,kk-1) - Press                        )  &
              & / ( xyr_Press(i,j,kk-1) - xyr_Press(i,j,kk-1) * 1.0e-1_DP )  &
              & * xyz_DelCompMass(i,j,kk)
          else
            LogPress =                                           &
              &   log( xyr_Press(i,j,kk) / xyr_Press(i,j,kk-1) ) &
              & / xyz_DelZ(i,j,kk) * xyr_FracSedDist(i,j,k)      &
              & + log( xyr_Press(i,j,kk-1) )
            Press = exp( LogPress )
            xyr_FracQMixFlux(i,j,k) =                          &
              &   ( xyr_Press(i,j,kk-1) - Press             )  &
              & / ( xyr_Press(i,j,kk-1) - xyr_Press(i,j,kk) )  &
              & * xyz_DelCompMass(i,j,kk)
          end if
          !-----
          xyr_FracQMixFlux(i,j,k) = xyr_FracQMixFlux(i,j,k) &
            & / ( 2.0_DP * DelTime )
        end do
      end do
    end do

    ! Add sign of sedimentation velocity.
    ! This is equivalent to mulplying -1.
    xyr_FracQMixFlux = sign( 1.0_DP, xyr_SedVel ) * xyr_FracQMixFlux


    xyr_QMixFlux = xyr_IntQMixFlux + xyr_FracQMixFlux


    !
    ! estimate dust mixing ratio at next time step
    !
    do k = 1, kmax
      xyz_DQMixDt(:,:,k) =                                  &
        &   ( xyr_QMixFlux(:,:,k) - xyr_QMixFlux(:,:,k-1) ) &
        & / ( xyr_Press   (:,:,k) - xyr_Press   (:,:,k-1) ) &
        & * Grav
    end do


    xyz_QMixA = xyz_QMix + xyz_DQMixDt * 2.0_DP * DelTime

    xyz_QMix  = xyz_QMixA


    if ( present ( xy_SurfGravSedFlux ) ) then
      xy_SurfGravSedFlux = xyr_QMixFlux(:,:,0)
    end if


    ! ヒストリデータ出力
    ! History data output
    !
    if ( SpcName == 'MarsH2OCloud' ) then
      call HistoryAutoPut( TimeN, 'MarsH2OCloudRadius', xyz_PartDia/2.0_DP )
    end if

  end subroutine GravSed

  !------------------------------------------------------------------------------------
  ! Sedimentation velocity is calculated by the use of a formula of Conrath (1975)
  !

  function aaa_SedVel(                                                   &
    & MolVisCoef, MeanFreePathRef, PressLambdaRef, PartDen, aaa_PartDia, &
    & aaa_Press                                                          &
    & ) result( aaa_Result )

    ! 物理定数設定
    ! Physical constants settings
    !
    use constants, only: &
      & Grav
                              ! $ g $ [m s-2].
                              ! 重力加速度.
                              ! Gravitational acceleration


    real(DP), intent(in) :: MolVisCoef
    real(DP), intent(in) :: MeanFreePathRef
    real(DP), intent(in) :: PressLambdaRef
    real(DP), intent(in) :: PartDen
    real(DP), intent(in) :: aaa_PartDia(:,:,:)
    real(DP), intent(in) :: aaa_Press(:,:,:)

    real(DP) :: aaa_Result(size(aaa_Press,1),size(aaa_Press,2),size(aaa_Press,3))

    !
    ! local variables
    !
    real(DP) :: aaa_MeanFreePath(size(aaa_Press,1),size(aaa_Press,2),size(aaa_Press,3))


    ! 実行文 ; Executable statement
    !

    aaa_MeanFreePath = MeanFreePathRef * ( PressLambdaRef / aaa_Press )
!!$    aaa_Result =                                                     &
!!$      & - PartDen * Grav * aaa_PartDia**2 / ( 18.0_DP * MolVisCoef ) &
!!$      & * ( 1.0_DP + 2.0_DP * aaa_MeanFreePath / aaa_PartDia )
    aaa_Result =                                                  &
      & - PartDen * Grav * aaa_PartDia / ( 18.0_DP * MolVisCoef ) &
      & * ( aaa_PartDia + 2.0_DP * aaa_MeanFreePath )


  end function aaa_SedVel

  !--------------------------------------------------------------------------------------
!!$  !
!!$  ! This routine works as an adjustment-type one. 
!!$  !
!!$
!!$  subroutine dust_borrowingfrombelow( damassn, gdmmrn, gdmassg, ijs, ije )
!!$
!!$    use matype
!!$    use maparam, only : imax, jmax, kmax
!!$
!!$    implicit none
!!$
!!$    real(dp)    , intent(in   ) :: damassn( imax, jmax, kmax )
!!$    real(dp)    , intent(inout) :: gdmmrn ( imax, jmax, kmax )
!!$    real(dp)    , intent(inout) :: gdmassg( imax, jmax )
!!$    integer(i4b), intent(in   ) :: ijs, ije
!!$
!!$
!!$    !
!!$    ! local variables
!!$    !
!!$    ! ddm     : dust deficit mass
!!$    !
!!$    real(dp)                :: ddm
!!$
!!$    integer(i4b)            :: ij, k
!!$    integer(i4b), parameter :: j = 1
!!$
!!$
!!$    !
!!$    ! borrowing
!!$    !
!!$    do k = 1, kmax-1
!!$      do ij = ijs, ije
!!$        ddm = -min( gdmmrn( ij, j, k ), 0.0d0 ) * damassn( ij, j, k )
!!$        gdmmrn( ij, j, k   ) = max( gdmmrn( ij, j, k ), 0.0d0 )
!!$        gdmmrn( ij, j, k+1 ) = gdmmrn( ij, j, k+1 ) &
!!$          - ddm / damassn( ij, j, k+1 )
!!$      end do
!!$    end do
!!$    k = kmax
!!$    do ij = ijs, ije
!!$      ddm = -min( gdmmrn( ij, j, k ), 0.0d0 ) * damassn( ij, j, k )
!!$      gdmmrn ( ij, j, k ) = max( gdmmrn( ij, j, k ), 0.0d0 )
!!$      gdmassg( ij, j )    = gdmassg( ij, j ) + ddm
!!$    end do
!!$
!!$
!!$  end subroutine dust_borrowingfrombelow
!!$
  !--------------------------------------------------------------------------------------

  subroutine GravSedInit

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

    ! 種別型パラメタ
    ! Kind type parameter
    !
    use dc_types, only: STDOUT ! 標準出力の装置番号. Unit number of standard output

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

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

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

    ! 座標データ設定
    ! Axes data settings
    !
    use axesset, only : &
      & AxnameX, &
      & AxnameY, &
      & AxnameZ, &
      & AxnameT


    ! メッセージ制御
    ! Message control
    !
    use mpi_messagecntl, only : DoesOutputMPIMessage


    ! 宣言文 ; Declaration statements
    !
    implicit none

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


    ! NAMELIST 変数群
    ! NAMELIST group name
    !
    namelist /grav_sed_nml/    &
      & RadiusMarsH2OCloud,    &
      & IceNumRatio,           &
      & DOD067ForMarsH2OCloud
!!$      & FlagConstDiffCoef,              &
!!$      & ConstDiffCoefM, ConstDiffCoefH, &
!!$!
!!$      & SquareVelMin, BulkRiNumMin,     &
!!$!
!!$      & MixLengthMax, TildeShMin, TildeSmMin, &
!!$      & VelDiffCoefMin, TempDiffCoefMin, QMixDiffCoefMin, &
!!$      & VelDiffCoefMax, TempDiffCoefMax, QMixDiffCoefMax, &
!!$!
!!$      & MYLv2ParamA1, MYLv2ParamB1, MYLv2ParamA2, MYLv2ParamB2, MYLv2ParamC1, &
!!$      & FlagCalcRiWithTv


    ! 実行文 ; Executable statement
    !

    if ( grav_sed_inited ) return


    ! デフォルト値の設定
    ! Default values settings
    !
    RadiusMarsH2OCloud    =  1.0e-6_DP
    IceNumRatio           = -1.0_DP
    DOD067ForMarsH2OCloud =  0.3_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 = grav_sed_nml, &           ! (out)
        & iostat = iostat_nml )           ! (out)
      close( unit_nml )

      call NmlutilMsg( iostat_nml, module_name ) ! (in)
      if ( iostat_nml == 0 .AND. DoesOutputMPIMessage() ) write( STDOUT, nml = grav_sed_nml )
    end if

    ! ヒストリデータ出力のためのへの変数登録
    ! Register of variables for history data output
    !
    call HistoryAutoAddVariable( 'MarsH2OCloudRadius', &
      & (/ AxnameX, AxnameY, AxnameZ, AxnameT /), &
      & 'Mars H2O Cloud Radius', &
      & 'm' )

!!$    do n = 1, ncmax
!!$      call HistoryAutoAddVariable( 'Surf'//trim(a_QMixName(n))//'GravSedFlux', &
!!$        & (/ 'lon ', 'lat ', 'time' /), &
!!$        & 'surface gravitational sedimentation flux of ' // trim(a_QMixName(n)), &
!!$        & 'kg m-2 s-1' )
!!$    end do


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

    grav_sed_inited = .true.

  end subroutine GravSedInit

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

end module grav_sed