Commit 02b0c809 authored by Lena Noack's avatar Lena Noack
Browse files

Several updates evaluation software

parent bc41363f
......@@ -52,6 +52,7 @@ function input = read_input(fpath)
input.DeltaVisc = value(strcmp("DeltaVisc", data));
input.use_melt = value(strcmp("use_melt", data));
input.Fe_Mantle = value(strcmp("Fe_Mantle", data));
if (isempty(input.Fe_Mantle)), input.Fe_Mantle=10;, endif
if (isempty(input.use_melt)), input.use_melt=0;, endif
if (input.use_melt==1)
input.m1s0 = value(strcmp("melt1s_0", data)); % K
......
......@@ -17,6 +17,7 @@ M_Ti = 47.867;
M_Fe = 55.845;
M_CO2 = M_C + 2*M_O;
M_H2O = 2*M_H + M_O;
M_P = 30.973762;
M_carb = M_C + 3*M_O; %molar mass of carbonate CO_2^3-
% partial pressures P_CO2 from previous degassing, how to integrate that?
......@@ -34,6 +35,19 @@ M_carb = M_C + 3*M_O; %molar mass of carbonate CO_2^3-
% hydrous = 1;
%endif
% Etna composition, all in weight fractions
%X_K2O = 0.0199;
%X_Na2O = 0.0345;
%X_CaO = 0.1093;
%X_MgO = 0.0576;
%X_FeO = 0.1024;
%X_Al2O3 = 0.1732;
%X_SiO2 = 0.4795;
%X_TiO2 = 0.0167;
%X_P2O5 = 0.0051;
%X_Fe2O3 = 0;
%X_H2O = X_H2O_wt; %molar fraction of water in the melt
% basaltic example composition, aparently in molar fractions, Saal et al 2012, D-6-1 sample
X_K2O = 0.00023;
X_Na2O = 0.0194;
......@@ -44,6 +58,7 @@ X_Al2O3 = 0.164;
X_SiO2 = 0.4957;
X_TiO2 = 0.0084;
X_Fe2O3 = 0;
X_P2O5 = 0;
% 1921 Kilauea tholeiitic basalt, Holloway 1992, in mass fraction:
%X_K2O = 0.0051;
......@@ -55,6 +70,7 @@ X_Fe2O3 = 0;
%X_Al2O3 = 0.0172;
%X_SiO2 = 0.4916;
%X_TiO2 = 0.1333;
%X_P2O5 = 0;
% basaltic example composition, all in weight fractions
%X_K2O = 0.0772;
......@@ -66,6 +82,7 @@ X_Fe2O3 = 0;
%X_SiO2 = 0.4989;
%X_TiO2 = 0.0089;
%X_Fe2O3 = 0;
%X_P2O5 = 0;
%% Dixon 1997, MORB Thol. Basalt, wt-fractions:
%X_K2O = 0.07/100;
......@@ -77,6 +94,7 @@ X_Fe2O3 = 0;
%X_SiO2 = 49.1/100;
%X_TiO2 = 0.74/100;
%X_Fe2O3 = 0;
%X_P2O5 = 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:
......@@ -85,9 +103,9 @@ X_Fe2O3 = 0;
%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_Fe2O3/(2*M_Fe+3*M_O) + X_P2O5/(2*M_P+5*M_O) + X_H2O/(2*M_H+M_O);
%
%X_H2O = X_H2O_wt/Sum_mol/(2*M_H+M_O);
%X_H2O = X_H2O/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);
......@@ -97,9 +115,10 @@ X_Fe2O3 = 0;
%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_P2O5 = X_P2O5/Sum_mol/(2*M_P+5*M_O);
% if already in molar fractions:
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);
Sum_mol = X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Fe2O3 + X_Al2O3 + X_SiO2 + X_TiO2 + X_P2O5 + 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;
......@@ -111,24 +130,25 @@ X_Fe2O3 = X_Fe2O3/Sum_mol;
X_Al2O3 = X_Al2O3/Sum_mol;
X_SiO2 = X_SiO2/Sum_mol;
X_TiO2 = X_TiO2/Sum_mol;
X_P2O5 = X_P2O5/Sum_mol;
% sum of all should be 1:
%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_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));
%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)
%%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_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);
......@@ -178,15 +198,17 @@ P_CO2 = (1.0-H2O_fl_mole)*P_bar;
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-%
X_H2O_mol = exp(ln_H2O)/100/Sum_mol/(2*M_H+M_O); % in wt-amount
% 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;
else
ln_CO_3_2min = X_AI*d_AI + (X_FeO+X_MgO)*d_FeMg + (X_Na2O+X_K2O)*d_NaK ... % in ppm
%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
ln_CO_3_2min = X_H2O_mol*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;
end
%else
%ln_CO_3_2min = 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;
%end
H2O_dissolved = exp(ln_H2O); % in wt-%
if (H2O_fl_mole==0), H2O_dissolved = 0;, end
......
......@@ -15,6 +15,7 @@ M_Ti = 47.867;
M_Fe = 55.845;
M_CO2 = M_C + 2*M_O;
M_H2O = 2*M_H + M_O;
M_P = 30.973762;
M_carb = M_C + 3*M_O; %molar mass of carbonate CO_2^3-
% partial pressures P_CO2 from previous degassing, how to integrate that?
......@@ -22,7 +23,7 @@ P_bar = 4000; %p_lid;
T_K = 1200+273.15; %T_lid;
H2O_fl_mole = 0.01; % pre-defined mole fraction of H2O (over H2O+CO2) in gas bubble, calc partial pressures P_CO2/P_H2O from that
hydrous = 0; % unhydrous okay for water contents below ~4 wt-%; in online calculator only anhydrous version used
hydrous = 1; % unhydrous okay for water contents below ~4 wt-%; in online calculator only anhydrous version used
%if (H2O_fl_mole<0.5) % when to be anhydrous? for now if more CO2 moles than H2O moles
% hydrous = 0;
......@@ -39,7 +40,8 @@ X_FeO = 0.1024;
X_Al2O3 = 0.1732;
X_SiO2 = 0.4795;
X_TiO2 = 0.0167;
X_H2O = 0.01; %molar fraction of water in the melt
X_P2O5 = 0.0051;
X_H2O = 0.0001;%0.04;%0.0001; %weight fraction of water in the melt
%% basaltic example composition, all in weight fractions
%X_K2O = 0.0772;
......@@ -68,7 +70,7 @@ X_H2O = 0.01; %molar fraction of water in the melt
%%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/(2*M_H+M_O);
+ X_P2O5/(2*M_P+5*M_O) + X_H2O/(2*M_H+M_O);
X_H2O = X_H2O/Sum_mol/(2*M_H+M_O);
X_K2O = X_K2O/Sum_mol/(2*M_K+M_O);
......@@ -79,6 +81,7 @@ 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_P2O5 = X_P2O5/Sum_mol/(2*M_P+5*M_O);
% sum of all should be 1:
%X_K2O + X_Na2O + X_CaO + X_MgO + X_FeO + X_Al2O3 + X_SiO2 + X_TiO2
......@@ -94,7 +97,7 @@ X_AI = X_Al2O3 / (X_CaO+X_K2O+X_Na2O);
if hydrous==0
% non-bridging oxygen divided by ixygen
% non-bridging oxygen divided by oxygen
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 );
else
......@@ -145,13 +148,17 @@ for H2O_fl_mole=0:0.01:1
ln_H2O = a_H2O*log(P_H2O) + b_H2O*NBO_O + B_H2O + C_H2O*P_bar/T_K; % in wt-%
X_H2O_mol = exp(ln_H2O)/100/Sum_mol/(2*M_H+M_O); % in mol-amount
X_H2O_mol = min(X_H2O_mol,X_H2O);
% 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;
% 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
ln_CO_3_2min = X_H2O_mol*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;
% else
% ln_CO_3_2min = 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;
% endif
% end
H2O_dissolved = exp(ln_H2O); % in wt-%
if (H2O_fl_mole==0), H2O_dissolved = 0;, end
......@@ -160,7 +167,7 @@ for H2O_fl_mole=0:0.01:1
if (H2O_fl_mole==1), CO2_dissolved = 0;, end
% fprintf('%f %f %f %f\n',P_bar,H2O_fl_mole,H2O_dissolved,CO2_dissolved)
CO2(i) = CO2_dissolved;
CO2(i) = CO32m_dissolved;%CO2_dissolved;
H2O(i) = H2O_dissolved;
i=i+1;
end
......
......@@ -11,8 +11,8 @@ p_lid = 100; % in bar
T_lid = 1300 + 273.15; % in K
% solubility calculated from Iacono-Marzian et al., 2012; H2O in wt-% and CO2 in wt-ppm
X_H2O = 0.01;%0.1;%0.4;%0.1; %weight percent of water in the melt
X_CO2 = 10000;%600;%4000; %600; %weight ppm of CO2 in the melt
X_H2O = 0.1;%0.1;%0.4;%0.1; %weight percent of water in the melt
X_CO2 = 600;%600;%4000; %600; %weight ppm of CO2 in the melt
if (X_H2O<4) % personal communication with Fabrice Gaillard
hydrous=0;
......@@ -29,7 +29,7 @@ i_max = 200;
%fprintf('Start with CO2_melt = %f wt-ppm and H2O_melt = %f wt-%%; p=%f, T=%f\n',...
% X_CO2,X_H2O,p_lid,T_lid)
ln_p = [-3:0.01:3];
ln_p = [-6:0.01:3];
CO2_diss = zeros(size(ln_p));
H2O_diss = zeros(size(ln_p));
......@@ -41,7 +41,9 @@ H2O_fl_mole = 0.5; % = P_H2O/P_lid
no_water = 0; % if water pressure goes to zero, set no_water to one, otherwise wrong oscillating solution
i=1;
fac=1;
do
H2O_dissolved_old=1.0;
CO2_dissolved_old=1.0;
while ((abs(H2O_dissolved-H2O_dissolved_old)>tol && abs(CO2_dissolved-CO2_dissolved_old)>tol) && (i<i_max))
H2O_dissolved_old = H2O_dissolved;
CO2_dissolved_old = CO2_dissolved;
[H2O_dissolved,CO2_dissolved] = solubility(p_lid,T_lid,H2O_fl_mole,X_H2O,hydrous);
......@@ -65,7 +67,7 @@ do
dX_H2O = max(0,X_H2O-H2O_dissolved);
dX_CO2 = max(0,X_CO2-CO2_dissolved);
i=i+1;
until ((abs(H2O_dissolved-H2O_dissolved_old)<tol && abs(CO2_dissolved-CO2_dissolved_old)<tol) || (i>i_max))
end
%if (i>i_max)
% fprintf('Warning, iteration stopped after %d iterations at p=%f bar!\n',i,p_lid)
......@@ -76,17 +78,17 @@ H2O_diss(p_i) = 1-dX_H2O/X_H2O;
end
graphics_toolkit gnuplot
%graphics_toolkit gnuplot
%f = figure('visible','off');
semilogx(10.^(ln_p),CO2_diss,'r-',10.^(ln_p),H2O_diss,'b-')
set (findall (gcf (), '-property', 'fontsize'), 'fontsize', 20);
set (findall (gcf (), '-property', 'linewidth'), 'linewidth', 5);
title(sprintf(strcat('Melt (ppm): X_{H2O}=',num2str(X_H2O*10000,"%4.0f"),'; X_{CO2}=',num2str(X_CO2,"%4.0f"))),'FontSize',30)
set (findall (gcf (), '-property', 'fontsize'), 'fontsize', 12);
set (findall (gcf (), '-property', 'linewidth'), 'linewidth', 2);
title(sprintf(strcat('Melt (ppm): X_{H2O}=',num2str(X_H2O*10000,"%4.0f"),'; X_{CO2}=',num2str(X_CO2,"%4.0f"))),'FontSize',14)
%title(sprintf(strcat('Melt: X_{H2O}=',num2str(X_H2O,"%4.3f"),'; X_{CO2}=',num2str(X_CO2/10000.0,"%4.3f"),' wt-%%')),'FontSize',24)
xlabel('Total pressure in bar','FontSize',28)
ylabel('Abundance of volatiles dissolved in melt','FontSize',28)
xlabel('Total pressure in bar','FontSize',14)
ylabel('Abundance of volatiles dissolved in melt','FontSize',14)
%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);
......
This diff is collapsed.
......@@ -585,7 +585,7 @@ module initialize
cell%ref_a(:,:,0) = cell%ref_a(:,:,1) + 0.5_dp*(cell%ref_a(:,:,1)-cell%ref_a(:,:,2))
!cell%k(:,:,0) = (13.6_dp+1.88e-5_dp*pF%D) * ((300.0_dp-pF%Ts)/pF%DeltaT/cell%ref_T(1,1,0))**0.58_dp / pF%k ! Tosi et al 2013, pv + per, profile A
!cell%k(:,:,0) = (3.48_dp+5.17e-6_dp*pF%D) * (300.0_dp/(cell%ref_T(1,1,0)*pF%DeltaT+pF%Ts))**0.31_dp / pF%k ! Tosi et al 2013, pv + per, profile B
cell%k(:,:,0) = (3.48_dp+1.171e-10_dp*pI%D*pF%D*pF%g*pF%rho) * &
cell%k(:,:,0) = (3.48_dp+1.171e-10_dp*pI%D*pF%D*pF%g*pF%rho) * & !ToDo Check, why 1.171e-10 instead of 1.2e-10 as used above???
& (300.0_dp/(cell%ref_T(1,1,0)*pF%DeltaT+pF%Ts))**0.31_dp / pF%k ! p-dependence derived from Tosi paper using Earth values
cell%cp(:,:,0) = cell%cp(:,:,1) + 0.5_dp*(cell%cp(:,:,1)-cell%cp(:,:,2))
cell%shear(:,:,0) = cell%shear(:,:,1) + 0.5_dp*(cell%shear(:,:,1)-cell%shear(:,:,2))
......
......@@ -965,7 +965,7 @@ program CHIC
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="&
if(rank.eq.0) write(*,'("TS ",i8," 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,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,Te_save*pF%DeltaT+pF%Ts,Te2*pF%DeltaT+pF%Ts,Nu,Nu_t,Nu_m,Nu_b,ViscCon,VD,Work,ErrVD*100_dp,0.0_dp
......@@ -1011,13 +1011,13 @@ program CHIC
! write(*,'(str_out,i3)') i
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,&
write(*,'("TS ",i8," 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," (",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,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,&
write(*,'("TS ",i8," 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," (",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,&
......@@ -1130,7 +1130,7 @@ program CHIC
if (pE%output_dim.eq.1) then
if (rank.eq.0) then
write(*,'("TS ",i6," t=",ES10.3," dt=",ES10.3," Tr=",F7.2," Tcr=",F7.2,&
write(*,'("TS ",i8," t=",ES10.3," dt=",ES10.3," Tr=",F7.2," Tcr=",F7.2,&
& " Tl=",F7.2," Tm=",F7.2," Tb=",F8.2," Tc=",F8.2," Dr=",ES9.3," Dcr=",ES9.3," Dl=",ES9.3,&
& " Dut=",ES9.3," Dct=",ES9.3," Deff=",ES9.3)') i,t/pF%time_yr,dt/pF%time_yr,&
& tet%Tr*pF%DeltaT+pF%Ts,tet%Tcr*pF%DeltaT+pF%Ts,tet%Tl*pF%DeltaT+pF%Ts,tet%Tm*pF%DeltaT+pF%Ts,&
......@@ -1190,7 +1190,7 @@ program CHIC
endif
else if (pE%output_dim.eq.0) then
if (rank.eq.0) then
write(*,'("TS ",i6," t=",ES10.3," dt=",ES10.3," Tr=",F7.4," Tcr=",F7.4," Tl=",F7.4," Tm=",F7.4,&
write(*,'("TS ",i8," t=",ES10.3," dt=",ES10.3," Tr=",F7.4," Tcr=",F7.4," Tl=",F7.4," Tm=",F7.4,&
& " Tb=",F7.4," Tc=",F7.4," Dr=",F8.5," Dcr=",F8.5," Dl=",F8.5," Dut=",F8.5," Dct=",F8.5," Deff=",F8.5)') &
& i,t,dt,tet%Tr,tet%Tcr,tet%Tl,tet%Tm,tet%Tb,tet%Tc,tet%Dreg,tet%Dcr,tet%Dl,tet%Delta_u,tet%Delta_c,tet%Deff
......@@ -1237,7 +1237,7 @@ program CHIC
! & i,t,dt,pr%Rp/1000.0_dp,pr%Rc/1000.0_dp,pr%Tb(pR%nc),pr%Tc_av,pr%Tm_av,pr%Tw_av,pr%rho_av,pr%Dl/1000.0_dp,&
! & maxval(pr%v(pr%nc+1:pr%nc+pr%nm))*365.25_dp*24.0_dp*3600.0_dp*100.0_dp
! write(51,'(8(x e23.15,x) )') t,dt,pr%Rp,pr%Tb(pR%nc),pr%Tc_av,pr%Tm_av,pr%Tw_av,pr%rho_av
write(*,'("TS ",i6," t=",ES10.3," dt=",ES10.3," Rp",F9.2," Rc",F9.2," Tc",F9.2,&
write(*,'("TS ",i8," t=",ES10.3," dt=",ES10.3," Rp",F9.2," Rc",F9.2," Tc",F9.2,&
& " TcA",F9.2," TmA",F9.2," TwA",F9.2," rhoA",F10.2," Dl",F9.2," vmax",F10.3," qs",F9.3," qc",F9.3)') &
& i,t,dt,pr%Rp/1000.0_dp,pr%Rc/1000.0_dp,pr%Tb(pR%nc),pr%Tc_av,pr%Tm_av,pr%Tw_av,pr%rho_av,pr%Dl/1000.0_dp,&
& maxval(pr%v(pr%nc+1:pr%nc+pr%nm))*365.25_dp*24.0_dp*3600.0_dp*100.0_dp,&
......@@ -1323,7 +1323,7 @@ program CHIC
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="&
if(rank.eq.0)write(*,'("TS ",i8," 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
if ((pN%debug.ge.1).and.(rank.eq.0)) write(*,*) "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++"
......@@ -1334,13 +1334,13 @@ program CHIC
! &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,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,&
if(rank.eq.0) write(*,'("TS ",i8," 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," (",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,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,&
if(rank.eq.0) write(*,'("TS ",i8," 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," (",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,&
......@@ -1950,7 +1950,7 @@ program CHIC
timeS2=CLOCK() !omp_get_wtime()
if ((pN%debug.ge.2).and.(rank.eq.0)) write(*,*) "--------------------------------------------------------"
if (pE%output_dim.eq.0) then
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,&
if(rank.eq.0) write(*,'("TS ",i8," 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,nb_inner,t,dt,Vrms,V0rms,V0rms_max,T_mean,Te_save,Te2,Nu,Nu_t,Nu_m,Nu_b,ViscCon,VD,Work,&
......@@ -1980,7 +1980,7 @@ program CHIC
& * 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,&
write(*,'("TS ",i8," 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," (",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,&
......@@ -2180,7 +2180,7 @@ program CHIC
!write(*,*) timeS1,timeS2,timeS2-timeS1
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=",&
if(rank.eq.0) write(*,'("TS ",i8," 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,nb_inner,t,dt,Vrms,V0rms,V0rms_max,T_mean,Te_save,Te2,Nu,Nu_t,Nu_m,Nu_b,ViscCon,VD,Work,&
&ErrVD*100_dp,timeS2-timeS1
......@@ -2211,13 +2211,13 @@ program CHIC
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,&
if(rank.eq.0) write(*,'("TS ",i8," 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," (",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,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,&
if(rank.eq.0) write(*,'("TS ",i8," 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," (",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,&
......@@ -2597,7 +2597,7 @@ program CHIC
if (pE%output_dim.eq.1) then
if (rank.eq.0) then
write(*,'("TS ",i6," t=",ES10.3," dt=",ES10.3," Tr=",F7.2," Tcr=",F7.2,&
write(*,'("TS ",i8," t=",ES10.3," dt=",ES10.3," Tr=",F7.2," Tcr=",F7.2,&
& " Tl=",F7.2," Tm=",F7.2," Tb=",F8.2," Tc=",F8.2," Dr=",ES9.3," Dcr=",ES9.3," Dl=",ES9.3,&
& " Dut=",ES9.3," Dct=",ES9.3," Deff=",ES9.3)') i,t/pF%time_yr,dt/pF%time_yr,&
& tet%Tr*pF%DeltaT+pF%Ts,tet%Tcr*pF%DeltaT+pF%Ts,tet%Tl*pF%DeltaT+pF%Ts,tet%Tm*pF%DeltaT+pF%Ts,&
......@@ -2648,7 +2648,7 @@ program CHIC
endif
else if (pE%output_dim.eq.0) then
if (rank.eq.0) then
write(*,'("TS ",i6," t=",ES10.3," dt=",ES10.3," Tr=",F7.4," Tcr=",F7.4," Tl=",F7.4," Tm=",F7.4,&
write(*,'("TS ",i8," t=",ES10.3," dt=",ES10.3," Tr=",F7.4," Tcr=",F7.4," Tl=",F7.4," Tm=",F7.4,&
& " Tb=",F7.4," Tc=",F7.4," Dr=",F8.5," Dcr=",F8.5," Dl=",F8.5," Dut=",F8.5," Dct=",F8.5," Deff=",F8.5)') &
& i,t,dt,tet%Tr,tet%Tcr,tet%Tl,tet%Tm,tet%Tb,tet%Tc,tet%Dreg,tet%Dcr,tet%Dl,tet%Delta_u,tet%Delta_c,tet%Deff
......@@ -2722,7 +2722,7 @@ program CHIC
call Run1DProf(pr,dt,t,pN%debug,pG%nr,pE,pI,pV,pT%Ttop,pT%Tbottom)
write(*,'("TS ",i6," t=",ES10.3," dt=",ES10.3," Rp",F9.2," Rc",F9.2," Tc",F9.2,&
write(*,'("TS ",i8," t=",ES10.3," dt=",ES10.3," Rp",F9.2," Rc",F9.2," Tc",F9.2,&
& " TcA",F9.2," TmA",F9.2," TwA",F9.2," rhoA",F10.2," Dl",F9.2," vmax",F10.3," qs",F9.3," qc",F9.3)') &
& i,t,dt,pr%Rp/1000.0_dp,pr%Rc/1000.0_dp,pr%Tb(pR%nc),pr%Tc_av,pr%Tm_av,pr%Tw_av,pr%rho_av,pr%Dl/1000.0_dp,&
& maxval(pr%v(pr%nc+1:pr%nc+pr%nm))*365.25_dp*24.0_dp*3600.0_dp*100.0_dp,&
......
......@@ -114,7 +114,7 @@ contains
end subroutine
subroutine melt_comp_p(field,cell,mesh,mc,part,pF,pM)
subroutine melt_comp_p(field,cell,mesh,mc,part,pF,pM,pV,pX)
type(variables_unknowns),intent(inout)::field ! out: if T>solidus in crust, then cut down to avoid numerical problems
type(cell_properties),intent(in)::cell
type(mesh_cp),intent(in)::mesh
......@@ -122,9 +122,11 @@ contains
type(particle),intent(inout)::part
type(param_F),intent(in)::pF
type(param_M),intent(in)::pM
type(param_V),intent(in)::pV
type(param_X),intent(in)::pX
integer::i,j,k,m
real(dp)::z,solidus,liquidus,solidus_ref
real(dp)::z,solidus,liquidus,solidus_ref,k_solid,k_melt
! change depletion, TODO later on: add change in heat sources etc. (maybe in extra function like dehydration?)
......@@ -184,6 +186,8 @@ if ((z.lt.pM%melt_depth).and.(field%T(i,j,k).gt.mc%solidus_ref(i,j,k)-solidus_hy
! part%df(m) = pM%max_depl-part%d(m)
! endif
if (pM%max_depl.ne.0.0_dp) part%df(m) = min( part%df(m),pM%max_depl-part%d(m) )
! in any case, df cannot be larger than 100%:
part%df(m) = min( part%df(m),1.0_dp )
if (pM%melt_extract.eq.1) part%df(m) = max(part%df(m),0.0_dp) ! no melt kept
......@@ -195,8 +199,32 @@ if ((z.lt.pM%melt_depth).and.(field%T(i,j,k).gt.mc%solidus_ref(i,j,k)-solidus_hy
! field%T(i,j,k) = solidus ! cut temp to solidus (though not yet crust solidus) to avoid unrealistic high temps in crust
!endif
part%d(m) = max(part%d(m) + part%df(m),0.0_dp) ! shoudln't go under 0.0 anyway, but just in case
if (pM%melt_extract.eq.1) part%d(m) = max(part%d(m) + part%df(m),0.0_dp) ! shoudln't go under 0.0 anyway, but just in case
if (pM%melt_extract.eq.0) part%d(m) = max(part%d(m),part%df(m)) ! only relevant if not extrating melt and looking only at what melt fraction
! one would have if considering melting (benchmarks, magma oceans etc.)
if (pM%melt_extract.eq.0) then ! add how melt would change local properties:
part%c(m) = 1.0_dp - (1.0_dp-pX%density_melt)*part%df(m) ! density reduced if molten; take surface density here; scales with ref_r with depth
! conductivity: following Golabek and Solomatov papers, an effective k can be estimated for magma oceans
! however, here we may have only local magma ponds, so cannot use that approach
! in Golabek et al 2014, Meteoritics, it is writeen that k_eff <= 10^6 W/mK; use here 1000 (since smaller ponds) as upper k limit
! needs to be refined in future studies!!!
k_solid = (3.48_dp+1.2e-10_dp*z*pF%D*pF%g*pF%rho) &
& * (300.0_dp/(field%T(i,j,k)*pF%DeltaT+pF%Ts))**0.31_dp / pF%k ! Tosi for pv/ppv actually
k_melt = 1000.0_dp / pF%k !ToDo improve; read-in k_melt; calculated depending on magma pond size...
! exp(ln(k_s)+DF*(ln(k_l)-ln(k_s)))
! -> DF=0: exp(ln(k_s))=k_s
! -> DF=1: exp(ln(k_s)+ln(k_l)-ln(k_s))=exp(ln(k_l))=k_l
part%k(m) = exp(log(k_solid) + part%df(m)*(log(k_melt)-log(k_solid)))
! reduce viscosity locally; note that lower-mantle-visc / visc prefactor is set in A_dif/A_dis; water is set directly in visc, so v is normally always one and can be changed here
! TODO: makes no sense, for df=0 visc is higher than solid visc...
! check instead Schaefer paper (what Patrick also uses in his model)
!part%v(m) = 1.0_dp*exp((1.0_dp-part%df(m))*(2.5_dp+((1-part%df(m))/part%df(m))^0.48_dp)) ! Formular from Golabek et al., 2014
part%v(m) = 1.0_dp*exp(-part%df(m)*(2.5_dp+(part%df(m)/(1.0_dp-part%df(m)))**0.48_dp)) ! Formular adapted from Golabek et al., 2014 -> v is 1 for df=0; 0.77 for df=0.1; 0.03 for df=0.8 and 0.002 for df=0.95
part%v(m) = max(part%v(m),pV%eta_min)
endif
else
! two-phase flow benchmark
! assume here that liquidus temperature increases as well as solidus temperature!!!
......
......@@ -1211,7 +1211,7 @@ contains
if (pr%Tmb(nc_old)-pr%Tb(nc_old).gt.0.0_dp) then
!write(*,*) "Temperature jump at CMB=",abs(pE%DTc)*(pr%Tmb(k)-pr%Tb(k)),"K",", Tm=",pr%Tmb(k),", Tad=",pr%Tb(k)
pr%DTc = abs(pE%DTc)*(pr%Tmb(nc_old)-pr%Tb(nc_old)) ! use nc_old instead of pr%nc, since otherwise if core index changed values can be core instead of mantle values
pr%Tb(k-1)=pr%Tb(k) + abs(pE%DTc)*(pr%Tmb(nc_old)-pr%Tb(nc_old)) ! if DTc=-1 then temperature at CMB is set to melting temperature of peridotite
pr%Tb(k-1)=pr%Tb(k) + pr%DTc ! if DTc=-1 then temperature at CMB is set to melting temperature of peridotite
!write(*,*) "TCMB=",abs(pE%DTc)*(pr%Tmb(nc_old)-pr%Tb(nc_old)),k-1,k,nc_old,pr%nc,&
! &", nc_old:",pr%Tmb(nc_old),pr%Tb(nc_old),", nc_new:",pr%Tmb(pr%nc),pr%Tb(pr%nc)
else
......@@ -1221,9 +1221,13 @@ contains
! add linear temp increase in TBL above CMB, overwrite temperature:
if (pE%DTc.ne.0) then
do i=k,lTBL_k-1 ! from first shell within TBL to last shell above CMB
pr%Tb(i) = pr%Tb(lTBL_k) + (pr%Tb(k-1)-pr%Tb(lTBL_k))/pE%Delta_c*(pr%rb(lTBL_k)-pr%rb(i)) ! k-1 -> 1st shell in core, is set as T_CMB
!pr%Tb(i) = pr%Tb(lTBL_k) + (pr%Tb(k-1)-pr%Tb(lTBL_k))/pE%Delta_c*(pr%rb(lTBL_k)-pr%rb(i)) ! k-1 -> 1st shell in core, is set as T_CMB
pr%Tb(i) = pr%Tb(lTBL_k) + (pr%Tb(k-1)-pr%Tb(lTBL_k))/(pr%rb(lTBL_k)-pr%rb(k-1))*(pr%rb(lTBL_k)-pr%rb(i))
enddo
endif
! write(*,*) " T prof:",pr%Tb(pr%nc),pr%Tb(pr%nc+1),pr%Tb(pr%nc+2),pr%Tb(pr%nc+3),pr%Tb(pr%nc+4)
! write(*,*) " Ind:",lTBL_k,nc_old,pr%nc+1,lTBL_k-1,pE%DTc,pr%DTc
! write(*,*) "Calc:",pr%Tb(k-1),pr%Tb(k),pr%Tmb(nc_old-1)-pr%Tb(nc_old-1),pr%Tmb(nc_old)-pr%Tb(nc_old)
endif
enddo
......
......@@ -1892,7 +1892,7 @@ contains
Ra_depl = pM%Ra_depl
endif
if (pM%update_solidus.eq.1) then
if (pM%update_solidus.eq.1) then ! ToDo: should depend on both cells with indices "k" and "k-1"!!!
! ToDo: no refR here? then also no phase transition influence?
BoussC = -pM%Ra_melt*mc%porosity(i,j,k) - Ra_depl*cell%d(i,j,k)*(1.0_dp-mc%porosity(i,j,k)) ! TODO Check: is also mutliplied below with Ra, correct? How is Ra_depl defined?
else ! depletion not correct sind no update soliuds -> use DeltaF here instead
......@@ -1905,10 +1905,8 @@ contains
! & + refRu*g_u*FPu*(Cu-C_ref)*Du*mesh%dVi(i,j,k-1) )*dVi ! prescribed chemical density via C or depletion via D
! BoussC = B*( refR*g_t*FPt*(C-C_ref)*(1.0_dp-D_t*pI%drho)*mesh%dVi(i,j,k) &
! & + refRu*g_u*FPu*(Cu-C_ref)*(1.0_dp-Du*pI%drho)*mesh%dVi(i,j,k-1) )*dVi ! prescribed chemical density via C or depletion via D
BoussC = B*( refR*g_t*FPt*(C-C_ref-pI%drho*D_t)*mesh%dVi(i,j,k) &
BoussC = B*( refR*g_t*FPt*(C-C_ref-pI%drho*D_t)*mesh%dVi(i,j,k) &
& + refRu*g_u*FPu*(Cu-C_ref-pI%drho*Du)*mesh%dVi(i,j,k-1) )*dVi ! prescribed chemical density via C or depletion via D
if (pN%NoMMSolveA) then
!BoussC = -B*pI%drho*((C-C_ref)*D_t*mesh%dVi(i,j,k)+(Cu-C_ref)*Du*mesh%dVi(i,j,k-1))*dVi
!BoussC = B*pI%drho*((C-C_ref)*(1.0_dp-D_t)*mesh%dVi(i,j,k)+(Cu-C_ref)*(1.0_dp-Du)*mesh%dVi(i,j,k-1))*dVi
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment