anchor.diff:  
           END DO
       END DO
  -
  -
  -!STOP
   END DO
   
   
  
positive.code:                 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
  
  
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
  