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:SpcNum) = a_MolFrIni(1:SpcNum) !地表面での平均分子量を決める MolWtMeanDry = MolWtDry * (1.0d0 - sum(a_MolFrIni)) MolWtMeanWet = dot_product(MolWtWet, a_MolFrIni) z_MolWtMean(RegZMin) = MolWtMeanDry + MolWtMeanWet !地表面での温度(RegZMin は, 高度 DelZ / 2 に相当) z_Temp = 1.0d-60 z_Temp(RegZMin) = TempSfc - Grav * z_MolWtMean(RegZMin) / CpDryMol * ( DelZ * 5.0d-1 ) !地表面での圧力(RegZMin は, 高度 DelZ / 2 に相当) z_Press = 1.0d-60 z_Press(RegZMin) = PressSfc *((TempSfc / z_Temp(RegZMin)) ** (- CpDryMol / GasRUniv)) !----------------------------------------------------------- ! (1) 乾燥断熱線に沿った温度を決める ! (2) 静水圧平衡から圧力を求める ! (3) (1),(2) の温度圧力に対して, とある相対湿度となるモル比を決める !----------------------------------------------------------- DtDz: do k = RegZMin, 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:SpcNum) = a_MolFrIni(1:SpcNum) !地表面での平均分子量を決める MolWtMeanDry = MolWtDry * (1.0d0 - sum(a_MolFrIni)) MolWtMeanWet = dot_product(MolWtWet, a_MolFrIni) z_MolWtMean(RegZMin) = MolWtMeanDry + MolWtMeanWet !地表面での温度(RegZMin は, 高度 DelZ / 2 に相当) z_Temp = 1.0d-60 z_Temp(RegZMin) = TempSfc - Grav * z_MolWtMean(RegZMin) / CpDryMol * ( DelZ * 5.0d-1 ) !地表面での圧力(RegZMin は, 高度 DelZ / 2 に相当) z_Press = 1.0d-60 z_Press(RegZMin) = PressSfc *((TempSfc / z_Temp(RegZMin)) ** (- CpDryMol / GasRUniv)) !----------------------------------------------------------- ! 断熱減率 dT/dz の計算. !----------------------------------------------------------- DtDz: do k = RegZMin, 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