Update atmosphere, correct CO2 calculation

Octave/make_movie.m

@@ -6,7 +6,7 @@ filenames = sprintf('%s%s', fpath, '/data_snap*')
 files = glob(filenames);
 for i=1:numel(files)
   [~, name] = fileparts(files{i});
-  eval(sprintf('plot_snap(30,4000,-1,0,1,"%s", "-ascii");', files{i}));
+  eval(sprintf('plot_snap(30,4000,-1,1,0,"%s", "-ascii");', files{i}));
 %  eval(sprintf('%s = load("%s", "-ascii");', name, files{i}));
 endfor
 

Octave/melt_output_for_speciation_lid.m

@@ -424,7 +424,8 @@ end
 %%%%%%%%%%%%%%%%%%%%%%
 mu_av_extr = H2_mole*mass_H2 + H2O_mole*mass_H2O + CO_mole*mass_CO + CO2_mole*mass_CO2;
 M_H2=zeros(1,nt);,M_H2O=zeros(1,nt);,M_CO=zeros(1,nt);,M_CO2=zeros(1,nt);
-x = X_H2O_mole~=0; %x = X_H2O_mole!=0; in octave
+x = X_H2O_mole~=0; % in octave
+%x = X_H2O_mole!=0; in Matlab
 M_H2(x)  = X_H2O_mantle/mass_H2O .* H2_mole(x) ./X_H2O_mole(x) * mass_H2 .* melt_vol(x) * 1e-6 * Chi_extr; % mass of H2 in kg
 M_H2O(x) = X_H2O_mantle/mass_H2O .* H2O_mole(x)./X_H2O_mole(x) * mass_H2O.* melt_vol(x) * 1e-6 * Chi_extr; % mass of H2O in kg
 M_CO(x)  = X_CO2_mantle/mass_CO2 .* CO_mole(x) ./X_CO2_mole(x) * mass_CO .* melt_vol(x) * 1e-6 * Chi_extr; % mass of H2 in kg

Octave/plot_snap.m

@@ -111,14 +111,14 @@ 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

src/atmosphere.f90

@@ -225,6 +225,12 @@ contains
 !            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)
 
@@ -242,12 +248,15 @@ contains
           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)
+!          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)=((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) * 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
+          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)
                
         end if