model1D.f90 77.3 KB
Newer Older
Lena Noack's avatar
Lena Noack committed
1
2
3
4
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
5
  use stixrudeMantleSpecies
Lena Noack's avatar
Lena Noack committed
6
7
8
9
10
11
12
13
14
  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
15
    real(dp),allocatable::gr(:),KT(:),KS(:),GS(:),S(:),ec(:)
Lena Noack's avatar
Lena Noack committed
16
17
18
19
    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
20
    real(dp)::DTc ! temperature jump (either fixed or determined via melting temperature)
Lena Noack's avatar
Lena Noack committed
21
22
23
24
25
    integer::nc,nm,nw
    real(dp)::Tc_av,Tm_av,Tw_av,rho_av ! output values, averaged value
  end type props1D

!
26
27
28
!  ----------- rb(nr) = Rp (planet surface radius), nr = nc+nm+nw
!       -      r(nr), rho(nr), v(nr), ...
!  ----------- rb(nr-1), Tb(nr-1)
Lena Noack's avatar
Lena Noack committed
29
30
!
!
31
!  rb(nc) = Rc, rb(nc+nm) = Rm
Lena Noack's avatar
Lena Noack committed
32
33
!
!
34
35
36
!  ----------- rb(1)
!       -      r(1) (core)
!  ----------- rb(0) = center planet, r=0
Lena Noack's avatar
Lena Noack committed
37
38
39
40
41
42
43
44
45
46
!

contains

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

47
  subroutine Init1DProf(pE,pF,pI,pM,pO,pT,pV,pr,debug,nr,tyr)
48
    type(param_E),intent(inout)::pE
Lena Noack's avatar
Lena Noack committed
49
50
51
52
53
54
55
56
    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
57
    real(dp),intent(in)::tyr
Lena Noack's avatar
Lena Noack committed
58

59
60
    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
61
    integer::i,k,kinv,lTBL_k,error,maxOI,iter,nc_old
62
    real(dp)::res(9),res_c(6)
Lena Noack's avatar
Lena Noack committed
63
64
65
66
67
68
69

    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))
70
    allocate(pr%gr(nr),pr%KT(nr),pr%KS(nr),pr%GS(nr),pr%S(nr),pr%ec(nr))
Lena Noack's avatar
Lena Noack committed
71
72
    allocate(pr%Vb(nr),pr%Vc(nr),pr%Lc(nr),pr%dr(nr))

73
74
75
    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

76
77
78
79
80
81
82
83
84
85
86
87
88
89
90

    ! 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


Lena Noack's avatar
Lena Noack committed
91
92
    M = pE%M_E * 5.9736e+24 ! Mass in g
    pi_ = acos(-1.0_dp)
93
    maxOI = 20 !100
Lena Noack's avatar
Lena Noack committed
94
95
96
97
98
99
100
101
102
103
104
105

    ! 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

106
107
108
109
110
111
112
113
    ! 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
114
115
116
117
118
119
120
121
    
    !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
122
123
124
    !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)
125
126
127
128
    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
129
130
    !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
131

Lena Noack's avatar
Lena Noack committed
132
    ! simplified scaling law for first radii, will be updated by Get1DStruct()
133

Lena Noack's avatar
Lena Noack committed
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
    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

164
165
166
    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...
167
168
169

    pr%DTc = 0.0_dp

Lena Noack's avatar
Lena Noack committed
170
171
172
173
174
175
176
    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

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

      p_GPa = pr%pb(k)*1.0e-9_dp
      T_K = pr%Tb(k)
Lena Noack's avatar
Lena Noack committed
181

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

Lena Noack's avatar
Lena Noack committed
184
185
      if (pr%rb(k).gt.pr%Rm) then
        pr%m(k) = 1 ! first assume water
186
        call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
Lena Noack's avatar
Lena Noack committed
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203

        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
204
        pr%m(k) = 3 ! first assume olivine
205
        call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
        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)
Lena Noack's avatar
Lena Noack committed
235
      else
236
237
238
239
240
        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)
Lena Noack's avatar
Lena Noack committed
241
      endif
242
243
244
245
246
247
248
      !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 )
249

250
251
252
253
254
255
      ! 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
256
257
258
259
260
261
262
        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)
Lena Noack's avatar
Lena Noack committed
263
      endif
264
265
266
267
268
269
      ! 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
270
271
272
273
274
275
276
277
278
279
280
281
      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
282
283
284
285
286
287
        ! 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
288
      endif
289
290
291
292
293
294
295
296
      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

297
298
      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)

Lena Noack's avatar
Lena Noack committed
299
300
    enddo

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

Lena Noack's avatar
Lena Noack committed
303
304
305
306
307
308
309
310
311
312
    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
313
    iter = 0
Lena Noack's avatar
Lena Noack committed
314
    if (debug.ge.2) write(*,*) "-----------------------------"
315
316
    do while ((abs(Tcore-Tcore_old).gt.0.0001_dp).and.(iter.le.maxOI))
      iter=iter+1
317
      nc_old=pr%nc
Lena Noack's avatar
Lena Noack committed
318
319
      ! get interior structure for T-profile depending on pE%H2O and pE%Fe
      if (debug.ge.1) write(*,*) "Get1DStruct for T-profile"
320
      if (Tcore_old.eq.0.0_dp) then !first call -> T-prof can be completely off, just few iterations
321
322
        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)
323
      else
324
325
        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)
326
      endif
Lena Noack's avatar
Lena Noack committed
327

328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
      ! 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
Lena Noack's avatar
Lena Noack committed
375
376
377
378
379
380
      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
381
382
    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)
Lena Noack's avatar
Lena Noack committed
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438


    !!!!!!!!!!!!!!!!!!!!!!!!!!
    ! 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
439
    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
Lena Noack's avatar
Lena Noack committed
440
441
442
443
444
445
446
447
448
449
450
451
452
453


    ! 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

454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
    ! 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

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

Lena Noack's avatar
Lena Noack committed
490
491
492
493
  end subroutine Init1DProf


  ! is called both from Init and running simulation to obtain changes in radii due to thermal contraction
494
495
  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)
Lena Noack's avatar
Lena Noack committed
496
    type(props1D),intent(inout)::pr
497
    real(dp),intent(in)::M_E,X_Fe,X_H2O,rho_c,rho_m,rho_w,Psurf,Ttop,X_FeM,X_H2OM_w,tyr
Lena Noack's avatar
Lena Noack committed
498
    real(dp),intent(inout)::xs(:)
499
    type(param_E),intent(inout)::pE
500
    type(param_I),intent(inout)::pI
Lena Noack's avatar
Lena Noack committed
501
    real(dp),intent(out)::R_E !radius in Earth radii
502
    integer,intent(in)::nr,debug,writeIS,maxOI,iterO
Lena Noack's avatar
Lena Noack committed
503

504
    integer::iter,i,k,kinv,lTBL_k,error,nc_old
505
506
    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
507
    real(dp)::HtU235,HtU238,HtTh232,HtK40
508
    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
Lena HPC's avatar
Lena HPC committed
509
    real(dp)::X_FeM_mol ! molar fraction of iron in mantle
Lena Noack's avatar
Lena Noack committed
510
    real(dp),allocatable::mass(:)
511
    real(dp)::res(9),res_c(6)
Lena Noack's avatar
Lena Noack committed
512
513
514
515
516
    character(4)::str_M
    character(7)::str_Fe,str_H2O

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

517
!    maxOI = 100 ! in input.txt?
Lena Noack's avatar
Lena Noack committed
518
519
520
521
522
    tolOI = 1.0e-8_dp
    error = 0

    M = M_E * 5.9736e+24 ! Mass in g
    pi_ = acos(-1.0_dp)  
523
    kB = 8.61733e-5_dp !eV/K - Boltzman constant
Lena HPC's avatar
Lena HPC committed
524
    !X_H2OM = 0.0_dp !mantle water content in wt-%, for now use constant value here
Lena Noack's avatar
Lena Noack committed
525
526
527
528
529

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

530
531
532
533
534
    !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
535
536
537
538
539
!    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
!
Lena HPC's avatar
Lena HPC committed
540
    X_H2OM = X_H2OM_w*100.0_dp ! from weight value to weight percent as needed here
541
    !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
Lena HPC's avatar
Lena HPC committed
542

543
544
    ! 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
545
546
547
548
    !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)
549
    if (Mc.lt.0.0_dp) then
550
      write(*,*) "Error, not enough iron in planet to have ",pI%X_FeM," as iron number!"
551
552
      Mc=0.0_dp
    endif
553
554
    Mm = M - Mc - Mw
    !Mm = Mm + 2.0_dp*X_FeM*55.845_dp/(X_FeM*63.08_dp+140.69_dp)*Mm
555
556
    

Lena Noack's avatar
Lena Noack committed
557
558
559
560
561
562
563
564
565
566
567
568
!    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
569
570
571
    Rp_old2 = 0.0_dp
    Rp_old3 = 0.0_dp
    Rp_old4 = 0.0_dp
572
    Tc_old = 0.0_dp
Lena Noack's avatar
Lena Noack committed
573

574
575
!    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))
Lena Noack's avatar
Lena Noack committed
576
    do while ((iter.le.maxOI).and.((abs(Rp_old-pr%Rp)/max(Rp_old,pr%Rp)).gt.tolOI))
577
578
579
      Rp_old4 = Rp_old3
      Rp_old3 = Rp_old2
      Rp_old2 = Rp_old
Lena Noack's avatar
Lena Noack committed
580
      Rp_old = pr%Rp  
581
      Tc_old = pr%Tb(pr%nc)
Lena Noack's avatar
Lena Noack committed
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
      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
610
        call MeltPhase(pr%Tb(nr),pr%Tmb(nr),pr%pb(nr),pr%m(nr),pr%xs(nr),pI%X_FeM,pI%X_M0)
Lena Noack's avatar
Lena Noack committed
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629

        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
630
        pr%S(nr) = 0.0_dp
631
        pr%ec(nr) = 0.0_dp
Lena Noack's avatar
Lena Noack committed
632
      else if (X_Fe.lt.1.0_dp) then
633
634
635
         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
636
          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
637
638
639
640
          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
641
          call MeltPhase(pr%Tb(nr),pr%Tmb(nr),pr%pb(nr),pr%m(nr),pr%xs(nr),pI%X_FeM,pI%X_M0)
642
643
644
645
646
          !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
Lena Noack's avatar
Lena Noack committed
647
648
649
650
651
652
653
654
        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) 
655
          pr%S(nr) = res(9)
Lena Noack's avatar
Lena Noack committed
656
        else
657
658
659
          pr%rho(nr) = 0.0_dp
          pr%alpha(nr) = 0.0_dp
          pr%Cp(nr) = 0.0_dp
Lena Noack's avatar
Lena Noack committed
660
661
662
663
          pr%gr(nr) = 0.0_dp
          pr%KT(nr) = 0.0_dp
          pr%KS(nr) = 0.0_dp
          pr%GS(nr) = 0.0_dp
664
          pr%S(nr) = 0.0_dp
Lena Noack's avatar
Lena Noack committed
665
        endif
Lena Noack's avatar
Lena Noack committed
666
667
668
669
670
671
672
673
674
675
        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
Lena Noack's avatar
Lena Noack committed
676
      else
677
        pr%m(nr) = 8 ! Fe or FeS
678
        call MeltPhase(pr%Tb(nr),pr%Tmb(nr),pr%pb(nr),pr%m(nr),pr%xs(nr),pI%X_FeM,pI%X_M0)
679
680
        !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
Lena Noack's avatar
Lena Noack committed
681
682
683
684
685
686
687
688
689
        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
690
          pr%S(nr) = 0.0_dp
Lena Noack's avatar
Lena Noack committed
691
        else
692
693
694
          pr%rho(nr) = 0.0_dp
          pr%alpha(nr) = 0.0_dp
          pr%Cp(nr) = 0.0_dp
Lena Noack's avatar
Lena Noack committed
695
696
697
698
          pr%gr(nr) = 0.0_dp
          pr%KT(nr) = 0.0_dp
          pr%KS(nr) = 0.0_dp
          pr%GS(nr) = 0.0_dp
699
          pr%S(nr) = 0.0_dp
Lena Noack's avatar
Lena Noack committed
700
        endif
701
        pr%ec(nr) = 0.0_dp
Lena Noack's avatar
Lena Noack committed
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
      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
725
          call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
Lena Noack's avatar
Lena Noack committed
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746

          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
747
          pr%S(k) = 0.0_dp
748
          pr%ec(k) = 0.0_dp
Lena Noack's avatar
Lena Noack committed
749
        else if (pr%rb(k).gt.pr%Rc) then !mantle
750
751
752
753
754
          if (pr%rb(k).gt.pr%Rm-pE%Dcr0) then 
            !!!!!!!!!
            ! crust !
            !!!!!!!!!
            pr%m(k) = 2
755
            call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
            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
Lena Noack's avatar
Lena Noack committed
791
792
793
794
            !!!!!!!!!!! 
            ! olivine !
            !!!!!!!!!!!
            pr%m(k) = 3
795
            call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
            !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
Lena Noack's avatar
Lena Noack committed
819
820
821
822
823
824
825
826
827
828
            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
829

830
831
832
833
834
          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
835
            call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
836
837
838
839
            !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)
Lena Noack's avatar
Lena Noack committed
840
841
842
843
844
845
846
847
            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)
848
              pr%S(k) = res(9)
Lena Noack's avatar
Lena Noack committed
849
            else
850
851
852
853
854
855
856
857
858
              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
Lena Noack's avatar
Lena Noack committed
859
860
861
862
863
864
865
866
            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
867

868
869
870
871
872
          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
873
            call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
            !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
Lena Noack's avatar
Lena Noack committed
891
892
893
894
              pr%gr(k) = 0.0_dp
              pr%KT(k) = 0.0_dp
              pr%KS(k) = 0.0_dp
              pr%GS(k) = 0.0_dp
895
              pr%S(k) = 0.0_dp
Lena Noack's avatar
Lena Noack committed
896
            endif
Lena Noack's avatar
Lena Noack committed
897
898
899
900
901
902
903
904
905
906
907
908
            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
Lena Noack's avatar
Lena Noack committed
909
910
911
912
          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 !
            !!!!!!!!!!!!!!
913
            pr%m(k) = 6
914
            call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
915
916
917
918
919
            !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)
Lena Noack's avatar
Lena Noack committed
920
921
922
923
924
925
926
927
            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)
928
              pr%S(k) = res(9)
Lena Noack's avatar
Lena Noack committed
929
            else
930
931
932
              pr%rho(k) = 0.0_dp
              pr%alpha(k) = 0.0_dp
              pr%Cp(k) = 0.0_dp
Lena Noack's avatar
Lena Noack committed
933
934
935
936
              pr%gr(k) = 0.0_dp
              pr%KT(k) = 0.0_dp
              pr%KS(k) = 0.0_dp
              pr%GS(k) = 0.0_dp
937
              pr%S(k) = 0.0_dp
Lena Noack's avatar
Lena Noack committed
938
            endif
939
940
941
942
943
            ! 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))
Lena Noack's avatar
Lena Noack committed
944
945
946
947
          else
            !!!!!!!!!!!!!!!!!!!
            ! post-perovskite !
            !!!!!!!!!!!!!!!!!!!
948
            pr%m(k) = 7
949
            call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
950
951
952
953
954
            !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)
Lena Noack's avatar
Lena Noack committed
955
956
957
958
959
960
961
962
            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)
963
              pr%S(k) = res(9)
Lena Noack's avatar
Lena Noack committed
964
            else
965
966
967
              pr%rho(k) = 0.0_dp
              pr%alpha(k) = 0.0_dp
              pr%Cp(k) = 0.0_dp
Lena Noack's avatar
Lena Noack committed
968
969
970
971
              pr%gr(k) = 0.0_dp
              pr%KT(k) = 0.0_dp
              pr%KS(k) = 0.0_dp
              pr%GS(k) = 0.0_dp
972
              pr%S(k) = 0.0_dp
Lena Noack's avatar
Lena Noack committed
973
            endif
974
975
976
            ! 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))
Lena Noack's avatar
Lena Noack committed
977
978
979
          endif
        else ! core
          ! ToDo: add treatment for xs
980
          pr%m(k) = 8
981
          call MeltPhase(pr%Tb(k),pr%Tmb(k),pr%pb(k),pr%m(k),pr%xs(k),pI%X_FeM,pI%X_M0)
982
          !call getProp(p_GPa,T,pr%rho(k),pr%alpha(k),pr%Cp(k),pr%k(k),debug,error,8)
983
          pr%k(k) = 50.0_dp ! not used
Lena Noack's avatar
Lena Noack committed
984
985
986
987
988
989
990
991
992
          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
993
            pr%S(k) = 0.0_dp
Lena Noack's avatar
Lena Noack committed
994
          else
995
996
997
998
999
1000
            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
For faster browsing, not all history is shown. View entire blame