Commit 8df9f47f authored by Lena Noack's avatar Lena Noack
Browse files

Start update Octave code to use in Matlab

parent 91161f01
......@@ -28,13 +28,13 @@ function melt_output_for_speciation(fpath=' ',all=1,redox="IW")
end
function retval = rdir (d)
## info about D
% info about D
d_info = dir (d);
## list of directories in D, not including . and ..
% list of directories in D, not including . and ..
names = {d_info.name};
subdirs = names([d_info.isdir]);
subdirs = subdirs(3:end);
## adjust filenames to include directory info
% adjust filenames to include directory info
for j = 1:numel (d_info)
d_info(j).name = fullfile (d, d_info(j).name);
endfor
......
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=0)
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)
if nargin<5
redox=0;
end
......@@ -180,7 +182,7 @@ switch redox_buffer
QIF_ratio_H2_H2O=log10(H2_mole_fraction_QIF./H2O_mole_fraction_QIF);
else
QIF_ratio_H2_H2O=0;
endif
end
% oxygen fugacity at QIF buffer
......@@ -255,7 +257,7 @@ switch redox_buffer
IW_ratio_H2_H2O=log10(H2_mole_fraction_IW./H2O_mole_fraction_IW);
else
IW_ratio_H2_H2O=0;
endif
end
% oxygen fugacity at IW buffer
......@@ -326,7 +328,7 @@ switch redox_buffer
QFM_ratio_H2_H2O=log10(H2_mole_fraction_QFM./H2O_mole_fraction_QFM);
else
QFM_ratio_H2_H2O=0;
endif
end
% oxygen fugacity at QFM buffer
......@@ -397,7 +399,7 @@ switch redox_buffer
NiNiO_ratio_H2_H2O=log10(H2_mole_fraction_NiNiO./H2O_mole_fraction_NiNiO);
else
NiNiO_ratio_H2_H2O=0;
endif
end
% oxygen fugacity at NiNiO buffer
......
......@@ -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=1200; %600;
res=5000; %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
......@@ -80,7 +80,8 @@ 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('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('D2D_',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=-10,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 (reduc==0)
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
......
......@@ -41,7 +41,7 @@ function [data,input,program_dir,fname,fpath] = read_snap(finput=' ')
data.m_vel_u=0.0;,data.m_vel_v=0.0;,data.m_vel_w=0.0;,data.mprod=0.0;,data.poros=0.0;
while(strncmp(name,"EF",2)==0)
%name
name
if(name == "nl"), bin_b=fread(data_snap,1,'int32');, nl=fread(data_snap,1,'int32'), bin_e=fread(data_snap,1,'int32');
elseif(name == "ny"), bin_b=fread(data_snap,1,'int32');, ny=fread(data_snap,1,'int32'), bin_e=fread(data_snap,1,'int32');
elseif(name == "nr")
......
......@@ -34,15 +34,26 @@ M_carb = M_C + 3*M_O; %molar mass of carbonate CO_2^3-
% hydrous = 1;
%endif
% basaltic example composition, aparently in molar fractions, Saal et al 2012, D-6-1 sample
X_K2O = 0.00023;
X_Na2O = 0.0194;
X_CaO = 0.1201;
X_MgO = 0.1134;
X_FeO = 0.08;
X_Al2O3 = 0.164;
X_SiO2 = 0.4957;
X_TiO2 = 0.0084;
% basaltic example composition, all in weight fractions
X_K2O = 0.0772;
X_Na2O = 0.0195;
X_CaO = 0.114;
X_MgO = 0.0575;
X_FeO = 0.0782;
X_Al2O3 = 0.1557;
X_SiO2 = 0.4989;
X_TiO2 = 0.0089;
%X_K2O = 0.0772;
%X_Na2O = 0.0195;
%X_CaO = 0.114;
%X_MgO = 0.0575;
%X_FeO = 0.0782;
%X_Al2O3 = 0.1557;
%X_SiO2 = 0.4989;
%X_TiO2 = 0.0089;
%% Dixon 1997, MORB Thol. Basalt, wt-fractions:
%X_K2O = 0.07/100;
......@@ -59,22 +70,35 @@ X_TiO2 = 0.0089;
%X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Al2O3 + X_SiO2 + X_TiO2
%%get molar fractions
Sum_mol = X_K2O/(2*M_K+M_O) + X_Na2O/(2*M_Na+M_O) + X_CaO/(M_Ca+M_O) + X_MgO/(M_Mg+M_O) ...
+ X_FeO/(M_Fe+M_O) + X_Al2O3/(2*M_Al+3*M_O) + X_SiO2/(M_Si+2*M_O) + X_TiO2/(M_Ti+2*M_O) ...
+ X_H2O_wt/(2*M_H+M_O);
%Sum_mol = X_K2O/(2*M_K+M_O) + X_Na2O/(2*M_Na+M_O) + X_CaO/(M_Ca+M_O) + X_MgO/(M_Mg+M_O) ...
% + X_FeO/(M_Fe+M_O) + X_Al2O3/(2*M_Al+3*M_O) + X_SiO2/(M_Si+2*M_O) + X_TiO2/(M_Ti+2*M_O) ...
% + X_H2O_wt/(2*M_H+M_O);
%X_H2O = X_H2O_wt/Sum_mol/(2*M_H+M_O);
%X_K2O = X_K2O/Sum_mol/(2*M_K+M_O);
%X_Na2O = X_Na2O/Sum_mol/(2*M_Na+M_O);
%X_CaO = X_CaO/Sum_mol/(M_Ca+M_O);
%X_MgO = X_MgO/Sum_mol/(M_Mg+M_O);
%X_FeO = X_FeO/Sum_mol/(M_Fe+M_O);
%X_Al2O3 = X_Al2O3/Sum_mol/(2*M_Al+3*M_O);
%X_SiO2 = X_SiO2/Sum_mol/(M_Si+2*M_O);
%X_TiO2 = X_TiO2/Sum_mol/(M_Ti+2*M_O);
% if already in molar fractions:
Sum_mol = X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Al2O3 + X_SiO2 + X_TiO2 + X_H2O_wt/(2*M_H+M_O);
X_H2O = X_H2O_wt/Sum_mol/(2*M_H+M_O);
X_K2O = X_K2O/Sum_mol/(2*M_K+M_O);
X_Na2O = X_Na2O/Sum_mol/(2*M_Na+M_O);
X_CaO = X_CaO/Sum_mol/(M_Ca+M_O);
X_MgO = X_MgO/Sum_mol/(M_Mg+M_O);
X_FeO = X_FeO/Sum_mol/(M_Fe+M_O);
X_Al2O3 = X_Al2O3/Sum_mol/(2*M_Al+3*M_O);
X_SiO2 = X_SiO2/Sum_mol/(M_Si+2*M_O);
X_TiO2 = X_TiO2/Sum_mol/(M_Ti+2*M_O);
X_K2O = X_K2O/Sum_mol;
X_Na2O = X_Na2O/Sum_mol;
X_CaO = X_CaO/Sum_mol;
X_MgO = X_MgO/Sum_mol;
X_FeO = X_FeO/Sum_mol;
X_Al2O3 = X_Al2O3/Sum_mol;
X_SiO2 = X_SiO2/Sum_mol;
X_TiO2 = X_TiO2/Sum_mol;
% sum of all should be 1:
%X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Al2O3 + X_SiO2 + X_TiO2
X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Al2O3 + X_SiO2 + X_TiO2;
%% M_melt normalized to one oxygen (???)
%M_melt = (X_K2O*(2*M_K+M_O) + X_Na2O*(2*M_Na+M_O) + X_CaO*(M_Ca+M_O) ...
......@@ -138,13 +162,13 @@ ln_CO_3_2min = X_H2O*d_H2O + X_AI*d_AI + (X_FeO+X_MgO)*d_FeMg + (X_Na2O+X_K2O)*d
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
end
H2O_dissolved = exp(ln_H2O); % in wt-%
if (H2O_fl_mole==0), H2O_dissolved = 0;, endif
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;, endif
if (H2O_fl_mole==1), CO2_dissolved = 0;, end
%fwm=36.594; % KILAUEA THOLEIITIC BASALT, HOLLOWAY 1992
%%CO2_dissolved = CO32m_dissolved*1.202656 / (1+CO32m_dissolved*0.202656)
......
......@@ -24,7 +24,7 @@ H2O_dissolved = 0;
CO2_dissolved = 0;
tol = 10e-10;
i_max = 100;
i_max = 200;
%fprintf('Start with CO2_melt = %f wt-ppm and H2O_melt = %f wt-%%; p=%f, T=%f\n',...
% X_CO2,X_H2O,p_lid,T_lid)
......@@ -58,8 +58,10 @@ do
no_water=1;
H2O_fl_mole = 0;%0.5*H2O_fl_mole; % neither CO2 nor H2O go out, so set to smaller value to see if at least some CO2 goes out?
endif
% if (i>90)
% fprintf('%d H-C-mole = %f CO2_diss = %f H2O_diss = %f CO2_deg = %f%% H2O_deg = %f%%\n',i,...
% H2O_fl_mole,CO2_dissolved,H2O_dissolved,dX_CO2/X_CO2*100,dX_H2O/X_H2O*100)
% end
dX_H2O = max(0,X_H2O-H2O_dissolved);
dX_CO2 = max(0,X_CO2-CO2_dissolved);
i=i+1;
......
......@@ -1017,7 +1017,7 @@ module initialize
& ((0.5_dp*(rc(k)+rc(k+1)))**3-(0.5_dp*(rc(k)+rc(k-1)))**3)
enddo
write(*,*) time/(1.0e+6_dp*pF%time_yr), Q*1.0e+7_dp*(pF%k*pF%DeltaT)*pF%D,"erg s^-1" ! *k*DT/(D^2*rho) [W/kg] * D^3*rho (mass mantle) [kg] -> W = 10^7 erg/s
!write(*,*) time/(1.0e+6_dp*pF%time_yr), Q*1.0e+7_dp*(pF%k*pF%DeltaT)*pF%D,"erg s^-1" ! *k*DT/(D^2*rho) [W/kg] * D^3*rho (mass mantle) [kg] -> W = 10^7 erg/s
pH%Q_eddy = Q
......
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