anchor.diff:              a(nt-1) = -xa
               b(nt-1) =  xa + xc
               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)
   
               ! 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
               xa = xc
               xc = kt * rpos(nt) / rdel(nt)
               a(nt) = -xa
  
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:             a(nt-1) = -xa
              b(nt-1) =  xa + xc + eta
              c(nt-1) = -xc
              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]
              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
  