summaryrefslogtreecommitdiff
path: root/Dragon/src/SPHMPO.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/SPHMPO.f')
-rw-r--r--Dragon/src/SPHMPO.f716
1 files changed, 716 insertions, 0 deletions
diff --git a/Dragon/src/SPHMPO.f b/Dragon/src/SPHMPO.f
new file mode 100644
index 0000000..64acbf8
--- /dev/null
+++ b/Dragon/src/SPHMPO.f
@@ -0,0 +1,716 @@
+*DECK SPHMPO
+ SUBROUTINE SPHMPO(IPMPO,IPMAC,ICAL,IMPX,HEQUI,HMASL,NMIL,NALBP,
+ 1 NGROUP,HEDIT,VOSAP,ILUPS,SPH,B2)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Extract a Macrolib corresponding to an elementary calculation in a
+* MPO file
+*
+*Copyright:
+* Copyright (C) 2022 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
+* IPMPO pointer to the MPO file.
+* ICAL index of the elementary calculation being considered.
+* 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.
+* NMIL number of mixtures in the elementary calculation.
+* NALBP number of physical albedos per energy group.
+* NGROUP number of energy groups in the elementary calculation.
+* HEDIT name of output group for a (multigroup mesh, output geometry)
+* couple (generally equal to 'output_0').
+* VOSAP mixture volumes in the MPO file.
+* ILUPS up-scattering removing flag (=1 to remove up-scattering from
+* output cross-sections).
+* B2 imposed buckling.
+*
+*Parameters: output
+* IPMAC pointer to the Macrolib (L_MACROLIB signature).
+* SPH SPH-factor set extracted from the MPO file.
+* B2 buckling recovered from MPO file.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+ USE hdf5_wrap
+ IMPLICIT NONE
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ TYPE(C_PTR) IPMPO,IPMAC
+ INTEGER ICAL,IMPX,NMIL,NALBP,NGROUP,ILUPS
+ REAL VOSAP(NMIL),SPH(NMIL+NALBP,NGROUP),B2
+ CHARACTER(LEN=80) HEQUI,HMASL
+ CHARACTER(LEN=12) HEDIT
+*----
+* LOCAL VARIABLES
+*----
+ INTEGER, PARAMETER::NSTATE=40
+ INTEGER ISTATE(NSTATE)
+ INTEGER I,NADDRXS,NREA,NISO,NISOM,ADDRZX,ADDRZI,IBM,ISO,ISOM,
+ & NBISO,NBYTE,RANK,TYPE,DIMSR(5),NL,ILOC,NLOC,IPROF,JOFS,NL1,
+ & NL2,IPRC,NPRC,IOF,IGR,JGR,ITRANC,NED,IFISS,IGMAX,IGMIN,
+ & IPOSDE,IL,IREA,NSURFD,IDF,NALBP2,TYPE2,TYPE4,NITMA
+ REAL FLOTT,DEN,ZIL,FF,CSCAT
+ LOGICAL LSPH,LMASL,LSTRD,LDIFF,LHFACT,LNEW
+ CHARACTER RECNAM*80,RECNA2*80,TEXT12*12,CM*2,HSMG*131
+ TYPE(C_PTR) JPMAC,KPMAC
+*----
+* ALLOCATABLE ARRAYS
+*----
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: REACTION,ISOTOPE,IDATAP,
+ & ADDRISO,LOCAD,FAG,ADR
+ INTEGER, ALLOCATABLE, DIMENSION(:,:,:) :: ADDRXS
+ REAL, ALLOCATABLE, DIMENSION(:) :: FLUXS,CONCEN,RDATAX,RVALO,
+ & FMASL,LAMBDAD,BETADF,INVEL,LAMB,VREAL
+ REAL, ALLOCATABLE, DIMENSION(:,:) :: NWT0,EFACT,CHID,SIGS0,TOTAL,
+ & DIFF,BETAR,INVELS,DISFAC,ALBP,VFLUX
+ REAL, ALLOCATABLE, DIMENSION(:,:,:) :: XS,SIGS,CHIRS,ALBP_ERM,
+ & SFLUX
+ REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: SS2D
+ CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: HEDI,HADF
+ CHARACTER(LEN=24), ALLOCATABLE, DIMENSION(:) :: TEXT24,NOMREA,
+ & NOMISO
+ LOGICAL, ALLOCATABLE, DIMENSION(:) :: LXS
+ CHARACTER(LEN=80), ALLOCATABLE, DIMENSION(:) :: LOCTYP,LOCKEY
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: IPOS,NJJM,IJJM
+ REAL, ALLOCATABLE, DIMENSION(:) :: SCAT,GAR
+*----
+* SCRATCH STORAGE ALLOCATION
+* SIGS0 P0 scattering cross sections.
+* TOTAL total cross sections.
+* DIFF diffusion coefficients.
+* FMASL heavy element mass.
+*----
+ ALLOCATE(IPOS(NMIL),NJJM(NMIL),IJJM(NMIL))
+ ALLOCATE(SIGS0(NMIL,NGROUP),TOTAL(NMIL,NGROUP),DIFF(NMIL,NGROUP),
+ & FMASL(NMIL))
+ SIGS0(:NMIL,:NGROUP)=0.0
+ TOTAL(:NMIL,:NGROUP)=0.0
+ DIFF(:NMIL,:NGROUP)=0.0
+ FMASL(:NMIL)=0.0
+*----
+* RECOVER MPO FILE CHARACTERISTICS
+*----
+ 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)//"NREA",NREA)
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"NISO",NISO)
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ADDRXS",ADDRXS)
+ 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)//"TRANSPROFILE",IDATAP)
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ADDRISO",ADDRISO)
+ NBISO=ADDRISO(SIZE(ADDRISO,1))
+ IF(NBISO.EQ.0) CALL XABORT('SPHMPO: NO CROSS SECTIONS.')
+ ALLOCATE(NOMREA(NREA),NOMISO(NBISO))
+ CALL hdf5_read_data(IPMPO,"contents/reactions/REACTIONAME",TEXT24)
+ DO I=1,NREA
+ NOMREA(I)=TEXT24(REACTION(I)+1)
+ ENDDO
+ DEALLOCATE(TEXT24)
+ CALL hdf5_read_data(IPMPO,"contents/isotopes/ISOTOPENAME",TEXT24)
+ DO I=1,NBISO
+ NOMISO(I)=TEXT24(ISOTOPE(I)+1)
+ ENDDO
+ DEALLOCATE(TEXT24)
+ IF(IMPX.GT.1) THEN
+ WRITE(6,'(/24H SPHMPO: reaction names:)')
+ DO I=1,NREA
+ WRITE(6,'(5X,7HNOMREA(,I3,2H)=,A)') I,TRIM(NOMREA(I))
+ ENDDO
+ WRITE(6,'(/23H SPHMPO: isotope names:)')
+ DO I=1,NBISO
+ WRITE(6,'(5X,7HNOMISO(,I3,2H)=,A)') I,TRIM(NOMISO(I))
+ ENDDO
+ WRITE(6,'(/34H SPHMPO: number of energy groups =,I4)') NGROUP
+ WRITE(6,'(30H SPHMPO: number of mixtures =,I4)') NMIL
+ WRITE(6,'(30H SPHMPO: number of reactions =,I4)') NREA
+ WRITE(6,'(30H SPHMPO: number of isotopes =,I4)') NBISO
+ ENDIF
+ WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0)') TRIM(HEDIT),ICAL-1
+ IF(.not.hdf5_group_exists(IPMPO,TRIM(RECNAM))) THEN
+ WRITE(HSMG,'(38HSPHMPO: missing elementary calculation,I5,
+ & 10H in group ,A,1H.)') ICAL,TRIM(HEDIT)
+ CALL XABORT(HSMG)
+ 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.
+*----
+ NL=0
+ DO I=1,NADDRXS
+ DO ISO=1,NISO
+ 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 SPHMPO: number of legendre orders =,I4)') NL
+ ENDIF
+*----
+* RECOVER GENERAL INFORMATION
+*----
+ LSTRD=(B2.EQ.0.0)
+ WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,8H/addons/)')
+ & TRIM(HEDIT),ICAL-1
+ CALL hdf5_info(IPMPO,TRIM(RECNAM)//"KEFF",RANK,TYPE,NBYTE,DIMSR)
+ IF(TYPE.NE.99) THEN
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"KEFF",FLOTT)
+ CALL LCMPUT(IPMAC,'K-EFFECTIVE',1,2,FLOTT)
+ IF(IMPX.GT.1) THEN
+ WRITE(6,'(22H SPHMPO: K-EFFECTIVE =,1P,E13.6)') FLOTT
+ ENDIF
+ ENDIF
+ CALL hdf5_info(IPMPO,TRIM(RECNAM)//"KINF",RANK,TYPE,NBYTE,DIMSR)
+ IF(TYPE.NE.99) THEN
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"KINF",FLOTT)
+ CALL LCMPUT(IPMAC,'K-INFINITY',1,2,FLOTT)
+ IF(IMPX.GT.1) THEN
+ WRITE(6,'(21H SPHMPO: K-INFINITY =,1P,E13.6)') FLOTT
+ ENDIF
+ ENDIF
+ CALL hdf5_info(IPMPO,TRIM(RECNAM)//"B2",RANK,TYPE,NBYTE,DIMSR)
+ IF(TYPE.NE.99) THEN
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"B2",B2)
+ LSTRD=(B2.EQ.0.0)
+ CALL LCMPUT(IPMAC,'B2 B1HOM',1,2,B2)
+ IF(IMPX.GT.1) THEN
+ WRITE(6,'(13H SPHMPO: B2 =,1P,E14.6)') B2
+ ENDIF
+ ENDIF
+*----
+* SET NSURFD
+*----
+ IDF=0
+ NSURFD=0
+ WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,7H/zone_0,
+ & 15H/discontinuity/)') TRIM(HEDIT),ICAL-1
+ LNEW=hdf5_group_exists(IPMPO,TRIM(RECNAM))
+ IF(LNEW) THEN
+* new specification
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"NSURF",NSURFD)
+ ELSE
+* old specification
+ WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,12H/flux/NSURF/)')
+ & TRIM(HEDIT),ICAL-1
+ CALL hdf5_info(IPMPO,TRIM(RECNAM),RANK,TYPE,NBYTE,DIMSR)
+ IF(TYPE.NE.99) CALL hdf5_read_data(IPMPO,TRIM(RECNAM),NSURFD)
+ ENDIF
+ IF(NSURFD.EQ.0) GO TO 10
+*----
+* RECOVER DISCONTINUITY FACTOR INFORMATION
+*----
+ IF(LNEW) THEN
+* new specification
+ ALLOCATE(SFLUX(NMIL,NGROUP**2,NSURFD),VFLUX(NMIL,NGROUP))
+ DO IBM=1,NMIL
+ WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,6H/zone_,I0,1H/)')
+ & TRIM(HEDIT),ICAL-1,IBM-1
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ZONEFLUX",VREAL)
+ VFLUX(IBM,:NGROUP)=VREAL(:NGROUP)/VOSAP(IBM)
+ DEALLOCATE(VREAL)
+ WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,6H/zone_,I0,
+ & 15H/discontinuity/)') TRIM(HEDIT),ICAL-1,IBM-1
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"NSURF",NITMA)
+ IF(NITMA.NE.NSURFD) THEN
+ WRITE(HSMG,'(32HSPHMPO: THE NUMBER OF SURFACES (,I5,
+ & 12H) IN MIXTURE,I5,31H IS DIFFERENT FROM THE NUMBER (,I5,
+ & 15H) IN MIXTURE 1.)') NITMA,IBM,NSURFD
+ CALL XABORT(HSMG)
+ ENDIF
+ CALL hdf5_info(IPMPO,TRIM(RECNAM)//"DFACTOR",RANK,TYPE2,
+ & NBYTE,DIMSR)
+ CALL hdf5_info(IPMPO,TRIM(RECNAM)//"DFACTORGxG",RANK,TYPE4,
+ & NBYTE,DIMSR)
+ IF(TYPE2.NE.99) THEN
+ IDF=3 ! discontinuity factor information
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"DFACTOR",DISFAC)
+ DO I=1,NSURFD
+ SFLUX(IBM,:NGROUP,I)=DISFAC(I,:NGROUP)
+ ENDDO
+ DEALLOCATE(DISFAC)
+ ELSE IF(TYPE4.NE.99) THEN
+ IDF=4 ! matrix discontinuity factor information
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"DFACTORGxG",DISFAC)
+ DO I=1,NSURFD
+ SFLUX(IBM,:NGROUP**2,I)=DISFAC(I,:NGROUP**2)
+ ENDDO
+ DEALLOCATE(DISFAC)
+ ELSE
+ CALL hdf5_list(IPMPO,TRIM(RECNAM))
+ CALL XABORT('SPHMPO: UNABLE TO SET TYPE OF DF.')
+ ENDIF
+ ENDDO
+ ELSE
+* old specification
+ ALLOCATE(SFLUX(NMIL,NGROUP,NSURFD),VFLUX(NMIL,NGROUP))
+ IDF=3 ! discontinuity factor information
+ CALL SPHMOL(IPMPO,ICAL,NMIL,NGROUP,NSURFD,HEDIT,VOSAP,SFLUX,
+ 1 VFLUX)
+ ENDIF
+*----
+* WRITE DISCONTINUITY FACTOR INFORMATION ON IPMAC
+*----
+ CALL LCMSIX(IPMAC,'ADF',1)
+ CALL LCMPUT(IPMAC,'NTYPE',1,1,NSURFD)
+ CALL LCMPUT(IPMAC,'AVG_FLUX',NMIL*NGROUP,2,VFLUX)
+ ALLOCATE(HADF(NSURFD))
+ IF((IDF.EQ.2).OR.(IDF.EQ.3)) THEN
+ DO I=1,NSURFD
+ WRITE(HADF(I),'(3HFD_,I5.5)') I
+ CALL LCMPUT(IPMAC,HADF(I),NMIL*NGROUP,2,SFLUX(1,1,I))
+ ENDDO
+ ELSE IF(IDF.EQ.4) THEN
+ DO I=1,NSURFD
+ WRITE(HADF(I),'(3HFD_,I5.5)') I
+ CALL LCMPUT(IPMAC,HADF(I),NMIL*NGROUP**2,2,SFLUX(1,1,I))
+ ENDDO
+ ENDIF
+ CALL LCMPTC(IPMAC,'HADF',8,NSURFD,HADF)
+ DEALLOCATE(VFLUX,SFLUX)
+ CALL LCMSIX(IPMAC,' ',2)
+*----
+* RECOVER ALBEDO INFORMATION
+*----
+ 10 WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,6H/flux/)')
+ & TRIM(HEDIT),ICAL-1
+ CALL hdf5_info(IPMPO,TRIM(RECNAM)//"ALBEDO",RANK,TYPE,NBYTE,
+ & DIMSR)
+ IF(TYPE.NE.99) THEN
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"NALBP",NALBP2)
+ IF(NALBP2.NE.NALBP) CALL XABORT('SPHMPO: INVALID NALBP(1).')
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ALBEDO",ALBP)
+ CALL LCMPUT(IPMAC,'ALBEDO',NALBP*NGROUP,2,ALBP)
+ DEALLOCATE(ALBP)
+ ENDIF
+ CALL hdf5_info(IPMPO,TRIM(RECNAM)//"ALBEDOGxG",RANK,TYPE,NBYTE,
+ & DIMSR)
+ IF(TYPE.NE.99) THEN
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"NALBP",NALBP2)
+ IF(NALBP2.NE.NALBP) CALL XABORT('SPHMPO: INVALID NALBP(2).')
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ALBEDOGxG",ALBP_ERM)
+ CALL LCMPUT(IPMAC,'ALBEDO',NALBP*NGROUP*NGROUP,2,ALBP_ERM)
+ DEALLOCATE(ALBP_ERM)
+ ENDIF
+*----
+* ALLOCATE MACROLIB WORKING ARRAYS.
+*----
+ ALLOCATE(LXS(NREA),NWT0(NMIL,NGROUP),EFACT(NMIL,NGROUP),
+ 1 SIGS(NMIL,NGROUP,NL),SS2D(NMIL,NGROUP,NGROUP,NL),
+ 2 XS(NMIL,NGROUP,NREA),FAG(NGROUP),ADR(NGROUP))
+ NWT0(:NMIL,:NGROUP)=0.0
+ EFACT(:NMIL,:NGROUP)=0.0
+ SIGS(:NMIL,:NGROUP,:NL)=0.0
+ SS2D(:NMIL,:NGROUP,:NGROUP,:NL)=0.0
+ XS(:NMIL,:NGROUP,:NREA)=0.0
+ LXS(:NREA)=.FALSE.
+ LDIFF=.FALSE.
+*----
+* ALLOCATE DELAYED NEUTRON WORKING ARRAYS.
+*----
+ WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,
+ & 24H/zone_0/kinetics/LAMBDAD)') TRIM(HEDIT),ICAL-1
+ NPRC=0
+ CALL hdf5_info(IPMPO,TRIM(RECNAM),RANK,TYPE,NBYTE,DIMSR)
+ IF(TYPE.NE.99) NPRC=DIMSR(1)
+ IF(IMPX.GT.1) THEN
+ WRITE(6,'(37H SPHMPO: number of precursor groups =,I4)') NPRC
+ ENDIF
+ ALLOCATE(LAMB(NPRC),CHIRS(NGROUP,NPRC,NMIL),BETAR(NPRC,NMIL),
+ & INVELS(NGROUP,NMIL))
+ LAMB(:NPRC)=0.0
+ CHIRS(:NGROUP,:NPRC,:NMIL)=0.0
+ BETAR(:NPRC,:NMIL)=0.0
+ INVELS(:NGROUP,:NMIL)=0.0
+*----
+* LOOP OVER MPO MIXTURES
+*----
+ DO IBM=1,NMIL
+ WRITE(RECNAM,'(8H/output/,A,9H/statept_,I0,6H/zone_,I0,1H/)')
+ & TRIM(HEDIT),ICAL-1,IBM-1
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ZONEFLUX",FLUXS)
+ DO I=1,NGROUP
+ NWT0(IBM,I)=NWT0(IBM,I)+FLUXS(I)
+ ENDDO
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ADDRZI",ADDRZI)
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"ADDRZX",ADDRZX)
+ NISOM=ADDRISO(ADDRZI+2)-ADDRISO(ADDRZI+1)
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"CONCENTRATION",CONCEN)
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"CROSSECTION",RDATAX)
+ DO ISOM=1,NISOM
+ ! loop over the isotopes present in the mix
+ NL1=ADDRXS(NREA+1,ISOM,ADDRZX+1)
+ NL2=ADDRXS(NREA+2,ISOM,ADDRZX+1)
+ IF((NL1.GT.NL).OR.(NL2.GT.NL)) THEN
+ CALL XABORT('SPHMPO: NL OVERFLOW.')
+ ENDIF
+ DEN=CONCEN(ISOM)
+ IF(DEN.NE.0.0) THEN
+ IF(IMPX.GT.3) THEN
+ ISO=ADDRISO(ADDRZI+1)+ISOM
+ WRITE(6,'(10H mixture=,I5,15H concentration(,A,2H)=,1P,
+ & E12.4)') IBM,TRIM(NOMISO(ISO)),DEN
+ ENDIF
+ DO IREA=1,NREA
+ IOF=ADDRXS(IREA,ISOM,ADDRZX+1)
+ IF(IOF.LT.0) CYCLE
+ IF(NOMREA(IREA).EQ.'Diffusion') THEN
+ DO IL=1,NL1
+ DO IGR=1,NGROUP
+ FLOTT=DEN*RDATAX(IOF+(IL-1)*NGROUP+IGR)
+ SIGS(IBM,IGR,IL)=SIGS(IBM,IGR,IL)+FLOTT
+ LXS(IREA)=LXS(IREA).OR.(FLOTT.NE.0.0)
+ ENDDO
+ ENDDO
+ ELSE IF(NOMREA(IREA).EQ.'Scattering') THEN
+ IPROF=ADDRXS(NREA+3,ISOM,ADDRZX+1)
+ DO IGR=1,NGROUP
+ FAG(IGR)=IDATAP(IPROF+IGR)+1
+ ADR(IGR)=IDATAP(IPROF+NGROUP+IGR)
+ ENDDO
+ ADR(NGROUP+1)=IDATAP(IPROF+1+2*NGROUP)
+ JOFS=0
+ DO IL=1,NL2
+ ZIL=REAL(2*IL-1)
+ DO IGR=1,NGROUP
+ DO JGR=FAG(IGR),FAG(IGR)+(ADR(IGR+1)-ADR(IGR))-1
+ IF(JGR.GT.NGROUP) CALL XABORT('SPHMPO: SS2D OVER'
+ & //'FLOW.')
+ FLOTT=DEN*RDATAX(IOF+JOFS+1)/ZIL
+ SS2D(IBM,JGR,IGR,IL)=SS2D(IBM,JGR,IGR,IL)+FLOTT ! JGR <-- IGR
+ JOFS=JOFS+1
+ LXS(IREA)=LXS(IREA).OR.(FLOTT.NE.0.0)
+ ENDDO
+ ENDDO
+ ENDDO
+ ELSE IF(NOMREA(IREA).EQ.'FissionSpectrum') THEN
+ DO IGR=1,NGROUP
+ XS(IBM,IGR,IREA)=RDATAX(IOF+IGR)
+ LXS(IREA)=LXS(IREA).OR.(RDATAX(IOF+IGR).NE.0.0)
+ ENDDO
+ ELSE
+ DO IGR=1,NGROUP
+ XS(IBM,IGR,IREA)=XS(IBM,IGR,IREA)+DEN*RDATAX(IOF+IGR)
+ LXS(IREA)=LXS(IREA).OR.(DEN*RDATAX(IOF+IGR).NE.0.0)
+ ENDDO
+ ENDIF
+ ENDDO ! end of loop over reactions
+ ENDIF
+ ENDDO ! end of loop over isotopes
+ DEALLOCATE(RDATAX,CONCEN,FLUXS)
+*
+ IF(hdf5_group_exists(IPMPO,TRIM(RECNAM)//"kinetics")) THEN
+ WRITE(RECNA2,'(A,9Hkinetics/)') TRIM(RECNAM)
+ CALL hdf5_list(IPMPO,TRIM(RECNA2))
+ CALL hdf5_read_data(IPMPO,TRIM(RECNA2)//"LAMBDAD",LAMBDAD)
+ IF(SIZE(LAMBDAD,1).NE.NPRC) CALL XABORT('SPHMPO: WRONG NPRC.')
+ CALL hdf5_read_data(IPMPO,TRIM(RECNA2)//"CHID",CHID)
+ CALL hdf5_read_data(IPMPO,TRIM(RECNA2)//"BETADF",BETADF)
+ CALL hdf5_read_data(IPMPO,TRIM(RECNA2)//"INVERSESPEED",INVEL)
+ CHIRS(:NGROUP,:NPRC,IBM)=CHID(:NGROUP,:NPRC)
+ BETAR(:NPRC,IBM)=BETADF(:NPRC)
+ INVELS(:NGROUP,IBM)=INVEL(:NGROUP)
+ DEALLOCATE(INVEL,BETADF,CHID,LAMBDAD)
+ ENDIF
+*
+* UP-SCATTERING CORRECTION OF THE MACROLIB.
+ IF(ILUPS.EQ.1) THEN
+ DO JGR=2,NGROUP
+ DO IGR=1,JGR-1 ! IGR < JGR
+ FF=NWT0(IBM,JGR)/NWT0(IBM,IGR)
+ CSCAT=SS2D(IBM,IGR,JGR,1) ! IGR < JGR
+ DO IL=1,NL
+ CSCAT=SS2D(IBM,IGR,JGR,IL)
+ SIGS(IBM,IGR,IL)=SIGS(IBM,IGR,IL)-CSCAT*FF
+ SIGS(IBM,JGR,IL)=SIGS(IBM,JGR,IL)-CSCAT
+ SS2D(IBM,JGR,IGR,IL)=SS2D(IBM,JGR,IGR,IL)-CSCAT*FF
+ SS2D(IBM,IGR,JGR,IL)=0.0
+ ENDDO
+ ENDDO
+ ENDDO
+ ENDIF
+*
+ IF(LSPH) 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.NGROUP) THEN
+ CALL XABORT('SPHMPO: INVALID NUMBER OF COMPONENTS FOR '
+ & //'SPH FACTORS')
+ ENDIF
+ DO IGR=1,NGROUP
+ SPH(IBM,IGR)=RVALO(LOCAD(ILOC)+IGR-1)
+ ENDDO
+ ENDIF
+ ENDDO
+ DEALLOCATE(LOCAD,RVALO)
+ ELSE
+ SPH(IBM,:NGROUP)=1.0
+ ENDIF
+ IF(LMASL) 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.'HEAVY_METAL_DENSITY').AND.
+ & (LOCKEY(ILOC).EQ.HMASL)) THEN
+ IF(LOCAD(ILOC+1)-LOCAD(ILOC).NE.1) THEN
+ CALL XABORT('SPHMPO: INVALID NUMBER OF COMPONENTS FOR '
+ & //'HEAVY_METAL_DENSITY')
+ ENDIF
+ FMASL(IBM)=RVALO(LOCAD(ILOC))
+ ENDIF
+ ENDDO
+ DEALLOCATE(RVALO)
+ ENDIF
+*----
+* RECOVER DIFFUSION COEFFICIENT INFORMATION
+*----
+ IF(hdf5_group_exists(IPMPO,TRIM(RECNAM)//"leakage")) THEN
+ CALL hdf5_info(IPMPO,TRIM(RECNAM)//"leakage/DIFFCOEF",RANK,
+ & TYPE,NBYTE,DIMSR)
+ IF(TYPE.NE.99) THEN
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"leakage/DIFFCOEF",
+ & VREAL)
+ DO IGR=1,NGROUP
+ DIFF(IBM,IGR)=VREAL(IGR)
+ ENDDO
+ DEALLOCATE(VREAL)
+ LDIFF=.TRUE.
+ LSTRD=.FALSE.
+ GO TO 20
+ ENDIF
+ CALL hdf5_info(IPMPO,TRIM(RECNAM)//"leakage/DB2",RANK,TYPE,
+ & NBYTE,DIMSR)
+ IF(TYPE.NE.99) THEN
+ CALL hdf5_read_data(IPMPO,TRIM(RECNAM)//"leakage/DB2",VREAL)
+ DO IGR=1,NGROUP
+ DIFF(IBM,IGR)=VREAL(IGR)/B2
+ ENDDO
+ DEALLOCATE(VREAL)
+ LDIFF=.TRUE.
+ LSTRD=.FALSE.
+ ENDIF
+ ENDIF
+ 20 CONTINUE
+ ENDDO ! end of loop over mixtures
+ IF(NALBP.GT.0) THEN
+ SPH(NMIL+1:NMIL+NALBP,:NGROUP)=1.0 ! assigned to albedo function
+ ENDIF
+ DEALLOCATE(ADDRISO,IDATAP,ISOTOPE,REACTION,ADDRXS)
+*----
+* IDENTIFY SPECIAL FLUX EDITS
+*----
+ ALLOCATE(HEDI(NREA))
+ NED=0
+ DO 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
+ NED=NED+1
+ IF(NOMREA(IREA).EQ.'Fission') THEN
+ HEDI(NED)='NFTOT'
+ ELSE
+ HEDI(NED)=NOMREA(IREA)(:8)
+ ENDIF
+ ENDDO
+*----
+* STORE MACROLIB.
+*----
+ CALL LCMPUT(IPMAC,'VOLUME',NMIL,2,VOSAP)
+ IF(LMASL) CALL LCMPUT(IPMAC,'MASL',NMIL,2,FMASL)
+ IFISS=0
+ ITRANC=0
+ LHFACT=.FALSE.
+ ALLOCATE(VREAL(NMIL))
+ JPMAC=LCMLID(IPMAC,'GROUP',NGROUP)
+ DO IGR=1,NGROUP
+ KPMAC=LCMDIL(JPMAC,IGR)
+ CALL LCMPUT(KPMAC,'FLUX-INTG',NMIL,2,NWT0(1,IGR))
+ IF(NPRC.GT.0) THEN
+ DO IBM=1,NMIL
+ VREAL(IBM)=INVELS(IGR,IBM)
+ ENDDO
+ CALL LCMPUT(KPMAC,'OVERV',NMIL,2,VREAL)
+ ENDIF
+ DO IREA=1,NREA
+ IF(.NOT.LXS(IREA)) CYCLE
+ IF(NOMREA(IREA).EQ.'Absorption') THEN
+ TOTAL(:,IGR)=TOTAL(:,IGR)+XS(:,IGR,IREA)
+ ELSE IF(NOMREA(IREA).EQ.'Nexcess') THEN
+* correct scattering XS with excess XS
+ SIGS0(:,IGR)=SIGS0(:,IGR)+XS(:,IGR,IREA)
+ CALL LCMPUT(KPMAC,'N2N',NMIL,2,XS(1,IGR,IREA))
+ ELSE IF(NOMREA(IREA).EQ.'Fission') THEN
+ CALL LCMPUT(KPMAC,'NFTOT',NMIL,2,XS(1,IGR,IREA))
+ ELSE IF(NOMREA(IREA).EQ.'FissionSpectrum') THEN
+ CALL LCMPUT(KPMAC,'CHI',NMIL,2,XS(1,IGR,IREA))
+ DO IPRC=1,NPRC
+ DO IBM=1,NMIL
+ VREAL(IBM)=CHIRS(IGR,IPRC,IBM)
+ ENDDO
+ WRITE(TEXT12,'(A3,I2.2)') 'CHI',IPRC
+ CALL LCMPUT(KPMAC,TEXT12,NMIL,2,VREAL)
+ ENDDO
+ ELSE IF(NOMREA(IREA).EQ.'NuFission') THEN
+ IFISS=1
+ CALL LCMPUT(KPMAC,'NUSIGF',NMIL,2,XS(1,IGR,IREA))
+ DO IPRC=1,NPRC
+ DO IBM=1,NMIL
+ VREAL(IBM)=XS(IBM,IGR,IREA)*BETAR(IPRC,IBM)
+ ENDDO
+ WRITE(TEXT12,'(A6,I2.2)') 'NUSIGF',IPRC
+ CALL LCMPUT(KPMAC,TEXT12,NMIL,2,VREAL)
+ ENDDO
+ ELSE IF(NOMREA(IREA).EQ.'CaptureEnergyCapture') THEN
+ LHFACT=.TRUE.
+ EFACT(:,IGR)=EFACT(:,IGR)+XS(:,IGR,IREA)
+ ELSE IF(NOMREA(IREA).EQ.'FissionEnergyFission') THEN
+ LHFACT=.TRUE.
+ EFACT(:,IGR)=EFACT(:,IGR)+XS(:,IGR,IREA)
+ ELSE IF(NOMREA(IREA).EQ.'TransportCorrection') THEN
+ ITRANC=2
+ CALL LCMPUT(KPMAC,'TRANC',NMIL,2,XS(1,IGR,IREA))
+ ELSE IF(NOMREA(IREA).EQ.'Diffusion') THEN
+ DO IL=1,NL
+ WRITE(CM,'(I2.2)') IL-1
+ IF(IL.EQ.1) THEN
+ DO IBM=1,NMIL
+ SIGS0(IBM,IGR)=SIGS0(IBM,IGR)+SIGS(IBM,IGR,IL)
+ TOTAL(IBM,IGR)=TOTAL(IBM,IGR)+SIGS(IBM,IGR,IL)
+ ENDDO
+ ELSE
+ CALL LCMPUT(KPMAC,'SIGS'//CM,NMIL,2,SIGS(1,IGR,IL))
+ ENDIF
+ ENDDO
+ ELSE IF(NOMREA(IREA).EQ.'Scattering') THEN
+ ALLOCATE(SCAT(NGROUP*NMIL),GAR(NMIL))
+ DO IL=1,NL
+ WRITE(CM,'(I2.2)') IL-1
+ IPOSDE=0
+ DO IBM=1,NMIL
+ IPOS(IBM)=IPOSDE+1
+ IGMIN=IGR
+ IGMAX=IGR
+ DO JGR=NGROUP,1,-1
+ IF(SS2D(IBM,IGR,JGR,IL).NE.0.0) THEN
+ IGMIN=MIN(IGMIN,JGR)
+ IGMAX=MAX(IGMAX,JGR)
+ ENDIF
+ ENDDO
+ IJJM(IBM)=IGMAX
+ NJJM(IBM)=IGMAX-IGMIN+1
+ DO JGR=IGMAX,IGMIN,-1
+ IPOSDE=IPOSDE+1
+ SCAT(IPOSDE)=SS2D(IBM,IGR,JGR,IL) ! IGR <-- JGR
+ ENDDO
+ GAR(IBM)=SCAT(IPOS(IBM)+IJJM(IBM)-IGR)
+ ENDDO
+ CALL LCMPUT(KPMAC,'SCAT'//CM,IPOSDE,2,SCAT)
+ CALL LCMPUT(KPMAC,'NJJS'//CM,NMIL,1,NJJM)
+ CALL LCMPUT(KPMAC,'IJJS'//CM,NMIL,1,IJJM)
+ CALL LCMPUT(KPMAC,'IPOS'//CM,NMIL,1,IPOS)
+ CALL LCMPUT(KPMAC,'SIGW'//CM,NMIL,2,GAR)
+ ENDDO
+ DEALLOCATE(GAR,SCAT)
+ ELSE
+ CALL LCMPUT(KPMAC,NOMREA(IREA),NMIL,2,XS(1,IGR,IREA))
+ ENDIF
+ ENDDO ! end of loop over reactions
+ IF(LSTRD) THEN
+ IF((ITRANC.EQ.0).AND.(NL.GT.1)) THEN
+* Apollo-type transport correction
+ VREAL(:)=TOTAL(:,IGR)-SIGS(:,IGR,2)
+ ELSE
+ VREAL(:)=TOTAL(:,IGR)
+ ENDIF
+ DO IBM=1,NMIL
+ DIFF(IBM,IGR)=1.0/(3.0*VREAL(IBM))
+ ENDDO
+ LDIFF=.TRUE.
+ ENDIF
+ IF((ITRANC.EQ.0).AND.(NL.GT.1)) THEN
+* Apollo-type transport correction
+ IF(IGR.EQ.NGROUP) ITRANC=2
+ CALL LCMPUT(KPMAC,'TRANC',NMIL,2,SIGS(1,IGR,2))
+ ENDIF
+ IF(LDIFF) CALL LCMPUT(KPMAC,'DIFF',NMIL,2,DIFF(1,IGR))
+ IF(LHFACT) CALL LCMPUT(KPMAC,'H-FACTOR',NMIL,2,EFACT(1,IGR))
+ ENDDO
+ DEALLOCATE(VREAL)
+*----
+* RELEASE MEMORY
+*----
+ DEALLOCATE(INVELS,BETAR,CHIRS,LAMB)
+ DEALLOCATE(ADR,FAG,XS,SS2D,SIGS,EFACT,NWT0,LXS)
+ DEALLOCATE(NOMISO,NOMREA)
+*----
+* SAVE SCATTERING P0 AND TOTAL CROSS SECTION INFO
+*----
+ DO IGR=1,NGROUP
+ KPMAC=LCMDIL(JPMAC,IGR)
+ CALL LCMPUT(KPMAC,'SIGS00',NMIL,2,SIGS0(1,IGR))
+ CALL LCMPUT(KPMAC,'NTOT0',NMIL,2,TOTAL(1,IGR))
+ ENDDO
+*----
+* WRITE STATE VECTOR
+*----
+ IF(IMPX.GT.1) THEN
+ WRITE(6,'(32H SPHMPO: fissile isotope index =,I4)') IFISS
+ WRITE(6,'(37H SPHMPO: transport correction index =,I4)') ITRANC
+ ENDIF
+ TEXT12='L_MACROLIB'
+ CALL LCMPTC(IPMAC,'SIGNATURE',12,TEXT12)
+ ISTATE(:NSTATE)=0
+ ISTATE(1)=NGROUP
+ ISTATE(2)=NMIL
+ ISTATE(3)=NL ! 1+scattering anisotropy
+ ISTATE(4)=IFISS
+ ISTATE(5)=NED
+ ISTATE(6)=ITRANC
+ ISTATE(8)=NALBP
+ IF(LDIFF) ISTATE(9)=1
+ ISTATE(12)=IDF ! ADF information
+ CALL LCMPUT(IPMAC,'STATE-VECTOR',NSTATE,1,ISTATE)
+ IF(NED.GT.0) CALL LCMPTC(IPMAC,'ADDXSNAME-P0',8,NED,HEDI)
+ DEALLOCATE(HEDI)
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(FMASL,DIFF,TOTAL,SIGS0)
+ DEALLOCATE(IJJM,NJJM,IPOS)
+ RETURN
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 11:33:42 +0000