Initial commit from Polytechnique Montreal

1 files changed, 501 insertions, 0 deletions

diff --git a/Dragon/src/EVOBLD.f b/Dragon/src/EVOBLD.f
new file mode 100644
index 0000000..76763ab
--- /dev/null
+++ b/Dragon/src/EVOBLD.f
@@ -0,0 +1,501 @@
+*DECK EVOBLD
+      SUBROUTINE EVOBLD(IMPX,INR,IGLOB,NBMIX,NBISO,NCOMB,ISONAM,YDPL,
+     1 VX,MILVO,JM,NVAR,NDFP,NSUPS,NREAC,NPAR,NFISS,XT,EPS1,EPS2,EXPMAX,
+     2 H1,ITYPE,IDIRAC,FIT,DELTA,ENERG,KPAR,BPAR,YIELD,IDR,RER,RRD,AWR,
+     3 FUELDN,SIG,VPH,VPHV,MIXPWR,VTOTD,IEVOLB,KFISS,KPF)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Perform flux normalization and call EVOSOL to solve the depletion
+* system for each depleting mixture between times XT(1) and XT(2).
+*
+*Copyright:
+* Copyright (C) 2002 Ecole Polytechnique de Montreal
+* This library is free software; you can redistribute it and/or
+* modify it under the terms of the GNU Lesser General Public
+* License as published by the Free Software Foundation; either
+* version 2.1 of the License, or (at your option) any later version
+*
+*Author(s): A. Hebert
+*
+*Parameters: input/output
+* IMPX    print flag (equal to zero for no print).
+* INR     type of flux normalization:
+*         =0: out-of-core depletion;
+*         =1: constant flux depletion;
+*         =2: constant fuel power depletion;
+*         =3: constant assembly power depletion.
+* IGLOB   out-of-fuel power in flux normalization. Compute the burnup:
+*         =-1: using the Serpent mode 0 empirical formula in the fuel;
+*         =0: using the power released in the fuel;
+*         =1: using the power released in the global geometry.
+* NBMIX   number of mixtures.
+* NBISO   number of isotopes/materials including non-depleting ones.
+* NCOMB   number of depleting mixtures.
+* ISONAM  alias name of isotopes.
+* YDPL    initial/final number density of isotope in the depletion
+*         chain. YDPL(NVAR+1,2,ICMB) is the stage burnup increment
+*         in region ICMB.
+* VX      volumes of the depleting mixtures.
+* MILVO   mixture index corresponding to each depleting mixture.
+* JM      position in isotope list of each nuclide of the depletion
+*         chain.
+* NVAR    number of depleting nuclides.
+* NDFP    number of direct fission products (fission fragments).
+* NSUPS   number of non-depleting isotopes producing energy.
+* NREAC   maximum number of depletion reactions.
+* NPAR    maximum number of parent nuclides in the depletion chain.
+* NFISS   number of fissile isotopes producing fission products.
+* XT      initial and final time (independent variable).
+* EPS1    required accuracy for the ode solver.
+* EPS2    required accuracy for constant power iterations.
+* EXPMAX  saturation limit. A nuclide is saturating if
+*         -ADPL(MU1(I))*(XT(2)-XT(1)).GT.EXPMAX. Suggested value:
+*         EXPMAX=80.0.
+* H1      guessed first stepsize.
+* ITYPE   type of ODE solution:
+*         =1 fifth-order Runge-Kutta method;
+*         =2 fourth-order Kaps-Rentrop method.
+* IDIRAC  saturation model flag (=1 to use Dirac function contributions
+*         in the saturating nuclide number densities).
+* FIT     flux normalization factor:
+*         n/cm**2/s if INR=1;
+*         MW/tonne of initial heavy elements if INR=2;
+*         W/cc of assembly volume if INR=3.
+* DELTA   burnup stage increments:
+*         DELTA(1): increment in fuel burnup for this stage;
+*         DELTA(2): increment in fuel neutron exposure for this stage;
+*         DELTA(3): target increment in fuel burnup for this stage.
+*         Cross section should be tabulated with respect to the sum
+*         of the DELTA(1) of all the previous stages.
+* ENERG   increment in fuel burnup for this stage in each mixture.
+* KPAR    position in chain of the parent nuclide and type of
+*         reaction.
+* BPAR    branching ratio for neutron induced reactions.
+* YIELD   mixture-dependent fission yields.
+* IDR     identifier for each depleting reaction.
+* RER     energy (Mev) per reaction. If RER(3,J)=0., the fission energy
+*         includes radiative capture energy. Neutrino energy is
+*         never included.
+* RRD     sum of radioactive decay constants in 10**-8/s.
+* AWR     mass of the nuclides in unit of neutron mass.
+* FUELDN  fuel initial density and mass.
+* SIG     initial/final microscopic depletion reaction rates for nuclide
+*         I in mixture IBM:
+*         SIG(I,1,IBM,:) fission reaction rate;
+*         SIG(I,2,IBM,:) gamma reaction rate;
+*         SIG(I,3,IBM,:) N2N reaction rate;
+*         cont...;
+*         SIG(I,NREAC,IBM,:) neutron-induced energy released;
+*         SIG(I,NREAC+1,IBM,:) decay energy released (10**-8 MeV/s).
+* VPH     initial/final integrated flux in fuel.
+* VPHV    initial/final integrated flux in each mixture.
+* MIXPWR  flags for mixtures to include in power normalization.
+* VTOTD   total fuel volume.
+* IEVOLB  flag making an isotope non-depleting:
+*         =0 the isotope is depleting;
+*         =1 to force an isotope to be non-depleting;
+*         =2 to force an isotope to be depleting;
+*         =3 to force an isotope to be at saturation
+* KFISS   position in chain of the fissile isotopes.
+* KPF     position in chain of the direct fission products (fission
+*         fragments).
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER IMPX,INR,IGLOB,NBMIX,NBISO,NCOMB,ISONAM(3,NBISO),
+     1 MILVO(NCOMB),JM(NBMIX,NVAR+NSUPS),NVAR,NDFP,NSUPS,NREAC,NPAR,
+     2 NFISS,ITYPE,IDIRAC,KPAR(NPAR,NVAR),IDR(NREAC,NVAR+NSUPS),
+     3 MIXPWR(NBMIX),IEVOLB(NVAR,NBMIX),KFISS(NFISS,NBMIX),
+     4 KPF(NDFP,NBMIX)
+      REAL YDPL(NVAR+1,2,NCOMB),VX(NBMIX),XT(2),EPS1,EPS2,EXPMAX,H1,FIT,
+     1 DELTA(3),ENERG(NBMIX),BPAR(NPAR,NVAR),YIELD(NFISS,NDFP,NBMIX),
+     2 RER(NREAC,NVAR+NSUPS),RRD(NVAR+NSUPS),AWR(NVAR),FUELDN(3),
+     3 SIG(NVAR+1,NREAC+1,NBMIX,2),VPH(2),VPHV(NBMIX,2)
+      DOUBLE PRECISION VTOTD
+*----
+*  LOCAL VARIABLES
+*----
+      CHARACTER TEXT8*8,HSMG*131
+      DOUBLE PRECISION GAR,GARD,XDRCST,EVJ,FITD,PHI2
+      LOGICAL LCOOL,LSIMPL
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: MU1,IMA,LP,CHAIN
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(MU1(NVAR+1),IMA(NVAR+1),LP(NVAR))
+*----
+*  CHECK IF ONLY THE HEAVY ISOTOPES ARE PRODUCING ENERGY. IN THIS CASE,
+*  SOME SIMPLIFICATIONS ARE POSSIBLE
+*----
+      LSIMPL=.TRUE.
+      DO 10 IS=1,NVAR
+      LSIMPL=LSIMPL.AND.(RER(3,IS).EQ.0.0)
+   10 CONTINUE
+*
+      EVJ=XDRCST('eV','J')*1.0E22
+      LCOOL=(INR.EQ.0)
+      IF(LCOOL) GO TO 410
+      IF(IMPX.GT.1) WRITE (6,640) VPH(1)/VTOTD,VPH(2)/VTOTD
+*----
+*  CONVERGE ON A FIXED FINAL POWER. SOLVE THE DEPLETION CHAIN WITHOUT
+*  THE FISSION PRODUCTS.
+*----
+      IF((INR.GE.2).AND.(EPS2.LT.10.0)) THEN
+         ITER=0
+  250    ITER=ITER+1
+         IF(ITER.GT.20) CALL XABORT('EVOBLD: UNABLE TO CONVERGE.')
+         DO 330 ICMB=1,NCOMB
+*        DETERMINE THE ROW AND COLUMN PROFILE OF THE ADPL MATRIX AND
+*        COMPUTE NVAR2, THE NUMBER OF DEPLETING NUCLIDES IN REGION ICMB.
+         IBM=MILVO(ICMB)
+         IF(IBM.EQ.0) GO TO 330
+         IF(IMPX.GT.3) WRITE(6,'(/34H EVOBLD: PROCESS DEPLETING MIXTURE,
+     1   I5,15H (REAL MIXTURE=,I5,2H).)') ICMB,IBM
+         NVAR2=0
+         NSUPL2=0
+         LP(:NVAR)=0
+         DO 270 IS=1,NVAR
+         KDRI=IDR(2,IS)
+         IF(KDRI.EQ.0) GO TO 270
+         IF((MOD(KDRI,100).NE.3).AND.(MOD(KDRI,100).NE.4)) GO TO 270
+         IF(JM(IBM,IS).GT.0) THEN
+            NVAR2=NVAR2+1
+            LP(IS)=NVAR2
+         ENDIF
+  270    CONTINUE
+         NSUPL2=NVAR2
+         DO 280 IS=1,NVAR
+         IF(LSIMPL.AND.(AWR(IS).LE.210.0)) GO TO 280
+         IF((JM(IBM,IS).GT.0).AND.(LP(IS).EQ.0)) THEN
+            NVAR2=NVAR2+1
+            LP(IS)=NVAR2
+         ENDIF
+  280    CONTINUE
+         IF(NVAR2.EQ.0) GO TO 330
+*        CHECK IF ONLY THE HEAVY ISOTOPES ARE PRODUCING ENERGY. IN
+*        THIS CASE, IT IS POSSIBLE TO AVOID THE SOLUTION FOR FISSION
+*        PRODUCTS.
+         NSUPF2=NVAR2-NSUPL2
+         IF(LSIMPL) NSUPF2=0
+         CALL EVOMU1(IMPX,NVAR,NREAC,LP,XT,LCOOL,NPAR,KPAR,RRD,
+     1   SIG(1,1,IBM,1),SIG(1,1,IBM,2),EXPMAX,IEVOLB(1,IBM),MU1,
+     2   IMA,MAXA)
+         MU1(NVAR2+1)=IMA(NVAR2)+NVAR2+1
+         IMA(NVAR2+1)=IMA(NVAR2)+NVAR2+1
+         MAXA=MAXA+10*(NVAR2+1)
+         NFISS2=0
+         DO 300 I=1,NFISS
+         IF(KFISS(I,IBM).EQ.0) GO TO 300
+         IF(LP(KFISS(I,IBM)).GT.0) NFISS2=NFISS2+1
+  300    CONTINUE
+         ALLOCATE(CHAIN(2*(NVAR2+1)))
+         DO 310 IS=1,NVAR
+         IF(LP(IS).GT.0) THEN
+            K=JM(IBM,IS)
+            CHAIN((LP(IS)-1)*2+1)=ISONAM(1,K)
+            CHAIN((LP(IS)-1)*2+2)=ISONAM(2,K)
+         ENDIF
+  310    CONTINUE
+         TEXT8='*POWER*'
+         READ(TEXT8,'(2A4)') (CHAIN(NVAR2*2+I0),I0=1,2)
+         CALL EVOSOL(IMPX,LCOOL,NVAR,NREAC,NDFP,NPAR,NFISS,XT,EPS1,
+     1   EXPMAX,H1,ITYPE,IDIRAC,RRD,KPAR,BPAR,KFISS(1,IBM),KPF(1,IBM),
+     2   YIELD(1,1,IBM),LP,IEVOLB(1,IBM),SIG(1,1,IBM,1),SIG(1,1,IBM,2),
+     3   NVAR2,NFISS2,NSUPF2,MU1,IMA,MAXA,YDPL(1,1,ICMB),CHAIN)
+*
+         DEALLOCATE(CHAIN)
+  330    CONTINUE
+         GAR=0.0D0
+         GARD=0.0D0
+         DO 360 IS=1,NVAR
+         DO 350 ICMB=1,NCOMB
+         IBM=MILVO(ICMB)
+         IF(IBM.EQ.0) GO TO 350
+         IF(MIXPWR(IBM).EQ.1) THEN
+            GAR=GAR+VX(IBM)*YDPL(IS,2,ICMB)*SIG(IS,NREAC,IBM,2)
+            GARD=GARD+VX(IBM)*YDPL(IS,2,ICMB)*SIG(IS,NREAC+1,IBM,2)
+         ENDIF
+  350    CONTINUE
+  360    CONTINUE
+         IF((IGLOB.EQ.1).OR.(INR.EQ.3)) THEN
+            DO 370 IBM=1,NBMIX
+            IF(MIXPWR(IBM).EQ.1) THEN
+               GAR=GAR+SIG(NVAR+1,NREAC,IBM,2)
+               GAR=GAR+SIG(NVAR+1,NREAC+1,IBM,2)
+            ENDIF
+  370       CONTINUE
+         ELSE
+            DO 380 ICMB=1,NCOMB
+            IBM=MILVO(ICMB)
+            IF(IBM.EQ.0) GO TO 380
+            IF(MIXPWR(IBM).EQ.1) THEN
+               GAR=GAR+SIG(NVAR+1,NREAC,IBM,2)
+               GAR=GAR+SIG(NVAR+1,NREAC+1,IBM,2)
+            ENDIF
+  380       CONTINUE
+         ENDIF
+         IF(GAR.EQ.0.0D0) CALL XABORT('EVOBLD: UNABLE TO NORMALIZE.')
+         IF(INR.EQ.3) THEN
+*           FITD IS THE DECAY POWER IN WATT PER CUBIC CENTIMETER.
+            FITD=(EVJ*GARD*FUELDN(3))/(FUELDN(1)*VTOTD)
+            IF(FITD.GT.FIT) THEN
+               WRITE(HSMG,'(35HEVOBLD: NEGATIVE FIT(1) FIT(DECAY)=,1P,
+     1         E11.4,12H FIT(INPUT)=,E11.4,1H.)') FITD,FIT
+               CALL XABORT(HSMG)
+            ENDIF
+            PHI2=(FIT-FITD)*FUELDN(1)*VPH(2)/(EVJ*GAR*FUELDN(3))
+         ELSE
+*           FITD IS THE DECAY POWER IN WATT PER GRAM.
+            FITD=(EVJ*GARD)/(FUELDN(1)*VTOTD)
+            IF(FITD.GT.FIT) THEN
+               WRITE(HSMG,'(35HEVOBLD: NEGATIVE FIT(2) FIT(DECAY)=,1P,
+     1         E11.4,12H FIT(INPUT)=,E11.4,1H.)') FITD,FIT
+               CALL XABORT(HSMG)
+            ENDIF
+            PHI2=(FIT-FITD)*FUELDN(1)*VPH(2)/(EVJ*GAR)
+         ENDIF
+         ERROR=REAL(ABS(PHI2-VPH(2)/VTOTD)/ABS(PHI2))
+         DO 400 IBM=1,NBMIX
+         VPHV(IBM,2)=VPHV(IBM,2)*REAL(PHI2*VTOTD)/VPH(2)
+         DO 395 IQ=1,NREAC
+         DO 390 IS=1,NVAR+1
+         SIG(IS,IQ,IBM,2)=SIG(IS,IQ,IBM,2)*REAL(PHI2*VTOTD)/VPH(2)
+  390    CONTINUE
+  395    CONTINUE
+  400    CONTINUE
+         VPH(2)=REAL(PHI2*VTOTD)
+         IF(IMPX.GT.3) THEN
+            WRITE (6,650) ITER,ERROR,VPH(1)/VTOTD,VPH(2)/VTOTD
+         ENDIF
+         IF(ERROR.LT.EPS2) THEN
+            IF(IMPX.GT.-1) WRITE(6,'(/29H EVOBLD: POWER CONVERGENCE IN,
+     1      I3,19H ITERATIONS. ERROR=,1P,E9.2,1H.)') ITER,ERROR
+            GO TO 410
+         ELSE
+            GO TO 250
+         ENDIF
+      ENDIF
+*----
+*  SOLVE THE COMPLETE DEPLETION CHAIN, INCLUDING FISSION PRODUCTS
+*----
+  410 DELTA(1)=0.0
+      ENERG(:NBMIX)=0.0
+      DO 500 ICMB=1,NCOMB
+*     DETERMINE THE ROW AND COLUMN PROFILE OF THE ADPL MATRIX AND
+*     COMPUTE NVAR2, THE NUMBER OF DEPLETING NUCLIDES IN REGION ICMB.
+      IBM=MILVO(ICMB)
+      IF(IBM.EQ.0) GO TO 500
+      IF(IMPX.GT.3) WRITE(6,'(/34H EVOBLD: PROCESS DEPLETING MIXTURE,
+     1 I5,15H (REAL MIXTURE=,I5,2H).)') ICMB,IBM
+      NVAR2=0
+      NSUPL2=0
+      LP(:NVAR)=0
+      DO 440 IS=1,NVAR
+      KDRI=IDR(2,IS)
+      IF(KDRI.EQ.0) GO TO 440
+      IF((MOD(KDRI,100).NE.3).AND.(MOD(KDRI,100).NE.4)) GO TO 440
+      IF(JM(IBM,IS).GT.0) THEN
+         NVAR2=NVAR2+1
+         LP(IS)=NVAR2
+      ENDIF
+  440 CONTINUE
+      NSUPL2=NVAR2
+      DO 450 IS=1,NVAR
+      IF((JM(IBM,IS).GT.0).AND.(LP(IS).EQ.0)) THEN
+         NVAR2=NVAR2+1
+         LP(IS)=NVAR2
+      ENDIF
+  450 CONTINUE
+      CALL EVOMU1(IMPX,NVAR,NREAC,LP,XT,LCOOL,NPAR,KPAR,RRD,
+     1 SIG(1,1,IBM,1),SIG(1,1,IBM,2),EXPMAX,IEVOLB(1,IBM),MU1,
+     2 IMA,MAXA)
+      MU1(NVAR2+1)=IMA(NVAR2)+NVAR2+1
+      IMA(NVAR2+1)=IMA(NVAR2)+NVAR2+1
+      MAXA=MAXA+10*(NVAR2+1)
+      NFISS2=0
+      DO 460 I=1,NFISS
+      IF(KFISS(I,IBM).EQ.0) GO TO 460
+      IF(LP(KFISS(I,IBM)).GT.0) NFISS2=NFISS2+1
+  460 CONTINUE
+      NSUPF2=NVAR2-NSUPL2
+      ALLOCATE(CHAIN(2*(NVAR2+1)))
+      DO 470 IS=1,NVAR
+      IF(LP(IS).GT.0) THEN
+         K=JM(IBM,IS)
+         CHAIN((LP(IS)-1)*2+1)=ISONAM(1,K)
+         CHAIN((LP(IS)-1)*2+2)=ISONAM(2,K)
+      ENDIF
+  470 CONTINUE
+      TEXT8='*POWER*'
+      READ(TEXT8,'(2A4)') (CHAIN(NVAR2*2+I0),I0=1,2)
+      CALL EVOSOL(IMPX,LCOOL,NVAR,NREAC,NDFP,NPAR,NFISS,XT,EPS1,
+     1 EXPMAX,H1,ITYPE,IDIRAC,RRD,KPAR,BPAR,KFISS(1,IBM),KPF(1,IBM),
+     2 YIELD(1,1,IBM),LP,IEVOLB(1,IBM),SIG(1,1,IBM,1),SIG(1,1,IBM,2),
+     3 NVAR2,NFISS2,NSUPF2,MU1,IMA,MAXA,YDPL(1,1,ICMB),CHAIN)
+*
+      DEALLOCATE(CHAIN)
+      IF(MIXPWR(IBM).EQ.1) THEN
+         DELTA(1)=DELTA(1)+YDPL(NVAR+1,2,ICMB)*VX(IBM)
+      ENDIF
+  500 CONTINUE
+*----
+*  BURNUP CALCULATION. TAKE THE CONTRIBUTION OBTAINED FROM THE ODE
+*  SOLVER AND ADD THE CONTRIBUTION FROM THE NON-DEPLETING ISOTOPES
+*  PRODUCING ENERGY
+*----
+      DO 510 IBM=1,NBMIX
+      IF(MIXPWR(IBM).EQ.1) THEN
+         GAR=0.5D0*(SIG(NVAR+1,NREAC,IBM,1)+SIG(NVAR+1,NREAC,IBM,2)
+     1             +SIG(NVAR+1,NREAC+1,IBM,1)+SIG(NVAR+1,NREAC+1,IBM,2))
+         ENERG(IBM)=ENERG(IBM)+REAL(GAR*(XT(2)-XT(1))*EVJ)
+      ENDIF
+  510 CONTINUE
+      DO 516 ICMB=1,NCOMB
+      IBM=MILVO(ICMB)
+      IF(IBM.EQ.0) GO TO 516
+      IF(MIXPWR(IBM).EQ.1) THEN
+         DO 515 IS=1,NVAR
+         GAR=0.5D0*(YDPL(IS,1,ICMB)*SIG(IS,NREAC,IBM,1)
+     1             +YDPL(IS,2,ICMB)*SIG(IS,NREAC,IBM,2)
+     1             +YDPL(IS,1,ICMB)*SIG(IS,NREAC+1,IBM,1)
+     1             +YDPL(IS,2,ICMB)*SIG(IS,NREAC+1,IBM,2))
+         ENERG(IBM)=ENERG(IBM)+REAL(GAR*(XT(2)-XT(1))*VX(IBM)*EVJ)
+  515    CONTINUE
+      ENDIF
+  516 CONTINUE
+      DELTA(3)=0.0
+      IF(IGLOB.LE.0) THEN
+         DO 520 ICMB=1,NCOMB
+         IBM=MILVO(ICMB)
+         IF(IBM.EQ.0) GO TO 520
+         IF(MIXPWR(IBM).EQ.1) THEN
+            GAR=0.5D0*(SIG(NVAR+1,NREAC,IBM,1)+SIG(NVAR+1,NREAC,IBM,2)
+     1                +SIG(NVAR+1,NREAC+1,IBM,1)
+     2                +SIG(NVAR+1,NREAC+1,IBM,2))
+            DELTA(1)=DELTA(1)+REAL(GAR)*(XT(2)-XT(1))
+            DELTA(3)=DELTA(3)+REAL(GAR)*(XT(2)-XT(1))
+         ENDIF
+  520    CONTINUE
+      ELSE IF(IGLOB.EQ.1) THEN
+         DO 530 IBM=1,NBMIX
+         IF(MIXPWR(IBM).EQ.1) THEN
+            GAR=0.5D0*(SIG(NVAR+1,NREAC,IBM,1)+SIG(NVAR+1,NREAC,IBM,2)
+     1                +SIG(NVAR+1,NREAC+1,IBM,1)
+     2                +SIG(NVAR+1,NREAC+1,IBM,2))
+            DELTA(1)=DELTA(1)+REAL(GAR)*(XT(2)-XT(1))
+            DELTA(3)=DELTA(3)+REAL(GAR)*(XT(2)-XT(1))
+         ENDIF
+  530    CONTINUE
+      ENDIF
+      DO 545 ICMB=1,NCOMB
+      IBM=MILVO(ICMB)
+      IF(IBM.EQ.0) GO TO 545
+      IF(MIXPWR(IBM).EQ.1) THEN
+         DO 540 IS=1,NVAR
+         DELTA(3)=DELTA(3)+0.5*(YDPL(IS,1,ICMB)*SIG(IS,NREAC,IBM,1)
+     1                 +YDPL(IS,2,ICMB)*SIG(IS,NREAC,IBM,2))*VX(IBM)
+     2                 *(XT(2)-XT(1))
+         DELTA(3)=DELTA(3)+0.5*(YDPL(IS,1,ICMB)*SIG(IS,NREAC+1,IBM,1)
+     1                 +YDPL(IS,2,ICMB)*SIG(IS,NREAC+1,IBM,2))*VX(IBM)
+     2                 *(XT(2)-XT(1))
+  540    CONTINUE
+      ENDIF
+  545 CONTINUE
+      IF(FUELDN(2) .EQ. 0.0) THEN
+        DELTA(1)=0.0
+        DELTA(3)=0.0
+      ELSE
+        DELTA(1)=DELTA(1)*REAL(EVJ)/FUELDN(2)
+        DELTA(3)=DELTA(3)*REAL(EVJ)/FUELDN(2)
+      ENDIF
+      DELTA(2)=0.5*(VPH(1)+VPH(2))*(XT(2)-XT(1))/REAL(VTOTD)
+      IF((.NOT.LCOOL).AND.(IMPX.GT.0)) THEN
+         IF(DELTA(1) .EQ. 0.0) THEN
+           WRITE (6,661) DELTA(2)
+         ELSE
+           WRITE (6,660) DELTA(1),DELTA(1)/8.64E-4,DELTA(2)
+           WRITE (6,665) DELTA(3),DELTA(3)/8.64E-4
+         ENDIF
+      ENDIF
+*----
+*  PRINT THE BEGINNING- AND END-OF-STAGE NORMALIZATION POWERS
+*----
+      IF(IMPX.GT.-1) THEN
+         DO 580 IP=1,2
+         DELTA1=0.0
+         DELTA2=0.0
+         IF(IGLOB.LE.0) THEN
+            DO 550 ICMB=1,NCOMB
+            IBM=MILVO(ICMB)
+            IF(IBM.EQ.0) GO TO 550
+            IF(MIXPWR(IBM).EQ.1) THEN
+               DELTA1=DELTA1+SIG(NVAR+1,NREAC,IBM,IP)
+               DELTA2=DELTA2+SIG(NVAR+1,NREAC+1,IBM,IP)
+            ENDIF
+  550       CONTINUE
+         ELSE IF(IGLOB.EQ.1) THEN
+            DO 560 IBM=1,NBMIX
+            IF(MIXPWR(IBM).EQ.1) THEN
+               DELTA1=DELTA1+SIG(NVAR+1,NREAC,IBM,IP)
+               DELTA2=DELTA2+SIG(NVAR+1,NREAC+1,IBM,IP)
+            ENDIF
+  560       CONTINUE
+         ENDIF
+         DO 575 ICMB=1,NCOMB
+         IBM=MILVO(ICMB)
+         IF(IBM.EQ.0) GO TO 575
+         IF(MIXPWR(IBM).EQ.1) THEN
+           DO 570 IS=1,NVAR
+           DELTA1=DELTA1+YDPL(IS,IP,ICMB)*SIG(IS,NREAC,IBM,IP)*VX(IBM)
+           DELTA2=DELTA2+YDPL(IS,IP,ICMB)*SIG(IS,NREAC+1,IBM,IP)*VX(IBM)
+  570      CONTINUE
+         ENDIF
+  575    CONTINUE
+         IF(FUELDN(2) .EQ. 0.0) THEN
+           DELTA1=0.0
+           DELTA2=0.0
+         ELSE
+           DELTA1=DELTA1*REAL(EVJ)/FUELDN(2)
+           DELTA2=DELTA2*REAL(EVJ)/FUELDN(2)
+           IF(IP.EQ.1) THEN
+              WRITE(6,680) 'BEGINNING-OF-STAGE',DELTA1
+              WRITE(6,690) 'BEGINNING-OF-STAGE',DELTA2
+           ELSE IF(IP.EQ.2) THEN
+              WRITE(6,680) 'END-OF-STAGE',DELTA1
+              WRITE(6,690) 'END-OF-STAGE',DELTA2
+           ENDIF
+         ENDIF
+  580    CONTINUE
+         WRITE(6,'(/48H NOTE: POWER MAY EXIBITS VARIATIONS OUTSIDE THE ,
+     1   52HBEGINNING- AND END-OF-STAGE VALUES DURING THE STAGE.)')
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(LP,IMA,MU1)
+      RETURN
+*
+  640 FORMAT(/53H EVOBLD: STARTING VALUES OF INITIAL/FINAL AVERAGED FL,
+     1 11HUX IN FUEL=,1P,2E12.4,15H E+13 N/CM**2/S)
+  650 FORMAT(/37H EVOBLD: ITERATION ON FINAL POWER NB.,I3,3X,6HERROR=,
+     1 1P,E12.4,3X,36HINITIAL/FINAL AVERAGED FLUX IN FUEL=,2E12.4,
+     2 15H E+13 N/CM**2/S)
+  660 FORMAT(/' EVOBLD: ',
+     1 'FUEL BURNUP INCREMENT DURING THIS STAGE =',1P,
+     1 E12.4,' MW*S**8/TONNE (',E12.4,' MW*DAY/TONNE).'/9X,
+     2 'NEUTRON EXPOSURE (FLUENCE) INCREMENT =',E12.4,' N/KB.')
+  661 FORMAT(/' EVOBLD: ',
+     1 'NEUTRON EXPOSURE (FLUENCE) INCREMENT DURING THIS STAGE =',1P,
+     2 E12.4,' N/KB.')
+  665 FORMAT(' EVOBLD: ',
+     1 'TARGET FUEL BURNUP INCREMENT DURING THIS STAGE =',1P,
+     2 E12.4,' MW*S**8/TONNE (',E12.4,' MW*DAY/TONNE).')
+  680 FORMAT(/34H EVOBLD: NEUTRON-INDUCED POWER AT ,A,2H =,1P,E12.4,
+     > 10H MW/TONNE.)
+  690 FORMAT(/24H EVOBLD: DECAY POWER AT ,A,2H =,1P,E12.4,10H MW/TONNE.)
+      END