model1D.f90

module model1D
  use parameters, only:param_E,param_F,param_I,param_M,param_O,param_T,param_V
  use bouchet ! author: Attilio Rivoldini, 2015
  use thermoElastic ! author: Attilio Rivoldini, 2015
  use stixrudeMantleSpecies
  use speciesAssemblage
  use precision
  use water
  use int_struct
  implicit none

  type props1D
    real(dp),allocatable::rb(:),Tb(:),Tmb(:),pb(:),gb(:) ! at boundaries
    real(dp),allocatable::r(:),eta(:),rho(:),k(:),alpha(:),Cp(:),H(:),v(:),vl(:),vr(:),d(:),m(:),xs(:),xw(:) ! in center of each layer
    real(dp),allocatable::gr(:),KT(:),KS(:),GS(:),S(:),ec(:)
    real(dp),allocatable::Vb(:),Vc(:),Lc(:),dr(:)
    real(dp)::rho_c,rho_m,rho_w ! averaged densities, needed for update of interior structure
    real(dp)::alpha_c,alpha_m,alpha_w,Cp_c,Cp_m,Cp_w,k_m,k_w,g_m,shear_m,gr_m
    real(dp)::Rc,Rm,Rp,Dl
    real(dp)::DTc ! temperature jump (either fixed or determined via melting temperature)
    integer::nc,nm,nw
    real(dp)::Tc_av,Tm_av,Tw_av,rho_av ! output values, averaged value
  end type props1D

!
!  ----------- rb(nr) = Rp (planet surface radius), nr = nc+nm+nw
!       -      r(nr), rho(nr), v(nr), ...
!  ----------- rb(nr-1), Tb(nr-1)
!
!
!  rb(nc) = Rc, rb(nc+nm) = Rm
!
!
!  ----------- rb(1)
!       -      r(1) (core)
!  ----------- rb(0) = center planet, r=0
!

contains

  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!                                                                                                                           !!
  !!                                                     Initialization                                                        !!
  !!                                                                                                                           !!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  subroutine Init1DProf(pE,pF,pI,pM,pO,pT,pV,pr,debug,nr,tyr)
    type(param_E),intent(inout)::pE
    type(param_F),intent(in)::pF
    type(param_I),intent(inout)::pI 
    type(param_M),intent(in)::pM
    type(param_O),intent(in)::pO
    type(param_T),intent(in)::pT
    type(param_V),intent(in)::pV
    type(props1D),intent(inout)::pr
    integer,intent(in)::debug,nr
    real(dp),intent(in)::tyr

    real(dp)::R_E,Mc,Mw,Mm,pi_,M,t,Tcore,Tcore_old,p_GPa,T_K,vol,vol_c,vol_m,vol_w,XmFe
    real(dp)::t_T,t_p,t_m
    integer::i,k,kinv,lTBL_k,error,maxOI,iter,nc_old
    real(dp)::res(9),res_c(6)

    if (pE%IntStruct.ge.0) &
    &write(*,*) "nr=",nr

    allocate(pr%rb(0:nr),pr%Tb(0:nr),pr%Tmb(0:nr),pr%gb(0:nr),pr%pb(0:nr))
    allocate(pr%r(nr),pr%eta(nr),pr%rho(nr),pr%k(nr),pr%alpha(nr),pr%Cp(nr),pr%xs(nr),pr%xw(nr))
    allocate(pr%H(nr),pr%v(nr),pr%vl(nr),pr%vr(0:nr),pr%d(nr),pr%m(nr))
    allocate(pr%gr(nr),pr%KT(nr),pr%KS(nr),pr%GS(nr),pr%S(nr),pr%ec(nr))
    allocate(pr%Vb(nr),pr%Vc(nr),pr%Lc(nr),pr%dr(nr))

    pr%Tmb=0.0_dp;pr%k=0.0_dp;pr%xs=0.0_dp;pr%xw=0.0_dp;pr%H=0.0_dp;pr%d=0.0_dp
    pr%v=0.0_dp;pr%vl=0.0_dp;pr%vr=0.0_dp;pr%eta=0.0_dp


    ! test profiles for p/T/KT/G along melting line for partition coefficients
    !t_p = 0.0_dp
    !t_T = 0.0_dp
    !t_m = 3.0_dp
    !do while (t_p.le.15_dp)
    !  call MeltPhase(100.0_dp,t_T,t_p*1e+9_dp,t_m,0.0_dp,0.1_dp)
    !  res = assemblageThermoElastic(t_p,t_T,[forsterite,fayalite],[0.9_dp,1.0_dp],error)
    !
    !  write(*,*) t_p, t_T, res(6), res(8)
    !
    !  t_p = t_p + 0.5_dp
    !enddo


    M = pE%M_E * 5.9736e+24 ! Mass in g
    pi_ = acos(-1.0_dp)
    maxOI = 20 !100

    ! note: Get1DStruct() uses dimensional(!) values, calculation Run1DProf() uses either dimensional or non-dimensional values

    ! simplified scaling law for first iteration, will be updapted by Get1DStruct()
    pr%rho_c = 11032.0_dp * pE%M_E**0.2612_dp
    pr%rho_m = 4457.0_dp * pE%M_E**0.1779_dp
    pr%rho_w = 2120.8_dp * pE%M_E**0.2007_dp

    Mc = pE%X_Fe * M
    Mw = pE%X_H2O * M
    Mm = M - Mc - Mw

    ! subtract iron in mantle
    !Mc = Mc - pI%X_FeM * Mm
    !if (Mc.lt.0.0_dp) then
    !  pI%X_FeM = pE%X_Fe ! all iron in mantle
    !  Mc = 0.0_dp
    !  if (debug.ge.1) write(*,*) "Error: Mc is negative, set Rc to 0 and Fe_Mantle to ",pI%X_FeM
    !endif
    !Mm = M - Mc - Mw
    
    !if (pI%X_FeM.gt.pE%X_Fe) pI%X_FeM=pE%X_Fe
    !Mm = Mm / (1.0_dp-pI%X_FeM)  ! hence pI%X_FeM is iron content of Mm
    !Mc = M - Mm - Mw
    !if (pI%X_FeM.eq.pE%X_Fe) Mc=0.0_dp

    ! X_FeM is iron number in mantle, assuming mantle composed of (1-X_FeM) Mg2SiO4 + X_FeM Fe2SiO4, how much iron needs to be extracted from core?
    ! molar masses (since Mc weight fraction of planet mass): Mg=24.305, Fe=55.845, Si=28.0855, O=15.9994
    !Mc = Mc - 2.0_dp*pI%X_FeM*55.845_dp/(pI%X_FeM*63.08_dp+140.69_dp)*Mm
    XmFe = 2.0_dp*pI%X_FeM*55.845_dp/(pI%X_FeM*63.08_dp+140.69_dp)
    Mc = ( pE%X_Fe*M - XmFe*(M-Mw) ) / ( 1.0_dp - XmFe)
    if (Mc.lt.0.0_dp) then
      write(*,*) "Error, not enough iron in planet to have ",pI%X_FeM," as iron number!"
      Mc=0.0_dp
    endif
    !Mm = Mm + 2.0_dp*pI%X_FeM*55.845_dp/(pI%X_FeM*63.08_dp+140.69_dp)*Mm
    Mm = M - Mc - Mw

    ! simplified scaling law for first radii, will be updated by Get1DStruct()

    pr%Rc = ( Mc*3.0_dp / (4.0_dp*pi_*pr%rho_c) )**(1.0_dp/3.0_dp)
    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)

    ! 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
    pr%nm = 0
    pr%nw = 0
    do k=1,nr
      pr%rb(k) = real(k)*pr%Rp/real(nr)
      pr%r(k) = 0.5_dp * ( pr%rb(k-1)+pr%rb(k) )
      if (pr%rb(k).le.pr%Rc) then
        pr%rho(k) = pr%rho_c
        pr%m(k) = 6 ! material: 1-water/ice, 2-crust, 3-silicates:olivine, 4-silicates:pv+mw, 5-silicates:ppv+mv, 6-FeS or pure Fe
        pr%nc = k ! update until not in core anymore -> correct number of grid points in core
      else if (pr%rb(k).le.pr%Rm) then
        pr%rho(k) = pr%rho_m
        pr%m(k) = 3
        pr%nm = k ! update until not in mantle anymore -> correct number of grid points in mantle
      else
        pr%rho(k) = pr%rho_w
        pr%m(k) = 1 
        pr%nw = k ! update until not in ocean anymore -> correct number of grid points in ocean
      endif
    enddo

    pr%gb(nr) = 6.67384e-11_dp * M / pr%Rp**2 ! surface gravity
    pr%pb(nr) = pE%Psurf!*pF%D*pF%rho*pF%g

    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...

    pr%DTc = 0.0_dp

    if (debug.gt.1) write(*,*) "Set preliminary temperature profile"

    ! set preliminary temperature profile 
    pr%Tb(nr) = pT%Ttop ! pT%Ttop*pF%DeltaT+pF%Ts
    do kinv = 1,nr
      k = nr+1-kinv !do k=nr,1

      if (debug.ge.1) write(*,*) "Ini T: k=",k,"(",nr,"),T=",pr%Tb(k)

      p_GPa = pr%pb(k)*1.0e-9_dp
      T_K = pr%Tb(k)

      !call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM)

      if (pr%rb(k).gt.pr%Rm) then
        pr%m(k) = 1 ! first assume water
        call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)

        if (pr%m(k).eq.1) then
          ! liquid water
          call getProp(pr%pb(k),min(pr%Tb(k),645.0_dp),pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,1) ! 1=water,2=ice VI, 3=VII, 4=X, 5=ol, 6=pv/mw, 7=ppv/mw, 8=Fe
          !if (debug.gt.2) write(*,*) k,pr%Tb(k),pr%rho(k) 
        else
          ! ice (here ice VI, ToDo: include Ih, III and V)
          if (pr%pb(k).lt.2.216e+9_dp) then !ice VI
            call getProp(pr%pb(k),pr%Tb(k),pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,2)
          else if (pr%pb(k).lt.47.0e+9_dp) then !ice VII after Schlager et al., 2004
            call getProp(pr%pb(k),pr%Tb(k),pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,3)
          else ! ice X
            call getProp(pr%pb(k),pr%Tb(k),pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,4)
          endif
        endif
        !call getProp(pr%pb(k),pr%Tb(k),pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,1)
      else if (pr%rb(k).gt.pr%Rc) then
        pr%m(k) = 3 ! first assume olivine
        call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
        if (pr%pb(k).le.(4.92_dp+5.4_dp*(1.0_dp-pI%X_FeM)+0.0025_dp*pr%Tb(k))*1.0e+9_dp) then
          res = assemblageThermoElastic(p_GPa,T_K,[forsterite,fayalite],[1.0_dp-pI%X_FeM,pI%X_FeM],error)
          pr%rho(k) = res(2)
          pr%alpha(k) = res(3)
          pr%Cp(k) = res(4)
        elseif (pr%pb(k).le.(6.2_dp+16.7_dp*(1.0_dp-pI%X_FeM)+0.005_dp*(pr%Tb(k)-2260.0_dp))*1.0e+9_dp) then
          res = assemblageThermoElastic(p_GPa,T_K,[MgWadsleyite,FeWadsleyite],[1.0_dp-pI%X_FeM,pI%X_FeM],error)
          pr%rho(k) = res(2)
          pr%alpha(k) = res(3)
          pr%Cp(k) = res(4)
        elseif (pr%pb(k).le.(22.9_dp-0.0025_dp*(pr%Tb(k)-2260.0_dp))*1.0e+9_dp) then 
          res = assemblageThermoElastic(p_GPa,T_K,[MgRingwoodite,FeRingwoodite],[1.0_dp-pI%X_FeM,pI%X_FeM],error)
          pr%rho(k) = res(2)
          pr%alpha(k) = res(3)
          pr%Cp(k) = res(4)
        else if (pr%pb(k).le.(109.54_dp+0.0058_dp*pr%Tb(k))*1.0e+9_dp) then 
          res = assemblageThermoElastic(p_GPa,T_K,[MgPerovskite,periclase,FePerovskite,wustite],&
                &[0.7135_dp*(1.0_dp-pI%X_FeM),0.2865_dp*(1.0_dp-pI%X_FeM),0.7135_dp*pI%X_FeM,0.2865_dp*pI%X_FeM],error)
          pr%rho(k) = res(2)
          pr%alpha(k) = res(3)
          pr%Cp(k) = res(4)
        else
          res = assemblageThermoElastic(p_GPa,T_K,[MgPostPerovskite,periclase,FePostPerovskite,wustite],&
                &[0.7135_dp*(1.0_dp-pI%X_FeM),0.2865_dp*(1.0_dp-pI%X_FeM),0.7135_dp*pI%X_FeM,0.2865_dp*pI%X_FeM],error)
          pr%rho(k) = res(2)
          pr%alpha(k) = res(3)
          pr%Cp(k) = res(4)
        endif
!        call getProp(pr%pb(k),pr%Tb(k),pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,5)
      else
        res_c = eos_core(p_GPa,T_K,error)
        pr%rho(k) = res_c(2)
        pr%alpha(k) = res_c(3)
        pr%Cp(k) = res_c(4)
!        call getProp(pr%pb(k),pr%Tb(k),pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,8)
      endif
      !if (pr%alpha(k).eq.0.0_dp) pr%alpha(k)=pr%alpha(k+1) ! can happen for initial profile, does not happen anymore later when profiles converge
      !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)

      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 )

      ! lithosphere plus upper TBL
      pr%Dl = pE%Dl + pE%Delta_u ! later: depending on viscosity profile and thus velocity profile 

      if ((pr%rb(k).le.pr%Rm).and.(pr%rb(k).gt.pr%Rm-pr%Dl)) pr%Tb(k-1) = pT%Ttop + (pr%Rm-pr%rb(k))/pr%Dl*(abs(pE%Tm)-pT%Ttop)
      ! at bottom of lithosphere
      if ((pr%rb(k).gt.pr%Rm-pr%Dl).and.(pr%rb(k-1).le.pr%Rm-pr%Dl)) then
        if (pE%Tm.lt.0.0_dp) then
          pE%Tm = -pr%Tmb(k) ! set temp at bottom of lid to melting temp at that depth (lid temp is adjusted in iteration below)
          pr%Tb(k-1) = abs(pE%Tm)
        else
          pr%Tb(k-1) = max(pE%Tm,pr%Tb(k-1))
        endif
        !write(*,*) "set upper mantle with Tm=",pE%Tm,pr%Tb(k),pr%Tmb(k)
      endif
      ! find radius where lower TBL starts
      if ((pr%rb(k).gt.pr%Rc+pE%Delta_c).and.(pr%rb(k-1).le.pr%Rc+pE%Delta_c)) then
        ! store k, temp will be set once CMB temp was calculated
        lTBL_k = k
      endif
      ! at CMB
      if ((pr%rb(k-1).le.pr%Rc).and.(pr%rb(k).gt.pr%Rc)) then
        ! add temperature jump at CMB
        if (pE%DTc.gt.0) then
          pr%Tb(k-1)=pr%Tb(k-1) + pE%DTc
        elseif (pE%DTc.lt.0) then
          ! scale CMB temp difference with -pE%DTc (negtive value: scale difference to silicte melting temperature at CMB with -pE%DTc)
          if (pr%Tmb(k)-pr%Tb(k).gt.0.0_dp) then
            pr%DTc = abs(pE%DTc)*(pr%Tmb(k)-pr%Tb(k))
            !write(*,*) "Temperature jump at CMB=",abs(pE%DTc)*(pr%Tmb(k)-pr%Tb(k)),"K",", Tm=",pr%Tmb(k),", Tad=",pr%Tb(k)
            pr%Tb(k-1)=pr%Tb(k) + abs(pE%DTc)*(pr%Tmb(k)-pr%Tb(k)) ! if DTc=-1 then temperature at CMB is set to melting temperature of peridotite
          endif
        endif
        ! 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
          enddo
        endif
      endif
      if (pr%Tb(k-1).gt.2.0_dp*pr%Tb(k)) 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)
        pr%alpha(k) = pr%alpha(k+1)
        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)

    enddo

    if (pr%nc.eq.0) pr%Tb(0)=pr%Tb(1) ! if no core

    if (debug.ge.1) write(*,*) "Tcore=",pr%Tb(pr%nc)

    if (pr%nm.ne.0) pr%nm = pr%nm - pr%nc ! only number of points in mantle
    if (pr%nw.ne.0) pr%nw = pr%nw - pr%nm - pr%nc ! only number of points in water-ice layer

    ! check number of gridpoints in each layer:
    if (debug.ge.1) write(*,*) "nr=",nr,", nc=",pr%nc,", nm=",pr%nm,", nw=",pr%nw

    Tcore = 1000.0_dp
    Tcore_old = 0.0_dp
    iter = 0
    if (debug.ge.2) write(*,*) "-----------------------------"
    do while ((abs(Tcore-Tcore_old).gt.0.0001_dp).and.(iter.le.maxOI))
      iter=iter+1
      nc_old=pr%nc
      ! get interior structure for T-profile depending on pE%H2O and pE%Fe
      if (debug.ge.1) write(*,*) "Get1DStruct for T-profile"
      if (Tcore_old.eq.0.0_dp) then !first call -> T-prof can be completely off, just few iterations
        call Get1DStruct(pr,pE%M_E,pE%X_Fe,pE%X_H2O,pr%rho_c,pr%rho_m,pr%rho_w,pE%Psurf,pT%Ttop,R_E,nr,debug,pr%xs,0,&
                        &pI%X_FeM,pI%X_H2OM,tyr,pE,pI,2,iter)
      else
        call Get1DStruct(pr,pE%M_E,pE%X_Fe,pE%X_H2O,pr%rho_c,pr%rho_m,pr%rho_w,pE%Psurf,pT%Ttop,R_E,nr,debug,pr%xs,0,&
                        &pI%X_FeM,pI%X_H2OM,tyr,pE,pI,maxOI,iter)
      endif

      ! when oscillations occur, nc can jump between two indices
!      nc_old = pr%nc+1 !NINT(real(nc_old+pr%nc)/2.0_dp) ! nearest integer

      ! Update temp now directly in Get1DStruct

!      !Update Temp Profile with new material properties
!      if (debug.ge.1) write(*,*) "Update Temp"
!      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
!        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) = 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 )
!        if (INAN(pr%Tb(k-1))) write(*,*) "Error NaN Tb,",pr%Tb(k),pr%alpha(k),pr%gb(k),pr%Cp(k),k
!        if ((pr%rb(k).le.pr%Rm).and.(pr%rb(k).gt.pr%Rm-pE%Dl)) pr%Tb(k-1) = pT%Ttop + (pr%Rm-pr%rb(k))/pE%Dl*(pE%Tm-pT%Ttop)
!        if ((pr%rb(k).gt.pr%Rm-pE%Dl).and.(pr%rb(k-1).le.pr%Rm-pE%Dl)) then
!          pr%Tb(k-1) = max(pE%Tm,pr%Tb(k-1))
!        endif
!        if ((pr%rb(k-1).le.pr%Rc).and.(pr%rb(k).gt.pr%Rc)) then
!          ! add temperature jump at CMB
!          if (pE%DTc.gt.0) then
!            pr%Tb(k-1)=pr%Tb(k-1) + pE%DTc
!          elseif (pE%DTc.lt.0) then
!            ! scale CMB temp difference with -pE%DTc (negative value: scale difference to silicte melting temperature at CMB with -pE%DTc)
!            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
!              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
!              write(*,*) "Error, mantle frozen at CMB, no temp jump applied",k-1,k,nc_old,pr%nc
!            endif
!          endif
!        endif
!
!      enddo
      Tcore_old = Tcore
      Tcore = pr%Tb(pr%nc)
      k = pr%nc+1
      !write(*,*) "Tcore=",Tcore,pr%Tb(k),exp(pr%alpha(k)*pr%gb(k)/pr%Cp(k)*(pr%rb(k)-pr%rb(k-1))),&
      !        &pr%alpha(k),pr%gb(k),pr%Cp(k),(pr%rb(k)-pr%rb(k-1)),pr%rb(k)
    enddo
    call Get1DStruct(pr,pE%M_E,pE%X_Fe,pE%X_H2O,pr%rho_c,pr%rho_m,pr%rho_w,pE%Psurf,pT%Ttop,R_E,nr,debug,pr%xs,1,&
                    &pI%X_FeM,pI%X_H2OM,tyr,pE,pI,maxOI,iter)


    !!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Heat sources in mantle !
    !!!!!!!!!!!!!!!!!!!!!!!!!!

    pr%H = 0.0_dp ! values in core and water-ice layer stays zero for now, mantle values changed in the following:

    t=4500000000.0_dp ! for re-scaling of present-day heat sources

    ! specific heating rates in W/kg
    pI%H0U238  = 9.46e-5_dp
    pI%H0U235  = 5.69e-4_dp
    pI%H0Th232 = 2.64e-5_dp
    pI%H0K40   = 2.92e-5_dp 
    ! half-life times in yr
    pI%t12U238  = 4.47e+9_dp
    pI%t12U235  = 7.04e+8_dp
    pI%t12Th232 = 1.40e+10_dp
    pI%t12K40   = 1.25e+9_dp

    ! get initial global concentration from enriched crustal value
    if (pE%H_from_surf.ge.1) then
      pI%C_U  = pI%C_U  / pM%Lambda
      pI%C_Th = pI%C_Th / pM%Lambda
      pI%C_K  = pI%C_K  / pM%Lambda
    endif
 
    if (pE%IntStruct.ge.0) &
    &write(*,*) "Bulk concentrations at 4.5 Gyr: U=",pI%C_U," ppb, Th=",pI%C_Th," ppb, K=",pI%C_K," ppm"

    pI%X0U238 = 0.9928_dp * pI%C_U * 1.0e-9 ! C in ppb
    pI%X0U235 = 0.0071_dp * pI%C_U * 1.0e-9 ! C in ppb
    pI%X0Th232 =            pI%C_Th * 1.0e-9 !C in ppb
    pI%X0K40 = 0.000128_dp* pI%C_K * 1.0e-6 ! C in ppm

    if (pE%IntStruct.ge.0) &
    &write(*,*) "Concentrations after 4.5 Gyr: U235=",pI%X0U235,",U238=",pI%X0U238,",Th:",pI%X0Th232,",K:",pI%X0K40

    if (pE%IntStruct.ge.0) &
    &write(*,*) "Time factors",exp( Log(2.0_dp)*t/pI%t12U238 ),exp( Log(2.0_dp)*t/pI%t12U235 ),&
             &exp( Log(2.0_dp)*t/pI%t12Th232 ),exp( Log(2.0_dp)*t/pI%t12K40 )
    if (pE%IntStruct.ge.0) &
    &write(*,*) "Exponents",( Log(2.0_dp)*t/pI%t12U238 ),( Log(2.0_dp)*t/pI%t12U235 ),&
             &( Log(2.0_dp)*t/pI%t12Th232 ),( Log(2.0_dp)*t/pI%t12K40 )

    ! set back concentration to values at t=0.0Gyr (values in input file are present-day values)
    pI%X0U238 = pI%X0U238 * exp( Log(2.0_dp)*t/pI%t12U238 )
    pI%X0U235 = pI%X0U235 * exp( Log(2.0_dp)*t/pI%t12U235 )
    pI%X0Th232 = pI%X0Th232 * exp( Log(2.0_dp)*t/pI%t12Th232 )
    pI%X0K40 = pI%X0K40 * exp( Log(2.0_dp)*t/pI%t12K40 )

    if (pE%IntStruct.ge.0) &
    &write(*,*) "Initial concentrations: U235=",pI%X0U235,",U238=",pI%X0U238,",Th:",pI%X0Th232,",K:",pI%X0K40

    ! H so far constant in mantle and crust, ToDo: adapt to different layers
    pr%H(pr%nc+1:pr%nc+pr%nm) = pI%X0U238*pI%H0U238 + pI%X0U235*pI%H0U235 + pI%X0Th232*pI%H0Th232 + pI%X0K40*pI%H0K40


    ! preliminary: set core temp via initial bottom temperature
!    pr%Tb(pr%nc) = pT%Tbottom
!    do kinv=1,pr%nc
!      k = pr%nc+1-kinv  ! from nc to 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 )
!    enddo
!
!    ! preliminary: set ocean temp to Ttop
!    do k=pr%nc+pr%nm+1,nr
!      pr%Tb(k) = pT%Ttop
!    enddo

    ! get average values for output file
    pr%Tc_av = 0.0_dp
    pr%Tm_av = 0.0_dp
    pr%Tw_av = 0.0_dp
    pr%rho_av = 0.0_dp
    vol = 0.0_dp
    vol_c = 0.0_dp
    vol_m = 0.0_dp
    vol_w = 0.0_dp

    do k=1,nr
      if (k.le.pr%nc) then
        pr%Tc_av = pr%Tc_av + 0.5_dp*(pr%Tb(k)+pr%Tb(k-1))*(4.0_dp/3.0_dp*pi_)*(pr%rb(k)**3-pr%rb(k-1)**3)
        vol_c = vol_c + (4.0_dp/3.0_dp*pi_)*(pr%rb(k)**3-pr%rb(k-1)**3)
      elseif (k.le.pr%nc+pr%nm) then
        pr%Tm_av = pr%Tm_av + 0.5_dp*(pr%Tb(k)+pr%Tb(k-1))*(4.0_dp/3.0_dp*pi_)*(pr%rb(k)**3-pr%rb(k-1)**3)
        vol_m = vol_m + (4.0_dp/3.0_dp*pi_)*(pr%rb(k)**3-pr%rb(k-1)**3)
      else
        pr%Tw_av = pr%Tw_av + 0.5_dp*(pr%Tb(k)+pr%Tb(k-1))*(4.0_dp/3.0_dp*pi_)*(pr%rb(k)**3-pr%rb(k-1)**3)
        vol_w = vol_w + (4.0_dp/3.0_dp*pi_)*(pr%rb(k)**3-pr%rb(k-1)**3)
      endif
      pr%rho_av = pr%rho_av + pr%rho(k)*(4.0_dp/3.0_dp*pi_)*(pr%rb(k)**3-pr%rb(k-1)**3)
      vol = vol + (4.0_dp/3.0_dp*pi_)*(pr%rb(k)**3-pr%rb(k-1)**3)
    enddo

    if (vol_c.eq.0.0_dp) vol_c = 1.0_dp
    if (vol_m.eq.0.0_dp) vol_m = 1.0_dp
    if (vol_w.eq.0.0_dp) vol_w = 1.0_dp

    pr%Tc_av = pr%Tc_av/vol_c
    pr%Tm_av = pr%Tm_av/vol_m
    pr%Tw_av = pr%Tw_av/vol_w
    pr%rho_av = pr%rho_av/vol

    !write(*,*) "Pressure at center:",pr%pb(0)*1.0e-9_dp," GPa"

  end subroutine Init1DProf


  ! is called both from Init and running simulation to obtain changes in radii due to thermal contraction
  subroutine Get1DStruct(pr,M_E,X_Fe,X_H2O,rho_c,rho_m,rho_w,Psurf,Ttop,R_E,nr,debug,xs,&
                        &writeIS,X_FeM,X_H2OM_w,tyr,pE,pI,maxOI,iterO)
    type(props1D),intent(inout)::pr
    real(dp),intent(in)::M_E,X_Fe,X_H2O,rho_c,rho_m,rho_w,Psurf,Ttop,X_FeM,X_H2OM_w,tyr
    real(dp),intent(inout)::xs(:)
    type(param_E),intent(inout)::pE
    type(param_I),intent(inout)::pI
    real(dp),intent(out)::R_E !radius in Earth radii
    integer,intent(in)::nr,debug,writeIS,maxOI,iterO

    integer::iter,i,k,kinv,lTBL_k,error,nc_old
    real(dp)::tolOI,M,Mc,Mm,Mw,Rp_old,Tc_old,dr,Vc,Vm,Vw,vol,pi_,p,p_GPa,T,kB,X_H2OM,XmFe,damp
    real(dp)::Rp_old2,Rp_old3,Rp_old4
    real(dp)::HtU235,HtU238,HtTh232,HtK40
    real(dp)::Cp_c,Cp_m,Cp_w,alpha_c,alpha_m,alpha_w,k_m,k_w,shear_m,gr_m,ec_log,gr_c,R_ICB,p_ICB
    real(dp)::X_FeM_mol ! molar fraction of iron in mantle
    real(dp),allocatable::mass(:)
    real(dp)::res(9),res_c(6)
    character(4)::str_M
    character(7)::str_Fe,str_H2O

!    type(phase)::forsterite,MgPerovskite,MgPostPerovskite,periclase

!    maxOI = 100 ! in input.txt?
    tolOI = 1.0e-8_dp
    error = 0

    M = M_E * 5.9736e+24 ! Mass in g
    pi_ = acos(-1.0_dp)  
    kB = 8.61733e-5_dp !eV/K - Boltzman constant
    !X_H2OM = 0.0_dp !mantle water content in wt-%, for now use constant value here

    Mc = X_Fe * M
    Mw = X_H2O * M
    Mm = M - Mc - Mw

    !write(*,*) X_Fe,X_H2O,X_FeM,Mc,Mw,Mm,Mm/(1.0_dp-X_FeM),M-Mw-Mm/(1.0_dp-X_FeM)

    ! subtract iron in mantle
!    Mc = Mc - X_FeM * Mm
!    Mm = M - Mc - Mw
!    Mm = Mm / (1.0_dp-X_FeM) 
!    Mc = M - Mm - Mw
!    if (X_FeM.eq.X_Fe) Mc=0.0_dp
!    if (Mc.lt.0.0_dp) Mc=0.0_dp ! rounding errors can lead to negative core radius
!
    X_H2OM = X_H2OM_w*100.0_dp ! from weight value to weight percent as needed here
    !X_FeM_mol = X_FeM * (X_FeM*203.7731_dp+(1.0_dp-X_FeM)*140.6931_dp)/203.7731_dp !conversion mass fraction to molar fraction for Fe2SiO4-Mg2SiO4 solid solution

    ! X_FeM is iron number in mantle, assuming mantle composed of (1-X_FeM) Mg2SiO4 + X_FeM Fe2SiO4, how much iron needs to be extracted from core?
    ! molar masses (since Mc weight fraction of planet mass): Mg=24.305, Fe=55.845, Si=28.0855, O=15.9994
    !Mc = Mc - 2.0_dp*X_FeM*55.845_dp/(X_FeM*63.08_dp+140.69_dp)*Mm
    XmFe = 2.0_dp*X_FeM*55.845_dp/(X_FeM*63.08_dp+140.69_dp) ! iron weight fraction in mantle
    X_FeM_mol = X_FeM ! input X_FeM is already molar fraction (iron number)
    Mc = ( X_Fe*M - XmFe*(M-Mw) ) / ( 1.0_dp - XmFe)
    if (Mc.lt.0.0_dp) then
      write(*,*) "Error, not enough iron in planet to have ",pI%X_FeM," as iron number!"
      Mc=0.0_dp
    endif
    Mm = M - Mc - Mw
    !Mm = Mm + 2.0_dp*X_FeM*55.845_dp/(X_FeM*63.08_dp+140.69_dp)*Mm
    

!    nc = nint(nr*Rc/Rp) ! get integer value for number of grid points in core
!    nm = nint(nr*(Rm-Rc)/Rp) ! get integer value for number of grid points in mantle
!    nw = nint(nr*(Rp-Rm)/Rp) ! get integer value for number of grid points in water-ice layer

    allocate( mass(0:nr) )

    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! Iteration
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    iter = 0 ! iteration: get new radii (changes with more accurate pressure/density profiles)
    Rp_old = 0.0_dp
    Rp_old2 = 0.0_dp
    Rp_old3 = 0.0_dp
    Rp_old4 = 0.0_dp
    Tc_old = 0.0_dp

!    do while ((iter.le.maxOI).and.((abs(Tc_old-pr%Tb(pr%nc))/max(Tc_old,pr%Tb(pr%nc))).gt.tolOI)&
!            &.and.((abs(Rp_old-pr%Rp)/max(Rp_old,pr%Rp)).gt.tolOI))
    do while ((iter.le.maxOI).and.((abs(Rp_old-pr%Rp)/max(Rp_old,pr%Rp)).gt.tolOI))
      Rp_old4 = Rp_old3
      Rp_old3 = Rp_old2
      Rp_old2 = Rp_old
      Rp_old = pr%Rp  
      Tc_old = pr%Tb(pr%nc)
      iter = iter + 1

      ! mass and gravity for given radii and densities
      mass(0) = 0.0_dp ! also at shell boundaries, i.e. nr is planet surface
      pr%gb(0) = 0.0_dp
      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
      enddo

      ! pressure (at layer boundaries) for given radii and densities
      do kinv = 1,nr
        k = nr-kinv !do k=nr-1,0

        dr = pr%rb(k+1)-pr%rb(k)
        pr%pb(k) = 0.5_dp*(pr%gb(k+1)+pr%gb(k)) * pr%rho(k+1) * dr + pr%pb(k+1) ! rho(k+1) is density in cell center of cell k between rb(k) and rb(k+1)
      enddo

      !!!!!!!!!!!!!!!!!!
      ! surface values !
      !!!!!!!!!!!!!!!!!!
      ! cell center values
      p = 0.5_dp*(pr%pb(nr)+pr%pb(nr-1))
      p_GPa = p*1.0e-9_dp
      T = 0.5_dp*(pr%Tb(nr)+pr%Tb(nr-1))

      if (X_H2O.gt.0.0_dp) then
        pr%m(nr) = 1 ! first assume water
        call MeltPhase(pr%Tb(nr),pr%Tmb(nr),pr%pb(nr),pr%m(nr),pr%xs(nr),pI%X_FeM,pI%X_M0)

        if (pr%m(nr).eq.1) then
          ! liquid water
          if (debug.gt.2) write(*,*) nr,p,T
          call getProp(p,min(T,645.0_dp),pr%rho(nr),pr%alpha(nr),pr%Cp(nr),pr%k(nr),debug,error,1) ! 1=water, 2=ice VI, 3=VII, 4=X, 5=ol, 6=pv/mw, 7=ppv/mw, 8=Fe
        else
          ! ice (here ice VI, ToDo: include Ih, III and V)
          if (p_GPa.lt.2.216_dp) then !ice VI
            call getProp(p,T,pr%rho(nr),pr%alpha(nr),pr%Cp(nr),pr%k(nr),debug,error,2)
          else if (p_GPa.lt.47.0_dp) then !ice VII after Schlager et al., 2004
            call getProp(p,T,pr%rho(nr),pr%alpha(nr),pr%Cp(nr),pr%k(nr),debug,error,3)
          else ! ice X
            call getProp(p,T,pr%rho(nr),pr%alpha(nr),pr%Cp(nr),pr%k(nr),debug,error,4)
          endif
        endif
        pr%gr(nr) = 0.0_dp
        pr%KT(nr) = 0.0_dp
        pr%KS(nr) = 0.0_dp
        pr%GS(nr) = 0.0_dp
        pr%S(nr) = 0.0_dp
        pr%ec(nr) = 0.0_dp
      else if (X_Fe.lt.1.0_dp) then
         if (pE%Dcr0.gt.0.0_dp) then
           ! simple mafic crust as mixture between diopside and anorthite, if more felsic then more quarz (not in St11) and albite instead
          pr%m(nr) = 2
          call MeltPhase(pr%Tb(nr),pr%Tmb(nr),pr%pb(nr),pr%m(nr),pr%xs(nr),pI%X_FeM,pI%X_M0) ! melt temp as for upper mantle for now
          pr%k(nr) = pI%k_cr
          res = assemblageThermoElastic(p_GPa,T,[diopside,anorthite],[0.5_dp,0.5_dp],error)
         else
          pr%m(nr) = 3 ! upper mantle olivine-rich
          call MeltPhase(pr%Tb(nr),pr%Tmb(nr),pr%pb(nr),pr%m(nr),pr%xs(nr),pI%X_FeM,pI%X_M0)
          !call getProp(p,T,pr%rho(nr),pr%alpha(nr),pr%Cp(nr),pr%k(nr),debug,error,5)  ! keep this call for thermal conductivity, only
          !res = eos(p_GPa,T,forsterite,error) ! ToDo: try actual crust composition??? Also below...
          pr%k(nr) = (2.47_dp + 0.33_dp*p_GPa)*(300.0_dp/T)**0.48_dp ! mantle conductivities all follow Tosi et al., 2013, PEPI
          res = assemblageThermoElastic(p_GPa,T,[forsterite,fayalite],[1.0_dp-X_FeM,X_FeM],error)
        endif
        if (error.eq.0) then
          pr%rho(nr) = res(2)
          pr%alpha(nr) = res(3)
          pr%Cp(nr) = res(4)
          pr%gr(nr) = res(5)
          pr%KT(nr) = res(6)
          pr%KS(nr) = res(7)
          pr%GS(nr) = res(8) 
          pr%S(nr) = res(9)
        else
          pr%rho(nr) = 0.0_dp
          pr%alpha(nr) = 0.0_dp
          pr%Cp(nr) = 0.0_dp
          pr%gr(nr) = 0.0_dp
          pr%KT(nr) = 0.0_dp
          pr%KS(nr) = 0.0_dp
          pr%GS(nr) = 0.0_dp
          pr%S(nr) = 0.0_dp
        endif
        if (pI%ec_model.eq.1) then
          pr%ec(nr) = 10.0_dp**2.69_dp*exp(-1.62_dp/(kB*T))  ! Xu dry profile
        elseif (pI%ec_model.eq.2) then
          pr%ec(nr) = 10.0_dp**4.73_dp*exp(-2.31_dp/(kB*T)) &
                  & + 10.0_dp**1.9_dp*X_H2OM*exp(-(0.92_dp-0.16_dp*X_H2OM**(1.0_dp/3.0_dp))/(kB*T))
        else
          pr%ec(nr) = 10.0_dp**4.73_dp*exp(-2.31_dp/(kB*T)) &
                  & + 10.0_dp**2.98_dp*X_FeM_mol*exp(-(1.71_dp-2.0_dp*X_FeM_mol**(1.0_dp/3.0_dp))/(kB*T)) &
                  & + 10.0_dp**1.9_dp*X_H2OM*exp(-(0.92_dp-0.16_dp*X_H2OM**(1.0_dp/3.0_dp))/(kB*T))
        endif
      else
        pr%m(nr) = 8 ! Fe or FeS
        call MeltPhase(pr%Tb(nr),pr%Tmb(nr),pr%pb(nr),pr%m(nr),pr%xs(nr),pI%X_FeM,pI%X_M0)
        !call getProp(p_GPa,T,pr%rho(nr),pr%alpha(nr),pr%Cp(nr),pr%k(nr),debug,error,8) ! keep this call for thermal conductivity, only
        pr%k(nr) = 50.0_dp ! not used
        res_c = eos_core(p_GPa,T,error)
        if (error.eq.0) then
          pr%rho(nr) = res_c(2)
          pr%alpha(nr) = res_c(3)
          pr%Cp(nr) = res_c(4)
          pr%gr(nr) = res_c(5)
          pr%KT(nr) = res_c(6)
          pr%KS(nr) = 0.0_dp
          pr%GS(nr) = 0.0_dp
          pr%S(nr) = 0.0_dp
        else
          pr%rho(nr) = 0.0_dp
          pr%alpha(nr) = 0.0_dp
          pr%Cp(nr) = 0.0_dp
          pr%gr(nr) = 0.0_dp
          pr%KT(nr) = 0.0_dp
          pr%KS(nr) = 0.0_dp
          pr%GS(nr) = 0.0_dp
          pr%S(nr) = 0.0_dp
        endif
        pr%ec(nr) = 0.0_dp
      endif


      if (debug.ge.2) then
      write(*,*) "         nr           rho                       alpha                     Cp                       k&
            &                       p_GPa                       T                          r"
      write(*,*) nr,pr%rho(nr),pr%alpha(nr),pr%Cp(nr),pr%k(nr),p_GPa,T,pr%r(nr)
      endif

      ! From 2nd shell downwards
      do kinv = 1,nr-1
        k = nr-kinv !do k=nr-1,1

        dr = pr%rb(k+1)-pr%rb(k)
        pr%pb(k) = 0.5_dp*(pr%gb(k+1)+pr%gb(k)) * pr%rho(k+1) * dr + pr%pb(k+1) 

        ! cell center values
        p = 0.5_dp*(pr%pb(k)+pr%pb(k-1))
        p_GPa = p*1.0e-9_dp
        T = 0.5_dp*(pr%Tb(k)+pr%Tb(k-1))

        if (pr%rb(k).gt.pr%Rm) then ! ice/ocean
          pr%m(k) = 1 ! first assume water
          call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)

          if (debug.ge.1) write(*,*) "Test1:",k,p,T,pr%rb(k),pr%m(k)

          if (pr%m(k).eq.1) then
            ! liquid water
            if (debug.gt.2) write(*,*) k,p,T,pr%rb(k) 
            call getProp(p,min(T,645.0_dp),pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,1) 
          else
            ! ice (here ice VI, ToDo: include Ih, III and V)
            if (p_GPa.lt.2.216_dp) then !ice VI
              call getProp(p,T,pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,2)
            else if (p_GPa.lt.47.0_dp) then !ice VII after Schlager et al., 2004
              call getProp(p,T,pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,3)
            else ! ice X
              call getProp(p,T,pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,4)
            endif
          endif
          pr%gr(k) = 0.0_dp
          pr%KT(k) = 0.0_dp
          pr%KS(k) = 0.0_dp
          pr%GS(k) = 0.0_dp
          pr%S(k) = 0.0_dp
          pr%ec(k) = 0.0_dp
        else if (pr%rb(k).gt.pr%Rc) then !mantle
          if (pr%rb(k).gt.pr%Rm-pE%Dcr0) then 
            !!!!!!!!!
            ! crust !
            !!!!!!!!!
            pr%m(k) = 2
            call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
            pr%k(k) = pI%k_cr
            res = assemblageThermoElastic(p_GPa,T,[diopside,anorthite],[0.5_dp,0.5_dp],error)
            if (error.eq.0) then
              pr%rho(k) = res(2)
              pr%alpha(k) = res(3)
              pr%Cp(k) = res(4)
              pr%gr(k) = res(5)
              pr%KT(k) = res(6)
              pr%KS(k) = res(7)
              pr%GS(k) = res(8)
              pr%S(k) = res(9)
            else
              pr%rho(k) = 0.0_dp
              pr%alpha(k) = 0.0_dp
              pr%Cp(k) = 0.0_dp
              pr%gr(k) = 0.0_dp
              pr%KT(k) = 0.0_dp
              pr%KS(k) = 0.0_dp
              pr%GS(k) = 0.0_dp
              pr%S(k) = 0.0_dp
            endif
            ! electrical conductivity same as for olivine mantle below for now
            if (pI%ec_model.eq.1) then
              pr%ec(k) = 10.0_dp**2.69_dp*exp(-1.62_dp/(kB*T)) ! Xu dry profile
            elseif (pI%ec_model.eq.2) then
              pr%ec(k) = 10.0_dp**4.73_dp*exp(-2.31_dp/(kB*T)) &  !electrical conductivity in upper mantle following Yoshino and Katsura 2013
                     & + 10.0_dp**1.9_dp*X_H2OM*exp(-(0.92_dp-0.16_dp*X_H2OM**(1.0_dp/3.0_dp))/(kB*T))
            else
              pr%ec(k) = 10.0_dp**4.73_dp*exp(-2.31_dp/(kB*T)) &  !electrical conductivity in upper mantle following Yoshino and Katsura 2013
                     & + 10.0_dp**2.98_dp*X_FeM_mol*exp(-(1.71_dp-2.0_dp*X_FeM_mol**(1.0_dp/3.0_dp))/(kB*T)) & !neglect above iron influence (seems to lead to different effect than in paper???)
                     & + 10.0_dp**1.9_dp*X_H2OM*exp(-(0.92_dp-0.16_dp*X_H2OM**(1.0_dp/3.0_dp))/(kB*T))
            endif
            

          elseif (pr%pb(k).le.(4.92_dp+5.4_dp*(1-X_FeM)+0.0025_dp*pr%Tb(k))*1.0e+9_dp) then ! from Doctoral Thesis of Monika Buske, University of Goettingen, Germany, 2006
            !!!!!!!!!!! 
            ! olivine !
            !!!!!!!!!!!
            pr%m(k) = 3
            call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
            !call getProp(p,T,pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,5) ! used only for pr%k
            !res = eos(p_GPa,T,forsterite,error)
            pr%k(k) = (2.47_dp + 0.33_dp*p_GPa)*(300.0_dp/T)**0.48_dp
            res = assemblageThermoElastic(p_GPa,T,[forsterite,fayalite],[1.0_dp-X_FeM,X_FeM],error)
            if (error.eq.0) then
              pr%rho(k) = res(2)
              pr%alpha(k) = res(3)
              pr%Cp(k) = res(4)
              pr%gr(k) = res(5)
              pr%KT(k) = res(6)
              pr%KS(k) = res(7)
              pr%GS(k) = res(8)
              pr%S(k) = res(9)
            else
              pr%rho(k) = 0.0_dp
              pr%alpha(k) = 0.0_dp
              pr%Cp(k) = 0.0_dp
              pr%gr(k) = 0.0_dp
              pr%KT(k) = 0.0_dp
              pr%KS(k) = 0.0_dp
              pr%GS(k) = 0.0_dp
              pr%S(k) = 0.0_dp
            endif
            if (pI%ec_model.eq.1) then
              pr%ec(k) = 10.0_dp**2.69_dp*exp(-1.62_dp/(kB*T)) ! Xu dry profile
            elseif (pI%ec_model.eq.2) then
              pr%ec(k) = 10.0_dp**4.73_dp*exp(-2.31_dp/(kB*T)) &  !electrical conductivity in upper mantle following Yoshino and Katsura 2013
                     & + 10.0_dp**1.9_dp*X_H2OM*exp(-(0.92_dp-0.16_dp*X_H2OM**(1.0_dp/3.0_dp))/(kB*T))
            else
              pr%ec(k) = 10.0_dp**4.73_dp*exp(-2.31_dp/(kB*T)) &  !electrical conductivity in upper mantle following Yoshino and Katsura 2013
                     & + 10.0_dp**2.98_dp*X_FeM_mol*exp(-(1.71_dp-2.0_dp*X_FeM_mol**(1.0_dp/3.0_dp))/(kB*T)) & !neglect above iron influence (seems to lead to different effect than in paper???)
                     & + 10.0_dp**1.9_dp*X_H2OM*exp(-(0.92_dp-0.16_dp*X_H2OM**(1.0_dp/3.0_dp))/(kB*T))
            endif

          elseif (pr%pb(k).le.(6.2_dp+16.7_dp*(1-X_FeM)+0.005_dp*(pr%Tb(k)-2260.0_dp))*1.0e+9_dp) then ! from Doctoral Thesis of Monika Buske, University of Goettingen, Germany, 2006
            !!!!!!!!!!!!!! 
            ! wadsleyite !
            !!!!!!!!!!!!!!
            pr%m(k) = 4
            call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
            !call getProp(p,T,pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,5) ! used only for pr%k
            !res = eos(p_GPa,T,forsterite,error)
            pr%k(k) = (3.81_dp + 0.34_dp*p_GPa)*(300.0_dp/T)**0.56_dp
            res = assemblageThermoElastic(p_GPa,T,[MgWadsleyite,FeWadsleyite],[1.0_dp-X_FeM,X_FeM],error)
            if (error.eq.0) then
              pr%rho(k) = res(2)
              pr%alpha(k) = res(3)
              pr%Cp(k) = res(4)
              pr%gr(k) = res(5)
              pr%KT(k) = res(6)
              pr%KS(k) = res(7)
              pr%GS(k) = res(8)
              pr%S(k) = res(9)
            else
              pr%rho(k) = 0.0_dp
              pr%alpha(k) = 0.0_dp
              pr%Cp(k) = 0.0_dp
              pr%gr(k) = 0.0_dp
              pr%KT(k) = 0.0_dp
              pr%KS(k) = 0.0_dp
              pr%GS(k) = 0.0_dp
              pr%S(k) = 0.0_dp
            endif
            if (pI%ec_model.eq.1) then
              pr%ec(k) = 10.0_dp**3.29_dp*exp(-1.29_dp/(kB*T)) !Xu profile
            elseif (pI%ec_model.eq.2) then
              pr%ec(k) = 399.0_dp*exp(-1.49_dp/(kB*T))+7.74_dp*X_H2OM*exp(-(0.68_dp-0.02_dp*X_H2OM**(1.0_dp/3.0_dp))/(kB*T)) ! Yoshino 2008
            else
              pr%ec(k) = 10.0_dp**2.46_dp*X_FeM_mol*exp(-(1.45_dp-2.0_dp*X_FeM_mol**(1.0_dp/3.0_dp))/(kB*T)) & !Yoshino and Katsura 2013
                     & + 10.0_dp**1.4_dp*X_H2OM*exp(-(0.83_dp-0.2_dp*X_H2OM**(1.0_dp/3.0_dp))/(kB*T))
            endif

          elseif (pr%pb(k).le.(22.9_dp-0.0025_dp*(pr%Tb(k)-2260.0_dp))*1.0e+9_dp) then ! from Doctoral Thesis of Monika Buske, University of Goettingen, Germany, 2006
            !!!!!!!!!!!!!!! 
            ! ringwoodite !
            !!!!!!!!!!!!!!!
            pr%m(k) = 5
            call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
            !call getProp(p,T,pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,5) ! used only for pr%k
            !res = eos(p_GPa,T,forsterite,error)
            pr%k(k) = (3.52_dp + 0.36_dp*p_GPa)*(300.0_dp/T)**0.61_dp
            res = assemblageThermoElastic(p_GPa,T,[MgRingwoodite,FeRingwoodite],[1.0_dp-X_FeM,X_FeM],error)
            if (error.eq.0) then
              pr%rho(k) = res(2)
              pr%alpha(k) = res(3)
              pr%Cp(k) = res(4)
              pr%gr(k) = res(5)
              pr%KT(k) = res(6)
              pr%KS(k) = res(7)
              pr%GS(k) = res(8)
              pr%S(k) = res(9)
            else
              pr%rho(k) = 0.0_dp
              pr%alpha(k) = 0.0_dp
              pr%Cp(k) = 0.0_dp
              pr%gr(k) = 0.0_dp
              pr%KT(k) = 0.0_dp
              pr%KS(k) = 0.0_dp
              pr%GS(k) = 0.0_dp
              pr%S(k) = 0.0_dp
            endif
            if (pI%ec_model.eq.1) then
              pr%ec(k) = 10.0_dp**2.92_dp*exp(-1.16_dp/(kB*T)) ! Xu dry profile
            elseif (pI%ec_model.eq.2) then
              pr%ec(k) = 838.0_dp*exp(-1.36_dp/(kB*T))+27.7_dp*X_H2OM*exp(-(1.12_dp-0.67_dp*X_H2OM**(1.0_dp/3.0_dp))/(kB*T)) !Yoshino 2008
            elseif (pI%ec_model.eq.3) then
!             pr%ec(k) = 10.0_dp**2.92_dp*X_FeM_mol*exp(-(1.36_dp-2.0_dp*X_FeM_mol**(1.0_dp/3.0_dp))/(kB*T)) &
              pr%ec(k) = 10042.0_dp*X_FeM_mol*exp(-(2.48_dp-2.25_dp*X_FeM_mol**(1.0_dp/3.0_dp))/(kB*T)) & !enthalpy and alpha from Yoshino&Katsura 2009
                     & + 10.0_dp**1.44_dp*X_H2OM*exp(-(1.12_dp-0.67_dp*X_H2OM**(1.0_dp/3.0_dp))/(kB*T))
            else
              pr%ec(k) = 10.0_dp**2.92_dp*X_FeM_mol*exp(-(1.36_dp-2.0_dp*X_FeM_mol**(1.0_dp/3.0_dp))/(kB*T)) &
                     & + 10.0_dp**1.44_dp*X_H2OM*exp(-(1.12_dp-0.67_dp*X_H2OM**(1.0_dp/3.0_dp))/(kB*T))
            endif
          else if (pr%pb(k).le.(109.54_dp+0.0058_dp*pr%Tb(k))*1.0e+9_dp) then ! formula for plot in Murakami et al., 2004 based on data from Sidorin
            !!!!!!!!!!!!!!
            ! perovskite !
            !!!!!!!!!!!!!!
            pr%m(k) = 6
            call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
            !call getProp(p,T,pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,6) ! used only for pr%k
            !res = assemblageThermoElastic(p_GPa,T,[MgPerovskite,periclase],[0.7135_dp,0.2865_dp],error)
            pr%k(k) = (3.48_dp + 0.12_dp*p_GPa)*(300.0_dp/T)**0.31_dp
            res = assemblageThermoElastic(p_GPa,T,[MgPerovskite,periclase,FePerovskite,wustite],&
                  &[0.7135_dp*(1.0_dp-X_FeM),0.2865_dp*(1.0_dp-X_FeM),0.7135_dp*X_FeM,0.2865_dp*X_FeM],error)
            if (error.eq.0) then
              pr%rho(k) = res(2)
              pr%alpha(k) = res(3)
              pr%Cp(k) = res(4)
              pr%gr(k) = res(5)
              pr%KT(k) = res(6)
              pr%KS(k) = res(7)
              pr%GS(k) = res(8)
              pr%S(k) = res(9)
            else
              pr%rho(k) = 0.0_dp
              pr%alpha(k) = 0.0_dp
              pr%Cp(k) = 0.0_dp
              pr%gr(k) = 0.0_dp
              pr%KT(k) = 0.0_dp
              pr%KS(k) = 0.0_dp
              pr%GS(k) = 0.0_dp
              pr%S(k) = 0.0_dp
            endif
            ! electrical conductivity in lower mantle following Xu 2000
            ! use for now volume fractions of pe (0.3149) and pv (0.6851) at zero pressure condition
            ! ToDo: have volume fractions stored in res vector or calculate electrical conductivity in EOS file
            pr%ec(k) = 0.6851_dp*10.0_dp**1.12_dp*exp(-0.62_dp/(kB*T)) &
                   & + 0.3149_dp*10.0_dp**2.69_dp*exp(-0.85_dp/(kB*T))
          else
            !!!!!!!!!!!!!!!!!!!
            ! post-perovskite !
            !!!!!!!!!!!!!!!!!!!
            pr%m(k) = 7
            call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
            !call getProp(p,T,pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,7)
            !res = assemblageThermoElastic(p_GPa,T,[MgPostPerovskite,periclase],[0.7135_dp,0.2865_dp],error)
            pr%k(k) = (3.48_dp + 0.12_dp*p_GPa)*(300.0_dp/T)**0.31_dp
            res = assemblageThermoElastic(p_GPa,T,[MgPostPerovskite,periclase,FePostPerovskite,wustite],&
                  &[0.7135_dp*(1.0_dp-X_FeM),0.2865_dp*(1.0_dp-X_FeM),0.7135_dp*X_FeM,0.2865_dp*X_FeM],error)
            if (error.eq.0) then
              pr%rho(k) = res(2)
              pr%alpha(k) = res(3)
              pr%Cp(k) = res(4)
              pr%gr(k) = res(5)
              pr%KT(k) = res(6)
              pr%KS(k) = res(7)
              pr%GS(k) = res(8)
              pr%S(k) = res(9)
            else
              pr%rho(k) = 0.0_dp
              pr%alpha(k) = 0.0_dp
              pr%Cp(k) = 0.0_dp
              pr%gr(k) = 0.0_dp
              pr%KT(k) = 0.0_dp
              pr%KS(k) = 0.0_dp
              pr%GS(k) = 0.0_dp
              pr%S(k) = 0.0_dp
            endif
            ! use for now pv values also for ppv
            pr%ec(k) = 0.6851_dp*10.0_dp**1.12_dp*exp(-0.62_dp/(kB*T)) &
                   & + 0.3149_dp*10.0_dp**2.69_dp*exp(-0.85_dp/(kB*T))
          endif
        else ! core
          ! ToDo: add treatment for xs
          pr%m(k) = 8
          call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
          !call getProp(p_GPa,T,pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,8)
          pr%k(k) = 50.0_dp ! not used
          res_c = eos_core(p_GPa,T,error)
          if (error.eq.0) then
            pr%rho(k) = res_c(2)
            pr%alpha(k) = res_c(3)
            pr%Cp(k) = res_c(4)
            pr%gr(k) = res_c(5)
            pr%KT(k) = res_c(6)
            pr%KS(k) = 0.0_dp
            pr%GS(k) = 0.0_dp
            pr%S(k) = 0.0_dp
          else
            if (debug.ge.1) write(*,*) "Error ",error," occured in core EOS for p=",p_GPa," and T=",T ! Error can occur in first iterations, but goes away when solution converges
            pr%rho(k) = pr%rho(k+1)!0.0_dp
            pr%alpha(k) = pr%alpha(k+1)!0.0_dp
            pr%Cp(k) = pr%Cp(k+1)!0.0_dp
            pr%gr(k) = pr%gr(k+1)!0.0_dp
            pr%KT(k) = pr%KT(k+1)!0.0_dp