Commit fe248efc authored by Lena Noack's avatar Lena Noack
Browse files

Update outgassing, now also for regional sphere, add H2O/CO2 comparision pressures

parent 6f70a019
This diff is collapsed.
function [H2_mole_fraction, H2O_mole_fraction, CO_mole_fraction, CO2_mole_fraction]=ortenzi_simultaneous_eq_C_O_H(temperature, P, H2O_mole_fraction, CO2_mole_fraction,redox_buffer)
function [H2_mole_fraction, H2O_mole_fraction, CO_mole_fraction, CO2_mole_fraction]=ortenzi_simultaneous_eq_C_O_H(temperature, P, H2O_mole_fraction, CO2_mole_fraction, redox_buffer, buffer)
if nargin<5
redox=0;
redox_buffer=0;
end
% redox_buffer -> which buffer to use
% buffer -> add +/- log10 values to buffer, i.e. IW+1 -> buffer = 1
%simultaneous equilibria Fegley 2013 pag. 458
......@@ -25,7 +27,9 @@ T=temperature;
buffer=0;
if nargin<6
buffer=0;
end
%Equation for the calculation of the delta G of formation of carbon monoxyde 2C (graphite) + 2H2 = CH4 (Fegley 2013)
......
......@@ -111,14 +111,13 @@ function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
endif
endif
% if(size(data.depl,1)>1), produce_plot('D2Dc_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),-data.depl,input,data,res,ar,50,'jet',min_z=-max_depl,max_z=0,full=fullV), endif
% if(size(data.depl,1)>1), produce_plot('D2D_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),data.depl,input,data,res,ar,20,'bgry',min_z=0,max_z=max_depl,full=fullV), endif
%
% if(size(data.depl,1)>1), produce_plot('D2Db_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),data.depl,input,data,res,ar,20,'hot',min_z=0,max_z=max_depl,full=fullV), endif
if(size(data.depl,1)>1), produce_plot('D2Dc_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),-data.depl,input,data,res,ar,50,'jet',min_z=-max_depl,max_z=0,full=fullV), endif
if(size(data.depl,1)>1), produce_plot('D2D_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),data.depl,input,data,res,ar,20,'bgry',min_z=0,max_z=max_depl,full=fullV), endif
if(size(data.depl,1)>1), produce_plot('D2Db_',sprintf(strcat("Depletion at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),data.depl,input,data,res,ar,20,'hot',min_z=0,max_z=max_depl,full=fullV), endif
if (reduc==0)
if(size(data.DeltaF,1)>1), produce_plot('DF2D_',"Melt fraction",data.DeltaF,input,data,res,ar,20,'jet',min_z=-1,max_z=-1,full=fullV), endif
% if(size(data.w,1)>1), produce_plot('W2D_',"Water content in ppm",-data.w,input,data,res,ar,50,'jet',min_z=-100,max_z=0,full=fullV), endif %input.H2O), endif
if(size(data.w,1)>1), produce_plot('W2D_',"Water content in ppm",-data.w,input,data,res,ar,50,'jet',min_z=-100,max_z=0,full=fullV), endif %input.H2O), endif
if(size(data.ref_a,1)>1), produce_plot('Al2D_',"Heat expansion",data.ref_a,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
if(size(data.rec,1)>1), produce_plot('Rec2D_',"Recycling",data.rec,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif %input.nr_ph
......
......@@ -41,7 +41,7 @@ contains
real(dp):: KI, KII, X_carbonate_melt , fO2
integer::i,j,k,m,jmin,jmax
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
real(dp)::M_melt,M_w_loss,M_w_degas
real(dp)::M_melt,M_w_loss,M_w_degas,M_c_loss,M_c_degas,vol
call atmos_outgassing(t,part,mesh,cell,mc,field,di,pF,pG,pM,pN,pP,pX,Psurf,M_melt)
if (pX%use_atmos.ge.2) call atmos_escaping(t,mc,pF,pM,mesh)
......@@ -55,31 +55,49 @@ contains
jmax = mesh%ny+1
endif
pi__ = acos(-1.0_dp)
FD=(pF%D**3)*pF%rho ! FD correspond to Final Dimensionalization; values in atmos_outgassing scaled up to 3D sphere, hence "*D^3"
M_w_loss = 0.0_dp
do i=0,mesh%nl+1
do j=jmin,jmax
do k=0,mesh%nr+1
M_w_loss = M_w_loss + (pI%mantle_w_H2O-cell%w(i,j,k))*mesh%dVi(i,j,k)*cell%ref_r(i,j,k)*(1.0e-6_dp)&
& * 4.0_dp/3.0_dp * (mesh%rmax**(3)-mesh%rmin**(3)) / (mesh%rmax**(2)-mesh%rmin**(2)) * pM%Chi_extr
M_c_loss = 0.0_dp
vol = 0.0_dp
do i=1,mesh%nl
do j=1,mesh%ny
do k=1,mesh%nr
! masses of H2O and CO2 in kg that were lost from mantle to atmosphere (only looking at extrusive melt here)
!M_w_loss = M_w_loss + (pI%mantle_w_H2O-cell%w(i,j,k))*mesh%dVi(i,j,k)*cell%ref_r(i,j,k)*(1.0e-6_dp)&
! & * 4.0_dp/3.0_dp * (mesh%rmax**(3)-mesh%rmin**(3)) / (mesh%rmax**(2)-mesh%rmin**(2)) * pM%Chi_extr * FD
!M_c_loss = M_c_loss + (pX%CO2_ini-cell%CO2(i,j,k))*mesh%dVi(i,j,k)*cell%ref_r(i,j,k)*(1.0e-6_dp)&
! & * 4.0_dp/3.0_dp * (mesh%rmax**(3)-mesh%rmin**(3)) / (mesh%rmax**(2)-mesh%rmin**(2)) * pM%Chi_extr * FD
! fraction of H2O and CO2 lost from mantle in ppm:
M_w_loss = M_w_loss + (pI%mantle_w_H2O-cell%w(i,j,k))*mesh%dVi(i,j,k)
M_c_loss = M_c_loss + (pX%CO2_ini-cell%CO2(i,j,k))*mesh%dVi(i,j,k)
vol = vol + mesh%dVi(i,j,k)
enddo
enddo
enddo
M_w_degas = mc%total_mass_H2O + mc%total_mass_H2/2.01588_dp*18.01488_dp ! what would be mass if only H2O, not H2 would be degassed
M_w_loss = M_w_loss / vol
M_c_loss = M_c_loss / vol
!M_w_degas = mc%total_mass_H2O + mc%total_mass_H2/2.01588_dp*18.01488_dp ! what would be mass if only H2O, not H2 would be degassed
!write(*,*) "Water taken out of mantle=",M_w_loss,", Water put in atmosphere=",M_w_degas
! Degassing pressure in bar
M_w_degas = M_w_loss * sum(cell%ref_r(1,1,:))/mesh%nr * (1.0e-6_dp) &
& * 4.0_dp/3.0_dp * pi__ * (mesh%rmax**(3)-mesh%rmin**(3)) &
& * pM%Chi_extr * FD * pI%gravity*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5) ! pressure in bar if only H2O would go out (not actual partial pressure, which depends on mean molecular weigth)
M_c_degas = M_c_loss * sum(cell%ref_r(1,1,:))/mesh%nr * (1.0e-6_dp) &
& * 4.0_dp/3.0_dp * pi__ * (mesh%rmax**(3)-mesh%rmin**(3)) &
& * pM%Chi_extr * FD * pI%gravity*pF%g/(4.0_dp*pi__*mesh%rmax**2*pF%D**2)*10.0_dp**(-5) ! pressure in bar if only CO2 would go out
! output per time step
open(unit=60,file="data_atmosphere.res",form='formatted',status='unknown',POSITION='APPEND')
pi__ = acos(-1.0_dp)
FD=(pF%D**3)*pF%rho ! FD correspond to Final Dimensionalization; values in atmos_outgassing scaled up to 3D sphere, hence "*D^3"
! 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) | mass Melt | mass H2O | &
if (t.eq.dt) write(60,*) ' Time(Gyr) | T_mantle(K) | mass Melt | 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 | &
& X_H2O_melt | X_CO2_melt '
& X_H2O_melt | X_CO2_melt | Deg_H2O(ppm) | Deg_CO2(ppm) | Deg_PH2O(bar) | Deg_PCO2(bar) '
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
......@@ -104,7 +122,7 @@ contains
! & ", flux in kg/yr=",M_tot*FD/dt*pF%time_yr
write(60,'(F13.8,F12.3,ES16.4,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,F12.6,F12.6)') &
& F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6)') &
& t/pF%time_yr*1.0e-9_dp,& ! time in Gyr
& T_mean*pF%DeltaT+pF%Ts,& ! average mantle temperature in K
& M_melt*FD,& ! maxval(part%df),& ! to see when melting starts
......@@ -127,7 +145,8 @@ contains
& (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_H2O,mc%frac_CO2,mc%frac_H2,mc%frac_CO,& ! mole fractions
& sum(cell%X_H2O(1:mesh%nl,1,1:mesh%nr))/(mesh%nl*mesh%nr), &
& sum(cell%X_CO2(1:mesh%nl,1,1:mesh%nr))/(mesh%nl*mesh%nr)
& sum(cell%X_CO2(1:mesh%nl,1,1:mesh%nr))/(mesh%nl*mesh%nr), &
& M_w_loss, M_c_loss, M_w_degas, M_c_degas
! else
!
......@@ -181,7 +200,7 @@ contains
type(param_X),intent(in)::pX
integer,intent(in)::di
real(dp)::KI, KII, X_carbonate_melt,T_K,Tm_surf,P_bar,Psurf,fO2,D,pi__
real(dp)::KI, KII, Kf, 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,vol
real(dp)::M_H2O,M_CO2,M_H2,M_CO,M_melt,vol_cell
integer::i,j,k,m,jmin,jmax,y
......@@ -210,31 +229,16 @@ contains
j = part%cell_y(m)
k = part%cell_r(m)
!ToDo - is it problem that particle calc. concentration of volatiles in melt, without actually havng a mass/volume
! themselves? Hence if 1 particle has tiny melt fraction, and one intermediate melt fraction, concentration in
! melt is linear average of both, but should be exponentially averaged?
! Hence better to extract water below for cell based on average DeltaF, and subtract water equally for all
! particles in cell (if possible)?
! But if particle has less water already, then need to extract more from other cell particles
! BUT: can't be problem, since with dH2O=1 this would not be a problem, but still jumps in DF occur
! H2O equation ! CHECK: better to use aggregate formulation instead of batch melting?
part%X_H2O(m) = part%w(m)/( pX%dH2O + part%df(m)*(1.0_dp-pX%dH2O) ) ! cell%X_H2O / part%X_H2O is concentration in melt
part%X_H2O(m) = part%X_H2O(m)*part%df(m)
! if (part%X_H2O(m)*part%df(m).gt.part%w(m)) then ! condition for not extracting more water than available
! part%X_H2O(m)= part%w(m)/part%df(m)
! end if
! should not be necessary for H2O, but to be on the safe side...
if (part%X_H2O(m).gt.part%w(m)) then ! condition for not extracting more water than available
part%X_H2O(m)= part%w(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)),cell%w(i,j,k)
part%w(m) = part%w(m) - part%X_H2O(m) !part%X_H2O(m)*part%df(m) (X_H2O already includes melt fraction now)
part%w(m) = part%w(m) - part%X_H2O(m) ! (X_H2O already includes melt fraction)
! CO2 equations (Holloway)
T_K= (field%T(i,j,k)*pF%DeltaT)+pF%Ts
......@@ -242,22 +246,19 @@ contains
KI= 40.07639_dp-((2.53932_dp*10.0_dp**(-2))*T_K)+(5.27096_dp*10.0_dp**(-6))*(T_K**2)+0.0267_dp*(P_bar-1.0_dp)/T_K
KII= -6.24763_dp-(282.56_dp/T_K)-(0.119242_dp*(P_bar-1000.0_dp)/T_K)
KI=10.0_dp**(KI) ! for the inverse of LOG10
KII= 10.0_dp**(KII)
fO2= 6.899_dp-(27714.0_dp/T_K)+0.05_dp*(P_bar-1.0_dp)/T_K + pX%fO2 ! instead fO2-field, where with change in O atoms we calculate with the same formula back how fO2 changes locally
fO2= 10.0_dp**(fO2)
X_carbonate_melt = (KI*KII*fO2)/(1.0_dp+(KI*KII*fO2))
! part%X_CO2(m)=((44.01_dp/pX%fwm)*X_carbonate_melt)/(1.0_dp+(1.0_dp-(44.01_dp/pX%fwm))* X_carbonate_melt) ! with 44.01 represent the MM of CO2 (12+2*16)
Kf = 10.0_dp**(KI+KII+fO2)
X_carbonate_melt = Kf/(1.0_dp+Kf)
part%X_CO2(m)=((44.01_dp/pX%fwm)*X_carbonate_melt)/(1.0_dp-(1.0_dp-(44.01_dp/pX%fwm))* X_carbonate_melt) ! with 44.01 represent the MM of CO2 (12+2*16)
part%X_CO2(m) = (part%X_CO2(m))*(1.0_dp-(1.0_dp-part%df(m))**(1.0_dp/pX%dCO2)) !/part%df(m) (df already in X_CO2 now) ! fractional melting
part%X_CO2(m) = part%X_CO2(m) * (1.0_dp-(1.0_dp-part%df(m))**(1.0_dp/pX%dCO2)) !/part%df(m) (df already in X_CO2 now) ! fractional melting
part%X_CO2(m) = part%X_CO2(m) * 10.0_dp**6
! if (part%X_CO2(m)*part%df(m).gt.part%CO2(m)) then
! part%X_CO2(m) = part%CO2(m)/part%df(m)
if (part%X_CO2(m).gt.part%CO2(m)) then
part%X_CO2(m) = part%CO2(m)
end if
part%CO2(m) = part%CO2(m) - part%X_CO2(m)!*part%df(m)
part%CO2(m) = part%CO2(m) - part%X_CO2(m)
end if
end do
......@@ -269,10 +270,12 @@ contains
! for cells !
!!!!!!!!!!!!!
!D=((4*pi__)/3)*((mesh%rmax**(3)-mesh%rmin**(3))/((mesh%rmax**(2)-mesh%rmin**(2)))) !D is for the shift from 2D to 3D calculations
! scale from regional 2D sphere to full 2D sphere
D=1.0_dp/pG%phi_factor
! calc for 2D sphere -> / (pi * (rmax^2-rmin^2))
! scale to 3D volume -> * (4pi/3 * (rmax^3-rmax^2))
D=4.0_dp/3.0_dp * (mesh%rmax**(3)-mesh%rmin**(3)) / (mesh%rmax**(2)-mesh%rmin**(2)) !D is for the shift from 2D to 3D calculations
! scale to 3D volume -> * (4pi/3 * (rmax^3-rmax^3))
D=D*4.0_dp/3.0_dp * (mesh%rmax**(3)-mesh%rmin**(3)) / (mesh%rmax**(2)-mesh%rmin**(2)) !D is for the shift from 2D to 3D calculations
mc%mass_H2O=0.0_dp
mc%mass_CO2=0.0_dp
......@@ -296,12 +299,6 @@ contains
mc%Concentration_CO2_frac = 0.0_dp
mc%max_CO2_frac= 0.0_dp
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! ! Test depletion via cells, not particles !
! cell%X_CO2 = 0.0_dp
! cell%X_H2O = 0.0_dp
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
do i=0,mesh%nl+1
do j=jmin,jmax
do k=0,mesh%nr+1
......@@ -316,40 +313,6 @@ contains
Tm_surf = T_K*exp(-pX%alpha_melt*pF%alpha/(pX%density_melt*pF%rho*pX%Cp_melt*pF%Cp)*(P_bar-Psurf)*10.0**5)
!write(*,*) i,j,k,"T=",T_K,",Tm_surf=",Tm_surf,", P=",P_bar
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! ! Test: overwrite particle values below, used fo degassing cell-averaged values, only !
!
! cell%X_H2O(i,j,k) = cell%w(i,j,k)/( pX%dH2O + mc%DeltaF(i,j,k)*(1.0_dp-pX%dH2O) ) ! cell%X_H2O / part%X_H2O is concentration in melt
!! if (part%X_H2O(m)*part%df(m).gt.part%w(m)) then ! condition for not extracting more water than available
!! 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)),cell%w(i,j,k)
!
! cell%w(i,j,k) = cell%w(i,j,k) - cell%X_H2O(i,j,k)*mc%DeltaF(i,j,k)
!
! ! CO2 equations (Holloway)
!
! KI= 40.07639_dp-((2.53932_dp*10.0_dp**(-2))*T_K)+(5.27096_dp*10.0_dp**(-6))*(T_K**2)+0.0267_dp*(P_bar-1.0_dp)/T_K
! KII= -6.24763_dp-(282.56_dp/T_K)-(0.119242_dp*(P_bar-1000.0_dp)/T_K)
! KI=10.0_dp**(KI) ! for the inverse of LOG10
! KII= 10.0_dp**(KII)
!
! fO2= 6.899_dp-(27714.0_dp/T_K)+0.05_dp*(P_bar-1.0_dp)/T_K + pX%fO2 ! instead fO2-field, where with change in O atoms we calculate with the same formula back how fO2 changes locally
! fO2= 10.0_dp**(fO2)
! X_carbonate_melt = (KI*KII*fO2)/(1.0_dp+(KI*KII*fO2))
! cell%X_CO2(i,j,k)=((44.01_dp/pX%fwm)*X_carbonate_melt)/(1.0_dp+(1.0_dp-(44.01_dp/pX%fwm))* X_carbonate_melt) ! with 44.01 represent the MM of CO2 (12+2*16)
! cell%X_CO2(i,j,k) = (cell%X_CO2(i,j,k))*(1.0_dp-(1.0_dp-mc%DeltaF(i,j,k))**(1.0_dp/pX%dCO2))/mc%DeltaF(i,j,k) ! fractional melting
! cell%X_CO2(i,j,k) = cell%X_CO2(i,j,k) * 10.0_dp**6
!! if (part%X_CO2(m)*part%df(m).gt.part%CO2(m)) then
!! part%X_CO2(m) = part%CO2(m)/part%df(m)
!! end if
! cell%CO2(i,j,k) = cell%CO2(i,j,k) - cell%X_CO2(i,j,k)*mc%DeltaF(i,j,k)
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! go from wt-ppm of H2O/CO2 in melt to mole fractions (ppm does not matter here since only ratio counts)
! 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) ! X_H2O and X_CO2 are amount of volatiles leaving rock and partitioning into melt
......@@ -375,35 +338,25 @@ contains
if (mole_frac_H2O.gt.0.0_dp) then
mc%mass_H2O = mc%mass_H2O + (cell%X_H2O(i,j,k)/M_H2O * frac_H2O/mole_frac_H2O * M_H2O)&
& * vol_cell*cell%ref_r(i,j,k)*(1.0e-6_dp) ! no *DeltaF needed here, since already multiplied above for particles
! & * mc%DeltaF(i,j,k)*vol_cell*cell%ref_r(i,j,k)*(1.0e-6_dp)
mc%mass_H2 = mc%mass_H2 + (cell%X_H2O(i,j,k)/M_H2O * frac_H2/mole_frac_H2O * M_H2)&
& * vol_cell*cell%ref_r(i,j,k)*(1.0e-6_dp)
! & * mc%DeltaF(i,j,k)*vol_cell*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 * frac_CO2/mole_frac_CO2 * M_CO2)&
& * vol_cell*cell%ref_r(i,j,k)*(1.0e-6_dp)
! & * mc%DeltaF(i,j,k)*vol_cell*cell%ref_r(i,j,k)*(1.0e-6_dp)
mc%mass_CO = mc%mass_CO + (cell%X_CO2(i,j,k)/M_CO2 * frac_CO/mole_frac_CO2 * M_CO)&
& * vol_cell*cell%ref_r(i,j,k)*(1.0e-6_dp)
! & * mc%DeltaF(i,j,k)*vol_cell*cell%ref_r(i,j,k)*(1.0e-6_dp)
end if
else
! calc mass H2O/CO2 each in kg without speciation model
mc%mass_H2O = mc%mass_H2O+cell%X_H2O(i,j,k)*vol_cell& !*mc%DeltaF(i,j,k)*vol_cell&
mc%mass_H2O = mc%mass_H2O+cell%X_H2O(i,j,k)*vol_cell&
&*cell%ref_r(i,j,k)*(1.0e-6_dp) ! partial weight
mc%mass_CO2 = mc%mass_CO2+cell%X_CO2(i,j,k)*vol_cell& !*mc%DeltaF(i,j,k)*vol_cell&
mc%mass_CO2 = mc%mass_CO2+cell%X_CO2(i,j,k)*vol_cell&
&*cell%ref_r(i,j,k)*(1.0e-6_dp) ! partial weight
end if
M_melt = M_melt + mc%DeltaF(i,j,k)*vol_cell*cell%ref_r(i,j,k)
! trace CO2 in melt
! if (cell%X_CO2(i,j,k).gt.mc%max_CO2_frac) then
! mc%max_CO2_frac = cell%X_CO2(i,j,k)
! mc%Concentration_CO2_melt = mc%Concentration_CO2_melt + (cell%X_CO2(i,j,k)*mc%DeltaF(i,j,k) * mesh%dVi(i,j,k))
! mc%Concentration_CO2_frac = mc%Concentration_CO2_frac + ((cell%X_CO2(i,j,k)) * mesh%dVi(i,j,k))
! end if
if (cell%X_CO2(i,j,k)/mc%DeltaF(i,j,k).gt.mc%max_CO2_frac) then
mc%max_CO2_frac = cell%X_CO2(i,j,k)/mc%DeltaF(i,j,k)
mc%Concentration_CO2_melt = mc%Concentration_CO2_melt + cell%X_CO2(i,j,k) * mesh%dVi(i,j,k)
......@@ -415,24 +368,6 @@ contains
enddo
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! ! Test depletion via cells, not particles, reset particle values here !
!
! do m=1,part%n
! i = part%cell_l(m)
! j = part%cell_y(m)
! k = part%cell_r(m)
!
! part%X_H2O(m) = cell%X_H2O(i,j,k)
! part%X_CO2(m) = cell%X_CO2(i,j,k)
! part%w(m) = cell%w(i,j,k)
! part%CO2(m) = cell%CO2(i,j,k)
! part%df(m) = mc%DeltaF(i,j,k)
! part%d(m) = cell%d(i,j,k)
! enddo
!
! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
mc%mass_H2O = mc%mass_H2O * D * pM%Chi_extr ! in nondimensional Kg, scaled from 2D to 3D mantle, only extrusive melt here
mc%mass_CO2 = mc%mass_CO2 * D * pM%Chi_extr
mc%mass_H2 = mc%mass_H2 * D * pM%Chi_extr
......
......@@ -1183,7 +1183,7 @@ contains
case("n_perm")
pM%n_perm = int(val)
if (rank.eq.0) write(*,*) var,"=",pM%n_perm
case("Chi_CO2")
case("Chi_CO2") ! not used anymore -> instead CO2_ini as initial mantle value
pM%Chi_CO2 = val
if (rank.eq.0) write(*,*) var,"=",pM%Chi_CO2
case("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