Commit 1d27a1a7 authored by Claire Guimond's avatar Claire Guimond
Browse files

Merge branch 'master' of gitlab-as.oma.be:lenano/CHIC

parents d231b10c 198a57ad
ifndef server
server=local
else
ifneq ($(server),local) #planeto, odin, FU linux pc
ifneq ($(server),local) #planeto, odin, FU linux pc, FU PC Doc
ifneq ($(server),hlrn)
ifneq ($(server),hpc)
ifneq ($(server),plato)
......
% t,minval(cell%d),&
% &depl_av,maxval(cell%d),&
% &sum(mc%oc_vol)/(mesh%lmax-mesh%lmin)/(mesh%rmax-mesh%rmin),&
% &sum(mc%fc_vol)/(mesh%lmax-mesh%lmin)/(mesh%rmax-mesh%rmin),&
% &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
function melt_output_for_speciation(fpath=' ',all=1,redox="IW")
if (all==1)
%// List all sub-directories under MainFolder recursively
D = rdir(fpath); %// List of all sub-directories
for k = 1:length(D)
currpath = D(k).name; %// Name of current directory
if strcmp('data_melt_profiles.res',substr(currpath,-22))
% substr(currpath,1,length(currpath)-23)
run_speciation(substr(currpath,1,length(currpath)-23),redox)
end
end
else
run_speciation(fpath,redox)
end
end
function retval = rdir (d)
## info about D
d_info = dir (d);
## 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
for j = 1:numel (d_info)
d_info(j).name = fullfile (d, d_info(j).name);
endfor
retval = d_info;
for i = 1:numel (subdirs)
## recurse and append results for subdirectories to the list
d_info = rdir (fullfile (d, subdirs{i}));
retval = [retval; d_info];
endfor
endfunction
function run_speciation(fpath=' ',redox="IW")
program_dir = cd(); % store octave path to be able to return
%cd('/home/noacklen/Arbeit/SOURCE/CHIC/Octave');
cd('C:/Users/medhe/GIT/CHIC/Octave');
load('gianluigi_H20_500ppm_CO2_1000ppm.mat')
X_H2O_mantle = 500; %ppm % Needs to fit to speciation file grom Gianluigi or wrong pressures will come out
X_CO2_mantle = 1000; %ppm
% load('gianluigi_no_H20_only_CO2.mat')
% X_H2O_mantle = 0; %ppm
% X_CO2_mantle = 1000; %ppm
% load('gianluigi_only_H20_no_CO2.mat')
% X_H2O_mantle = 500; %ppm
% X_CO2_mantle = 0; %ppm
cd(program_dir);
if (fpath==' ')
[fname, fpath, fltidx] = uigetfile('data_melt_profiles.res');
else
fname = 'data_melt_profiles.res';
endif
cd(fpath);
data = load('-ascii',fname);
% for lid thickness
M = dlmread('data_temp_profs.res');
% output from screenlog:
% pF%Cp= 1270.4042841639362 , pF%k= 5.1123222586372110 , pF%g= 9.7488315239551309 , pF%rho= 4364.0674927388454 , pF%D= 3231564.5380912004 , pF%alpha= 2.0900769532248606E-005
% Rp= 6456645.6000000006 , Rc= 3225581.0619088002 , Cp_m= 1262.8001493999734 , Cp_c= 755.60358584686617
% Core density is 12816.270439083673 , mantle density is 4471.8398376656078 , average density is 5401.7509948414781
%rho = 4471.8398376656078;
% Ref values (D,t_yr,DeltaT,Ts,eta_ref):
% 3231564.5380912004 2.7865260888694197E-012 5720.0000000000000 280.00000000000000 1.0000000000000000E+021
%Rp= 6456645.6000000006;
%Rc= 3225581.0619088002;
%t_yr= 2.7865260888694197E-012;
%g= 9.7488315239551309
filenames = sprintf('%s/%s', fpath, 'output*');
files = glob(filenames);
if (size(files,1)==0)
fid = fopen('screenlog.0');
else
fid = fopen(files{1});
end
tline = fgetl(fid);
while ischar(tline)
if startsWith(tline,' pF%Cp=')
[start,remain] = strtok(tline,',');
[start,remain] = strtok(remain,',');
[g,remain] = strtok(remain,',');
g = g(7:end);
g = str2num(g);
end
if startsWith(tline,' Rp=')
% disp(tline)
[Rp_s,remain] = strtok(tline,',');
Rc_s = strtok(remain,',');
Rp_s = Rp_s(5:end);
Rc_s = Rc_s(5:end);
Rp = str2num(Rp_s)
Rc = str2num(Rc_s)
tline = fgetl(fid);
[Crho,remain] = strtok(tline,',');
Mrho = strtok(remain,',');
Mrho = Mrho(20:end);
rho = str2num(Mrho);
end
if startsWith(tline,' Ref values ')
% disp(tline)
tline = fgetl(fid);
% disp(tline)
arr = textscan(tline, '%f%f%f%f%f', 1);
arr2 = cell2mat(arr);
t_yr = arr2(2); %str2double(arr(2))
DeltaT = arr2(3)
Ts = arr2(4)
end
if startsWith(tline,' nl =')
% disp(tline)
[start,remain] = strtok(tline,',');
[start,remain] = strtok(remain,',');
[start,remain] = strtok(remain,'=');
nr_s = remain(5:end);
nr = str2num(nr_s);
end
tline = fgetl(fid);
end
fclose(fid);
% continue here with temp_profs, derivative. find lid thickness with temp
mantle_vol = 4/3*pi*(Rp^3-Rc^3);
nt = size(data,1);
melt_vol = 0; %! grow in interation to arrays as needed for true number of time steps
melt_rate = 0;%zeros(2,1);
time = data(:,1)/t_yr;
%%%%%%%%%%%%%%%%%%%%
% temp profiles
nr_sh = nr+2+1 % first line is time, plus two boundary cells
if (mod(nr_sh,4)==0)
nr_res = nr_sh;
else
nr_res = nr_sh-mod(nr_sh,4)+4;
end
%nr_res
A = reshape(M',nr_res,[]); % written in lines of 4, last line filled up with zeros
%size(A)
time_prof=A(1,:);
temp_prof=A(2:nr_sh+1,:);
%nr_sh = length(temp_prof(:,1))
r = linspace(Rc,Rp,nr_sh)';
% index where lithsophere ends, hence where temp profile changes from adiabatic to conductive
[maxv ind_l] = max((temp_prof(:,1)*DeltaT+Ts)-DeltaT/(Rp-Rc)*(Rp-r));
D_lid = (Rp-r(ind_l))/1000 % lid thickness in km
p_lid = D_lid*1000*g*rho % Pa
T_lid = temp_prof(ind_l,1) % temp at lid base in K
M_H = 1.0079;
M_C = 12.0107;
M_O = 15.9994;
M_CO2 = M_C + 2*M_O;
M_H2O = 2*M_H + M_O;
%X_H2O_mantle = 500; % ppm -> wt or mole?
%X_CO2_mantle = 1000; % ppm
% solubility calculated from Iacono-Marzian et al., 2012; H2O in wt-% and CO2 in wt-ppm
%X_H2O = 10 %weight percent of water in the melt
%X_CO2 = 10000 %weight ppm of CO2 in the melt
%hydrous = 0;
%H2O_fl_mole = 0.5; % = P_H2O/P_lid
%[H2O_dissolved,CO2_dissolved] = solubility(p_lid*10^-5,T_lid,H2O_fl_mole,X_H2O,hydrous);
%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
%dX_H2O_mole = dX_H2O/M_H2O/10^2 / (dX_H2O/M_H2O/10^2 + dX_CO2/M_CO2/10^6)
%dX_CO2_mole = dX_CO2/M_CO2/10^6 / (dX_H2O/M_H2O/10^2 + dX_CO2/M_CO2/10^6)
%H2O_fl_mole = dX_H2O_mole / (dX_H2O_mole+dX_CO2_mole)
%%%%%%%%%%%%%%%%%%%%
%time_plot(1) = time(1);
%time_restart = 0;
%j=1;
%for i=1:nt
% if time(i)>time_restart
% melt_vol(j) = data(i,12)*mantle_vol*rho;
% time_plot(j)=time(i);
% if (j==1)
% melt_rate(j) = melt_vol(j)/rho/time_plot(j);
% else
% melt_rate(j) = melt_vol(j)/rho/(time_plot(j)-time_plot(j-1));
% endif
% j=j+1;
% time_restart = time(i);
% end
%end
time_plot(1) = time(nt);
time_restart = time(nt);
j=1;
for i_flip=1:nt
i = nt+1-i_flip; % go backwarts from last time step to beginning
if time(i)<time_restart
j=j+1;
melt_vol(j) = data(i,12)*mantle_vol*rho;
time_plot(j)=time(i);
melt_rate(j) = melt_vol(j)/rho/(time_plot(j-1)-time_plot(j));
time_restart = time(i);
end
end
time_plot = flip(time_plot);
melt_rate = flip(melt_rate);
melt_vol = flip(melt_vol);
%for i=2:nt
% melt_vol(i) = data(i,12)*mantle_vol*rho;
% if (i==1)
% melt_rate(i) = melt_vol(i)/rho/time(i);
% else
% melt_rate(i) = melt_vol(i)/rho/(time(i)-time(i-1));
% endif
%end
%figure
%plot(time,melt_rate*10^-9)
Chi_extr = 0.1; % extrusive volcanism
mass_H2O = 2.0*1.00794 + 15.999;
mass_CO2 = 12.0107 + 2.0*15.999;
mass_H2 = 2.0*1.00794;
mass_CO = 12.0107 + 15.999;
X_H2O_mole = X_H2O_mantle/mass_H2O / (X_H2O_mantle/mass_H2O+X_CO2_mantle/mass_CO2); % mole amount of X_H2O and X_CO2
X_CO2_mole = X_CO2_mantle/mass_CO2 / (X_H2O_mantle/mass_H2O+X_CO2_mantle/mass_CO2); % sum of both is 1
% 1600 K is last (100th) entry in temp array
%%% IW
ind = 0;
for i=2:length(T)
if T(i-1)<1600 && T(i)>=1600
ind=i;
end
end
%ToDo: add linear interpolation case between different redox states (e.g. QFM->IW or IW->QFM)
if strcmp(redox,'QIF')
H2_mole = H2_mole_fraction_QIF(ind);
H2O_mole = H2O_mole_fraction_QIF(ind);
CO_mole = CO_mole_fraction_QIF(ind);
CO2_mole = CO2_mole_fraction_QIF(ind);
elseif strcmp(redox,'IW')
H2_mole = H2_mole_fraction_IW(ind);
H2O_mole = H2O_mole_fraction_IW(ind);
CO_mole = CO_mole_fraction_IW(ind);
CO2_mole = CO2_mole_fraction_IW(ind);
elseif strcmp(redox,'QFM')
H2_mole = H2_mole_fraction_QFM(ind);
H2O_mole = H2O_mole_fraction_QFM(ind);
CO_mole = CO_mole_fraction_QFM(ind);
CO2_mole = CO2_mole_fraction_QFM(ind);
elseif strcmp(redox,'NiNiO')
H2_mole = H2_mole_fraction_NiNiO(ind);
H2O_mole = H2O_mole_fraction_NiNiO(ind);
CO_mole = CO_mole_fraction_NiNiO(ind);
CO2_mole = CO2_mole_fraction_NiNiO(ind);
else
fprintf('ERROR - choose one of the following redox states: QIF-IW-QFM-NiNiO')
endif
%H2_mole
%H2O_mole
%CO_mole
%CO2_mole
mu_av = H2_mole*mass_H2 + H2O_mole*mass_H2O + CO_mole*mass_CO + CO2_mole*mass_CO2;
M_H2 = X_H2O_mantle/mass_H2O * H2_mole/X_H2O_mole * mass_H2 * melt_vol * 1e-6 * Chi_extr; % mass of H2 in kg
%P_H2 = M_H2 * g / (4*pi*Rp^2) * 1e-5;
M_H2O = X_H2O_mantle/mass_H2O * H2O_mole/X_H2O_mole * mass_H2O * melt_vol * 1e-6 * Chi_extr; % mass of H2O in kg
%P_H2O = M_H2O * g / (4*pi*Rp^2) * 1e-5;
M_CO = X_CO2_mantle/mass_CO2 * CO_mole/X_CO2_mole * mass_CO * melt_vol * 1e-6 * Chi_extr; % mass of H2 in kg
%P_CO = M_CO * g / (4*pi*Rp^2) * 1e-5;
M_CO2 = X_CO2_mantle/mass_CO2 * CO2_mole/X_CO2_mole * mass_CO2 * melt_vol * 1e-6 * Chi_extr; % mass of H2O in kg
%P_CO2 = M_CO2 * g / (4*pi*Rp^2) * 1e-5;
M_tot = M_H2 + M_H2O + M_CO + M_CO2; % total mass in kg
% partial pressure of species "i" depends on mole fraction, not mass fraction!
P_H2 = M_tot * H2_mole * g / ( 4*pi*Rp^2 ) * 1e-5;
P_H2O = M_tot * H2O_mole * g / ( 4*pi*Rp^2 ) * 1e-5;
P_CO = M_tot * CO_mole * g / ( 4*pi*Rp^2 ) * 1e-5;
P_CO2 = M_tot * CO2_mole * g / ( 4*pi*Rp^2 ) * 1e-5;
%P_tot_H2 = sum(P_H2);
%P_tot_H2O = sum(P_H2O);
%P_tot_CO = sum(P_CO);
%P_tot_CO2 = sum(P_CO2);
P_tot_H2(1) = P_H2(1);
P_tot_H2O(1) = P_H2O(1);
P_tot_CO(1) = P_CO(1);
P_tot_CO2(1) = P_CO2(1);
%time_plot(1) = time(1);
%time_restart = 0;
%j=2;
for i=2:length(P_H2)
% if time(i)>time_restart
P_tot_H2(i)=P_tot_H2(i-1)+P_H2(i);
P_tot_H2O(i)=P_tot_H2O(i-1)+P_H2O(i);
P_tot_CO(i)=P_tot_CO(i-1)+P_CO(i);
P_tot_CO2(i)=P_tot_CO2(i-1)+P_CO2(i);
% time_plot(j)=time(i);
% j=j+1;
% time_restart = time(i);
% end
end
P_tot_H2(nt)
P_tot_H2O(nt)
P_tot_CO(nt)
P_tot_CO2(nt)
%f = figure('visible','off');
%plot(time_plot*1e-9,P_tot_H2)
%hold on
%plot(time_plot*1e-9,P_tot_H2O)
%plot(time_plot*1e-9,P_tot_CO)
%plot(time_plot*1e-9,P_tot_CO2)
%hold off
%legend('H2','H2O','CO','CO2')
%xlabel('Time in Gyr')
%ylabel('Partial pressure in bar')
%mkdir 'plots';
%print(sprintf('plots/Profiles_Outgassing.png',f));
% write all output files in one main folder
path_write = 'C:\OwnCloud_FUB\Literature_Projects\Eigene_Projekte\Project_Ortenzi_Dorn\OctaveResultsNew2\';
%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
%name_file = substr(fpath,64); % /scratch-2/PublishedProjects/Project_CaroDorn/Project_CaroDorn_
%[name1,name2] = strtok(name_file,'/'); % Review /New_Case10_Visc10/New_Run_1.5M_LowAcc/RC0809_FS0.5_MS1.0
%[name2,name3] = strtok(name2,'/'); % New_Case10_Visc10 /New_Run_1.5M_LowAcc/RC0809_FS0.5_MS1.0
%[name3,name4] = strtok(name3,'/'); % New_Run_1.5M_LowAcc /RC0809_FS0.5_MS1.0
%[name4,residual] = strtok(name4,'/'); % RC0809_FS0.5_MS1.0
%name2=name2(9:end); % delete "New_Case"
%name3=name3(9:end); % delete "New_Run_"
%new_name = sprintf('Case%s--%s--%s.txt',name2,name3,name4);
name1=' ';,name2=' ';,name3=' ';,name4=' ';,new_name = 'Test_old';
cd(path_write)
fid=fopen(new_name,'w');
fprintf(fid,'time (Gyr) P_H2 (bar) P_H2O (bar) P_CO (bar) P_CO2(bar) frac_H2 frac_H2O frac_CO frac_CO2\n')
time_print = 1e+7; % 450 lines output
fprintf(fid,"%3.8f %3.8f %3.8f %3.8f %3.8f %3.8f %3.8f %3.8f %3.8f\n",time_plot(1)*1e-9,P_tot_H2(1),P_tot_H2O(1),P_tot_CO(1),P_tot_CO2(1),H2_mole,H2O_mole,CO_mole,CO2_mole)
np = length(time_plot);
for i=2:np
if time_plot(i)>time_print
fprintf(fid,"%3.8f %3.8f %3.8f %3.8f %3.8f %3.8f %3.8f %3.8f %3.8f\n",time_plot(i)*1e-9,P_tot_H2(i),P_tot_H2O(i),P_tot_CO(i),P_tot_CO2(i),H2_mole,H2O_mole,CO_mole,CO2_mole)
time_print = time_print + 1e+7;
end
end
fclose(fid);
fid=fopen('all_results.txt','a');
fprintf(fid,"%s %s %s %3.8f %3.8f %3.8f %3.8f %3.8f %3.8f %3.8f %3.8f %3.8f %3.8f %3.8f %3.8f %3.8f\n",...
name2,name3,name4,Rp,Rc,g,rho,P_tot_H2(np),P_tot_H2O(np),P_tot_CO(np),P_tot_CO2(np),H2_mole,H2O_mole,CO_mole,CO2_mole,mu_av)
fclose(fid);
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
endif
retval = strncmp(s, prefix, n);
endfunction
This diff is collapsed.
%simultaneous equilibria Fegley 2013 pag. 458
% to demonstrate that XCH4 is the same with all the equilibria. If yes, I
% can use the equilibrium without oxygen fugacity to calcalute the C-O-H
% system (verified. using only T and fO2 it was possible to calculate: 1) the wt% of CO, CO2 (from CO +1/2O2= CO2) and the wt% of H2 and H2O (2H2 + O2 = 2H2O). 2) the wt% of CH4 with the others 4 equilibrium (also direct and inverse reaction)
initial_H2O_mole_fraction=0.5;
initial_CO2_mole_fraction=0.5;
%T=Temperature in Kelvin, P=Pressure is in Bar
%The parameters are for the QFM buffer. It is possible to change the delta
%value of the buffer.
% file (fO2_ferric_ferrous ratio calculation).
T=1200;
P=1;
buffer=0;
%Equation for the calculation of the delta G of formation of carbon monoxyde 2C (graphite) + 2H2 = CH4 (Fegley 2013)
%P is 1 bar, T is in Kelvin. Range of T from 298 to 2500 K
% the equilibrium constant for reactions of ideal gases is constant at any
% given temperature and is unaffected by the total pressure ( Fegley 2013,
% pag 446). However, the equilibrium partial pressure of reactants and
% products are affected by the total pressure. Especially (P, As, Sb, Bi,
% S, and Se)
% 2C (graphite) + 2H2 = CH4
% delta G = A + BTLogT +CT in J/mol (A, B, C are parameters calculated in Fegley 2013)
% results checked with JANAF tables
Gf_CH4=-77437+22.5098.*T.*log10(T)+29.5967.*T;
%Equation for the calculation of the delta G of formation of carbon monoxyde 2C (graphite) + O2 = 2CO (Fegley 2013. pag 437). Delta G of formation is the change in the Gibbs energy due to the formation of one mole of a substance starting from the constituents (pure elements are always zero).
%P is 1 bar, T is in Kelvin. Range of T from 298 to 2500 K
% the equilibrium constant for reactions of ideal gases is constant at any
% given temperature and is unaffected by the total pressure ( Fegley 2013,
% pag 446). However, the equilibrium partial pressure of reactants and
% products are affected by the total pressure. Especially (P, As, Sb, Bi,
% S, and Se)
% 2C (graphite) + O2 = 2CO
% delta G = A + BTLogT +CT in J/mol (A, B, C are parameters calculated in Fegley 2013)
% results checked with JANAF tables
g_co=-214104+25.2183.*T.*log10(T)-262.1545.*T;
Gf_CO=g_co./2;
%Equation for the calculation of the delta G of formation of carbon dioxyde
%C (graphite) + O2 = CO2 (Fegley 2013) . Delta G of formation is the change in the Gibbs energy due to the formation of one mole of a substance starting from the constituents (pure elements are always zero).
%P is 1 bar, T is in Kelvin. Range of T from 298 to 2500 K
% the equilibrium constant for reactions of ideal gases is constant at any
% given temperature and is unaffected by the total pressure ( Fegley 2013,
% pag 446). However, the equilibrium partial pressure of reactants and
% products are affected by the total pressure. Especially (P, As, Sb, Bi,
% S, and Se)
% C (graphite) + O2 = CO2
% delta G = A + BTLogT +CT in J/mol (A, B, C are parameters calculated in Fegley 2013)
% results checked with JANAF tables
Gf_CO2=-392647+4.5855.*T.*log10(T)-16.9762.*T;
%Equation for the calculation of the delta G of formation of water 2H2 + O2 = 2H2o (Fegley 2013)
%P is 1 bar, T is in Kelvin. Range of T from 298 to 2500 K
% the equilibrium constant for reactions of ideal gases is constant at any
% given temperature and is unaffected by the total pressure ( Fegley 2013,
% pag 446). However, the equilibrium partial pressure of reactants and
% products are affected by the total pressure. Especially (P, As, Sb, Bi,
% S, and Se)
% 2H2 + O2 = 2H2O
% delta G = A + BTLogT +CT in J/mol (A, B, C are parameters calculated in Fegley 2013)
% results checked with JANAF tables
g_h2o=-483095+25.3687.*T.*log10(T)+21.9563.*T;
Gf_H2O=g_h2o./2;
% Calculation of the species mole per cent of the C-O-H system (Fegley 2013)
%1) 2H2 +O2 = 2H2O
% Gf_H2O * 2 because there are 2 moles of water
% K1 = equilibrium constant for reaction 1
K1=exp(((-Gf_H2O).*2)./(8.314.*T));
% 2) CO + 1/2O2 = CO2
%oxygen partial pressure (i.e an oxygen fugacity) for the reaction. Going
%to low ox fugacity the level of CO2 decrease.
%the delta G formations come from the functions in deltaG_formation folder
deltaG_reaction_CO_CO2=Gf_CO2-Gf_CO;
% K2 = equilibrium constant for reaction 2
K2=exp(-deltaG_reaction_CO_CO2./(8.314.*T));
%equation and parameters from Holloway 1992
%
%the equation for calculating the equilibria constant fits to better than +-0.01 in units of logfO2 for the range 1000< T <1700°C (Holloway 1992)
redox_buffer = input('Enter a number for a Redox buffer (1 for QIF, 2 for IW, 3 for QFM, 4 NiNiO) : ');
switch redox_buffer
case 1
LogfO2_QIF=((-29673./T)+7.679+((0.05*(P-1)./T))+buffer);
O2_fugacity_QIF=10.^LogfO2_QIF;
% oxygen fugacity at QIF buffer
ratio_H2O_H2_QIF=sqrt((O2_fugacity_QIF).*K1);
H2_mole_fraction_QIF=(initial_H2O_mole_fraction./(ratio_H2O_H2_QIF+1)); % (initial H2O mole fraction is from the solubility model)
H2_wt_percent_QIF=H2_mole_fraction_QIF.*100;
H2O_mole_fraction_QIF=initial_H2O_mole_fraction-H2_mole_fraction_QIF;
H2O_wt_percent_QIF=H2O_mole_fraction_QIF.*100;