summaryrefslogtreecommitdiff
path: root/Dragon/src/EVODRV.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/EVODRV.f')
-rw-r--r--Dragon/src/EVODRV.f1034
1 files changed, 1034 insertions, 0 deletions
diff --git a/Dragon/src/EVODRV.f b/Dragon/src/EVODRV.f
new file mode 100644
index 0000000..5015c51
--- /dev/null
+++ b/Dragon/src/EVODRV.f
@@ -0,0 +1,1034 @@
+*DECK EVODRV
+ SUBROUTINE EVODRV(IPDEPL,IPLIB,INDREC,IMPX,NBISO,NGROUP,NBMIX,
+ 1 ISONAM,ISONRF,MIX,DEN,IEVOL,ISTYP,VX,NDEPL,NSUPS,NREAC,NCOMB,
+ 2 EPS1,EPS2,EXPMAX,H1,ITYPE,INR,IEXTR,IGLOB,ISAT,IDIRAC,ITIXS,
+ 3 IFLMAC,IYLMIX,FIT,ISAVE,ISET,IDEPL,XTI,XTF,XT,LMACRO,FLUMIX,
+ 4 IPICK,MIXBRN,MIXPWR)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Isotopic depletion calculation main driver.
+*
+*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
+* IPDEPL pointer to the depletion history (L_BURNUP signature).
+* IPLIB pointer to the lattice microscopic cross section library
+* (L_LIBRARY signature).
+* INDREC beginning of depletion flag (=1: beginning; =2: else).
+* IMPX print flag (equal to zero for no print).
+* NBISO number of isotopes present in the calculation domain.
+* NGROUP number of energy groups.
+* NBMIX number of mixtures.
+* ISONAM alias name of isotopes.
+* ISONRF library name of isotopes.
+* MIX mix number of each isotope (can be zero).
+* DEN density of each isotope.
+* IEVOL non-depleting flag (=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).
+* ISTYP isotope type (=1 not fissile nor fission product; =2: fissile;
+* =3: fission product).
+* VX volume occupied by each mixture.
+* NDEPL number of isotopes in the depletion chain.
+* NSUPS number of non-depleting isotopes producing energy.
+* NREAC maximum number of depletion reactions.
+* NCOMB number of depleting mixtures.
+* 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.
+* INR type of flux normalization:
+* =0: out-of-core depletion;
+* =1: constant flux depletion;
+* =2: constant fuel power depletion;
+* =3: constant assembly power depletion.
+* IEXTR flux extrapolation flag (=0: no extrapolation; =1: perform
+* linear extrapolation; =2: perform parabolic extrapolation).
+* 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.
+* ISAT initial saturation flag (=1 to save initial saturated number
+* densities).
+* IDIRAC saturation model flag (=1 to use Dirac function contributions
+* in the saturating nuclide number densities).
+* ITIXS flag for time-dependent cross sections (=0/1: on/off).
+* IFLMAC 0/1/2 flag to recover fluxes from L_FLUX/L_MACROLIB/L_POWER.
+* IYLMIX 0/1 flag to recover fission yield data from DEPL-CHAIN/PYIELD
+* data.
+* 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.
+* ISAVE save flag:
+* =-1: do not save the last flux calculation in the depletion
+* table;
+* .GE.0 save the last flux calculation in the depletion
+* table at time XTI.
+* ISET set flag:
+* =-1: do not set the number densities to a selected time;
+* .GE.0 set the number densities to time XTF of the depletion
+* table.
+* IDEPL depletion flag:
+* =0: do not perform a depletion calculation
+* =1: perform a depletion calculation.
+* XTI initial save time (save the last flux calculation in the
+* depletion table at time XTI).
+* XTF final set time (recover the number densities from the
+* depletion table at time XTF and modify the internal library).
+* XT time variable (independent variable) for the depletion
+* calculation.
+* XT(1) initial time
+* XT(2) final time
+* LMACRO macrolib building flag (=.true. to compute the embedded
+* macrolib).
+* FLUMIX average fluxes in mixtures.
+* IPICK burnup recovery flag:
+* =0: do not recover the burnup in a CLE-2000 variable
+* =1: recover the burnup in a CLE-2000 variable.
+* MIXBRN flags for mixtures to burn.
+* MIXPWR flags for mixtures to include in power normalization.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ TYPE(C_PTR) IPDEPL,IPLIB
+ INTEGER INDREC,IMPX,NBISO,NGROUP,NBMIX,ISONAM(3,NBISO),
+ 1 ISONRF(3,NBISO),MIX(NBISO),IEVOL(NBISO),ISTYP(NBISO),NDEPL,
+ 2 NSUPS,NREAC,NCOMB,ITYPE,INR,IEXTR,IGLOB,ISAT,IDIRAC,ITIXS,
+ 3 IFLMAC,IYLMIX,ISAVE,ISET,IDEPL,IPICK,MIXBRN(NBMIX),
+ 4 MIXPWR(NBMIX)
+ REAL DEN(NBISO),VX(NBMIX),EPS1,EPS2,EXPMAX,H1,FIT,XTI,XTF,
+ 1 XT(2),FLUMIX(NGROUP,NBMIX)
+ LOGICAL LMACRO
+*----
+* LOCAL VARIABLES
+*----
+ PARAMETER (IUNOUT=6,NSTATE=40)
+ TYPE(C_PTR) KPLIB
+ CHARACTER TEXT12*12,HSMG*131
+ LOGICAL LCOOL
+ INTEGER IDIM(NSTATE),IPAR(NSTATE)
+ REAL DELTA(3),RPAR(5),BRNWIR(2),TMPDAY(3),VPH(2),
+ 1 FUELDN(3),DELTAT(2,2),TIMEP(2,3)
+ DOUBLE PRECISION T(3),WEI,DPD,XDRCST,AVCON,VTOTD,VPHINI,DBLLIR
+*----
+* ALLOCATABLE ARRAYS
+*----
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: MILVO,ISOCMB,NFISS2,
+ 1 NDFP2,HREAC,IPIFI,IZAE
+ INTEGER, ALLOCATABLE, DIMENSION(:,:) :: JM,INADPL,IEVOLB,KFISS,
+ 1 KPAR,IDR,KPF
+ REAL, ALLOCATABLE, DIMENSION(:) :: ENERG,RRD,AWR,PYIELD,TIMES
+ REAL, ALLOCATABLE, DIMENSION(:,:) :: BPAR,RER,VPHV
+ REAL, ALLOCATABLE, DIMENSION(:,:,:) :: YDPL,YIELD,YIELD2
+ REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: SIG
+ LOGICAL, ALLOCATABLE, DIMENSION(:) :: MASK,MASKL
+ TYPE(C_PTR), ALLOCATABLE, DIMENSION(:) :: IPISO
+*----
+* SCRATCH STORAGE ALLOCATION
+* MILVO mixture index corresponding to each depleting mixture.
+*----
+ ALLOCATE(JM(NBMIX,NDEPL),MILVO(NCOMB),ISOCMB(NBISO),
+ 1 INADPL(3,NDEPL),IEVOLB(NDEPL,NBMIX))
+ ALLOCATE(SIG(NDEPL-NSUPS+1,NREAC+1,NBMIX,1-IEXTR:2),
+ 1 VPHV(NBMIX,1-IEXTR:2),ENERG(NBMIX),AWR(NDEPL),IZAE(NDEPL),
+ 2 YDPL(NDEPL-NSUPS+1,2,NCOMB))
+ ALLOCATE(MASK(NBMIX),MASKL(NGROUP))
+ ALLOCATE(IPISO(NBISO))
+*----
+* INITIALIZE DATA.
+*----
+ AVCON=XDRCST('Neutron mass','amu')/
+ 1 (1.0D-24*XDRCST('Avogadro','N/moles'))
+*
+ IF(INDREC.EQ.1) THEN
+* BEGINNING OF DEPLETION.
+ NTIM=0
+ ELSE IF(INDREC.EQ.2) THEN
+ CALL LCMLEN(IPDEPL,'DEPL-TIMES',NTIM,ITYLCM)
+ ENDIF
+ IF(NTIM+2.GE.10000) CALL XABORT('EVODRV: No more than 9999'//
+ >'burnup steps permitted.')
+ ALLOCATE(TIMES(NTIM+2))
+ TIMES(:NTIM+2)=0.0
+ IF(NTIM.EQ.0) THEN
+ CALL LCMPUT(IPDEPL,'DEPL-TIMES',1,2,TIMES)
+ ELSE
+ CALL LCMGET(IPDEPL,'DEPL-TIMES',TIMES)
+ ENDIF
+*----
+* RECOVER DEPLETION CHAIN INFO FROM LCM
+*----
+ NVAR=NDEPL-NSUPS
+ CALL LCMLEN(IPLIB,'DEPL-CHAIN',LENGT,ITYLCM)
+ IF (LENGT.EQ.0) CALL XABORT('EVODRV: DEPLETION CHAIN MISSING.')
+ CALL LCMSIX(IPLIB,'DEPL-CHAIN',1)
+ CALL LCMGET(IPLIB,'STATE-VECTOR',IDIM)
+ IF((NDEPL.NE.IDIM(1)).OR.(NSUPS.NE.IDIM(7)))
+ 1 CALL XABORT('EVODRV: INCONSISTENT NDEPL OR NSUPS.')
+ IF(NVAR.EQ.0) CALL XABORT('EVODRV: NO DEPLETING ISOTOPES')
+ NFISS=IDIM(2)
+ NDFP=IDIM(3)
+ NPAR=IDIM(9)
+ NBESP=MAX(1,IDIM(10))
+ ALLOCATE(KPAR(NDEPL,NPAR),HREAC(2*NREAC),IDR(NREAC,NDEPL))
+ ALLOCATE(BPAR(NDEPL,NPAR),YIELD2(NBESP,NFISS,NDFP),
+ 1 RER(NREAC,NDEPL),RRD(NDEPL))
+ CALL LCMGET(IPLIB,'ISOTOPESDEPL',INADPL)
+ IF(IMPX.GT.1) THEN
+ WRITE(IUNOUT,'(/38HEVODRV: DEPLETING ISOTOPES FROM CHAIN:)')
+ WRITE(IUNOUT,'(1X,3A4,2X,3A4,2X,3A4,2X,3A4,2X,3A4,2X,3A4,2X,
+ 1 3A4,2X,3A4,2X,3A4,2X,3A4,2X,3A4,2X,3A4)') INADPL(:3,:NVAR)
+ ENDIF
+ CALL LCMGET(IPLIB,'PRODUCE-REAC',KPAR)
+ CALL LCMGET(IPLIB,'PRODUCE-RATE',BPAR)
+ CALL LCMGET(IPLIB,'DEPLETE-IDEN',HREAC)
+ CALL LCMGET(IPLIB,'DEPLETE-REAC',IDR)
+ CALL LCMGET(IPLIB,'DEPLETE-ENER',RER)
+ CALL LCMGET(IPLIB,'DEPLETE-DECA',RRD)
+ CALL LCMGET(IPLIB,'CHARGEWEIGHT',IZAE)
+ IF(NFISS*NDFP.GT.0) CALL LCMGET(IPLIB,'FISSIONYIELD',YIELD2)
+ CALL LCMSIX(IPLIB,' ',2)
+*----
+* SET THE LCM MICROLIB ISOTOPEWISE DIRECTORIES.
+*----
+ CALL LIBIPS(IPLIB,NBISO,IPISO)
+*----
+* DETECT THE DEPLETING ISOTOPES AND MIXTURES IN THE MICROLIB.
+*----
+ ICOMB=0
+ JM(:NBMIX,:NDEPL)=0
+ IEVOLB(:NDEPL,:NBMIX)=1
+ ALLOCATE(NFISS2(NBMIX),NDFP2(NBMIX))
+ NFISS2(:NBMIX)=0
+ NDFP2(:NBMIX)=0
+ AWR(:NDEPL)=0.0
+ DO 30 ISOT=1,NBISO
+ IBM=MIX(ISOT)
+ IF(IBM.EQ.0) GO TO 30
+ KPLIB=IPISO(ISOT) ! set ISOT-th isotope
+ DO 20 INUCL=1,NDEPL
+ IF((ISONRF(1,ISOT).EQ.INADPL(1,INUCL)).AND.(ISONRF(2,ISOT).EQ.
+ 1 INADPL(2,INUCL))) THEN
+ IF(JM(IBM,INUCL).GT.0) GO TO 20
+ IF(C_ASSOCIATED(KPLIB)) THEN
+ CALL LCMLEN(KPLIB,'AWR',ILONG,ITYLCM)
+ IF(ILONG.EQ.1) CALL LCMGET(KPLIB,'AWR',AWR(INUCL))
+ ENDIF
+ IF(INUCL.GT.NVAR) THEN
+ IF(JM(IBM,INUCL).EQ.0) THEN
+ IEVOLB(INUCL,IBM)=1
+ JM(IBM,INUCL)=-ISOT
+ ENDIF
+ ELSE IF(IEVOL(ISOT).EQ.1) THEN
+ IF(JM(IBM,INUCL).EQ.0) THEN
+ IEVOLB(INUCL,IBM)=1
+ JM(IBM,INUCL)=-ISOT
+ ENDIF
+ IF(ISTYP(ISOT).EQ.2) THEN
+ NFISS2(IBM)=NFISS2(IBM)+1
+ JM(IBM,INUCL)=ISOT
+ ENDIF
+ ELSE
+ IF(ISTYP(ISOT).EQ.2) THEN
+ NFISS2(IBM)=NFISS2(IBM)+1
+ ELSE IF(ISTYP(ISOT).EQ.3) THEN
+ NDFP2(IBM)=NDFP2(IBM)+1
+ ENDIF
+ IEVOLB(INUCL,IBM)=IEVOL(ISOT)
+ JM(IBM,INUCL)=ISOT
+ IF(IEVOL(ISOT).EQ.1) GO TO 30
+ DO 10 J=1,ICOMB
+ IF(IBM.EQ.MILVO(J)) GO TO 30
+ 10 CONTINUE
+ ICOMB=ICOMB+1
+ MILVO(ICOMB)=IBM
+ ENDIF
+ GO TO 30
+ ENDIF
+ 20 CONTINUE
+ 30 CONTINUE
+ IF(ICOMB.NE.NCOMB) THEN
+ WRITE(HSMG,'(38HEVODRV: INVALID VALUE OF NCOMB (ICOMB=,I5,
+ 1 7H NCOMB=,I5,2H).)') ICOMB,NCOMB
+ CALL XABORT(HSMG)
+ ENDIF
+ DO 35 J=1,ICOMB
+ IBM=MILVO(J)
+ IF(MIXBRN(IBM).EQ.0) MILVO(J)=0
+ 35 CONTINUE
+ IF(IYLMIX.EQ.1) THEN
+ NFISS=MAXVAL(NFISS2)
+ NDFP=MAXVAL(NDFP2)
+ ENDIF
+ DEALLOCATE(NDFP2,NFISS2)
+ IF(IMPX.GT.0) WRITE(IUNOUT,500) NFISS,NDFP
+*----
+* SET KFISS, KPF AND YIELD
+*----
+ ALLOCATE(KFISS(NFISS,NBMIX),KPF(NDFP,NBMIX),
+ 1 YIELD(NFISS,NDFP,NBMIX))
+ KFISS(:NFISS,:NBMIX)=0
+ KPF(:NDFP,:NBMIX)=0
+ IF(IYLMIX.EQ.0) THEN
+* Use fission yield data from 'DEPL-CHAIN'
+ DO 40 IS=1,NVAR
+ KDRI=IDR(2,IS)/100
+ IF((KDRI.GT.0).AND.(MOD(IDR(2,IS),100).EQ.4)) THEN
+ IF(KDRI.GT.NFISS) CALL XABORT('EVODRV: INVALID NFISS.')
+ KFISS(KDRI,:NBMIX)=IS
+ ELSE IF((KDRI.GT.0).AND.(MOD(IDR(2,IS),100).EQ.5)) THEN
+ IF(KDRI.GT.NDFP) CALL XABORT('EVODRV: INVALID NDFP.')
+ KPF(KDRI,:NBMIX)=IS
+ ENDIF
+ 40 CONTINUE
+ DO IDFP=1,NDFP
+ DO IFISS=1,NFISS
+ YIELD(IFISS,IDFP,:NBMIX)=YIELD2(NBESP,IFISS,IDFP)
+ ENDDO
+ ENDDO
+ ELSE IF(IYLMIX.EQ.1) THEN
+* Use isotopic PIFI/PYIELD fission yield data
+ YIELD(:NFISS,:NDFP,:NBMIX)=0.0
+ DO 65 IBM=1,NBMIX
+ IFISS=0
+ IDFP=0
+ DO 50 IS=1,NVAR
+ ISOT=ABS(JM(IBM,IS))
+ IF(ISOT.EQ.0) GO TO 50
+ IF(ISTYP(ISOT).EQ.2) THEN
+ IFISS=IFISS+1
+ IF(IFISS.GT.NFISS) CALL XABORT('EVODRV: NFISS OVERFLOW.')
+ KFISS(IFISS,IBM)=IS
+ IF(IDR(2,IS).EQ.0) IDR(2,IS)=4
+ ELSE IF(ISTYP(ISOT).EQ.3) THEN
+ IDFP=IDFP+1
+ IF(IDFP.GT.NDFP) CALL XABORT('EVODRV: NDFP OVERFLOW.')
+ KPF(IDFP,IBM)=IS
+ IF(IDR(2,IS).EQ.0) IDR(2,IS)=5
+ ENDIF
+ 50 CONTINUE
+ DO 60 IS=1,NVAR
+ ISOT=JM(IBM,IS)
+ IF(ISOT.LE.0) GO TO 60
+ KPLIB=IPISO(ISOT) ! set ISOT-th isotope
+ IF(.NOT.C_ASSOCIATED(KPLIB)) THEN
+ WRITE(HSMG,'(17HEVODRV: ISOTOPE '',3A4,16H'' IS NOT AVAILAB,
+ > 22HLE IN THE MICROLIB(1).)') (ISONAM(I0,ISOT),I0=1,3)
+ CALL XABORT(HSMG)
+ ENDIF
+ CALL LCMLEN(KPLIB,'PIFI',NDFI,ITYLCM)
+ IF(NDFI.GT.0) THEN
+ ALLOCATE(IPIFI(NDFI),PYIELD(NDFI))
+ CALL LCMGET(KPLIB,'PIFI',IPIFI)
+ CALL LCMGET(KPLIB,'PYIELD',PYIELD)
+ IDFP=0
+ DO I=1,NDFP
+ IF(KPF(I,IBM).EQ.IS) THEN
+ IDFP=I
+ EXIT
+ ENDIF
+ ENDDO
+ IF(IDFP.EQ.0) THEN
+ WRITE(IUNOUT,510) 'FISSION PRODUCT',IS,(KPF(I,IBM),
+ 1 I=1,NDFP)
+ WRITE(HSMG,'(39HEVODRV: UNABLE TO FIND FP INDEX FOR ISO,
+ 1 5HTOPE ,3A4,5H (1).)')(ISONAM(I0,ISOT),I0=1,3)
+ CALL XABORT(HSMG)
+ ENDIF
+ DO 55 I=1,NDFI
+ IF(IPIFI(I).EQ.0) GO TO 55
+ DO JST=1,NVAR
+ IF(ABS(JM(IBM,JST)).EQ.IPIFI(I)) THEN
+ IFISS=0
+ DO J=1,NFISS
+ IF(KFISS(J,IBM).EQ.JST) THEN
+ IFISS=J
+ EXIT
+ ENDIF
+ ENDDO
+ IF(IFISS.EQ.0) THEN
+ WRITE(IUNOUT,510) 'FISSILE ISOTOPE',JST,
+ 1 (KFISS(J,IBM),J=1,NFISS)
+ CALL XABORT('EVODRV: UNABLE TO FIND FISSILE I'
+ 1 //'SOTOPE INDEX')
+ ENDIF
+ YIELD(IFISS,IDFP,IBM)=PYIELD(I)
+ GO TO 55
+ ENDIF
+ ENDDO
+ WRITE(HSMG,'(39HEVODRV: UNABLE TO FIND FP INDEX FOR ISO,
+ 1 5HTOPE ,3A4,5H (2).)')(ISONAM(I0,ISOT),I0=1,3)
+ CALL XABORT(HSMG)
+ 55 CONTINUE
+ DEALLOCATE(PYIELD,IPIFI)
+ ENDIF
+ 60 CONTINUE
+ 65 CONTINUE
+ ELSE
+ CALL XABORT('EVODRV: INVALID VALUE OF FLAG IYLMIX.')
+ ENDIF
+ DEALLOCATE(YIELD2)
+*----
+* COMPUTE THE INITIAL INTEGRATED FLUX
+*----
+ VPHINI=0.0D0
+ DO 85 IU=1,NGROUP
+ DO 80 ICMB=1,NCOMB
+ IBM=MILVO(ICMB)
+ IF(IBM.GT.0) VPHINI=VPHINI+FLUMIX(IU,IBM)*VX(IBM)
+ 80 CONTINUE
+ 85 CONTINUE
+*----
+* CHECK IF PERTURBATION XS OR STANDARD XS.
+*----
+ NXSPER=1
+ CALL LCMLEN(IPLIB,'TIMESPER',LENGTH,ITYLCM)
+ IF((LENGTH.GE.2).AND.(LENGTH.LE.6)) THEN
+ CALL LCMGET(IPLIB,'TIMESPER',TIMEP)
+ DELTAT(1,1)=TIMEP(1,1)
+ DELTAT(2,1)=TIMEP(2,1)
+ TMPREF=DELTAT(1,1)
+ NXSPER=2
+ IF(ITIXS.EQ.0) THEN
+ DO 90 IP=1,2
+ DELTAT(1,IP)=1.0
+ DELTAT(2,IP)=XT(IP)/8.64E-4-TMPREF
+ 90 CONTINUE
+ ELSE
+ XREF=XT(1)/8.64E-4
+ DO 100 IP=1,2
+ DELTAT(1,IP)=1.0
+ DELTAT(2,IP)=XREF-TMPREF
+ 100 CONTINUE
+ ENDIF
+ ELSE
+ DO 110 IP=1,2
+ DELTAT(1,IP)=1.0
+ DELTAT(2,IP)=0.0
+ 110 CONTINUE
+ ENDIF
+*----
+* COMPUTE AND SAVE THE INITIAL MASS OF HEAVY ELEMENTS
+*----
+ VTOTD=0.0D0
+ DO ICMB=1,NCOMB
+ IBM=MILVO(ICMB)
+ IF(IBM.EQ.0) CYCLE
+ IF(MIXPWR(IBM).GT.0) VTOTD=VTOTD+DBLE(VX(IBM))
+ ENDDO
+ IF(INDREC.EQ.1) THEN
+*
+* COMPUTE THE GLOBAL HEAVY-ELEMENT MASS FOR ISOTOPES IN MIXPWR.
+ FUELDN(1)=0.0
+ IF(IGLOB.LE.0) THEN
+ DO 120 ICMB=1,NCOMB
+ IBM=MILVO(ICMB)
+ IF(IBM.EQ.0) GO TO 120
+ IF(MIXPWR(IBM).GT.0) THEN
+ DO 115 IS=1,NBISO
+ IF((MIX(IS).NE.IBM).OR.(IEVOL(IS).EQ.1)) GO TO 115
+ KPLIB=IPISO(IS) ! set IS-th isotope
+ IF(C_ASSOCIATED(KPLIB)) THEN
+ AWRGAR=0.0
+ CALL LCMLEN(KPLIB,'AWR',ILONG,ITYLCM)
+ IF(ILONG.EQ.1) CALL LCMGET(KPLIB,'AWR',AWRGAR)
+ IF((AWRGAR.GT.210.0).OR.(ISTYP(IS).EQ.2)) THEN
+ FUELDN(1)=FUELDN(1)+AWRGAR*DEN(IS)*VX(IBM)
+ ENDIF
+ ENDIF
+ 115 CONTINUE
+ ENDIF
+ 120 CONTINUE
+ ELSE IF(IGLOB.EQ.1) THEN
+ DO 125 IBM=1,NBMIX
+ IF(MIXPWR(IBM).GT.0) THEN
+ DO 124 IS=1,NBISO
+ IF(MIX(IS).NE.IBM) GO TO 124
+ KPLIB=IPISO(IS) ! set IS-th isotope
+ IF(C_ASSOCIATED(KPLIB)) THEN
+ AWRGAR=0.0
+ CALL LCMLEN(KPLIB,'AWR',ILONG,ITYLCM)
+ IF(ILONG.EQ.1) CALL LCMGET(KPLIB,'AWR',AWRGAR)
+ IF((AWRGAR.GT.210.0).OR.(ISTYP(IS).EQ.2)) THEN
+ FUELDN(1)=FUELDN(1)+AWRGAR*DEN(IS)*VX(IBM)
+ ENDIF
+ ENDIF
+ 124 CONTINUE
+ ENDIF
+ 125 CONTINUE
+ ENDIF
+ IF(FUELDN(1).EQ.0) THEN
+ IF(INR.LE.1) THEN
+ FUELDN(1)=0.0
+ ELSE
+ CALL XABORT('EVODRV: Burnup at fixed power without '//
+ 1 'heavy fissile isotopes is forbidden')
+ ENDIF
+ ELSE
+ FUELDN(1)=FUELDN(1)*REAL(AVCON/VTOTD)
+ ENDIF
+ FUELDN(2)=FUELDN(1)*REAL(VTOTD)
+ VASSMB=0.0
+ DO 130 IBM=1,NBMIX
+ IF(MIXPWR(IBM).EQ.1) THEN
+ VASSMB=VASSMB+VX(IBM)
+ ENDIF
+ 130 CONTINUE
+ FUELDN(3)=FUELDN(2)/VASSMB
+ CALL LCMPUT(IPDEPL,'FUELDEN-INIT',3,2,FUELDN)
+*
+* COMPUTE THE HEAVY-ELEMENT MASS PER MIXTURE.
+ DO 150 IBM=1,NBMIX
+ ENERG(IBM)=0.0
+ DO 140 IS=1,NBISO
+ IF((MIX(IS).NE.IBM).OR.(IEVOL(IS).EQ.1)) GO TO 140
+ KPLIB=IPISO(IS) ! set IS-th isotope
+ IF(C_ASSOCIATED(KPLIB)) THEN
+ AWRGAR=0.0
+ CALL LCMLEN(KPLIB,'AWR',ILONG,ITYLCM)
+ IF(ILONG.EQ.1) CALL LCMGET(KPLIB,'AWR',AWRGAR)
+ IF((AWRGAR.GT.210.0).OR.(ISTYP(IS).EQ.2)) THEN
+ ENERG(IBM)=ENERG(IBM)+AWRGAR*DEN(IS)*VX(IBM)
+ ENDIF
+ ENDIF
+ 140 CONTINUE
+ ENERG(IBM)=ENERG(IBM)*REAL(AVCON)
+ 150 CONTINUE
+ CALL LCMPUT(IPDEPL,'FUELDEN-MIX',NBMIX,2,ENERG)
+*
+* COMPUTE THE TOTAL MASS PER MIXTURE.
+ ENERG(:NBMIX)=0.0
+ DO 170 IS=1,NBISO
+ IF(DEN(IS).EQ.0.0) GO TO 170
+ KPLIB=IPISO(IS) ! set IS-th isotope
+ AWRGAR=0.0
+ IF(C_ASSOCIATED(KPLIB)) THEN
+ CALL LCMLEN(KPLIB,'AWR',ILONG,ITYLCM)
+ IF(ILONG.EQ.1) CALL LCMGET(KPLIB,'AWR',AWRGAR)
+ ENDIF
+ DO 160 IBM=1,NBMIX
+ IF(MIX(IS).EQ.IBM) THEN
+ ENERG(IBM)=ENERG(IBM)+REAL(AWRGAR*DEN(IS)*VX(IBM)*AVCON)
+ ENDIF
+ 160 CONTINUE
+ 170 CONTINUE
+ CALL LCMPUT(IPDEPL,'WEIGHT-MIX',NBMIX,2,ENERG)
+ ELSE
+ CALL LCMGET(IPDEPL,'FUELDEN-INIT',FUELDN)
+ ENDIF
+ IF(IMPX.GT.0) THEN
+ WRITE (6,610) FUELDN(1),VTOTD,FUELDN(2),FUELDN(3)
+ ENDIF
+*----
+* SAVE THE LAST FLUX CALCULATION SET POINT IN THE DEPLETION TABLE.
+* CROSS-SECTION PERTURBATION IS ENABLED IF ISAVE=0 AND NXSPER=2.
+*----
+ IF(ISAVE.EQ.0) THEN
+ DO 200 IP=1,NXSPER
+ ITIM=0
+ IF(NTIM.GT.0) THEN
+ DO 180 I=1,NTIM
+ IF(ABS(TIMES(I)-XT(IP)).LE.1.0E-4*XT(IP)) ITIM=I
+ 180 CONTINUE
+ ENDIF
+ IF(ITIM.EQ.0) THEN
+ IF(NTIM.GT.0) THEN
+ IF(XT(IP).LT.TIMES(NTIM)) CALL XABORT('EVODRV: INVALID X1.')
+ ENDIF
+ NTIM=NTIM+1
+ TIMES(NTIM)=XT(IP)
+ CALL LCMPUT(IPDEPL,'DEPL-TIMES',NTIM,2,TIMES)
+ ITIM=NTIM
+ ENDIF
+ WRITE(TEXT12,'(8HDEPL-DAT,I4.4)') ITIM
+ IF(IMPX.GT.0) WRITE(IUNOUT,530) XT(IP),XT(IP)/8.64E-4,TEXT12
+ CALL LCMSIX(IPDEPL,TEXT12,1)
+ CALL LCMPUT(IPDEPL,'ISOTOPESDENS',NBISO,2,DEN)
+*----
+* COMPUTE, NORMALIZE AND SAVE THE MICROSCOPIC REACTION RATES.
+*----
+ CALL EVOSIG(IMPX,INR,IGLOB,NGROUP,NBMIX,NBISO,NCOMB,ISONAM,
+ 1 IPISO,DEN,FLUMIX,VX,MILVO,JM,NVAR,NSUPS,NREAC,HREAC,IDR,
+ 2 RER,RRD,FIT,AWR,IZAE,FUELDN,NXSPER,DELTAT(1,IP),MIXPWR,PFACT,
+ 3 SIG(1,1,1,IP),VPHV(1,IP))
+ NLENGT=(NVAR+1)*(NREAC+1)*NBMIX
+ CALL LCMPUT(IPDEPL,'MICRO-RATES',NLENGT,2,SIG(1,1,1,IP))
+ CALL LCMPUT(IPDEPL,'INT-FLUX',NBMIX,2,VPHV(1,IP))
+ VPHD=0.0
+ DO 190 ICMB=1,NCOMB
+ IBM=MILVO(ICMB)
+ IF(IBM.GT.0) VPHD=VPHD+VPHV(IBM,IP)
+ 190 CONTINUE
+ VPH(IP)=VPHD
+ IF(INR.NE.0) THEN
+ FNORM=VPH(IP)/REAL(VPHINI)
+ IF(INR.EQ.3) CALL LCMPUT(IPDEPL,'FORM-POWER',1,2,PFACT)
+ ELSE
+ FNORM=0.0
+ ENDIF
+ CALL LCMPUT(IPDEPL,'FLUX-NORM',1,2,FNORM)
+ IF((INDREC.EQ.1).AND.(IP.EQ.1)) THEN
+ BRNWIR(1)=0.0
+ BRNWIR(2)=0.0
+ CALL LCMPUT(IPDEPL,'BURNUP-IRRAD',2,2,BRNWIR)
+ ENDIF
+ CALL LCMSIX(IPDEPL,' ',2)
+ 200 CONTINUE
+ ELSE IF(ISAVE.EQ.1) THEN
+ IP=1
+ ITIM=0
+ IF(NTIM.GT.0) THEN
+ DO 210 I=1,NTIM
+ IF(ABS(TIMES(I)-XTI).LE.1.0E-4*XTI) ITIM=I
+ 210 CONTINUE
+ ENDIF
+ IF(ITIM.EQ.0) THEN
+ IF(NTIM.GT.0) THEN
+ IF(XTI.LT.TIMES(NTIM)) CALL XABORT('EVODRV: INVALID X1.')
+ ENDIF
+ NTIM=NTIM+1
+ TIMES(NTIM)=XTI
+ CALL LCMPUT(IPDEPL,'DEPL-TIMES',NTIM,2,TIMES)
+ ITIM=NTIM
+ ENDIF
+ WRITE(TEXT12,'(8HDEPL-DAT,I4.4)') ITIM
+ IF(IMPX.GT.0) WRITE(IUNOUT,530) XTI,XTI/8.64E-4,TEXT12
+ CALL LCMSIX(IPDEPL,TEXT12,1)
+ CALL LCMPUT(IPDEPL,'ISOTOPESDENS',NBISO,2,DEN)
+*----
+* COMPUTE, NORMALIZE AND SAVE THE MICROSCOPIC REACTION RATES.
+*----
+ CALL EVOSIG(IMPX,INR,IGLOB,NGROUP,NBMIX,NBISO,NCOMB,ISONAM,
+ 1 IPISO,DEN,FLUMIX,VX,MILVO,JM,NVAR,NSUPS,NREAC,HREAC,IDR,
+ 2 RER,RRD,FIT,AWR,IZAE,FUELDN,NXSPER,DELTAT(1,IP),MIXPWR,PFACT,
+ 3 SIG(1,1,1,IP),VPHV(1,IP))
+ NLENGT=(NVAR+1)*(NREAC+1)*NBMIX
+ CALL LCMPUT(IPDEPL,'MICRO-RATES',NLENGT,2,SIG(1,1,1,IP))
+ CALL LCMPUT(IPDEPL,'INT-FLUX',NBMIX,2,VPHV(1,IP))
+ VPHD=0.0
+ DO 220 ICMB=1,NCOMB
+ IBM=MILVO(ICMB)
+ IF(IBM.GT.0) VPHD=VPHD+VPHV(IBM,IP)
+ 220 CONTINUE
+ VPH(IP)=VPHD
+ IF(INR.NE.0) THEN
+ FNORM=VPH(IP)/REAL(VPHINI)
+ IF(INR.EQ.3) CALL LCMPUT(IPDEPL,'FORM-POWER',1,2,PFACT)
+ ELSE
+ FNORM=0.0
+ ENDIF
+ CALL LCMPUT(IPDEPL,'FLUX-NORM',1,2,FNORM)
+ IF((INDREC.EQ.1).AND.(IP.EQ.1)) THEN
+ BRNWIR(1)=0.0
+ BRNWIR(2)=0.0
+ CALL LCMPUT(IPDEPL,'BURNUP-IRRAD',2,2,BRNWIR)
+ ENDIF
+ CALL LCMSIX(IPDEPL,' ',2)
+ ENDIF
+*
+ IF(IDEPL.EQ.1) THEN
+*----
+* PERFORM A DEPLETION CALCULATION BETWEEN TIMES XT(1) AND XT(2).
+*----
+ IF(IMPX.GT.0) WRITE(IUNOUT,600) XT(1),XT(2)
+ LCOOL=.TRUE.
+ DO 300 IP=1,2
+ ITIM=0
+ DO 230 I=1,NTIM
+ IF(ABS(TIMES(I)-XT(IP)).LE.1.0E-4*XT(IP)) ITIM=I
+ 230 CONTINUE
+ IF(ITIM.GT.0) THEN
+ WRITE(TEXT12,'(8HDEPL-DAT,I4.4)') ITIM
+ IF(IMPX.GT.0) WRITE(IUNOUT,520) XT(IP),XT(IP)/8.64E-4,TEXT12
+ CALL LCMSIX(IPDEPL,TEXT12,1)
+ CALL LCMGET(IPDEPL,'ISOTOPESDENS',DEN)
+ IF(IP.EQ.1) CALL LCMGET(IPDEPL,'BURNUP-IRRAD',BRNWIR)
+ CALL LCMLEN(IPDEPL,'MICRO-RATES',LENGT,ITYLCM)
+ IF((LENGT.GT.0).OR.(IP.EQ.1)) THEN
+ CALL LCMGET(IPDEPL,'MICRO-RATES',SIG(1,1,1,IP))
+ CALL LCMGET(IPDEPL,'INT-FLUX',VPHV(1,IP))
+ ENDIF
+ IF((IP.EQ.1).AND.(INR.EQ.3)) THEN
+ CALL LCMGET(IPDEPL,'FORM-POWER',PFACT)
+ ENDIF
+ CALL LCMSIX(IPDEPL,' ',2)
+ ELSE
+ IF(IP.EQ.1) CALL XABORT('EVODRV: NO DEPLETION DATA STORED.')
+ IF((IEXTR.GE.1).AND.(NTIM.GE.2).AND.(INR.NE.0)) THEN
+* PERFORM MICRO REACTION RATE EXTRAPOLATION.
+ ITIM=0
+ DO 240 I=1,NTIM
+ IF(ABS(TIMES(I)-XT(1)).LE.1.0E-4*XT(1)) ITIM=I
+ 240 CONTINUE
+ IF(ITIM.EQ.0) CALL XABORT('EVODRV: TABLE LOOKUP FAILURE.')
+ NLENGT=(NVAR+1)*(NREAC+1)*NBMIX
+ DO IEX=1,MIN(ITIM-IEXTR+1,IEXTR)
+ WRITE(TEXT12,'(8HDEPL-DAT,I4.4)') ITIM-IEX
+ CALL LCMSIX(IPDEPL,TEXT12,1)
+ CALL LCMLEN(IPDEPL,'MICRO-RATES',LENGT,ITYLCM)
+ IF(LENGT.NE.NLENGT) THEN
+ CALL XABORT('EVODRV: MICRO-RATES OVERFLOW.')
+ ENDIF
+ CALL LCMGET(IPDEPL,'MICRO-RATES',SIG(1,1,1,1-IEX))
+ CALL LCMGET(IPDEPL,'INT-FLUX',VPHV(1,1-IEX))
+ CALL LCMSIX(IPDEPL,' ',2)
+ ENDDO
+ N=0
+ IF((IEXTR.GE.2).AND.(ITIM.GE.3)) THEN
+ ! parabolic extrapolation
+ N=3
+ T(1)=TIMES(ITIM-2)
+ T(2)=TIMES(ITIM-1)
+ T(3)=TIMES(ITIM)
+ IF(T(1).GE.T(3)) CALL XABORT ('EVODRV: T(1).GE.T(3).')
+ IF(IMPX.GT.0) WRITE(IUNOUT,'(/21H EVODRV: PARABOLIC EX,
+ 1 51HTRAPOLAPOLATION OF MICRO REACTION RATES AT END-OF-S,
+ 2 5HTAGE.)')
+ ELSE IF((IEXTR.GE.1).AND.(ITIM.GE.2)) THEN
+ ! linear extrapolation
+ N=2
+ T(1)=TIMES(ITIM-1)
+ T(2)=TIMES(ITIM)
+ IF(T(1).GE.T(2)) CALL XABORT ('EVODRV: T(1).GE.T(2).')
+ IF(IMPX.GT.0) WRITE(IUNOUT,'(/23H EVODRV: LINEAR EXTRAPO,
+ 1 51HLAPOLATION OF MICRO REACTION RATES AT END-OF-STAGE.)')
+ ELSE IF(ITIM.EQ.1) THEN
+ N=1
+ T(1)=TIMES(ITIM)
+ IF(IMPX.GT.0) WRITE(IUNOUT,'(/21H EVODRV: NO EXTRAPOLA,
+ 1 49HPOLATION OF MICRO REACTION RATES AT END-OF-STAGE.)')
+ ENDIF
+ DPD=1.0D0 ! perform Lagrange extrapolation
+ DO I=1,N
+ DPD=(XT(2)-T(I))*DPD
+ ENDDO
+ VPHV(:NBMIX,2)=0.0
+ SIG(:NVAR+1,:NREAC+1,:NBMIX,2)=0.0
+ DO I=1,N
+ WEI=DPD/(XT(2)-T(I))
+ DO J=1,N
+ IF(J.EQ.I) CYCLE
+ IF(T(I).EQ.T(J)) CALL XABORT('EVODRV: DIVIDE CHECK.')
+ WEI=WEI/(T(I)-T(J))
+ ENDDO
+ DO IBM=1,NBMIX
+ VPHV(IBM,2)=VPHV(IBM,2)+REAL(WEI)*VPHV(IBM,I-N+1)
+ SIG(:NVAR+1,:NREAC+1,IBM,2)=SIG(:NVAR+1,:NREAC+1,IBM,2)+
+ 1 REAL(WEI)*SIG(:NVAR+1,:NREAC+1,IBM,I-N+1)
+ ENDDO
+ ENDDO
+ ELSE
+ DO IBM=1,NBMIX
+ VPHV(IBM,2)=VPHV(IBM,1)
+ SIG(:NVAR+1,:NREAC+1,IBM,2)=SIG(:NVAR+1,:NREAC+1,IBM,1)
+ ENDDO
+ IF(IMPX.GT.0) WRITE(IUNOUT,'(/23H EVODRV: USE BEGINNING-,
+ 1 46HOF-STAGE MICRO REACTION RATES AT END-OF-STAGE.)')
+ ENDIF
+ ENDIF
+ VPHD=0.0
+ DO 250 ICMB=1,NCOMB
+ IBM=MILVO(ICMB)
+ IF(IBM.GT.0) VPHD=VPHD+VPHV(IBM,IP)
+ 250 CONTINUE
+ VPH(IP)=VPHD
+*
+ DO 270 ICMB=1,NCOMB
+ IBM=MILVO(ICMB)
+ IF(IBM.EQ.0) GO TO 270
+ DO 260 IS=1,NVAR
+ IF(JM(IBM,IS).GT.0) THEN
+ YDPL(IS,IP,ICMB)=DEN(JM(IBM,IS))
+ ELSE
+ YDPL(IS,IP,ICMB)=0.0
+ ENDIF
+ 260 CONTINUE
+ 270 CONTINUE
+*
+ IF(INR.NE.0) THEN
+* CHECK FOR OUT-OF-CORE DEPLETION.
+ DO 290 IBM=1,NBMIX
+ DO 280 IU=1,NGROUP
+ LCOOL=LCOOL.AND.(FLUMIX(IU,IBM).EQ.0.)
+ 280 CONTINUE
+ 290 CONTINUE
+ ENDIF
+ 300 CONTINUE
+*
+ IF(LCOOL.AND.(FIT.NE.0.0)) CALL XABORT('EVODRV: NEUTRON FLUX I'
+ 1 //'S ZERO. UNABLE TO NORMALIZE.')
+ IF(LCOOL.AND.(IMPX.GT.1)) THEN
+ WRITE (IUNOUT,'(/31H EVODRV: OUT-OF-CORE DEPLETION.)')
+ ELSE IF(IMPX.GT.1) THEN
+ WRITE (IUNOUT,'(/27H EVODRV: IN-CORE DEPLETION.)')
+ IF((FUELDN(3).GT.0.0).AND.(INR.EQ.3)) THEN
+ WRITE(IUNOUT,'(/31H EVODRV: FUEL POWER NORMALISATI,
+ 1 10HON FACTOR=,1P,E12.5,29H MW/TONNE. OUT-OF-FUEL POWER ,
+ 2 12HFORM FACTOR=,E12.5)') FIT/FUELDN(3)/PFACT,PFACT
+ ENDIF
+ ENDIF
+ INR2=INR
+ IF(LCOOL) INR2=0
+*----
+* PERFORM THE DEPLETION CALCULATION
+*----
+ CALL EVOBLD(IMPX,INR2,IGLOB,NBMIX,NBISO,NCOMB,ISONAM,YDPL,VX,
+ 1 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,
+ 3 AWR,FUELDN,SIG(1,1,1,1),VPH,VPHV(1,1),MIXPWR,VTOTD,IEVOLB,
+ 4 KFISS,KPF)
+*----
+* SAVE THE INITIAL SATURATED NUMBER DENSITIES IN THE DEPLETION TABLE
+*----
+ IF((ISAVE.GE.0).AND.(ISAT.EQ.1)) THEN
+ DO 315 ICMB=1,NCOMB
+ IBM=MILVO(ICMB)
+ IF(IBM.EQ.0) GO TO 315
+ DO 310 IS=1,NVAR
+ IF(JM(IBM,IS).GT.0) DEN(JM(IBM,IS))=YDPL(IS,1,ICMB)
+ 310 CONTINUE
+ 315 CONTINUE
+ ITIM=0
+ DO 320 I=1,NTIM
+ IF(ABS(TIMES(I)-XT(1)).LE.1.0E-4*XT(1)) ITIM=I
+ 320 CONTINUE
+ IF(ITIM.EQ.0) CALL XABORT('EVODRV: MISSING TIME ENTRY.')
+ WRITE(TEXT12,'(8HDEPL-DAT,I4.4)') ITIM
+ IF(IMPX.GT.0) WRITE(IUNOUT,530) XT(1),XT(1)/8.64E-4,TEXT12
+ CALL LCMSIX(IPDEPL,TEXT12,1)
+ CALL LCMPUT(IPDEPL,'ISOTOPESDENS',NBISO,2,DEN)
+ CALL LCMSIX(IPDEPL,' ',2)
+ ENDIF
+*----
+* SAVE THE DEPLETION CALCULATION RESULT IN THE DEPLETION TABLE
+*----
+ DO 335 ICMB=1,NCOMB
+ IBM=MILVO(ICMB)
+ IF(IBM.EQ.0) GO TO 335
+ DO 330 IS=1,NVAR
+ IF(JM(IBM,IS).GT.0) DEN(JM(IBM,IS))=YDPL(IS,2,ICMB)
+ 330 CONTINUE
+ 335 CONTINUE
+ ITIM=0
+ DO 340 I=1,NTIM
+ IF(ABS(TIMES(I)-XT(2)).LE.1.0E-4*XT(2)) ITIM=I
+ 340 CONTINUE
+ IF(ITIM.EQ.0) THEN
+ IF(XT(2).LT.TIMES(NTIM)) CALL XABORT('EVODRV: INVALID X2')
+ NTIM=NTIM+1
+ TIMES(NTIM)=XT(2)
+ CALL LCMPUT(IPDEPL,'DEPL-TIMES',NTIM,2,TIMES)
+ ITIM=NTIM
+ ENDIF
+ WRITE(TEXT12,'(8HDEPL-DAT,I4.4)') ITIM
+ IF(IMPX.GT.0) WRITE(IUNOUT,530) XT(2),XT(2)/8.64E-4,TEXT12
+ CALL LCMSIX(IPDEPL,TEXT12,1)
+ CALL LCMPUT(IPDEPL,'ISOTOPESDENS',NBISO,2,DEN)
+ NLENGT=(NVAR+1)*(NREAC+1)*NBMIX
+ CALL LCMPUT(IPDEPL,'MICRO-RATES',NLENGT,2,SIG(1,1,1,2))
+ CALL LCMPUT(IPDEPL,'INT-FLUX',NBMIX,2,VPHV(1,2))
+ CALL LCMPUT(IPDEPL,'ENERG-MIX',NBMIX,2,ENERG)
+* We use DELTA(3) instead of DELTA(1) in order to avoid different
+* base points in multi-D tables.
+ BRNWIR(1)=BRNWIR(1)+DELTA(3)/8.64E-4
+ BRNWIR(2)=BRNWIR(2)+DELTA(2)
+ CALL LCMPUT(IPDEPL,'BURNUP-IRRAD',2,2,BRNWIR)
+ IF(IMPX.GE.1) WRITE(IUNOUT,580) XT(2)/8.64E-4,
+ > BRNWIR(1),BRNWIR(2)
+ CALL LCMSIX(IPDEPL,' ',2)
+ ENDIF
+*----
+* RELEASE THE ALLOCATED MEMORY
+*----
+ DEALLOCATE(IDR,HREAC,KPAR)
+ DEALLOCATE(RRD,RER,YIELD,BPAR)
+ DEALLOCATE(KPF,KFISS)
+*----
+* USE THE RESULT OF A DEPLETION CALCULATION IN THE FOLLOWING RUN
+*----
+ IF(ISET.GE.0) THEN
+ ITIM=0
+ DO 350 I=1,NTIM
+ IF(ABS(TIMES(I)-XTF).LE.1.0E-4*XTF) ITIM=I
+ 350 CONTINUE
+ IF(ITIM.EQ.0) CALL XABORT('EVODRV: NO DEPLETION DATA STORED.')
+ WRITE(TEXT12,'(8HDEPL-DAT,I4.4)') ITIM
+ IF(IMPX.GT.0) WRITE(IUNOUT,520) XTF,XTF/8.64E-4,TEXT12
+ CALL LCMSIX(IPDEPL,TEXT12,1)
+ CALL LCMGET(IPDEPL,'ISOTOPESDENS',DEN)
+ CALL LCMGET(IPDEPL,'BURNUP-IRRAD',BRNWIR)
+ CALL LCMSIX(IPDEPL,' ',2)
+ IF(IMPX.GT.1) THEN
+ WRITE(IUNOUT,550) XTF,XTF/8.64E-4
+ DO 370 ICMB=1,NCOMB
+ IMIXC=MILVO(ICMB)
+ IF(IMIXC.EQ.0) GO TO 370
+ NISOCC=0
+ DO 360 ISOT=1,NBISO
+ IMIXI=MIX(ISOT)
+ IF(IMIXI.EQ.IMIXC) THEN
+ NISOCC=NISOCC+1
+ ISOCMB(NISOCC)=ISOT
+ ENDIF
+ 360 CONTINUE
+ WRITE(IUNOUT,560) IMIXC
+ WRITE(IUNOUT,570) ((ISONAM(I0,ISOCMB(I)),I0=1,2),
+ 1 DEN(ISOCMB(I)),I=1,NISOCC)
+ 370 CONTINUE
+ ENDIF
+ ENDIF
+*----
+* RECOVER THE BURNUP AND SAVE IT IN A CLE-2000 VARIABLE
+*----
+ IF(IPICK.EQ.1) THEN
+ ITIM=0
+ DO 375 I=1,NTIM
+ IF(ABS(TIMES(I)-XTF).LE.1.0E-4*XTF) ITIM=I
+ 375 CONTINUE
+ IF(ITIM.EQ.0) CALL XABORT('EVODRV: NO DEPLETION DATA STORED.')
+ WRITE(TEXT12,'(8HDEPL-DAT,I4.4)') ITIM
+ IF(IMPX.GT.0) WRITE(IUNOUT,520) XTF,XTF/8.64E-4,TEXT12
+ CALL LCMSIX(IPDEPL,TEXT12,1)
+ CALL LCMGET(IPDEPL,'BURNUP-IRRAD',BRNWIR)
+ CALL LCMSIX(IPDEPL,' ',2)
+ CALL REDGET(ITYPLU,INTLIR,REALIR,TEXT12,DBLLIR)
+ IF(ITYPLU.NE.-2) CALL XABORT('EVODRV: OUTPUT REAL EXPECTED.')
+ ITYPLU=2
+ REALIR=BRNWIR(1)
+ IF(IMPX.GT.2) WRITE(IUNOUT,540) REALIR
+ CALL REDPUT(ITYPLU,INTLIR,REALIR,TEXT12,DBLLIR)
+ CALL REDGET(ITYPLU,INTLIR,REALIR,TEXT12,DBLLIR)
+ IF((ITYPLU.NE.3).OR.(TEXT12.NE.';')) THEN
+ CALL XABORT('EVODRV: ; CHARACTER EXPECTED.')
+ ENDIF
+ ENDIF
+*
+ IF((IDEPL.EQ.1).OR.(ISET.EQ.1).AND.LMACRO) THEN
+* COMPUTE THE NEW DEPLETED MACROSCOPIC CROSS SECTIONS.
+ MASKL(:NGROUP)=.TRUE.
+ MASK(:NBMIX)=.FALSE.
+ DO 380 ICOMB=1,NCOMB
+ IBM=MILVO(ICOMB)
+ IF(IBM.GT.0) MASK(IBM)=.TRUE.
+ 380 CONTINUE
+*
+ ITSTMP=2
+ TMPDAY(1)=XT(2)/8.64E-4
+ TMPDAY(2)=BRNWIR(1)
+ TMPDAY(3)=BRNWIR(2)
+ CALL LIBMIX(IPLIB,NBMIX,NGROUP,NBISO,ISONAM,MIX,DEN,MASK,MASKL,
+ 1 ITSTMP,TMPDAY)
+ ENDIF
+*----
+* STORE THE GENERAL DEPLETION RELATED PARAMETERS
+*----
+ CALL LCMPUT(IPDEPL,'VOLUME-MIX',NBMIX,2,VX)
+ CALL LCMPUT(IPDEPL,'DEPLETE-MIX',NVAR*NBMIX,1,JM)
+ CALL LCMPUT(IPDEPL,'MIXTURESBurn',NBMIX,1,MIXBRN)
+ CALL LCMPUT(IPDEPL,'MIXTURESPowr',NBMIX,1,MIXPWR)
+ IPAR(:NSTATE)=0
+ IPAR(1)=ITYPE
+ IPAR(2)=INR
+ IPAR(3)=NTIM
+ IPAR(4)=NBISO
+ IPAR(5)=NCOMB
+ IPAR(6)=NREAC
+ IPAR(7)=NVAR
+ IPAR(8)=NBMIX
+ IPAR(9)=IEXTR
+ IPAR(10)=IGLOB
+ IPAR(11)=ISAT
+ IPAR(12)=IDIRAC
+ IPAR(13)=ITIXS
+ IPAR(14)=IFLMAC
+ IPAR(15)=IYLMIX
+ CALL LCMPUT(IPDEPL,'STATE-VECTOR',NSTATE,1,IPAR)
+ RPAR(1)=EPS1
+ RPAR(2)=EPS2
+ RPAR(3)=EXPMAX
+ RPAR(4)=H1
+ RPAR(5)=FIT
+ CALL LCMPUT(IPDEPL,'EVOLUTION-R',5,2,RPAR)
+ IF((IMPX.GT.1).OR.((IMPX.GT.0).AND.(INDREC.EQ.1))) THEN
+ WRITE(IUNOUT,590) IMPX,ITYPE,INR,NTIM,NBISO,NCOMB,NREAC,
+ 1 NVAR,NBMIX,IEXTR
+ WRITE(IUNOUT,595) IGLOB,ISAT,IDIRAC,ITIXS,IFLMAC,IYLMIX
+ ENDIF
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(TIMES)
+ DEALLOCATE(IPISO)
+ DEALLOCATE(MASKL,MASK)
+ DEALLOCATE(YDPL,IZAE,AWR,ENERG,VPHV,SIG)
+ DEALLOCATE(IEVOLB,INADPL,ISOCMB,MILVO,JM)
+ RETURN
+*
+ 500 FORMAT(/36H EVODRV: NUMBER OF FISSILE ISOTOPES=,I4/9X,8HNUMBER O,
+ 1 19HF FISSION PRODUCTS=,I4)
+ 510 FORMAT(/24H EVODRV: UNABLE TO FIND ,A15,6H INDEX,I5,6H AMONG,10I5
+ 1 /(56X,10I5))
+ 520 FORMAT(/44H EVODRV: RECOVER INFORMATION RELATED TO TIME,1P,E12.4,
+ 1 8H E+8 S (,E12.4,32H DAY) FROM LCM DIRECTORY NAMED ',A12,2H'.)
+ 530 FORMAT(/41H EVODRV: SAVE INFORMATION RELATED TO TIME,1P,E12.4,
+ 1 8H E+8 S (,E12.4,30H DAY) ON LCM DIRECTORY NAMED ',A12,2H'.)
+ 540 FORMAT(/21H EVODRV: PICK BURNUP=,1P,E12.4,10H MWd/tonne)
+ 550 FORMAT(/' EVODRV: NUMBER DENSITIES PER ISOTOPE AT TIME',1P,
+ 1 E12.4,' E+8 S (',E12.4,' DAY)')
+ 560 FORMAT(/' ISOTOPIC DENSITIES AFTER BURNUP FOR MIXTURE = ',I5,
+ 1 22H (10**24 PARTICLES/CC))
+ 570 FORMAT(1P,5(4X,2A4,':',E12.4))
+ 580 FORMAT(' -> FINAL BURNUP AT TIME = ',1P,E14.6,' DAYS'/
+ > ' FUEL BURNUP = ',E14.6,' MW*D/TONNE'/
+ > ' NEUTRON EXPOSURE = ',E14.6,' N/KB')
+ 590 FORMAT(/8H OPTIONS/8H -------/
+ 1 7H IMPX ,I8,30H (0=NO PRINT/1=SHORT/2=MORE)/
+ 2 7H ITYPE ,I8,31H (1=CASH-KARP/2=KAPS-RENTROP)/
+ 3 7H INR ,I8,47H (0=OUT-OF-CORE/1=CONSTANT FLUX/2=CONSTANT PO,
+ 4 4HWER)/
+ 5 7H NTIM ,I8,35H (NUMBER OF DEPLETION SET POINTS)/
+ 6 7H NBISO ,I8,29H (TOTAL NUMBER OF ISOTOPES)/
+ 7 7H NCOMB ,I8,33H (NUMBER OF DEPLETING MIXTURES)/
+ 8 7H NREAC ,I8,34H (NUMBER OF DEPLETING REACTIONS)/
+ 9 7H NVAR ,I8,33H (NUMBER OF DEPLETING ISOTOPES)/
+ 1 7H NBMIX ,I8,23H (NUMBER OF MIXTURES)/
+ 2 7H IEXTR ,I8,47H (FLUX EXTRAPOLATION: 0=NONE/1=LINEAR/2=PARAB,
+ 3 5HOLIC))
+ 595 FORMAT(
+ 1 7H IGLOB ,I8,47H (-1=SERPENT EDEPMODE-0 FORMULA/0=COMPUTE BUR,
+ 2 44HNUP IN FUEL/1=COMPUTE BURNUP IN GLOBAL CELL)/
+ 3 7H ISAT ,I8,47H (0/1=DO NOT/DO SAVE SATURATED INITIAL NUMBER,
+ 4 11H DENSITIES)/
+ 5 7H IDIRAC,I8,47H (0/1=DO NOT/DO USE DIRAC FUNCTION CONTRIBUTI,
+ 6 34HONS IN SATURATED NUMBER DENSITIES)/
+ 7 7H ITIXS ,I8,38H (0/1=TIME-DEPENDENT XS FLAG ON/OFF)/
+ 8 7H IFLMAC,I8,47H (0/1/2=RECOVER FLUX FROM L_FLUX/L_MACROLIB/L,
+ 9 7H_POWER)/
+ 1 7H IYLMIX,I8,47H (0/1=RECOVER FISSION YIELD DATA FROM DEPL-CH,
+ 2 16HAIN/PYIELD DATA))
+ 600 FORMAT(/54H EVODRV: SOLUTION OF A DEPLETION SYSTEM BETWEEN TIMES ,
+ 1 1P,E12.4,4H AND,E12.4,6H E+8 S)
+ 610 FORMAT(/' EVODRV: FUEL INITIAL DENSITY = ',1P,E14.6,' G/CC'/
+ > ' FUEL TOTAL VOLUME = ',E14.6,' CC'/
+ > ' FUEL INITIAL MASS = ',E14.6,' G'/
+ > ' FUEL INITIAL MASS/CELL VOL = ',E14.6,' G/CC')
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 15:16:10 +0000