! ECCM の再計算 ! ! 内部で利用する式は 中島 (1998) を参照 program main use gt4_history implicit none real(8), allocatable :: Temp(:) real(8), allocatable :: Temp_low(:) real(8), allocatable :: Temp_high(:) real(8), allocatable :: VTemp(:) real(8), allocatable :: VTemp_low(:) real(8), allocatable :: VTemp_high(:) real(8), allocatable :: dtdz(:) real(8), allocatable :: dtdz_low(:) real(8), allocatable :: dtdz_high(:) real(8), allocatable :: Vdtdz(:) real(8), allocatable :: Vdtdz_low(:) real(8), allocatable :: Vdtdz_high(:) real(8), allocatable :: vstab(:) real(8), allocatable :: vstab_low(:) real(8), allocatable :: vstab_high(:) real(8), allocatable :: Mol(:,:) !モル数 (mol) real(8), allocatable :: Mol_low(:,:) real(8), allocatable :: Mol_high(:,:) real(8), allocatable :: MolFr(:,:) !モル分率 real(8), allocatable :: MolFr_low(:,:) real(8), allocatable :: MolFr_high(:,:) real(8), allocatable :: MolWtMean(:) !平均分子量 real(8), allocatable :: MolWtMean_low(:) real(8), allocatable :: MolWtMean_high(:) real(8), allocatable :: Pres(:) real(8), allocatable :: MolWt(:) !分子量 (kg/mol) real(8), allocatable :: Cp(:) !比熱 (J/K mol) real(8), allocatable :: Fr(:) real(8) :: VapPres !分圧 real(8) :: SatVapPres !飽和蒸気圧 real(8) :: MolGas !気相の全モル数 real(8) :: MolWtDry !非凝結成分の平均分子量 real(8) :: MolWtWet !凝結成分の平均分子量 real(8) :: A real(8) :: B real(8) :: dp real(8) :: CpMean real(8) :: CpDry real(8) :: VdtdzDry real(8), parameter :: GasR = 8.31d0 !気体定数 (J/k mol) real(8), parameter :: LatentHeat = 40.66d3 !潜熱 (J/mol) real(8), parameter :: Grav = 23.2d0 !木星重力 real(8), parameter :: PresIni = 2.0d6 real(8), parameter :: PresEnd = 1.0d4 real(8), parameter :: TempIni = 4.5d2 integer, parameter :: NStep = 2000 integer, parameter :: e = 3 integer, parameter :: v = 3 integer, parameter :: s = 4 integer :: i character(40) :: OutputFile !!! !!!配列の初期化 !!! !出力ファイル OutputFile = "eccm.nc" !配列の割り当て allocate(Pres( Nstep ), & & vstab( Nstep ), vstab_low( Nstep ), vstab_high( Nstep ), & & Temp( Nstep ), Temp_low( Nstep ), Temp_high( Nstep ), & & VTemp( Nstep ), VTemp_low( Nstep ), VTemp_high( Nstep ), & & dtdz( NStep ), dtdz_low( NStep ), dtdz_high( NStep ), & & vdtdz( NStep ), vdtdz_low( NStep ), vdtdz_high( NStep ), & & Mol( Nstep, s ), Mol_low( Nstep, s ), Mol_high( Nstep, s ), & & MolFr( Nstep, s ), MolFr_low( Nstep, s ), MolFr_high( Nstep, s ), & & MolWtMean( NStep), MolWtMean_low( NStep), MolWtMean_high( NStep), & & MolWt( s ), Cp( s ), Fr( s ) ) !値を初期化する Temp = TempIni; Temp_low = TempIni; Temp_high = TempIni vTemp = TempIni; vTemp_low = TempIni; vTemp_high = TempIni dtdz = 0.0d0; dtdz_low = 0.0d0; dtdz_high = 0.0d0 vdtdz = 0.0d0; vdtdz_low = 0.0d0; vdtdz_high = 0.0d0 vstab = 0.0d0; vstab_low = 0.0d0; vstab_high = 0.0d0 Mol = 0.0d0; Mol_low = 0.0d0; Mol_high = 0.0d0; MolFr = 0.0d0; MolFr_low = 0.0d0; MolFr_high = 0.0d0; MolWtMean = 0.0d0; MolWtMean_low = 0.0d0; MolWtMean_high = 0.0d0; !分子量 (単位は kg/mol) MolWt(1) = 2.0d-3 !水素 MolWt(2) = 4.0d-3 !ヘリウム MolWt(3) = 18.0d-3 !水(水蒸気) MolWt(4) = 18.0d-3 !水(氷) !比熱 Cp = 0.0d0 Cp(1) = 25.0d0 Cp(2) = 20.786d0 Cp(3) = 35.0d0 !初期モル数 Mol(:,1) = 5.0d-1 !水素 Mol(:,2) = 9.75d-2 !ヘリウム ! Mol(:,3) = 8.5d-6 !水(水蒸気) 0.01 x solar ! Mol(:,3) = 8.5d-5 !水(水蒸気) 0.1 x solar ! Mol(:,3) = 8.5d-4 !水(水蒸気) 1 x solar Mol(:,3) = 8.5d-3 !水(水蒸気) 10 x solar ! Mol(:,3) = 8.5d-2 !水(水蒸気) 100 x solar Mol(:,4) = 0.0d0 !水(氷) Mol_low = Mol; Mol_high = Mol !!! !!!圧力刻みの設定 !!! Pres = PresIni do i = 2, NStep Pres(i) = PresIni & & * 10.0d0**( ( i - 1 ) * dlog10( PresEnd / PresIni ) & & / ( NStep - 1 ) ) end do !!! !!!NetCDF ファイルの用意 !!! call gt4f90io_init !!! !!!安定度の計算 !!! call lapserate(Temp, dtdz, Mol, MolFr, MolWtMean, 1) call lapserate(Temp_low, dtdz_low, Mol_low, MolFr_low, MolWtMean_low, 2) call lapserate(Temp_high, dtdz_high, Mol_high, MolFr_high, MolWtMean_high, 3) call stability() !!! !!!仮温度への変換 !!! !!!平均分子量, モル分率は既に求まっている !乾燥気塊の分子量. 最上層でのモル数を利用 MolWtDry = MolWtMean( NStep ) !乾燥気塊の比熱. 非凝縮性成分で考える Fr(1:v-1) = Mol(NStep, 1:v-1) / sum( Mol(NStep, 1:v-1) ) CpDry = dot_product( Fr(1:v-1), Cp(1:v-1) ) do i = 1, NStep !仮温度への変換 vTemp(i) & & = Temp(i) * MolWtDry / MolWtMean(i) vTemp_low(i) & & = Temp_low(i) * MolWtDry / MolWtMean(i) vTemp_high(i) & & = Temp_high(i) * MolWtDry / MolWtMean(i) end do !乾燥気塊の仮温度減率 VdtdzDry = - Grav * MolWtDry / CpDry do i = 1, NStep - 1 !圧力差 dp = Pres(i+1) - Pres(i) !仮温度減率 vdtdz(i) = - ( VTemp(i+1) - VTemp(i) ) & & * ( Pres(i) * Grav * MolWtMean(i) )/ ( GasR * Temp(i) * dp ) vdtdz_low(i) = - ( VTemp_low(i+1) - VTemp_low(i) ) & & * ( Pres(i) * Grav * MolWtMean(i) ) & & / ( GasR * Temp_low(i) * dp ) vdtdz_high(i) = - ( VTemp_high(i+1) - VTemp_high(i)) & & * ( Pres(i) * Grav * MolWtMean(i) ) & & / ( GasR * Temp_high(i) * dp ) end do vdtdz(NStep) = vdtdz(NStep - 1) vdtdz_low(NStep) = vdtdz_low(NStep - 1) vdtdz_high(NStep) = vdtdz_high(NStep - 1) !静的安定度 do i = 1, NStep - 1 vstab(i) = Grav * ( vdtdz(i) - VdtdzDry ) / VTemp(i) vstab_low(i) & & = Grav * ( vdtdz_low(i) - VdtdzDry ) / VTemp_low(i) vstab_high(i) & & = Grav * ( vdtdz_high(i) - VdtdzDry ) / VTemp_high(i) end do vstab(NStep) = vstab(NStep - 1) vstab_low(NStep) = vstab_low(NStep - 1) vstab_high(NStep) = vstab_high(NStep - 1) !!! !!!単位の換算 (m --> km), 符号の付け換え !!! dtdz = dtdz * (- 1.0d3 ) dtdz_low = dtdz_low * ( - 1.0d3 ) dtdz_high = dtdz_high * ( - 1.0d3 ) vdtdz = vdtdz * (- 1.0d3 ) vdtdz_low = vdtdz_low * ( - 1.0d3 ) vdtdz_high = vdtdz_high * ( - 1.0d3 ) !!! !!!物理量の書き出し !!! call HistoryPut('mol', Mol) call HistoryPut('temp', temp) call HistoryPut('vtemp', vtemp) call HistoryPut('dtdz', dtdz) call HistoryPut('vdtdz', vdtdz) call HistoryPut('vstab', vstab) call HistoryPut('temp_low', Temp_low) call HistoryPut('vtemp_low', vTemp_low) call HistoryPut('dtdz_low', dtdz_low) call HistoryPut('vdtdz_low', vdtdz_low) call HistoryPut('vstab_low', vstab_low) call HistoryPut('temp_high', Temp_high) call HistoryPut('vtemp_high', vTemp_high) call HistoryPut('dtdz_high', dtdz_high) call HistoryPut('vdtdz_high', vdtdz_high) call HistoryPut('vstab_high', vstab_high) !!! !!!NetCDF ファイルを閉じる !!! call HistoryClose contains subroutine Antoine(temp, svap) implicit none real(8), intent(in) :: temp !温度 real(8), intent(out) :: svap !飽和蒸気圧 real(8) :: svap_log !以下は水での値 real(8), parameter :: AA = 8.184254d0 real(8), parameter :: AB = 1791.3d0 real(8), parameter :: AC = 238.1d0 !Antoine の式 svap_log = (AA - (AB / (AC + Temp - 273.15D0) ) ) & & * dlog(10.0D0) + dlog(133.322D0) svap = dexp(svap_log) end subroutine Antoine subroutine Amp(temp, svap) implicit none real(8), intent(in) :: temp !温度 real(8), intent(out) :: svap !飽和蒸気圧 real(8) :: svap_log !以下は水での値 real(8), parameter :: A = -5631.1206d0 real(8), parameter :: B = -8.363602d0 real(8), parameter :: C = 8.2312d0 real(8), parameter :: D = -0.03861449d0 real(8), parameter :: E = 2.77494d-5 svap_log = & & A / temp & & + B & & + (C * dlog(temp)) & & + (D * temp) & & + (E * (temp**2)) & & + dlog(1.0D-1) svap = dexp(svap_log) end subroutine Amp subroutine gt4f90io_init implicit none real(8) :: spc(s) real(8) :: elm(s) real(8) :: sizes(20) integer :: i integer :: num character(1000) :: dataset write(dataset, *) "; PresIni=", pres(1), & & "; PresEnd=", pres(NStep), & & "; TempIni=", temp(1), & & "; step=", NStep num = len_trim(dataset) call HistoryCreate(& & file = OutputFile, & & title = "Equilibrium Calculation by ECCM", & & source = dataset(:num+1), & & institution = "GFD_Dennou_Club oboro Project", & & dims = (/'pres', 'spc ', 'elm ', 'size'/), & & dimsizes = (/NStep, s, s, 20/), & & longnames = (/"Pressure ", "chemical species", & & "element ", "size of name "/), & & units = (/"Pa", "1 ", "1 ", "1 "/), & & origin = 0.0, & & interval = 0.0) spc = (/(i, i = 1, s)/) call HistoryPut('spc', real(spc, 4)) elm = (/(i, i = 1, s)/) call HistoryPut('elm', real(elm, 4)) call HistoryPut('pres', real(pres, 4)) sizes = (/(i, i = 1, 20)/) call HistoryPut('size', real(sizes, 4)) !変数定義, 温度 call HistoryAddVariable( & & varname="temp", & & dims=(/'pres'/), & & longname="Temperature", & & units='K', & & xtype="double" ) !変数定義, 温度 call HistoryAddVariable( & & varname="temp_low", & & dims=(/'pres'/), & & longname="Temperature", & & units='K', & & xtype="double" ) !変数定義, 温度 call HistoryAddVariable( & & varname="temp_high", & & dims=(/'pres'/), & & longname="Temperature", & & units='K', & & xtype="double" ) !変数定義, 仮温度 call HistoryAddVariable( & & varname="vtemp", & & dims=(/'pres'/), & & longname="Temperature", & & units='K', & & xtype="double" ) !変数定義, 仮温度 call HistoryAddVariable( & & varname="vtemp_low", & & dims=(/'pres'/), & & longname="Temperature", & & units='K', & & xtype="double" ) !変数定義, 仮温度 call HistoryAddVariable( & & varname="vtemp_high", & & dims=(/'pres'/), & & longname="Temperature", & & units='K', & & xtype="double" ) !変数定義, モル call HistoryAddVariable( & & varname="mol", & & dims=(/'pres', 'spc '/), & & longname="Chemical Species Abundance", & & units='mol', & & xtype="double" ) !変数定義, 温度減率 call HistoryAddVariable( & & varname="dtdz", & & dims=(/'pres'/), & & longname="Lapse Rate", & & units='K/km', & & xtype="double" ) !変数定義, 温度減率 call HistoryAddVariable( & & varname="dtdz_low", & & dims=(/'pres'/), & & longname="Lapse Rate", & & units='K/km', & & xtype="double" ) !変数定義, 温度減率 call HistoryAddVariable( & & varname="dtdz_high", & & dims=(/'pres'/), & & longname="Lapse Rate", & & units='K/km', & & xtype="double" ) !変数定義, 温度減率 call HistoryAddVariable( & & varname="vdtdz", & & dims=(/'pres'/), & & longname="Lapse Rate", & & units='K/km', & & xtype="double" ) !変数定義, 温度減率 call HistoryAddVariable( & & varname="vdtdz_low", & & dims=(/'pres'/), & & longname="Lapse Rate", & & units='K/km', & & xtype="double" ) !変数定義, 温度減率 call HistoryAddVariable( & & varname="vdtdz_high", & & dims=(/'pres'/), & & longname="Lapse Rate", & & units='K/km', & & xtype="double" ) !変数定義, 温度減率 call HistoryAddVariable( & & varname="vstab", & & dims=(/'pres'/), & & longname="Stability", & & units='s^-2', & & xtype="double" ) !変数定義, 温度減率 call HistoryAddVariable( & & varname="vstab_low", & & dims=(/'pres'/), & & longname="Stability", & & units='s^-2', & & xtype="double" ) !変数定義, 温度減率 call HistoryAddVariable( & & varname="vstab_high", & & dims=(/'pres'/), & & longname="Stability", & & units='s^-2', & & xtype="double" ) end subroutine gt4f90io_init !!! !!!断熱減率の計算 !!! subroutine lapserate(Temp_anal, dtdz_anal, Mol_anal, MolFr_anal, MolWtMean_anal, case) implicit none real(8), intent(inout) :: Temp_anal(:) real(8), intent(out) :: dtdz_anal(:) real(8), intent(inout) :: Mol_anal(:,:) real(8), intent(out) :: MolFr_anal(:,:) real(8), intent(out) :: MolWtMean_anal(:) integer, intent(in) :: case integer :: i do i = 1, NStep - 1 !圧力変化 dp = Pres(i+1) - Pres(i) !前のステップでの気相の全モル数 MolGas = sum(Mol_anal(i, 1:v), 1) !前のステップでのモル分率 MolFr_anal(i, 1:v) = Mol_anal(i, 1:v) / MolGas !前のステップでのモル分率を用いて現在の蒸気圧を計算 VapPres = MolFr_anal(i, e) * Pres(i) !飽和蒸気圧 call Amp( Temp_anal(i), SatVapPres ) if ( VapPres < SatVapPres ) then !平均分子量 (モル数は前のステップと変わらない) MolWtMean_anal(i) = dot_product( MolFr_anal(i, 1:v), MolWt(1:v) ) !平均比熱の計算 (モル数は前のステップと変わらない) CpMean = dot_product( MolFr_anal(i, 1:v), Cp(1:v) ) !乾燥断熱減率を計算 dtdz_anal(i) = - Grav * MolWtMean_anal(i) / CpMean else !飽和蒸気圧と圧力から現在のモル数を計算 Mol_anal(i,e) = SatVapPres * ( MolGas - Mol_anal(i,e) ) & & / ( Pres(i) - SatVapPres) !現在の気相の全モル数 MolGas = sum(Mol_anal(i, 1:v), 1) !現在のモル分率 MolFr_anal(i, 1:v) = Mol_anal(i,1:v) / MolGas !平均分子量 MolWtMean_anal(i) = dot_product( MolFr_anal(i, 1:v), MolWt(1:v) ) !平均比熱の計算 CpMean = dot_product( MolFr_anal(i, 1:v), Cp(1:v) ) !係数 A = LatentHeat * Molfr(i,e) / ( GasR * Temp_anal(i) ) B = ( LatentHeat ** 2.0d0 ) * MolFr_anal(i,e) & & / ( CpMean * GasR * ( Temp_anal(i) ** 2.0d0 ) ) if (case == 1) then !断熱減率 (一般形) dtdz_anal(i) & & = - Grav * MolWtMean_anal(i) * ( 1 + A ) & & / ( CpMean * ( 1 + B ) ) elseif (case == 2) then !断熱減率 (凝結成分が少ない近似) dtdz_anal(i) & & = - Grav * MolWtMean_anal(i) * ( 1 + A - B ) / CpMean elseif (case == 3) then !断熱減率 (凝結成分が多い近似) dtdz_anal(i) & & = - Grav * MolWtMean_anal(i) * A / ( CpMean * B ) end if end if Temp_anal(i+1) = Temp_anal(i) & & + dtdz_anal(i) * dp * ( - GasR ) * Temp_anal(i) & & / ( Pres(i) * Grav * MolWtMean_anal(i) ) end do !決まらない部分を与える dtdz_anal( NStep ) = dtdz_anal( NStep - 1 ) Mol_anal( NStep, : ) = Mol_anal( NStep - 1, : ) MolFr_anal( NStep, : ) = MolFr_anal( NStep - 1, : ) MolWtMean_anal( NStep ) = MolWtMean_anal( NStep - 1 ) end subroutine lapserate end program main