anchor.diff:  !    To update thermal parameters
   !
   
  -USE sdata, ONLY: mtem, cden, ftem, tin, xyz, cflow, nyy, nxx, nzz, cf, ent, heatf, &
  -                 ystag, tfm, nt, rpos, rdel, rf, rg, rc, pi, zdel, dia, ystag
  +USE sdata, ONLY: nnod, mtem, cden, ftem, tin, ix, iy, iz, nxx, nyy, cflow, cf,  &
  +                 ent, heatf, tfm, nt, rpos, rdel, rf, rg, rc, pi, zdel, dia
   
   IMPLICIT NONE
   
   REAL(DP), DIMENSION(:), INTENT(IN) :: xpline    ! Linear Power Density (W/cm)
   
  -INTEGER :: i, j, k, n
  +INTEGER :: i, n
   REAL(DP), DIMENSION(nt+1) :: a, b, c, d
   REAL(DP) :: hs, Hg = 1.d4, kt, kt1, kt2
   REAL(DP) :: alp = 0.7_DP
  
positive.code: !    To perform fuel pin thermal transient
  !
  
  USE sdata, ONLY: mtem, cden, ftem, tin, xyz, cflow, nyy, nzz, nxx, cf, ent, heatf, nnod, &
                   ystag, tfm, nt, rpos, rdel, rf, rg, rc, farea, dia, pi, zdel, ystag
  
  IMPLICIT NONE
  
  REAL(DP), DIMENSION(:), INTENT(IN) :: xpline    ! Linear Power Density (W/cm)
  REAL(DP), INTENT(IN) :: h                       ! Time step
  
  INTEGER :: i, j, k, n
  REAL(DP), DIMENSION(nt+1) :: a, b, c, d
  REAL(DP) :: hs, hg = 1.d4, kt , kt1, kt2          ! coolant heat transfer coef., gap heat transfer coef, and thermal conductivity
  REAL(DP) :: alpha = 0.7_DP
  
negative.code: 
  CALL getent(tin, enti)
  
  DO n = 1, nnod
     cpline = heatf(n) * pi * dia  + cf * xpline(n) * 100._DP          ! Coolant Linear power densisty (W/m)
     IF (iz(n) == 1) THEN                                              ! For most bootom channel
        ent(n) = enti + 0.5_DP * cpline * zdel(iz(n)) * 0.01_DP / cflow    ! Calculate coolant enthalpy
        CALL gettd(ent(n), mtem(n), cden(n), Pr, kv, tcon)             ! Get corresponding temp and density
        entm(ix(n),iy(n)) = 2._DP * ent(n) - enti                      ! Extrapolate enthalpy at node boundary
     ELSE
        ent(n) = entm(ix(n),iy(n)) + 0.5_DP * cpline * zdel(iz(n)) * 0.01_DP / cflow
        CALL gettd(ent(n), mtem(n), cden(n), Pr, kv, tcon)
        entm(ix(n),iy(n)) = 2._DP * ent(n) - entm(ix(n),iy(n))
     END IF
  
     hs = geths(cden(n), Pr, kv, tcon)
     pdens = (1. - cf) * 100._DP * xpline(n) / (pi * rf**2)        ! Fuel pin Power Density (W/m3)
  
     ! Calculate tridiagonal matrix: a, b, c and source: d
     ! For nt=1 [FUEL CENTERLINE]
     kt1 = getkf(tfm(n,1))                                         ! Get thermal conductivity
     kt2 = getkf(tfm(n,2))
     kt  = 2._DP * kt1 * kt2 / (kt1 + kt2)
     xc  = kt * rpos(1) / rdel(1)
     b(1) =  xc
     c(1) = -xc
     d(1) = pdens * 0.5_DP * rpos(1)**2
  
     DO i = 2, nt-2
        kt1 = kt2
        kt2 = getkf(tfm(n,i+1))
        kt  = 2._DP * kt1 * kt2 / (kt1 + kt2)
        xa = xc
        xc = kt * rpos(i) / rdel(i)
        a(i) = -xa
        b(i) =  xa + xc
        c(i) = -xc
        d(i) = pdens * 0.5_DP * (rpos(i)**2 - rpos(i-1)**2)
     END DO
  
     ! For nt-1 [FUEL-GAP INTERFACE]
     xa = xc
     xc = rg * Hg
     a(nt-1) = -xa
     b(nt-1) =  xa + xc
     c(nt-1) = -xc
     d(nt-1) = pdens * 0.5_DP * (rf**2 - rpos(nt-2)**2)
  
     ! For nt [GAP-CLADDING INTERFACE]
     kt1 = getkc(tfm(n,nt))
     kt2 = getkc(tfm(n,nt+1))
     kt  = 2._DP * kt1 * kt2 / (kt1 + kt2)     ! For cladding
     xa = xc
     xc = kt * rpos(nt) / rdel(nt)
     a(nt) = -xa
     b(nt) =  xa + xc
     c(nt) = -xc
     d(nt) = 0.
  
     ! For nt+1  [CLADDING-COOLANT INTERFACE]
     xa = xc
     a(nt+1) = -xa
     b(nt+1) =  xa + hs * rc
     d(nt+1) = rc * hs * mtem(n)
  
     ! Solve tridiagonal matrix
     CALL TridiaSolve(a, b, c, d, tfm(n, :))
  
     ! Get lumped fuel temp
     ftem(n) = (1.-alp) * tfm(n, 1) + alp * tfm(n, nt-1)
  
     ! Calculate heat flux
     heatf(n) = hs * (tfm(n, nt+1) - mtem(n))
  END DO
  
  
  