| Class | ECCM |
| In: |
moist/eccm.f90
|
| Subroutine : | |||
| a_MolFrIni(1:SpcNum) : | real(8), intent(in)
| ||
| Humidity : | real(8), intent(in)
| ||
| z_Temp(DimZMin:DimZMax) : | real(8), intent(out)
| ||
| z_Press(DimZMin:DimZMax) : | real(8), intent(out)
| ||
| z_MolWtMean(DimZMin:DimZMax) : | real(8), intent(out)
| ||
| za_MolFr(DimZMin:DimZMax, 1:SpcNum) : | real(8), intent(out)
|
* 乾燥断熱温度減率に沿った温度、圧力を求め、それらに対して指定された相対湿度となるように凝結成分のモル比を決める * 比熱は乾燥気塊のもので代表させる * 流体の方程式において比熱は乾燥成分のもので代表させているため * 大気の平均分子量には湿潤成分の分子量を効果 * 流体の方程式において, 湿潤成分の分子量は考慮しているため
subroutine ECCM_Dry( a_MolFrIni, Humidity, z_Temp, z_Press, z_MolWtMean, za_MolFr )
!
!== 概要
! * 乾燥断熱温度減率に沿った温度、圧力を求め、それらに対して指定された相対湿度となるように凝結成分のモル比を決める
! * 比熱は乾燥気塊のもので代表させる
! * 流体の方程式において比熱は乾燥成分のもので代表させているため
! * 大気の平均分子量には湿潤成分の分子量を効果
! * 流体の方程式において, 湿潤成分の分子量は考慮しているため
!
!暗黙の型宣言禁止
implicit none
real(8), intent(in) :: a_MolFrIni(1:SpcNum) !下部境界でのモル比
real(8), intent(in) :: Humidity !相対湿度 ( Humidity <= 1.0 )
real(8), intent(out):: z_Temp(DimZMin:DimZMax) !温度
real(8), intent(out):: z_Press(DimZMin:DimZMax)!圧力
real(8), intent(out):: z_MolWtMean(DimZMin:DimZMax)
!平均分子量
real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum)
!モル分率
real(8) :: MolWtMeanDry !乾燥気塊の分子量の作業変数
real(8) :: MolWtMeanWet !湿潤気塊の分子量の作業変数
real(8) :: SatPress !飽和蒸気圧
real(8) :: VapPress !蒸気圧
real(8) :: DelMolFr
integer :: k, s
!-------------------------------------------------------------
! 配列の初期化
!-------------------------------------------------------------
!地表面の分子量を決める
za_MolFr(RegZMin+1, 1:SpcNum) = a_MolFrIni(1:SpcNum)
!地表面での平均分子量を決める
MolWtMeanDry = MolWtDry * (1.0d0 - sum(a_MolFrIni))
MolWtMeanWet = dot_product(MolWtWet, a_MolFrIni)
z_MolWtMean(RegZMin+1) = MolWtMeanDry + MolWtMeanWet
!地表面での温度(RegZMin+1 は, 高度 DelZ / 2 に相当)
z_Temp = 1.0d-60
z_Temp(RegZMin+1) = TempSfc - Grav * z_MolWtMean(RegZMin+1) / CpDryMol * ( DelZ * 5.0d-1 )
!地表面での圧力(RegZMin は, 高度 DelZ / 2 に相当)
z_Press = 1.0d-60
z_Press(RegZMin+1) = PressSfc *((TempSfc / z_Temp(RegZMin+1)) ** (- CpDryMol / GasRUniv))
!-----------------------------------------------------------
! (1) 乾燥断熱線に沿った温度を決める
! (2) 静水圧平衡から圧力を求める
! (3) (1),(2) の温度圧力に対して, とある相対湿度となるモル比を決める
!-----------------------------------------------------------
DtDz: do k = RegZMin+1, DimZMax-1
!(1)乾燥断熱線に沿って k+1 での温度を計算
z_Temp(k+1) = z_Temp(k) - Grav * z_MolWtMean(k) / CpDryMol * DelZ
!念為
if (z_Temp(k+1) <= 0.0d0 ) z_Temp(k+1) = z_Temp(k)
!(2)圧力を静水圧平衡から計算
z_Press(k+1) = z_Press(k) * ((z_Temp(k) / z_Temp(k+1)) ** (- CpDryMol / GasRUniv))
!(3)モル比の計算
! まずはモル比は変化しないものとしてモル比を与える
! 飽和蒸気圧と平衡定数との平衡条件の前に適用しておく
za_MolFr(k+1,:) = za_MolFr(k,:)
do s = 1, CondNum
!飽和蒸気圧
SatPress = SvapPress( SpcWetID(IdxCC(s)), z_Temp(k+1) )
!元々のモル分率を用いて現在の蒸気圧を計算
VapPress = za_MolFr(k,IdxCG(s)) * z_Press(k+1)
!飽和蒸気圧と圧力から現在のモル比を計算
if ( VapPress > SatPress ) then
za_MolFr(k+1,IdxCG(s)) = max(SatPress * Humidity / z_Press(k+1), 1.0d-16)
end if
end do
!NH4SH の平衡条件
if ( IdxNH3 /= 0 ) then
DelMolFr = max ( DelMolFrNH4SH( z_Temp(k+1), z_Press(k+1), za_MolFr(k+1,IdxNH3), za_MolFr(k+1,IdxH2S), Humidity ), 0.0d0 )
za_MolFr(k+1,IdxNH3) = za_MolFr(k+1,IdxNH3) - DelMolFr
za_MolFr(k+1,IdxH2S) = za_MolFr(k+1,IdxH2S) - DelMolFr
end if
!------------------------------------------------------------
!温度勾配を計算
!------------------------------------------------------------
!地表面での平均分子量を決める
MolWtMeanDry = MolWtDry * (1.0d0 - sum(za_MolFr(k+1, 1:SpcNum)))
MolWtMeanWet = dot_product(MolWtWet(1:SpcNum), za_MolFr(k+1, 1:SpcNum))
z_MolWtMean(k+1) = MolWtMeanDry + MolWtMeanWet
end do DtDz
end subroutine ECCM_Dry
| Subroutine : | |
| a_MolFrIni(1:SpcNum) : | real(8), intent(in) |
| Humidity : | real(8), intent(in) |
| z_Temp(DimZMin:DimZMax) : | real(8), intent(in) |
| z_Press(DimZMin:DimZMax) : | real(8), intent(in) |
| za_MolFr(DimZMin:DimZMax, 1:SpcNum) : | real(8), intent(out) |
与えられた温度に対し, 気塊が断熱的に上昇した時に実現される モル比のプロファイルを求める
subroutine ECCM_MolFr( a_MolFrIni, Humidity, z_Temp, z_Press, za_MolFr )
!
! 与えられた温度に対し, 気塊が断熱的に上昇した時に実現される
! モル比のプロファイルを求める
!
!暗黙の型宣言禁止
implicit none
real(8), intent(in) :: a_MolFrIni(1:SpcNum)
real(8), intent(in) :: Humidity
real(8), intent(in) :: z_Temp(DimZMin:DimZMax)
real(8), intent(in) :: z_Press(DimZMin:DimZMax)
real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum)
real(8) :: DelMolFr
integer :: k, s
!-----------------------------------------------------------
! 配列の初期化
!-----------------------------------------------------------
do s = 1, SpcNum
za_MolFr(:,s) = a_MolFrIni(s)
end do
!-----------------------------------------------------------
! 断熱減率 dT/dz の計算.
!-----------------------------------------------------------
do k = RegZMin, DimZMax
za_MolFr(k,:) = za_MolFr(k-1,:)
!------------------------------------------------------------
!NH4SH 以外の化学種の平衡条件
!------------------------------------------------------------
do s = 1, CondNum
!モル比を求める
!モル比は前のステップでのモル比を超えることはない
za_MolFr(k,IdxCG(s)) = min( za_MolFr(k-1,IdxCG(s)), SvapPress( SpcWetID(IdxCC(s)), z_Temp(k) ) * Humidity / z_Press(k) )
end do
!------------------------------------------------------------
!NH4SH の平衡条件
!------------------------------------------------------------
if ( IdxNH3 /= 0 ) then
!モル比の変化.
!とりあえず NH4SH に対する飽和比は 1.0 とする(手抜き...).
DelMolFr = max ( DelMolFrNH4SH( z_Temp(k), z_Press(k), za_MolFr(k,IdxNH3), za_MolFr(k,IdxH2S), Humidity ), 0.0d0 )
za_MolFr(k,IdxNH3) = za_MolFr(k,IdxNH3) - DelMolFr
za_MolFr(k,IdxH2S) = za_MolFr(k,IdxH2S) - DelMolFr
end if
end do
end subroutine ECCM_MolFr
| Subroutine : | |||||
| xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax) : | real(8), intent(in) | ||||
| xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax) : | real(8), intent(in) | ||||
| xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum) : | real(8), intent(in)
& xz_Stab, xz_StabTemp, xz_StabMolWt ) |
subroutine ECCM_Stab( xz_PotTemp, xz_Exner, xza_MixRt )
! & xz_Stab, xz_StabTemp, xz_StabMolWt )
use gridset, only: DimXMin, DimXMax, DimZMin, DimZMax, SpcNum !
use basicset,only: MolWtDry, MolWtWet, CpDry, Grav, xz_ExnerBasicZ, xz_PotTempBasicZ, xz_EffMolWtBasicZ, xza_MixRtBasicZ
use average, only: xz_avr_xr
use differentiate_center2, only: xr_dz_xz
implicit none
real(8), intent(in) :: xz_PotTemp(DimXMin:DimXMax,DimZMin:DimZMax)
real(8), intent(in) :: xz_Exner(DimXMin:DimXMax, DimZMin:DimZMax)
real(8), intent(in) :: xza_MixRt(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
! real(8), intent(out) :: xz_Stab(DimXMin:DimXMax, DimZMin:DimZMax)
! real(8), intent(out) :: xz_StabTemp(DimXMin:DimXMax, DimZMin:DimZMax)
! real(8), intent(out) :: xz_StabMolWt(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: xz_Stab(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: xz_StabTemp(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: xz_StabMolWt(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: xza_MolFrAll(DimXMin:DimXMax,DimZMin:DimZMax,SpcNum)
real(8) :: xz_TempAll(DimXMin:DimXMax, DimZMin:DimZMax)
real(8) :: xz_MolWtWet(DimXMin:DimXMax, DimZMin:DimZMax)
integer :: i, k, s
xz_TempAll = (xz_PotTemp + xz_PotTempBasicZ) * (xz_Exner + xz_ExnerBasicZ)
do s = 1, SpcNum
xza_MolFrAll(:,:,s) = (xza_MixRt(:,:,s) + xza_MixRtBasicZ(:,:,s)) * MolWtDry / MolWtWet(s)
end do
do k = DimZMin, DimZMax
do i = DimXMin, DimXMax
xz_MolWtWet(i,k) = dot_product( MolWtWet(1:GasNum), xza_MolFrAll(i,k,1:GasNum) )
end do
end do
xz_StabTemp = Grav / xz_TempAll * ( xz_avr_xr( xr_dz_xz( xz_TempAll ) ) + Grav * xz_EffMolWtBasicZ / CpDry )
xz_StabMolWt = - Grav * xz_avr_xr( xr_dz_xz( xz_MolWtWet ) ) / ( MolWtDry * xz_EffMolWtBasicZ )
xz_Stab = xz_StabTemp + xz_StabMolWt
! xz_Stab = &
! & Grav / xz_TempAll &
! & * ( xz_avr_xr( xr_dz_xz( xz_TempAll ) ) &
! & + Grav * xz_EffMolWtBasicZ / CpDry ) &
! & - Grav * xz_avr_xr( xr_dz_xz( xz_MolWtWet ) ) &
! & / ( MolWtDry * xz_EffMolWtBasicZ )
call StoreStabTemp( xz_StabTemp )
call StoreStabMolWt( xz_StabMolWt )
where (xz_Stab < 1.0d-7)
xz_Stab = 1.0d-7
end where
end subroutine ECCM_Stab
| Subroutine : | |||
| a_MolFrIni(1:SpcNum) : | real(8), intent(in)
| ||
| Humidity : | real(8), intent(in)
| ||
| z_Temp(DimZMin:DimZMax) : | real(8), intent(out)
| ||
| z_Press(DimZMin:DimZMax) : | real(8), intent(out)
| ||
| z_MolWtMean(DimZMin:DimZMax) : | real(8), intent(out)
| ||
| za_MolFr(DimZMin:DimZMax, 1:SpcNum) : | real(8), intent(out)
|
* 湿潤断熱温度減率に沿った温度、圧力を求め、それらに対して指定された相対湿度となるように凝結成分のモル比を決める * 比熱は乾燥気塊のもので代表させる * 流体の方程式において比熱は乾燥成分のもので代表させているため * 大気の平均分子量には湿潤成分の分子量を効果 * 流体の方程式において, 湿潤成分の分子量は考慮しているため
subroutine ECCM_Wet( a_MolFrIni, Humidity, z_Temp, z_Press, z_MolWtMean, za_MolFr )
!
!== 概要
! * 湿潤断熱温度減率に沿った温度、圧力を求め、それらに対して指定された相対湿度となるように凝結成分のモル比を決める
! * 比熱は乾燥気塊のもので代表させる
! * 流体の方程式において比熱は乾燥成分のもので代表させているため
! * 大気の平均分子量には湿潤成分の分子量を効果
! * 流体の方程式において, 湿潤成分の分子量は考慮しているため
!
!暗黙の型宣言禁止
implicit none
real(8), intent(in) :: a_MolFrIni(1:SpcNum) !下部境界でのモル比
real(8), intent(in) :: Humidity !相対湿度 ( Humidity <= 1.0 )
real(8), intent(out):: z_Temp(DimZMin:DimZMax) !温度
real(8), intent(out):: z_Press(DimZMin:DimZMax)!圧力
real(8), intent(out):: z_MolWtMean(DimZMin:DimZMax)
!平均分子量
real(8), intent(out):: za_MolFr(DimZMin:DimZMax, 1:SpcNum)
!モル分率
real(8) :: MolWtMeanDry !乾燥気塊の分子量の作業変数
real(8) :: MolWtMeanWet !湿潤気塊の分子量の作業変数
real(8) :: SatPress !飽和蒸気圧
real(8) :: VapPress !蒸気圧
real(8) :: DelMolFr
real(8) :: a_MolFr(SpcNum) !モル比の作業配列
integer :: k, s
real(8) :: Temp1, Press1, DTempDZ1
real(8) :: Temp2, Press2, DTempDZ2
real(8) :: Temp3, Press3, DTempDZ3
real(8) :: Temp4, Press4, DTempDZ4
real(8) :: DTempDZ
!-------------------------------------------------------------
! 配列の初期化
!-------------------------------------------------------------
!地表面の分子量を決める
za_MolFr(RegZMin+1, 1:SpcNum) = a_MolFrIni(1:SpcNum)
!地表面での平均分子量を決める
MolWtMeanDry = MolWtDry * (1.0d0 - sum(a_MolFrIni))
MolWtMeanWet = dot_product(MolWtWet, a_MolFrIni)
z_MolWtMean(RegZMin+1) = MolWtMeanDry + MolWtMeanWet
!地表面での温度(RegZMin は, 高度 DelZ / 2 に相当)
z_Temp = 1.0d-60
z_Temp(RegZMin+1) = TempSfc - Grav * z_MolWtMean(RegZMin+1) / CpDryMol * ( DelZ * 5.0d-1 )
!地表面での圧力(RegZMin は, 高度 DelZ / 2 に相当)
z_Press = 1.0d-60
z_Press(RegZMin+1) = PressSfc *((TempSfc / z_Temp(RegZMin+1)) ** (- CpDryMol / GasRUniv))
!-----------------------------------------------------------
! 断熱減率 dT/dz の計算.
!-----------------------------------------------------------
DtDz: do k = RegZMin+1, DimZMax-1
!初期化
za_MolFr(k+1,:) = za_MolFr(k,:)
!----------------------------------------------------
!湿潤断熱減率をルンゲクッタ法を用いて計算
!----------------------------------------------------
! (0) 高度 k での値を作業配列に保管
Temp1 = z_Temp(k)
Press1 = z_Press(k)
a_MolFr = za_MolFr(k,:)
! (1) 高度 k での値を用いて温度変化を計算
call ECCM_DTempDZ( Temp1, Press1, a_MolFr, DTempDZ1 )
! (2) (1) で求めた値を用いて, 高度 k + Δk/2 での値を用いて温度変化を計算
! このとき, 分子量は変化しないものとする.
Temp2 = Temp1 + DTempDZ1 * DelZ * 5.0d-1
Press2 = Press1 * ((Temp1 / Temp3) ** (Grav * MolWtDry / (GasRUniv * DTempDZ2)))
call ECCM_DTempDZ( Temp2, Press2, a_MolFr, DTempDZ2 )
! (3) (2) で求めた値を用いて, 高度 k + Δk/2 での値を用いて温度変化を計算
! このとき, 分子量は変化しないものとする.
Temp3 = Temp1 + DTempDZ2 * DelZ * 5.0d-1
Press3 = Press1 * ((Temp1 / Temp3) ** (Grav * MolWtDry / (GasRUniv * DTempDZ2)))
call ECCM_DTempDZ( Temp3, Press3, a_MolFr, DTempDZ3 )
! (4) (3) で求めた値を用いて, 高度 k + Δk での値を用いて温度変化を計算
! このとき, 分子量は変化しないものとする.
Temp4 = Temp1 + DTempDZ3 * DelZ
Press4 = Press1 * ((Temp1 / Temp4) ** (Grav * MolWtDry / (GasRUniv * DTempDZ3)))
call ECCM_DTempDZ( Temp4, Press4, a_MolFr, DTempDZ4 )
! (5) 最終的な傾きを求める
DTempDZ = (DTempDZ1 + DTempDZ2 * 2.0d0 + DTempDZ3 * 2.0d0 + DTempDZ4) / 6.0d0
!----------------------------------------------------
!得られた温度減率より温度と圧力を決める
!----------------------------------------------------
!温度を計算
z_Temp(k+1) = z_Temp(k) + DTempDz * DelZ
!為念
if(z_Temp(k+1) < 0.0d0) z_Temp(k+1) = z_Temp(k)
!圧力を静水圧平衡から計算
z_Press(k+1) = z_Press(k) * ( ( z_Temp(k) / z_Temp(k+1)) ** (Grav * MolWtDry / ( DTempDZ * GasRUniv ) ) )
!----------------------------------------------------
!モル比の計算
!----------------------------------------------------
do s = 1, CondNum
!飽和蒸気圧
SatPress = SvapPress( SpcWetID(IdxCC(s)), z_Temp(k+1) )
!元々のモル分率を用いて現在の蒸気圧を計算
VapPress = za_MolFr(k,IdxCG(s)) * z_Press(k+1)
!飽和蒸気圧と圧力から現在のモル比を計算
if ( VapPress > SatPress ) then
za_MolFr(k+1,IdxCG(s)) = max(SatPress * Humidity / z_Press(k+1), 1.0d-16)
end if
end do
!NH4SH の平衡条件
if ( IdxNH3 /= 0 ) then
DelMolFr = max ( DelMolFrNH4SH( z_Temp(k+1), z_Press(k+1), za_MolFr(k+1,IdxNH3), za_MolFr(k+1,IdxH2S), Humidity ), 0.0d0 )
za_MolFr(k+1,IdxNH3) = za_MolFr(k+1,IdxNH3) - DelMolFr
za_MolFr(k+1,IdxH2S) = za_MolFr(k+1,IdxH2S) - DelMolFr
end if
!------------------------------------------------------------
!地表面での平均分子量を決める
!------------------------------------------------------------
MolWtMeanDry = MolWtDry * (1.0d0 - sum(za_MolFr(k+1, 1:SpcNum)))
MolWtMeanWet = dot_product(MolWtWet(1:SpcNum), za_MolFr(k+1, 1:SpcNum))
z_MolWtMean(k+1) = MolWtMeanDry + MolWtMeanWet
end do DtDz
end subroutine ECCM_Wet