Skip to content
GitLab
Commit dd7adfad authored by Lena Noack's avatar Lena Noack
Browse files

Some updates

parent 1d715dc6
Expand all Hide whitespace changes
Inline Side-by-side
Octave/make_movie.m
View file @ dd7adfad
%function make_movie(start=0,fpath=' ',vals=' ') for normal plots
function make_movie(start=0,fpath='D:\Curta_PhaseTrans_Proxima\V5_Case3\',vals='TWHDvVRMt')
%function make_movie(start=0,fpath='D:\Curta_PhaseTrans_Proxima\V5_Case3\',vals='TWHDvVRMt')
function make_movie(start=0,fpath='/scratch-3/Simulations_New/Project_Venus/Ref_MOprofs_Highres_PhaseTrans_NoChemBuoy/',vals='TWHDvVRMt')
%function make_movie(start=0,fpath='/scratch-3/Simulations_New/Project_Venus/Ref_MOprofs_Highres_PhaseTrans_HotSurf1000K/',vals='TWHDvVRMt')
if fpath==' '
[fname, fpath, fltidx] = uigetfile("data_char_val_ts.res");
......
Octave/plot_snap.m
View file @ dd7adfad
......@@ -4,8 +4,8 @@
function plot_snap(max_depl=30,max_T=0,fullV=-1,part=0,reduc=1,fin=' ',print_plot='TWHDvVRMtOo')%' ')
% if reduc=1, then plot only T2D, V2D, Vi2D, D2D, Mradp2D
%addpath("/home/noacklen/Arbeit/SOURCE/CHIC/Octave")
addpath("C:/Arbeit/Git_Source/CHIC/Octave")
addpath("/home/noacklen/Arbeit/SOURCE/CHIC/Octave")
%addpath("C:/Arbeit/Git_Source/CHIC/Octave")
%addpath("C:/Arbeit/NextCloud_FUB/CHIC/Code_Backup_March2020/Octave")
[data,input,program_dir,fname,fpath] = read_snap(fin); % read data from snapshot
......
Python/Code_Nick_Wogan/Volcano-Speciation-master/Explanation_of_model.ipynb
View file @ dd7adfad
......@@ -508,7 +508,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.4"
"version": "3.7.3"
}
},
"nbformat": 4,
......
Python/Code_Nick_Wogan/Volcano-Speciation-master/Several_plots.ipynb
View file @ dd7adfad
......@@ -815,7 +815,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.4"
"version": "3.7.3"
}
},
"nbformat": 4,
......
src/heating.f90
View file @ dd7adfad
......@@ -19,11 +19,11 @@ contains
real(dp)::HtU238,HtU235,HtTh232,HtK40
if (pI%C_U.eq.0.0_dp) then
radiogenic_heating = H0*exp(-lambda*t)
radiogenic_heating = H0*exp(-lambda*t) ! H0 is in nondimensional W/kg; hence radiogenic_heating is in ND W/kg
else
! calculate from C_U, C_Th and C_K, lambda is 1 (cooling) or 0 (constant heat sources)
HtU238 = pI%H0U238*Exp(-Log(2.0_dp)*t*lambda/pI%t12U238)
HtU238 = pI%H0U238*Exp(-Log(2.0_dp)*t*lambda/pI%t12U238) ! in ND W/kg
HtU235 = pI%H0U235*Exp(-Log(2.0_dp)*t*lambda/pI%t12U235)
HtTh232 = pI%H0Th232*Exp(-Log(2.0_dp)*t*lambda/pI%t12Th232)
HtK40 = pI%H0K40*Exp(-Log(2.0_dp)*t*lambda/pI%t12K40)
......
src/initialize.f90
View file @ dd7adfad
......@@ -51,13 +51,14 @@ module initialize
end subroutine initialize_pressure_strR
subroutine initialize_cell(cell,field,mesh,pE,pF,pI,pT,pV,pX,rho_av,max_depl,rc,compress,&
subroutine initialize_cell(cell,field,mesh,pE,pF,pH,pI,pT,pV,pX,rho_av,max_depl,rc,compress,&
& Di,DiT,Gr,T0,rmin,TinclTref,Tini,Tbot,shear,DeltaVisc,DeltaTc,MOprofs)
type(cell_properties),intent(inout)::cell
type(variables_unknowns),intent(inout)::field
type(mesh_cp),intent(in)::mesh
type(param_E),intent(in)::pE
type(param_F),intent(in)::pF
type(param_H),intent(inout)::pH
type(param_I),intent(inout)::pI
type(param_T),intent(inout)::pT
type(param_V),intent(in)::pV
......@@ -70,11 +71,12 @@ module initialize
real(dp),intent(inout)::Tini,Tbot ! adapt Tbot depending on initial upper mantle temperature and read-in adiabatic profile
real(dp),intent(in)::DeltaVisc,DeltaTc
integer::k,kinv,kp,nl,ny,nr,jmin,jmax,n_t,n_d
integer::i,j,k,kinv,kp,nl,ny,nr,jmin,jmax,n_t,n_d,ind
real(dp)::z,Cp_count,Core_count,vol,T_CMB
real(dp),allocatable::vals_t(:,:),vals_d(:,:),k_e(:)
integer,allocatable::mat(:)
real(dp)::gc,pi_,dr,g,rho,rb,Vav,check
real(dp)::M_MgFeO, M_MgFeSiO3
pi_ = acos(-1.0_dp)
......@@ -91,6 +93,21 @@ module initialize
endif
cell%c = pI%rho ! composition density (e.g. by melt)
! mass fractions of Mg/Fe O and Mg/Fe Si O3 in lower mantle for phonon conductivity formula
M_MgFeO = 0.2865_dp*(1.0_dp-pI%X_FeM)+0.3526_dp*pI%X_FeM
M_MgFeSiO3 = 0.7135_dp*(1.0_dp-pI%X_FeM)+0.6474_dp*pI%X_FeM
ind=1
if (size(pI%image,1)>1) then
do i=1,nl
do k=1,nr
cell%c(i,1,nr+1-k) = real(pI%image(ind))
ind=ind+1
enddo
enddo
endif
!cell%w = pI%mantle_w_H2O
cell%wr = pI%mantle_w_r ! water exponent for viscosity
cell%CO2 = pX%CO2_ini
......@@ -480,6 +497,11 @@ module initialize
deallocate( vals_t )
! tidal heating needs planet-averaged density, if zero as input, then set here:
if (pH%rho_av.eq.0.0_dp) then
pH%rho_av = rho_av
endif
!profiles produced by 1D CHIC model
else if (pI%read_profs.ge.4) then
n_t = 1000
......@@ -572,8 +594,13 @@ module initialize
cell%ref_T(:,:,nr+1) = vals_t(n_t,5)
cell%ref_a(:,:,nr+1) = vals_t(n_t,8)
cell%cp(:,:,nr+1) = vals_t(n_t,7)
! cell%k(:,:,nr+1) = 13.6_dp * ((300.0_dp-pF%Ts)/pF%DeltaT/vals_t(n_t,5))**0.58_dp / pF%k ! Tosi et al 2013, pv + per, profile A
cell%k(:,:,nr+1) = 3.48_dp * (300.0_dp/(vals_t(n_t,5)*pF%DeltaT+pF%Ts))**0.31_dp / pF%k ! Tosi et al 2013, pv + per, profile B
! cond. at surface -> use Tosi model, since Stamenkovic only for lower mantle!
!if (pI%k_prof.eq.1) then
! cell%k(:,:,nr+1) = 13.6_dp * ((300.0_dp-pF%Ts)/pF%DeltaT/vals_t(n_t,5))**0.58_dp / pF%k ! Tosi et al 2013, pv + per, profile A
cell%k(:,:,nr+1) = 3.48_dp * (300.0_dp/(vals_t(n_t,5)*pF%DeltaT+pF%Ts))**0.31_dp / pF%k ! Tosi et al 2013, pv + per, profile B
!elseif (pI%k_prof.eq.2) then
! cell%k(:,:,nr+1) = 8.5e-11_dp*(vals_t(n_t,5)*pF%DeltaT+pF%Ts)**3 / pF%k ! Stamenkovic et al., 2011; Hofmeister, 1999 for radiative conductivity
!endif
cell%bulk(:,:,nr+1) = vals_t(n_t,10)
cell%shear(:,:,nr+1) = vals_t(n_t,12)
if (pI%read_profs.ge.5) k_e(nr+1) = vals_t(n_t,13)
......@@ -584,40 +611,57 @@ module initialize
z = 0.5_dp*(rc(nr)+rc(nr+1)) - rc(k)
do kp=1,n_t-1
if ((vals_t(kp,4).le.rc(k)).and.(vals_t(kp+1,4).ge.rc(k))) then
if (pI%read_profs.ge.5) mat(k)=vals_t(kp+1,14)
cell%ref_g(:,:,k) = vals_t(kp+1,1)-(vals_t(kp+1,1)-vals_t(kp,1))*(vals_t(kp+1,4)-rc(k))/(vals_t(kp+1,4)-vals_t(kp,4))
cell%ref_p(:,:,k) = vals_t(kp+1,2)-(vals_t(kp+1,2)-vals_t(kp,2))*(vals_t(kp+1,4)-rc(k))/(vals_t(kp+1,4)-vals_t(kp,4))
cell%ref_r(:,:,k) = vals_t(kp+1,3)-(vals_t(kp+1,3)-vals_t(kp,3))*(vals_t(kp+1,4)-rc(k))/(vals_t(kp+1,4)-vals_t(kp,4))
cell%ref_T(:,:,k) = vals_t(kp+1,5)-(vals_t(kp+1,5)-vals_t(kp,5))*(vals_t(kp+1,4)-rc(k))/(vals_t(kp+1,4)-vals_t(kp,4))
cell%ref_a(:,:,k) = vals_t(kp+1,8)-(vals_t(kp+1,8)-vals_t(kp,8))*(vals_t(kp+1,4)-rc(k))/(vals_t(kp+1,4)-vals_t(kp,4))
cell%cp(:,:,k) = vals_t(kp+1,7)-(vals_t(kp+1,7)-vals_t(kp,7))*(vals_t(kp+1,4)-rc(k))/(vals_t(kp+1,4)-vals_t(kp,4))
! cell%k(:,:,k) = (13.6_dp+1.88e-5_dp*z*pF%D) * ((300.0_dp-pF%Ts)/pF%DeltaT/cell%ref_T(1,1,k))**0.58_dp / pF%k
! cell%k(:,:,k) = (3.48_dp+5.17e-6_dp*z*pF%D) * (300.0_dp/(cell%ref_T(1,1,k)*pF%DeltaT+pF%Ts))**0.31_dp / pF%k
cell%k(:,:,k) = (3.48_dp+1.2e-10_dp*vals_t(kp+1,2)*pF%D*pF%g*pF%rho) &
& * (300.0_dp/(cell%ref_T(1,1,k)*pF%DeltaT+pF%Ts))**0.31_dp / pF%k
if ((pI%k_prof.eq.1).or.((pI%k_prof.eq.2).and.(mat(k).le.5))) then ! for Tosi model & for upper mantle Stamenkovic model
! cell%k(:,:,k) = (13.6_dp+1.88e-5_dp*z*pF%D) * ((300.0_dp-pF%Ts)/pF%DeltaT/cell%ref_T(1,1,k))**0.58_dp / pF%k
! cell%k(:,:,k) = (3.48_dp+5.17e-6_dp*z*pF%D) * (300.0_dp/(cell%ref_T(1,1,k)*pF%DeltaT+pF%Ts))**0.31_dp / pF%k
cell%k(:,:,k) = (3.48_dp+1.2e-10_dp*vals_t(kp+1,2)*pF%D*pF%g*pF%rho) &
& * (300.0_dp/(cell%ref_T(1,1,k)*pF%DeltaT+pF%Ts))**0.31_dp / pF%k
elseif (pI%k_prof.eq.2) then
! only used for lower mantle, otherwise Tosi formula is used
cell%k(:,:,k) = ( M_MgFeO*7.01_dp*(cell%ref_r(1,1,k)*pF%rho/3300.0_dp)**1.94_dp&
& *(2000.0_dp/(cell%ref_T(1,1,k)*pF%DeltaT+pF%Ts)) & ! exp 1.94 comes from: 3*0.758+2*0-1/3 in infinite pressure limit (0.758=K_inf/2-1/6)
&+ M_MgFeSiO3*0.71_dp*(cell%ref_r(1,1,k)*pF%rho/3870.0_dp)**3.12_dp&
& *(2000.0_dp/(cell%ref_T(1,1,k)*pF%DeltaT+pF%Ts)) & ! exp 3.12 comes from: 3*1.15+2*0-1/3 in infinite pressure limit
&+ 8.5e-11_dp*(cell%ref_T(1,1,k)*pF%DeltaT+pF%Ts)**3) / pF%k ! Stamenkovic et al., 2011; Hofmeister, 1999 for radiative conductivity
endif
cell%bulk(:,:,k)=vals_t(kp+1,10)-(vals_t(kp+1,10)-vals_t(kp,10))*(vals_t(kp+1,4)-rc(k))/(vals_t(kp+1,4)-vals_t(kp,4))
cell%shear(:,:,k)=vals_t(kp+1,12)-(vals_t(kp+1,12)-vals_t(kp,12))*(vals_t(kp+1,4)-rc(k))/(vals_t(kp+1,4)-vals_t(kp,4))
if (pI%read_profs.ge.5) k_e(k)=vals_t(kp+1,13)-(vals_t(kp+1,13)-vals_t(kp,13))*(vals_t(kp+1,4)-rc(k))/(vals_t(kp+1,4)-vals_t(kp,4))
if (pI%read_profs.ge.5) mat(k)=vals_t(kp+1,14)
endif
enddo
enddo
! boundary values at CMB might be wrong since rho etc interpolated
! between core and mantle values, take next mantle value instead:
if (pI%read_profs.ge.5) mat(0) = mat(1)
cell%ref_g(:,:,0) = cell%ref_g(:,:,1) + 0.5_dp*(cell%ref_g(:,:,1)-cell%ref_g(:,:,2))
cell%ref_p(:,:,0) = cell%ref_p(:,:,1) + 0.5_dp*(cell%ref_p(:,:,1)-cell%ref_p(:,:,2))
cell%ref_T(:,:,0) = cell%ref_T(:,:,1) + 0.5_dp*(cell%ref_T(:,:,1)-cell%ref_T(:,:,2))
cell%ref_r(:,:,0) = cell%ref_r(:,:,1) + 0.5_dp*(cell%ref_r(:,:,1)-cell%ref_r(:,:,2))
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) * & !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
if ((pI%k_prof.eq.1).or.((pI%k_prof.eq.2).and.(mat(0).le.5))) then
!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) * & !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
elseif (pI%k_prof.eq.2) then
cell%k(:,:,0) = ( M_MgFeO*7.01_dp*(cell%ref_r(1,1,0)*pF%rho/3300.0_dp)**1.94_dp&
& *(2000.0_dp/(cell%ref_T(1,1,0)*pF%DeltaT+pF%Ts)) &
&+ M_MgFeSiO3*0.71_dp*(cell%ref_r(1,1,0)*pF%rho/3870.0_dp)**3.12_dp&
& *(2000.0_dp/(cell%ref_T(1,1,0)*pF%DeltaT+pF%Ts)) &
&+ 8.5e-11_dp*(cell%ref_T(1,1,0)*pF%DeltaT+pF%Ts)**3 ) / pF%k ! Stamenkovic et al., 2011; Hofmeister, 1999 for radiative conductivity
endif
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))
cell%bulk(:,:,0) = cell%bulk(:,:,1) + 0.5_dp*(cell%bulk(:,:,1)-cell%bulk(:,:,2))
if (pI%read_profs.ge.5) k_e(0) = k_e(1) + 0.5_dp*(k_e(1)-k_e(2))
if (pI%read_profs.ge.5) mat(0) = mat(1)
write(*,*) " Radius (km) Ref. Temp. (K) Density (kg/m^3) &
& Alpha (1e-5 1/K) Heat capacity (J/kgK) Therm. Cond. (W/mK) El. Cond. (log10 S/M) Layer"
......@@ -741,7 +785,7 @@ module initialize
nr = mesh%nr
! specific heating rates in W/kg
! specific heating rates in nondimensional W/kg
pI%H0U238 = 9.46e-5 * pF%D**2/(pF%DeltaT * pF%k) * pF%rho
pI%H0U235 = 5.69e-4 * pF%D**2/(pF%DeltaT * pF%k) * pF%rho
pI%H0Th232 = 2.64e-5 * pF%D**2/(pF%DeltaT * pF%k) * pF%rho
......@@ -991,7 +1035,8 @@ module initialize
& ((0.5_dp*(rc(k)+rc(k+1)))**3-(0.5_dp*(rc(k)+rc(k-1)))**3)
! Q here multiplied with dimensionless ref_r and radii, to account for variations
! with depth -> in main.f90 Q_eddy is then calculated by multiplying with D, rho_ref etc.
! only question: should it be time rho or not? Want to have pW/kg???
! only question: should it be time rho or not? Want to have pW/kg???
! -> vec_magn_H here stored in nondimensional W/kg, but still needs to be multiplied with local ND rho in thermal.f90
enddo
if (pH%use_magn_H.ge.4.0_dp) Q = Q/2.0_dp ! actually is angle-dependent, but divide here by 2 to get correct mantle-averaged Q
pH%Q_eddy = Q
......
src/input.f90
View file @ dd7adfad
......@@ -111,6 +111,7 @@ contains
pI%k=val
case("eta_ref") ! use reference viscosity to calculate Rayleigh number (Ra is not used then)
pF%eta_ref=val
pI%eta_ref=val
case("OS_SuL") ! additional surface layer on top of domain
if (val.gt.0.0_dp) then
pG%rmax = pG%rmax+val
......@@ -148,6 +149,7 @@ contains
pG%rmax=val
case("eta_ref") ! use reference viscosity to calculate Rayleigh number (Ra is not used then)
pF%eta_ref=val
pI%eta_ref=val
case("crust_ini")
pC%crust_ini=val
case("add_crust")
......@@ -414,6 +416,7 @@ contains
pI%rho = pI%rho / pF%rho
pI%kappa = pI%kappa / pF%kappa
if (pI%C_U.eq.0.0_dp) pI%H0 = pI%H0 * pF%D**2 / (pF%DeltaT * pF%k) ! otherwise already non-dimensional enrichment/depletion factor
! pI%H0 already multiplied with pF%rho above -> H0 is now in nondimensionalized W/kg
pI%k = pI%k / pF%k
pV%R = pV%R / pF%R
......@@ -470,7 +473,7 @@ contains
case("k")
if (rank.eq.0) write(*,*) var,"=",val," -> ",pI%k
case("eta_ref")
if (rank.eq.0) write(*,*) var,"=",val," -> ",pF%eta_ref
if (rank.eq.0) write(*,*) var,"=",val," -> ",pI%eta_ref
case("C_U")
if (rank.eq.0) write(*,*) var,"=",val," -> ",pI%C_U
case("C_Th")
......@@ -846,6 +849,12 @@ contains
case("MOprofs")
pI%MOprofs = int(val)
if (rank.eq.0) write(*,*) var,"=",pI%MOprofs
case("image") ! file to be read-in
allocate(pI%image(pG%nl*pG%nr))
OPEN(1,FILE='Image.txt',form='formatted',status='unknown')
READ(1,*) pI%image(:)
CLOSE(1)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! H: Heating parameters !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
......@@ -1315,7 +1324,9 @@ contains
if (rank.eq.0) write(*,*) var,"=",pV%rheol_vec
case("DeltaVisc")
pV%DeltaVisc = val
pI%eta_ref = val
if (rank.eq.0) write(*,*) var,"=",pV%DeltaVisc
if (rank.eq.0) write(*,*) "eta_ref -> ",pI%eta_ref
case("DeltaVisc_LM")
pV%DeltaVisc_LM = val
if (rank.eq.0) write(*,*) var,"=",pV%DeltaVisc_LM
......@@ -3138,6 +3149,9 @@ contains
case("X_H2OM") ! values between 0 and 100, in wt-% (so far only for electrical conductivity profile in model1d.f90)
pI%X_H2OM = val*0.01_dp
if (rank.eq.0) write(*,*) var,"=",pI%X_H2OM
case("k_prof")
pI%k_prof = int(val)
if (rank.eq.0) write(*,*) var,"=",pI%k_prof
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! X: concentrations and outgassing !!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
......@@ -3246,7 +3260,9 @@ contains
if (rank.eq.0) write(*,*) var,"=",val," -> ",pT%s_top
case("ExpDeltaVisc")
pV%DeltaVisc = 10.0_dp**val
pI%eta_ref = pV%DeltaVisc
if (rank.eq.0) write(*,*) var,"=",pV%DeltaVisc
if (rank.eq.0) write(*,*) "eta_ref -> ",pI%eta_ref
case("DeltaTc")
pT%DeltaTc=val/pF%DeltaT
if (rank.eq.0) write(*,*) var,"=",val," -> ",pT%DeltaTc
......
src/main.f90
View file @ dd7adfad
......@@ -116,6 +116,7 @@ program CHIC
&vol_all,FS_min,FS_max,geom_d,relax,delta_pr,delta_pr_old,relax_v,divvel,div,div_save,To_save,T_lid_m,T_upm
integer::magma_local,magma_global
real(dp)::magma !time when first post-accretional magma ocean appears
real(dp)::pi_
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,heat_av
character(len=10):: str_i
......@@ -142,6 +143,7 @@ program CHIC
big_number=10D99
time1=CLOCK() !omp_get_wtime()
depl_av=0.0_dp
pi_ = acos(-1.0_dp)
call SetTimeParallel()
......@@ -282,18 +284,25 @@ program CHIC
if (pE%IntStruct.lt.0) then
if (pE%IntStruct.eq.-1) then ! for melt paper, mat prop and core heat flux study
allocate( Masses(11) , XiFe(6) , XiH2O(4) )
Masses(1) = 0.8_dp !0.1_dp
Masses(2) = 0.9_dp !0.2_dp !0.5_dp
Masses(3) = 1.0_dp !0.4_dp !1.0_dp
Masses(4) = 1.1_dp !0.7_dp !1.5_dp
Masses(5) = 1.2_dp !1.0_dp !2.0_dp
Masses(6) = 1.3_dp !1.5_dp !2.5_dp
Masses(7) = 1.4_dp !2.0_dp !3.0_dp
Masses(8) = 1.5_dp !2.5_dp !3.5_dp
Masses(9) = 1.6_dp !3.0_dp !4.0_dp
Masses(10) = 1.8_dp !3.5_dp !4.5_dp
Masses(11) = 2.0_dp !4.0_dp !5.0_dp
allocate( Masses(10) , XiFe(7) , XiH2O(1) )
Masses(1) = 5.5_dp
Masses(2) = 6.0_dp
Masses(3) = 6.5_dp
Masses(4) = 7.0_dp
Masses(5) = 7.5_dp
Masses(6) = 8.0_dp
Masses(7) = 8.5_dp
Masses(8) = 9.0_dp
Masses(9) = 9.5_dp
Masses(10) = 10.0_dp
!Masses(1) = 2.5_dp !0.1_dp
!Masses(2) = 3.0_dp !0.2_dp !0.5_dp
!Masses(3) = 3.5_dp !0.4_dp !1.0_dp
!Masses(4) = 4.0_dp !0.7_dp !1.5_dp
!Masses(5) = 4.5_dp !1.0_dp !2.0_dp
!Masses(6) = 5.0_dp !1.5_dp !2.5_dp
!Masses(1) = 1.7_dp
!Masses(2) = 1.9_dp
! Masses(1) = 0.1_dp
! Masses(2) = 0.2_dp !0.5_dp
! Masses(3) = 0.4_dp !1.0_dp
......@@ -306,24 +315,24 @@ program CHIC
! Masses(10) = 3.5_dp !4.5_dp
! Masses(11) = 4.0_dp !5.0_dp
XiFe(1) = 0.20_dp !0.15_dp
XiFe(2) = 0.30_dp !0.25_dp
XiFe(3) = 0.40_dp !0.35_dp
XiFe(4) = 0.50_dp !0.45_dp
XiFe(5) = 0.60_dp !0.55_dp
XiFe(6) = 0.70_dp !0.65_dp
XiFe(7) = 0.80_dp !0.75_dp
XiFe(1) = 0.15_dp
XiFe(2) = 0.25_dp
XiFe(3) = 0.35_dp
XiFe(4) = 0.45_dp
XiFe(5) = 0.55_dp
XiFe(6)= 0.65_dp
XiFe(7)= 0.75_dp
! XiH2O array is used here not for water but for iron content in mantle; this is iron number and not weight fraction (e.g. Fe# 0.4 -> iron mantle weight fraction=0.1346) -> X_Fe-X_FeM*(1-X_Fe)>0!!!
XiH2O(1) = 0.0_dp !0.05_dp
XiH2O(2) = 0.05_dp
XiH2O(3) = 0.1_dp !0.15_dp
XiH2O(4) = 0.15_dp !0.2_dp
XiH2O(5) = 0.2_dp
do ind_k=1,11
do ind_j=1,7
do ind_i=1,5
XiH2O(1) = 0.1_dp !0.0_dp !0.05_dp
!XiH2O(2) = 0.05_dp
!XiH2O(3) = 0.1_dp !0.15_dp
!XiH2O(4) = 0.15_dp !0.2_dp
!XiH2O(5) = 0.2_dp
do ind_k=1,10 !13
do ind_j=1,7 !14!7
do ind_i=1,1 !5
pE%M_E = Masses(ind_k)
pE%X_Fe = XiFe(ind_j)
pI%X_FeM = XiH2O(ind_i)
......@@ -528,7 +537,7 @@ program CHIC
endif
! first call needed for adiabatic temperature profile
call initialize_cell(cell,field,mesh,pE,pF,pI,pT,pV,pX,pH%rho_av,pM%max_depl,mesh%rc,pI%compress,pI%Di,pI%DiT,pI%Gr,pT%T0,&
call initialize_cell(cell,field,mesh,pE,pF,pH,pI,pT,pV,pX,pH%rho_av,pM%max_depl,mesh%rc,pI%compress,pI%Di,pI%DiT,pI%Gr,pT%T0,&
&pG%save_rmin,pI%TinclTref,pT%Tini,pT%Tbottom,pH%shear_modulus,pV%DeltaVisc,pT%DeltaTc,pI%MOprofs)
call initialize_temperature(jmin_2D,field,mesh%rc,field%T,pT%Tbottom,pT%Ttop,pG%inner_bound,pG%outer_bound,&
......@@ -537,7 +546,7 @@ program CHIC
! second call for re-setting temperature if needed with profs.res
if (pI%read_profs.ge.6) then
call initialize_cell(cell,field,mesh,pE,pF,pI,pT,pV,pX,pH%rho_av,pM%max_depl,mesh%rc,pI%compress,pI%Di,pI%DiT,pI%Gr,pT%T0,&
call initialize_cell(cell,field,mesh,pE,pF,pH,pI,pT,pV,pX,pH%rho_av,pM%max_depl,mesh%rc,pI%compress,pI%Di,pI%DiT,pI%Gr,pT%T0,&
&pG%save_rmin,pI%TinclTref,pT%Tini,pT%Tbottom,pH%shear_modulus,pV%DeltaVisc,pT%DeltaTc,pI%MOprofs)
endif
......@@ -1017,12 +1026,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,heat_m,heat_av,0.0_dp
else
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")') &
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, H_tidal=",ES10.3," 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*pI%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,&
&heat_m,heat_av,0.0_dp
&ViscCon,pH%Q_eddy*(pF%k*pF%DeltaT)/(pF%D**2*pF%rho*pI%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,heat_m,heat_av,&
&pH%tidalH_global_2/pI%rho/pF%rho*1.0e12_dp,0.0_dp
endif
endif
endif
......@@ -1342,11 +1352,12 @@ program CHIC
&ViscCon,ErrVD*100_dp,heat_m,heat_av,0.0_dp
else
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")') &
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," Q_eddy=",ES10.3," pW/kg, "," H_av=",F10.7," (",F10.7,") pW/kg, H_tidal=",ES10.3,&
&" 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*pI%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,&
&heat_m,heat_av,0.0_dp
&ViscCon,pH%Q_eddy*(pF%k*pF%DeltaT)/(pF%D**2*pF%rho*pI%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,heat_m,&
&heat_av,pH%tidalH_global_2/pI%rho/pF%rho*1.0e12_dp,0.0_dp
endif
endif
......@@ -2221,11 +2232,12 @@ program CHIC
&ViscCon,ErrVD*100_dp,heat_m,heat_av,timeS2-timeS1
else
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")') &
&" qs=",F8.4,", qc=",F8.4,", VC=",ES10.3," Q_eddy=",ES10.3," pW/kg, "," H_av=",F10.7," (",F10.7,") pW/kg, H_tidal=",ES10.3,&
&" 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*pI%rho*4/3*acos(-1.0_dp)*(pI%Rp**3-pI%Rc**3))*1.0e+12,heat_m,heat_av,&
&timeS2-timeS1
&pH%tidalH_global_2/pI%rho/pF%rho*1.0e12_dp,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")') &
......@@ -3139,6 +3151,7 @@ program CHIC
! ToDo pG%dim.eq.1
endif
deallocate( melt_local )
if (size(pI%image,1)>1) deallocate( pI%image )
if (rank.eq.0) write(*,*) " "
if (rank.eq.0) write(*,*) "+++ Simulation terminated correctly (deallocation was successful) +++"
......
src/model1D.f90
View file @ dd7adfad
......@@ -90,7 +90,7 @@ contains
M = pE%M_E * 5.9736e+24 ! Mass in g
pi_ = acos(-1.0_dp)
maxOI = 20 !100
maxOI = 100 !20 !100
! note: Get1DStruct() uses dimensional(!) values, calculation Run1DProf() uses either dimensional or non-dimensional values
......@@ -99,9 +99,9 @@ contains
pr%rho_m = 4457.0_dp * pE%M_E**0.1779_dp
pr%rho_w = 2120.8_dp * pE%M_E**0.2007_dp
pr%alpha_m = pI%alpha ! only needed if EOS does not work in first iteration of EOS
pr%Cp_m = pI%Cp
pr%rho_m = pI%rho
pr%alpha_m = 2.0e-5_dp !pI%alpha ! only needed if EOS does not work in first iteration of EOS
pr%Cp_m = 1200.0_dp !pI%Cp
!pr%rho_m = pI%rho
Mc = pE%X_Fe * M
Mw = pE%X_H2O * M
......@@ -157,6 +157,9 @@ contains
pr%Rm = ( Mm*3.0_dp / (4.0_dp*pi_*pr%rho_m) + pr%Rc**3 )**(1.0_dp/3.0_dp)
pr%Rp = ( Mw*3.0_dp / (4.0_dp*pi_*pr%rho_w) + pr%Rm**3 )**(1.0_dp/3.0_dp)
!write(*,*) M, Mc, Mm, pr%Rc, pr%Rm, pr%Rp
!write(*,*) pE%M_E,pr%rho_c,pr%rho_m,pr%rho_w
! set-up preliminary arrays for radius and density, needed by Get1DStruct(), later on: use values from last time-step
pr%rb(0) = 0.0_dp
pr%nc = 0
......@@ -267,10 +270,16 @@ contains
!if (pr%Cp(k).eq.0.0_dp) pr%Cp(k)=pr%Cp(k+1)
!if (pr%rho(k).eq.0.0_dp) pr%rho(k)=pr%rho(k+1)
if (INAN(pr%alpha(k))) then
pr%alpha(k) = pr%alpha_m
pr%Cp(k) = pr%Cp_m
pr%rho(k) = pr%rho_m
if (INAN(pr%alpha(k)).or.(pr%alpha(k).lt.0.0_dp)) then ! ToDo: only true if no water
if (k.lt.nr) then
pr%alpha(k) = pr%alpha(k+1)
pr%Cp(k) = pr%Cp(k+1)
pr%rho(k) = pr%rho(k+1)
else
pr%alpha(k) = pr%alpha_m
pr%Cp(k) = pr%Cp_m
pr%rho(k) = pr%rho_m
endif
endif
pr%gb(k-1) = pr%gb(k) - (4.0_dp*pi_*6.67384e-11_dp-2.0_dp*pr%gb(k)/pr%rb(k))*(pr%rb(k)-pr%rb(k-1))
......@@ -316,7 +325,7 @@ contains
enddo
endif
endif
if (pr%Tb(k-1).gt.2.0_dp*pr%Tb(k)) then
if ((pr%Tb(k-1).gt.2.0_dp*pr%Tb(k)).and.(k.lt.nr)) then
if (debug.gt.1) write(*,*) "Reset T from ",pr%Tb(k-1)," to ",2.0_dp*pr%Tb(k)
pr%Tb(k-1)=2.0_dp*pr%Tb(k) ! just ini T prof, if wrong thermod. prop then T prof can increase to infinity
pr%rho(k) = pr%rho(k+1)
......@@ -324,8 +333,9 @@ contains
pr%Cp(k) = pr%Cp(k+1)
endif
if (debug.ge.2) write(*,*) "r=",pr%rb(k-1)/1000.0_dp,pr%Tb(k-1),pr%alpha(k),pr%gb(k),pr%rb(k)-pr%rb(k-1),pr%Cp(k),pr%pb(k)
if (debug.ge.2) then
write(*,*) k,"r=",pr%rb(k-1)/1000.0_dp,pr%Tb(k-1),pr%alpha(k),pr%gb(k),pr%rb(k)-pr%rb(k-1),pr%Cp(k),pr%pb(k)*1.0e-9_dp
endif
enddo
if (pr%nc.eq.0) pr%Tb(0)=pr%Tb(1) ! if no core
......@@ -546,7 +556,7 @@ contains
! type(phase)::forsterite,MgPerovskite,MgPostPerovskite,periclase
! maxOI = 100 ! in input.txt?
tolOI = 1.0e-8_dp
tolOI = 1.0e-9_dp !1.0e-8_dp
error = 0
M = M_E * 5.9736e+24 ! Mass in g
......@@ -601,7 +611,9 @@ contains
endif
Mm = M - Mc - Mw
!Mm = Mm + 2.0_dp*X_FeM*55.845_dp/(X_FeM*63.08_dp+140.69_dp)*Mm
!write(*,*) "Masses: ",M,Mc,Mw
! store updated values:
pI%X_FeM = X_FeM
pI%X_H2OM = X_H2OM_w
......@@ -1061,6 +1073,19 @@ contains
pr%ec(k)