diff options
Diffstat (limited to 'Donjon/src/MCRLIB.f')
| -rw-r--r-- | Donjon/src/MCRLIB.f | 856 |
1 files changed, 856 insertions, 0 deletions
diff --git a/Donjon/src/MCRLIB.f b/Donjon/src/MCRLIB.f new file mode 100644 index 0000000..ab94d3f --- /dev/null +++ b/Donjon/src/MCRLIB.f @@ -0,0 +1,856 @@ +*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 |