Skip to content
GitLab
Commit 3d7155c5 authored by Lena Noack's avatar Lena Noack
Browse files

Update mantle volatiles output

parent 268c4cd3
Hide whitespace changes
Inline Side-by-side
Octave/plot_snap.m
View file @ 3d7155c5
......@@ -25,7 +25,7 @@ function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
mkdir "plots";
graphics_toolkit gnuplot
res=5000; %1200; %600;
res=1200; %600;
ar=(max(max(data.p_xx_dim))-min(min(data.p_xx_dim)))/(max(max(data.p_zz_dim))-min(min(data.p_zz_dim)));
time = data.time/data.t_yr/1000000.0
......@@ -319,4 +319,4 @@ function produce_plot(name,title_plot,field,input,data,res,ar,nr_contour,cmap,mi
clear f;
end
\ No newline at end of file
end
src/atmosphere.f90
View file @ 3d7155c5
......@@ -53,10 +53,9 @@ contains
if (pX%gas_spec.ge.1) then
! 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_H2O | frac_CO2 | frac_H2 | frac_CO'
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_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
......@@ -187,6 +186,10 @@ contains
if (part%X_H2O(m)*part%df(m).gt.part%w(m)) then ! condition for having more water after than before the calculation
part%X_H2O(m)= part%w(m)/part%df(m)
end if
! if ((i.eq.1).and.(j.eq.1).and.(k.eq.mesh%nr+1)) &
! &write(*,*) part%w(m),part%X_H2O(m),part%w(m) - (part%X_H2O(m)*part%df(m)),field%w(i,j,k)
part%w(m) = part%w(m) - (part%X_H2O(m)*part%df(m))
! CO2 equations (Holloway)
......
src/main.f90
View file @ 3d7155c5
......@@ -117,7 +117,7 @@ program CHIC
integer::magma_local,magma_global
real(dp)::magma !time when first post-accretional magma ocean appears
real(dp)::save_lmin,save_lmax,save_ymin,save_ymax,save_rmin,save_rmax
real(dp)::depl_av,water_av,DeltaF_av,vol,vol_tot,heat_m
real(dp)::depl_av,water_av,carbon_av,DeltaF_av,vol,vol_tot,heat_m
character(len=10):: str_i
real(dp),allocatable::T_profile(:),v_profile(:),visc_profile(:),T_center(:),save_rightBnd(:,:),nmes_ch(:,:),nmes_me(:,:),&
&nmes_li(:,:),nmes_ch1(:),nmes_me1(:),nmes_li1(:),Tprof(:),Tprof0(:),r(:),OcTProf(:,:)
......@@ -1007,6 +1007,7 @@ program CHIC
depl_av = 0.0_dp
water_av = 0.0_dp
carbon_av = 0.0_dp
DeltaF_av = 0.0_dp
vol_tot = 0.0_dp
do ind_i = 1,mesh%nl
......@@ -1015,6 +1016,7 @@ program CHIC
vol = 4.0_dp/3.0_dp * acos(-1.0_dp) * (mesh%rb(ind_k+1)**3 - mesh%rb(ind_k)**3)
depl_av = depl_av + cell%d(ind_i,ind_j,ind_k)*vol
water_av = water_av + cell%w(ind_i,ind_j,ind_k)*vol
carbon_av = carbon_av + cell%CO2(ind_i,ind_j,ind_k)*vol
DeltaF_av = DeltaF_av + mc%DeltaF(ind_i,ind_j,ind_k)*vol
vol_tot = vol_tot + vol
enddo
......@@ -1022,6 +1024,7 @@ program CHIC
enddo
depl_av = depl_av/vol_tot
water_av = water_av/vol_tot
carbon_av = carbon_av/vol_tot
DeltaF_av = DeltaF_av/vol_tot
if (pM%add_crust.ne.0) then
......@@ -1045,7 +1048,7 @@ program CHIC
&sum(mc%cc_vol)/(mesh%lmax-mesh%lmin)/(mesh%rmax-mesh%rmin),&
&mc%mass_H2O,mc%total_mass_H2O,water_av,&
&minval(mc%DeltaF),DeltaF_av,maxval(mc%DeltaF),&
&field%T(1,1,ind_k),mc%DeltaF(1,1,ind_k),field%w(1,1,ind_k+1),0.0_dp,0.0_dp,0.0_dp
&field%T(1,1,ind_k),mc%DeltaF(1,1,ind_k),field%w(1,1,ind_k+1),carbon_av,minval(cell%w),minval(cell%CO2)
endif
endif
......@@ -1653,7 +1656,6 @@ program CHIC
!!!!!!!!!!!!!!! END INNER ITERATION !!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! if (pW%OS_SuL.gt.0.0_dp) then
! do ind_k = 0,mesh%nr+1
! if ((mesh%rmax-mesh%rc(ind_k)).le.pW%OS_SuL) then ! in ocean
......@@ -2169,6 +2171,7 @@ program CHIC
depl_av = 0.0_dp
water_av = 0.0_dp
carbon_av = 0.0_dp
DeltaF_av = 0.0_dp
vol_tot = 0.0_dp
do ind_i = 1,mesh%nl
......@@ -2177,6 +2180,7 @@ program CHIC
vol = 4.0_dp/3.0_dp * acos(-1.0_dp) * (mesh%rb(ind_k+1)**3 - mesh%rb(ind_k)**3)
depl_av = depl_av + cell%d(ind_i,ind_j,ind_k)*vol
water_av = water_av + cell%w(ind_i,ind_j,ind_k)*vol
carbon_av = carbon_av + cell%CO2(ind_i,ind_j,ind_k)*vol
DeltaF_av = DeltaF_av + mc%DeltaF(ind_i,ind_j,ind_k)*vol
vol_tot = vol_tot + vol
enddo
......@@ -2184,6 +2188,7 @@ program CHIC
enddo
depl_av = depl_av/vol_tot
water_av = water_av/vol_tot
carbon_av = carbon_av/vol_tot
DeltaF_av = DeltaF_av/vol_tot
if (pM%add_crust.ne.0) then
......@@ -2210,10 +2215,9 @@ program CHIC
&sum(mc%cc_vol)/(mesh%lmax-mesh%lmin)/(mesh%rmax-mesh%rmin),&
&mc%mass_H2O,mc%total_mass_H2O,water_av,&
&minval(mc%DeltaF),DeltaF_av,maxval(mc%DeltaF),&
&field%T(1,1,ind_k),mc%DeltaF(1,1,ind_k),field%w(1,1,ind_k+1),0.0_dp,0.0_dp,0.0_dp
&field%T(1,1,ind_k),mc%DeltaF(1,1,ind_k),field%w(1,1,ind_k+1),carbon_av,minval(cell%w),minval(cell%CO2) ! attention: field%w is radial velocity, not water content
endif
!write(*,*) minval(mc%DeltaF),sum(mc%DeltaF)/(mesh%nl+2)/(mesh%nr+2),&
! & sum(mc%DeltaF(1:mesh%nl,1,1:mesh%nr))/real(mesh%nl*mesh%nr),maxval(mc%DeltaF)
......@@ -2221,6 +2225,7 @@ program CHIC
depl_av = 0.0_dp
water_av = 0.0_dp
carbon_av = 0.0_dp
DeltaF_av = 0.0_dp
vol_tot = 0.0_dp
do ind_i = 1,mesh%nl
......@@ -2229,6 +2234,7 @@ program CHIC
vol = 4.0_dp/3.0_dp * acos(-1.0_dp) * (mesh%rb(ind_k+1)**3 - mesh%rb(ind_k)**3)
depl_av = depl_av + cell%d(ind_i,ind_j,ind_k)*vol
water_av = water_av + cell%w(ind_i,ind_j,ind_k)*vol
carbon_av = carbon_av + cell%CO2(ind_i,ind_j,ind_k)*vol
DeltaF_av = DeltaF_av + mc%DeltaF(ind_i,ind_j,ind_k)*vol
vol_tot = vol_tot + vol
enddo
......@@ -2236,6 +2242,7 @@ program CHIC
enddo
depl_av = depl_av/vol_tot
water_av = water_av/vol_tot
carbon_av = carbon_av/vol_tot
DeltaF_av = DeltaF_av/vol_tot
nmes_me(nmes_iter,1) = t
......@@ -2262,9 +2269,9 @@ program CHIC
nmes_me(nmes_iter,14) = field%T(1,1,ind_k)
nmes_me(nmes_iter,15) = mc%DeltaF(1,1,ind_k)
nmes_me(nmes_iter,16) = field%w(1,1,ind_k+1)
nmes_me(nmes_iter,17) = 0.0_dp
nmes_me(nmes_iter,18) = 0.0_dp
nmes_me(nmes_iter,19) = 0.0_dp
nmes_me(nmes_iter,17) = carbon_av !0.0_dp
nmes_me(nmes_iter,18) = minval(cell%w) !0.0_dp
nmes_me(nmes_iter,19) = minval(cell%CO2) !0.0_dp
!nmes_me(nmes_iter,20)= mc%mass_CO2
!nmes_me(nmes_iter,21) = mc%total_mass_CO2
endif
......
src/melt.f90
View file @ 3d7155c5
......@@ -162,6 +162,9 @@ contains
part%df(m) = 0.0_dp
endif
if ((i.eq.1).and.(j.eq.1).and.(k.eq.mesh%nr+1)) &
& write(*,*) "Melt p 2:",part%df(m),part%w(m)
! if ((part%df(m).gt.(pM%max_depl-part%d(m))).and.(pM%max_depl.ne.0.0_dp)) then ! more melt than defined in pM%max_depl is not allowed
! part%df(m) = pM%max_depl-part%d(m)
! endif
......@@ -225,10 +228,12 @@ contains
endif
do m=1,part%n
! if (part%d(m).lt.dehydr_melt) then ! immediate dehydration when melt occurs if dehydr_melt = depl_max; no dehydration at all if dehydr_melt = 0
if ((dehydr_melt.gt.0.0_dp).and.(part%d(m).gt.dehydr_melt)) then ! dehydr_melt is melt in percent where all water is dehydrated, e.g. 1-5%
part%w(m) = 0.0_dp ! Complete dehydration
endif
enddo
vol_H2O = 0.0_dp
......
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