anchor.diff:              a(nt-1) = -xa
               b(nt-1) =  xa + xc + eta
               c(nt-1) = -xc
  -            d(nt-1) = pdens * 0.5 * (rf**2 - rpos(nt-2)**2) &
  +            d(nt-1) = pdens * 0.5d0 * (rf**2 - rpos(nt-2)**2) &
                       + eta * tfm(xyz(i,j,k),nt-1)
   
               ! For nt [GAP-CLADDING INTERFACE]
  -            tem = 0.5 * (tfm(xyz(i,j,k),nt) + tfm(xyz(i,j,k),nt+1))
  -            kt = getkc(tem)      ! For cladding
  +            kt1 = getkc(tfm(xyz(i,j,k),nt))
  +            kt2 = getkc(tfm(xyz(i,j,k),nt+1))
  +            kt  = 2.d0 * kt1 * kt2 / (kt1 + kt2)     ! For cladding
               cp = getcpc(tfm(xyz(i,j,k),nt))
               eta = cdens * cp * (rpos(nt)**2 - rg**2) / (2. * h)
               xa = xc
  
positive.code: 
  INTEGER :: i, j, k, n
  DOUBLE PRECISION, DIMENSION(nt+1) :: a, b, c, d
  DOUBLE PRECISION :: hs, hg = 1.e4, kt           ! coolant heat transfer coef., gap heat transfer coef, and thermal conductivity
  DOUBLE PRECISION :: alpha = 0.7
  DOUBLE PRECISION :: xa, xc, tem
  DOUBLE PRECISION :: fdens = 10.412e3            ! UO2 density (kg/m3)
  DOUBLE PRECISION :: cdens = 6.6e3               ! Cladding density (kg/m3)
  DOUBLE PRECISION :: cp                          ! Specific heat capacity
  
negative.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
  