Commit 632bcef5 authored by Lena Noack's avatar Lena Noack
Browse files

Bugfix degassing results / Bugfix redistribute particles for empty cells

parent 3d7155c5
......@@ -64,6 +64,7 @@ function data = get_data(data,input)
if(size(data.depl,1)>1), data.depl = cut_field(data.depl,nl,ny,nr)*100.0;, endif % depletion in percent
if(size(data.DeltaF,1)>1), data.DeltaF = cut_field(data.DeltaF,nl,ny,nr)*100.0;, endif % melt fraction of last time step in percent
if(size(data.cr_age,1)>1), data.cr_age = cut_field(data.cr_age,nl,1,1)*data.D^2.0/data.kappa/3600.0/24.0/365.0/1000000.0 ;, endif %2D matrix, no radial dependence; Age in Myr
if(size(data.oc_vol_tot,1)>1), data.oc_vol_tot = cut_field(data.oc_vol_tot,nl,1,1)*data.D^2;, endif % vol in 2D multiplied with D^2, in 3D would need to be D^3
if(size(data.entropy,1)>1), data.entropy = cut_field(data.entropy,nl,ny,nr)*data.rho*data.gsurf*data.D^4.0/data.DT;, endif % entropy
if(size(data.CO2,1)>1), data.CO2 = cut_field(data.CO2,nl,ny,nr);, endif % in ppm (?)
if(size(data.ref_r,1)>1), data.ref_r = cut_field(data.ref_r,nl,ny,nr)*data.rho;, endif
......
......@@ -42,17 +42,32 @@ function melt_output_for_speciation_lid(fpath,all,redox)
if isfolder(currpath)
cd(currpath);
if exist('data_melt_profiles.res')
currpath
run_speciation(currpath,redox)
% run_speciation(sprintf('%s/%s',D(k).folder,D(k).name),redox)
else
% subfolders
D2 = dirr(currpath); %// List of all sub-directories
for k = 1:length(D2)
currpath = sprintf('%s/%s',D2(k).folder,D2(k).name)
for j = 1:length(D2)
currpath = sprintf('%s/%s',D2(j).folder,D2(j).name);
if isfolder(currpath)
cd(currpath);
if exist('data_melt_profiles.res')
currpath
run_speciation(currpath,redox)
else
% subsubfolders
D3 = dirr(currpath); %// List of all sub-directories
for i = 1:length(D3)
currpath = sprintf('%s/%s',D3(i).folder,D3(i).name);
if isfolder(currpath)
cd(currpath);
if exist('data_melt_profiles.res')
currpath
run_speciation(currpath,redox)
end
end
end
end
end
end
......@@ -129,13 +144,16 @@ files = dirr(filenames);%glob(filenames);
if (size(files,1)==0)
fid = fopen('screenlog.0');
path = 'screenlog.0';
else
path1 = files.folder;
path2 = files.name;
path = sprintf('%s/%s',path1(1,1:end),path2(1,1:end))
path = sprintf('%s/%s',path1(1,1:end),path2(1,1:end));
fid = fopen(path);
end
%path
% 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
......@@ -573,7 +591,8 @@ dr_r1=dr_r(end);
% write all output files in one main folder
%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 = '/home/noacklen/ownCloud/Literature_Projects/Eigene_Projekte/Project_Ortenzi_Dorn/OctaveResultsNew3/';
path_write = '/scratch-2/PublishedProjects/Project_CaroDorn/Results_New/';
%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
......@@ -597,7 +616,7 @@ new_name = sprintf('Case%s--%s--%s--%s.txt',name2,name3,name4,redox);
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) ', ...
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_av frac_H2O_av frac_CO_av frac_CO2_av 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 ',...
'M_H2_in_m M_H2O_in_m M_CO_in_m M_CO2_in_m M_H2_extr M_H2O_extr M_CO_extr M_CO2_extr ',...
'Mt_H2_in_l(kg/yr) Mt_H2O_in_l Mt_CO_in_l Mt_CO2_in_l Mt_H2_in_c Mt_H2O_in_c Mt_CO_in_c Mt_CO2_in_c ',...
......@@ -648,7 +667,7 @@ fclose(fid);
% before executing this function
if exist('all_results.txt')==0
fid=fopen('all_results.txt','w');
fprintf(fid,'Case : time(Gyr) | P_atm(bar) | P_H2 | P_H2O | P_CO | P_CO2 | mu_av | M_tot(kg) | d_atm(km) | d_atm(Rp) | crust(km) | d_lith(km) | M_H2_l | M_H2O_l | M_CO_l | M_CO2_l | M_H2_c | M_H2O_c | M_CO_c | M_CO2_c | M_H2_m | M_H2O_m | M_CO_m | M_CO2_m | M_H2 | M_H2O | M_CO | M_CO2 | H2_mole_frac_av | H2O_mole_frac_av | CO_mole_frac_av | CO2_mole_frac_av\n');
fprintf(fid,'Case: time(Gyr) P_atm(bar) P_H2 P_H2O P_CO P_CO2 mu_av M_tot(kg) d_atm(km) d_atm(Rp) crust(km) d_lith(km) M_H2_l M_H2O_l M_CO_l M_CO2_l M_H2_c M_H2O_c M_CO_c M_CO2_c M_H2_m M_H2O_m M_CO_m M_CO2_m M_H2 M_H2O M_CO M_CO2 H2_mole_frac_av H2O_mole_frac_av CO_mole_frac_av CO2_mole_frac_av\n');
fclose(fid);
end
......@@ -755,4 +774,4 @@ end
deg_CO2 = dX_CO2; % wt-ppm
deg_H2O = dX_H2O*10^4; % wt-% -> wt-ppm
end
\ No newline at end of file
end
No preview for this file type
......@@ -61,17 +61,18 @@ function plot_profile_melt(D=1,atmos=0)
%%%%%%%%%%%%%%
if (atmos>0)
%Time (Gyr) T_mantle(K) maxF massH2O totalH2O massCO2 totalCO2 PH2O(bar) PCO2(bar) max_CO2_frac esc_H2O_(kg/Gyr) esc_CO2_(kg/Gyr) tot_H2O_atm tot_CO2_atm PH2O_atm_end PCO2_atm_end EGL_H2O
% data_atmos = load("data_atmosphere.res");
% 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_deg | X_CO2_deg
data_atmos = dlmread("data_atmosphere.res", "", 4, 0);
data_atmos(1,3:17) = 0.0;
data_atmos(1,3:end) = 0.0;
endif
%%%%%%%%%%%%%
% Melt data %
%%%%%%%%%%%%%
% 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
% add_crust!=0: t,min_depl,depl_av,max_depl,OC_av,FC_av,CC_av,mass_H2O,total_mass_H2O,water_av,min DeltaF,DeltaF_av,max DeltaF,cr_av_OC,cr_av_FC,cr_av_CC,oc_vol_tot,fc_vol_tot,cc_vol_tot
% add_crust==0: t,min_depl,depl_av,max_depl,OC_av,FC_av,CC_av,mass_H2O,total_mass_H2O,water_av,min DeltaF,DeltaF_av,max DeltaF,T at 1.5*s_bot,DeltaF at 1.5*s,w at 1.5*s_top,0.0_dp,0.0_dp,0.0_dp
......@@ -81,6 +82,35 @@ function plot_profile_melt(D=1,atmos=0)
if (atmos>0)
data_melt(:,1) = data_melt(:,1)*max(data_atmos(:,1))/max(data_melt(:,1));
endif
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define index for restart %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
idx_m = ones(size(data_melt(:,1)));
idx_a = ones(size(data_atmos(:,1)));
min_time = data_melt(end,1);
for i=length(data_melt(:,1))-1:-1:1
if min_time > data_melt(i,1)
min_time = data_melt(i,1);
else
idx_m(i)=0; % skip this one, since restart leads to overwriting these output lines
end
end
idx_m = logical(idx_m); % create index vector with true/false entries
min_time = data_atmos(end,1);
for i=length(data_atmos(:,1))-1:-1:1
if min_time > data_atmos(i,1)
min_time = data_atmos(i,1);
else
idx_a(i)=0; % skip this one, since restart leads to overwriting these output lines
end
end
idx_a = logical(idx_a); % create index vector with true/false entries
%%%%%%%%%%%%%%%%%
% Plot profiles %
......@@ -90,62 +120,158 @@ function plot_profile_melt(D=1,atmos=0)
f = figure('visible','off');
max_time=max(data_atmos(:,1));
dt = [data_atmos(1);data_atmos(2:end,1)-data_atmos(1:end-1,1)];
dtm = [data_melt(1);data_melt(2:end,1)-data_melt(1:end-1,1)];
if (atmos>0)
% subplot(3,2,1)
subplot(2,2,1)
plot(data_atmos(:,1),data_atmos(:,2),"linewidth",3);
subplot(4,2,1)
% subplot(2,2,1)
% plot(data_atmos(:,1),data_atmos(:,2),"linewidth",3);
% xlabel ("Time (Gyr)");
% ylabel ("T_{mean} (K)");
% axis([0,max_time,min(data_atmos(:,2)),max(data_atmos(:,2))]);
semilogy(data_atmos(idx_a,1),data_atmos(idx_a,21),"linewidth",3);
hold on
semilogy(data_atmos(idx_a,1),data_atmos(idx_a,22),"linewidth",3);
hold off
xlabel ("Time (Gyr)");
ylabel ("T_{mean} (K)");
axis([0,5,min(data_atmos(:,2)),max(data_atmos(:,2))]);
ylabel ("X_{H2O}/X_{CO2} degassed");
% axis([0,max_time,0,0.1]);
subplot(4,2,2)
% plot(data(:,1)/1000,data(:,3),"linewidth",3);
% hold on
% plot(data(:,1)/1000,data(:,4),'r',"linewidth",3);
% hold off
% xlabel ("Time (Gyr)");
% ylabel ("Vm (k), V0 (r)");
%% axis([0,max_time,min(data(:,6)),max(data(:,4))]);
plot(data_atmos(idx_a,1),data_atmos(idx_a,4),"linewidth",3);
hold on
plot(data_atmos(idx_a,1),data_atmos(idx_a,8),"linewidth",3);
plot(data_atmos(idx_a,1),data_atmos(idx_a,6),"linewidth",3);
plot(data_atmos(idx_a,1),data_atmos(idx_a,10),"linewidth",3);
hold off
xlabel ("Time (Gyr)")
ylabel ("Mass H2O/H2/CO2/CO")
legend ("H2O","H2","CO2","CO")
% subplot(3,2,2)
subplot(2,2,4)
plot(data_atmos(:,1),data_atmos(:,8),"linewidth",3);
subplot(4,2,3)
% subplot(2,2,4)
plot(data_atmos(idx_a,1),data_atmos(idx_a,12),"linewidth",3);
hold on
plot(data_atmos(idx_a,1),data_atmos(idx_a,14),"linewidth",3);
plot(data_atmos(idx_a,1),data_atmos(idx_a,13),"linewidth",3);
plot(data_atmos(idx_a,1),data_atmos(idx_a,15),"linewidth",3);
hold off
xlabel ("Time (Gyr)");
ylabel ("P_{H2O} (bar)");
axis([0,5,0,max(data_atmos(:,8))]);
ylabel ("P_{H2O}/P_{H2}/P_{CO2}/P_{CO}");
% axis([0,max_time,0,max([max(data_atmos(:,12)),max(data_atmos(:,14))])]);
subplot(4,2,4)
% plot(data_atmos(:,1),data_atmos(:,13),"linewidth",3);
% hold on
% plot(data_atmos(:,1),data_atmos(:,15),'r',"linewidth",3);
% hold off
% xlabel ("Time (Gyr)");
% ylabel ("P_{CO2} (bar) / CO");
% axis([0,max_time,0,max([max(data_atmos(:,13)),max(data_atmos(:,15))])]);
% if only H2O, how much would it be (to compare shape of degassing):
P_H2O_total = 4/3*(Rp^3-Rc^3)*rho*((input.H2O-data_melt(:,10))/1000000)*grav/(4*Rp^2)/100000*input.Chi_extr;
subplot(4,2,4)
plot(data_melt(idx_m,1),P_H2O_total(idx_m),'b:',"linewidth",3);
xlabel ("Time (Gyr)");
ylabel ("P deg H2O tot");
subplot(2,2,3)
% subplot(3,2,3)
plot(data_atmos(:,1),data_atmos(:,9),"linewidth",3);
% subplot(2,2,3)
subplot(4,2,5)
% plot(data_atmos(:,1),data_atmos(:,7),"linewidth",3);
% xlabel ("Time (Gyr)");
% ylabel ("P_{CO2} (bar)");
% axis([0,max_time,0,max(data_atmos(:,7))]);
% plot(data_melt(:,1),data_melt(:,11),"linewidth",3);
% hold on;
idx_m2 = idx_m;
idx_m2(1:10)=0;
plot(data_melt(idx_m2,1),data_melt(idx_m2,12)./dtm(idx_m2),"k","linewidth",3);
% plot(data_melt(:,1),data_melt(:,13),"r--","linewidth",3);
% hold off;
% legend("Min","Mean","Max","location","northwest");
% axis ([0, max_time, 0, 0.3]);
% axis ([0, max_time, 0, 5]);
xlabel ("Time (Gyr)");
ylabel ("P_{CO2} (bar)");
axis([0,5,0,max(data_atmos(:,9))]);
ylabel ("Mean melt fraction");
% subplot(3,2,4)
% plot(data_atmos(:,1),data_atmos(:,17),"linewidth",3);
subplot(4,2,6)
% plot(data_atmos(:,1),data_atmos(:,16),"linewidth",3);
% xlabel ("Time (Gyr)");
% ylabel ("EGL_{H2O} (m)");
% axis([0,max_time,0,max(data_atmos(:,16))]);
% plot(data_atmos(:,1),dt,"linewidth",3);
% hold on
% plot(data_melt(:,1),dtm,'r',"linewidth",3);
% hold off
% xlabel ("Time (Gyr)");
% ylabel ("Time step");
plot(data_atmos(idx_a,1),data_atmos(idx_a,17),"linewidth",3);
hold on
plot(data_atmos(idx_a,1),data_atmos(idx_a,19),"linewidth",3);
plot(data_atmos(idx_a,1),data_atmos(idx_a,18),"linewidth",3);
plot(data_atmos(idx_a,1),data_atmos(idx_a,20),"linewidth",3);
hold off
xlabel ("Time (Gyr)");
ylabel ("frac H2O/H2/CO2/CO");
% subplot(3,2,5)
subplot(4,2,7)
plot(data_melt(idx_m,1),data_melt(idx_m,3),"linewidth",3);
% plot(data_melt(:,1),data_melt(:,2),"linewidth",3);
% hold on;
hold on
% plot(data_melt(:,1),data_melt(:,3),"k-.","linewidth",3);
% plot(data_melt(:,1),data_melt(:,4),"r--","linewidth",3);
% hold off;
plot(data_melt(idx_m,1),cumsum(data_melt(idx_m,12)),"m:","linewidth",3);
hold off
% legend("Min","Mean","Max","location","northwest");
% axis ([0, max(data_melt(:,1)), 0, 0.3]);
% xlabel ("Time (Gyr)");
% ylabel ("Melt depletion");
legend("Mean D","Sum DF","location","southeast");
% axis ([0, max_time, 0, 0.3]);
xlabel ("Time (Gyr)");
ylabel ("Melt depletion");
subplot(2,2,2)
% subplot(3,2,6)
plot(data_melt(:,1),data_melt(:,10),"linewidth",3);
% subplot(2,2,2)
subplot(4,2,8)
% plot(data_melt(:,1),data_melt(:,10),"linewidth",3);
% xlabel ("Time (Gyr)");
% ylabel ("Av. mantle water (ppm)");
% % axis([0,max_time,min(data_melt(:,10)),max(data_melt(:,10))]);
% min(data_melt(:,10))
% max(data_melt(:,10))
plot(data_atmos(idx_a,1),data_atmos(idx_a,4)./dt(idx_a)/18.01488,"linewidth",3);
hold on
plot(data_atmos(idx_a,1),data_atmos(idx_a,8)./dt(idx_a)/2.01588,"linewidth",3);
plot(data_atmos(idx_a,1),data_atmos(idx_a,6)./dt(idx_a)/44.0087,"linewidth",3);
plot(data_atmos(idx_a,1),data_atmos(idx_a,10)./dt(idx_a)/28.0097,"linewidth",3);
hold off
xlabel ("Time (Gyr)");
ylabel ("Av. mantle water (ppm)");
axis([0,5,min(data_melt(:,10)),max(data_melt(:,10))]);
ylabel ("Mole/dt H2O/H2/CO2/CO");
else
plot(data_melt(:,1),data_melt(:,10),"linewidth",3);
plot(data_melt(idx_m,1),data_melt(idx_m,10),"linewidth",3);
xlabel ("Time (Gyr)");
ylabel ("Av. mantle water (ppm)");
% axis([0,5,min(data_melt(:,10)),max(data_melt(:,10))]);
axis([0,max_time,min(data_melt(idx_m,10)),max(data_melt(idx_m,10))]);
endif
print(sprintf('plots/Profiles_Melt.png',f));
cd(program_dir)
\ No newline at end of file
% cd(program_dir)
\ No newline at end of file
......@@ -25,6 +25,9 @@ function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
mkdir "plots";
graphics_toolkit gnuplot
size(data.depl)
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)));
......@@ -45,13 +48,42 @@ function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
min_T=min(min(data.T));
if(max_T==0),max_T=max(max(data.T));, endif
if (reduc<=1)
if (reduc<=2)
% useful colormaps are for example hot, copper, jet (also viridis, but load extra?)
produce_plot('T2D_',"Temperature in K",data.T,input,data,res,ar,50,'jet',min_z=min_T,max_z=max_T,full=fullV)
produce_plot('Vi2D_',"Viscosity in log_{10} Pa\ s",log10(data.visc),input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)
produce_plot('V2D_',"Velocity in cm/yr",data.vel,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)
% data.cr_age only surface values, need to be set everywhere where there is crust
% data.oc_vol_tot is total volume per lateral/perpendicular cell, need to get depth of crust
nl = size(data.T,1);
nr = size(data.T,3);
cr_depth = zeros(nl,1);
cr_vol = data.oc_vol_tot;
dl = data.p_lat(2)-data.p_lat(1)
dr = data.p_rad(2)-data.p_rad(1)
for i=1:nl % in 3D would need do add for-loop for size(data.T,2)
for k=nr:-1:1
vol = 2*pi*dl*((data.p_rad(k)+dr/2)^2-(data.p_rad(k)-dr/2)^2); % vol of cell in 2D
% vol = 4/3*pi*dl*((data.p_rad(k)+dr/2)^3-(data.p_rad(k)-dr/2)^3); % vol of cell
if (cr_vol(i)>0);%vol)
cr_vol(i) = cr_vol(i)-vol;
cr_depth(i) = cr_depth(i) + 1; % number of cells from top downwards filled as crust
else
break
end
end
end
fprintf('Crust is between %d and %d cells thick\n',min(cr_depth),max(cr_depth))
crust = zeros(size(data.T));
for i=1:size(data.T,1) % in 3D would need do add for-loop for size(data.T,2)
crust(i,1,nr-cr_depth(i):nr)=max(data.cr_age)-data.cr_age(i);
end
produce_plot('Cr2D_',sprintf(strcat("Crust age at t=",num2str(data.time/data.t_yr/1000000.0,"%04.0f"),"Myr")),crust,input,data,res,ar,50,'bgry',min_z=-1,max_z=-1,full=fullV)
endif
if (reduc==0)
if(part==1)
produce_plot('OP2D_',"Tracer",data.output_tracer,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV,partX=A(:,1),partY=A(:,2),partV=A(:,21))
%% produce_plot('PD2D_',"Tracer",data.output_tracer,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV,partX=A(:,1),partY=A(:,2),partV=A(:,17))
......@@ -60,7 +92,6 @@ function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
produce_plot('O2D_',"Tracer",data.output_tracer,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)
endif
if (reduc==0)
%min_z=-max(abs(min(min(data.vl))),max(max(data.vl)));
%max_z=max(abs(min(min(data.vl))),max(max(data.vl)));
produce_plot('Vlat2D_',"Lateral velocity in cm/yr",data.vl,input,data,res,ar,10,'jet',min_z=-1,max_z=-1,full=fullV)
......@@ -75,15 +106,20 @@ function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
if (reduc==0)
produce_plot('De2D_',"Density",data.ref_r,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV)
produce_plot('Mrad2D_',"Radial mass transport in kg/(m^2 yr)",data.vr.*data.ref_r/100,input,data,res,ar,10,'jet',min_z=-1,max_z=-1,full=fullV)
elseif (reduc==1)
produce_plot('Mradp2D_',"Total radial mass transport in kg/(m^2 yr)",abs(data.vr.*data.ref_r/100),input,data,res,ar,10,'hot',min_z=-1,max_z=-1,full=fullV)
%elseif (reduc==1)
%produce_plot('Mradp2D_',"Total radial mass transport in kg/(m^2 yr)",abs(data.vr.*data.ref_r/100),input,data,res,ar,10,'hot',min_z=-1,max_z=-1,full=fullV)
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,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(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.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
% if(size(data.rec,1)>1), produce_plot('Rec2D_',"Recycling",data.rec,input,data,res,ar,50,'jet',min_z=3,max_z=7,full=fullV), endif %input.nr_ph
......@@ -91,14 +127,11 @@ function plot_snap(max_depl=0,max_T=0,fullV=-1,part=0,reduc=0,fin=' ')
if(size(data.cond,1)>1), produce_plot('K2D_',"Conductivity",data.cond,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
if(size(data.heat,1)>1), produce_plot('H2D_',"Heat souce redistribution",data.heat,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
%if(size(data.visc_diff,1)>1), produce_plot('ViDiff2D_',"Viscosity Diff in log",log10(data.visc_diff),input,data,res,ar,50,'jet',min_z=-1,max_z=-1,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,'hot',min_z=0,max_z=max_depl,full=fullV), endif
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.CO2,1)>1), produce_plot('CD2D_',"CO2 content",data.CO2,input,data,res,ar,50,'jet',min_z=0,max_z=500,full=fullV), endif %input.CO2), endif
if(size(data.U238,1)>1), produce_plot('HS_U2382D_',"U238 in ppm",data.U238,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
if(size(data.U235,1)>1), produce_plot('HS_U2352D_',"U235 in ppm",data.U235,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
if(size(data.Th232,1)>1), produce_plot('HS_Th2322D_',"Th232 in ppm",data.Th232,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
if(size(data.K40,1)>1), produce_plot('HS_K402D_',"K40 in ppm",data.K40,input,data,res,ar,50,'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.entropy,1)>1), produce_plot('En2D_',"Entropy",data.entropy,input,data,res,ar,50,'jet',min_z=-1,max_z=-1,full=fullV), endif
% lateral temperature variations
......@@ -208,7 +241,7 @@ function produce_plot(name,title_plot,field,input,data,res,ar,nr_contour,cmap,mi
max_z = max(max(field));
endif
font_title=round(30*res/600);
font_title=round(20*res/600);%(30*res/600);
font_axes_title=round(24*res/600);
font_axes=round(20*res/600);
......@@ -270,6 +303,7 @@ function produce_plot(name,title_plot,field,input,data,res,ar,nr_contour,cmap,mi
% contourf(data.p_xx_dim,data.p_zz_dim,dummy,nr_contour);
[C,h] = contourf(data.p_xx_dim,data.p_zz_dim,dummy,nr_contour);
set(h,'LineColor','none')
set(gca,'visible','off')
endif
% surface(data.p_xx_dim,data.p_zz_dim,dummy,'EdgeColor','none'); % works both for spherical and cartesian plots
......@@ -299,8 +333,8 @@ function produce_plot(name,title_plot,field,input,data,res,ar,nr_contour,cmap,mi
axis ([min(min(data.p_xx_dim)), max(max(data.p_xx_dim)), min(min(data.p_zz_dim)), max(max(data.p_zz_dim))]); % min_z max_z dazu?
dimen="";
if (data.D>1.0), dimen=" [km]";, endif
xlabel (sprintf(strcat("Lateral direction",dimen)),'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf','FontSize',font_axes_title);%dejavu/DejaVuSans.ttf','FontSize',font_axes_title);
ylabel (sprintf(strcat("Radius",dimen)),'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf','FontSize',font_axes_title);%dejavu/DejaVuSans.ttf','FontSize',font_axes_title);
% xlabel (sprintf(strcat("Lateral direction",dimen)),'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf','FontSize',font_axes_title);%dejavu/DejaVuSans.ttf','FontSize',font_axes_title);
% ylabel (sprintf(strcat("Radius",dimen)),'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Regular.ttf','FontSize',font_axes_title);%dejavu/DejaVuSans.ttf','FontSize',font_axes_title);
title (title_plot,'FontName','/usr/share/fonts/truetype/liberation/LiberationSansNarrow-Bold.ttf','FontSize',font_title);%dejavu/DejaVuSans-Bold.ttf','FontSize',font_title);
filename=sprintf(strcat('plots/',name,num2str(data.iter,"%07.0f"),'.png'),f);
......
......@@ -44,6 +44,7 @@ function input = read_input(fpath)
input.phi_factor = value(strcmp("phi_factor", data));
input.CO2 = value(strcmp("CO2_ini", data)); % in ppm
input.H2O = value(strcmp("mantle_w_H2O", data)); % in ppm
input.Chi_extr = value(strcmp("Chi_extr", data));
input.res = value(strcmp("resolution", data)); % resolution in m
input.read_profs = value(strcmp("read_profs", data));
input.nr_ph = value(strcmp("PhT_nr", data));
......
......@@ -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")
......@@ -91,7 +91,11 @@ function [data,input,program_dir,fname,fpath] = read_snap(finput=' ')
elseif(name == "dl"), hnl=nl+2;,hny=ny+2*n3D;,hnr=nr+2;,data.dl_dz=read_matrix(hnl,hny,hnr,data_snap);
elseif(name == "df"), hnl=nl+2;,hny=ny+2*n3D;,hnr=nr+2;,data.DeltaF=read_matrix(hnl,hny,hnr,data_snap);
elseif(name == "sp"), hnl=nl+2;,hny=ny+2*n3D;,hnr=nr+2;,data.serp=read_matrix(hnl,hny,hnr,data_snap);
elseif(name == "tH"), bin_b=fread(data_snap,1,'int32');, data.total_mass_H2O=fread(data_snap,1,'double');, bin_e=fread(data_snap,1,'int32');
elseif(name == "tH"), bin_b=fread(data_snap,1,'int32');, data.total_mass_H2O=fread(data_snap,1,'double');, bin_e=fread(data_snap,1,'int32'); % old version, now M1-M4
elseif(name == "M1"), bin_b=fread(data_snap,1,'int32');, data.total_mass_H2O=fread(data_snap,1,'double');, bin_e=fread(data_snap,1,'int32');
elseif(name == "M2"), bin_b=fread(data_snap,1,'int32');, data.total_mass_CO2=fread(data_snap,1,'double');, bin_e=fread(data_snap,1,'int32');
elseif(name == "M3"), bin_b=fread(data_snap,1,'int32');, data.total_mass_H2=fread(data_snap,1,'double');, bin_e=fread(data_snap,1,'int32');
elseif(name == "M4"), bin_b=fread(data_snap,1,'int32');, data.total_mass_CO=fread(data_snap,1,'double');, bin_e=fread(data_snap,1,'int32');
elseif(name == "oc"), hnl=nl+2;,hny=ny+2*n3D;,hnr=1;,data.oc_vol_tot=read_matrix(hnl,hny,hnr,data_snap); %2D matrix, no radial dependence
elseif(name == "fe"), hnl=nl+2;,hny=ny+2*n3D;,hnr=1;,data.fc_vol_tot=read_matrix(hnl,hny,hnr,data_snap); %2D matrix, no radial dependence
elseif(name == "cc"), hnl=nl+2;,hny=ny+2*n3D;,hnr=1;,data.cc_vol_tot=read_matrix(hnl,hny,hnr,data_snap); %2D matrix, no radial dependence
......
This diff is collapsed.
......@@ -51,15 +51,15 @@ O2_fugacity_IW=10.0**LogfO2_IW
LogfO2_CI=2.4976_dp+(-9.8605e+3_dp/T)+(-17.0701e+6_dp/(T**2.0_dp))+(7.522e+9_dp/(T**3.0_dp))+(-1.0404e+12_dp/(T**4.0_dp))
LogfO2_CV=9.0621_dp+(-31.193e+3_dp/T)+(5.1092e+6_dp/(T**2.0_dp))+(-1.8475e+9_dp/(T**3.0_dp))+(0.2e+12_dp/(T**4.0_dp))
LogfO2_H=5.0743_dp+(-22.906e+3_dp/T)+(-5.661e+6_dp/(T**2.0_dp))+(2.0634e+9_dp/(T**3.0_dp))+(-0.2618e+12_dp/(T**4.0_dp))
LogfO2_EH=4.9495_dp+(-24.024e+3_dp/T)+(-4.6236e+6_dp/(T**2.0_dp))+(1.7177e+9_dp/(T**3.0_dp))+(-0.2332e+12_dp/(T**4.0_dp))
LogfO2_Eucrite=5.4856_dp+(-25.127e+3_dp/T)+(-3.658e+6_dp/(T**2.0_dp))+(1.3014e+9_dp/(T**3.0_dp))+(-0.165e+12_dp/(T**4.0_dp))
!LogfO2_CI=2.4976_dp+(-9.8605e+3_dp/T)+(-17.0701e+6_dp/(T**2.0_dp))+(7.522e+9_dp/(T**3.0_dp))+(-1.0404e+12_dp/(T**4.0_dp))
!
!LogfO2_CV=9.0621_dp+(-31.193e+3_dp/T)+(5.1092e+6_dp/(T**2.0_dp))+(-1.8475e+9_dp/(T**3.0_dp))+(0.2e+12_dp/(T**4.0_dp))
!
!LogfO2_H=5.0743_dp+(-22.906e+3_dp/T)+(-5.661e+6_dp/(T**2.0_dp))+(2.0634e+9_dp/(T**3.0_dp))+(-0.2618e+12_dp/(T**4.0_dp))
!
!LogfO2_EH=4.9495_dp+(-24.024e+3_dp/T)+(-4.6236e+6_dp/(T**2.0_dp))+(1.7177e+9_dp/(T**3.0_dp))+(-0.2332e+12_dp/(T**4.0_dp))
!
!LogfO2_Eucrite=5.4856_dp+(-25.127e+3_dp/T)+(-3.658e+6_dp/(T**2.0_dp))+(1.3014e+9_dp/(T**3.0_dp))+(-0.165e+12_dp/(T**4.0_dp))
......
......@@ -485,7 +485,7 @@ module initialize
enddo
close(70)
do k=1,n_t
do k=2,n_t
if (k.lt.n_t) then
! either core or mantle
vol = 4.0_dp/3.0_dp*acos(-1.0_dp) * (vals_t(n_t+1-k,4)**3-vals_t(n_t-k,4)**3)
......
......@@ -117,11 +117,12 @@ 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,carbon_av,DeltaF_av,vol,vol_tot,heat_m
real(dp)::depl_av,depl_av_old,depl_p_av_old,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(:,:)
real(dp),allocatable::dl_dummy(:),Masses(:),XiFe(:),XiH2O(:)
real(dp)::frac_CO2,frac_CO,frac_H2O,frac_H2
integer,allocatable::melt_local(:,:,:)
real(dp)::Layer_0,Layer_1,LayerOld,LayerNext
! logical::finished,finished_OI,finished_II,outputVisc
......@@ -140,6 +141,7 @@ program CHIC
t = 0
big_number=10D99
time1=CLOCK() !omp_get_wtime()
depl_av=0.0_dp
call SetTimeParallel()
......@@ -692,6 +694,26 @@ program CHIC
mc%mass_CO2 = 0.0_dp
mc%total_mass_CO2 = 0.0_dp
! example output speciation code to see if limitations occur
! call spec_COH_Fegley2013(0.5_dp,0.5_dp,pX%fO2,1.0_dp,1600.0_dp,frac_CO2,frac_CO,frac_H2O,frac_H2)
! write(*,*) "0.5 0.5 ",pX%fO2," 1 1600: ",frac_H2O,"/",frac_CO2,"/",frac_H2,"/",frac_CO