anchor.diff:  
       ! Modify removal xsec
       sigrp = sigr    ! Save sigr to sigrp
  -    IF (i > 1) THEN
  -      DO g = 1, ng
  -         DO n = 1, nnod
  -           omeg(n,g) = LOG(f0(n,g) / ft(n,g)) / tstep1
  -           sigr(n,g) = sigr(n,g) + 1. / (velo(g) * tstep1) + omeg(n,g) / velo(g)
  -         END DO
  -      END DO
  -    ELSE
  -       DO g = 1, ng
  -          DO n = 1, nnod
  -             omeg(n,g) = 0.
  -             sigr(n,g) = sigr(n,g) + 1. / (velo(g) * tstep1)
  -          END DO
  +    DO g = 1, ng
  +       DO n = 1, nnod
  +          sigr(n,g) = sigr(n,g) + 1. / (velo(g) * tstep1) + omeg(n,g) / velo(g)
          END DO
  -    END IF
  +    END DO
   
       ! Save the previous flux and DN precusor density
       ft = f0
  
positive.code:     step = step + 1
      t1 = t2
      t2 = REAL(i)*tstep1
  
      ! Rod bank changes
      DO n = 1, nb
  
negative.code:   tbeta = tbeta + iBeta(j)
  END DO
  
  ! Calculate reactivity
  CALL react(af, sigr, rho)
  
  