Commit 51b0d1b6 authored by Lena Noack's avatar Lena Noack
Browse files

Update atmosphere to include solubility treatment of gases (application for Venus)

parent 8440b503
......@@ -111,13 +111,13 @@ 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.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
......
......@@ -47,8 +47,8 @@ do
[H2O_dissolved,CO2_dissolved] = solubility(p_lid,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
if (H2O_fl_mole==1), dX_CO2=0;, CO2_dissolved=CO2_dissolved_old;, endif
if (H2O_fl_mole==0), dX_H2O=0;, H2O_dissolved=H2O_dissolved_old;, endif
if (H2O_fl_mole==1), dX_CO2=0;, CO2_dissolved=CO2_dissolved_old;, end
if (H2O_fl_mole==0), dX_H2O=0;, H2O_dissolved=H2O_dissolved_old;, end
if (dX_H2O+dX_CO2>0)&&(no_water==0)
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);
......@@ -57,7 +57,7 @@ do
else
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
end
% 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)
......@@ -78,7 +78,7 @@ end
graphics_toolkit gnuplot
f = figure('visible','off');
%f = figure('visible','off');
semilogx(10.^(ln_p),CO2_diss,'r-',10.^(ln_p),H2O_diss,'b-')
set (findall (gcf (), '-property', 'fontsize'), 'fontsize', 20);
......@@ -88,8 +88,8 @@ title(sprintf(strcat('Melt (ppm): X_{H2O}=',num2str(X_H2O*10000,"%4.0f"),'; X_{C
xlabel('Total pressure in bar','FontSize',28)
ylabel('Abundance of volatiles dissolved in melt','FontSize',28)
filename=sprintf(strcat('C:/OwnCloud_FUB/Literature_Projects/Eigene_Projekte/Project_Ortenzi_Dorn/Solubility/',...
'Solubility_H2O_',num2str(X_H2O*10000,'%4.0f'),'_CO2_',num2str(X_CO2,'%4.0f'),'.png'),f);
print(filename,'-dpng','-S800,800');
clear f;
%filename=sprintf(strcat('C:/OwnCloud_FUB/Literature_Projects/Eigene_Projekte/Project_Ortenzi_Dorn/Solubility/',...
% 'Solubility_H2O_',num2str(X_H2O*10000,'%4.0f'),'_CO2_',num2str(X_CO2,'%4.0f'),'.png'),f);
%print(filename,'-dpng','-S800,800');
%
%clear f;
This diff is collapsed.
......@@ -946,6 +946,9 @@ contains
case("melt_extract")
pM%melt_extract=int(val)
if (rank.eq.0) write(*,*) var,"=",pM%melt_extract
case("DeltaF_min") ! in %
pM%DeltaF_min=val/100.0_dp
if (rank.eq.0) write(*,*) var,"=",pM%DeltaF_min
case("add_crust_depth")
pM%add_crust_depth=val/pF%D
if (rank.eq.0) write(*,*) var,"=",val," -> ",pM%add_crust_depth
......@@ -3174,6 +3177,27 @@ contains
case("density_melt")
pX%density_melt = val / pF%rho
if (rank.eq.0) write(*,*) var,"=",pX%density_melt
case("solubility")
pX%solubility = int(val)
if (rank.eq.0) write(*,*) var,"=",pX%solubility
case("mass_H2O")
pX%mass_H2O = val/(pF%D**3 * pF%rho) ! from kg to nondim
if (rank.eq.0) write(*,*) var,"=",val," -> ",pX%mass_H2O
case("mass_H2")
pX%mass_H2 = val/(pF%D**3 * pF%rho) ! from kg to nondim
if (rank.eq.0) write(*,*) var,"=",val," -> ",pX%mass_H2
case("mass_CO2")
pX%mass_CO2 = val/(pF%D**3 * pF%rho) ! from kg to nondim
if (rank.eq.0) write(*,*) var,"=",val," -> ",pX%mass_CO2
case("mass_CO")
pX%mass_CO = val/(pF%D**3 * pF%rho) ! from kg to nondim
if (rank.eq.0) write(*,*) var,"=",val," -> ",pX%mass_CO
case("mass_N2")
pX%mass_N2 = val/(pF%D**3 * pF%rho) ! from kg to nondim
if (rank.eq.0) write(*,*) var,"=",val," -> ",pX%mass_N2
case("tau_ex_H2O")
pX%tau_ex_H2O = val*pF%time_yr ! residence time of water in atmosphere in years -> applied to H2O and H2 in the atmosphere
if (rank.eq.0) write(*,*) var,"=",val," -> ",pX%tau_ex_H2O
end select
read(42,*) var,equal,val
end do
......@@ -3287,6 +3311,9 @@ contains
case("CO2_ini") ! initial CO2 concentration in ppm
pX%CO2_ini = val
if (rank.eq.0) write(*,*) var,"=",pX%CO2_ini
case("mass_CO2")
pX%mass_CO2 = val/(pF%D**3 * pF%rho) ! from kg to nondim
if (rank.eq.0) write(*,*) var,"=",val," -> ",pX%mass_CO2
case("fO2") ! oxygen fugacity
pX%fO2 = val
if (rank.eq.0) write(*,*) var,"=",pX%fO2
......@@ -3306,6 +3333,9 @@ contains
pM%melt_depth=val
endif
if (rank.eq.0) write(*,*) var,"=",val," -> ",pM%melt_depth
case("DeltaF_min") ! in %
pM%DeltaF_min=val/100.0_dp
if (rank.eq.0) write(*,*) var,"=",pM%DeltaF_min
end select
read(44,*) var,equal,val
end do
......
......@@ -117,7 +117,7 @@ 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,depl_av_old,depl_p_av_old,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,heat_av
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(:,:)
......@@ -686,13 +686,11 @@ program CHIC
&pN%p_average,pN%p_average2,pN%pr_time,pP%moveRef) ! transfer particle properties to cells
if (pI%compress.ge.4) call UpdatePressure(jmin_2D,cell%ref_p,cell%ref_g,cell%ref_r,mesh%rc,pE%Gr,pE%Psurf,pF) ! Get real pressure profile depending on actual density and gravity fields
if (pX%use_atmos.ne.0) call init_atmos(pE,pF,pI,pX,mc,mesh%rmax)
endif
mc%mass_H2O = 0.0_dp
mc%total_mass_H2O = 0.0_dp
mc%mass_CO2 = 0.0_dp
mc%total_mass_CO2 = 0.0_dp
! example output speciation code to see if limitations occur
......@@ -961,17 +959,28 @@ program CHIC
if ((pN%debug.ge.1).and.(rank.eq.0)) write(*,*) "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
else
if(rank.eq.0) then
heat_m = 0.0_dp
vol_tot = 0.0_dp
if (pX%dU238.eq.0.0_dp) then
heat_m = radiogenic_heating(pI%H0,pI%lambda,t,pI)*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12
heat_av = heat_m
else
! initial value as it should be on all cells
do ind_i=1,mesh%nl
do ind_k=1,mesh%nr
cell%h(ind_i,1,ind_k)=radiogenic_nuclei(cell%X_U238(ind_i,1,ind_k),cell%X_U235(ind_i,1,ind_k),&
&cell%X_Th232(ind_i,1,ind_k),cell%X_K40(ind_i,1,ind_k),t,pI)
heat_m = heat_m + cell%h(ind_i,1,ind_k) * mesh%dVi(ind_i,1,ind_k)*cell%ref_r(ind_i,1,ind_k)
vol_tot = vol_tot + mesh%dVi(ind_i,1,ind_k)*cell%ref_r(ind_i,1,ind_k)
enddo
enddo
heat_m = sum(cell%h(1:mesh%nl,1,1:mesh%nr))/(mesh%nl*mesh%nr)*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12 ! ToDo: vol-averaged instead of radius-averaged?
!heat_m = sum(cell%h(1:mesh%nl,1,1:mesh%nr))/(mesh%nl*mesh%nr)*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12 ! ToDo: vol-averaged instead of radius-averaged?
heat_m = heat_m/vol_tot*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12
heat_av = radiogenic_nuclei(pI%C_U*0.9928e-9_dp*exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12U238),&
&pI%C_U*0.0071e-9_dp*exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12U235),&
&pI%C_Th*1.0e-9_dp *exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12Th232),&
&pI%C_K*0.000128e-6_dp*exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12K40),t,pI)&
& * pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12
endif
! str_out1 = "D(VD,W)="
! str_out2 = "%, "
......@@ -989,16 +998,16 @@ program CHIC
if (pH%use_magn_H.lt.2.0_dp) then
write(*,'("TS ",i6," IO",i3," II",i3," t=",ES9.2," dt=",ES9.2," v=",F10.5," v0=",F8.5,"/",F8.5," T=",F8.2,&
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," D(VD,W)=",F6.2,"%, "," H_av=",F10.7," pW/kg, ",F6.2,"s")') &
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," D(VD,W)=",F6.2,"%, "," H_av=",F10.7," (",F10.7,") pW/kg, ",F6.2,"s")') &
&i,i_outer,i_inner+i_inner_v,t/pF%time_yr,dt/pF%time_yr,Vrms*(100.0_dp*pF%D*pF%time_yr),V0rms*(100.0_dp*pF%D*pF%time_yr),&
&V0rms_max*(100.0_dp*pF%D*pF%time_yr),T_mean*pF%DeltaT+pF%Ts,Nu_t*pF%k*pF%DeltaT/pF%D,0.0_dp,&!Nu_b*pF%k*pF%DeltaT/pF%D,&
&ViscCon,ErrVD*100_dp,heat_m,0.0_dp
&ViscCon,ErrVD*100_dp,heat_m,heat_av,0.0_dp
else
write(*,'("TS ",i6," IO",i3," II",i3," t=",ES9.2," dt=",ES9.2," v=",F10.5," v0=",F8.5,"/",F8.5," T=",F8.2,&
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," Q_eddy=",ES10.3," pW/kg, "," H_av=",F10.7," pW/kg, ",F6.2,"s")') &
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," Q_eddy=",ES10.3," pW/kg, "," H_av=",F10.7," (",F10.7,") pW/kg, ",F6.2,"s")') &
&i,i_outer,i_inner+i_inner_v,t/pF%time_yr,dt/pF%time_yr,Vrms*(100.0_dp*pF%D*pF%time_yr),V0rms*(100.0_dp*pF%D*pF%time_yr),&
&V0rms_max*(100.0_dp*pF%D*pF%time_yr),T_mean*pF%DeltaT+pF%Ts,Nu_t*pF%k*pF%DeltaT/pF%D,0.0_dp,&!Nu_b*pF%k*pF%DeltaT/pF%D,&
&ViscCon,pH%Q_eddy*(pF%k*pF%DeltaT)/(pF%D**2*pF%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,heat_m,0.0_dp
&ViscCon,pH%Q_eddy*(pF%k*pF%DeltaT)/(pF%D**2*pF%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,heat_m,heat_av,0.0_dp
endif
endif
endif
......@@ -1275,23 +1284,34 @@ program CHIC
endif
ViscCon = maxval(eta%CV)/minval(eta%CV)
heat_m = 0.0_dp
vol_tot = 0.0_dp
if (pX%dU238.eq.0.0_dp) then
heat_m = radiogenic_heating(pI%H0,pI%lambda,t,pI)*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12
heat_av = heat_m
else
do ind_i=1,mesh%nl
do ind_k=1,mesh%nr
cell%h(ind_i,1,ind_k)=radiogenic_nuclei(cell%X_U238(ind_i,1,ind_k),cell%X_U235(ind_i,1,ind_k),&
&cell%X_Th232(ind_i,1,ind_k),cell%X_K40(ind_i,1,ind_k),t,pI)
heat_m = heat_m + cell%h(ind_i,1,ind_k) * mesh%dVi(ind_i,1,ind_k)*cell%ref_r(ind_i,1,ind_k)
vol_tot = vol_tot + mesh%dVi(ind_i,1,ind_k)*cell%ref_r(ind_i,1,ind_k)
enddo
enddo
heat_m = sum(cell%h(1:mesh%nl,1,1:mesh%nr))/(mesh%nl*mesh%nr)*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12 ! ToDo: vol-averaged instead of radius-averaged?
!heat_m = sum(cell%h(1:mesh%nl,1,1:mesh%nr))/(mesh%nl*mesh%nr)*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12 ! ToDo: vol-averaged instead of radius-averaged?
heat_m = heat_m/vol_tot*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12
heat_av = radiogenic_nuclei(pI%C_U*0.9928e-9_dp*exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12U238),&
&pI%C_U*0.0071e-9_dp*exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12U235),&
&pI%C_Th*1.0e-9_dp *exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12Th232),&
&pI%C_K*0.000128e-6_dp*exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12K40),t,pI)&
& * pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12
endif
if (pE%output_dim.eq.0) then
if ((pN%debug.ge.1).and.(rank.eq.0)) write(*,*) "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
if(rank.eq.0)write(*,'("TS ",i6," IO",i3," II",i3," t=",F9.6," dt=",ES9.2," v=",F12.4," v0=",F8.3,"/",F8.3," T=",F8.5," Te="&
&,F7.4,"/",F7.4," Nu=",F8.4,"/",F8.4,"/",F8.4,"/",F8.4,", VC=",ES10.3," VD=",F7.4," W=",F7.4," E",F6.2,"% ",F6.2,"s")') &
&i,i_outer,i_inner+i_inner_v,t,dt,Vrms,V0rms,V0rms_max,T_mean,Te_save,Te2,Nu,Nu_t,Nu_m,Nu_b,ViscCon,VD,Work,ErrVD*100_dp,0.0_dp
&i,i_outer,i_inner+i_inner_v,t,dt,Vrms,V0rms,V0rms_max,T_mean,Te_save,Te2,Nu,Nu_t,Nu_m,Nu_b,ViscCon,VD,Work,ErrVD*100_dp,0.0_dp
if ((pN%debug.ge.1).and.(rank.eq.0)) write(*,*) "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
else
! if(rank.eq.0) write(*,'("TS ",i6," IO",i3," II",i3," t=",ES9.2," dt=",ES9.2," v=",F10.5," v0=",F8.5,"/",F8.5," T=",F8.2,&
......@@ -1301,16 +1321,16 @@ program CHIC
! &ViscCon,ErrVD*100_dp,0.0_dp
if (pH%use_magn_H.lt.2.0_dp) then
if(rank.eq.0) write(*,'("TS ",i6," IO",i3," II",i3," t=",ES9.2," dt=",ES9.2," v=",F10.5," v0=",F8.5,"/",F8.5," T=",F8.2,&
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," D(VD,W)=",F6.2,"%, "," H_av=",F10.7," pW/kg, ",F6.2,"s")') &
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," D(VD,W)=",F6.2,"%, "," H_av=",F10.7," (",F10.7,") pW/kg, ",F6.2,"s")') &
&i,i_outer,i_inner+i_inner_v,t/pF%time_yr,dt/pF%time_yr,Vrms*(100.0_dp*pF%D*pF%time_yr),V0rms*(100.0_dp*pF%D*pF%time_yr),&
&V0rms_max*(100.0_dp*pF%D*pF%time_yr),T_mean*pF%DeltaT+pF%Ts,Nu_t*pF%k*pF%DeltaT/pF%D,0.0_dp,&!Nu_b*pF%k*pF%DeltaT/pF%D,&
&ViscCon,ErrVD*100_dp,heat_m,0.0_dp
&ViscCon,ErrVD*100_dp,heat_m,heat_av,0.0_dp
else
if(rank.eq.0) write(*,'("TS ",i6," IO",i3," II",i3," t=",ES9.2," dt=",ES9.2," v=",F10.5," v0=",F8.5,"/",F8.5," T=",F8.2,&
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," Q_eddy=",ES10.3," pW/kg, "," H_av=",F10.7," pW/kg, ",F6.2,"s")') &
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," Q_eddy=",ES10.3," pW/kg, "," H_av=",F10.7," (",F10.7,") pW/kg, ",F6.2,"s")') &
&i,i_outer,i_inner+i_inner_v,t/pF%time_yr,dt/pF%time_yr,Vrms*(100.0_dp*pF%D*pF%time_yr),V0rms*(100.0_dp*pF%D*pF%time_yr),&
&V0rms_max*(100.0_dp*pF%D*pF%time_yr),T_mean*pF%DeltaT+pF%Ts,Nu_t*pF%k*pF%DeltaT/pF%D,0.0_dp,&!Nu_b*pF%k*pF%DeltaT/pF%D,&
&ViscCon,pH%Q_eddy*(pF%k*pF%DeltaT)/(pF%D**2*pF%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,heat_m,0.0_dp
&ViscCon,pH%Q_eddy*(pF%k*pF%DeltaT)/(pF%D**2*pF%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,heat_m,heat_av,0.0_dp
endif
endif
......@@ -1917,23 +1937,34 @@ program CHIC
&ErrVD*100_dp,timeS2-timeS1
else
if(rank.eq.0) then
heat_m = 0.0_dp
vol_tot = 0.0_dp
if (pX%dU238.eq.0.0_dp) then
heat_m = radiogenic_heating(pI%H0,pI%lambda,t,pI)*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12
heat_av = heat_m
else
do ind_i=1,mesh%nl
do ind_k=1,mesh%nr
cell%h(ind_i,1,ind_k)=radiogenic_nuclei(cell%X_U238(ind_i,1,ind_k),cell%X_U235(ind_i,1,ind_k),&
&cell%X_Th232(ind_i,1,ind_k),cell%X_K40(ind_i,1,ind_k),t,pI)
heat_m = heat_m + cell%h(ind_i,1,ind_k) * mesh%dVi(ind_i,1,ind_k)*cell%ref_r(ind_i,1,ind_k)
vol_tot = vol_tot + mesh%dVi(ind_i,1,ind_k)*cell%ref_r(ind_i,1,ind_k)
enddo
enddo
heat_m = sum(cell%h(1:mesh%nl,1,1:mesh%nr))/(mesh%nl*mesh%nr)*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12
!heat_m = sum(cell%h(1:mesh%nl,1,1:mesh%nr))/(mesh%nl*mesh%nr)*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12
heat_m = heat_m/vol_tot*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12
heat_av = radiogenic_nuclei(pI%C_U*0.9928e-9_dp*exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12U238),&
&pI%C_U*0.0071e-9_dp*exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12U235),&
&pI%C_Th*1.0e-9_dp *exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12Th232),&
&pI%C_K*0.000128e-6_dp*exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12K40),t,pI)&
& * pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12
endif
write(*,'("TS ",i6," IO",i3," II",i3," t=",ES9.2," dt=",ES9.2," v=",F10.5," v0=",F8.5,"/",F8.5," T=",F8.2,&
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," D(VD,W)=",F6.2,"% "," H_av=",F10.7," pW/kg, ",F6.2,"s")') i,i_outer,&
&i_inner+i_inner_v,t/pF%time_yr,dt/pF%time_yr,Vrms*(100.0_dp*pF%D*pF%time_yr),V0rms*(100.0_dp*pF%D*pF%time_yr),&
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," D(VD,W)=",F6.2,"% "," H_av=",F10.7," (",F10.7,") pW/kg, ",F6.2,"s")') i,&
&i_outer,i_inner+i_inner_v,t/pF%time_yr,dt/pF%time_yr,Vrms*(100.0_dp*pF%D*pF%time_yr),V0rms*(100.0_dp*pF%D*pF%time_yr),&
&V0rms_max*(100.0_dp*pF%D*pF%time_yr),T_mean*pF%DeltaT+pF%Ts,Nu_t*pF%k*pF%DeltaT/pF%D,Nu_b*pF%k*pF%DeltaT/pF%D,&
&ViscCon,ErrVD*100_dp,heat_m,timeS2-timeS1
&ViscCon,ErrVD*100_dp,heat_m,heat_av,timeS2-timeS1
endif
endif
if ((pN%debug.ge.2).and.(rank.eq.0)) write(*,*) "--------------------------------------------------------"
......@@ -2093,6 +2124,8 @@ program CHIC
if (pX%use_atmos.ge.1) then
call get_atmos(t,dt,pF,pG,pI,pM,pN,pP,pX,field,depl_av,T_mean,part,mesh,cell,mc,pG%dim,pE%Psurf)
pI%Psurf = pE%Psurf * (pF%D*pF%rho*pF%g) * pF%D**2/(pF%kappa*pF%eta_ref) ! from lithostatic pressure to stress as used e.g. for pore pressure in plasticity.f90
pM%Psurf = pE%Psurf ! used in melt.f90 -> no need to use pE as additional input parameter then
endif
......@@ -2133,30 +2166,42 @@ program CHIC
if ((pN%debug.ge.1).and.(rank.eq.0)) write(*,*) "+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
else
if(rank.eq.0) then
heat_m = 0.0_dp
vol_tot = 0.0_dp
if (pX%dU238.eq.0.0_dp) then
heat_m = radiogenic_heating(pI%H0,pI%lambda,t,pI)*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12 ! H = H'*k*DeltaT/D² -> W/m³ -> divided by rho -> W/kg
heat_av = heat_m
else
do ind_i=1,mesh%nl
do ind_k=1,mesh%nr
cell%h(ind_i,1,ind_k)=radiogenic_nuclei(cell%X_U238(ind_i,1,ind_k),cell%X_U235(ind_i,1,ind_k),&
&cell%X_Th232(ind_i,1,ind_k),cell%X_K40(ind_i,1,ind_k),t,pI)
heat_m = heat_m + cell%h(ind_i,1,ind_k) * mesh%dVi(ind_i,1,ind_k)*cell%ref_r(ind_i,1,ind_k)
vol_tot = vol_tot + mesh%dVi(ind_i,1,ind_k)*cell%ref_r(ind_i,1,ind_k)
enddo
enddo
heat_m = sum(cell%h(1:mesh%nl,1,1:mesh%nr))/(mesh%nl*mesh%nr)*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12
!heat_m = sum(cell%h(1:mesh%nl,1,1:mesh%nr))/(mesh%nl*mesh%nr)*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12
heat_m = heat_m/vol_tot*pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12
heat_av = radiogenic_nuclei(pI%C_U*0.9928e-9_dp*exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12U238),&
&pI%C_U*0.0071e-9_dp*exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12U235),&
&pI%C_Th*1.0e-9_dp *exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12Th232),&
&pI%C_K*0.000128e-6_dp*exp(Log(2.0_dp)*(4.5e+9*pF%time_yr)/pI%t12K40),t,pI)&
& * pF%k*pF%DeltaT/pF%D**2/pF%rho*1.0e+12
endif
if (pH%use_magn_H.lt.2.0_dp) then
if(rank.eq.0) write(*,'("TS ",i6," IO",i3," II",i3," t=",ES9.2," dt=",ES9.2," v=",F10.5," v0=",F8.5,"/",F8.5," T=",F8.2,&
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," D(VD,W)=",F6.2,"%, "," H_av=",F10.7," pW/kg, ",F6.2,"s")') &
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," D(VD,W)=",F6.2,"%, "," H_av=",F10.7," (",F10.7,") pW/kg, ",F6.2,"s")') &
&i,i_outer,i_inner+i_inner_v,t/pF%time_yr,dt/pF%time_yr,Vrms*(100.0_dp*pF%D*pF%time_yr),V0rms*(100.0_dp*pF%D*pF%time_yr),&
&V0rms_max*(100.0_dp*pF%D*pF%time_yr),T_mean*pF%DeltaT+pF%Ts,Nu_t*pF%k*pF%DeltaT/pF%D,0.0_dp,&!Nu_b*pF%k*pF%DeltaT/pF%D,&
&ViscCon,ErrVD*100_dp,heat_m,timeS2-timeS1
&ViscCon,ErrVD*100_dp,heat_m,heat_av,timeS2-timeS1
else
if(rank.eq.0) write(*,'("TS ",i6," IO",i3," II",i3," t=",ES9.2," dt=",ES9.2," v=",F10.5," v0=",F8.5,"/",F8.5," T=",F8.2,&
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," Q_eddy=",ES10.3," pW/kg, "," H_av=",F10.7," pW/kg, ",F6.2,"s")') &
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," Q_eddy=",ES10.3," pW/kg, "," H_av=",F10.7," (",F10.7,") pW/kg, ",F6.2,"s")') &
&i,i_outer,i_inner+i_inner_v,t/pF%time_yr,dt/pF%time_yr,Vrms*(100.0_dp*pF%D*pF%time_yr),V0rms*(100.0_dp*pF%D*pF%time_yr),&
&V0rms_max*(100.0_dp*pF%D*pF%time_yr),T_mean*pF%DeltaT+pF%Ts,Nu_t*pF%k*pF%DeltaT/pF%D,0.0_dp,&!Nu_b*pF%k*pF%DeltaT/pF%D,&
&ViscCon,pH%Q_eddy*(pF%k*pF%DeltaT)/(pF%D**2*pF%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,heat_m,timeS2-timeS1
&ViscCon,pH%Q_eddy*(pF%k*pF%DeltaT)/(pF%D**2*pF%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,heat_m,heat_av,&
&timeS2-timeS1
endif
! write(*,'("TS ",i6," IO",i3," II",i3," t=",ES9.2," dt=",ES9.2," v=",F10.5," v0=",F8.5,"/",F8.5," T=",F8.2,&
! &" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," D(VD,W)=",F6.2,"% "," H_av=",F10.7," pW/kg, ",F6.2,"s")') &
......
......@@ -156,6 +156,7 @@ contains
!write(*,*) "Melting occurred",i,k,field%T(i,j,k),solidus,mc%liquidus(i,j,k),part%df(m),part%d(m)
!write(*,*) "Solidus=",solidus,"=",mc%solidus_ref(i,j,k),"+",part%s(m),"+",&
! &-solidus_hydrated(pM%sol_hyd,pM%sol_hyd_exp,part%w(m))
else if (pM%melt_extract.ne.1) then
part%df(m) = (field%T(i,j,k)-solidus_ref)/(mc%liquidus(i,j,k)-solidus_ref) !mc%solidus_ref(i,j,k)) (see explanation above for hydrated case) ! CHECK if T-solidus_ref or T-solidus ?
else
......@@ -166,6 +167,8 @@ contains
part%df(m) = 0.0_dp
endif
if ((pM%DeltaF_min.gt.0.0_dp).and.(part%df(m).lt.pM%DeltaF_min)) part%df(m) = 0.0_dp ! extract melt only above specific threshold
!if ((i.eq.1).and.(j.eq.1).and.(k.eq.mesh%nr)) &
! & write(*,*) "Melt p 2:",part%df(m),part%w(m)
......@@ -649,7 +652,7 @@ contains
!! !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! call in main.f90 before particles are set
subroutine set_init_crust(mc,cell,mesh,pC,pI,pM,pX)
type(melt_properties),intent(inout)::mc
type(cell_properties),intent(inout)::cell
......@@ -689,16 +692,26 @@ contains
if (cell%fld_f) cell%fc(:,:,k) = pC%cont_OC_fr_coh
if (cell%fld_a) cell%ref_a(:,:,k) = pC%cont_OC_alpha
if (mesh%rb(k).gt.(mesh%rmax-pC%crust_ini)) then ! whole cell in crust
mc%oc_vol_tot(:,:) = mc%oc_vol_tot(:,:) + mesh%dVi(:,:,k)
mc%oc_vol_tot(:,:) = mc%oc_vol_tot(:,:) + mesh%dVi(:,:,k)*cell%ref_r(1,1,k)
else
mc%oc_vol_tot(:,:) = mc%oc_vol_tot(:,:) + mesh%dVi(:,:,k)*(mesh%rb(k+1)-(mesh%rmax-pC%crust_ini))/mesh%dr_vec(k)
mc%oc_vol_tot(:,:) = mc%oc_vol_tot(:,:) + &
&mesh%dVi(:,:,k)*cell%ref_r(1,1,k)*(mesh%rb(k+1)-(mesh%rmax-pC%crust_ini))/mesh%dr_vec(k)
endif
else if (mesh%rb(k+1).ge.(mesh%rmax-pC%crust_ini)) then !upper part of cell still in crust
! main part of cell in mantle, do not change properties of cell
mc%oc_vol_tot(:,:) = mc%oc_vol_tot(:,:) + mesh%dVi(:,:,k)*(mesh%rb(k+1)-(mesh%rmax-pC%crust_ini))/mesh%dr_vec(k)
mc%oc_vol_tot(:,:) = mc%oc_vol_tot(:,:) + &
&mesh%dVi(:,:,k)*cell%ref_r(1,1,k)*(mesh%rb(k+1)-(mesh%rmax-pC%crust_ini))/mesh%dr_vec(k)
endif
enddo
! Variante 1 -> trace average crust ppm of trace elements
!mc%X_U238_add = cell%X_U238(1,1,1)*pM%Lambda
!mc%X_U235_add = cell%X_U235(1,1,1)*pM%Lambda
!mc%X_Th232_add = cell%X_Th232(1,1,1)*pM%Lambda
!mc%X_K40_add = cell%X_K40(1,1,1)*pM%Lambda
write(*,*) "+++Crust+++",pM%Lambda,cell%X_U238(1,1,1),cell%X_U238(1,1,mesh%nr)
end subroutine
......@@ -738,10 +751,11 @@ contains
do i=1,nl
do j=1,ny
crust_old = mc%oc_vol_tot(i,j) + mc%fc_vol_tot(i,j) + mc%cc_vol_tot(i,j)
!mc%X_U238_add(i,j) = mc%X_U238_add(i,j) * crust_old * cell%ref_r(i,j,nr) ! from ppm fraction back to volume of heat sources / mass
!mc%X_U235_add(i,j) = mc%X_U235_add(i,j) * crust_old * cell%ref_r(i,j,nr)
!mc%X_Th232_add(i,j) = mc%X_Th232_add(i,j) * crust_old * cell%ref_r(i,j,nr)
!mc%X_K40_add(i,j) = mc%X_K40_add(i,j) * crust_old * cell%ref_r(i,j,nr)
! Version 1
!mc%X_U238_add(i,j) = mc%X_U238_add(i,j) * crust_old !* cell%ref_r(i,j,nr) ! from ppm fraction back to volume of heat sources / mass
!mc%X_U235_add(i,j) = mc%X_U235_add(i,j) * crust_old !* cell%ref_r(i,j,nr)
!mc%X_Th232_add(i,j) = mc%X_Th232_add(i,j) * crust_old !* cell%ref_r(i,j,nr)
!mc%X_K40_add(i,j) = mc%X_K40_add(i,j) * crust_old !* cell%ref_r(i,j,nr)
do k=1,nr
if (mc%DeltaF(i,j,k).ne.0.0_dp) then ! bring new crust to surface
mc%cr_age(i,j) = t ! always assume that newest crust goes (partly) to the surface
......@@ -810,22 +824,22 @@ contains
if ((pX%dU238.ne.0.0_dp).and.(cell%r(i,j,k).lt.0.5_dp).and.(mc%DeltaF(i,j,k).gt.0.0_dp)) then ! is (in average) primordial mantle material or wedge material (should be 0 if primordial mantle)
mc%X_U238_add(i,j) = mc%X_U238_add(i,j) + cell%ref_r(i,j,k) * mesh%dVi(i,j,k) * cell%X_U238(i,j,k) * &
mc%X_U238_add(i,j) = mc%X_U238_add(i,j) + compr_fac * mesh%dVi(i,j,k) * cell%X_U238(i,j,k) * &
& (1.0_dp - (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dU238))
mc%X_U235_add(i,j) = mc%X_U235_add(i,j) + cell%ref_r(i,j,k) * mesh%dVi(i,j,k) * cell%X_U235(i,j,k) * &
mc%X_U235_add(i,j) = mc%X_U235_add(i,j) + compr_fac * mesh%dVi(i,j,k) * cell%X_U235(i,j,k) * &
& (1.0_dp - (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dU235))
mc%X_Th232_add(i,j) = mc%X_Th232_add(i,j) + cell%ref_r(i,j,k) * mesh%dVi(i,j,k) * cell%X_Th232(i,j,k) * &
mc%X_Th232_add(i,j) = mc%X_Th232_add(i,j) + compr_fac * mesh%dVi(i,j,k) * cell%X_Th232(i,j,k) * &
& (1.0_dp - (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dTh232))
mc%X_K40_add(i,j) = mc%X_K40_add(i,j) + cell%ref_r(i,j,k) * mesh%dVi(i,j,k) * cell%X_K40(i,j,k) * &
mc%X_K40_add(i,j) = mc%X_K40_add(i,j) + compr_fac * mesh%dVi(i,j,k) * cell%X_K40(i,j,k) * &
& (1.0_dp - (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dK40))
elseif ((pX%dU238.ne.0.0_dp).and.(mc%DeltaF(i,j,k).gt.0.0_dp)) then ! part of cell in crust, hence cannot use cell radiogenic heat sources -> use cell value below
mc%X_U238_add(i,j) = mc%X_U238_add(i,j) + cell%ref_r(i,j,k) * mesh%dVi(i,j,k) * cell%X_U238(i,j,k-1) * &
mc%X_U238_add(i,j) = mc%X_U238_add(i,j) + compr_fac * mesh%dVi(i,j,k) * cell%X_U238(i,j,k-1) * &
& (1.0_dp - (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dU238))
mc%X_U235_add(i,j) = mc%X_U235_add(i,j) + cell%ref_r(i,j,k) * mesh%dVi(i,j,k) * cell%X_U235(i,j,k-1) * &
mc%X_U235_add(i,j) = mc%X_U235_add(i,j) + compr_fac * mesh%dVi(i,j,k) * cell%X_U235(i,j,k-1) * &
& (1.0_dp - (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dU235))
mc%X_Th232_add(i,j) = mc%X_Th232_add(i,j) + cell%ref_r(i,j,k) * mesh%dVi(i,j,k) * cell%X_Th232(i,j,k-1) * &
mc%X_Th232_add(i,j) = mc%X_Th232_add(i,j) + compr_fac * mesh%dVi(i,j,k) * cell%X_Th232(i,j,k-1) * &
& (1.0_dp - (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dTh232))
mc%X_K40_add(i,j) = mc%X_K40_add(i,j) + cell%ref_r(i,j,k) * mesh%dVi(i,j,k) * cell%X_K40(i,j,k-1) * &
mc%X_K40_add(i,j) = mc%X_K40_add(i,j) + compr_fac * mesh%dVi(i,j,k) * cell%X_K40(i,j,k-1) * &
& (1.0_dp - (1.0_dp - mc%DeltaF(i,j,k))**(1.0_dp/pX%dK40))
endif
......@@ -886,8 +900,8 @@ contains
!mc%X_U235_add(i,j) = ( mc%X_U235_add(i,j) + X_U235_old*crust_old ) / crust_new
!mc%X_Th232_add(i,j) = ( mc%X_Th232_add(i,j) + X_Th232_old*crust_old ) / crust_new
!mc%X_K40_add(i,j) = ( mc%X_K40_add(i,j) + X_K40_old*crust_old ) / crust_new
! Version 2 -> add same amount everywhere, hence older crust more enriched then newer crust
mc%X_U238_add(i,j) = mc%X_U238_add(i,j) / crust_new ! what will be added to particles in new crust
! Version 1 new & 2 -> add same amount everywhere, hence older crust more enriched then newer crust
mc%X_U238_add(i,j) = mc%X_U238_add(i,j) / crust_new ! what will be added to particles in new crust ; in non-dimensional kg/m³
mc%X_U235_add(i,j) = mc%X_U235_add(i,j) / crust_new
mc%X_Th232_add(i,j) = mc%X_Th232_add(i,j) / crust_new
mc%X_K40_add(i,j) = mc%X_K40_add(i,j) / crust_new
......@@ -1350,13 +1364,14 @@ contains
k = part%cell_r(m)
! if recycling value has been reset above, then all particles in the cell get the same, new information
if (cell%r(i,j,k).lt.0.0_dp) then
!if (cell%r(i,j,k).lt.0.0_dp) then
if (part%pos_r(m).gt.mesh%rmax-depth_r_OC(i,j)) then
!if (m.eq.19) write(*,*) "Particle new crust: ",m,i,j,k,part%r(m),cell%r(i,j,k)
if (cell%r(i,j,k).lt.-15.0_dp) then ! residual of molten material, filled up with surrounding material
part%r(m) = 0.0_dp
else
part%r(m) = -cell%r(i,j,k)
endif
!if (cell%r(i,j,k).lt.-15.0_dp) then ! residual of molten material, filled up with surrounding material
! part%r(m) = 0.0_dp
!else
part%r(m) = 7.0_dp !-cell%r(i,j,k)
!endif
!if (m.eq.19) write(*,*) "Particle new crust: ",m,i,j,k,part%r(m),cell%r(i,j,k)
if (cell%fld_c) part%c(m) = cell%c(i,j,k)
if (cell%fld_t) part%t(m) = cell%ref_T(i,j,k)
......@@ -1374,16 +1389,16 @@ contains
if (cell%fld_f) part%fa(m) = cell%fa(i,j,k)
if (cell%fld_f) part%fc(m) = cell%fc(i,j,k)
if (cell%fld_a) part%a(m) = cell%ref_a(i,j,k)
endif
else !endif
!else if ((part%df(m).gt.0.0_dp)) then !.and.(part%r(m).lt.0.5_dp)) then
! ! then local depletion in mantle, reduce particle values
if (pX%dU238.ne.0.0_dp) then
! V1: X_new = X_old * (1-Xl) = X_old * (1-dF)^(1/D)
! V2: X_new (C_new) = X_old * (1-dF)^(1/D) / (1-dF)
!part%X_U238(m) = part%X_U238(m) * (1.0_dp - part%df(m))**(1.0_dp/pX%dU238) !/ (1.0_dp - part%df(m)) ! /(1-dF) to fill up whole cell again with same concentration
!part%X_U235(m) = part%X_U235(m) * (1.0_dp - part%df(m))**(1.0_dp/pX%dU235) !/ (1.0_dp - part%df(m))
!part%X_Th232(m) = part%X_Th232(m) * (1.0_dp - part%df(m))**(1.0_dp/pX%dTh232) !/ (1.0_dp - part%df(m))
!part%X_K40(m) = part%X_K40(m) * (1.0_dp - part%df(m))**(1.0_dp/pX%dK40) !/ (1.0_dp - part%df(m))
!part%X_U238(m) = part%X_U238(m) * (1.0_dp - part%df(m))**(1.0_dp/pX%dU238) !/ (1.0_dp-mc%DeltaF(i,j,k)) !/ (1.0_dp - part%df(m)) ! /(1-dF) to fill up whole cell again with same concentration
!part%X_U235(m) = part%X_U235(m) * (1.0_dp - part%df(m))**(1.0_dp/pX%dU235) !/ (1.0_dp-mc%DeltaF(i,j,k)) !/ (1.0_dp - part%df(m))
!part%X_Th232(m) = part%X_Th232(m) * (1.0_dp - part%df(m))**(1.0_dp/pX%dTh232) !/ (1.0_dp-mc%DeltaF(i,j,k)) !/ (1.0_dp - part%df(m))
!part%X_K40(m) = part%X_K40(m) * (1.0_dp - part%df(m))**(1.0_dp/pX%dK40) !/ (1.0_dp-mc%DeltaF(i,j,k)) !/ (1.0_dp - part%df(m))
part%X_U238(m) = max(0.0_dp,part%X_U238(m) *(1.0_dp-mc%DeltaF(i,j,k))**(1.0_dp/pX%dU238) / (1.0_dp-mc%DeltaF(i,j,k))) ! /(1-dF) to fill up whole cell again with same concentration
part%X_U235(m) = max(0.0_dp,part%X_U235(m) *(1.0_dp-mc%DeltaF(i,j,k))**(1.0_dp/pX%dU235) / (1.0_dp-mc%DeltaF(i,j,k)))
part%X_Th232(m)= max(0.0_dp,part%X_Th232(m)*(1.0_dp-mc%DeltaF(i,j,k))**(1.0_dp/pX%dTh232) / (1.0_dp-mc%DeltaF(i,j,k)))
......@@ -1392,7 +1407,7 @@ contains
! test with D=1 shows: V1 leads to too strong reduction in sum(H), V2 leads to equal error but increased sum(H) -> truth somewhere in the middle???
endif
!endif
endif
! 0 - mantle
! 1 - weak zone
......@@ -1407,12 +1422,13 @@ contains
! 10 - CC
if (part%r(m).gt.6.0_dp) then !if (part%r(m).ge.7.0_dp) then
compr_fac = cell%ref_r(i,j,k)
! all crust particles get heat source enrichment, not just the new produced ones (where part%r(m) would be below zero)
! Version 2 -> add to crust (different values old and new crust)
if (cell%fld_x) part%X_U238(m) = part%X_U238(m) + mc%X_U238_add(i,j) / cell%ref_r(i,j,k) !cell%X_U238(i,j,k)
if (cell%fld_x) part%X_U235(m) = part%X_U235(m) + mc%X_U235_add(i,j) / cell%ref_r(i,j,k) !cell%X_U235(i,j,k)
if (cell%fld_x) part%X_Th232(m) = part%X_Th232(m) + mc%X_Th232_add(i,j) / cell%ref_r(i,j,k) !cell%X_Th232(i,j,k)
if (cell%fld_x) part%X_K40(m) = part%X_K40(m) + mc%X_K40_add(i,j) / cell%ref_r(i,j,k) !cell%X_K40(i,j,k)
if (cell%fld_x) part%X_U238(m) = part%X_U238(m) + mc%X_U238_add(i,j) !/ compr_fac !cell%X_U238(i,j,k) ! do not divide again by density since was already divided by melt mass above
if (cell%fld_x) part%X_U235(m) = part%X_U235(m) + mc%X_U235_add(i,j) !/ compr_fac !cell%X_U235(i,j,k)
if (cell%fld_x) part%X_Th232(m) = part%X_Th232(m) + mc%X_Th232_add(i,j) !/ compr_fac !cell%X_Th232(i,j,k)