! ECCM の再計算 ! ! 内部で利用する式は 中島 (1998) を参照 program main use gt4_history implicit none real(8), allocatable :: Temp(:) real(8), allocatable :: Temp_approx(:) real(8), allocatable :: VTemp(:) real(8), allocatable :: VTemp_approx(:) real(8), allocatable :: dtdz(:) real(8), allocatable :: dtdz_approx(:) real(8), allocatable :: Vdtdz(:) real(8), allocatable :: Vdtdz_approx(:) real(8), allocatable :: vstab(:) real(8), allocatable :: vstab_approx(:) real(8), allocatable :: Pres(:) real(8), allocatable :: Mol(:,:) !モル数 (mol) real(8), allocatable :: MolWt(:) !分子量 (kg/mol) real(8), allocatable :: MolFr(:,:) !モル分率 real(8), allocatable :: Cp(:) !比熱 (J/K mol) real(8), allocatable :: MolWtMean(:) !平均分子量 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 !!! !!!配列の初期化 !!! allocate(Pres( Nstep ), & & vstab( Nstep ), vstab_approx( Nstep ), & & Temp( Nstep ), Temp_approx( Nstep ), & & VTemp( Nstep ), VTemp_approx( Nstep ), & & dtdz( NStep ), dtdz_approx( NStep ), & & vdtdz( NStep ), vdtdz_approx( NStep ), & & Mol( Nstep, s ), MolWt( s ), Cp( s ), & & MolFr( Nstep, s ), MolWtMean( NStep), Fr( s ) ) !温度 Temp = TempIni Temp_approx = TempIni dtdz = 0.0d0 dtdz_approx = 0.0d0 MolFr = 0.0d0 !出力ファイル OutputFile = "eccm.nc" !分子量 (単位は 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 !水(氷) !!! !!!圧力刻みの設定 !!! Pres = PresIni do i = 2, NStep Pres(i) = PresIni & & * 10.0d0**( ( i - 1 ) * dlog10( PresEnd / PresIni ) & & / ( NStep - 1 ) ) end do !!! !!!NetCDF ファイルの用意 !!! call gt4f90io_init !!! !!!断熱減率の計算 !!! do i = 1, NStep - 1 !圧力変化 dp = Pres(i+1) - Pres(i) !前のステップでの気相の全モル数 MolGas = sum(Mol(i, 1:v), 1) !前のステップでのモル分率 MolFr(i, 1:v) = Mol(i, 1:v) / MolGas !前のステップでのモル分率を用いて現在の蒸気圧を計算 VapPres = MolFr(i, e) * Pres(i) !飽和蒸気圧 call Amp( Temp(i), SatVapPres ) if ( VapPres < SatVapPres ) then !平均分子量 (モル数は前のステップと変わらない) MolWtMean(i) = dot_product( MolFr(i, 1:v), MolWt(1:v) ) !平均比熱の計算 (モル数は前のステップと変わらない) CpMean = dot_product( MolFr(i, 1:v), Cp(1:v) ) !乾燥断熱減率を計算 dtdz(i) = - Grav * MolWtMean(i) / CpMean dtdz_approx(i) = - Grav * MolWtMean(i) / CpMean else !飽和蒸気圧と圧力から現在のモル数を計算 Mol(i,e) = SatVapPres * ( MolGas - Mol(i,e) ) & & / ( Pres(i) - SatVapPres) !現在の気相の全モル数 MolGas = sum(Mol(i, 1:v), 1) !現在のモル分率 MolFr(i, 1:v) = Mol(i,1:v) / MolGas !平均分子量 MolWtMean(i) = dot_product( MolFr(i, 1:v), MolWt(1:v) ) !平均比熱の計算 CpMean = dot_product( MolFr(i, 1:v), Cp(1:v) ) !係数 A = LatentHeat * Molfr(i,e) / ( GasR * Temp(i) ) B = ( LatentHeat ** 2.0d0 ) * MolFr(i,e) & & / ( CpMean * GasR * ( Temp(i) ** 2.0d0 ) ) write(*,*) A write(*,*) B !断熱減率 dtdz_approx(i) = - Grav * MolWtMean(i) * ( 1 + A - B ) / CpMean dtdz(i) = - Grav * MolWtMean(i) * ( 1 + A ) / ( CpMean * ( 1 + B ) ) end if Temp(i+1) = Temp(i) + dtdz(i) * dp * ( - GasR ) * Temp(i) & & / ( Pres(i) * Grav * MolWtMean(i) ) Temp_approx(i+1) = Temp_approx(i) & & + dtdz_approx(i) * dp * ( - GasR ) * Temp(i) & & / ( Pres(i) * Grav * MolWtMean(i) ) end do !決まらない部分を与える dtdz( NStep ) = dtdz( NStep - 1 ) dtdz_approx( NStep ) = dtdz_approx( NStep - 1 ) MolFr( NStep, : ) = MolFr( NStep - 1, : ) MolWtMean( NStep ) = MolWtMean( NStep - 1 ) !!! !!!仮温度への変換 !!! !!!平均分子量, モル分率は既に求まっている !乾燥気塊の分子量. 最上層でのモル数を利用 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) * ( MolWtMean(i) / MolWtDry ) ! vTemp_approx(i) = Temp_approx(i) * ( MolWtMean(i) / MolWtDry ) vTemp(i) = Temp(i) * MolWtDry / MolWtMean(i) vTemp_approx(i) = Temp_approx(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_approx(i) = - ( VTemp_approx(i+1) - VTemp_approx(i) ) & & * ( Pres(i) * Grav * MolWtMean(i) )/ ( GasR * Temp(i) * dp ) end do vdtdz(NStep) = vdtdz(NStep - 1) vdtdz_approx(NStep) = vdtdz_approx(NStep - 1) !静的安定度 do i = 1, NStep - 1 vstab(i) = Grav * ( vdtdz(i) - VdtdzDry ) / VTemp(i) vstab_approx(i) = Grav * ( vdtdz_approx(i) - VdtdzDry ) / VTemp(i) write(*,*) vdtdz(i), VdtdzDry, vstab(i) end do vstab(NStep) = vstab(NStep - 1) vstab_approx(NStep) = vstab_approx(NStep - 1) !!! !!!単位の換算 (m --> km) !!! dtdz = dtdz * (- 1.0d3 ) dtdz_approx = dtdz_approx * ( - 1.0d3 ) vdtdz = vdtdz * (- 1.0d3 ) vdtdz_approx = vdtdz_approx * ( - 1.0d3 ) !!! !!!物理量の書き出し !!! call HistoryPut('mol', Mol) call HistoryPut('temp', temp) call HistoryPut('temp_approx', Temp_approx) call HistoryPut('vtemp', vtemp) call HistoryPut('vtemp_approx', vTemp_approx) call HistoryPut('dtdz', dtdz) call HistoryPut('dtdz_approx', dtdz_approx) call HistoryPut('vdtdz', vdtdz) call HistoryPut('vdtdz_approx', vdtdz_approx) call HistoryPut('vstab', vstab) call HistoryPut('vstab_approx', vstab_approx) !!! !!!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_approx", & & dims=(/'pres'/), & & longname="Temperature", & & units='K', & & xtype="double" ) !変数定義, 仮温度 call HistoryAddVariable( & & varname="vtemp", & & dims=(/'pres'/), & & longname="Temperature", & & units='K', & & xtype="double" ) !変数定義, 仮温度 call HistoryAddVariable( & & varname="vtemp_approx", & & 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_approx", & & 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_approx", & & 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_approx", & & dims=(/'pres'/), & & longname="Stability", & & units='s^-2', & & xtype="double" ) end subroutine gt4f90io_init end program main