Initial commit from Polytechnique Montreal

1 files changed, 351 insertions, 0 deletions

diff --git a/Dragon/src/TONSPH.f b/Dragon/src/TONSPH.f
new file mode 100644
index 0000000..812ac62
--- /dev/null
+++ b/Dragon/src/TONSPH.f
@@ -0,0 +1,351 @@
+*DECK TONSPH
+      SUBROUTINE TONSPH(IPLIB,IPTRK,IFTRAK,NREG,NUN,NBM,NBISO,ISONAM,
+     1 MAT,VOL,KEYFLX,CDOOR,INRS,LEAKSW,IMPX,DEN,MIX,LSHI,ITRANC,
+     2 IPHASE,NGRO,IGRMIN,IGRMAX,NBNRS,TITR,SIGT2,SIGT3,SN,SPH,ICPIJ,
+     3 TK3,TK4)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* SPH equivalence procedure over the self-shielded cross sections. Use
+* all the standard solution doors of Dragon.
+*
+*Copyright:
+* Copyright (C) 2017 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
+* IPLIB   pointer to the internal microscopic cross section library
+*         (L_LIBRARY signature).
+* IPTRK   pointer to the tracking (L_TRACK signature).
+* IFTRAK  file unit number used to store the tracks.
+* NREG    number of regions.
+* NUN     number of unknowns per energy group.
+* NBM     number of mixtures in the internal library.
+* NBISO   number of isotopes.
+* ISONAM  alias name of isotopes in IPLIB.
+* 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.
+* INRS    index of the resonant isotope under consideration.
+* 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).
+* LSHI    resonant region number associated with each isotope.
+*         Infinite dilution will be assumed if LSHI(i)=0.
+* ITRANC  type of transport correction.
+* IPHASE  type of flux solution (=1 use a native flux solution door;
+*         =2 use collision probabilities).
+* NGRO    number of energy groups.
+* IGRMIN  first group where the self-shielding is applied.
+* IGRMAX  most thermal group where the self-shielding is applied.
+* NBNRS   number of totally correlated fuel regions. NBNRS=max(IRES).
+* TITR    title.
+* SIGT2   total macroscopic cross sections.
+* SIGT3   transport correction.
+* SN      computed dilution cross section in each energy group of
+*         each isotope.
+*
+*Parameters: output
+* SPH     SPH factors.
+* ICPIJ   number of flux solution door calls.
+*
+*Parameters: input/output
+* TK3     cpu time to compute system matrices.
+* TK4     cpu time to compute fluxes.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+      USE DOORS_MOD
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPLIB,IPTRK
+      INTEGER IFTRAK,NREG,NUN,NBM,NBISO,ISONAM(3,NBISO),MAT(NREG),
+     1 KEYFLX(NREG),INRS,IMPX,MIX(NBISO),LSHI(NBISO),ITRANC,IPHASE,
+     2 NGRO,IGRMIN,IGRMAX,NBNRS,ICPIJ
+      REAL VOL(NREG),DEN(NBISO),SIGT2(0:NBM,NGRO),SIGT3(0:NBM,NGRO),
+     1 SN(NGRO,NBISO),SPH(NBM,NGRO),TK3,TK4
+      LOGICAL LEAKSW
+      CHARACTER CDOOR*12,TITR*72,HNAMIS*12
+*----
+*  LOCAL VARIABLES
+*----
+      PARAMETER(NSTATE=40)
+      TYPE(C_PTR) JPLIB,KPLIB,IPMACR,IPSOU
+      LOGICAL LHOMOG,LPROB,LEXAC,LOGDO,REBFLG
+      CHARACTER TEXT12*12
+      INTEGER NALBP,ISTATE(NSTATE)
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: IRES,ISONR,NPSYS
+      REAL, ALLOCATABLE, DIMENSION(:) :: VST,SIGTXS,SIGS0X,SIGG,FLNEW
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: FLUX,SIG0,SIG1,SIG3,TOTAL,
+     1 SIGS0,TRANC,PHGAR,SUNKNO,FUNKNO
+      LOGICAL, ALLOCATABLE, DIMENSION(:) :: MASKI
+      TYPE(C_PTR), ALLOCATABLE, DIMENSION(:) :: IPISO
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(IRES(NBM),ISONR(NBISO),NPSYS(NGRO))
+      ALLOCATE(VST(NBNRS),SIGTXS(0:NBM),SIGS0X(0:NBM),SIGG(0:NBM),
+     1 FLNEW(NBNRS),SUNKNO(NUN,NGRO),FUNKNO(NUN,NGRO),FLUX(NBM,NGRO),
+     2 SIG0(NBM,NGRO),SIG1(NBM,NGRO),SIG3(NBM,NGRO),TOTAL(NGRO,NBNRS),
+     3 SIGS0(NGRO,NBNRS),TRANC(NGRO,NBNRS),PHGAR(NGRO,NBNRS))
+      ALLOCATE(MASKI(NBISO))
+      ALLOCATE(IPISO(NBISO))
+*----
+*  FIND THE RESONANT MIXTURE NUMBERS AND THE CORRELATED ISOTOPES
+*  ASSOCIATED WITH REGION INRS
+*----
+      IRES(:NBM)=0
+      ISONR(:NBISO)=0
+      IRS=0
+      TEXT12=' '
+      DO 30 IBM=1,NBM
+      LOGDO=.FALSE.
+      DO 10 I=1,NREG
+      LOGDO=LOGDO.OR.(MAT(I).EQ.IBM)
+   10 CONTINUE
+      IF(.NOT.LOGDO) GO TO 30
+      DO 20 ISO=1,NBISO
+      IF((MIX(ISO).EQ.IBM).AND.(LSHI(ISO).EQ.INRS)) THEN
+         WRITE(HNAMIS,'(3A4)') (ISONAM(I0,ISO),I0=1,3)
+         IF(HNAMIS.NE.TEXT12) THEN
+           IRS=IRS+1
+           TEXT12=HNAMIS
+         ENDIF
+         ISONR(ISO)=IRS
+         IRES(IBM)=IRS
+      ENDIF
+   20 CONTINUE
+   30 CONTINUE
+      IF(IRS.NE.NBNRS) CALL XABORT('TONSPH: INVALID VALUE OF NBNRS.')
+*----
+*  SET THE LCM MICROLIB ISOTOPEWISE DIRECTORIES.
+*----
+      CALL LIBIPS(IPLIB,NBISO,IPISO)
+*----
+*  UNLOAD MICROSCOPIC X-S FROM LCM TO SCRATCH STORAGE.
+*----
+      DO 40 ISO=1,NBISO
+      IRS=ISONR(ISO)
+      IF(IRS.GT.0) THEN
+        KPLIB=IPISO(ISO) ! set ISO-th isotope
+        CALL LCMGET(KPLIB,'NTOT0',TOTAL(1,IRS))
+        CALL LCMGET(KPLIB,'SIGS00',SIGS0(1,IRS))
+        DO IGRP=IGRMIN,IGRMAX
+*         Compute a ST flux for the homogeneous equivalent medium.
+          PHGAR(IGRP,IRS)=MAX(0.0,SN(IGRP,ISO)/(SN(IGRP,ISO)+
+     1    (TOTAL(IGRP,IRS)-SIGS0(IGRP,IRS))))
+        ENDDO
+        IF(ITRANC.NE.0) CALL LCMGET(KPLIB,'TRANC',TRANC(1,IRS))
+      ENDIF
+   40 CONTINUE
+*----
+*  COMPUTE THE MERGED VOLUMES.
+*----
+      NALBP=0
+      LHOMOG=.TRUE.
+      VST(:NBNRS)=0.0
+      DO 50 I=1,NREG
+      IBM=MAT(I)
+      IF(IBM.EQ.0) GO TO 50
+      IND=IRES(IBM)
+      IF(IND.EQ.0) THEN
+         LHOMOG=.FALSE.
+      ELSE
+         VST(IND)=VST(IND)+VOL(I)
+      ENDIF
+   50 CONTINUE
+      IF(LHOMOG.AND.(NBNRS.EQ.1)) GO TO 260
+      IF(IMPX.GE.3) WRITE(6,'(37H TONSPH: SPH FACTOR CALCULATION (NBNR,
+     1 2HS=,I5,1H)/)') NBNRS
+*----
+*  SET THE MIXTURE-DEPENDENT MACROSCOPIC XS.
+*----
+      FUNKNO(:NUN,:NGRO)=0.0
+      SUNKNO(:NUN,:NGRO)=0.0
+      NPSYS(:NGRO)=0
+      CALL LCMSIX(IPLIB,'SHIBA',1)
+      JPLIB=LCMLID(IPLIB,'GROUP',NGRO)
+      DO 110 IGRP=IGRMIN,IGRMAX
+        NPSYS(IGRP)=IGRP
+*
+*       COMPUTE THE LIGHT AND RESONANT COMPONENTS OF THE MACROSCOPIC
+*       CROSS SECTIONS IN EACH RESONANT MIXTURE.
+        DO 70 IBM=1,NBM
+        SIG0(IBM,IGRP)=0.0
+        SIG1(IBM,IGRP)=0.0
+        SIG3(IBM,IGRP)=SIGT3(IBM,IGRP)
+   70   CONTINUE
+        DO 80 ISO=1,NBISO
+        IRS=ISONR(ISO)
+        IF(IRS.GT.0) THEN
+          IBM=MIX(ISO)
+          FLUX(IBM,IGRP)=PHGAR(IGRP,IRS)
+          SIGT2(IBM,IGRP)=SIGT2(IBM,IGRP)-TOTAL(IGRP,IRS)*DEN(ISO)
+          SIG0(IBM,IGRP)=TOTAL(IGRP,IRS)*DEN(ISO)
+          SIG1(IBM,IGRP)=SIGS0(IGRP,IRS)*DEN(ISO)
+          IF(ITRANC.NE.0) THEN
+            SIG3(IBM,IGRP)=SIGT3(IBM,IGRP)-TRANC(IGRP,IRS)*DEN(ISO)
+          ENDIF
+        ENDIF
+   80   CONTINUE
+        IF(IMPX.GE.10) THEN
+          WRITE (6,400) IGRP,(SIG0(I,IGRP),I=1,NBM)
+          WRITE (6,410) IGRP,(SIG1(I,IGRP),I=1,NBM)
+          WRITE (6,420) IGRP,(SIGT2(I,IGRP),I=1,NBM)
+          WRITE (6,430) IGRP,(FLUX(I,IGRP),I=1,NBM)
+        ENDIF
+*----
+*  COMPUTE THE SOURCES.
+*----
+        SIGG(0)=0.0
+        DO 90 IBM=1,NBM
+        SIGG(IBM)=SIGT2(IBM,IGRP)
+        IF(IRES(IBM).GT.0) THEN
+          SIGG(IBM)=SIGG(IBM)+FLUX(IBM,IGRP)*(SIG1(IBM,IGRP)-
+     >    SIG0(IBM,IGRP))
+          IF(.NOT.LHOMOG) SIGG(IBM)=SIGG(IBM)-FLUX(IBM,IGRP)*
+     >    SIGT2(IBM,IGRP)
+        ENDIF
+   90   CONTINUE
+        SUNKNO(:NUN,IGRP)=0.0
+        CALL DOORS(CDOOR,IPTRK,NBM,0,NUN,SIGG,SUNKNO(1,IGRP))
+*
+        IF(NPSYS(IGRP).NE.0) THEN
+          ICPIJ=ICPIJ+1
+          SIGTXS(0)=0.0
+          SIGS0X(0)=0.0
+          DO 100 IBM=1,NBM
+          IND=IRES(IBM)
+          IF((ITRANC.NE.0).AND.(IND.EQ.0)) THEN
+            SIGTXS(IBM)=SIGT2(IBM,IGRP)-SIG3(IBM,IGRP)
+          ELSE
+            SIGTXS(IBM)=SIGT2(IBM,IGRP)
+          ENDIF
+          IF(IND.EQ.0) THEN
+*           REMOVE TRANSPORT CORRECTION.
+            IF(ITRANC.NE.0) THEN
+              SIGS0X(IBM)=-SIG3(IBM,IGRP)
+            ELSE
+              SIGS0X(IBM)=0.0
+            ENDIF
+          ELSE
+*           BELL ACCELERATION.
+            SIGTXS(IBM)=SIGTXS(IBM)+SIG0(IBM,IGRP)
+            SIGS0X(IBM)=SIGTXS(IBM)
+            IF(LHOMOG) SIGS0X(IBM)=SIGS0X(IBM)-SIGT2(IBM,IGRP)
+          ENDIF
+  100     CONTINUE
+          KPLIB=LCMDIL(JPLIB,IGRP)
+          CALL LCMPUT(KPLIB,'DRAGON-TXSC',NBM+1,2,SIGTXS)
+          CALL LCMPUT(KPLIB,'DRAGON-S0XSC',NBM+1,2,SIGS0X)
+        ENDIF
+  110 CONTINUE
+*----
+*  SOLVE FOR THE FLUX USING DIRECT SELF-SHIELDED CROSS SECTIONS
+*----
+      CALL KDRCPU(TKA)
+      ISTRM=1
+      NANI=1
+      NW=0
+      KNORM=1
+      IMPY=MAX(0,IMPX-3)
+      IF(IPHASE.EQ.1) THEN
+*        USE A NATIVE DOOR.
+         CALL DOORAV(CDOOR,JPLIB,NPSYS,IPTRK,IFTRAK,IMPY,NGRO,NREG,
+     1   NBM,NANI,NW,MAT,VOL,KNORM,LEAKSW,TITR,NALBP,ISTRM)
+      ELSE IF(IPHASE.EQ.2) THEN
+*        USE A COLLISION PROBABILITY DOOR.
+         IPIJK=1
+         ITPIJ=1
+         CALL DOORPV(CDOOR,JPLIB,NPSYS,IPTRK,IFTRAK,IMPY,NGRO,NREG,
+     1   NBM,NANI,MAT,VOL,KNORM,IPIJK,LEAKSW,ITPIJ,.FALSE.,TITR,NALBP)
+      ENDIF
+      CALL KDRCPU(TKB)
+      TK3=TK3+(TKB-TKA)
+      CALL KDRCPU(TKA)
+      IDIR=0
+      LEXAC=.FALSE.
+      IPMACR=C_NULL_PTR
+      IPSOU=C_NULL_PTR
+      REBFLG=.FALSE.
+      CALL DOORFV(CDOOR,JPLIB,NPSYS,IPTRK,IFTRAK,IMPX,NGRO,NBM,IDIR,
+     1 NREG,NUN,IPHASE,LEXAC,MAT,VOL,KEYFLX,TITR,SUNKNO,FUNKNO,IPMACR,
+     2 IPSOU,REBFLG)
+      CALL LCMSIX(IPLIB,' ',2)
+      TK4=TK4+(TKB-TKA)
+*----
+*  HOMOGENIZE THE FLUX
+*----
+      DO 150 IGRP=IGRMIN,IGRMAX
+      IF(NPSYS(IGRP).NE.0) THEN
+        FLNEW(:NBNRS)=0.0
+        DO 120 I=1,NREG
+        IF(MAT(I).EQ.0) GO TO 120
+        IND=IRES(MAT(I))
+        IF(IND.GT.0) FLNEW(IND)=FLNEW(IND)+FUNKNO(KEYFLX(I),IGRP)*VOL(I)
+  120   CONTINUE
+        DO 130 IND=1,NBNRS
+        FLNEW(IND)=FLNEW(IND)/VST(IND)
+  130   CONTINUE
+*----
+*  SPH FACTOR CONTROL
+*----
+        DO 140 IBM=1,NBM
+        IND=IRES(IBM)
+        IF(IND.GT.0) THEN
+          SPHNEW=PHGAR(IGRP,IND)/FLNEW(IND)
+          LPROB=(SPHNEW.LE.0.).OR.(SPHNEW.GT.1.).OR.(FLNEW(IND).LT.0.05)
+          IF(LPROB) SPHNEW=1.0
+          SPH(IBM,IGRP)=SPHNEW
+        ENDIF
+  140   CONTINUE
+      ENDIF
+  150 CONTINUE
+*----
+*  SPH CORRECTION OF THE MICROLIB
+*----
+      CALL LCMGET(IPLIB,'STATE-VECTOR',ISTATE)
+      NL=ISTATE(4)
+      NED=ISTATE(13)
+      NDEL=ISTATE(19)
+      DO 160 ISO=1,NBISO
+      MASKI(ISO)=(ISONR(ISO).GT.0)
+  160 CONTINUE
+      CALL TONCMI(IPLIB,IMPX,NBM,NBISO,NGRO,NL,NED,NDEL,MASKI,SPH)
+      IF(IMPX.GT.3) THEN
+        DO 170 IGRP=IGRMIN,IGRMAX
+        WRITE (6,440) IGRP,(SPH(IBM,IGRP),IBM=1,NBM)
+  170   CONTINUE
+      ENDIF
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+  260 DEALLOCATE(IPISO)
+      DEALLOCATE(MASKI)
+      DEALLOCATE(PHGAR,TRANC,SIGS0,TOTAL,SIG3,SIG1,SIG0,FLUX,FUNKNO,
+     1 SUNKNO,FLNEW,SIGG,SIGS0X,SIGTXS,VST)
+      DEALLOCATE(NPSYS,ISONR,IRES)
+      RETURN
+  400 FORMAT(/51H TOTAL MACROSCOPIC CROSS SECTIONS OF THE RESONANT M,
+     1 31HATERIALS IN EACH MIXTURE (GROUP,I5,2H):/(1X,1P,11E11.3))
+  410 FORMAT(/51H SCATTERING MACROSCOPIC CROSS SECTIONS OF THE OTHER,
+     1 33H MATERIALS IN EACH MIXTURE (GROUP,I5,2H):/(1X,1P,11E11.3))
+  420 FORMAT(/51H TOTAL MACROSCOPIC CROSS SECTIONS OF THE OTHER MATE,
+     1 28HRIALS IN EACH MIXTURE (GROUP,I5,2H):/(1X,1P,11E11.3))
+  430 FORMAT(/19H TABSN3 FLUX (GROUP,I5,2H):/(1X,1P,11E11.3))
+  440 FORMAT(/19H SPH FACTORS (GROUP,I5,2H):/(1X,1P,11E11.3))
+      END