!--------------------------------------------------------------------- ! Copyright (C) GFD Dennou Club, 2005. All rights reserved. !--------------------------------------------------------------------- != Subroutine ECCM ! ! * Developer: SUGIYAMA Ko-ichiro ! * Version: $Id: eccm.f90,v 1.1.2.1 2006/01/13 06:10:20 sugiyama Exp $ ! * Tag Name: $Name: $ ! * Change History: ! !== Overview ! !断熱的に上昇する気塊の温度減率を計算し, 静水圧平衡から圧力を決める ! !== Error Handling ! !== Known Bugs ! !== Note ! ! * 比熱は定数, 平均分子量は配列 ! * オイラースキームだと精度が足りないので, ルンゲクッタスキームを利用. ! !== Future Plans ! module ECCM !モジュール読み込み use gridset, only: DimZMin, &! 配列の Z 方向の下限 & DimZMax, &! 配列の Z 方向の上限 & RegZMin, &! & SpcNum, &! & DelZ ! use basicset,only: MolWtDry, &! & MolWtWet, &! & CpDryMol, &! & SpcWetID, &! & SpcWetSymbol, &! & TempSfc, &! & PressSfc, &! & Grav ! use chemcalc, only: SvapPress, &! & LatentHeat ! use ChemData, only: GasRUniv, &! & ChemData_OneSpcID !暗黙の型宣言禁止 implicit none !属性の指定 private !関数の公開 public ECCM_Init public ECCM_Temp_Press public ECCM_MolFr !変数定義 integer :: LoopNum = 0 integer :: GasNum(10) = 0 integer :: CloudNum(10) = 0 integer :: NH3Num = 0 integer :: H2SNum = 0 save LoopNum, GasNum, CloudNum save NH3Num, H2SNum contains subroutine ECCM_Init( ) !暗黙の型宣言禁止 implicit none !変数定義 integer :: s integer :: n1 !----------------------------------------------------------- ! 雲粒と気体の ID の組を作る !----------------------------------------------------------- !初期化 LoopNum = 0 !化学種の中から雲粒を作るものを選び, その配列添え字と分子量を保管. SelectCloud: do s = 1, SpcNum ! NH4SH については無視 if ( trim(SpcWetSymbol(s)) == 'NH4SH-s-Cloud' ) then cycle SelectCloud end if !'Cloud' という文字列が含まれるものの個数を数える n1 = index(SpcWetSymbol(s), '-Cloud' ) if (n1 /= 0) then LoopNum = LoopNum + 1 CloudNum(LoopNum)= s GasNum(LoopNum) = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID(SpcWetSymbol(s)(1:n1-3) // '-g')) end if end do SelectCloud !----------------------------------------------------------- ! 硫化アンモニウム, およびアンモニアと硫化水素の ID を取得 !----------------------------------------------------------- NH3Num = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('NH3-g')) H2SNum = minloc(SpcWetID, 1, SpcWetID == ChemData_OneSpcID('H2S-g')) !----------------------------------------------------------- ! 確認 !----------------------------------------------------------- if ( LoopNum == 0 ) then write(*,*) "MoistAdjust: CloudNum = 0, please comment out of MoistAdjust" stop end if write(*,*) "ECCM_Init, LoopNum: ", LoopNum write(*,*) "ECCM_Init, CloudNum: ", CloudNum write(*,*) "ECCM_Init, GasNum: ", GasNum write(*,*) "ECCM_Init, NH3Num: ", NH3Num write(*,*) "ECCM_Init, H2SNum: ", H2SNum end subroutine ECCM_Init !!!------------------------------------------------------------------------------!!! subroutine ECCM_Temp_Press( MolFrIni, z_Temp, z_Press ) !暗黙の型宣言禁止 implicit none real(8), intent(in) :: MolFrIni(1:SpcNum) real(8), intent(out):: z_Temp(DimZMin:DimZMax) real(8), intent(out):: z_Press(DimZMin:DimZMax) real(8) :: z_MolFr(DimZMin:DimZMax, 0:SpcNum) !モル分率 real(8) :: z_DTempDZ(DimZMin:DimZMax) real(8) :: z_DTempDZDry(DimZMin:DimZMax) real(8) :: MolWt( 0:SpcNum) real(8) :: z_MolWtMean(DimZMin:DimZMax) real(8) :: z_dM(DimZMin:DimZMax) real(8) :: z_dp(DimZMin:DimZMax) real(8) :: z_dMdz(DimZMin:DimZMax) real(8) :: z_Stab(DimZMin:DimZMax) real(8) :: ReactHeat real(8) :: Heat(SpcNum) real(8) :: SatPress real(8) :: VapPress !分圧 real(8) :: A, B real(8) :: Solution(2) real(8) :: PPress(2) real(8) :: Kp real(8) :: DelPress integer :: k, s !------------------------------------------------------------- ! 配列の初期化 !------------------------------------------------------------- !初期モル比 z_MolFr = 0.0d0 z_MolFr(RegZMin, 1:SpcNum) = MolFrIni(1:SpcNum) z_MolFr(RegZMin, 0) = 1.0d0 - sum(MolFrIni) z_MolFr(RegZMin-1, 1:SpcNum) = MolFrIni(1:SpcNum) z_MolFr(RegZMin-1, 0) = 1.0d0 - sum(MolFrIni) !地表面での温度(RegZMin は, 高度 DelZ / 2 に相当) z_Temp = 0.0d0 z_Temp(RegZMin) = TempSfc - Grav * MolWtDry / CpDryMol * DelZ * 5.0d-1 !地表面での圧力(RegZMin は, 高度 DelZ / 2 に相当) z_Press = 0.0d0 z_Press(RegZMin) = PressSfc *((TempSfc / z_Temp(RegZMin)) & & ** ( - CpDryMol / GasRUniv )) !分子量の初期化 MolWt(0) = MolWtDry MolWt(1:SpcNum) = MolWtWet(1:SpcNum) z_MolWtMean = 0.0d0 !----------------------------------------------------------- ! 断熱減率 dT/dz の計算. !----------------------------------------------------------- do k = RegZMin, DimZMax-1 !初期化 Heat = 0.0d0 ReactHeat = 0.0d0 z_MolFr(k,:) = z_MolFr(k-1,:) !------------------------------------------------------------ !NH4SH 以外の化学種の平衡条件 !------------------------------------------------------------ do s = 1, LoopNum !飽和蒸気圧 SatPress = SvapPress( SpcWetID(CloudNum(s)), z_Temp(k) ) !潜熱 Heat(GasNum(s)) = LatentHeat( SpcWetID(CloudNum(s)), z_Temp(k) ) !前のステップでのモル分率を用いて現在の蒸気圧を計算 VapPress = z_MolFr(k-1, GasNum(s)) * z_Press(k) !飽和蒸気圧から凝結の有無を決める if ( VapPress < SatPress ) then !モル比 z_MolFr(k,GasNum(s)) = z_MolFr(k-1,GasNum(s)) !凝結していないので潜熱なし. Heat(GasNum(s)) = 0.0d0 else !飽和蒸気圧と圧力から現在のモル比を計算 z_MolFr(k,GasNum(s)) = max(SatPress / z_Press(k), 1.0d-16) end if write(*,*) 'ECCM, Heat: ', Heat(GasNum(s)) end do !------------------------------------------------------------ !NH4SH の平衡条件 !------------------------------------------------------------ !アンモニアと硫化水素の分圧 PPress(1) = z_MolFr(k, NH3Num) * z_Press(k) PPress(2) = z_MolFr(k, H2SNum) * z_Press(k) !圧平衡定数 Kp = 61.781d0 - 10834.0d0 / z_Temp(k) - dlog(1.0d2) !気圧変化を二次方程式の解として求める. Solution(1) = 5.0d-1 * (sum(PPress) & & + dsqrt( (PPress(1) - PPress(2))**2.0d0 & & + 4.0d0 * dexp(Kp)) ) Solution(2) = 5.0d-1 * (sum(PPress) & & - dsqrt( (PPress(1) - PPress(2))**2.0d0 & & + 4.0d0 * dexp(Kp)) ) !成立条件のチェック if ( minval(Solution, 1, Solution > 0.0d0) > minval(PPress, 1) ) then DelPress = 0.0d0 ! z_MolFr(k,NH3Num) = z_MolFr(k,NH3Num) ! z_MolFr(k,H2SNum) = z_MolFr(k,H2SNum) else DelPress = minval(Solution, 1, Solution > 0.0d0) ! z_MolFr(k,NH3Num) = z_MolFr(k,NH3Num) - DelPress / z_Press(k) ! z_MolFr(k,H2SNum) = z_MolFr(k,H2SNum) - DelPress / z_Press(k) end if !成立条件のチェック ! 変化量はゼロ以上でなければならない ! DelPress = minval( Solution, Solution > 0.0d0 ) ! DelPress = min( DelPress, minval(PPress) ) !反応熱解放 [K / kg] ! ReactHeat = 10834.0d0 * GasRUniv * DelPress / CpDryMol !------------------------------------------------------------ !次のステップの温度, 圧力を求める !------------------------------------------------------------ !モル比 z_MolFr(k,0) = 1.0d0 - sum(z_MolFr(k,1:SpcNum)) !平均分子量 z_MolWtMean(k) = dot_product(MolWt, z_MolFr(k, :)) !係数. 温度 Temp(i) で評価 A = dot_product( Heat(1:SpcNum), z_MolFr(k,1:SpcNum)) & & / ( GasRUniv * z_Temp(k) ) B = dot_product(( Heat(1:SpcNum) ** 2.0d0), z_MolFr(k,1:SpcNum)) & & / ( CpDryMol * GasRUniv * ( z_Temp(k) ** 2.0d0 ) ) !断熱温度減率 z_DTempDZ(k) = - Grav * z_MolWtMean(k) * ( 1.0d0 + A ) & & / ( CpDryMol * ( 1.0d0 + B ) ) & & + ReactHeat / DelZ !温度を計算 z_Temp(k+1) = z_Temp(k) + z_DTempDz(k) * DelZ !圧力を静水圧平衡から計算 z_Press(k+1) = z_Press(k) * ( ( z_Temp(k) / z_Temp(k+1)) & & ** (Grav * z_MolWtMean(k) / (GasRUniv * z_DTempDZ(k)) ) ) write(*,*) 'ECCM, MolFr: ', z_MolFr(k,:) write(*,*) 'ECCM, Temp: ', z_Temp(k) write(*,*) 'ECCM, Press: ', z_Press(k) write(*,*) 'ECCM, DTempDZ: ', z_DTempDZ(k) !モル比 (RegZMax でも値を持つようにするため) z_MolFr(k+1,:) = z_MolFr(k,:) z_MolFr(k+1,0) = 1.0d0 - sum( z_MolFr(k,1:SpcNum) ) end do !----------------------------------------------------------- ! 安定度の計算 !----------------------------------------------------------- !分子量, 圧力差 do k = RegZMin+1, DimZMax z_dM(k) = ( z_MolWtMean(k+1) - z_MolWtMean(k-1) ) * 5.0d-1 z_dp(k) = ( z_Press(k+1) - z_Press(k-1) ) * 5.0d-1 end do !乾燥断熱温度減率 do k = RegZMin+1, DimZMax z_DTempDZDry(k) = - Grav * z_MolWtMean(k) / CpDryMol end do do k = RegZMin, DimZMax z_dMdz(k) = - z_dM(k) * z_Press(k) * Grav * z_MolWtMean(k) & & / ( z_dp(k) * GasRUniv * z_Temp(k) ) z_Stab(k) = Grav * ( z_DTempDz(k) - z_DTempDzDry(k) ) / z_Temp(k) & & - Grav * z_dMdz(k) / z_MolWtMean(k) end do where (z_Stab < 1.0d-7) z_Stab = 1.0d-7 end where end subroutine ECCM_Temp_Press !!!------------------------------------------------------------------------------!!! subroutine ECCM_MolFr( MolFrIni, Humidity, z_Temp, z_Press, z_MolFr ) !暗黙の型宣言禁止 implicit none real(8), intent(in):: 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):: z_MolFr(DimZMin:DimZMax, 1:SpcNum) real(8) :: SatPress real(8) :: Solution(2) real(8) :: PPress(2) real(8) :: Kp real(8) :: DelPress integer :: k, s !----------------------------------------------------------- ! 配列の初期化 !----------------------------------------------------------- !初期モル比 z_MolFr = 0.0d0 z_MolFr(RegZMin, 1:SpcNum) = MolFrIni(1:SpcNum) z_MolFr(RegZMin-1, 1:SpcNum) = MolFrIni(1:SpcNum) !----------------------------------------------------------- ! 断熱減率 dT/dz の計算. !----------------------------------------------------------- do k = RegZMin, DimZMax z_MolFr(k,:) = z_MolFr(k-1,:) !------------------------------------------------------------ !NH4SH 以外の化学種の平衡条件 !------------------------------------------------------------ do s = 1, LoopNum !飽和蒸気圧 SatPress = SvapPress( SpcWetID(CloudNum(s)), z_Temp(k) ) !飽和蒸気圧と圧力からモル比を計算 ! 前のステップでの値よりも大きくならないようにする z_MolFr(k,GasNum(s)) = & & min( SatPress * Humidity / z_Press(k), z_MolFr(k-1,GasNum(s)) ) end do !------------------------------------------------------------ !NH4SH の平衡条件 !------------------------------------------------------------ !アンモニアと硫化水素の分圧 PPress(1) = z_MolFr(k, NH3Num) * z_Press(k) PPress(2) = z_MolFr(k, H2SNum) * z_Press(k) !圧平衡定数 Kp = 61.781d0 - 10834.0d0 / z_Temp(k) - dlog(1.0d2) !気圧変化を二次方程式の解として求める. Solution(1) = 5.0d-1 * (sum(PPress) & & + dsqrt( (PPress(1) - PPress(2))**2.0d0 & & + 4.0d0 * dexp(Kp)) ) Solution(2) = 5.0d-1 * (sum(PPress) & & - dsqrt( (PPress(1) - PPress(2))**2.0d0 & & + 4.0d0 * dexp(Kp)) ) !成立条件のチェック ! NH4SH の昇華を考えないので, 変化量は正でなければならない. if ( minval(Solution, 1, Solution > 0.0d0) > minval(PPress, 1) ) then DelPress = 0.0d0 z_MolFr(k,NH3Num) = z_MolFr(k,NH3Num) z_MolFr(k,H2SNum) = z_MolFr(k,H2SNum) else DelPress = minval(Solution, 1, Solution > 0.0d0) z_MolFr(k,NH3Num) = z_MolFr(k,NH3Num) - DelPress * Humidity / z_Press(k) z_MolFr(k,H2SNum) = z_MolFr(k,H2SNum) - DelPress * Humidity / z_Press(k) end if end do end subroutine ECCM_MolFr end module ECCM