anchor.diff:  
   USE sdata, ONLY: nnod, ftem, mtem, cden, bcon, bpos, npow, pow, ppow, nout, nac,  &
                    zdel, node_nf, ix, iy, iz, th_err, node_nf, ix, iy, iz, th_niter
  -USE nodal, ONLY: nodal_coup4, outer4, PowDis
  +USE nodal, ONLY: nodal_coup4, outer4th, PowDis
   USE InpOutp, ONLY: XS_updt, ounit
   
   IMPLICIT NONE
   
  +INTEGER, INTENT(IN), OPTIONAL :: ind    ! if iteration reaching th_iter: ind = 0 => stop
   REAL, DIMENSION(nnod) :: pline
   REAL, DIMENSION(nnod) :: otem
   INTEGER :: n, niter
  
positive.code: DO k = 1, nzz
      DO j = 1, nyy
          DO i = ystag(j)%smin, ystag(j)%smax
          ! WRITE(ounit,*) cflow , node_nf(i,j), i, j
  
              cpline = heatf(xyz(i,j,k)) * pi * dia  &
                     + cf * xpline(xyz(i,j,k)) * 100.                             ! Coolant Linear power densisty (W/m)
  
negative.code:             IF (k == 1) THEN                                                    ! Calculate coolant enthalpy and
                  ent(xyz(i,j,k)) = enti + cpline * zdel(k) * 0.01 / cflow        ! corresponding temp and density
                  CALL gettd(0.5 * (enti + ent(xyz(i,j,k))), &
                             mtem(xyz(i,j,k)), cden(xyz(i,j,k)), Pr, kv, tc)          ! Get corresponding temp and density
              ELSE
                  ent(xyz(i,j,k)) = ent(xyz(i,j,k-1)) &
                                  + cpline * zdel(k) * 0.01 / cflow
                  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, tc)           ! Get corresponding temp and density
              END IF
  
              hs = geths(cden(xyz(i,j,k)))
  