Skip to content
GitLab
Commit d231b10c authored by Claire Guimond's avatar Claire Guimond
Browse files

updated atmosphere.f90

parent 17b440bf
Hide whitespace changes
Inline Side-by-side
Octave/make_movie.m
View file @ d231b10c
function make_movie(fpath)
[fname, fpath, fltidx] = uigetfile("data_char_val_ts.res");
filenames = sprintf('%s%s', fpath, '/data_snap*')
......
Octave/plot_snap.m
View file @ d231b10c
......@@ -285,7 +285,7 @@ function produce_plot(name,title_plot,field,input,data,res,ar,nr_contour,cmap,mi
else
colormap(cmap)
endif
caxis([min_z max_z]);
%caxis([min_z max_z]);
colorbar
FN = findall(f,'-property','FontName');
......
src/atmosphere.f90
View file @ d231b10c
......@@ -40,6 +40,7 @@ contains
real(dp)::pi__,FD ! FD correspond to Final Dimensionalization
real(dp):: KI, KII, X_carbonate_melt , fO2
integer::i,j,k,m
real(dp)::M_tot,X_mol_H2,X_mol_H2O,X_mol_CO,X_mol_CO2,sum_moles,frac_tot_H2,frac_tot_H2O,frac_tot_CO,frac_tot_CO2
call atmos_outgassing(t,part,mesh,cell,mc,field,di,pF,pG,pM,pN,pP,pX,Psurf)
if (pX%use_atmos.ge.2) call atmos_escaping(t,mc,pF,pM,mesh)
......@@ -54,7 +55,27 @@ contains
! no escape included, but different species outgassed
if (t.eq.dt) write(60,*) ' Time (Gyr) | T_mantle(K) | maxF | mass H2O |&
& total H2O | mass CO2 | total CO2 | mass H2 | total H2 | mass CO | total CO |&
& PH2O(bar) | PCO2(bar) | PH2(bar) | PCO(bar) | EGL_H2O | frac_CO2 | frac_CO | frac_H2O | frac_H2'
& PH2O(bar) | PCO2(bar) | PH2(bar) | PCO(bar) | EGL_H2O | frac_H2O | frac_CO2 | frac_H2 | frac_CO'
M_tot = mc%total_mass_H2O + mc%total_mass_CO2 + mc%total_mass_H2 + mc%total_mass_CO
X_mol_H2 = mc%total_mass_H2*1000.0_dp / 2.01588_dp ! number of moles
X_mol_H2O = mc%total_mass_H2O*1000.0_dp / 18.01488_dp
X_mol_CO = mc%total_mass_CO*1000.0_dp / 28.0097_dp
X_mol_CO2 = mc%total_mass_CO2*1000.0_dp / 44.0087_dp
sum_moles = X_mol_H2+X_mol_H2O+X_mol_CO+X_mol_CO2
if (sum_moles.gt.0.0_dp) then
frac_tot_H2 = X_mol_H2 / sum_moles
frac_tot_H2O = X_mol_H2O / sum_moles
frac_tot_CO = X_mol_CO / sum_moles
frac_tot_CO2 = X_mol_CO2 / sum_moles
else
frac_tot_H2 = 0.0_dp
frac_tot_H2O = 0.0_dp
frac_tot_CO = 0.0_dp
frac_tot_CO2 = 0.0_dp
endif
write(60,'(F13.8,F12.3,F13.8,ES16.4,ES16.4,ES16.4,ES16.4,ES16.4,ES16.4,ES16.4,ES16.4,&
& F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6)') &
......@@ -69,12 +90,16 @@ contains
& mc%total_mass_H2*FD,&
& mc%mass_CO*FD,&
& mc%total_mass_CO*FD,&
& mc%total_mass_H2O*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed H2O accumulating over time in bar
& mc%total_mass_CO2*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed CO2 accumulating over time in bar
& mc%total_mass_H2*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed H2 accumulating over time in bar
& mc%total_mass_CO*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed CO accumulating over time in bar
! & mc%total_mass_H2O*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed H2O accumulating over time in bar
! & mc%total_mass_CO2*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed CO2 accumulating over time in bar
! & mc%total_mass_H2*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed H2 accumulating over time in bar
! & mc%total_mass_CO*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed CO accumulating over time in bar
& M_tot * frac_tot_H2O*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed H2O accumulating over time in bar
& M_tot * frac_tot_CO2*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed CO2 accumulating over time in bar
& M_tot * frac_tot_H2 *pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed H2 accumulating over time in bar
& M_tot * frac_tot_CO *pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*FD,& ! pressure outgassed CO accumulating over time in bar
& (mc%total_mass_H2O*FD - mc%escape_H2O)/(4.0_dp*pi__*mesh%rmax**2*pF%D**2*1000.0_dp),& ! EGL H2O in m, assuming rho_water=1000 !X_atm_H2O*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5)*1e5_dp/(pF%g*1000.0_dp) ! EGL H2O
& mc%frac_CO2,mc%frac_CO,mc%frac_H2O,mc%frac_H2
& mc%frac_H2O,mc%frac_CO2,mc%frac_H2,mc%frac_CO ! mole fractions
else
......@@ -129,7 +154,7 @@ contains
integer,intent(in)::di
real(dp)::KI, KII, X_carbonate_melt,T_K,Tm_surf,P_bar,Psurf,fO2,D,pi__
real(dp)::frac_CO2,frac_CO,frac_H2O,frac_H2,mole_frac_H2O,mole_frac_CO2
real(dp)::frac_CO2,frac_CO,frac_H2O,frac_H2,mole_frac_H2O,mole_frac_CO2,vol
real(dp)::M_H2O,M_CO2,M_H2,M_CO
integer::i,j,k,m,jmin,jmax,y
pi__ = acos(-1.0_dp) !3.1415
......@@ -201,6 +226,12 @@ contains
mc%mass_H2=0.0_dp
mc%mass_CO=0.0_dp
mc%frac_CO = 0.0_dp
mc%frac_CO2 = 0.0_dp
mc%frac_H2 = 0.0_dp
mc%frac_H2O = 0.0_dp
vol = 0.0_dp
M_H2O = 2.0_dp*1.00794_dp + 15.999_dp
M_CO2 = 12.0107_dp + 2.0_dp*15.999_dp
M_H2 = 2.0_dp*1.00794_dp
......@@ -227,22 +258,27 @@ contains
! sum of mole_frac_H2O and mole_frac_CO2 needs to be one
mole_frac_H2O = cell%X_H2O(i,j,k)/M_H2O/(cell%X_H2O(i,j,k)/M_H2O+cell%X_CO2(i,j,k)/M_CO2)
mole_frac_CO2 = cell%X_CO2(i,j,k)/M_CO2/(cell%X_H2O(i,j,k)/M_H2O+cell%X_CO2(i,j,k)/M_CO2)
call spec_COH_Fegley2013(mole_frac_H2O,mole_frac_CO2,pX%fO2,Psurf,Tm_surf,mc%frac_CO2,&
&mc%frac_CO,mc%frac_H2O,mc%frac_H2)
call spec_COH_Fegley2013(mole_frac_H2O,mole_frac_CO2,pX%fO2,Psurf,Tm_surf,frac_CO2,&
&frac_CO,frac_H2O,frac_H2)
! frac_CO2+frac_CO+frac_H2O+frac_H2 is equal to one
mc%frac_CO = mc%frac_CO + frac_CO * mesh%dVi(i,j,k)
mc%frac_CO2 = mc%frac_CO2 + frac_CO2 * mesh%dVi(i,j,k)
mc%frac_H2 = mc%frac_H2 + frac_H2 * mesh%dVi(i,j,k)
mc%frac_H2O = mc%frac_H2O + frac_H2O * mesh%dVi(i,j,k)
vol = vol + mesh%dVi(i,j,k)
! calc mass H2O/CO2/H2/CO each in kg (when dimensionalized, here still dimensionless value)
! first calculate mole of H2O and CO2, use frac_X to get moles of each species, multiply with molar masses to get final gas masses
if (mole_frac_H2O.gt.0.0_dp) then
mc%mass_H2O = mc%mass_H2O + (cell%X_H2O(i,j,k)/M_H2O * mc%frac_H2O/mole_frac_H2O * M_H2O)&
mc%mass_H2O = mc%mass_H2O + (cell%X_H2O(i,j,k)/M_H2O * frac_H2O/mole_frac_H2O * M_H2O)&
& * mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)*cell%ref_r(i,j,k)*(1.0e-6_dp)
mc%mass_H2 = mc%mass_H2 + (cell%X_H2O(i,j,k)/M_H2O * mc%frac_H2/mole_frac_H2O * M_H2)&
mc%mass_H2 = mc%mass_H2 + (cell%X_H2O(i,j,k)/M_H2O * frac_H2/mole_frac_H2O * M_H2)&
& * mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)*cell%ref_r(i,j,k)*(1.0e-6_dp)
end if
if (mole_frac_CO2.gt.0.0_dp) then
mc%mass_CO2 = mc%mass_CO2 + (cell%X_CO2(i,j,k)/M_CO2 * mc%frac_CO2/mole_frac_CO2 * M_CO2)&
mc%mass_CO2 = mc%mass_CO2 + (cell%X_CO2(i,j,k)/M_CO2 * frac_CO2/mole_frac_CO2 * M_CO2)&
& * mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)*cell%ref_r(i,j,k)*(1.0e-6_dp)
mc%mass_CO = mc%mass_CO + (cell%X_CO2(i,j,k)/M_CO2 * mc%frac_CO/mole_frac_CO2 * M_CO)&
mc%mass_CO = mc%mass_CO + (cell%X_CO2(i,j,k)/M_CO2 * frac_CO/mole_frac_CO2 * M_CO)&
& * mc%DeltaF(i,j,k)*mesh%dVi(i,j,k)*cell%ref_r(i,j,k)*(1.0e-6_dp)
end if
else
......@@ -263,10 +299,15 @@ contains
enddo
enddo
enddo
mc%mass_H2O = mc%mass_H2O * D ! in Kg
mc%mass_H2O = mc%mass_H2O * D ! in nondimensional Kg
mc%mass_CO2 = mc%mass_CO2 * D
mc%mass_H2 = mc%mass_H2 * D
mc%mass_CO = mc%mass_CO * D
mc%frac_CO = mc%frac_CO / vol
mc%frac_CO2 = mc%frac_CO2 / vol
mc%frac_H2 = mc%frac_H2 / vol
mc%frac_H2O = mc%frac_H2O / vol
mc%total_mass_H2O = mc%total_mass_H2O + mc%mass_H2O*pM%Chi_extr
mc%total_mass_CO2 = mc%total_mass_CO2 + mc%mass_CO2*pM%Chi_extr
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment