Commit 44403c63 authored by Lena Noack's avatar Lena Noack
Browse files

Update interior structure model if first iteration guess is completely off

parent ba5997ba
......@@ -99,6 +99,10 @@ 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
Mc = pE%X_Fe * M
Mw = pE%X_H2O * M
Mm = M - Mc - Mw
......@@ -179,6 +183,8 @@ contains
pr%gb(nr) = 6.67384e-11_dp * M / pr%Rp**2 ! surface gravity
pr%pb(nr) = pE%Psurf!*pF%D*pF%rho*pF%g
!write(*,*) pr%gb(nr),pr%pb(nr),M,pr%Rp
pr%xs = 0.0_dp ! light elements in core
pr%xs(1:pr%nc) = pE%XiS0 ! Light elements in core: in beginning equally distributed (later then via particles different)
pr%xw(:) = 0.0_dp ! Water in core/mantle, for ocean: set 0 since it is water anyway...
......@@ -261,6 +267,12 @@ 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
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))
pr%pb(k-1) = pr%pb(k) + pr%gb(k)*pr%rho(k)*(pr%rb(k)-pr%rb(k-1))
pr%Tb(k-1) = pr%Tb(k) * exp(pr%alpha(k)*pr%gb(k)/pr%Cp(k)*(pr%rb(k)-pr%rb(k-1))) ! T = T_t exp( alpha_t g / Cp_t * dr )
......@@ -629,7 +641,9 @@ contains
do k=1,nr
mass(k) = 4.0_dp*pi_*(pr%rb(k)**3-pr%rb(k-1)**3)*pr%rho(k)/3.0_dp + mass(k-1) ! density (pr%rho) is already average shell density
pr%gb(k) = 6.67384e-11_dp * mass(k) / pr%rb(k)**2
!write(*,*) k, mass(k), pr%gb(k), pr%rb(k), pr%rb(k-1), pr%rho(k), mass(k-1)
enddo
!write(*,*) iter,pr%gb(nr),mass(nr),pr%rb(nr),pr%rho(nr)
! pressure (at layer boundaries) for given radii and densities
do kinv = 1,nr
......@@ -1212,16 +1226,9 @@ contains
do kinv = 1,nr
k = nr+1-kinv !do k=nr,1
if (pr%alpha(k).eq.0.0_dp) then
! EOS had a problem, zero values
if (k.ne.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
pr%alpha(k)=pI%alpha
pr%Cp(k)=pI%Cp
pr%rho(k)=pI%rho
endif
pr%pb(k-1) = pr%pb(k) + pr%gb(k)*pr%rho(k)*(pr%rb(k)-pr%rb(k-1))
pr%Tb(k-1) = damp * pr%Tb(k) * exp(pr%alpha(k)*pr%gb(k)/pr%Cp(k)*(pr%rb(k)-pr%rb(k-1))) & ! T = T_t exp(alpha_t g / Cp_t * dr )
......
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