*DECK EVODRV
      SUBROUTINE EVODRV(IPDEPL,IPLIB,INDREC,IMPX,NBISO,NGROUP,NBMIX,
     1 NREGIO,ISONAM,ISONRF,MIX,DEN,IEVOL,MAT,VOL,NDEPL,NSUPS,NREAC,
     2 MCOMB,EPS1,EPS2,EXPMAX,H1,ITYPE,INR,IEXTR,IGLOB,ISAT,IDIRAC,
     3 ITIXS,FIT,ISAVE,ISET,IDEPL,XTI,XTF,XT,LMACRO,FLUXES,MASK,MASKL,
     4 JM,MILVO,YDPL,VX,ISOCMB,INADPL,SIG,VPHV,ENERG,AWR,
     5 MIXBRN,MIXPWR,MILPW)
*
*-----------------------------------------------------------------------
*
*Purpose:
* isotopic depletion calculation.
*
*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.
* NREGIO  number of volumes.
* 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).
* MAT     index-number of the mixture type assigned to each volume.
* VOL     volumes.
* NDEPL   number of isotopes in the depletion chain.
* NSUPS   number of non-depleting isotopes producing energy.
* NREAC   maximum number of depletion reactions.
* MCOMB   maximum 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
*         extrapolation).
* IGLOB   out-of-fuel power in flux normalization:
*         =0: compute the burnup using the power released in the fuel;
*         =1: compute the burnup 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).
* 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(1)   initial time (independent variable) for the depletion
*         calculation.
* XT(2)   final time (independent variable) for the depletion
*         calculation.
* LMACRO  macrolib building flag (=.true. to compute the embedded
*         macrolib).
* FLUXES  resulting average fluxes obtained using the input internal
*         library.
* MIXBRN  flags for mixtures to burn.
* MIXPWR  flags for mixtures to include in power normalization.
*
*Parameters: scratch
* MASK
* MASKL
* JM
* YDPL
* VX
* ISOCMB
* INADPL
* SIG
* VPHV
* ENERG
* AWR
* MILVO   mixture index corresponding to each depleting mixture.
* MILPW   flag for power evaluation for each depleting mixture.
*
*-----------------------------------------------------------------------
*
*----
*  SUBROUTINE ARGUMENTS
*----
      INTEGER IPDEPL,IPLIB,INDREC,IMPX,NBISO,NGROUP,NBMIX,NREGIO,
     1 ISONAM(3,NBISO),ISONRF(3,NBISO),MIX(NBISO),IEVOL(NBISO),
     2 MAT(NREGIO),NDEPL,NSUPS,NREAC,MCOMB,ITYPE,INR,IEXTR,IGLOB,
     3 ISAT,IDIRAC,ITIXS,ISAVE,ISET,IDEPL,JM(NBMIX,NDEPL),MILVO(MCOMB),
     4 ISOCMB(NBISO),INADPL(3,NDEPL)
      REAL DEN(NBISO),VOL(NREGIO),EPS1,EPS2,EXPMAX,H1,FIT,XTI,XTF,
     1 XT(2),FLUXES(NREGIO,NGROUP),YDPL(NDEPL-NSUPS+1,2,MCOMB),
     2 VX(NBMIX),SIG(NDEPL-NSUPS+1,NREAC+1,NBMIX,0:2),VPHV(NBMIX,0:2),
     3 ENERG(NBMIX),AWR(NDEPL)
      LOGICAL LMACRO,MASK(NBMIX),MASKL(NGROUP)
      INTEGER MIXBRN(NBMIX),MIXPWR(NBMIX),MILPW(MCOMB)
*----
*  LOCAL VARIABLES
*----
      PARAMETER (IUNOUT=6,NSTATE=40,MAXTIM=9999)
      CHARACTER TEXT12*12
      LOGICAL LCOOL
      INTEGER IDIM(NSTATE),IPAR(NSTATE)
      REAL TIMES(MAXTIM),DELTA(3),RPAR(5),BRNWIR(2),TMPDAY(3),VPH(2),
     1 FUELDN(3),DELTAT(2,2),TIMEP(2,3)
      DOUBLE PRECISION XDRCST,AVCON,VTOTD
      COMMON BASE(1)
*----
*  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.EQ.0) THEN
        CALL XDRSET(TIMES,MAXTIM,0.0)
        CALL LCMPUT(IPDEPL,'DEPL-TIMES',1,2,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)
      CALL SETARA(BASE,NDEPL*NPAR,IKPAR)
      CALL SETARA(BASE,NDEPL*NPAR,IBPAR)
      CALL SETARA(BASE,NDFP*NFISS,IYIELD)
      CALL SETARA(BASE,2*NREAC,IHREAC)
      CALL SETARA(BASE,NREAC*NDEPL,IIDR)
      CALL SETARA(BASE,NREAC*NDEPL,IRER)
      CALL SETARA(BASE,NDEPL,IRRD)
      CALL LCMGET(IPLIB,'ISOTOPESDEPL',INADPL)
      CALL LCMGET(IPLIB,'PRODUCE-REAC',BASE(IKPAR))
      CALL LCMGET(IPLIB,'PRODUCE-RATE',BASE(IBPAR))
      IF(NDFP.GT.0) CALL LCMGET(IPLIB,'FISSIONYIELD',BASE(IYIELD))
      CALL LCMGET(IPLIB,'DEPLETE-IDEN',BASE(IHREAC))
      CALL LCMGET(IPLIB,'DEPLETE-REAC',BASE(IIDR))
      CALL LCMGET(IPLIB,'DEPLETE-ENER',BASE(IRER))
      CALL LCMGET(IPLIB,'DEPLETE-DECA',BASE(IRRD))
      CALL LCMSIX(IPLIB,' ',2)
*----
*  DETECT THE DEPLETING ISOTOPES AND MIXTURES IN THE MICROLIB.
*----
      ICOMB=0
      CALL XDISET(JM,NBMIX*NDEPL,0)
      DO 30 ISOT=1,NBISO
      IBM=MIX(ISOT)
      IF(IBM.EQ.0) GO TO 30
      IF(MIXBRN(IBM) .EQ. 0) GO TO 30
      DO 20 INUCL=1,NDEPL
      IF((ISONRF(1,ISOT).EQ.INADPL(1,INUCL)).AND.(ISONRF(2,ISOT).EQ.
     1 INADPL(2,INUCL)).AND.(ISONRF(3,ISOT).EQ.INADPL(3,INUCL))) THEN
         IF(INUCL.GT.NVAR) THEN
            JM(IBM,INUCL)=-ISOT
         ELSE IF(IEVOL(ISOT).EQ.1) THEN
            JM(IBM,INUCL)=-ISOT
         ELSE IF(IEVOL(ISOT).EQ.0) THEN
            JM(IBM,INUCL)=ISOT
            DO 10 J=1,ICOMB
            IF(IBM.EQ.MILVO(J)) GO TO 30
   10       CONTINUE
            ICOMB=ICOMB+1
            MILVO(ICOMB)=IBM
            IF(MIXPWR(IBM) .EQ. 1) THEN
              MILPW(ICOMB)=1
            ELSE
              MILPW(ICOMB)=0
            ENDIF
         ENDIF
         GO TO 30
      ENDIF
   20 CONTINUE
   30 CONTINUE
      IF(ICOMB.GT.MCOMB) CALL XABORT('EVODRV: INVALID VALUE OF MCOMB.')
      NCOMB=ICOMB
*----
*  COMPUTE THE VOLUMES OF THE MIXTURES
*----
      CALL XDRSET(VX,NBMIX,0.0)
      DO 40 I=1,NREGIO
      IBM=MAT(I)
      IF(IBM.GT.0) VX(IBM)=VX(IBM)+VOL(I)
   40 CONTINUE
*----
*  COMPUTE THE INITIAL INTEGRATED FLUX
*----
      VPHINI=0.0
      DO ICMB=1,NCOMB
        IBM=MILVO(ICMB)
        IF(MILPW(ICMB) .EQ. 1) THEN
          IF(VX(IBM).NE.0.0) THEN
            DO I=1,NREGIO
              IF(MAT(I).EQ.IBM) THEN
                DO IU=1,NGROUP
                  VPHINI=VPHINI+FLUXES(I,IU)*VOL(I)
                ENDDO
              ENDIF
            ENDDO
          ENDIF
        ENDIF
      ENDDO
*----
*  RECOVER THE NUCLIDE MASS RATIOS.
*----
      DO 80 INUCL=1,NDEPL
      AWR(INUCL)=0.0
      DO 80 IBM=1,NBMIX
      K=ABS(JM(IBM,INUCL))
      IF(K.EQ.0) GO TO 80
      WRITE(TEXT12,'(3A4)') (ISONAM(I0,K),I0=1,3)
      CALL LCMSIX(IPLIB,TEXT12,1)
      CALL LCMGET(IPLIB,'AWR',AWR(INUCL))
      CALL LCMSIX(IPLIB,' ',2)
   80 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
   90       DELTAT(2,IP)=XT(IP)/8.64E-4-TMPREF
         ELSE
            XREF=XT(1)/8.64E-4
            DO 100 IP=1,2
            DELTAT(1,IP)=1.0
  100       DELTAT(2,IP)=XREF-TMPREF
         ENDIF
      ELSE
         DO 110 IP=1,2
         DELTAT(1,IP)=1.0
  110    DELTAT(2,IP)=0.0
      ENDIF
*----
*  COMPUTE AND SAVE THE INITIAL MASS OF HEAVY ELEMENTS
*----
      VTOTD=0.0D0
      DO ICMB=1,NCOMB
        IBM=MILVO(ICMB)
        IF(MILPW(ICMB) .GT. 0) THEN
          VTOTD=VTOTD+DBLE(VX(IBM))
        ENDIF
      ENDDO
      IF(INDREC.EQ.1) THEN
*
*        COMPUTE THE GLOBAL HEAVY-ELEMENT MASS.
*        for isotopes in MIXPWR
        FUELDN(1)=0.0
        DO ICMB=1,NCOMB
          IBM=MILVO(ICMB)
          IF(MILPW(ICMB) .GT. 0) THEN
            DO IS=1,NDEPL
              IF((JM(IBM,IS).GT.0).AND.(AWR(IS).GT.210.0)) THEN
                FUELDN(1)=FUELDN(1)+AWR(IS)*DEN(JM(IBM,IS))*VX(IBM)
              ENDIF
            ENDDO
          ENDIF
        ENDDO
        FUELDN(1)=FUELDN(1)*AVCON/VTOTD
        FUELDN(2)=FUELDN(1)*VTOTD
        VASSMB=0.0
        DO IBM=1,NBMIX
          IF(MIXPWR(IBM) .EQ. 1) THEN
            VASSMB=VASSMB+VX(IBM)
          ENDIF
        ENDDO
        FUELDN(3)=FUELDN(2)/VASSMB
        CALL LCMPUT(IPDEPL,'FUELDEN-INIT',3,2,FUELDN)
*
*        COMPUTE THE HEAVY-ELEMENT MASS PER MIXTURE.
        DO IBM=1,NBMIX
          ENERG(IBM)=0.0
          DO IS=1,NDEPL
            IF((JM(IBM,IS).GT.0).AND.(AWR(IS).GT.210.0)) THEN
              ENERG(IBM)=ENERG(IBM)+AWR(IS)*DEN(JM(IBM,IS))*VX(IBM)
            ENDIF
          ENDDO
          ENERG(IBM)=ENERG(IBM)*AVCON
        ENDDO
        CALL LCMPUT(IPDEPL,'FUELDEN-MIX',NBMIX,2,ENERG)
*
*        COMPUTE THE TOTAL MASS PER MIXTURE.
        CALL XDRSET(ENERG,NBMIX,0.0)
        DO IS=1,NBISO
          IF(DEN(IS).NE.0.0) THEN 
            WRITE(TEXT12,'(3A4)') (ISONAM(I0,IS),I0=1,3)
            CALL LCMSIX(IPLIB,TEXT12,1)
            CALL LCMGET(IPLIB,'AWR',AWRGAR)
            CALL LCMSIX(IPLIB,' ',2)
            DO IBM=1,NBMIX
              IF(MIX(IS).EQ.IBM) THEN
                ENERG(IBM)=ENERG(IBM)+AWRGAR*DEN(IS)*VX(IBM)*AVCON
              ENDIF
            ENDDO
          ENDIF
        ENDDO
        CALL LCMPUT(IPDEPL,'WEIGHT-MIX',NBMIX,2,ENERG)
      ELSE
        CALL LCMGET(IPDEPL,'FUELDEN-INIT',FUELDN)
      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
          CALL XDRSET(TIMES,MAXTIM,0.0)
          CALL LCMGET(IPDEPL,'DEPL-TIMES',TIMES)
          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
          IF(NTIM.GT.MAXTIM) CALL XABORT('EVODRV: INVALID MAXTIM(1).')
          TIMES(NTIM)=XT(IP)
          CALL LCMPUT(IPDEPL,'DEPL-TIMES',NTIM,2,TIMES)
          ITIM=NTIM
        ENDIF
        WRITE(TEXT12,'(8HDEPL-DAT,I4.4)') ITIM
        IF(ABS(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 SETARA(BASE,NFISS,IKFISS)
        CALL SETARA(BASE,NGROUP*NBMIX,IFLX)
        CALL SETARA(BASE,NGROUP*(NREAC-1),IXSREC)
        CALL EVOSIG(IPLIB,IMPX,INR,IGLOB,NGROUP,NREGIO,NBMIX,NBISO,
     1  NCOMB,ISONAM,DEN,MAT,VOL,FLUXES,VX,MILVO,JM,NVAR,NSUPS,NREAC,
     2  NFISS,BASE(IHREAC),BASE(IIDR),BASE(IRER),BASE(IRRD),FIT,FUELDN,
     3  NXSPER,DELTAT(1,IP),PFACT,SIG(1,1,1,IP),VPHV(1,IP),BASE(IKFISS),
     4  BASE(IFLX),BASE(IXSREC),VTOTD,MIXPWR,MILPW)
        CALL RLSARA(BASE(IXSREC))
        CALL RLSARA(BASE(IFLX))
        CALL RLSARA(BASE(IKFISS))
        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 ICMB=1,NCOMB
          IF(MILPW(ICMB). EQ. 1) VPHD=VPHD+VPHV(MILVO(ICMB),IP)
        ENDDO
        VPH(IP)=VPHD
        IF(INR.NE.0) THEN
           FNORM=VPH(IP)/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
          CALL XDRSET(TIMES,MAXTIM,0.0)
          CALL LCMGET(IPDEPL,'DEPL-TIMES',TIMES)
          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
          IF(NTIM.GT.MAXTIM) CALL XABORT('EVODRV: INVALID MAXTIM(2).')
          TIMES(NTIM)=XTI
          CALL LCMPUT(IPDEPL,'DEPL-TIMES',NTIM,2,TIMES)
          ITIM=NTIM
        ENDIF
        WRITE(TEXT12,'(8HDEPL-DAT,I4.4)') ITIM
        IF(ABS(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 SETARA(BASE,NFISS,IKFISS)
        CALL SETARA(BASE,NGROUP*NBMIX,IFLX)
        CALL SETARA(BASE,NGROUP*(NREAC-1),IXSREC)
        CALL EVOSIG(IPLIB,IMPX,INR,IGLOB,NGROUP,NREGIO,NBMIX,NBISO,
     1  NCOMB,ISONAM,DEN,MAT,VOL,FLUXES,VX,MILVO,JM,NVAR,NSUPS,NREAC,
     2  NFISS,BASE(IHREAC),BASE(IIDR),BASE(IRER),BASE(IRRD),FIT,FUELDN,
     3  NXSPER,DELTAT(1,IP),PFACT,SIG(1,1,1,IP),VPHV(1,IP),BASE(IKFISS),
     4  BASE(IFLX),BASE(IXSREC),VTOTD,MIXPWR,MILPW)
        CALL RLSARA(BASE(IXSREC))
        CALL RLSARA(BASE(IFLX))
        CALL RLSARA(BASE(IKFISS))
        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 ICMB=1,NCOMB
          IF(MILPW(ICMB). EQ. 1) VPHD=VPHD+VPHV(MILVO(ICMB),IP)
        ENDDO
        VPH(IP)=VPHD
        IF(INR.NE.0) THEN
           FNORM=VPH(IP)/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).
*----
        CALL XDRSET(TIMES,MAXTIM,0.0)
        CALL LCMGET(IPDEPL,'DEPL-TIMES',TIMES)
        IF(ABS(IMPX).GT.1) 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(ABS(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.EQ.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) THEN
                CALL XABORT('EVODRV: TABLE LOOKUP FAILURE.')
             ELSE IF(ITIM.GE.2) THEN
                T1=TIMES(ITIM-1)
                T2=TIMES(ITIM)
                IF(T1.GE.T2) CALL XABORT ('EVODRV: ALGORITHM FAILURE.')
                FACT=(XT(2)-T2)/(T2-T1)
                WRITE(TEXT12,'(8HDEPL-DAT,I4.4)') ITIM-1
                CALL LCMSIX(IPDEPL,TEXT12,1)
                CALL LCMGET(IPDEPL,'MICRO-RATES',SIG(1,1,1,0))
                CALL LCMGET(IPDEPL,'INT-FLUX',VPHV(1,0))
                CALL LCMSIX(IPDEPL,' ',2)
                DO 250 IBM=1,NBMIX
                VPHV(IBM,2)=VPHV(IBM,1)+FACT*(VPHV(IBM,1)-VPHV(IBM,0))
                DO 250 IQ=1,NREAC+1
                DO 250 IS=1,NVAR+1
  250           SIG(IS,IQ,IBM,2)=SIG(IS,IQ,IBM,1)+
     1             FACT*(SIG(IS,IQ,IBM,1)-SIG(IS,IQ,IBM,0))
                IF(ABS(IMPX).GT.0) 
     >          WRITE(IUNOUT,'(/20H EVODRV: USE EXTRAPO,
     1          43HLATED MICRO REACTION RATES AT END-OF-STAGE.)')
             ENDIF
          ELSE
             DO 260 IBM=1,NBMIX
             VPHV(IBM,2)=VPHV(IBM,1)
             DO 260 IQ=1,NREAC+1
             DO 260 IS=1,NVAR+1
  260        SIG(IS,IQ,IBM,2)=SIG(IS,IQ,IBM,1)
             IF(ABS(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 ICMB=1,NCOMB
          IF(MILPW(ICMB). EQ. 1) VPHD=VPHD+VPHV(MILVO(ICMB),IP)
        ENDDO
        VPH(IP)=VPHD
*
        DO 280 ICMB=1,NCOMB
        IBM=MILVO(ICMB)
        DO 280 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
  280   CONTINUE
*
        IF(INR.NE.0) THEN
*          CHECK FOR OUT-OF-CORE DEPLETION.
           DO 290 I=1,NREGIO
           DO 290 IU=1,NGROUP
           LCOOL=LCOOL.AND.(FLUXES(I,IU).EQ.0.)
  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.(ABS(IMPX).GT.1)) THEN
           WRITE (IUNOUT,'(/31H EVODRV: OUT-OF-CORE DEPLETION.)')
        ELSE IF(ABS(IMPX).GT.1) THEN
           WRITE (IUNOUT,'(/27H EVODRV: IN-CORE DEPLETION.)')
           IF(INR.EQ.3) WRITE(IUNOUT,'(/27H EVODRV: FUEL POWER NORMALI,
     1     14HSATION FACTOR=,1P,E12.5,29H MW/TONNE. OUT-OF-FUEL POWER ,
     2     12HFORM FACTOR=,E12.5)') FIT/FUELDN(3)/PFACT,PFACT
        ENDIF
        INR2=INR
        IF(LCOOL) INR2=0
*----
*  PERFORM THE DEPLETION CALCULATION
*----
        CALL SETARA(BASE,NVAR+1,IMU1)
        CALL SETARA(BASE,NVAR+1,IIMA)
        CALL SETARA(BASE,NVAR,ILP)
        CALL SETARA(BASE,NFISS,IKFISS)
        CALL SETARA(BASE,NDFP,IKPF)
        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,BASE(IKPAR),BASE(IBPAR),
     3  BASE(IYIELD),BASE(IIDR),BASE(IRER),BASE(IRRD),AWR,FUELDN,
     4  SIG(1,1,1,1),VPH,BASE(IMU1),BASE(IIMA),BASE(ILP),BASE(IKFISS),
     5  BASE(IKPF),MIXPWR,MILPW,VTOTD)
        CALL RLSARA(BASE(IKPF))
        CALL RLSARA(BASE(IKFISS))
        CALL RLSARA(BASE(ILP))
        CALL RLSARA(BASE(IIMA))
        CALL RLSARA(BASE(IMU1))
*----
*  SAVE THE INITIAL SATURATED NUMBER DENSITIES IN THE DEPLETION TABLE
*----
        IF((ISAVE.GE.0).AND.(ISAT.EQ.1)) THEN
          DO 310 ICMB=1,NCOMB
            IBM=MILVO(ICMB)
            DO 310 IS=1,NVAR
              IF(JM(IBM,IS).GT.0) DEN(JM(IBM,IS))=YDPL(IS,1,ICMB)
  310     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(ABS(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 330 ICMB=1,NCOMB
          IBM=MILVO(ICMB)
          DO 330 IS=1,NVAR
            IF(JM(IBM,IS).GT.0) DEN(JM(IBM,IS))=YDPL(IS,2,ICMB)
  330   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
          IF(NTIM.GT.MAXTIM) CALL XABORT('EVODRV: INVALID MAXTIM(3).')
          TIMES(NTIM)=XT(2)
          CALL LCMPUT(IPDEPL,'DEPL-TIMES',NTIM,2,TIMES)
          ITIM=NTIM
        ENDIF
        WRITE(TEXT12,'(8HDEPL-DAT,I4.4)') ITIM
        IF(ABS(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(ABS(IMPX).GE.1) WRITE(IUNOUT,580) XT(2)/8.64E-4,
     >                BRNWIR(1),BRNWIR(2)
        CALL LCMSIX(IPDEPL,' ',2)
      ENDIF
*----
*  RELEASE THE ALLOCATED MEMORY
*----
      CALL RLSARA(BASE(IRRD))
      CALL RLSARA(BASE(IRER))
      CALL RLSARA(BASE(IIDR))
      CALL RLSARA(BASE(IHREAC))
      CALL RLSARA(BASE(IYIELD))
      CALL RLSARA(BASE(IBPAR))
      CALL RLSARA(BASE(IKPAR))
*----
*  USE THE RESULT OF A DEPLETION CALCULATION IN THE FOLLOWING RUN
*----
      IF(ISET.GE.0) THEN
         CALL XDRSET(TIMES,MAXTIM,0.0)
         CALL LCMGET(IPDEPL,'DEPL-TIMES',TIMES)
         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(ABS(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(ABS(IMPX).GT.1) THEN
           WRITE(IUNOUT,550) XTF,XTF/8.64E-4
           DO 370 ICMB=1,NCOMB
             IMIXC=MILVO(ICMB)
             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
*
      IF((IDEPL.EQ.1).OR.(ISET.EQ.1).AND.LMACRO) THEN
*       COMPUTE THE NEW DEPLETED MACROSCOPIC CROSS SECTIONS.
        CALL XDLSET(MASKL,NGROUP,.TRUE.)
        CALL XDLSET(MASK,NBMIX,.FALSE.)
        DO 380 ICOMB=1,NCOMB
 380    MASK(MILVO(ICOMB))=.TRUE.
*
        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)
      CALL XDISET(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
      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((ABS(IMPX).GT.1).OR.
     >   ((ABS(IMPX).GT.0).AND.(INDREC.EQ.1))) THEN
         WRITE(IUNOUT,590) IMPX,ITYPE,INR,NTIM,NBISO,NCOMB,NREAC,
     1   NVAR,NBMIX,IEXTR,IGLOB,ISAT,IDIRAC,ITIXS
      ENDIF
      RETURN
*
  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(/54H EVODRV: *** WARNING *** FLUX MAY BE INCONSISTENT WITH,
     1 16H CROSS SECTIONS.)
  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 = ',I3,
     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  ,I6,30H   (0=NO PRINT/1=SHORT/2=MORE)/
     2 7H ITYPE ,I6,31H   (1=CASH-KARP/2=KAPS-RENTROP)/
     3 7H INR   ,I6,47H   (0=OUT-OF-CORE/1=CONSTANT FLUX/2=CONSTANT PO,
     4 4HWER)/
     5 7H NTIM  ,I6,35H   (NUMBER OF DEPLETION SET POINTS)/
     6 7H NBISO ,I6,29H   (TOTAL NUMBER OF ISOTOPES)/
     7 7H NCOMB ,I6,33H   (NUMBER OF DEPLETING MIXTURES)/
     8 7H NREAC ,I6,34H   (NUMBER OF DEPLETING REACTIONS)/
     9 7H NVAR  ,I6,33H   (NUMBER OF DEPLETING ISOTOPES)/
     1 7H NBMIX ,I6,23H   (NUMBER OF MIXTURES)/
     2 7H IEXTR ,I6,47H   (0=NO FLUX EXTRAPOLATION/1=FLUX EXTRAPOLATIO,
     3 2HN)/
     4 7H IGLOB ,I6,47H   (0=COMPUTE BURNUP IN FUEL/1=COMPUTE BURNUP I,
     5 14HN GLOBAL CELL)/
     6 7H ISAT  ,I6,47H   (0/1=DO NOT/DO SAVE SATURATED INITIAL NUMBER,
     7 11H DENSITIES)/
     8 7H IDIRAC,I6,47H   (0/1=DO NOT/DO USE DIRAC FUNCTION CONTRIBUTI,
     9 34HONS IN SATURATED NUMBER DENSITIES)/
     1 7H ITIXS ,I6,38H   (0/1=TIME-DEPENDENT XS FLAG ON/OFF))
  600 FORMAT(/54H EVODRV: SOLUTION OF A DEPLETION SYSTEM BETWEEN TIMES ,
     1 1P,E12.4,4H AND,E12.4,6H E+8 S)
      END
