summaryrefslogtreecommitdiff
path: root/Dragon/src/AUTFLU.f
diff options
context:
space:
mode:
Diffstat (limited to 'Dragon/src/AUTFLU.f')
-rw-r--r--Dragon/src/AUTFLU.f379
1 files changed, 379 insertions, 0 deletions
diff --git a/Dragon/src/AUTFLU.f b/Dragon/src/AUTFLU.f
new file mode 100644
index 0000000..47fdeed
--- /dev/null
+++ b/Dragon/src/AUTFLU.f
@@ -0,0 +1,379 @@
+*DECK AUTFLU
+ SUBROUTINE AUTFLU(IPTRK,IPLIB,IPLI0,IFTRAK,NREG,NUN,NBMIX,NBISO,
+ 1 NIRES,MAT,VOL,KEYFLX,CDOOR,LEAKSW,IMPX,DEN,MIX,IAPT,IPHASE,NGRP,
+ 2 IGRMIN,IGRRES,IGRMAX,DIL,TITR,IALTER,DELI,LBIN,NBIN,EBIN,MAXTRA,
+ 3 ISEED,ITRANC,UUU,FUNKNO,SIGT,SIGS,SIGS1,SIGF,SIGGAR,MASKG)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the flux in the Autolib fine groups using the Autosecol
+* method.
+*
+*Copyright:
+* Copyright (C) 2023 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
+* IPTRK pointer to the tracking (L_TRACK signature).
+* IPLIB pointer to the internal microscopic cross section library
+* with subgroups (L_LIBRARY signature).
+* IPLI0 pointer to the internal microscopic cross section library
+* builded by the self-shielding module.
+* IFTRAK file unit number used to store the tracks.
+* NREG number of regions.
+* NUN number of unknowns per energy group and band.
+* NBMIX number of mixtures in the internal library.
+* NBISO number of isotopes.
+* NIRES number of correlated resonant isotopes.
+* MAT index-number of the mixture type assigned to each volume.
+* VOL volumes.
+* KEYFLX pointers of fluxes in unknown vector.
+* CDOOR name of the geometry/solution operator.
+* LEAKSW leakage flag (LEAKSW=.true. if neutron leakage through
+* external boundary is present).
+* IMPX print flag (equal to zero for no print).
+* DEN density of each isotope.
+* MIX mix number of each isotope (can be zero).
+* IAPT resonant isotope index associated with isotope I. Mixed
+* moderator if IAPT(I)=NIRES+1. Out-of-fuel isotope if
+* IAPT(I)=0.
+* IPHASE type of flux solution (=1 use a native flux solution door;
+* =2 use collision probabilities).
+* NGRP number of energy groups.
+* IGRMIN first group where the self-shielding is applied.
+* IGRRES first resolved energy group.
+* IGRMAX most thermal group where the self-shielding is applied.
+* DIL microscopic dilution cross section of each isotope.
+* TITR title.
+* IALTER type of elastic slowing-down kernel (=0: use exact kernel;
+* =1: use an approximate kernel for the resonant isotopes).
+* DELI elementary lethargy width used by the elastic kernel.
+* LBIN total number of fine energy groups in the Autolib.
+* NBIN number of fine energy groups in each coarse energy group.
+* EBIN energy limits of the Autolib fine groups.
+* MAXTRA maximum number of down-scattering terms.
+* ISEED the seed for the generation of random numbers in the
+* unresolved energy domain.
+* ITRANC type of transport correction.
+*
+*Parameters: output
+* UUU lethargy limits of the Autolib fine groups.
+* FUNKNO flux in the Autolib fine groups.
+* SIGT total microscopic x-s.
+* SIGS P0 scattering microscopic x-s.
+* SIGS1 P1 scattering microscopic x-s.
+* SIGF nu*fission microscopic x-s.
+* SIGGAR macroscopic x-s of the non-resonant isotopes in each mixture:
+* (*,*,*,1) total; (*,*,*,2) transport correction;
+* (*,*,*,3) P0 scattering.
+* MASKG energy group mask pointing on self-shielded groups.
+*
+*-----------------------------------------------------------------------
+*
+ USE GANLIB
+*----
+* SUBROUTINE ARGUMENTS
+*----
+ TYPE(C_PTR) IPTRK,IPLIB,IPLI0
+ INTEGER IFTRAK,NREG,NUN,NBMIX,NBISO,NIRES,MAT(NREG),KEYFLX(NREG),
+ 1 IMPX,MIX(NBISO),IAPT(NBISO),IPHASE,NGRP,IGRMIN,IGRRES,IGRMAX,
+ 2 IALTER,LBIN,NBIN(NGRP),MAXTRA,ISEED,ITRANC
+ REAL VOL(NREG),DEN(NBISO),DIL(NBISO),DELI,EBIN(LBIN+1),
+ 1 UUU(LBIN+1),FUNKNO(NUN,LBIN),SIGT(LBIN,NBISO),SIGS(LBIN,NBISO),
+ 2 SIGS1(LBIN,NBISO),SIGF(LBIN,NBISO),SIGGAR(NBMIX,0:NIRES,NGRP,3)
+ LOGICAL LEAKSW,MASKG(NGRP)
+ CHARACTER CDOOR*12,TITR*72
+*----
+* LOCAL VARIABLES
+*----
+ TYPE(C_PTR) IPP,KPLIB,KPLI0,IPSYS
+ LOGICAL LLIB
+ DOUBLE PRECISION DUUU
+ CHARACTER HNAMIS*12,HNABIS*12,HSMG*131
+*----
+* ALLOCATABLE ARRAYS
+*----
+ TYPE(C_PTR), ALLOCATABLE, DIMENSION(:) :: IPISO1,IPISO2
+ INTEGER, ALLOCATABLE, DIMENSION(:) :: NJJ,IJJ,NEXT,III
+ REAL, ALLOCATABLE, DIMENSION(:) :: GA1,DELTAU,SGAR,PRI
+ REAL, ALLOCATABLE, DIMENSION(:,:) :: SIGINF,SCAT,GA2,CONC
+ LOGICAL, ALLOCATABLE, DIMENSION(:) :: MASKH
+*----
+* SCRATCH STORAGE ALLOCATION
+*----
+ ALLOCATE(GA1(NGRP),GA2(NGRP,NGRP),CONC(NBMIX,NBISO),DELTAU(NGRP),
+ 1 MASKH(NGRP),SIGINF(NGRP,3))
+ ALLOCATE(IPISO1(NBISO),IPISO2(NBISO))
+ ALLOCATE(NJJ(NGRP),IJJ(NGRP),SGAR(NGRP**2),SCAT(NGRP,NGRP))
+ ALLOCATE(PRI(MAXTRA),NEXT(NBISO),III(NBISO+1))
+ NEXT(:NBISO)=0
+*
+ CALL KDRCPU(TK1)
+ SIGGAR(:NBMIX,0:NIRES,:NGRP,:3)=0.0
+ CALL LIBIPS(IPLIB,NBISO,IPISO1)
+ CALL LIBIPS(IPLI0,NBISO,IPISO2)
+ CALL LCMGET(IPLI0,'DELTAU',DELTAU)
+ III(1)=1
+ DO 120 ISO=1,NBISO
+ SIGT(:LBIN,ISO)=0.0
+ SIGS(:LBIN,ISO)=0.0
+ SIGS1(:LBIN,ISO)=0.0
+ SIGF(:LBIN,ISO)=0.0
+ IF(IMPX.GT.1) WRITE(6,'(/32H AUTFLU: RECOVER XS FOR ISOTOPE=,I5)')
+ 1 ISO
+ IBM=MIX(ISO)
+ DO 10 NRE=1,NREG
+ IF(MAT(NRE).EQ.IBM) GO TO 20
+ 10 CONTINUE
+ III(ISO+1)=III(ISO)
+ GO TO 120
+ 20 KPLIB=IPISO1(ISO) ! infinite dilution isotope
+ KPLI0=IPISO2(ISO) ! self-shielded isotope
+ CALL LCMGTC(KPLIB,'ALIAS',12,HNAMIS)
+ IRES=IAPT(ISO)
+ JRES=IRES
+ IF(IRES.EQ.NIRES+1) JRES=0
+ DENN=DEN(ISO)
+ CALL LCMGET(KPLIB,'AWR',AWR)
+ CALL LCMLEN(KPLIB,'BIN-NFS',LENGT,ITYLCM)
+ IF((IRES.GT.0).AND.(IRES.LE.NIRES).AND.(LENGT.GT.0)) THEN
+ ! resonant isotope
+ IF(IMPX.GT.2) WRITE(6,'(26H AUTFLU: PROCESS AUTOLIB '',A12,
+ 1 1H'')') HNAMIS
+ IF(LBIN.EQ.0) CALL XABORT('AUTFLU: MISSING AUTOLIB DATA.')
+ DUUU=0.0D0
+ UUU(1)=0.0
+ DO 40 IGR=1,LBIN
+ DUUU=DUUU+LOG(EBIN(IGR)/EBIN(IGR+1))
+ UUU(IGR+1)=REAL(DUUU)
+ 40 CONTINUE
+ ! recover unresolved xs values
+ LLL=0
+ IF(IGRRES.GT.IGRMIN) THEN
+ CALL LCMLEN(KPLIB,'NUSIGF',N10,ITYLCM)
+ CALL LCMLEN(KPLIB,'SIGS00',N12,ITYLCM)
+ SIGINF(:NGRP,:3)=0.0
+ CALL LCMGET(KPLIB,'NTOT0',SIGINF(1,1))
+ IF(N10.GT.0) CALL LCMGET(KPLIB,'NUSIGF',SIGINF(1,2))
+ IF(N12.GT.0) CALL LCMGET(KPLIB,'SIGS00',SIGINF(1,3))
+ CALL AUTTAB(KPLIB,HNAMIS,IGRMIN,IGRRES,NGRP,LBIN,NBIN,UUU,
+ 1 ISEED,SIGINF,LLL,SIGT(1,ISO),SIGS(1,ISO),SIGF(1,ISO))
+ ENDIF
+ ! recover resolved xs values
+ IF(LLL.LT.LBIN) THEN
+ CALL LCMLEN(KPLIB,'BIN-NTOT0',LENG,ITYLCM)
+ IF(LENG.GT.LBIN) CALL XABORT('AUTFLU: LBIN OVERFLOW.')
+ CALL LCMGET(KPLIB,'BIN-NTOT0',SIGT(LLL+1,ISO))
+ CALL LCMGET(KPLIB,'BIN-SIGS00',SIGS(LLL+1,ISO))
+ CALL LCMLEN(KPLIB,'BIN-NUSIGF',ILEN,ITYLCM)
+ IF(ILEN.GT.0) CALL LCMGET(KPLIB,'BIN-NUSIGF',SIGF(LLL+1,ISO))
+ ENDIF
+*----
+* ELASTIC SCATTERING INFORMATION.
+*----
+ MAXIII=MAXTRA-III(ISO)+1
+ CALL LIBPRI(MAXIII,DELI,AWR,IALTER,0,NEXT(ISO),PRI(III(ISO)))
+ ELSE
+ ! infinite dilution isotope
+ IF(C_ASSOCIATED(KPLI0)) THEN
+ CALL LCMLEN(KPLI0,'NTOT0',ILEN0,ITYLCM)
+ IF(ILEN0.NE.0) THEN
+ LLIB=.FALSE.
+ IPP=KPLI0
+ ELSE
+ LLIB=.TRUE.
+ IPP=IPISO1(ISO) ! set ISO-th isotope
+ ENDIF
+ ELSE
+ LLIB=.TRUE.
+ IPP=IPISO1(ISO) ! set ISO-th isotope
+ ENDIF
+ IF(LLIB.AND.(.NOT.C_ASSOCIATED(IPP))) THEN
+ WRITE(HSMG,'(18H AUTFLU: ISOTOPE '',A12,7H'' (ISO=,I8,5H) IS ,
+ 1 39HNOT AVAILABLE IN THE ORIGINAL MICROLIB.)') HNAMIS,ISO
+ CALL XABORT(HSMG)
+ ENDIF
+ IF((.NOT.LLIB).AND.(IMPX.GT.2)) THEN
+ CALL LCMGTC(KPLI0,'ALIAS',12,HNABIS)
+ WRITE(6,'(26H AUTFLU: RECOVER ISOTOPE '',A12,11H'' FROM THE ,
+ 1 12HNEW LIBRARY.)') HNABIS
+ ELSE IF(LLIB.AND.(IMPX.GT.2)) THEN
+ WRITE(6,'(26H AUTFLU: RECOVER ISOTOPE '',A12,11H'' FROM THE ,
+ 1 12HOLD LIBRARY.)') HNAMIS
+ ENDIF
+ NSCAR=0
+ SCAT(:NGRP,:NGRP)=0.0
+ CALL LCMLEN(IPP,'NTOT0',N13,ITYLCM)
+ IF(N13.EQ.0) THEN
+ CALL LCMLIB(IPP)
+ CALL XABORT('AUTFLU: NO INFINITE-DILUTION TOTAL XS.')
+ ELSE IF(N13.GT.LBIN) THEN
+ CALL XABORT('AUTFLU: LBIN OVERFLOW.')
+ ELSE IF(N13.NE.NGRP) THEN
+ CALL XABORT('AUTFLU: INVALID X-SECTIONS.')
+ ENDIF
+ CALL LCMGET(IPP,'NTOT0',SIGT(1,ISO))
+ CALL LCMLEN(IPP,'NUSIGF',N10,ITYLCM)
+ IF(N10.GT.0) CALL LCMGET(IPP,'NUSIGF',SIGF(1,ISO))
+ CALL LCMLEN(IPP,'SIGS00',N12,ITYLCM)
+ IF(N12.GT.0) THEN
+ CALL LCMGET(IPP,'SIGS00',SIGS(1,ISO))
+ CALL LCMLEN(IPP,'SIGS01',N14,ITYLCM)
+ IF(N14.GT.0) THEN
+ CALL LCMGET(IPP,'SIGS01',SIGS1(1,ISO))
+ ELSE
+ SIGS1(:NGRP,ISO)=0.0
+ ENDIF
+ ELSE
+ CALL LCMGET(IPP,'SCAT00',SGAR)
+ CALL LCMGET(IPP,'NJJS00',NJJ)
+ CALL LCMGET(IPP,'IJJS00',IJJ)
+ IGAR1=0
+ DO IG2=1,NGRP
+ DO IG1=IJJ(IG2),IJJ(IG2)-NJJ(IG2)+1,-1
+ IGAR1=IGAR1+1
+ SCAT(IG1,IG2)=SGAR(IGAR1)
+ ENDDO
+ ENDDO
+ DO IG1=1,NGRP
+ SUMSC=0.0D0
+ DO IG2=1,NGRP
+ SUMSC=SUMSC+SCAT(IG1,IG2)
+ ENDDO
+ SIGS(IG1,ISO)=REAL(SUMSC)
+ ENDDO
+ CALL LCMLEN(IPP,'SCAT01',N87,ITYLCM)
+ IF(N87.GT.0) THEN
+ CALL LCMGET(IPP,'SCAT01',SGAR)
+ CALL LCMGET(IPP,'NJJS01',NJJ)
+ CALL LCMGET(IPP,'IJJS01',IJJ)
+ IGAR1=0
+ DO IG2=1,NGRP
+ DO IG1=IJJ(IG2),IJJ(IG2)-NJJ(IG2)+1,-1
+ IGAR1=IGAR1+1
+ SCAT(IG1,IG2)=SGAR(IGAR1)
+ ENDDO
+ ENDDO
+ DO IG1=1,NGRP
+ SUMSC=0.0D0
+ DO IG2=1,NGRP
+ SUMSC=SUMSC+SCAT(IG1,IG2)
+ ENDDO
+ SIGS1(IG1,ISO)=REAL(SUMSC)
+ ENDDO
+ ENDIF
+ ENDIF
+ ! compute SIGGAR used by SPH equivalence
+ IF((DENN.NE.0.0).AND.(IBM.NE.0).AND.(JRES.EQ.0)) THEN
+ DO 70 IG1=1,NGRP
+ SIGGAR(IBM,JRES,IG1,1)=SIGGAR(IBM,JRES,IG1,1)+DENN*
+ 1 SIGT(IG1,ISO)
+ SIGGAR(IBM,JRES,IG1,3)=SIGGAR(IBM,JRES,IG1,3)+DENN*
+ 1 SIGS(IG1,ISO)
+ 70 CONTINUE
+ CALL LCMLEN(IPP,'TRANC',LENGT,ITYLCM)
+ IF(LENGT.GT.0) THEN
+ CALL LCMGET(IPP,'TRANC',GA1)
+ DO 80 IG1=1,NGRP
+ SIGGAR(IBM,JRES,IG1,2)=SIGGAR(IBM,JRES,IG1,2)+DENN*GA1(IG1)
+ 80 CONTINUE
+ ENDIF
+ ENDIF
+ ! expand xs of non-resonant isotopes
+ CALL AUTPRD(NGRP,LBIN,NBIN,SIGS(1,ISO))
+ CALL AUTPRD(NGRP,LBIN,NBIN,SIGT(1,ISO))
+ CALL AUTPRD(NGRP,LBIN,NBIN,SIGF(1,ISO))
+ CALL AUTPRD(NGRP,LBIN,NBIN,SIGS1(1,ISO))
+ ENDIF
+ III(ISO+1)=III(ISO)+NEXT(ISO)
+*----
+* CROSS SECTION EDITION.
+*----
+ IF(IMPX.GT.7) THEN
+ CALL LCMGTC(KPLIB,'ALIAS',12,HNAMIS)
+ WRITE(6,540) HNAMIS,(K,SIGS(K,ISO),SIGT(K,ISO),SIGF(K,ISO),
+ 1 SIGS1(K,ISO),K=1,LBIN)
+ I2=III(ISO+1)-1
+ WRITE(6,550) III(ISO),I2,(PRI(K),K=III(ISO),I2)
+ ENDIF
+ 120 CONTINUE
+ CALL KDRCPU(TK2)
+ IF(IMPX.GT.1) WRITE(6,'(/36H AUTFLU: CPU TIME SPENT TO RECOVER A,
+ 1 33HUTOLIB AND INFINITE-DILUTION XS =,F8.1,8H SECOND./)') TK2-TK1
+*
+ CALL KDRCPU(TK1)
+ TK4=0.0
+ TK5=0.0
+ ICPIJ=0
+*----
+* SET THE NUMBER DENSITIES.
+*----
+ CONC(:NBMIX,:NBISO)=0.0
+ DO 130 ISO=1,NBISO
+ IBM=MIX(ISO)
+ IF(IBM.LE.0) CYCLE
+ CONC(IBM,ISO)=DEN(ISO)
+ 130 CONTINUE
+*----
+* DETERMINE WHICH GROUPS ARE SELF-SHIELDED.
+*----
+ MASKG(:NGRP)=.FALSE.
+ MASKG(IGRMIN:IGRMAX)=.TRUE.
+*----
+* COMPUTE THE NEUTRON FLUX.
+*----
+ CALL KDRCPU(TKA)
+ CALL KDRCPU(TKA)
+ DO 330 IG1=1,NGRP
+ ICPIJ=ICPIJ+NBIN(IG1)
+ 330 CONTINUE
+ CALL LCMOP(IPSYS,'**SYSTEM**',0,1,0)
+ KNORM=1
+ IMPY=MAX(0,IMPX-3)
+ IF(IPHASE.EQ.1) THEN
+* USE A NATIVE DOOR.
+ CALL AUTIT2(IPTRK,IFTRAK,IPSYS,MAXTRA,KNORM,NUN,LBIN,NREG,
+ 1 NBMIX,NBISO,MAT,VOL,KEYFLX,NIRES,IAPT,CDOOR,LEAKSW,TITR,IMPY,
+ 2 CONC,SIGS,SIGT,SIGS1,DIL,PRI,UUU,DELI,ITRANC,NEXT,III,FUNKNO)
+ ELSE IF(IPHASE.EQ.2) THEN
+* USE A COLLISION PROBABILITY DOOR.
+ CALL AUTIT1(IPTRK,IFTRAK,IPSYS,MAXTRA,KNORM,LBIN,NREG,
+ 1 NBMIX,NBISO,MAT,VOL,NIRES,IAPT,CDOOR,LEAKSW,TITR,IMPY,CONC,
+ 2 SIGS,SIGT,SIGS1,DIL,PRI,UUU,DELI,ITRANC,NEXT,III,FUNKNO)
+ ENDIF
+ CALL LCMCL(IPSYS,2)
+ CALL KDRCPU(TKB)
+ TK4=TK4+(TKB-TKA)
+* ***************************************************************
+ CALL LCMVAL(IPLI0,' ')
+*----
+* RESET MASKG FOR SPH CALCULATION IN SMALL LETHARGY WIDTH GROUPS.
+*----
+ DO 380 IG1=1,NGRP
+ IF(MASKG(IG1)) THEN
+ IF(DELTAU(IG1).GT.0.1) MASKG(IG1)=.TRUE.
+ ENDIF
+ 380 CONTINUE
+ CALL KDRCPU(TK2)
+ IF(IMPX.GT.1) WRITE(6,'(/34H AUTFLU: CPU TIME SPENT TO COMPUTE,
+ 1 31H SELF-SHIELDED REACTION RATES =,F8.1,19H SECOND, INCLUDING:
+ 2 /9X,F8.1,36H SECOND TO SOLVE FOR AUTOSECOL FLUX;/9X,7HNUMBER ,
+ 3 25HOF ASSEMBLY DOORS CALLS =,I5,1H.)') TK2-TK1,TK4,ICPIJ
+*----
+* SCRATCH STORAGE DEALLOCATION
+*----
+ DEALLOCATE(III,NEXT,PRI)
+ DEALLOCATE(SCAT,SGAR,IJJ,NJJ)
+ DEALLOCATE(IPISO2,IPISO1)
+ DEALLOCATE(SIGINF,MASKH,DELTAU,CONC,GA2,GA1)
+ RETURN
+ 540 FORMAT(/18H AUTFLU: ISOTOPE ',A12,1H'/12X,4HSIGS,16X,4HSIGT,
+ 1 16X,4HSIGF,16X,5HSIGS1/(I5,1P,4E20.7))
+ 550 FORMAT(/29H AUTFLU: SCATTERING ELEMENTS:,2I10/(1P,10E13.5))
+ END
generated by cgit v1.2.3 (git 2.46.0) at 2026-02-16 17:31:50 +0000