Update heat sources redistribution code 1D and 2D

Octave/call_gas_speciation.m

-buffer = -2:0.1:4;
+%buffer = -2:0.1:4;
+%
+%for i=1:length(buffer)
+%[H2(i), H2O(i), CO(i), CO2(i)]=ortenzi_simultaneous_eq_C_O_H(1000, 1, 0.5, 0.5, 2, buffer(i));
+%end
 
-for i=1:length(buffer)
-[H2(i), H2O(i), CO(i), CO2(i)]=ortenzi_simultaneous_eq_C_O_H(1800, 1, 0.5, 0.5, 2, buffer(i));
+T = 1000:100:2000;
+
+for i=1:length(T)
+%for i=1:length(buffer)
+[H2(i), H2O(i), CO(i), CO2(i)]=ortenzi_simultaneous_eq_C_O_H(T(i), 1, 0.5, 0.5, 2, 0);
 end
 
-  f = figure('visible','off');
+
+%  f = figure('visible','off');
+figure
   
-plot(buffer,H2*100,'r','linewidth',2)
+plot(T,H2*100,'r','linewidth',2)
+%plot(buffer,H2*100,'r','linewidth',2)
 hold on
-plot(buffer,H2O*100,'b','linewidth',2)
-plot(buffer,CO*100,'r--','linewidth',2)
-plot(buffer,CO2*100,'b--','linewidth',2)
+plot(T,H2O*100,'b','linewidth',2)
+plot(T,CO*100,'r--','linewidth',2)
+plot(T,CO2*100,'b--','linewidth',2)
 hold off
-xlabel('From reduced to oxidized','fontsize',16)
+xlabel('Temperature in K','fontsize',16)
+%xlabel('From reduced to oxidized','fontsize',16)
 %xlabel('Redox variation from IW in log_{10} steps','fontsize',16)
 ylabel('Abundance in %','fontsize',16)
 %ylabel('Molar gas fraction','fontsize',16)
@@ -23,5 +34,5 @@ text(0,3,'50% H_2O, 50% CO_2')
 %text(0,0.07,'1800 K, 1 bar')
 %text(0,0.04,'50% H_2O, 50% CO_2')
 
- filename='/home/noacklen/ownCloud/Conferences/International/2019/09_HearingFWF_Vienna/speciation.png';
- print(filename,'-dpng',f);
\ No newline at end of file
+% filename='/home/noacklen/ownCloud/Conferences/International/2019/09_HearingFWF_Vienna/speciation.png';
+% print(filename,'-dpng',f);
\ No newline at end of file

Octave/read_therm_evol_dim0.m

@@ -44,11 +44,11 @@ function read_therm_evol_dim0()
   plot(t_Gyr,data(:,7));
   plot(t_Gyr,data(:,6));
   plot(t_Gyr,data(:,5));
-  plot(t_Gyr,data(:,3));
+  plot(t_Gyr,data(:,4));
   hold off
   xlabel("Time (Gyr)");
   ylabel("Temperature (K)");
-  legend('CMB','Tb','Tm','Tl','Tr')
+  legend('CMB','Tb','Tm','Tl','Tcr')
 
   subplot(2,2,2)
   plot(t_Gyr,q_s);

Octave/solubility.m

@@ -35,15 +35,26 @@ M_carb = M_C + 3*M_O; %molar mass of carbonate CO_2^3-
 %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;
-
+%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;
+%X_Fe2O3 = 0;
+
+% 1921 Kilauea tholeiitic basalt, Holloway 1992, in mass fraction:
+X_K2O = 0.0051;
+X_Na2O = 0.0215; 
+X_CaO = 0.1093;
+X_MgO = 0.1041;
+X_FeO = 0.0968;
+X_Fe2O3 = 0.0172;
+X_Al2O3 = 0.0172;
+X_SiO2 = 0.4916;
+X_TiO2 = 0.1333;
 
 % basaltic example composition, all in weight fractions
 %X_K2O = 0.0772;
@@ -54,6 +65,7 @@ X_TiO2 = 0.0084;
 %X_Al2O3 = 0.1557;
 %X_SiO2 = 0.4989;
 %X_TiO2 = 0.0089;
+%X_Fe2O3 = 0;
 
 %% Dixon 1997, MORB Thol. Basalt, wt-fractions:
 %X_K2O = 0.07/100;
@@ -64,59 +76,70 @@ X_TiO2 = 0.0084;
 %X_Al2O3 = 16.4/100;
 %X_SiO2 = 49.1/100;
 %X_TiO2 = 0.74/100;
-%%X_K2O = 1 - X_Na2O - X_CaO - X_MgO - X_FeO - X_Al2O3 - X_SiO2 - X_TiO2;
+%X_Fe2O3 = 0;
+%%X_K2O = 1 - X_Na2O - X_CaO - X_MgO - X_FeO - X_Fe2O3 - X_Al2O3 - X_SiO2 - X_TiO2;
 
 % 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_Fe2O3 + 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);
+%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_Fe2O3/(2*M_Fe+3*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);
+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_Fe2O3 = X_Fe2O3/Sum_mol/(2*M_Fe+3*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;
-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_mol = X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Fe2O3 + 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;
+%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_Fe2O3 = X_Fe2O3/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_Fe2O3 + 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) ...
-%       + 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))
+M_melt = (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_Fe2O3*(2*M_Fe+3*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_MgO*(M_Mg+M_O) + X_FeO*(M_Fe+M_O) + X_Fe2O3*(2*M_Fe+3*M_O)/3 ...
 %%       + X_Al2O3*(2*M_Al+3*M_O)/3 + X_SiO2*(M_Si+2*M_O)/2 + X_TiO2*(M_Ti+2*M_O)/2)
 %%M_melt = M_melt / M_O
 
+%fwm = (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_Fe2O3*(2*M_Fe+3*M_O)/3 ...
+%       + X_Al2O3*(2*M_Al+3*M_O)/3 + X_SiO2*(M_Si+2*M_O)/2 + X_TiO2*(M_Ti+2*M_O)/2)
+fwm = (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_Fe2O3*(2*M_Fe+3*M_O)/3 ...
+       + X_Al2O3*(2*M_Al+3*M_O)/3 + X_SiO2*(M_Si+2*M_O)/2 + X_TiO2*(M_Ti+2*M_O)/2)
+
 X_AI = X_Al2O3 / (X_CaO+X_K2O+X_Na2O);
 
 
 if hydrous==0
   % non-bridging oxygen divided by ixygen
   NBO = 2 * (X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO - X_Al2O3);
-  NBO_O = NBO / ( 2 * X_SiO2 + 2 * X_TiO2 + 3 * X_Al2O3 + X_MgO + X_FeO + X_CaO + X_Na2O + X_K2O );
+  NBO_O = NBO / ( 2 * X_SiO2 + 2 * X_TiO2 + 3 * X_Al2O3 + X_MgO + X_FeO + 3*X_Fe2O3 + X_CaO + X_Na2O + X_K2O );
 else
   NBO = 2 * (X_H2O + X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO - X_Al2O3);
-  NBO_O = NBO / ( 2 * X_SiO2 + 2 * X_TiO2 + 3 * X_Al2O3 + X_MgO + X_FeO + X_CaO + X_Na2O + X_K2O + X_H2O );
+  NBO_O = NBO / ( 2 * X_SiO2 + 2 * X_TiO2 + 3 * X_Al2O3 + X_MgO + X_FeO + 3*X_Fe2O3 + X_CaO + X_Na2O + X_K2O + X_H2O );
 end
 
 if hydrous==1
@@ -156,6 +179,7 @@ P_H2O = H2O_fl_mole*P_bar;
 
 ln_H2O = a_H2O*log(P_H2O) + b_H2O*NBO_O + B_H2O + C_H2O*P_bar/T_K;     % in wt-%
 
+% in formula here Fe2O3 missing since not in original paper
 if hydrous==1
 ln_CO_3_2min = X_H2O*d_H2O + 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;

src/input.f90

@@ -1009,7 +1009,10 @@ contains
           if (rank.eq.0) write(*,*) var,"=",pM%update_solidus
         case("melt_depth") ! in GPa, typically at 8 to 12 or so...
           if (pE%DimValues.eq.1) then
-            pM%melt_depth=(val*1.0e9_dp-pE%Psurf)/(pI%rho*pI%gravity)  ! (melt_pressure - Psurf) / rho * g 
+            ! 1D code, there used as actual pressure, not depth! Also surface pressure taken care of in therm_evol_help_f90!
+            !pM%melt_depth=(val*1.0e9_dp-pE%Psurf)/(pI%rho*pI%gravity)  ! (melt_pressure - Psurf) / rho * g 
+            ! TODO: ADAPT IN PATRICK'S CODE, TOO
+            pM%melt_depth=val*1.0e9_dp
             if (pM%melt_depth.lt.0.0_dp) pM%melt_depth = 0.0_dp
           else if (pF%D.ne.1.0_dp) then
             pM%melt_depth=val*1.0e9_dp/(pF%D*pF%rho*pF%g) - pE%Psurf

src/main.f90

@@ -591,14 +591,14 @@ program CHIC
 !    if (pI%C_U.ne.0.0_dp) call initialize_heat(pC,pI,pF,pG,pX,pE%H_from_surf,pM%Lambda,cell) ! needs to be before crust_init
     if (pH%use_magn_H.ne.0.0_dp) call initialize_magn_heat(pF,pH,mesh%rc,mesh%nr,pI%prof_nm,0.0_dp,cell) 
 
-    if (pM%use_melt.eq.1) then
+!    if (pM%use_melt.eq.1) then ! always load melt curves since cut_T_solidus otherwise doesn't work
       call init_melt_curves(mesh,mc,pF,pM,pG%dim,cell,pI%X_FeM)
       if (pG%dim.eq.2) then
         mc%DeltaF(0:mesh%nl+1,1,0:mesh%nr+1) = 0.0_dp
       else if (pG%dim.eq.3) then
         mc%DeltaF(0:mesh%nl+1,0:mesh%ny+1,0:mesh%nr+1) = 0.0_dp
       endif
-    endif
+!    endif
 
     if ((pT%linear.ge.7).or.(pT%CutSolidus)) call cut_T_solidus(field,mc)
 
@@ -671,7 +671,7 @@ program CHIC
       endif
     endif
 
-    if ((pP%use_part.eq.1).and.(pW%weak_angle.eq.0.0_dp).and.(pW%BenchSubdFreeS.eq.0)) &
+    if ((pP%use_part.eq.1).and.(pW%weak_angle.eq.0.0_dp).and.(pW%BenchSubdFreeS.eq.0).and.(pD%nr_ph.ne.0)) &
       &call get_phases(part,cell,jmin_2D,field%T,field%w,pD,pF,pI,mesh%rmax,pE%Psurf,pI%compress,pN%debug,t)
     if (pD%nr_ph.gt.0.and.pD%set_tracers.eq.2) then
       ! trace material exchange -> set to initial material value, will not be updated later on
@@ -850,7 +850,7 @@ program CHIC
     endif  
    endif
   else ! 1D output
-    open(unit=41,file="data_not_used.res",form='formatted',status='unknown')
+    open(unit=41,file="data_heat_sources.res",form='formatted',status='unknown')
     open(unit=51,file="data_char_val_ts.res",form='formatted',status='unknown')
     open(unit=53,file="data_char_val_ts2.res",form='formatted',status='unknown')
     open(unit=43,file="data_ini_profs.res",form='formatted',status='unknown')
@@ -964,6 +964,13 @@ program CHIC
       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
       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)
+          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?
       endif
 !      str_out1 = "D(VD,W)="
@@ -1271,6 +1278,12 @@ program CHIC
       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
       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)
+          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?
       endif
 
@@ -1342,6 +1355,7 @@ program CHIC
     field%T_vOld = field%T_old ! T^{n-2}
     field%T_old = field%T ! temperature of last time step, T^{n-1}
 !    field%T_iter = field%T_old ! temperature of last iteration, for the first iteration it is the temperature of the last time step
+
     field%p_old = field%p
     cell%ref_p_old = cell%ref_p
     if (cell%fld_c) cell%C_vOld = cell%C_old
@@ -1686,6 +1700,7 @@ program CHIC
         if(rank.eq.0) call print_time(pN%pr_time,timeD1,timeD2,"Post-Exch v,p")
       enddo ! inner iteration, stokes equation
 
+
       !if ((pN%debug.ge.2).and.(rank.eq.0)) write(*,*) "--------------------------------------------------------"
 
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -1727,10 +1742,10 @@ program CHIC
       if (pS%minStrRate.gt.0.0_dp) field%strain_rate = max(field%strain_rate,pS%minStrRate)
       call Visc_Diss(field,mesh,jmin_2D,eta%CV,pG%geom,pS%str_pre_fac,pI%compress,pN%NoVD)
 
-
       ! Determine T and C for next time step
       if ((pN%debug.ge.2).and.(rank.eq.0)) write(*,*) "Energy and compositional solver"
 
+
       finished_II=0
       call exchange_int(finished_II,0)
       i_inner=0
@@ -1905,6 +1920,12 @@ program CHIC
              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
              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)
+                 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
              endif
 
@@ -1977,7 +1998,7 @@ program CHIC
       call part_cell4part(part,mesh,cell,mc,pG%dim,pN%debug) ! get new cells 
       if(rank.eq.0) call print_time(pN%pr_time,timeD1,timeD2,"Part/Post-C4P")
       !if ((abs(pI%compress).ge.3).or.(pD%nr_ph.ne.0)) &
-      if ((pW%weak_angle.eq.0.0_dp).and.(pW%BenchSubdFreeS.eq.0)) &
+      if ((pW%weak_angle.eq.0.0_dp).and.(pW%BenchSubdFreeS.eq.0).and.(pD%nr_ph.ne.0)) &
       & call get_phases(part,cell,jmin_2D,field%T,field%w,pD,pF,pI,mesh%rmax,pE%Psurf,pI%compress,pN%debug,t) !SetMatProp(part,jmin_2D,cell%r,mesh%rc,field%T,pE%Psurf,pD,pF,pI%compress)
       call part_part2cell(part,mc,mesh,cell,pG%geom,pN%debug,pG%periodic,pN%average,pN%p_average,pN%p_average2,&
                          &pN%pr_time,pP%moveRef) ! transfer particle properties to cells 
@@ -2014,6 +2035,7 @@ program CHIC
       if (pM%add_crust.eq.1) then
         call add_crust(mesh,cell,mc,part,t,dt,pC,pX)
         if (pM%add_crust_depth.gt.0.0_dp) call distr_crust(mesh,mc,pM%add_crust_depth,pM%add_crust_incr)
+
         !call move_crust(mc,mesh,field,dt)
         call set_new_crust(mc,mesh,cell,pC,pI,pW,pX,pE%Dreg,pI%k_r,part,t)
       endif
@@ -2027,6 +2049,9 @@ program CHIC
       if (pS%strain_max.ne.0.0_dp) call strain_weakening(field,mesh,part,pC,pI,pS,pW,pG%dim,dt) ! change the friction angle part%fa with strain
 
       if (pM%Rt.ne.0.0_dp) call melt_depletion(part,field,cell,mesh,mc,pM,pI,dt) ! depletion via specific function, was not determined in melt_comp_p
+
+      call part_part2cell(part,mc,mesh,cell,pG%geom,pN%debug,pG%periodic,pN%average,pN%p_average,pN%p_average2,&
+                           &pN%pr_time,pP%moveRef)
     endif
 
 !    if (abs(pI%compress).ge.3) then 
@@ -2111,6 +2136,12 @@ program CHIC
       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
       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)
+          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
       endif
 

src/particles.f90

@@ -1792,6 +1792,7 @@ contains
       if (cell%fld_x) cell%X_K40(0,:,:) = cell%X_K40(nl,:,:)
       if (cell%fld_x) cell%X_K40(nl+1,:,:) = cell%X_K40(1,:,:)
 
+      ! Todo: below ny instead of nl!!!
       if (ny.gt.1) then
         if (cell%fld_c) cell%c(:,0,:) = cell%c(:,nl,:)
         if (cell%fld_c) cell%c(:,nl+1,:) = cell%c(:,1,:)
@@ -1922,6 +1923,7 @@ contains
       if (cell%fld_x) cell%X_K40(0,:,0:) = cell%X_K40(1,:,0:)
       if (cell%fld_x) cell%X_K40(nl+1,:,0:) = cell%X_K40(nl,:,0:)
 
+      ! todo below: ny instead of nl
       if (ny.gt.1) then
         if (cell%fld_c) cell%c(:,0,:) = cell%c(:,1,:)
         if (cell%fld_c) cell%c(:,nl+1,:) = cell%c(:,nl,:)
@@ -2017,7 +2019,7 @@ contains
     integer,intent(in)::cell_l,cell_y,cell_r,geom,di,average
     real(dp),intent(inout)::u_tot,v_tot,w_tot
 
-    integer::i,j,k,i_,j_,k_,jmin
+    integer::i,j,k,i_,j_,k_,k__,jmin
     real(dp)::Ws_tot,weight
 
     if (di.eq.2) then
@@ -2158,6 +2160,8 @@ contains
           endif
           k_=k+cell_r
           if (k_.gt.mesh%nr+1) k_=mesh%nr+1
+          k__ = k_
+          if (k__.lt.1) k__=1 ! for vel w vector
           if (geom.eq.0) then ! diagonal of a cell divided by distance
             if (di.eq.2) then
               weight = sqrt(mesh%dl_vec(i_)**2+mesh%dr_vec(k_)**2) / sqrt((pos_r-mesh%rb(k_))**2 + (pos_l-mesh%lc(i_))**2)
@@ -2173,11 +2177,11 @@ contains
 !            weight = sqrt((mesh%dl_vec(i_)*mesh%rc(k_))**2+mesh%dr_vec(k_)**2) / sqrt((pos_r-mesh%rc(k_))**2 + ((pos_l-mesh%lc(i_))*mesh%rc(k_))**2)
           endif
           if (average.eq.1) then
-            w_tot = w_tot+weight/field%w(i_,j_,k_)
+            w_tot = w_tot+weight/field%w(i_,j_,k__)
           elseif (average.eq.2) then
-            w_tot = w_tot+field%w(i_,j_,k_)*weight
+            w_tot = w_tot+field%w(i_,j_,k__)*weight
           elseif (average.eq.3) then
-            w_tot = w_tot*field%w(i_,j_,k_)**weight
+            w_tot = w_tot*field%w(i_,j_,k__)**weight
           else
             if (di.eq.2) then
               j_=1

src/therm_evol.f90

@@ -515,7 +515,6 @@ contains
     call grav_New(g,pI%Rp,pI%Rc,Rcr,Rr,pI%Rp,pI%rho_c,pI%rho,pI%rho_cr,pI%rho_r,pE%Gr,dum1,dum2,dum3)
     write(*,*) "S: ",pI%Rp,g,g*pF%g
 
-
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     ! core freezing initial values !
     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -587,6 +586,21 @@ contains
                 & pI%rho*pI%gravity,debug,pE%Pc*pF%rho*pF%g*pF%D*1.0e-9_dp, pF%rho*pF%g*pF%D/(pF%R*pF%DeltaT) )
     write(*,*) "eta_um=",eta,", eta_c=",eta_c,", D_eta=",D_eta
 
+    !!!!!!!!!!!!!!!!!!!!!
+    ! mantle properties !
+    !!!!!!!!!!!!!!!!!!!!!
+
+    ! TODO: ADAPT IN PATRICK'S CODE, TOO
+    tet%rho_b = pI%rho
+    tet%alpha_b = pI%alpha
+    tet%Cp_b = pI%Cp
+    tet%k_b = pI%k
+    tet%rho_um = pI%rho
+    tet%alpha_um = pI%alpha
+    tet%Cp_um = pI%Cp
+    tet%k_um = pI%k
+
+
     !!!!!!!!!!!!!!!!!
     !! ocean layer !!
     !!!!!!!!!!!!!!!!!
@@ -921,7 +935,7 @@ contains
     logical::finished,convection
     real(dp)::dummyT
     real(dp)::A_t,A_b,V_tb
-    integer::k_min,iter_sub
+    integer::k_min,iter_sub,layer
     real(dp)::dt_s
 
     if (debug.gt.0) write(*,*) "Param. model at t=",t
@@ -988,10 +1002,42 @@ contains
 
     do while (.not.(finished)) 
 
-      call getProp_ND(pE%Pc,Tc,rho_b,alpha_b,Cp_b,k_b,0,error,pF,7)
-      call getProp_ND(tet%pOMB,Tm,rho_um,alpha_um,Cp_um,k_um,0,error,pF,7)
-
-      tet%rho_um = rho_um
+      ! TODO: ADAPT IN PATRICK'S CODE, TOO
+      alpha_b = tet%alpha_b
+      rho_b = tet%rho_b
+      Cp_b = tet%Cp_b
+      k_b = tet%k_b
+      alpha_um = tet%alpha_um
+      rho_um = tet%rho_um
+      Cp_um = tet%Cp_um
+      k_um = tet%k_um
+
+      layer = 5
+      if (pE%Pc.gt.(22.9_dp-0.0025_dp*(Tc-2260.0_dp))*1.0e+9_dp) layer = 6
+      if (pE%Pc.gt.(109.54_dp+0.0058_dp*Tc)*1.0e+9_dp) layer = 7
+      call getProp_ND(pE%Pc,Tc,tet%rho_b,tet%alpha_b,tet%Cp_b,tet%k_b,0,error,pF,layer) 
+      layer = 5
+      if (tet%pOMB.gt.(22.9_dp-0.0025_dp*(Tm-2260.0_dp))*1.0e+9_dp) layer = 6
+      if (tet%pOMB.gt.(109.54_dp+0.0058_dp*Tm)*1.0e+9_dp) layer = 7
+      call getProp_ND(tet%pOMB,Tm,tet%rho_um,tet%alpha_um,tet%Cp_um,tet%k_um,0,error,pF,layer)
+ 
+      if (debug.gt.2) write(*,*) "EOS c:",pE%Pc,Tc,tet%rho_b,tet%alpha_b,tet%Cp_b,tet%k_b
+      if (debug.gt.2) write(*,*) "EOS m:",tet%pOMB,Tm,tet%rho_um,tet%alpha_um,tet%Cp_um,tet%k_um
+      
+      if (tet%rho_b.lt.0) then
+        if (debug.gt.1) write(*,*) "----- Error with EOS at CMB!!! Use values from input.txt instead -----"
+        tet%rho_b = rho_b
+        tet%alpha_b = alpha_b
+        tet%Cp_b = Cp_b
+        tet%k_b = k_b
+      endif
+      if (tet%rho_um.lt.0) then
+        if (debug.gt.1) write(*,*) "----- Error with EOS at upper mantle!!! Use values from input.txt instead -----"
+        tet%rho_um = rho_um
+        tet%alpha_um = alpha_um
+        tet%Cp_um = Cp_um
+        tet%k_um = k_um
+      endif
 
       Tb_last = Tb
       Dcr_Del = tet%Dcr_Del
@@ -1026,7 +1072,6 @@ contains
         DeltaTi = abs(Tm - Ts_surf) + abs(Tc - Tb) !Tm - Ts_surf + Tc - Tb !abs(Tm - Ts) + abs(Tc - Tb)
         if (debug.gt.1) write(*,*) "T: ", DeltaT,":", Tcr, Tl, Tm, Tb, Tc
  
-        if (debug.gt.1) write(*,*) "Ra: ",pI%Ra, pI%alpha, pI%rho, pI%gravity, DeltaT, (Rl-pI%Rc)**3, (pI%kappa*eta )
         !call getRa( Ra,Rai,Rai_crit,Rai_crit_c,pI,eta,D_eta,DeltaT,DeltaTi,Rl,pI%Rc,pE%Ra_crit,debug )
         Ra = get_Ra( pI%Ra, pI%alpha, pI%rho, pI%gravity, DeltaT, Rl, pI%Rc, pI%kappa, eta )
         Rai= get_Ra( pI%Ra, pI%alpha, pI%rho, pI%gravity, DeltaTi, Rl, pI%Rc, pI%kappa, eta ) ! use pI%Rp or Rl? With Rl (e.g. Rueckriemen), results look smoother and convection slightly longer
@@ -1039,6 +1084,7 @@ contains
         !  Rai_crit_c = 11.74_dp*(log(D_eta))**4
         !endif
         Rai_crit_c = Rai_crit + max(0.0_dp,11.74_dp*(log(D_eta))**4)
+        if (debug.gt.1) write(*,*) "Ra: ",pI%Ra, pI%alpha, pI%rho, pI%gravity,DeltaT,(Rl-pI%Rc)**3,(pI%kappa*eta),Ra,Rai,Rai_crit
 
         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         ! Get new TBL, check if now stagnant !
@@ -1065,15 +1111,15 @@ contains
           Delta_c_max = pE%Delta_c_max 
         endif
 
-        kappa_b = k_b / (Cp_b*rho_b)
+        kappa_b = tet%k_b / (tet%Cp_b*tet%rho_b)
 
         Delta_ct_old = Delta_ct
         Delta_ct = min((pI%kappa*eta_c*Rai_crit_c/( pI%Ra*pI%alpha*pI%rho*gc*DeltaTi))**pE%beta,Delta_c_max) ! pI%Ra used for re-dimensionalization of factors... 
         !Delta_ct = min((pI%kappa*eta_c*Rai_crit_c/( pI%Ra*pI%alpha*pI%rho*gc*abs(Tc-Tb)))**pE%beta,Delta_c_max) ! this is the version from Morschhauser et al., above leads to thinner Delta_ct
 
-        !Delta_ct = (kappa_b*eta_c*Rai_crit_c/( pI%Ra*alpha_b*rho_b*gc*abs(Tc-Tb)))**pE%beta ! for Mercury, this leads to way too high lower TBL, need to be cut, unphysical...
-        !write(*,*) eta_c*pF%eta_ref,Delta_ct*pF%D,pI%Ra*alpha_b*rho_b*gc*abs(Tc-Tb)/(kappa_b*eta_c),Rai_crit_c
-        !write(*,*) pI%Ra,alpha_b,rho_b,gc,abs(Tc-Tb),kappa_b
+        !Delta_ct = (kappa_b*eta_c*Rai_crit_c/( pI%Ra*tet%alpha_b*tet%rho_b*gc*abs(Tc-Tb)))**pE%beta ! for Mercury, this leads to way too high lower TBL, need to be cut, unphysical...
+        !write(*,*) eta_c*pF%eta_ref,Delta_ct*pF%D,pI%Ra*tet%alpha_b*tet%rho_b*gc*abs(Tc-Tb)/(kappa_b*eta_c),Rai_crit_c
+        !write(*,*) pI%Ra,tet%alpha_b,tet%rho_b,gc,abs(Tc-Tb),kappa_b
 
 !TS     51 t= 2.531E+07 dt= 1.000E+06 Tr= 290.00 Tcr= 563.95 Tl=1576.09 Tm=1963.36 Tb= 2046.04 Tc= 6337.92 Dr=0.000E+00 Dcr=5.000E+03 Dl=2.503E+04 Dut=7.798E+03 Dct=5.392E+06 Deff=5.610E+05
 !  8.728771783613812E+019   24661.2566165471       4.846486599710517E-012   72.6896796768864
@@ -1112,8 +1158,10 @@ contains
 
       !ql = pI%k * (Tm-Tl)/Delta_ut 
       !qc = pI%k * (Tc-Tb)/Delta_ct 
-      ql = k_um * (Tm-Tl)/Delta_ut
-      qc = k_b  * (Tc-Tb)/Delta_ct
+      ql = tet%k_um * (Tm-Tl)/Delta_ut
+!      write(*,*) "ql:",ql,k_um,Tm,Tl,Delta_ut
+
+      qc = tet%k_b  * (Tc-Tb)/Delta_ct
 
       !!!!!!!!!!!!!
       ! Melt zone !
@@ -1122,7 +1170,7 @@ contains
       ! ToDo: include melt in iteration, but which temperature to use? from last time step? But new Dl?
       if ((convection).and.(pM%use_melt.ne.0)) then ! if still convecting
         call getMelt(Va,ma,dma_dTm,XavU238i,XavU235i,XavTh232i,XavK40i,mz,mT,mTs,mTl,tet,pE,pM,pI,Rl,Dl,Delta_ut,Delta_ct,&
-                &Tc,Tl,debug,Tprof,r,tet%pOMB,pE%DimValues,rho_um)
+                &Tc,Tl,debug,Tprof,r,tet%pOMB,pE%DimValues,tet%rho_um)
       else
         allocate( mz(1),mT(1),mTs(1),mTl(1) )
         Va = 0.0_dp
@@ -1142,12 +1190,12 @@ contains
 
       if (convection) then ! if still convecting
         if (Dcr.eq.0.0_dp) then
-          derivTrRl = -Q_m*Rl/(3.0_dp*k_um) - (Tl-Ts_surf-Q_m*(pI%Rp**2-Rl**2)/(6.0_dp*k_um))&
+          derivTrRl = -Q_m*Rl/(3.0_dp*tet%k_um) - (Tl-Ts_surf-Q_m*(pI%Rp**2-Rl**2)/(6.0_dp*tet%k_um))&
                                                 &/(Rl**2.0_dp*(1.0_dp/Rl-1.0_dp/pI%Rp))
 !          derivTrRl = -Q_m*Rl/(3.0_dp*pI%k_l) - (Tl-Ts_surf-Q_m*(pI%Rp**2-Rl**2)/(6.0_dp*pI%k_l))&
 !                                                &/(Rl**2.0_dp*(1.0_dp/Rl-1.0_dp/pI%Rp))
         else
-          derivTrRl = -Q_m*Rl/(3.0_dp*k_um) - (Tl-Tcr-Q_m*(Rcr**2-Rl**2)/(6.0_dp*k_um))&
+          derivTrRl = -Q_m*Rl/(3.0_dp*tet%k_um) - (Tl-Tcr-Q_m*(Rcr**2-Rl**2)/(6.0_dp*tet%k_um))&
                                                 &/(Rl**2.0_dp*(1.0_dp/Rl-1.0_dp/Rcr))
 !          derivTrRl = -Q_m*Rl/(3.0_dp*pI%k_l) - (Tl-Tcr-Q_m*(Rcr**2-Rl**2)/(6.0_dp*pI%k_l))&
 !                                                &/(Rl**2.0_dp*(1.0_dp/Rl-1.0_dp/Rcr))
@@ -1183,8 +1231,8 @@ contains
         if (debug.gt.2) write(*,*) ql,pI%rho_cr,pM%LatH,pI%Cp_cr,Tm-Tl,dDcrdt,pI%k,derivTrRl
 
         if (pE%PT.eq.0) then