anchor.diff:  END DO
   
   END SUBROUTINE rod_eject
  +
  +
  +SUBROUTINE trod_eject()
  +
  +!
  +! Purpose:
  +!    To perform rod ejection simulation
  +!
  +
  +USE sdata, ONLY: ng, nnod, sigr, nuf, sigs, f0, &
  +                 iBeta, lamb, nf, nout, nac, &
  +                 ttot, tdiv, tstep1, tstep2, &
  +                 ix, iy, iz, zdel,Ke, &
  +                 bcon, ftem, mtem, cden, &
  +                 fbpos, bpos, tmove, bspeed, mdir, &
  +                 nout, nac, nb, npow, ppow, pow, node_nf
  +USE InpOutp, ONLY: XS_updt, bther, ounit
  +USE nodal, ONLY: nodal_coup4, outer4, powdis
  +USE th, ONLY: th_iter, th_trans3
  +!, par_ave_f, par_max, par_ave
  +
  +IMPLICIT NONE
  +
  +REAL, DIMENSION(nnod,ng) ::  dsigr, dnuf
  +REAL, DIMENSION(nnod,ng,ng) :: dsigs
  +
  +REAL, DIMENSION(nnod,ng) ::  osigr, onuf
  +REAL, DIMENSION(nnod,ng,ng) :: osigs
  +
  +REAL, DIMENSION(nnod) :: pline       ! Linear power density
  +REAL, DIMENSION(nnod,ng) :: fl
  +
  +REAL :: A, rho
  +
  +REAL :: t1, t2
  +REAL :: xppow
  +INTEGER :: n, i, j, imax, step
  +REAL, PARAMETER :: hp = 0.0001 ! Point Kinetetic Time step
  +LOGICAL :: stime
  +
  +
  +! Allocate node power distribution npow
  +ALLOCATE(npow(nnod))
  +
  +! Calculate forward flux at t=0
  +CALL th_iter(0)
  +fl = f0
  +
  +!Initial amplitude function
  +amp = 1.
  +
  +! Save old sigr, nuf and sigs
  +osigr = sigr; onuf = nuf; osigs = sigs
  +dsigr = 0.; dnuf = 0.; dsigs = 0.
  +
  +WRITE(ounit, *)
  +WRITE(ounit, *) " TRANSIENT RESULTS :"
  +WRITE(ounit, *)
  +WRITE(ounit, *) " Step  Time(s)  React.($)    Power (W)    CR Bank Pos. (1-end)"
  +WRITE(ounit, *) "--------------------------------------------------------------"
  +WRITE(ounit,'(I4, F10.3, F10.4, ES15.4, 12F9.2)') 0, 0., 0., amp*ppow*0.01, (bpos(n), n = 1, nb)
  +
  +! Adjoint flux at t=0
  +CALL adj_calc()
  +
  +! Calculate intgral kinet parameters at t = 0
  +CALL kinet_par(dsigr, dnuf, dsigs, fl, A, rho)
  +
  +! Calculate Initial precursor density
  +DO j = 1, nf
  +    C(j) = beta(j) / (A * lamb(j))   ! See Eq. 6-32 D&H
  +END DO
  +
  +! Start transient calculation
  +step = 0
  +t2 = 0.
  +imax = CEILING(tdiv/tstep1)
  +stime = .FALSE.
  +
  +! First Time Step
  +DO i = 1, imax
  +
  +    step = step + 1
  +    t1 = t2
  +    t2 = t1 + tstep1
  +
  +    IF (t2 > tdiv) THEN
  +        t2 = tdiv
  +        stime = .TRUE.
  +    END IF
  +
  +    ! Rod bank changes
  +    DO n = 1, nb
  +        IF (mdir(n) == 1) THEN   ! If CR moving down
  +            IF (t2-tmove(n) > 1.d-5 .AND. fbpos(n)-bpos(n) < 1.d-5) THEN
  +                bpos(n) = bpos(n) - tstep1 *  bspeed(n)
  +                IF (bpos(n) < fbpos(n)) bpos(n) = fbpos(n)  ! If bpos exceed, set to fbpos
  +            END IF
  +        ELSE IF (mdir(n) == 2) THEN ! If CR moving up
  +            IF (t2-tmove(n) > 1.d-5 .AND. fbpos(n)-bpos(n) > 1.d-5) THEN
  +                bpos(n) = bpos(n) + tstep1 *  bspeed(n)
  +                IF (bpos(n) > fbpos(n)) bpos(n) = fbpos(n)  ! If bpos exceed, set to fbpos
  +            END IF
  +        ELSE
  +            CONTINUE
  +        END IF
  +     END DO
  +
  +    ! Calculate xsec after pertubation
  +    CALL XS_updt(bcon, ftem, mtem, cden, bpos)
  +
  +    ! Calculate shape function
  +    CALL nodal_coup4()
  +    CALL outer4(0)
  +
  +    ! Calculate xsec changes after rod is ejected
  +    dsigr = sigr - osigr
  +    dnuf = nuf - onuf
  +    dsigs = sigs - osigs
  +
  +    ! Calculate intgral kinet parameters
  +    CALL kinet_par(dsigr, dnuf, dsigs, f0, A, rho)
  +
  +    tbeta = 0.
  +    DO j = 1, nf
  +        tbeta = tbeta + beta(j)
  +    END DO
  +
  +    !Calculate amplitude function
  +    CALL point(t1, t2, hp, rho, A, beta, amp)
  +
  +    ! Calculate node power distribution
  +    CALL powdis(npow)
  +
  +    ! Power change
  +    xppow = ppow * amp * 0.01
  +
  +    ! Calculate linear power density for each nodes (W/cm)
  +    DO n = 1, nnod
  +        pline(n) = npow(n) * pow * xppow &
  +                 / (node_nf(ix(n),iy(n)) * zdel(iz(n)))     ! Linear power density (W/cm)
  +    END DO
  +
  +    ! TH transient
  +    CALL th_trans3(pline,tstep1)
  +
  +    WRITE(ounit,'(I4, F10.3, F10.4, ES15.4, 12F9.2)') step, t2, rho/tbeta, xppow, (bpos(n), n = 1, nb)
  +
  +    IF (stime) EXIT
  +
  +    IF (step>1000) THEN
  +        WRITE(ounit,*) 'TOO SMALL TIME STEPS. STOPPING'
  +        STOP
  +    END IF
  +
  +END DO
  +
  +
  +! Second Time Step
  +imax = CEILING((ttot-tdiv)/tstep2)
  +stime = .FALSE.
  +
  +DO i = 1, imax
  +
  +    step = step + 1
  +    t1 = t2
  +    t2 = t1 + tstep2
  +
  +    IF (t2 > ttot) THEN
  +        t2 = ttot
  +        stime = .TRUE.
  +    END IF
  +
  +    ! Rod bank changes
  +    DO n = 1, nb
  +        IF (mdir(n) == 1) THEN   ! If CR moving down
  +            IF (t2-tmove(n) > 1.d-5 .AND. fbpos(n)-bpos(n) < 1.d-5) THEN
  +                bpos(n) = bpos(n) - tstep2 *  bspeed(n)
  +                IF (bpos(n) < fbpos(n)) bpos(n) = fbpos(n)  ! If bpos exceed, set to fbpos
  +            END IF
  +        ELSE IF (mdir(n) == 2) THEN ! If CR moving up
  +            IF (t2-tmove(n) > 1.d-5 .AND. fbpos(n)-bpos(n) > 1.d-5) THEN
  +                bpos(n) = bpos(n) + tstep2 *  bspeed(n)
  +                IF (bpos(n) > fbpos(n)) bpos(n) = fbpos(n)  ! If bpos exceed, set to fbpos
  +            END IF
  +        ELSE
  +            CONTINUE
  +        END IF
  +     END DO
  +
  +    ! Calculate xsec after pertubation
  +    CALL XS_updt(bcon, ftem, mtem, cden, bpos)
  +
  +    ! Calculate shape function
  +    CALL nodal_coup4()
  +    CALL outer4(0)
  +
  +    ! Calculate xsec changes after rod is ejected
  +    dsigr = sigr - osigr
  +    dnuf = nuf - onuf
  +    dsigs = sigs - osigs
  +
  +    ! Calculate intgral kinet parameters
  +    CALL kinet_par(dsigr, dnuf, dsigs, f0, A, rho)
  +
  +    tbeta = 0.
  +    DO j = 1, nf
  +        tbeta = tbeta + beta(j)
  +    END DO
  +
  +    !Calculate amplitude function
  +    CALL point(t1, t2, hp, rho, A, beta, amp)
  +
  +    ! Calculate node power distribution
  +    CALL powdis(npow)
  +
  +    ! Power change
  +    xppow = ppow * amp * 0.01
  +
  +    ! Calculate linear power density for each nodes (W/cm)
  +    DO n = 1, nnod
  +        pline(n) = npow(n) * pow * xppow  &
  +                 / (node_nf(ix(n),iy(n)) * zdel(iz(n)))     ! Linear power density (W/cm)
  +    END DO
  +
  +    ! TH transient
  +    CALL th_trans3(pline,tstep2)
  +
  +    WRITE(ounit,'(I4, F10.3, F10.4, ES15.4, 12F9.2)') step, t2, rho/tbeta, xppow, (bpos(n), n = 1, nb)
  +
  +    IF (stime) EXIT
  +
  +    IF (step>1000) THEN
  +        WRITE(ounit,*) 'TOO SMALL TIME STEPS. STOPPING'
  +        STOP
  +    END IF
  +
  +END DO
  +
  +END SUBROUTINE trod_eject
  +
  +
  +SUBROUTINE trod_eject2()
  +
  +!
  +! Purpose:
  +!    To perform rod ejection simulation
  +!
  +
  +USE sdata, ONLY: ng, nnod, sigr, nuf, sigs, f0, &
  +                 iBeta, lamb, nf, nout, nac, &
  +                 ttot, tdiv, tstep1, tstep2, &
  +                 ix, iy, iz, zdel,Ke, &
  +                 bcon, ftem, mtem, cden, &
  +                 fbpos, bpos, tmove, bspeed, mdir, &
  +                 nout, nac, nb, npow, ppow, pow, node_nf
  +USE InpOutp, ONLY: XS_updt, bther, ounit
  +USE nodal, ONLY: nodal_coup4, outer4, powdis
  +USE th, ONLY: th_iter, th_trans4
  +!, par_ave_f, par_max, par_ave
  +
  +IMPLICIT NONE
  +
  +REAL, DIMENSION(nnod,ng) ::  dsigr, dnuf
  +REAL, DIMENSION(nnod,ng,ng) :: dsigs
  +
  +REAL, DIMENSION(nnod,ng) ::  osigr, onuf
  +REAL, DIMENSION(nnod,ng,ng) :: osigs
  +
  +REAL, DIMENSION(nnod) :: pline       ! Linear power density
  +REAL, DIMENSION(nnod,ng) :: fl
  +
  +REAL, DIMENSION(nnod) :: oftem, omtem, ocden
  +
  +REAL :: A, rho
  +
  +REAL :: t1, t2
  +REAL :: xppow
  +INTEGER :: n, i, j, imax, step
  +REAL, PARAMETER :: hp = 0.0001 ! Point Kinetetic Time step
  +LOGICAL :: stime
  +
  +
  +! Allocate node power distribution npow
  +ALLOCATE(npow(nnod))
  +
  +! Calculate forward flux at t=0
  +CALL th_iter(0)
  +fl = f0
  +
  +! Calculate node power distribution
  +CALL powdis(npow)
  +
  +!Initial amplitude function
  +amp = 1.
  +
  +! Save old sigr, nuf and sigs
  +osigr = sigr; onuf = nuf; osigs = sigs
  +dsigr = 0.; dnuf = 0.; dsigs = 0.
  +
  +WRITE(ounit, *)
  +WRITE(ounit, *) " TRANSIENT RESULTS :"
  +WRITE(ounit, *)
  +WRITE(ounit, *) " Step  Time(s)  React.($)    Power (W)    CR Bank Pos. (1-end)"
  +WRITE(ounit, *) "--------------------------------------------------------------"
  +WRITE(ounit,'(I4, F10.3, F10.4, ES15.4, 12F9.2)') 0, 0., 0., amp*ppow*0.01, (bpos(n), n = 1, nb)
  +
  +! Adjoint flux at t=0
  +CALL adj_calc()
  +
  +! Calculate intgral kinet parameters at t = 0
  +CALL kinet_par(dsigr, dnuf, dsigs, fl, A, rho)
  +
  +! Calculate Initial precursor density
  +DO j = 1, nf
  +    C(j) = beta(j) / (A * lamb(j))   ! See Eq. 6-32 D&H
  +END DO
  +
  +! Start transient calculation
  +step = 0
  +t2 = 0.
  +imax = CEILING(tdiv/tstep1)
  +stime = .FALSE.
  +
  +! First Time Step
  +DO i = 1, imax
  +
  +    step = step + 1
  +    t1 = t2
  +    t2 = t1 + tstep1
  +
  +    IF (t2 > tdiv) THEN
  +        t2 = tdiv
  +        stime = .TRUE.
  +    END IF
  +
  +    ! Rod bank changes
  +    DO n = 1, nb
  +        IF (mdir(n) == 1) THEN   ! If CR moving down
  +            IF (t2-tmove(n) > 1.d-5 .AND. fbpos(n)-bpos(n) < 1.d-5) THEN
  +                bpos(n) = bpos(n) - tstep1 *  bspeed(n)
  +                IF (bpos(n) < fbpos(n)) bpos(n) = fbpos(n)  ! If bpos exceed, set to fbpos
  +            END IF
  +        ELSE IF (mdir(n) == 2) THEN ! If CR moving up
  +            IF (t2-tmove(n) > 1.d-5 .AND. fbpos(n)-bpos(n) > 1.d-5) THEN
  +                bpos(n) = bpos(n) + tstep1 *  bspeed(n)
  +                IF (bpos(n) > fbpos(n)) bpos(n) = fbpos(n)  ! If bpos exceed, set to fbpos
  +            END IF
  +        ELSE
  +            CONTINUE
  +        END IF
  +     END DO
  +
  +     ! Power change
  +     xppow = ppow * amp * 0.01
  +
  +     ! Calculate linear power density for each nodes (W/cm)
  +     DO n = 1, nnod
  +         pline(n) = npow(n) * pow * xppow &
  +                  / (node_nf(ix(n),iy(n)) * zdel(iz(n)))     ! Linear power density (W/cm)
  +     END DO
  +
  +     !Save old th paramaters
  +     oftem = ftem; ocden = cden; omtem = mtem
  +
  +     ! TH transient
  +     CALL th_trans4(pline,tstep1)
  +
  +     ! Calculate xsec after pertubation
  +     CALL XS_updt(bcon, 0.5*(ftem+oftem), 0.5*(mtem+omtem), 0.5*(cden+ocden), bpos)
  +
  +    ! Calculate shape function
  +    CALL nodal_coup4()
  +    CALL outer4(0)
  +
  +    ! Calculate node power distribution
  +    CALL powdis(npow)
  +
  +    ! Calculate xsec changes after rod is ejected
  +    dsigr = sigr - osigr
  +    dnuf = nuf - onuf
  +    dsigs = sigs - osigs
  +
  +    ! Calculate intgral kinet parameters
  +    CALL kinet_par(dsigr, dnuf, dsigs, f0, A, rho)
  +
  +    tbeta = 0.
  +    DO j = 1, nf
  +        tbeta = tbeta + beta(j)
  +    END DO
  +
  +    !Calculate amplitude function
  +    CALL point(t1, t2, hp, rho, A, beta, amp)
  +
  +    WRITE(ounit,'(I4, F10.3, F10.4, ES15.4, 12F9.2)') step, t2, rho/tbeta, xppow, (bpos(n), n = 1, nb)
  +
  +    IF (stime) EXIT
  +
  +    IF (step>1000) THEN
  +        WRITE(ounit,*) 'TOO SMALL TIME STEPS. STOPPING'
  +        STOP
  +    END IF
  +
  +END DO
  +
  +
  +! Second Time Step
  +imax = CEILING((ttot-tdiv)/tstep2)
  +stime = .FALSE.
  +
  +DO i = 1, imax
  +
  +    step = step + 1
  +    t1 = t2
  +    t2 = t1 + tstep2
  +
  +    IF (t2 > ttot) THEN
  +        t2 = ttot
  +        stime = .TRUE.
  +    END IF
  +
  +    ! Rod bank changes
  +    DO n = 1, nb
  +        IF (mdir(n) == 1) THEN   ! If CR moving down
  +            IF (t2-tmove(n) > 1.d-5 .AND. fbpos(n)-bpos(n) < 1.d-5) THEN
  +                bpos(n) = bpos(n) - tstep2 *  bspeed(n)
  +                IF (bpos(n) < fbpos(n)) bpos(n) = fbpos(n)  ! If bpos exceed, set to fbpos
  +            END IF
  +        ELSE IF (mdir(n) == 2) THEN ! If CR moving up
  +            IF (t2-tmove(n) > 1.d-5 .AND. fbpos(n)-bpos(n) > 1.d-5) THEN
  +                bpos(n) = bpos(n) + tstep2 *  bspeed(n)
  +                IF (bpos(n) > fbpos(n)) bpos(n) = fbpos(n)  ! If bpos exceed, set to fbpos
  +            END IF
  +        ELSE
  +            CONTINUE
  +        END IF
  +     END DO
  +
  +     ! Power change
  +     xppow = ppow * amp * 0.01
  +
  +     ! Calculate linear power density for each nodes (W/cm)
  +     DO n = 1, nnod
  +         pline(n) = npow(n) * pow * xppow &
  +                  / (node_nf(ix(n),iy(n)) * zdel(iz(n)))     ! Linear power density (W/cm)
  +     END DO
  +
  +     !Save old th paramaters
  +     oftem = ftem; ocden = cden; omtem = mtem
  +
  +     ! TH transient
  +     CALL th_trans4(pline,tstep2)
  +
  +     ! Calculate xsec after pertubation
  +     CALL XS_updt(bcon, 0.5*(ftem+oftem), 0.5*(mtem+omtem), 0.5*(cden+ocden), bpos)
  +
  +    ! Calculate shape function
  +    CALL nodal_coup4()
  +    CALL outer4(0)
  +
  +    ! Calculate node power distribution
  +    CALL powdis(npow)
  +
  +    ! Calculate xsec changes after rod is ejected
  +    dsigr = sigr - osigr
  +    dnuf = nuf - onuf
  +    dsigs = sigs - osigs
  +
  +    ! Calculate intgral kinet parameters
  +    CALL kinet_par(dsigr, dnuf, dsigs, f0, A, rho)
  +
  +    tbeta = 0.
  +    DO j = 1, nf
  +        tbeta = tbeta + beta(j)
  +    END DO
  +
  +    !Calculate amplitude function
  +    CALL point(t1, t2, hp, rho, A, beta, amp)
  +
  +    WRITE(ounit,'(I4, F10.3, F10.4, ES15.4, 12F9.2)') step, t2, rho/tbeta, xppow, (bpos(n), n = 1, nb)
  +
  +    IF (stime) EXIT
  +
  +    IF (step>1000) THEN
  +        WRITE(ounit,*) 'TOO SMALL TIME STEPS. STOPPING'
  +        STOP
  +    END IF
  +
  +END DO
  +
  +END SUBROUTINE trod_eject2
   
   
   SUBROUTINE point(ti, tf, h, xrho, xA, xbet, xamp)
  
positive.code: af = f0       ! Save adjoint flux to af
  
  ALLOCATE(beta(nf), C(nf))
  
  END SUBROUTINE adj_calc
  
  
negative.code: END DO
  
  DO g = 1, ng
        vdum5 = 0.
      DO h = 1, ng
          DO n = 1, nnod
              vdum5(n) = vdum5(n) + dsigs(n,h,g) * sf(n,h)
  