subroutine GravSed( SpcName, xyr_Press, xyz_VirTemp, xyr_Height, xyz_QMix, xy_SurfGravSedFlux )
! ヒストリデータ出力
! History data output
!
use gtool_historyauto, only: HistoryAutoPut
! 組成に関わる配列の設定
! Settings of array for atmospheric composition
!
use composition, only: a_QMixName, CompositionInqIndex
! 時刻管理
! Time control
!
use timeset, only: TimeN, TimesetClockStart, TimesetClockStop, DelTime ! $ \Delta t $ [s]
! 物理定数設定
! Physical constants settings
!
use constants, only: Grav, GasRDry
! $ R $ [J kg-1 K-1].
! 乾燥大気の気体定数.
! Gas constant of air
character(*), intent(in ) :: SpcName
real(DP) , intent(in ) :: xyr_Press (0:imax-1, 1:jmax, 0:kmax)
real(DP) , intent(in ) :: xyz_VirTemp(0:imax-1, 1:jmax, 1: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) :: PartDia
real(DP) :: MeanFreePath
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
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
PartDia = 2.0_DP * 1.0e-6_DP
xyr_SedVel = aaa_SedVel( MolVisCoef, MeanFreePathRef, PressLambdaRef, PartDen, 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
xyr_SedVel = aaa_SedVel( MolVisCoef, MeanFreePathRef, PressLambdaRef, PartDen, 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
!
!!$ call HistoryAutoPut( TimeN, 'Surf'//trim(a_QMixName(n))//'GravSedFlux', &
!!$ & xyrf_QMixFlux(:,:,0,n) )
end subroutine GravSed