!--------------------------------------------------------------------- ! Copyright (C) GFD Dennou Club, 2005. All rights reserved. !--------------------------------------------------------------------- != Subroutine ECCM2 ! ! * 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 ! ! * 比熱は定数, 平均分子量は配列 ! * 積分は 1 次精度であるが, 十分にステップ数を細かくすれば OK. ! * 水のみの場合でルンゲクッタを用いた場合と値は同じ. ! !== Future Plans ! subroutine eccm2( MolFrIni) use gridset, only: RegZMin, RegZMax, SpcNum, DelZ, f_Z use dataset, only: MolWtDry, MolWtWet, & & CpDryMol, CpWetMol,SpcWetID, & & TempSfc, PressSfc, Grav use chemdata, only: GasRUniv use chemcalc, only: SvapPress, LatentHeat, ChemPotOneSpc, ChemPotNH4SH, & & EntropyOneSpc, EntropyNH4SH implicit none real(8), intent(in) :: MolFrIni(1:SpcNum) real(8) :: Temp( RegZMin:RegZMax+1 ) real(8) :: Press( RegZMin:RegZMax+1 ) real(8) :: DTempDz( RegZMin:RegZMax ) real(8) :: DTempDzDry( RegZMin:RegZMax ) real(8) :: MolFr( RegZMin-1:RegZMax+1, 0:SpcNum ) !モル分率 real(8) :: Cp(0:SpcNum) real(8) :: CpMean( RegZMin-1:RegZMax+1) real(8) :: MolWt( 0:SpcNum) real(8) :: MolWtMean( RegZMin-1:RegZMax+1) real(8) :: dM( RegZMin:RegZMax) real(8) :: dp( RegZMin:RegZMax) real(8) :: dMdz( RegZMin:RegZMax) real(8) :: Stab( RegZMin:RegZMax) real(8) :: LatentHeatAll(SpcNum ) real(8) :: SvapPressAll(SpcNum) !飽和蒸気圧 real(8) :: VapPress !分圧 real(8) :: A, B real(8) :: VapPressNH3 real(8) :: VapPressH2S real(8) :: MolFrNH3 real(8) :: MolFrH2S real(8) :: num real(8) :: del, Q integer :: i, s integer, parameter :: ID_NH3 = 7 integer, parameter :: ID_H2S = 9 integer, parameter :: ID_NH4SH = 9 integer, parameter :: Num_NH3 = 2 integer, parameter :: Num_H2S = 3 ! write(*,*) 'Sub ECCM2, MolFrIni: ', MolFrIni !!! !!!配列の初期化 !!! !初期モル比 MolFr = 0.0d0 MolFr(RegZMin, 1:SpcNum) = MolFrIni(1:SpcNum) MolFr(RegZMin, 0) = 1.0d0 - sum(MolFrIni) MolFr(RegZMin-1, 1:SpcNum) = MolFrIni(1:SpcNum) MolFr(RegZMin-1, 0) = 1.0d0 - sum(MolFrIni) !地表面での温度・圧力 Press = 0.0d0 Temp = 0.0d0 Press(RegZMin) = PressSfc Temp(RegZMin) = TempSfc Cp(0) = CpDryMol Cp(1:SpcNum) = CpWetMol ! CpMean = 0.0d0 MolWt(0) = MolWtDry MolWt(1:SpcNum) = MolWtWet MolWtMean = 0.0d0 num = 0.0d0 !!! !!!断熱減率 dT/dz の計算. !!! do i = RegZMin, RegZMax+1 write(*,*) i, Temp(i), Press(i) !初期化 Q = 0.0d0 !飽和蒸気圧 SvapPressAll(1:SpcNum) = SvapPress( Temp(i) ) !潜熱 LatentHeatAll(1:SpcNum) = LatentHeat( Temp(i) ) !------------------------------------------------------------ !NH4SH 以外の化学種の平衡条件 !------------------------------------------------------------ do s = 1, SpcNum if (s /= Num_H2S ) then !前のステップでのモル分率を用いて現在の蒸気圧を計算 VapPress = MolFr(i-1, s) * Press(i) !飽和蒸気圧から凝結の有無を決める if ( VapPress < SvapPressAll(s) ) then !モル比 MolFr(i,s) = MolFr(i-1,s) !凝結していないので潜熱なし. LatentHeatAll(s) = 0.0d0 else !飽和蒸気圧と圧力から現在のモル数を計算 MolFr(i,s) = max(SvapPressAll(s) / Press(i), 1.0d-16) end if end if end do !------------------------------------------------------------ !NH4SH の平衡条件 !------------------------------------------------------------ num = 0.0d0 Adjast: do del = 1.0d-9 * num MolFrNH3 = MolFr(i, Num_NH3) - del MolFrH2S = MolFr(i, Num_H2S) - del VapPressNH3 = MolFrNH3 * Press(i) VapPressH2S = MolFrH2S * Press(i) if (MolFrH2S <= 1.0d-16) then MolFr(i,Num_H2S) = max(MolFrH2S, 1.0d-16) LatentHeatAll(Num_H2S) = 0.0d0 exit Adjast end if if (MolFrNH3 <= 1.0d-16) then MolFr(i,Num_NH3) = max(MolFrNH3, 1.0d-16) LatentHeatAll(Num_NH3) = 0.0d0 exit Adjast end if if ( ChemPotNH4SH(ID_NH3, ID_H2S, Temp(i)) & & >= ChemPotOneSpc(ID_NH3, Temp(i), VapPressNH3) & & + ChemPotOneSpc(ID_H2S, Temp(i), VapPressH2S) ) then MolFr(i,Num_NH3) = max(MolFrNH3, 1.0d-16) MolFr(i,Num_H2S) = max(MolFrH2S, 1.0d-16) !反応熱解放 (これで OK だと思うが.... ) del = MolFr(i-1,Num_H2S) - MolFr(i,Num_H2S) Q = 10834.0d0 * GasRUniv * del / DelZ / CpDryMol LatentHeatAll(Num_NH3) = 0.0d0 LatentHeatAll(Num_H2S) = 0.0d0 if (num == 0.0d0) then LatentHeatAll(Num_NH3) = 0.0d0 end if exit Adjast end if num = num + 1.0d0 end do Adjast !------------------------------------------------------------ !次のステップの温度, 圧力を求める !------------------------------------------------------------ !モル比 MolFr(i,0) = 1.0d0 - sum(MolFr(i,1:SpcNum)) !平均分子量 MolWtMean(i) = dot_product(MolWt, MolFr(i,:)) !係数. 温度 Temp(i) で評価 A = dot_product(LatentHeatAll(1:SpcNum), MolFr(i,1:SpcNum)) & & / ( GasRUniv * Temp(i) ) B = dot_product((LatentHeatAll(1:SpcNum) ** 2.0d0), MolFr(i,1:SpcNum)) & & / ( CpDryMol * GasRUniv * ( Temp(i) ** 2.0d0 ) ) !断熱温度減率 DTempDz(i) = - Grav * MolWtMean(i) * ( 1.0d0 + A ) & & / ( CpDryMol * ( 1.0d0 + B ) ) & & + Q DTempDzDry(i) = - Grav * MolWtMean(i) / CpDryMol !温度を計算 Temp(i+1) = Temp(i) + DTempDz(i) * DelZ !圧力を静水圧平衡から計算 Press(i+1) = & & Press(i) - ( Press(i) * Grav * DelZ * MolWtMean(i)) & & / ( GasRUniv * Temp(i) ) !モル比 (RegZMax でも値を持つようにするため) MolFr(i+1,:) = MolFr(i,:) MolFr(i+1,0) = 1.0d0 - sum(MolFr(i,1:SpcNum)) end do !!! !!!安定度の計算 !!! !分子量, 圧力差 do i = RegZMin+1, RegZMax dM(i) = ( MolWtMean(i + 1) - MolWtMean(i - 1) ) * 5.0d-1 dp(i) = ( Press(i + 1) - Press(i - 1) ) * 5.0d-1 end do dM(RegZMin) = dM(RegZMin+1) dp(RegZMin) = dp(RegZMin+1) do i = RegZMin, RegZMax dMdz(i) = - dM(i) * Press(i) * Grav * MolWtMean(i) & & / ( dp(i) * GasRUniv * Temp(i) ) Stab(i) = Grav * ( DTempDz(i) - DTempDzDry(i) ) / Temp(i) & & - Grav * dMdz(i) / MolWtMean(i) write(*,*) "dM, dp ", dM(i), dp(i) write(*,*) "dMdz ", dMdz(i) write(*,*) "Deldtdzd ", DTempDz(i), DTempDzDry(i),DTempDz(i)-DTempDzDry(i) write(*,*) "Stab ", Stab(i) end do where (Stab < 1.0d-7) Stab = 1.0d-7 end where call HistoryFile() contains subroutine HistoryFile() use gt4_history use fileset, only: & & exptitle, & & expsrc, & & expinst implicit none real :: ID(0:SpcNum) ID(0) = 1.0 ID(1:SpcNum) = real(SpcWetID, 4) !----------------------------------------------------------- ! ヒストリー作成 !----------------------------------------------------------- call HistoryCreate( & & file = 'arare-eccm.nc', & & title = exptitle, & & source = expsrc, & & institution = expinst, & & dims=(/'p','s'/), & & dimsizes=(/RegZMax+1, SpcNum+1/), & & longnames=(/'Pressure ', & & 'Species Num '/), & & units=(/'Pa','1 '/), origin=0.0, & & interval=0.0 ) !----------------------------------------------------------- ! 変数出力 !----------------------------------------------------------- call HistoryPut('p', real(Press(RegZMin:RegZMax),4) ) call HistoryPut('s', ID) call HistoryAddAttr( 'p', 'positive', 'down' ) ! !圧力 ! call HistoryAddVariable( & ! & varname='Press', dims=(/'p'/), & ! & longname='Pressure', & ! & units='Pa', & ! & xtype='double' ) !温度 call HistoryAddVariable( & & varname='temp', dims=(/'p'/), & & longname='Temperature', & & units='K', & & xtype='double' ) !モル比 call HistoryAddVariable( & & varname='molfr', dims=(/'p','s'/), & & longname='Mole Fraction', & & units='1', & & xtype='double' ) !温度減率 call HistoryAddVariable( & & varname='dtdz', dims=(/'p'/), & & longname='Adiabatic Lapse Rate', & & units='K/km', & & xtype='double' ) !安定度 call HistoryAddVariable( & & varname='stab', dims=(/'p'/), & & longname='Static Stability', & & units='s^-2', & & xtype='double' ) ! call HistoryPut( 'Press', Press(RegZMin:RegZMax)) call HistoryPut( 'temp', Temp(RegZMin:RegZMax)) call HistoryPut( 'molfr', MolFr(RegZMin:RegZMax,:)) call HistoryPut( 'dtdz', - DTempDz(RegZMin:RegZMax) * 1.0d3) call HistoryPut( 'stab', Stab(RegZMin:RegZMax)) call HistoryClose end subroutine HistoryFile end subroutine eccm2