anchor.diff:  
   IF (afrad == 1) CALL AsmFlux(f0, 1.e0)
   
  -! CALL par_ave_f(ftem, tf)
  -! CALL par_ave(mtem, tm)
  -
  -! CALL par_max(ftem, maxtf)
  -! CALL par_max(tfm(:,1), maxfcl)
  -! CALL par_max(mtem, maxtm)
  -CALL getfq(npow)
  -
  -! Write Output
  -WRITE(ounit,*)
  -WRITE(ounit, 5001) tf, tf-273.15
  -WRITE(ounit, 5002)  maxfcl, maxfcl-273.15
  -WRITE(ounit, 5003) tm, tm-273.15
  -WRITE(ounit, 5004) maxtm, maxtm-273.15
  -
  -5001 FORMAT(2X, 'AVERAGE DOPPLER TEMPERATURE     : ', F7.1, ' K (', F7.1, ' C)')
  -5002 FORMAT(2X, 'MAX FUEL CENTERLINE TEMPERATURE : ', F7.1, ' K (', F7.1, ' C)')
  -5003 FORMAT(2X, 'AVERAGE MODERATOR TEMPERATURE   : ', F7.1, ' K (', F7.1, ' C)')
  -5004 FORMAT(2X, 'MAXIMUM MODERATOR TEMPERATURE   : ', F7.1, ' K (', F7.1, ' C)')
  -
  -! bpos(8) = 228.
  -
  -! CALL XS_updt(bcon, ftem, mtem, cden, bpos)
  -
  -! CALL th_iter()
  -! CALL outer4(1,1000,20)
  -
   
   END SUBROUTINE cbsearcht
   
   
   
  +
  +
  +
   END MODULE th
  
positive.code: 
  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 InpOutp, ONLY: XS_updt, ounit
  
  IMPLICIT NONE
  
  REAL, DIMENSION(nnod) :: pline
  REAL, DIMENSION(nnod) :: otem
  INTEGER :: n, niter
  
negative.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
          !IF (node_nf(i,j) < 2) STOP
  
              cpline = heatf(xyz(i,j,k)) * pi * dia  &
                     + cf * xpline(xyz(i,j,k)) * 100.                             ! Coolant Linear power densisty (W/m)
  