Update Interior Structure model to include TBL for initial temperature profile

src/model1D.f90

@@ -58,7 +58,7 @@ contains
 
     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::k,kinv,error,maxOI,iter,nc_old
+    integer::i,k,kinv,lTBL_k,error,maxOI,iter,nc_old
     real(dp)::res(9),res_c(6)
 
     if (pE%IntStruct.ge.0) &
@@ -247,8 +247,12 @@ contains
       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 ((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*(abs(pE%Tm)-pT%Ttop)
-      if ((pr%rb(k).gt.pr%Rm-pE%Dl).and.(pr%rb(k-1).le.pr%Rm-pE%Dl)) then
+      ! 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)
@@ -257,6 +261,12 @@ contains
         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
@@ -269,6 +279,12 @@ contains
             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)
@@ -366,8 +382,6 @@ contains
                     &pI%X_FeM,pI%X_H2OM,tyr,pE,pI,maxOI,iter)
 
 
-    pr%Dl = pE%Dl ! later: depending on viscosity profile and thus velocity profile
-
     !!!!!!!!!!!!!!!!!!!!!!!!!!
     ! Heat sources in mantle !
     !!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -487,7 +501,7 @@ contains
     real(dp),intent(out)::R_E !radius in Earth radii
     integer,intent(in)::nr,debug,writeIS,maxOI,iterO
 
-    integer::iter,k,kinv,error,nc_old
+    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
@@ -1152,6 +1166,7 @@ contains
 
       ! Test here already update temp profile:
       nc_old = pr%nc+1
+      lTBL_k = pr%nc+2 ! TBL thickness if problem to identify TBL below
       do kinv = 1,nr
         k = nr+1-kinv !do k=nr,1
         if (pr%alpha(k).eq.0.0_dp) then
@@ -1170,8 +1185,9 @@ contains
         pr%Tb(k-1) = damp * pr%Tb(k) * exp(pr%alpha(k)*pr%gb(k)/pr%Cp(k)*(pr%rb(k)-pr%rb(k-1))) &   ! T = T_t exp(alpha_t g / Cp_t * dr )
                    & + (1.0_dp-damp) * pr%Tb(k-1)
         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) = Ttop + (pr%Rm-pr%rb(k))/pE%Dl*(abs(pE%Tm)-Ttop)
-        if ((pr%rb(k).gt.pr%Rm-pE%Dl).and.(pr%rb(k-1).le.pr%Rm-pE%Dl)) then
+        if ((pr%rb(k).le.pr%Rm).and.(pr%rb(k).gt.pr%Rm-pr%Dl)) pr%Tb(k-1) = Ttop + (pr%Rm-pr%rb(k))/pr%Dl*(abs(pE%Tm)-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)
@@ -1180,7 +1196,12 @@ contains
           endif
           !write(*,*) "set upper mantle with Tm=",pE%Tm
         endif
-        !if ((pr%rb(k-1).le.pr%Rc).and.(pr%rb(k).gt.pr%Rc)) then
+        ! 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 (k.eq.pr%nc+1) then
           ! add temperature jump at CMB
           if (pE%DTc.gt.0) then
@@ -1197,6 +1218,12 @@ contains
               write(*,*) "Error, molten frozen at CMB, no temp jump applied",k-1,k,nc_old,pr%nc,pr%Tmb(nc_old),pr%Tb(nc_old),pE%Tm
             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
 
       enddo