path: root/Donjon/src/MCRLIB.f

*DECK MCRLIB
      SUBROUTINE MCRLIB(MAXNIS,MAXISO,IPLIB,IPMPO,IACCS,NMIX,NGRP,LADFM,
     1 IMPX,HEQUI,HMASL,NCAL,HEDIT,ITER,MY1,MY2,NBISO,TERP,NISO,LISO,
     2 HISO,CONC,ITODO,MIXC,LRES,LPURE,ILUPS,B2,VTOT,YLDS,DECAYC)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Build the Microlib by scanning the NCAL elementary calculations in
* a MPO file and weighting them with TERP factors.
*
*Copyright:
* Copyright (C) 2022 Ecole Polytechnique de Montreal
*
*Author(s): 
* A. Hebert
*
*Parameters: input
* MAXNIS  maximum value of NISO(I) in user data.
* MAXISO  maximum allocated space for output Microlib TOC information.
* IPLIB   address of the output Microlib LCM object.
* IPMPO   pointer to the MPO file.
* IACCS   =0 Microlib is created; =1 ... is updated.
* NMIX    maximum number of material mixtures in the Microlib.
* NGRP    number of energy groups.
* LADFM   type of discontinuity factors (.true.: diagonal; .false.: GxG).
* IMPX    print parameter (equal to zero for no print).
* HEQUI   keyword of SPH-factor set to be recovered.
* HMASL   keyword of MASL data set to be recovered.
* NCAL    number of elementary calculations in the MPO file.
* HEDIT   name of output group for a (multigroup mesh, output geometry)
*         couple (generally equal to 'output_0').
* ITER    completion flag (=0: compute the macrolib).
* MY1     number of fissile isotopes including macroscopic sets.
* MY2     number of fission fragment.
* NBISO   number of particularized isotopes in the MPO file.
* TERP    interpolation factors.
* NISO    number of user-selected isotopes.
* LISO    type of treatment (=.true.: ALL; =.false.: ONLY).
* HISO    name of the user-selected isotopes.
* CONC    user-defined number density of the user-selected isotopes. A
*         value of -99.99 is set to indicate that the MPO file value is
*         used.
* ITODO   non-depletion mask (=1 to force a user-selected isotope to be
*         non-depleting)
* MIXC    mixture index in the MPO file corresponding to each Microlib
*         mixture. Equal to zero if a Microlib mixture is not updated.
* LRES    =.true. if the interpolation is done without updating isotopic 
*         densities
* LPURE   =.true. if the interpolation is a pure linear interpolation 
*         with TERP factors.
* ILUPS   up-scattering removing flag (=1 to remove up-scattering from
*         output cross-sections).
* B2      buckling
* VTOT    volume of updated core.
* YLDS    fission yields.
* DECAYC  radioactive decay constants.
*
*-----------------------------------------------------------------------
*
      USE GANLIB
      USE hdf5_wrap
      IMPLICIT NONE
*----
*  SUBROUTINE ARGUMENTS
*----
      TYPE(C_PTR) IPLIB,IPMPO
      INTEGER MAXNIS,MAXISO,IACCS,NMIX,NGRP,IMPX,NCAL,ITER,MY1,MY2,
     > NBISO,NISO(NMIX),ITODO(NMIX,MAXNIS),MIXC(NMIX),ILUPS
      REAL TERP(NCAL,NMIX),CONC(NMIX,MAXNIS),B2
      DOUBLE PRECISION VTOT,YLDS(MY1,MY2),DECAYC(NBISO)
      LOGICAL LADFM,LISO(NMIX),LRES,LPURE
      CHARACTER(LEN=80) HEQUI,HMASL
      CHARACTER(LEN=12) HEDIT
      CHARACTER(LEN=8) HISO(NMIX,NBISO)
*----
*  LOCAL VARIABLES
*----
      INTEGER, PARAMETER::IOUT=6
      INTEGER, PARAMETER::MAXREA=50
      INTEGER, PARAMETER::NSTATE=40
      INTEGER, PARAMETER::MAXFRD=4
      TYPE(C_PTR) JPLIB,KPLIB
      REAL FACT0, WEIGHT, DEN
      INTEGER I, J, I0, IBM, IBMOLD, ICAL, IED2, IFISS, IGR, ILONG, IDF,
     > IPRC, IREA, IREAB, IREAF, ISO, ITRANC, ITSTMP, ITYLCM, IY1, IY2,
     > JSO, KSO, KSO1, LMY1, LSO, MAXMIX, NBISO2, NBISO2I, NBS1, NCALS,
     > NED2, NL, NMIL, NPAR, NPRC, NREA, NSURFD, NISOF, NISOP, NISOS,
     > RANK, NBYTE, TYPE, DIMSR(5), ILOC, NADDRXS, NLOC, NMGF, ID, ID_E,
     > ID_G, NENERG, NGEOME, ADDRZI, ISOM, NISOM, IGRC, NALBP, NALBP2
      CHARACTER RECNAM*80,RECNA2*80,TEXT8*8,TEXT12*12,HSMG*131,
     > HVECT2(MAXREA)*8,HRESID*8
      INTEGER ISTATE(NSTATE),INAME(2),IHRES(2)
      REAL TMPDAY(3)
      LOGICAL LUSER,LSTRD,LSPH,LMASL,LALBG
*----
*  ALLOCATABLE ARRAYS
*----
      INTEGER, ALLOCATABLE, DIMENSION(:) :: IMIX2,ISONA,ISOMI,ITOD2,
     > ISTY1,ISTY2,IPIFI,IMICR,ITOD1,JJSO,IPYMIX,LOCAD,REACTION,ISOTOPE,
     > ADDRISO,IGYELD,IADRY,DIMS_MPO
      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: HUSE2,HNAM2,OUPUTID
      INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: ADDRXS
      REAL, ALLOCATABLE, DIMENSION(:) :: DENS2,DENS3,VOL2,VOLMI2,SPH,
     > ENER,VOSAP,CONCE,TAUXFI,NWT0,FLUXS,DENIS,GAR1,GAR2,LAMB,BETAR,
     > INVELS,BETARB,INVELSB,DECAY2,RVALO
      REAL, ALLOCATABLE, DIMENSION(:,:) :: DENS1,FACT,CHIRS,CHIRSB,
     > TAUXGF
      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: XS,SIGS,DENS0,FLUX,YLDS2
      REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: SS2D
      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: YLDSM
      LOGICAL, ALLOCATABLE, DIMENSION(:) :: LXS,MASK,MASKL
      CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: HPYNAM
      CHARACTER(LEN=24), ALLOCATABLE, DIMENSION(:) :: TEXT24,NOMREA,
     > NOMISO
      CHARACTER(LEN=80), ALLOCATABLE, DIMENSION(:) :: LOCTYP,LOCKEY
*----
*  RECOVER MPO FILE CHARACTERISTICS
*----
      I=0
      CALL MPOTOC(IPMPO,HEDIT,0,NREA,I0,NMIL,NPAR,NLOC,NISOF,NISOP,
     > NISOS,NCALS,I,NSURFD,NALBP,NPRC)
      IF(NBISO.NE.I0) CALL XABORT('MCRLIB: INVALID VALUE OF NBISO.')
      IF(NGRP.NE.I) CALL XABORT('MCRLIB: INVALID VALUE OF NGRP.')
      IF(NREA+2.GT.MAXREA) CALL XABORT('MCRLIB: MAXREA OVERFLOW')
*----
*  SCRATCH STORAGE ALLOCATION
*----
      ALLOCATE(IMIX2(MAXISO),ITOD2(MAXISO),ISTY1(MAXISO),ISTY2(MAXISO),
     > HUSE2(3,MAXISO),HNAM2(3,MAXISO))
      ALLOCATE(DENS2(MAXISO),DENS3(MAXISO),VOL2(MAXISO),VOLMI2(NMIX),
     > FLUX(NMIX,NGRP,2),SPH(NGRP))
*----
*  MICROLIB INITIALIZATION
*----
      VOLMI2(:NMIX)=0.0
      DENS2(:MAXISO)=0.0
      VOL2(:MAXISO)=0.0
      IMIX2(:MAXISO)=0
      ITOD2(:MAXISO)=0
      ISTY2(:MAXISO)=0
      IF(IACCS.EQ.0) THEN
        IF(LRES) CALL XABORT('MCRLIB: RES OPTION IS INVALID.')
        NBISO2=0
        NED2=0
        TEXT12='L_LIBRARY'
        CALL LCMPTC(IPLIB,'SIGNATURE',12,TEXT12)
      ELSE
        CALL LCMGET(IPLIB,'STATE-VECTOR',ISTATE)
        IF(ISTATE(1).NE.NMIX) CALL XABORT('MCRLIB: INVALID NUMBER OF '
     1  //'MATERIAL MIXTURES IN THE MICROLIB.')
        IF(ISTATE(3).NE.NGRP) CALL XABORT('MCRLIB: INVALID NUMBER OF '
     1  //'ENERGY GROUPS IN THE MICROLIB.')
        NBISO2=ISTATE(2)
        IF(NBISO2.GT.MAXISO) CALL XABORT('MCRLIB: MAXISO OVERFLOW(1).')
        NED2=ISTATE(13)
        IF(NED2.GT.MAXREA) CALL XABORT('MCRLIB: MAXREA OVERFLOW.')
        CALL LCMLEN(IPLIB,'MIXTURESVOL',ILONG,ITYLCM)
        IF(ILONG.GT.0) THEN
          CALL LCMGET(IPLIB,'MIXTURESVOL',VOLMI2)
        ELSE
          VOLMI2(:NMIX)=0.0
        ENDIF
        CALL LCMGET(IPLIB,'ISOTOPESUSED',HUSE2)
        CALL LCMGET(IPLIB,'ISOTOPERNAME',HNAM2)
        CALL LCMGET(IPLIB,'ISOTOPESDENS',DENS2)
        CALL LCMGET(IPLIB,'ISOTOPESVOL',VOL2)
        CALL LCMGET(IPLIB,'ISOTOPESMIX',IMIX2)
        CALL LCMGET(IPLIB,'ISOTOPESTODO',ITOD2)
        CALL LCMGET(IPLIB,'ISOTOPESTYPE',ISTY2)
        IF(NED2.GT.0) CALL LCMGTC(IPLIB,'ADDXSNAME-P0',8,NED2,HVECT2)
      ENDIF
*----
*  SET EQUIVALENCE AND HEAVY DENSITY FLAGS.
*----
      LSPH=.FALSE.
      LMASL=.FALSE.
      NLOC=0
      IF(hdf5_group_exists(IPMPO,"/local_values/")) THEN
        CALL hdf5_read_data(IPMPO,"/local_values/LOCVALTYPE",LOCTYP)
        CALL hdf5_read_data(IPMPO,"/local_values/LOCVALNAME",LOCKEY)
        NLOC=SIZE(LOCTYP,1)
        DO ILOC=1,NLOC
          LSPH=LSPH.OR.((LOCTYP(ILOC).EQ.'EQUI').AND.
     >                  (LOCKEY(ILOC).EQ.HEQUI))
          LMASL=LMASL.OR.((LOCTYP(ILOC).EQ.'HEAVY_METAL_DENSITY').AND.
     >                    (LOCKEY(ILOC).EQ.HMASL))
        ENDDO
      ENDIF
*----
*  FIND SCATTERING ANISOTROPY.
*----
      WRITE(RECNAM,'(8H/output/,A,6H/info/)') TRIM(HEDIT)
      CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"NADDRXS",NADDRXS)
      CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ADDRXS",ADDRXS)
      NL=0
      DO I=1,NADDRXS
        DO ISO=1,NBISO
          NL=MAX(NL,ADDRXS(NREA+1,ISO,I))
          NL=MAX(NL,ADDRXS(NREA+2,ISO,I))
        ENDDO
      ENDDO
      IF(IMPX.GT.1) THEN
        WRITE(6,'(37H MCRLIB: number of legendre orders  =,I4)') NL
      ENDIF
*----
*  SET ENERGY MESH AND ZONE VOLUMES
*----
      CALL hdf5_read_data(IPMPO,"/energymesh/NENERGYMESH",NENERG)
      CALL hdf5_read_data(IPMPO,"/geometry/NGEOMETRY",NGEOME)
      CALL hdf5_read_data(IPMPO,"/output/OUPUTID",OUPUTID)
      READ(HEDIT,'(7X,I2)') ID
      ID_G=0
      ID_E=0
      DO I=1,NGEOME
        DO J=1,NENERG
          IF(OUPUTID(J,I).EQ.ID) THEN
            ID_G=I-1
            ID_E=J-1
            GO TO 10
          ENDIF
        ENDDO
      ENDDO
      CALL XABORT('MCRLIB: no ID found in /output/OUPUTID.')
   10 WRITE(RECNAM,'(23H/energymesh/energymesh_,I0)') ID_E
      CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"/ENERGY",ENER)
      IF(SIZE(ENER,1)-1.NE.NGRP) CALL XABORT('MCRLIB: INVALID NGRP VAL'
     > //'UE.')
      DO IGR=1,NGRP+1
        ENER(IGR)=ENER(IGR)/1.0E-6
      ENDDO
      WRITE(RECNAM,'(19H/geometry/geometry_,I0,1H/)') ID_G
      CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ZONEVOLUME",VOSAP)
      CALL LCMPUT(IPLIB,'ENERGY',NGRP+1,2,ENER)
      DO IGR=1,NGRP
        ENER(IGR)=LOG(ENER(IGR)/ENER(IGR+1))
      ENDDO
      CALL LCMPUT(IPLIB,'DELTAU',NGRP,2,ENER)
      DEALLOCATE(ENER)
*----
*  RECOVER INFORMATION ON REACTIONS AND ISOTOPE NAMES
*----
      IREAB=0
      IREAF=0
      WRITE(RECNAM,'(8H/output/,A,6H/info/)') TRIM(HEDIT)
      CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"REACTION",REACTION)
      CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ISOTOPE",ISOTOPE)
      CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ADDRISO",ADDRISO)
      NBISO=ADDRISO(SIZE(ADDRISO,1))
      IF(NBISO.EQ.0) CALL XABORT('MCRLIB: NO CROSS SECTIONS.')
      ALLOCATE(NOMREA(NREA+2),NOMISO(NBISO))
      IF(NREA.GT.0) THEN
        CALL hdf5_read_data(IPMPO,"/contents/reactions/REACTIONAME",
     >  TEXT24)
        DO I=1,NREA
          NOMREA(I)=TEXT24(REACTION(I)+1)
        ENDDO
        DEALLOCATE(TEXT24,REACTION)
        DO IREA=1,NREA
          IF(NOMREA(IREA).EQ.'Absorption') THEN
            IREAB=IREA
            EXIT
          ENDIF
        ENDDO
        DO IREA=1,NREA
          IF(NOMREA(IREA).EQ.'NuFission') THEN
            IREAF=IREA
            EXIT
          ENDIF
        ENDDO
      ENDIF
      NOMREA(NREA+1)='Total'
      NOMREA(NREA+2)='Leakage'
      NREA=NREA+2
      CALL hdf5_read_data(IPMPO,"/contents/isotopes/ISOTOPENAME",TEXT24)
      DO I=1,NBISO
        NOMISO(I)=TEXT24(ISOTOPE(I)+1)
      ENDDO
      DEALLOCATE(TEXT24,ADDRISO,ISOTOPE)
      IF(IMPX.GT.1) THEN
        WRITE(6,'(/24H MCRLIB: reaction names:)')
        DO I=1,NREA
          WRITE(6,'(5X,7HNOMREA(,I3,2H)=,A)') I,TRIM(NOMREA(I))
        ENDDO
        WRITE(6,'(/23H MCRLIB: isotope names:)')
        DO I=1,NBISO
          WRITE(6,'(5X,7HNOMISO(,I3,2H)=,A)') I,TRIM(NOMISO(I))
        ENDDO
      ENDIF
*----
*  LOOP OVER MPO MIXTURES TO COMPUTE DENS0(NMIL,NCAL,NBISO)
*----
      ALLOCATE(DENS0(NMIL,NCAL,NBISO))
      DENS0(:NMIL,:NCAL,:NBISO)=0.0
      DO 30 IBMOLD=1,NMIL
      DO ICAL=1,NCAL
        DO IBM=1,NMIX
          IF((TERP(ICAL,IBM).NE.0.0).AND.(MIXC(IBM).EQ.IBMOLD)) GO TO 15
        ENDDO
        CYCLE
   15   WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,6H/zone_,I0,1H/)')
     >  TRIM(HEDIT),ICAL-1,IBMOLD-1
        IF(NBISO.GT.0) THEN
          CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"CONCENTRATION",CONCE)
          DO 20 ISO=1,NBISO
            DENS0(IBMOLD,ICAL,ISO)=CONCE(ISO)
   20     CONTINUE
          DEALLOCATE(CONCE)
        ENDIF
      ENDDO
   30 CONTINUE
*----
*  LOOP OVER MICROLIB MIXTURES
*----
      YLDS(:MY1,:MY2)=0.0D0
      DECAYC(:NBISO)=0.0D0
      VTOT=0.0D0
      DO 40 IBM=1,NMIX
      IBMOLD=MIXC(IBM)
      IF(IBMOLD.NE.0) VTOT=VTOT+VOSAP(IBMOLD)
   40 CONTINUE
      ALLOCATE(JJSO(NBISO),YLDSM(MY1,MY2),ITOD1(NBISO))
      ALLOCATE(TAUXFI(NBISO),NWT0(NGRP),SIGS(NGRP,NL,NBISO),
     > SS2D(NGRP,NGRP,NL,NBISO),XS(NGRP,NREA,NBISO))
      ALLOCATE(LXS(NREA))
      ALLOCATE(CHIRS(NGRP,NPRC),BETAR(NPRC),INVELS(NGRP))
      CHIRS(:NGRP,:NPRC)=0.0
      BETAR(:NPRC)=0.0
      INVELS(:NGRP)=0.0
      ALLOCATE(BETARB(NPRC),INVELSB(NGRP))
      ALLOCATE(DENS1(NBISO,NCAL),FACT(NBISO,NCAL))
      JPLIB=LCMLID(IPLIB,'ISOTOPESLIST',NBISO*NMIX)
*
      DO 180 IBM=1,NMIX
      IBMOLD=MIXC(IBM)
      IF(IBMOLD.EQ.0) GO TO 180
      IF(NISO(IBM).GT.MAXNIS) CALL XABORT('MCRLIB: MAXNIS OVERFLOW.')
      VOLMI2(IBM)=VOSAP(IBMOLD)
*----
*  RECOVER ITOD1(NBISO) INDICES.
*----
      ITOD1(:NBISO)=0
      DO 50 ISO=1,NBISO ! MPO file isotope
        DO KSO=1,NISO(IBM) ! user-selected isotope
          IF(NOMISO(ISO).EQ.HISO(IBM,KSO)) THEN
            ITOD1(ISO)=ITODO(IBM,KSO)
            GO TO 50
          ENDIF
        ENDDO
   50 CONTINUE
*----
*  COMPUTE THE NUMBER DENSITIES OF EACH ELEMENTARY CALCULATION.
*----
      DENS1(:NBISO,:NCAL)=0.0
      DENS3(:NBISO)=0.0
      DO ICAL=1,NCAL
        WEIGHT=TERP(ICAL,IBM)
        IF(WEIGHT.EQ.0.0) CYCLE
        DO ISO=1,NBISO
          LUSER=.FALSE.
          KSO1=0
          DO KSO=1,NISO(IBM) ! user-selected isotope
            IF(NOMISO(ISO).EQ.HISO(IBM,KSO)) THEN
              KSO1=KSO
              LUSER=(CONC(IBM,KSO1).NE.-99.99)
              GO TO 60
            ENDIF
          ENDDO
   60     IF(LUSER) THEN
            DENS1(ISO,ICAL)=CONC(IBM,KSO1)
            CYCLE
          ENDIF
          IF(.NOT.LISO(IBM)) CYCLE
          DENS1(ISO,ICAL)=DENS0(IBMOLD,ICAL,ISO)
        ENDDO
        DO ISO=1,NBISO
          DENS3(ISO)=DENS3(ISO)+WEIGHT*DENS1(ISO,ICAL)
        ENDDO
      ENDDO
      FACT(:NBISO,:NCAL)=1.0
      IF(.NOT.LPURE) THEN
        DO ICAL=1,NCAL
          IF(TERP(ICAL,IBM).EQ.0.0) CYCLE
          DO ISO=1,NBISO
            IF(DENS3(ISO).GT.DENS1(ISO,ICAL)*1.0E-9) THEN
              FACT(ISO,ICAL)=DENS1(ISO,ICAL)/DENS3(ISO)
            ENDIF
          ENDDO
        ENDDO
      ENDIF
*----
*  INITIALIZE WORKING ARRAYS.
*----
      TAUXFI(:NBISO)=0.0
      NWT0(:NGRP)=0.0
      SIGS(:NGRP,:NL,:NBISO)=0.0
      SS2D(:NGRP,:NGRP,:NL,:NBISO)=0.0
      XS(:NGRP,:NREA,:NBISO)=0.0
      LXS(:NREA)=.FALSE.
      YLDSM(:MY1,:MY2)=0.0D0
*----
*  MAIN LOOP OVER ELEMENTARY CALCULATIONS
*----
      TEXT12='*MAC*RES'
      READ(TEXT12,'(2A4)') IHRES(1),IHRES(2)
      LSTRD=.FALSE.
      DO 80 ICAL=1,NCAL
      WEIGHT=TERP(ICAL,IBM)
      IF(WEIGHT.EQ.0.0) GO TO 80
*----
*  SELECT THE HDF5 GROUP CORRESPONDING TO ICAL AND IBMOLD
*----
      WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,6H/zone_,I0,1H/)')
     > TRIM(HEDIT),ICAL-1,IBMOLD-1
      NMGF=0
      IF(hdf5_group_exists(IPMPO,TRIM(RECNAM)//"yields")) THEN
        CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"yields/NMGF",NMGF)
        CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"yields/YIELDGROUP",
     >  IGYELD)
      ENDIF
      ALLOCATE(TAUXGF(NMGF,NBISO))
*----
*  RECOVER INFORMATION FROM caldir GROUP.
*----
      WRITE(RECNA2,'(A,9Hkinetics/)') TRIM(RECNAM)
      CALL hdf5_info(IPMPO,TRIM(RECNA2)//"LAMBDA",RANK,TYPE,NBYTE,DIMSR)
      NPRC=0
      IF(TYPE.NE.99) THEN
        NPRC=DIMSR(1)
        CALL hdf5_read_data(IPMPO,TRIM(RECNA2)//"LAMBDA",LAMB)
        CALL hdf5_read_data(IPMPO,TRIM(RECNA2)//"CHIDA",CHIRSB)
        CALL hdf5_read_data(IPMPO,TRIM(RECNA2)//"BETADA",BETARB)
        CALL hdf5_read_data(IPMPO,TRIM(RECNA2)//"INVELA",INVELSB)
      ENDIF
*----
*  RECOVER SPH FACTORS.
*----
      SPH(:NGRP)=1.0
      IF(HEQUI.NE.' ') THEN
        CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"LOCALVALUE",RVALO)
        CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"LOCALVALADDR",LOCAD)
        DO ILOC=1,NLOC
          IF((LOCTYP(ILOC).EQ.'EQUI').AND.(LOCKEY(ILOC).EQ.HEQUI))
     >    THEN
            IF(LOCAD(ILOC+1)-LOCAD(ILOC).NE.NGRP) THEN
              CALL XABORT('MCRLIB: INVALID NUMBER OF COMPONENTS FOR '
     >        //'SPH FACTORS')
            ENDIF
            DO IGR=1,NGRP
              SPH(IGR)=RVALO(LOCAD(ILOC)+IGR-1)
            ENDDO
          ENDIF
        ENDDO
        DEALLOCATE(LOCAD,RVALO)
      ENDIF
*----
*  RECOVER FLUXES.
*----
      CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ZONEFLUX",FLUXS)
      DO I=1,NGRP
        NWT0(I)=NWT0(I)+WEIGHT*FLUXS(I)/SPH(I)
      ENDDO
*----
*  RECOVER MICROSCOPIC CROSS SECTIONS.
*----
      DO ISO=1,NBISO
        FACT0=FACT(ISO,ICAL)
        DEN=DENS0(IBMOLD,ICAL,ISO)
        CALL MCRSX2(IPMPO,HEDIT,RECNAM,NREA,NGRP,NMGF,NL,ISO,NOMREA,
     1  NOMISO(ISO),DEN,FACT0,WEIGHT,SPH,FLUXS,IREAB,IREAF,LPURE,IGYELD,
     2  LXS,XS(1,1,ISO),SIGS(1,1,ISO),SS2D(1,1,1,ISO),TAUXFI(ISO),
     3  TAUXGF(1,ISO))
      ENDDO
      IF(NMGF.GT.0) DEALLOCATE(IGYELD)
      DEALLOCATE(FLUXS)
*
      IF(NPRC.GT.0) THEN
        DO IGR=1,NGRP
          INVELS(IGR)=INVELS(IGR)+SPH(IGR)*WEIGHT*INVELSB(IGR)
          DO IPRC=1,NPRC
            CHIRS(IGR,IPRC)=CHIRS(IGR,IPRC)+WEIGHT*CHIRSB(IGR,IPRC)
          ENDDO
        ENDDO
        DO IPRC=1,NPRC
          BETAR(IPRC)=BETAR(IPRC)+WEIGHT*BETARB(IPRC)
        ENDDO
      ENDIF
*----
*  COMPUTE DEPLETION CHAIN DATA
*----
      IF(NISOF*NISOP.GT.0) THEN
        IF(NMGF.EQ.0) CALL XABORT('MCRLIB: INVALID NMGF.')
        WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,6H/zone_,I0,1H/)')
     >  TRIM(HEDIT),ICAL-1,IBMOLD-1
        CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"yields/NISF",I0)
        IF(I0.NE.NISOF) CALL XABORT('MCRLIB: INVALID NISOF.')
        CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"yields/NISP",I0)
        IF(I0.NE.NISOP) CALL XABORT('MCRLIB: INVALID NISOP.')
        CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"yields/YIELD",YLDS2)
        CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"yields/ADDRY",
     >  DIMS_MPO)
        CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ADDRZI",ADDRZI)
        WRITE(RECNAM,'(8H/output/,A,6H/info/)') TRIM(HEDIT)
        CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ADDRISO",ADDRISO)
        NISOM=ADDRISO(ADDRZI+2)-ADDRISO(ADDRZI+1)
        DO IY1=1,NISOF
          ISO=0
          DO ISOM=1,NISOM
            IF(DIMS_MPO(ISOM).EQ.IY1) THEN
              ISO=ADDRISO(ADDRZI+1)+ISOM
              EXIT
            ENDIF
          ENDDO
          IF(ISO.EQ.0) CALL XABORT('MCRLIB: UNABLE TO FIND ISO.')
          DEN=0.0
          DO IGRC=1,NMGF
            DEN=DEN+TAUXGF(IGRC,ISO)
            DO IY2=1,NISOP
              YLDSM(IY1,IY2)=YLDSM(IY1,IY2)+WEIGHT*TAUXGF(IGRC,ISO)*
     >        YLDS2(IY1,IY2,IGRC)
              YLDS(IY1,IY2)=YLDS(IY1,IY2)+WEIGHT*TAUXGF(IGRC,ISO)*
     >        YLDS2(IY1,IY2,IGRC)*VOLMI2(IBM)/VTOT
            ENDDO
          ENDDO
          IF(DEN.EQ.0.0) CYCLE
          DO IY2=1,NISOP
            YLDSM(IY1,IY2)=YLDSM(IY1,IY2)/DEN
            YLDS(IY1,IY2)=YLDS(IY1,IY2)/DEN
          ENDDO
        ENDDO
        DEALLOCATE(ADDRISO,DIMS_MPO,YLDS2)
      ENDIF
      CALL hdf5_info(IPMPO,"/contents/isotopes/DECAYCONST",RANK,TYPE,
     1 NBYTE,DIMSR)
      IF(TYPE.NE.99) THEN
        CALL hdf5_read_data(IPMPO,"/contents/isotopes/DECAYCONST",
     >  DECAY2)
        DO ISO=1,NBISO
          DECAYC(ISO)=DECAYC(ISO)+WEIGHT*DECAY2(ISO)*VOLMI2(IBM)/VTOT
        ENDDO
        DEALLOCATE(DECAY2)
      ENDIF
      DEALLOCATE(TAUXGF)
   80 CONTINUE ! end of loop over elementary calculations.
*----
*  IDENTIFY SPECIAL FLUX EDITS
*----
      DO 100 IREA=1,NREA
        IF((NOMREA(IREA).EQ.'Total').or.
     &     (NOMREA(IREA).EQ.'Absorption').or.
     &     (NOMREA(IREA).EQ.'CaptureEnergyCapture').or.
     &     (NOMREA(IREA).EQ.'Diffusion').or.
     &     (NOMREA(IREA).EQ.'FissionEnergyFission').or.
     &     (NOMREA(IREA).EQ.'FissionSpectrum').or.
     &     (NOMREA(IREA).EQ.'NuFission').or.
     &     (NOMREA(IREA).EQ.'Scattering')) CYCLE
        DO 90 IED2=1,NED2
        IF(HVECT2(IED2).EQ.NOMREA(IREA)(:8)) GO TO 100
        IF(HVECT2(IED2).EQ.'NFTOT') GO TO 100
   90   CONTINUE
        NED2=NED2+1
        IF(NED2.GT.MAXREA) CALL XABORT('MCRLIB: MAXREA OVERFLOW.')
        IF(NOMREA(IREA).EQ.'Fission') THEN
          HVECT2(NED2)='NFTOT'
        ELSE
          HVECT2(NED2)=NOMREA(IREA)(:8)
        ENDIF
  100 CONTINUE
*----
*  SET FLAG LSTRD
*----
      LSTRD=.TRUE.
      DO IREA=1,NREA
        IF(NOMREA(IREA).EQ.'Leakage') THEN
          IF(LXS(IREA)) LSTRD=.FALSE.
          EXIT
        ENDIF
      ENDDO
*----
*  SET IADRY FOR MIXTURE IBMOLD
*----
      ALLOCATE(IADRY(NBISO))
      IADRY(:NBISO)=0
      DO ICAL=NCAL,1,-1
      IF(TERP(ICAL,IBM).EQ.0.0) CYCLE
        WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,6H/zone_,I0,1H/)')
     1  TRIM(HEDIT),ICAL-1,IBMOLD-1
        IF((hdf5_group_exists(IPMPO,TRIM(RECNAM)//"yields")).AND.
     1    (NISOP.GT.0)) THEN
          CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"yields/ADDRY",
     1    DIMS_MPO)
          CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ADDRZI",ADDRZI)
          WRITE(RECNAM,'(8H/output/,A,6H/info/)') TRIM(HEDIT)
          CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ADDRISO",ADDRISO)
          NISOM=ADDRISO(ADDRZI+2)-ADDRISO(ADDRZI+1)
          DO ISOM=1,NISOM
            ISO=ADDRISO(ADDRZI+1)+ISOM
            IADRY(ISO)=DIMS_MPO(ISOM)
          ENDDO
          DEALLOCATE(ADDRISO,DIMS_MPO)
        ENDIF
        EXIT
      ENDDO
*----
*  SAVE CROSS SECTIONS IN MICROLIB FOR MIXTURE IBM
*----
      ISTY1(:NBISO)=0
      JJSO(:NBISO)=0
      NBISO2I=NBISO2
      HRESID=' '
      DO ISO=1,NBISO
        READ(NOMISO(ISO),'(2A4)') INAME(:2)
        CALL SCRFND(MAXISO,NBISO2I,NBISO2,INAME,IBM,HRESID,HUSE2,
     1  HNAM2,IMIX2,JJSO(ISO))
        KPLIB=LCMDIL(JPLIB,JJSO(ISO)) ! step up isot JJSO(ISO)
        CALL MCRISO(KPLIB,NREA,NGRP,NL,NPRC,NOMREA,NWT0,XS(1,1,ISO),
     1  SIGS(1,1,ISO),SS2D(1,1,1,ISO),TAUXFI(ISO),LXS,LAMB,CHIRS,
     2  BETAR,INVELS,INAME,LSTRD,LPURE,ILUPS,ITRANC,IFISS)
        IF(MY1*MY2.GT.0) CALL MCRNDF(IMPX,NBISO,ISO,IBM,NOMISO,KPLIB,
     1  MY1,MY2,YLDSM,IADRY,ISTY1(ISO))
      ENDDO
      DEALLOCATE(IADRY)
*----
*  SET NUMBER DENSITIES AND VOLUMES IN OUTPUT MICROLIB
*----
      IF(LRES) THEN
*       -- Number densities are left unchanged except if they are
*       -- listed in HISO array.
        DO 110 KSO=1,NISO(IBM) ! user-selected isotope
          DO JSO=1,NBISO2 ! microlib isotope
            IF(IMIX2(JSO).NE.IBM) CYCLE
            WRITE(TEXT8,'(2A4)') HUSE2(1,JSO),HUSE2(2,JSO)
            IF(HISO(IBM,KSO).EQ.TEXT8) THEN
              ITOD2(JSO)=ITODO(IBM,KSO)
              IF(CONC(IBM,KSO).EQ.-99.99) THEN
*               -- Only number densities of isotopes set with "MICR" and
*               -- "*" keywords are interpolated
                DENS2(JSO)=0.0
                DO ISO=1,NBISO ! MPO file isotope
                  IF(JJSO(ISO).EQ.JSO) DENS2(JSO)=DENS2(JSO)+DENS3(ISO)
                ENDDO
              ELSE IF(CONC(IBM,KSO).NE.-99.99) THEN
*               -- Number densities of isotopes set with "MICR" and
*               -- fixed value are forced to this value
                DENS2(JSO)=CONC(IBM,KSO)
              ENDIF
              GO TO 110
            ENDIF
          ENDDO
          WRITE(HSMG,'(31HMCRLIB: UNABLE TO FIND ISOTOPE ,A8,6H IN MI,
     1    5HXTURE,I8,1H.)') HISO(IBM,KSO),IBM
          CALL XABORT(HSMG)
  110   CONTINUE
      ELSE
*       -- Number densities are interpolated or not according to
*       -- ALL/ONLY option
        DO JSO=1,NBISO2 ! microlib isotope
          WRITE(TEXT8,'(2A4)') HUSE2(1,JSO),HUSE2(2,JSO)
          IF(IBM.EQ.IMIX2(JSO)) THEN
            DO ISO=1,NBISO ! MPO file isotope
              IF(NOMISO(ISO).EQ.TEXT8) THEN
                DENS2(JSO)=0.0
                VOL2(JSO)=0.0
                CYCLE
              ENDIF
            ENDDO
          ENDIF
        ENDDO
        DO 130 ISO=1,NBISO ! MPO file isotope
        IF(.NOT.LISO(IBM)) THEN
*         --ONLY option
          DO KSO=1,NISO(IBM) ! user-selected isotope
            IF(NOMISO(ISO).EQ.HISO(IBM,KSO)) GO TO 120
          ENDDO
          GO TO 130
        ENDIF
  120   JSO=JJSO(ISO)
        IF(JSO.GT.0) THEN
          ITOD2(JSO)=ITOD1(ISO)
          ISTY2(JSO)=ISTY1(ISO)
          DENS2(JSO)=DENS2(JSO)+DENS3(ISO)
          VOL2(JSO)=VOL2(JSO)+VOSAP(IBMOLD)
        ENDIF
  130   CONTINUE
      ENDIF
*----
*  SET PIFI INFORMATION
*----
      ALLOCATE(IMICR(NBISO))
      IMICR(:NBISO)=0
      NBS1=0
      DO 140 JSO=1,NBISO2 ! microlib isotope
      IF(IMIX2(JSO).EQ.IBM) THEN
        NBS1=NBS1+1
        IF(NBS1.GT.NBISO) CALL XABORT('MCRLIB: NBISO OVERFLOW.')
        IMICR(NBS1)=JSO
      ENDIF
  140 CONTINUE
      DO 170 ISO=1,NBS1 ! MPO file isotope
      JSO=IMICR(ISO)
      KPLIB=LCMDIL(JPLIB,JSO) ! step up isot JSO
      CALL LCMLEN(KPLIB,'PYIELD',LMY1,ITYLCM)
      IF(LMY1.GT.0) THEN
        ALLOCATE(HPYNAM(LMY1),IPYMIX(LMY1),IPIFI(LMY1))
        IPIFI(:LMY1)=0
        CALL LCMGTC(KPLIB,'PYNAM',8,LMY1,HPYNAM)
        CALL LCMGET(KPLIB,'PYMIX',IPYMIX)
        DO 160 IY1=1,LMY1
        IF(HPYNAM(IY1).NE.' ') THEN
          DO 150 KSO=1,NBS1
          LSO=IMICR(KSO)
          WRITE(TEXT8,'(2A4)') HUSE2(:2,LSO)
          IF((HPYNAM(IY1).EQ.TEXT8).AND.(IPYMIX(IY1).EQ.IMIX2(LSO)))
     >    THEN
            IPIFI(IY1)=LSO
            GO TO 160
          ENDIF
  150     CONTINUE
          IF(IPIFI(IY1).EQ.0) THEN
            WRITE(HSMG,'(40HMCRLIB: FAILURE TO FIND FISSILE ISOTOPE ,
     1      A12,25H AMONG MICROLIB ISOTOPES.)') HPYNAM(IY1)
            CALL XABORT(HSMG)
          ENDIF
        ENDIF
  160   CONTINUE
        CALL LCMPUT(KPLIB,'PIFI',LMY1,1,IPIFI)
        DEALLOCATE(IPIFI,IPYMIX,HPYNAM)
      ENDIF
  170 CONTINUE
      DEALLOCATE(IMICR)
  180 CONTINUE ! end of loop over microlib mixtures.
*----
*  RELEASE MEMORY
*----
      DEALLOCATE(FACT,DENS1)
      IF(NPRC.GT.0) DEALLOCATE(INVELSB,BETARB,CHIRSB,INVELS,BETAR,
     > CHIRS,LAMB)
      DEALLOCATE(LXS,XS,SS2D,SIGS,NWT0,TAUXFI)
      DEALLOCATE(ITOD1,YLDSM)
      DEALLOCATE(JJSO,DENS0,NOMISO)
*----
*  MICROLIB FINALIZATION
*----
      IF(.NOT.LRES) THEN
        ISTATE(:NSTATE)=0
        ISTATE(1)=NMIX
        ISTATE(2)=NBISO2
        ISTATE(3)=NGRP
        ISTATE(4)=NL
        ISTATE(5)=ITRANC
        ISTATE(7)=1
        IF(ITER.EQ.3) ISTATE(12)=NMIX
        ISTATE(13)=NED2
        ISTATE(14)=NMIX
        ISTATE(18)=1
        ISTATE(19)=NPRC
        ISTATE(20)=MY1
        ISTATE(22)=MAXISO/NMIX
        IF(NBISO2.EQ.0) CALL XABORT('MCRLIB: NBISO2=0.')
        CALL LCMPUT(IPLIB,'STATE-VECTOR',NSTATE,1,ISTATE)
        CALL LCMPUT(IPLIB,'MIXTURESVOL',NMIX,2,VOLMI2)
        CALL LCMPUT(IPLIB,'ISOTOPESUSED',3*NBISO2,3,HUSE2)
        CALL LCMPUT(IPLIB,'ISOTOPERNAME',3*NBISO2,3,HNAM2)
        CALL LCMPUT(IPLIB,'ISOTOPESDENS',NBISO2,2,DENS2)
        CALL LCMPUT(IPLIB,'ISOTOPESMIX',NBISO2,1,IMIX2)
        CALL LCMPUT(IPLIB,'ISOTOPESVOL',NBISO2,2,VOL2)
        IF(NED2.GT.0) CALL LCMPTC(IPLIB,'ADDXSNAME-P0',8,NED2,HVECT2)
        CALL LCMPUT(IPLIB,'ISOTOPESTODO',NBISO2,1,ITOD2)
        CALL LCMPUT(IPLIB,'ISOTOPESTYPE',NBISO2,1,ISTY2)
      ELSE IF(LRES.AND.(NBISO.GT.0)) THEN
        CALL LCMPUT(IPLIB,'ISOTOPESDENS',NBISO2,2,DENS2)
        CALL LCMPUT(IPLIB,'ISOTOPESVOL',NBISO2,2,VOL2)
      ENDIF
      IF(IMPX.GT.5) CALL LCMLIB(IPLIB)
*----
*  COMPUTE THE MACROSCOPIC X-SECTIONS
*----
      IF((ITER.NE.0).AND.(ITER.NE.3)) GO TO 280
      CALL LCMGET(IPLIB,'STATE-VECTOR',ISTATE)
      MAXMIX=ISTATE(1)
      IF(MAXMIX.NE.NMIX) CALL XABORT('MCRLIB: INVALID NMIX.')
      NBISO=ISTATE(2)
      ALLOCATE(MASK(MAXMIX),MASKL(NGRP))
      ALLOCATE(ISONA(3*NBISO),ISOMI(NBISO),DENIS(NBISO))
      CALL LCMGET(IPLIB,'ISOTOPESUSED',ISONA)
      CALL LCMGET(IPLIB,'ISOTOPESMIX',ISOMI)
      CALL LCMGET(IPLIB,'ISOTOPESDENS',DENIS)
      MASK(:MAXMIX)=.TRUE.
      MASKL(:NGRP)=.TRUE.
      ITSTMP=0
      TMPDAY(1)=0.0
      TMPDAY(2)=0.0
      TMPDAY(3)=0.0
      CALL LCMLEN(IPLIB,'MACROLIB',ILONG,ITYLCM)
      IF(ILONG.NE.0) CALL LCMDEL(IPLIB,'MACROLIB')
      CALL LIBMIX(IPLIB,MAXMIX,NGRP,NBISO,ISONA,ISOMI,DENIS,MASK,MASKL,
     > ITSTMP,TMPDAY)
      DEALLOCATE(MASKL,MASK)
      DEALLOCATE(DENIS,ISOMI,ISONA)
*----
*  INCLUDE LEAKAGE IN THE MACROLIB (USED ONLY FOR NON-REGRESSION TESTS)
*----
      IF(B2.NE.0.0) THEN
        IF(IMPX.GT.0) WRITE(IOUT,'(/31H MCRLIB: INCLUDE LEAKAGE IN THE,
     >  14H MACROLIB (B2=,1P,E12.5,2H).)') B2
        CALL LCMSIX(IPLIB,'MACROLIB',1)
        JPLIB=LCMGID(IPLIB,'GROUP')
        ALLOCATE(GAR1(NMIX),GAR2(NMIX))
        DO 270 IGR=1,NGRP
          KPLIB=LCMGIL(JPLIB,IGR)
          CALL LCMGET(KPLIB,'NTOT0',GAR1)
          CALL LCMGET(KPLIB,'DIFF',GAR2)
          DO 260 IBM=1,NMIX
            IF(MIXC(IBM).NE.0) GAR1(IBM)=GAR1(IBM)+B2*GAR2(IBM)
  260     CONTINUE
          CALL LCMPUT(KPLIB,'NTOT0',NMIX,2,GAR1)
  270   CONTINUE
        DEALLOCATE(GAR2,GAR1)
        CALL LCMSIX(IPLIB,' ',2)
      ENDIF
*----
*  PROCESS ADF and physical albedos (if required)
*----
  280 LALBG=.TRUE.
      IDF=0
      IF(NALBP.GT.0) THEN
        WRITE(RECNAM,'(8H/output/,A,16H/statept_0/flux/)') TRIM(HEDIT)
        CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"NALBP",NALBP2)
        IF(NALBP2.NE.NALBP) CALL XABORT('MCRLIB: INVALID NALBP.')
        CALL hdf5_info(IPMPO,TRIM(RECNAM)//"ALBEDOGxG",RANK,TYPE,NBYTE,
     &  DIMSR)
        IF(TYPE.NE.99) LALBG=.FALSE.
      ENDIF
      CALL LCMSIX(IPLIB,'MACROLIB',1)
      CALL MCRAGF(IPLIB,IPMPO,IACCS,NMIL,NMIX,NGRP,NALBP,LALBG,LADFM,
     1 IMPX,NCAL,TERP,MIXC,NSURFD,HEDIT,VOSAP,VOLMI2,IDF)
      IF(NSURFD.GT.0) THEN
        CALL LCMGET(IPLIB,'STATE-VECTOR',ISTATE)
        ISTATE(12)=IDF ! ADF information
        CALL LCMPUT(IPLIB,'STATE-VECTOR',NSTATE,1,ISTATE)
      ENDIF
      CALL LCMSIX(IPLIB,' ',2)
      IF(NSURFD.GT.0) THEN
        CALL LCMGET(IPLIB,'STATE-VECTOR',ISTATE)
        ISTATE(24)=IDF ! ADF information
        CALL LCMPUT(IPLIB,'STATE-VECTOR',NSTATE,1,ISTATE)
      ENDIF
      IACCS=1
*----
*  SCRATCH STORAGE DEALLOCATION
*----
      DEALLOCATE(VOSAP)
      DEALLOCATE(SPH,FLUX,VOLMI2,VOL2,DENS3,DENS2)
      DEALLOCATE(HNAM2,HUSE2,ISTY2,ISTY1,ITOD2,IMIX2)
      RETURN
      END