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

Update initial profiles after MO; update melt temp in model 1D

parent 6be7933a
......@@ -357,6 +357,7 @@ T_surf_melt(1) = T_melt_um(1);%T_intr_lid(1)*exp(-alpha/(rho*Cp)*p_intr_lid(1));
%%% Init volatile arrays
deg_CO2_atm = zeros(size(crust));, deg_H2O_atm = zeros(size(crust));
deg_CO2_lid = zeros(size(crust));, deg_H2O_lid = zeros(size(crust));
deg_CO2_mid = zeros(size(crust));, deg_H2O_mid = zeros(size(crust));
deg_CO2_cru = zeros(size(crust));, deg_H2O_cru = zeros(size(crust));
......@@ -368,6 +369,7 @@ H2_mole_lid = zeros(size(crust));, H2O_mole_lid = zeros(size(crust));, CO_mole_l
H2_mole_mid = zeros(size(crust));, H2O_mole_mid = zeros(size(crust));, CO_mole_mid = zeros(size(crust));, CO2_mole_mid = zeros(size(crust));
H2_mole_cru = zeros(size(crust));, H2O_mole_cru = zeros(size(crust));, CO_mole_cru = zeros(size(crust));, CO2_mole_cru = zeros(size(crust));
X_H2O_mole = zeros(size(crust));, X_CO2_mole = zeros(size(crust));
X_H2O_mole_atm = zeros(size(crust));, X_CO2_mole_atm = zeros(size(crust));
X_H2O_mole_lid = zeros(size(crust));, X_CO2_mole_lid = zeros(size(crust));
X_H2O_mole_mid = zeros(size(crust));, X_CO2_mole_mid = zeros(size(crust));
X_H2O_mole_cru = zeros(size(crust));, X_CO2_mole_cru = zeros(size(crust));
......@@ -522,12 +524,29 @@ for i=1:nt
Count_melt_av = Count_melt_av + melt_vol(i);
end
% mole amount of H2O and CO2 for extrusive outgassing
% mole amount of H2O and CO2 in melt
X_H2O_mole(i) = C_H2O_melt(i)/mass_H2O / (C_H2O_melt(i)/mass_H2O+C_CO2_melt(i)/mass_CO2); % mole amount of X_H2O and X_CO2
X_CO2_mole(i) = C_CO2_melt(i)/mass_CO2 / (C_H2O_melt(i)/mass_H2O+C_CO2_melt(i)/mass_CO2); % sum of both is 1
[H2_mole(i), H2O_mole(i), CO_mole(i), CO2_mole(i)]=ortenzi_simultaneous_eq_C_O_H(T_surf_melt(i), 1, X_H2O_mole(i), X_CO2_mole(i), redox_buffer);
% extrusive outgassing
% Old version without surface solubility depending on P_H2O/P_CO2; reproduced results in Ortenzi et al., submitted, first version:
[H2_mole(i), H2O_mole(i), CO_mole(i), CO2_mole(i)]=ortenzi_simultaneous_eq_C_O_H(T_surf_melt(i), 1, X_H2O_mole(i), X_CO2_mole(i), redox_buffer);
% New version for revision of Ortenzi et al.:
%[deg_CO2_atm(i),deg_H2O_atm(i)] = get_sol_atm(X_H2O_mole(i),X_CO2_mole(i),1,T_surf_melt(i),0,0);%p_surf(i)*1e-5,T_surf_melt(i),P_CO2(i),P_H2O(i)); %still needs to be adapted to time-dependent pressure
%if (deg_H2O_atm(i)+deg_CO2_atm(i)==0)
% X_H2O_mole_atm(i) = 0;, X_CO2_mole_atm(i) = 0;, H2_mole(i) = 0;, H2O_mole(i) = 0; CO_mole(i) = 0;, CO2_mole(i) = 0;
%else
% X_H2O_mole_atm(i) = deg_H2O_atm(i)/mass_H2O / (deg_H2O_atm(i)/mass_H2O+deg_CO2_atm(i)/mass_CO2); % mole amount of X_H2O and X_CO2
% X_CO2_mole_atm(i) = deg_CO2_atm(i)/mass_CO2 / (deg_H2O_atm(i)/mass_H2O+deg_CO2_atm(i)/mass_CO2); % sum of both is 1
% [H2_mole(i), H2O_mole(i), CO_mole(i), CO2_mole(i)]=ortenzi_simultaneous_eq_C_O_H(T_surf_melt(i), 1, X_H2O_mole_atm(i), X_CO2_mole_atm(i), redox_buffer); % ToDo: psurf(i)
%end
% surface pressures change over time; needs different implementation than
% so far used -> stored in new file melt_output_for_speciation_lid_sol.m!
if (Chi_intr_lid>0)
% intrusive below lithosphere
[deg_CO2_lid(i),deg_H2O_lid(i)] = get_sol(X_H2O_mole(i),X_CO2_mole(i),p_intr_lid(i)*1e-5,T_intr_lid(i));
......@@ -574,10 +593,6 @@ end
%fprintf('At end: Crust is %f km thick, lithosphere is %f km thick\n',crust(nt)/1000,p_intr_lid(nt)/rho/g/1000)
%%%%%%%%%%%%%%%%%%%%%%
%ToDo: add linear interpolation case between different redox states (e.g. QFM->IW or IW->QFM)
%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%
% extrusive volcanism
%%%%%%%%%%%%%%%%%%%%%%
......@@ -835,7 +850,8 @@ dr_r1=dr_r(nt);
%path_write = 'C:\OwnCloud_FUB\Literature_Projects\Eigene_Projekte\Project_Ortenzi_Dorn\OctaveResultsNew2\';
%path_write = '/home/noacklen/ownCloud/Literature_Projects/Eigene_Projekte/Project_Ortenzi_Dorn/OctaveResultsNew3/';
path_write = '/scratch-2/PublishedProjects/Project_CaroDorn/Results_New/';
%path_write = '/scratch-2/PublishedProjects/Project_CaroDorn/Results_New/';
path_write = '/scratch-2/PublishedProjects/Project_CaroDorn/Results_Revision/';
%path_write = sprintf('/home/noacklen/ownCloud/Literature_Projects/Eigene_Projekte/Project_Ortenzi_Dorn/OctaveResultsNew/%s/',redox);
%path_write = sprintf('/home/noacklen/Arbeit/Ortenzi_Dorn_Paper/%s/',redox); % to avoid update conflicts with NextCloud
......@@ -857,7 +873,7 @@ new_name = sprintf('Case%s--%s--%s--%s--%dCO2--%dH2O.txt',name2,name3,name4,redo
cd(path_write);
fid=fopen(new_name,'w');
fprintf(fid,['time(Gyr) p_int_lid(bar) p_int_cru p_int_mid d_cr(km) d_lith(km) P_H2(bar) P_H2O P_CO P_CO2 mu_av M_tot_acc(kg) ', ...
'frac_H2_atm frac_H2O_atm frac_CO_atm frac_CO2_atm M_H2_in_l(kg) M_H2O_in_l M_CO_in_l M_CO2_in_l M_H2_in_c M_H2O_in_c M_CO_in_c M_CO2_in_c ',...
......@@ -1011,7 +1027,9 @@ function retval = startsWith(s, prefix)
retval = strncmp(s, prefix, n);
end
function [deg_CO2,deg_H2O] = get_sol(X_H2O_mantle,X_CO2_mantle,p_lid,T_lid);
% solubility in gas bubbles in lithosphere/crust - depends on lithostatic
% pressure but not any atmosphere partial pressures
function [deg_CO2,deg_H2O] = get_sol(X_H2O_mantle,X_CO2_mantle,p_lid,T_lid)
% solubility calculated from Iacono-Marzian et al., 2012; H2O in wt-% and CO2 in wt-ppm
X_H2O = X_H2O_mantle*10^-4; %ppm->weight percent of water in the melt
X_CO2 = X_CO2_mantle; %weight ppm of CO2 in the melt
......@@ -1075,3 +1093,28 @@ deg_CO2 = dX_CO2; % wt-ppm
deg_H2O = dX_H2O*10^4; % wt-% -> wt-ppm
end
% solubility at surface; depends on total atmospheric pressure but also
% P_H2O and P_CO2 partial pressures (H2 and CO assumed to have solubility % of 0)
function [deg_CO2,deg_H2O] = get_sol_atm(X_H2O_mantle,X_CO2_mantle,p_lid,T_lid,P_CO2,P_H2O)
% solubility calculated from Iacono-Marzian et al., 2012; H2O in wt-% and CO2 in wt-ppm
X_H2O = X_H2O_mantle*10^-4; %ppm->weight percent of water in the melt
X_CO2 = X_CO2_mantle; %weight ppm of CO2 in the melt
if (X_H2O<4) % personal communication with Fabrice Gaillard
hydrous=0;
else
hydrous=1;
end
H2O_fl_mole = P_H2O/p_lid; % H2O fraction in atmosphere
CO2_fl_mole = P_CO2/p_lid; % CO2 fraction in atmosphere
[H2O_dissolved,CO2_dissolved] = solubility(p_lid,T_lid,H2O_fl_mole,X_H2O,hydrous,CO2_fl_mole);
dX_H2O = max(0,X_H2O-H2O_dissolved); % what would go into gas bubble in wt-%
dX_CO2 = max(0,X_CO2-CO2_dissolved); % what would go into gas bubble in wt-ppm
deg_CO2 = dX_CO2; % wt-ppm
deg_H2O = dX_H2O*10^4; % wt-% -> wt-ppm
end
......@@ -4,6 +4,7 @@
function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
% if reduc=1, then plot only T2D, V2D, Vi2D, D2D, Mradp2D
addpath("/home/noacklen/Arbeit/SOURCE/CHIC/Octave")
[data,input,program_dir,fname,fpath] = read_snap(fin); % read data from snapshot
% cd(program_dir) % return to octave folder
......@@ -166,6 +167,10 @@ end
% plot min/mean/max profiles for field x
function x_profile = plot_x_profile(field,data,input,name,title_plot,orient,dimen,log=0,melt=0)
% ATTENTION: melt profiles and pressure profile not correct, since mantle-averaged density is used
% INSTEAD: should read in profs.res and take pressure profile from there
% (not for melt temp though, since calculated differently in convection simulationes) !!!
nl=data.nl; % for profiles independent of periodic or reflective boundary conditions
ny=data.ny;
nr=data.nr;
......@@ -200,6 +205,8 @@ function x_profile = plot_x_profile(field,data,input,name,title_plot,orient,dime
until k>nr
if (log==1), x_profile = log10(x_profile);, endif
% problem: rho is mantle-averaged, hence if strongly compressible mantle, then this shifts buoyant melt radius slightly
f = figure('visible','off');
......@@ -211,7 +218,6 @@ function x_profile = plot_x_profile(field,data,input,name,title_plot,orient,dime
h4=plot(m_profile(:,1),r,":g");
h5=plot(m_profile(:,2),r,"--g");
endif
hold off
set(h,'LineWidth',5)
set(h2,'LineWidth',10)
set(h3,'LineWidth',5)
......@@ -219,6 +225,14 @@ function x_profile = plot_x_profile(field,data,input,name,title_plot,orient,dime
set(h4,'LineWidth',5)
set(h5,'LineWidth',5)
endif
if (melt==1)
radius_melt_buoyant = max(r)-input.md*10.0^6/(input.rho*input.g);
plot([min(x_profile(:,1)),max(x_profile(:,3))],[radius_melt_buoyant,radius_melt_buoyant],'k:')
axis([min(x_profile(:,1)),3000,max(r)-3*input.md*10.0^6/(input.rho*input.g),max(r)])
endif
hold off
% legend(sprintf(strcat(name,' min')),sprintf(strcat(name," mean")),sprintf(strcat(name," max")), "location", orient);
% legend("boxoff");
......
......@@ -64,6 +64,7 @@ function input = read_input(fpath)
input.m1l2 = value(strcmp("melt1l_2", data));
input.m1l3 = value(strcmp("melt1l_3", data));
input.mz = value(strcmp("melt_z_change", data)); % GPa; in code non-dim with init rho,g to get depth, better to compare with actual pressure?
input.md = value(strcmp("melt_depth", data)); % until which pressure melt is buoyant
input.m2s0 = value(strcmp("melt2s_0", data));
input.m2s1 = value(strcmp("melt2s_1", data));
input.m2s2 = value(strcmp("melt2s_2", data));
......@@ -90,5 +91,58 @@ function input = read_input(fpath)
if (isempty(input.set_tracers)), input.set_tracers=1;, endif
input.DT = input.Tbot-input.Ttop;
if (input.read_profs>0) % read-in profiles and use mantle-averaged values for density, g etc.:
% take non-dim values from screenlog.0 if used:
fid = fopen('screenlog.0');
tline = fgetl(fid);
while ischar(tline)
if startsWith(tline,' pF%Cp=')
[Cpline,remain] = strtok(tline,',');
input.Cp = str2num(Cpline(8:end));
[kline,remain] = strtok(remain,',');
input.k = str2num(kline(7:end));
[gline,remain] = strtok(remain,',');
input.g = str2num(gline(7:end));
[rholine,remain] = strtok(remain,',');
input.rho = str2num(rholine(9:end));
[Dline,aline] = strtok(remain,',');
input.D = str2num(Dline(7:end));
input.alpha = str2num(aline(12:end));
end
tline = fgetl(fid);
end
fclose(fid);
cd(program_dir)
\ No newline at end of file
input.Cp,input.k,input.g,input.rho,input.D,input.alpha
end
% read-in profile if profs>=5
if (input.read_profs>0)
%fid = fopen('profs.res');
%data=textread('profs.res','%f');
%data = csvread('profs.res');
%data = fread(fid);
%fclose(fid);
% data = reshape(date,1000,size(data)/1000)
data = dlmread('profs.res');
size(data)
% data(1,:)
end
cd(program_dir)
end
% Check if a string starts with a given prefix
% Returns 1 if s starts with prefix, 0 else
function retval = startsWith(s, prefix)
n = length(prefix);
if n == 0 % Empty prefix
retval = 1; % Every string starts with empty prefix
return
end
retval = strncmp(s, prefix, n);
end
\ No newline at end of file
% solubility calculated from Iacono-Marzian et al., 2012
% valied only for T in 900-1500 °C and p in 0 - 10 000 bar (1 GPa)
function [H2O_dissolved,CO2_dissolved] = solubility(P_bar,T_K,H2O_fl_mole,X_H2O_wt,hydrous) % X_H2O here in wt-fraction, not molar!!!
function [H2O_dissolved,CO2_dissolved] = solubility(P_bar,T_K,H2O_fl_mole,X_H2O_wt,hydrous,CO2_fl_mole) % X_H2O here in wt-fraction, not molar!!!
if (nargin<7)
CO2_fl_mole = 1.0-H2O_fl_mole; % old version of solubility did not use CO2_fl_mole; used only H2O-CO2-pressures; now additional species possible
end
% molar masses, g/mol
M_H = 1.0079;
......@@ -194,7 +198,7 @@ else
end
% partial pressures from H2O_fl_mole (water mole fraction); P_i = mole_i*m_av*g/A
P_CO2 = (1.0-H2O_fl_mole)*P_bar;
P_CO2 = CO2_fl_mole*P_bar; %(1.0-H2O_fl_mole)*P_bar;
P_H2O = H2O_fl_mole*P_bar;
ln_H2O = a_H2O*log(P_H2O) + b_H2O*NBO_O + B_H2O + C_H2O*P_bar/T_K; % in wt-%
......@@ -203,7 +207,8 @@ X_H2O_mol = exp(ln_H2O)/100/Sum_mol/(2*M_H+M_O); % in wt-amount
% in formula here Fe2O3 missing since not in original paper
%if hydrous==1
%ln_CO_3_2min = X_H2O*d_H2O + X_AI*d_AI + (X_FeO+X_MgO)*d_FeMg + (X_Na2O+X_K2O)*d_NaK ... % in ppm
ln_CO_3_2min = X_H2O_mol*d_H2O + X_AI*d_AI + (X_FeO+X_MgO)*d_FeMg + (X_Na2O+X_K2O)*d_NaK ... % in ppm
%ln_CO_3_2min = X_H2O_mol*d_H2O + X_AI*d_AI + (X_FeO+X_MgO)*d_FeMg + (X_Na2O+X_K2O)*d_NaK ... % in ppm
ln_CO_3_2min = min(X_H2O,X_H2O_mol)*d_H2O + X_AI*d_AI + (X_FeO+X_MgO)*d_FeMg + (X_Na2O+X_K2O)*d_NaK ... % in ppm
+ a_CO2*log(P_CO2) + b_CO2*NBO_O + B_CO2 + C_CO2*P_bar/T_K;
%else
%ln_CO_3_2min = X_AI*d_AI + (X_FeO+X_MgO)*d_FeMg + (X_Na2O+X_K2O)*d_NaK ... % in ppm
......@@ -214,7 +219,7 @@ H2O_dissolved = exp(ln_H2O); % in wt-%
if (H2O_fl_mole==0), H2O_dissolved = 0;, end
CO32m_dissolved = exp(ln_CO_3_2min); % in wt-ppm
CO2_dissolved = CO32m_dissolved/M_carb*M_CO2; % in wt-ppm
if (H2O_fl_mole==1), CO2_dissolved = 0;, end
if ((H2O_fl_mole==1) || (CO2_fl_mole==0)), CO2_dissolved = 0;, end
%fwm=36.594; % KILAUEA THOLEIITIC BASALT, HOLLOWAY 1992
%%CO2_dissolved = CO32m_dissolved*1.202656 / (1+CO32m_dissolved*0.202656)
......@@ -231,4 +236,4 @@ if (H2O_fl_mole==1), CO2_dissolved = 0;, end
%
%%X_CO23m_melt = X_CO23m_melt - CO32m_dissolved
end %function
\ No newline at end of file
end %function
......@@ -71,6 +71,9 @@ contains
! Psurf (in bar) in code has to be non-dimensional:
pE%Psurf = pE%Psurf / (pF%g*pF%rho*pF%D*(10.0_dp**(-5)))
mc%escape_H2O = 0.0_dp
mc%escape_CO2 = 0.0_dp
end subroutine
......@@ -135,7 +138,7 @@ contains
! & * 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_w_loss = M_w_loss + (pI%mantle_w_H2O-cell%w(i,j,k))*mesh%dVi(i,j,k) !actually needs to take into account cell%ref_r, since wt-ppm are used
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
......@@ -343,7 +346,8 @@ 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)
if (pX%solubility.eq.1) then
call solubility_calc(mc%P_CO2,mc%P_H2O,Tm_surf,part%X_H2O(m)/part%df(m),part%X_CO2(m)/part%df(m),&
call solubility_calc(mc%P_N2+mc%P_CO2+mc%P_H2O+mc%P_CO2+mc%P_H2O,& ! psurf total
&mc%P_CO2,mc%P_H2O,Tm_surf,part%X_H2O(m)/part%df(m),part%X_CO2(m)/part%df(m),&
&H2O_diss_melt,CO2_diss_melt,H2O_gas,CO2_gas)
!write(*,*) m,part%df(m),part%X_H2O(m),part%X_CO2(m),H2O_diss_melt,CO2_diss_melt,H2O_gas,CO2_gas
......@@ -579,14 +583,14 @@ contains
end subroutine atmos_escaping
subroutine solubility_calc(p_CO2,p_H2O,T_volc,H2O_melt,CO2_melt,H2O_diss_melt,CO2_diss_melt,H2O_gas,CO2_gas)
subroutine solubility_calc(p_atm,p_CO2,p_H2O,T_volc,H2O_melt,CO2_melt,H2O_diss_melt,CO2_diss_melt,H2O_gas,CO2_gas)
! input: pressures of already existing atmosphere in bar
! lava temperature close to solidus temperature in K
! H2O and CO2 in ppm in melt
! output: ppm of H2O and CO2 available for degassing (not dissolved in melt)
real(dp),intent(in)::p_CO2,p_H2O,T_volc,H2O_melt,CO2_melt
real(dp),intent(in)::p_atm,p_CO2,p_H2O,T_volc,H2O_melt,CO2_melt
real(dp),intent(out)::H2O_diss_melt,CO2_diss_melt,H2O_gas,CO2_gas
real(dp)::H2O_wt,H2O_fl_mole
real(dp)::H2O_wt
integer::hydrous
real(dp)::H2O_dissolved,CO2_dissolved
......@@ -598,8 +602,7 @@ contains
hydrous=1
endif
H2O_fl_mole = p_H2O/(p_CO2+p_H2O)
call solubility(p_CO2+p_H2O,T_volc,H2O_fl_mole,H2O_wt,hydrous,H2O_dissolved,CO2_dissolved) ! return values H2O in wt-% and CO2 in wt-ppm
call solubility(p_atm,T_volc,p_H2O/p_atm,p_CO2/p_atm,H2O_wt,hydrous,H2O_dissolved,CO2_dissolved) ! return values H2O in wt-% and CO2 in wt-ppm
H2O_dissolved = H2O_dissolved*10000.0_dp ! back to wt-ppm
H2O_diss_melt = min(H2O_dissolved,H2O_melt)
......@@ -613,13 +616,13 @@ contains
! solubility calculated from Iacono-Marzian et al., 2012
! valied only for T in 900-1500 °C and p in 0 - 10 000 bar (1 GPa)
subroutine solubility(P_bar,T_K,H2O_fl_mole,X_H2O_wt,hydrous,H2O_dissolved,CO2_dissolved) ! here P_bar only P_CO2+P_H2O, neglect other species
real(dp),intent(in)::P_bar,T_K,H2O_fl_mole,X_H2O_wt
subroutine solubility(P_bar,T_K,H2O_fl_mole,CO2_fl_mole,X_H2O_wt,hydrous,H2O_dissolved,CO2_dissolved) ! here P_bar only P_CO2+P_H2O, neglect other species
real(dp),intent(in)::P_bar,T_K,H2O_fl_mole,CO2_fl_mole,X_H2O_wt
real(dp),intent(out)::H2O_dissolved,CO2_dissolved
real(dp)::M_H,M_C,M_O,M_Na,M_Mg,M_Al,M_Si,M_K,M_Ca,M_Ti,M_Fe,M_CO2,M_H2O,M_carb
real(dp)::X_H2O,X_K2O,X_Na2O,X_CaO,X_MgO,X_FeO,X_Al2O3,X_SiO2,X_TiO2
real(dp)::M_melt,X_AI,Sum_mol,NBO,NBO_O
real(dp)::M_melt,X_AI,Sum_mol,NBO,NBO_O,X_H2O_mol
real(dp)::d_H2O,d_AI,d_FeMg,d_NaK,a_CO2,b_CO2,C__CO2,B__CO2,a_H2O,b_H2O,B__H2O,C__H2O
real(dp)::P_CO2,P_H2O,ln_H2O,ln_CO_3_2min,CO32m_dissolved,X_CO23m_melt,fwm,Mfac
integer,intent(in)::hydrous
......@@ -758,24 +761,26 @@ contains
! partial pressures from H2O_fl_mole (water mole fraction); P_i = mole_i*m_av*g/A
P_CO2 = (1.0-H2O_fl_mole)*P_bar
P_CO2 = CO2_fl_mole*P_bar
P_H2O = H2O_fl_mole*P_bar
ln_H2O = a_H2O*log(P_H2O) + b_H2O*NBO_O + B__H2O + C__H2O*P_bar/T_K ! in wt-%
X_H2O_mol = exp(ln_H2O)/100/Sum_mol/(2*M_H+M_O) ! in mole-amount
if (hydrous.eq.1) then
ln_CO_3_2min = X_H2O*d_H2O + X_AI*d_AI + (X_FeO+X_MgO)*d_FeMg + (X_Na2O+X_K2O)*d_NaK & ! in ppm
! should be same formulation for hydrous and anhydrous, since d_H2O known in both cases
!if (hydrous.eq.1) then
ln_CO_3_2min = min(X_H2O,X_H2O_mol)*d_H2O + X_AI*d_AI + (X_FeO+X_MgO)*d_FeMg + (X_Na2O+X_K2O)*d_NaK & ! in ppm
& + a_CO2*log(P_CO2) + b_CO2*NBO_O + B__CO2 + C__CO2*P_bar/T_K
else
ln_CO_3_2min = X_AI*d_AI + (X_FeO+X_MgO)*d_FeMg + (X_Na2O+X_K2O)*d_NaK & ! in ppm
& + a_CO2*log(P_CO2) + b_CO2*NBO_O + B__CO2 + C__CO2*P_bar/T_K
endif
!else
! ln_CO_3_2min = X_AI*d_AI + (X_FeO+X_MgO)*d_FeMg + (X_Na2O+X_K2O)*d_NaK & ! in ppm
! & + a_CO2*log(P_CO2) + b_CO2*NBO_O + B__CO2 + C__CO2*P_bar/T_K
!endif
H2O_dissolved = exp(ln_H2O) ! in wt-%
if (H2O_fl_mole.eq.0) H2O_dissolved = 0.0_dp
CO32m_dissolved = exp(ln_CO_3_2min) ! in wt-ppm
CO2_dissolved = CO32m_dissolved/M_carb*M_CO2 ! in wt-ppm
if (H2O_fl_mole.eq.1) CO2_dissolved = 0.0_dp
if (CO2_fl_mole.eq.0) CO2_dissolved = 0.0_dp
!fwm=36.594; % KILAUEA THOLEIITIC BASALT, HOLLOWAY 1992
!!CO2_dissolved = CO32m_dissolved*1.202656 / (1+CO32m_dissolved*0.202656)
......
......@@ -51,10 +51,11 @@ module initialize
end subroutine initialize_pressure_strR
subroutine initialize_cell(cell,field,pE,pF,pI,pT,pV,pX,rho_av,max_depl,rc,compress,&
& Di,DiT,Gr,T0,rmin,TinclTref,Tini,Tbot,shear,DeltaVisc,DeltaTc)
subroutine initialize_cell(cell,field,mesh,pE,pF,pI,pT,pV,pX,rho_av,max_depl,rc,compress,&
& Di,DiT,Gr,T0,rmin,TinclTref,Tini,Tbot,shear,DeltaVisc,DeltaTc,MOprofs)
type(cell_properties),intent(inout)::cell
type(variables_unknowns),intent(inout)::field
type(mesh_cp),intent(in)::mesh
type(param_E),intent(in)::pE
type(param_F),intent(in)::pF
type(param_I),intent(inout)::pI
......@@ -64,7 +65,7 @@ module initialize
real(dp),intent(inout)::rho_av ! for TidHeat>0
real(dp),intent(in)::max_depl,compress,Di,DiT,Gr,T0,rmin
real(dp),intent(inout)::shear ! set average shear for TidHeat
integer,intent(in)::TinclTref
integer,intent(in)::TinclTref,MOprofs
real(dp),intent(in)::rc(0:)
real(dp),intent(inout)::Tini,Tbot ! adapt Tbot depending on initial upper mantle temperature and read-in adiabatic profile
real(dp),intent(in)::DeltaVisc,DeltaTc
......@@ -73,7 +74,7 @@ module initialize
real(dp)::z,Cp_count,Core_count,vol,T_CMB
real(dp),allocatable::vals_t(:,:),vals_d(:,:),k_e(:)
integer,allocatable::mat(:)
real(dp)::gc,pi_,dr,g,rho,rb
real(dp)::gc,pi_,dr,g,rho,rb,Vav,check
pi_ = acos(-1.0_dp)
......@@ -90,9 +91,10 @@ module initialize
endif
cell%c = pI%rho ! composition density (e.g. by melt)
cell%w = pI%mantle_w_H2O
!cell%w = pI%mantle_w_H2O
cell%wr = pI%mantle_w_r ! water exponent for viscosity
cell%CO2 = pX%CO2_ini
cell%X_H2O = 0.0_dp
cell%X_CO2 = 0.0_dp
cell%fa = pI%mantle_fr_angl ! friction angle
cell%fc = pI%mantle_fr_coh ! cohesion
......@@ -111,6 +113,29 @@ module initialize
cell%shear = shear
cell%ref_p = 0.0_dp
if (MOprofs.eq.0) then
cell%w = pI%mantle_w_H2O
cell%CO2 = pX%CO2_ini
else
! magma ocean should be enriched in upper layers at end of magma ocean and depleted at CMB
! works only for spherical annulus, no for box or cylinder:
Vav = dot_product(mesh%rc(1:nr),mesh%rc(1:nr))/real(nr) !(mesh%rc(1)**2+mesh%rc(nr)**2)/2.0_dp ! needed to scale incompatible elements per volume
do k=0,nr+1
! simple version: volume of water linear from surface to CMB
! ToDo: actually should be divided by mesh%rc(k)**2 * cell%ref_r -> but ref_r only defined later; leave for now...
! maybe do this in main routine?
cell%w(:,:,k) = (mesh%rc(k)-mesh%rmin)*2.0_dp*pI%mantle_w_H2O * Vav/mesh%rc(k)**2 ! hence sum_k cell%w(k)*V(k) = m_w_H2O*Vm
cell%CO2(:,:,k) = (mesh%rc(k)-mesh%rmin)*2.0_dp*pX%CO2_ini * Vav/mesh%rc(k)**2
enddo
! test if works:
!write(*,*) "homogeneous:",pI%mantle_w_H2O*sum(mesh%dVi(1,1,1:nr)),&
! &",heteorogeneous:",dot_product(cell%w(1,1,1:nr),mesh%dVi(1,1,1:nr))
check = dot_product(mesh%dVi(1,1,1:nr),cell%w(1,1,1:nr))/(sum(mesh%dVi(1,1,1:nr))*pI%mantle_w_H2O)
cell%w = cell%w/check
check = dot_product(mesh%dVi(1,1,1:nr),cell%CO2(1,1,1:nr))/(sum(mesh%dVi(1,1,1:nr))*pX%CO2_ini)
cell%CO2 = cell%CO2/check
endif
if ((pI%k_max.ne.0.0_dp).and.(cell%fld_k)) then
do k=0,nr+1
......@@ -699,18 +724,22 @@ module initialize
end subroutine initialize_cell
subroutine initialize_heat(pC,pI,pF,pG,pX,H_from_surf,Crust_Lambda,cell)
subroutine initialize_heat(mesh,pC,pI,pF,pG,pX,H_from_surf,Crust_Lambda,cell,MOprofs)
type(mesh_cp),intent(in)::mesh
type(param_C),intent(in)::pC
type(param_I),intent(inout)::pI
type(param_F),intent(in)::pF
type(param_G),intent(in)::pG
type(param_X),intent(in)::pX
integer,intent(in)::H_from_surf
integer,intent(in)::H_from_surf,MOprofs
real(dp),intent(in)::Crust_Lambda
type(cell_properties),intent(inout)::cell
real(dp)::t ! 4.5 Gyr
real(dp)::Vc,Vm
real(dp)::Vc,Vm,Vav,check
integer::k,nr
nr = mesh%nr
! specific heating rates in W/kg
pI%H0U238 = 9.46e-5 * pF%D**2/(pF%DeltaT * pF%k) * pF%rho
......@@ -757,14 +786,52 @@ module initialize
if (pX%dU238.gt.0.0_dp) then
! initial mantle radioactive heat sources, crust values are these values multiplied with Crust_Lambda and are set in melt.f90 in init crust routine
! Volume actually 4/3 pi (R*D)³, but later volumes divided by each oder, so take just non-dim radii here vor valume comparison
! Volume actually 4/3 pi (R*D)³, but later volumes divided by each oder, so take just non-dim radii here for volume comparison
! ToDo: should it not be r^2 since 2D spherical annulus?
Vm = (pG%rmax-pC%crust_ini)**3 - pG%rmin**3
Vc = pG%rmax**3 - (pG%rmax-pC%crust_ini)**3
cell%X_U238 = pI%X0U238*(Vm+Vc)/(Vm+Crust_Lambda*Vc)
cell%X_U235 = pI%X0U235*(Vm+Vc)/(Vm+Crust_Lambda*Vc)
cell%X_Th232 = pI%X0Th232*(Vm+Vc)/(Vm+Crust_Lambda*Vc)
cell%X_K40 = pI%X0K40*(Vm+Vc)/(Vm+Crust_Lambda*Vc)
!cell%X_U238 = pI%X0U238*(Vm+Vc)/(Vm+Crust_Lambda*Vc)
!cell%X_U235 = pI%X0U235*(Vm+Vc)/(Vm+Crust_Lambda*Vc)
!cell%X_Th232 = pI%X0Th232*(Vm+Vc)/(Vm+Crust_Lambda*Vc)
!cell%X_K40 = pI%X0K40*(Vm+Vc)/(Vm+Crust_Lambda*Vc)
! overwrite pI%X0's with reformulated mantle value; set crut value based on pI%X0's; not needed anymore apart after that
pI%X0U238 = pI%X0U238*(Vm+Vc)/(Vm+Crust_Lambda*Vc)
pI%X0U235 = pI%X0U235*(Vm+Vc)/(Vm+Crust_Lambda*Vc)
pI%X0Th232 = pI%X0Th232*(Vm+Vc)/(Vm+Crust_Lambda*Vc)
pI%X0K40 = pI%X0K40*(Vm+Vc)/(Vm+Crust_Lambda*Vc)
if (MOprofs.eq.0) then
cell%X_U238 = pI%X0U238
cell%X_U235 = pI%X0U235
cell%X_Th232= pI%X0Th232
cell%X_K40 = pI%X0K40
else
! magma ocean should be enriched in upper layers at end of magma ocean and depleted at CMB
! works only for spherical annulus, no for box or cylinder:
Vav = dot_product(mesh%rc(1:mesh%nr),mesh%rc(1:mesh%nr))/real(mesh%nr)
!incompatible elements per volume
do k=0,mesh%nr+1
cell%X_U238(:,:,k) = (mesh%rc(k)-mesh%rmin)*2.0_dp*pI%X0U238 * Vav/(mesh%rc(k)**2*cell%ref_r(1,1,k))
cell%X_U235(:,:,k) = (mesh%rc(k)-mesh%rmin)*2.0_dp*pI%X0U235 * Vav/(mesh%rc(k)**2*cell%ref_r(1,1,k))
cell%X_Th232(:,:,k)= (mesh%rc(k)-mesh%rmin)*2.0_dp*pI%X0Th232* Vav/(mesh%rc(k)**2*cell%ref_r(1,1,k))
cell%X_K40(:,:,k) = (mesh%rc(k)-mesh%rmin)*2.0_dp*pI%X0K40 * Vav/(mesh%rc(k)**2*cell%ref_r(1,1,k))
!write(*,*) k,cell%X_U238(1,1,k),cell%X_U235(1,1,k),cell%X_Th232(1,1,k),cell%X_K40(1,1,k)
enddo
if (pC%crust_ini.eq.0.0_dp) then
! if not then rescaling is done in melt.f90
check = 0.0_dp
do k=1,mesh%nr
check = check + mesh%dVi(1,1,k)*cell%ref_r(1,1,k)*cell%X_U238(1,1,k) ! scaling factor, same for all radiogenic species
enddo
check = check/(sum(mesh%dVi(1,1,1:nr))*pI%X0U238)
cell%X_U238 = cell%X_U238/check
cell%X_U235 = cell%X_U235/check
cell%X_Th232= cell%X_Th232/check
cell%X_K40 = cell%X_K40/check
endif
endif
!write(*,*) "X_U238=",cell%X_U238(1,1,1),pI%X0U238,Vm,Vc,Crust_Lambda
endif
......@@ -781,6 +848,8 @@ module initialize
real(dp)::heat_t,heat_b,Q
real(dp),allocatable::vals_t(:,:),time_star(:)
pH%Q_eddy=0.0_dp
if (time.eq.0.0_dp) then
! initialize
if (pH%use_magn_H.eq.1.0_dp) then
......
......@@ -843,6 +843,9 @@ contains
case("ec_model")
pI%ec_model = int(val)
if (rank.eq.0) write(*,*) var,"=",pI%ec_model
case("MOprofs")
pI%MOprofs = int(val)
if (rank.eq.0) write(*,*) var,"=",pI%MOprofs
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! H: Heating parameters !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
......
......@@ -528,8 +528,8 @@ program CHIC
endif
! first call needed for adiabatic temperature profile
call initialize_cell(cell,field,pE,pF,pI,pT,pV,pX,pH%rho_av,pM%max_depl,mesh%rc,pI%compress,pI%Di,pI%DiT,pI%Gr,pT%T0,&
&pG%save_rmin,pI%TinclTref,pT%Tini,pT%Tbottom,pH%shear_modulus,pV%DeltaVisc,pT%DeltaTc)
call initialize_cell(cell,field,mesh,pE,pF,pI,pT,pV,pX,pH%rho_av,pM%max_depl,mesh%rc,pI%compress,pI%Di,pI%DiT,pI%Gr,pT%T0,&
&pG%save_rmin,pI%TinclTref,pT%Tini,pT%Tbottom,pH%shear_modulus,pV%DeltaVisc,pT%DeltaTc,pI%MOprofs)
call initialize_temperature(jmin_2D,field,mesh%rc,field%T,pT%Tbottom,pT%Ttop,pG%inner_bound,pG%outer_bound,&
& pT%linear,pG%save_rmin,pG%save_rmax,pT%Tini,pT%s_bot,pT%s_top,pG%geom,pG%dim,pI%DiT,LBOUND(cell%cp,3),&
......@@ -537,8 +537,8 @@ program CHIC
! second call for re-setting temperature if needed with profs.res
if (pI%read_profs.ge.6) then
call initialize_cell(cell,field,pE,pF,pI,pT,pV,pX,pH%rho_av,pM%max_depl,mesh%rc,pI%compress,pI%Di,pI%DiT,pI%Gr,pT%T0,&
&pG%save_rmin,pI%TinclTref,pT%Tini,pT%Tbottom,pH%shear_modulus,pV%DeltaVisc,pT%DeltaTc)
call initialize_cell(cell,field,mesh,pE,pF,pI,pT,pV,pX,pH%rho_av,pM%max_depl,mesh%rc,pI%compress,pI%Di,pI%DiT,pI%Gr,pT%T0,&
&pG%save_rmin,pI%TinclTref,pT%Tini,pT%Tbottom,pH%shear_modulus,pV%DeltaVisc,pT%DeltaTc,pI%MOprofs)