!= 重力沈降過程
!
!= Gravitational sedimentation process
!
! Authors::   Yoshiyuki O. TAKAHASHI
! Version::   $Id: grav_sed.f90,v 1.3 2014/05/07 09:39:24 murashin Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. 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 :: GravSedAdjust
  public :: GravSedTransport1D
  public :: GravSedCalcSedVel
  public :: GravSedCalcSedVel1D
  public :: GravSedInit


  ! 公開変数
  ! Public variables
  !

  ! 非公開変数
  ! Private variables
  !

  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 2014/05/07 09:39:24 murashin 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 GravSedAdjust(              &
    & xyr_Press, xyr_Height,             & ! (in   )
    & xyz_PartDen, xyz_PartRad,          & ! (in   )
    & xyz_QMix,                          & ! (inout)
    & xyr_GravSedMassFlux                & ! (out  ) optional
    & )

    ! 時刻管理
    ! 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

    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(in   ) :: xyz_PartDen(0:imax-1, 1:jmax, 1:kmax)
    real(DP)    , intent(in   ) :: xyz_PartRad(0:imax-1, 1:jmax, 1:kmax)
    real(DP)    , intent(inout) :: xyz_QMix   (0:imax-1, 1:jmax, 1:kmax)
    real(DP)    , intent(out  ), optional :: xyr_GravSedMassFlux(0:imax-1, 1:jmax, 0:kmax)


    !
    ! local variables
    !
    real(DP) :: xyr_PartRad(0:imax-1, 1:jmax, 0:kmax)
    real(DP) :: xyr_PartDen(0:imax-1, 1:jmax, 0:kmax)

    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

    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
    !
    xyr_PartRad(:,:,0:kmax-1) = xyz_PartRad(:,:,1:kmax)
    xyr_PartRad(:,:,kmax    ) = xyz_PartRad(:,:,kmax)
    xyr_PartDen(:,:,0:kmax-1) = xyz_PartDen(:,:,1:kmax)
    xyr_PartDen(:,:,kmax    ) = xyr_PartDen(:,:,kmax)

    call GravSedCalcSedVel(          &
      & imax, jmax, kmax+1,          & ! (in )
      & xyr_Press,                   & ! (in )
      & xyr_PartDen, xyr_PartRad,    & ! (in )
      & xyr_SedVel                   & ! (out)
      & )


    ! 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( xyr_GravSedMassFlux ) ) then
      xyr_GravSedMassFlux = xyr_QMixFlux
    end if


  end subroutine GravSedAdjust

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

  subroutine GravSedTransport1D( &
    & TimeStep,                                     & ! (in )
    & kmax, r_Press, r_Height, r_SedVel, z_DelMass, & ! (in )
    & z_DelMassA,                                   & ! (out)
    & r_GravSedMassFlux                             & ! (out)
    & )


    real(DP), intent(in ) :: TimeStep
    integer , intent(in ) :: kmax
    real(DP), intent(in ) :: r_Press          (0:kmax)
    real(DP), intent(in ) :: r_Height         (0:kmax)
    real(DP), intent(in ) :: r_SedVel         (0:kmax)
    real(DP), intent(in ) :: z_DelMass        (1:kmax)
    real(DP), intent(out) :: z_DelMassA       (1:kmax)
    real(DP), intent(out) :: r_GravSedMassFlux(0:kmax)


    ! Local variables
    !
    real(dp) :: SedVel

    real(dp) :: r_APHeight    (0:kmax)
    real(dp) :: r_DPHeight    (0:kmax)
    integer  :: r_DPIndex     (0:kmax)
    real(dp) :: r_DPPress     (0:kmax)
    real(dp) :: PressBot
    real(dp) :: PressTop
    real(dp) :: HeightBot
    real(dp) :: HeightTop
    real(dp) :: DPHeight
    real(dp) :: ScaleHeight

    integer  :: k
    integer  :: kk


    ! Assumption:
    !   * Sedimentation velocity has linear subgrid distribution.
    !   * Sedimentation velocity is evaluated at center of arrival point
    !     and departure point.
    !
    ! w = ( w_{k} - w_{k-1} ) / ( z_{k} - z_{k-1} ) * ( z_{mp} - z_{k-1} ) + w_{k-1}
    !   = ( w_{k} - w_{k-1} ) / ( z_{k} - z_{k-1} ) * ( ( z_{ap} + z_{dp} ) / 2 - z_{k-1} ) + w_{k-1}
    !   = ( w_{k} - w_{k-1} ) / ( z_{k} - z_{k-1} ) * ( z_{ap} / 2 - z_{k-1} ) + w_{k-1}
    !   + ( w_{k} - w_{k-1} ) / ( z_{k} - z_{k-1} ) * z_{dp} / 2
    !
    ! z_{dp} = z_{ap} - w * dt
    !
    ! w = ( w_{k} - w_{k-1} ) / ( z_{k} - z_{k-1} ) * ( z_{ap} / 2 - z_{k-1} ) + w_{k-1}
    !   + ( w_{k} - w_{k-1} ) / ( z_{k} - z_{k-1} ) * ( z_{ap} - w * dt ) / 2
    !
    ! ( 1 + ( w_{k} - w_{k-1} ) / ( z_{k} - z_{k-1} ) * dt / 2 ) * w
    !   = ( w_{k} - w_{k-1} ) / ( z_{k} - z_{k-1} ) * ( z_{ap} / 2 - z_{k-1} ) + w_{k-1}
    !   + ( w_{k} - w_{k-1} ) / ( z_{k} - z_{k-1} ) * z_{ap} / 2
    !   = ( w_{k} - w_{k-1} ) / ( z_{k} - z_{k-1} ) * ( z_{ap} - z_{k-1} ) + w_{k-1}
    !   = ( w_{k} - w_{k-1} ) / ( z_{k} - z_{k-1} ) * ( z_{ap} - z_{k-1} )
    !             + w_{k-1}   / ( z_{k} - z_{k-1} ) * ( z_{k}  - z_{k-1} )
    !   = ( ( w_{k} - w_{k-1} ) * ( z_{ap} - z_{k-1} )
    !               + w_{k-1}   * ( z_{k}  - z_{k-1} ) )
    !     / ( z_{k} - z_{k-1} )
    !   = ( w_{k} * ( z_{ap} - z_{k-1} ) - w_{k-1} * ( z_{ap} - z_{k-1} )
    !                                    + w_{k-1} * ( z_{k}  - z_{k-1} ) )
    !     / ( z_{k} - z_{k-1} )
    !   = ( w_{k} * ( z_{ap} - z_{k-1} ) - w_{k-1} * z_{ap}
    !                                    + w_{k-1} * z_{k}  )
    !     / ( z_{k} - z_{k-1} )
    !   = ( w_{k} * ( z_{ap} - z_{k-1} ) - w_{k-1} * ( z_{ap} - z_{k} ) )
    !     / ( z_{k} - z_{k-1} )
    !   = ( w_{k} * ( z_{ap} - z_{k-1} ) + w_{k-1} * ( z_{k} - z_{ap} ) )
    !     / ( z_{k} - z_{k-1} )


    r_APHeight = r_Height
    do k = 0, kmax-1
      ! MEMO: r_SedVel is negative, i.e., downward.
      !
      ! Simple method
!!$      SedVel = r_SedVel(k)
      !
      ! Linear method
      !     ( 1 + ( w_{k} - w_{k-1} ) / ( z_{k} - z_{k-1} ) * dt / 2 ) * w
      !   = ( w_{k} * ( z_{ap} - z_{k-1} ) + w_{k-1} * ( z_{k} - z_{ap} ) )
      !     / ( z_{k} - z_{k-1} )
      SedVel = &
        &   (   r_SedVel(k+1) * ( r_APHeight(k  ) - r_Height  (k  ) )   &
        &     + r_SedVel(k  ) * ( r_Height  (k+1) - r_APHeight(k  ) ) ) &
        & / ( r_Height(k+1) - r_Height(k) )                             &
        & / ( 1.0d0 +   ( r_SedVel(k+1) - r_SedVel(k) )                 &
        &             / ( r_Height(k+1) - r_Height(k) ) * TimeStep / 2.0d0 )
      !
      r_DPHeight(k) = r_APHeight(k) - SedVel * TimeStep
      r_DPHeight(k) = min( r_DPHeight(k), r_Height(kmax) )
    end do
    k = kmax
    r_DPHeight(k) = r_Height(kmax)

    ! search index
    !   Departure point is located in a r_DPIndex'th layer.
    !   Here, r_DPIndex is layer index, which ranges from 1 to kmax.
    do k = 0, kmax-1
      do kk = k+1, kmax
        if ( r_DPHeight(k) <= r_Height(kk) ) exit
      end do
      r_DPIndex(k) = kk
    end do
    k = kmax
    r_DPIndex(k) = -1

    ! p_{k} = p_{k-1} * exp( - dz / H )
    ! ln(p_{k}/p_{k-1}) = - dz / H
    ! H = - dz / ln(p_{k}/p_{k-1})
    !   = - (z_{k}-z_{k-1}) / ln(p_{k}/p_{k-1})
    do k = 0, kmax-1
      kk = r_DPIndex(k)
      !
      if ( r_Press(kk) <= 0.0d0 ) then
        PressTop = r_Press(kk-1) * 1.0d-1
      else
        PressTop = r_Press(kk  )
      end if
      PressBot  = r_Press (kk-1)
      HeightTop = r_Height(kk  )
      HeightBot = r_Height(kk-1)
      DPHeight  = r_DPHeight(k)
      !
      ScaleHeight = - ( HeightTop - HeightBot ) / log( PressTop / PressBot )
      r_DPPress(k) = PressBot * exp( - ( DPHeight - HeightBot ) / ScaleHeight )
    end do
    k = kmax
    r_DPPress(k) = r_Press(kmax)
    !
    do k = 1, kmax
      z_DelMassA(k) = 0.0d0
      !
      if ( k < kmax ) then
        kk = r_DPIndex(k  )
        z_DelMassA(k) = z_DelMassA(k) &
          & + z_DelMass(kk) &
          &   * ( r_Press(kk-1) - r_DPPress(k) ) / ( r_Press(kk-1) - r_Press(kk) )
      end if
      !
      do kk = r_DPIndex(k  )-1, r_DPIndex(k-1)+1, -1
        z_DelMassA(k) = z_DelMassA(k) + z_DelMass(kk)
      end do
      !
      kk = r_DPIndex(k-1)
      z_DelMassA(k) = z_DelMassA(k) &
        & + z_DelMass(kk) &
        &   * ( r_DPPress(k-1) - r_Press(kk) ) / ( r_Press(kk-1) - r_Press(kk) )
    end do

    ! mass flux at layer interface
    r_GravSedMassFlux = 0.0_DP
    do k = 0, kmax-1
      if ( k < kmax ) then
        kk = r_DPIndex(k  )
        r_GravSedMassFlux(k) = r_GravSedMassFlux(k) &
          & + z_DelMass(kk) &
          &   * ( r_Press(kk-1) - r_DPPress(k) ) / ( r_Press(kk-1) - r_Press(kk) )
      end if
      do kk = r_DPIndex(k)-1, k+1, -1
        r_GravSedMassFlux(k) = r_GravSedMassFlux(k) + z_DelMass(kk)
      end do
    end do
    k = kmax
    r_GravSedMassFlux(kmax) = 0.0_DP

    ! sedimentation flux is negative / downward
    r_GravSedMassFlux = - r_GravSedMassFlux


  end subroutine GravSedTransport1D

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

  subroutine GravSedCalcSedVel(          &
    & imax, jmax, kmax,                  &
    & xyz_Press,                         & ! (in )
    & xyz_PartDen, xyz_PartRad,          & ! (in )
    & xyz_SedVel                         & ! (out)
    & )


    integer     , intent(in   ) :: imax
    integer     , intent(in   ) :: jmax
    integer     , intent(in   ) :: kmax
    real(DP)    , intent(in   ) :: xyz_Press  (0:imax-1, 1:jmax, 1:kmax)
    real(DP)    , intent(in   ) :: xyz_PartDen(0:imax-1, 1:jmax, 1:kmax)
    real(DP)    , intent(in   ) :: xyz_PartRad(0:imax-1, 1:jmax, 1:kmax)
    real(DP)    , intent(out  ) :: xyz_SedVel (0:imax-1, 1:jmax, 1:kmax)


    !
    ! local variables
    !
    real(DP) :: MolVisCoef

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



    real(DP) :: xyz_PartDia    (0:imax-1, 1:jmax, 1:kmax)

    real(DP) :: PlanetLonFromVE
    real(DP) :: xyr_DOD067     (0:imax-1, 1:jmax, 0: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


    !
    ! 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


    xyz_PartDia = 2.0_DP * xyz_PartRad

    xyz_SedVel = &
      & aaa_SedVel( &
      &   MolVisCoef, MeanFreePathRef, PressLambdaRef, xyz_PartDen, xyz_PartDia, &
      &   max( xyz_Press, 1.0e-20_DP )                                       &
      & )


  end subroutine GravSedCalcSedVel

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

  subroutine GravSedCalcSedVel1D(          &
    & kmax, z_Press, z_PartDen, z_PartRad, & ! (in )
    & z_SedVel                             & ! (out)
    & )

    integer     , intent(in   ) :: kmax
    real(DP)    , intent(in   ) :: z_Press  (1:kmax)
    real(DP)    , intent(in   ) :: z_PartDen(1:kmax)
    real(DP)    , intent(in   ) :: z_PartRad(1:kmax)
    real(DP)    , intent(out  ) :: z_SedVel (1:kmax)


    !
    ! local variables
    !
    real(DP) :: z_PartDia(1:kmax)
    real(DP) :: MolVisCoef

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

    integer  :: k


    ! 実行文 ; Executable statement
    !

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

    !
    ! 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


    z_PartDia = 2.0_DP * z_PartRad

    z_SedVel = &
      & a_SedVel( &
      &   MolVisCoef, MeanFreePathRef, PressLambdaRef, z_PartDen, z_PartDia, &
      &   max( z_Press, 1.0e-20_DP )                                         &
      & )


  end subroutine GravSedCalcSedVel1D

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

  function aaa_SedVel(                                                       &
    & MolVisCoef, MeanFreePathRef, PressLambdaRef, aaa_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) :: aaa_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 =                                                      &
      & - aaa_PartDen * Grav * aaa_PartDia / ( 18.0_DP * MolVisCoef ) &
      & * ( aaa_PartDia + 2.0_DP * aaa_MeanFreePath )


  end function aaa_SedVel

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

  function a_SedVel(                                                     &
    & MolVisCoef, MeanFreePathRef, PressLambdaRef, a_PartDen, a_PartDia, &
    & a_Press                                                            &
    & ) result( a_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) :: a_PartDen(:)
    real(DP), intent(in) :: a_PartDia(:)
    real(DP), intent(in) :: a_Press  (:)

    real(DP) :: a_Result(size(a_Press,1))

    !
    ! local variables
    !
    real(DP) :: a_MeanFreePath(size(a_Press,1))


    ! 実行文 ; Executable statement
    !

    a_MeanFreePath = MeanFreePathRef * ( PressLambdaRef / a_Press )
!!$    aaa_Result =                                                     &
!!$      & - PartDen * Grav * aaa_PartDia**2 / ( 18.0_DP * MolVisCoef ) &
!!$      & * ( 1.0_DP + 2.0_DP * aaa_MeanFreePath / aaa_PartDia )
    a_Result =                                                    &
      & - a_PartDen * Grav * a_PartDia / ( 18.0_DP * MolVisCoef ) &
      & * ( a_PartDia + 2.0_DP * a_MeanFreePath )


  end function a_SedVel

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

  function aaa_SedVelStokesSimple(              &
    & aaa_MolVisCoef, aaa_PartDen, aaa_PartDia  &
    & ) result( aaa_Result )

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


    real(DP), intent(in) :: aaa_MolVisCoef(:,:,:)
    real(DP), intent(in) :: aaa_PartDen   (:,:,:)
    real(DP), intent(in) :: aaa_PartDia   (:,:,:)

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

    !
    ! local variables
    !


    ! 実行文 ; Executable statement
    !

    aaa_Result =                                                      &
      & - aaa_PartDen * Grav * aaa_PartDia**2 / ( 18.0_DP * aaa_MolVisCoef )


  end function aaa_SedVelStokesSimple

  !--------------------------------------------------------------------------------------
!!$  !
!!$  ! 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

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


    ! 宣言文 ; 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/    &
!!$      & IceNumRatio,           &



    ! 実行文 ; Executable statement
    !

    if ( grav_sed_inited ) return


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


    ! 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 ) write( STDOUT, nml = grav_sed_nml )
!!$    end if



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

    grav_sed_inited = .true.

  end subroutine GravSedInit

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

end module grav_sed
