anchor.diff:  IMPLICIT NONE
   
   REAL, INTENT(IN) :: ent
  -REAL, INTENT(OUT) :: t, rho, prx, kvx
  -REAL :: t1, rho1, ent1, kv1, pr1
  -REAL :: t2, rho2, ent2, kv2, pr2
  +REAL, INTENT(OUT) :: t, rho, prx, kvx, tcx
  +REAL :: t1, rho1, ent1, kv1, pr1, tc1
  +REAL :: t2, rho2, ent2, kv2, pr2, tc2
   REAL :: ratx
   
  -IF ((ent < 857955) .OR. (ent > 1612217)) THEN
  +IF ((ent < 858341.5) .OR. (ent > 1624307.1)) THEN
       WRITE(ounit,*) '  Enthalpy. : ', ent
       WRITE(ounit,*) '  ERROR : ENTHALPY IS OUT OF THE RANGE OF DATA IN THE STEAM TABLE'
       WRITE(ounit,*) '  CHECK INPUT MASS FLOW RATE OR POWER'
  
positive.code:                   + cf * xpline(xyz(i,j,k)) * 100.                              ! Coolant Linear power densisty (W/m)
              vol   = farea * zdel(k) * 0.01
              IF (k == 1) THEN                                                    ! Calculate coolant enthalpy
                  eta = 0.5 * (enti + entp(xyz(i,j,k)))
                  ent(xyz(i,j,k)) = (cpline * zdel(k) * 0.01 * h &
                                  + (cflow * h - 0.5 * mdens * vol) * enti + mdens * vol * eta) &
                                  / (0.5 * mdens * vol + cflow * h)
                  CALL gettd(0.5 * (enti + ent(xyz(i,j,k))), mtem(xyz(i,j,k)), &
                            cden(xyz(i,j,k)), Pr, kv)                             ! Get corresponding temp and density
              ELSE
                  eta = 0.5 * (entp(xyz(i,j,k-1)) + entp(xyz(i,j,k)))
                  ent(xyz(i,j,k)) = (cpline * zdel(k) * 0.01 * h &
                                  + (cflow * h - 0.5 * mdens * vol) &
                                  * ent(xyz(i,j,k-1)) + mdens * vol * eta) &
                                  / (0.5 * mdens * vol + cflow * h)
                  CALL gettd(0.5 * (ent(xyz(i,j,k-1)) + ent(xyz(i,j,k))), &
                             mtem(xyz(i,j,k)), cden(xyz(i,j,k)), Pr, kv)          ! Get corresponding temp and density
              END IF
  
  
  
negative.code:             pdens = (1. - cf) * 100. * xpline(xyz(i,j,k)) / (pi * rf**2)                ! Fuel pin Power Density (W/m3)
  
              ! Calculate tridiagonal matrix: a, b, c and source: d
              ! For nt=1 [FUEL CENTERLINE]
              tem = 0.5 * (tfm(xyz(i,j,k),1) + tfm(xyz(i,j,k),2))                        ! Average temp. to get thermal conductivity
              kt = getkf(tem)                                                            ! Get thermal conductivity
              cp = getcpf(tfm(xyz(i,j,k),1))                                                           ! Get specific heat capacity
              eta = fdens * cp * rpos(1)**2 / (2. * h)
              xc  = kt * rpos(1) / rdel(1)
              b(1) =  xc + eta
              c(1) = -xc
              d(1) = pdens * 0.5 * rpos(1)**2 + eta * tfm(xyz(i,j,k),1)
  
              DO n = 2, nt-2
                  tem = 0.5 * (tfm(xyz(i,j,k),n) + tfm(xyz(i,j,k),n+1))
                  kt = getkf(tem)
                  cp = getcpf(tfm(xyz(i,j,k),n))
                  eta = fdens * cp * (rpos(n)**2 - rpos(n-1)**2) / (2. * h)
                  xa = xc
                  xc = kt * rpos(n) / rdel(n)
                  a(n) = -xa
                  b(n) =  xa + xc + eta
                  c(n) = -xc
                  d(n) = pdens * 0.5 * (rpos(n)**2 - rpos(n-1)**2) &
                       + eta * tfm(xyz(i,j,k),n)
              END DO
  
              ! For nt-1 [FUEL SURFACE]
              cp = getcpf(tfm(xyz(i,j,k),nt-1))
              eta = fdens * cp * (rf**2 - rpos(nt-2)**2) / (2. * h)
              xa = xc
              xc = hg * rpos(nt)   ! This is position of inner clad
              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) &
                      + eta * tfm(xyz(i,j,k),nt-1)
  
              ! For nt [INNER CLAD]
              tem = 0.5 * (tfm(xyz(i,j,k),nt) + tfm(xyz(i,j,k),nt+1))
              kt = getkc(tem)      ! For cladding
              cp = getcpc(tfm(xyz(i,j,k),nt))
              eta = cdens * cp * (rpos(nt+1)**2 - rpos(nt)**2) / (2. * h)
              xa = xc
              xc = kt * rpos(nt+1) / rdel(nt)
              a(nt) = -xa
              b(nt) =  xa + xc + eta
              c(nt) = -xc
              d(nt) = eta * tfm(xyz(i,j,k),nt)
  
              ! For nt+1  [OUTER CLAD]
              cp = getcpc(tfm(xyz(i,j,k),nt+1))
              eta = cdens * cp * (rpos(nt+2)**2 - rpos(nt+1)**2) / (2. * h)
              xa = xc
              xc = hs * rpos(nt+2)
              a(nt+1) = -xa
              b(nt+1) =  xa + xc + eta
              d(nt+1) = hs * rpos(nt+2) * mtem(xyz(i,j,k)) &
                      + eta * tfm(xyz(i,j,k),nt+1)
  
              ! Solve tridiagonal matrix
              CALL TridiaSolve(a, b, c, d, tfm(xyz(i,j,k), :))
  