Initial commit from Polytechnique Montreal

1 files changed, 458 insertions, 0 deletions

diff --git a/Dragon/src/USSRSE.f b/Dragon/src/USSRSE.f
new file mode 100644
index 0000000..aa0f60e
--- /dev/null
+++ b/Dragon/src/USSRSE.f
@@ -0,0 +1,458 @@
+*DECK USSRSE
+      SUBROUTINE USSRSE(IPTRK,IPLIB,IPLI0,IFTRAK,NREG,NUN,NBMIX,NBISO,
+     1 NIRES,NL,NED,NDEL,ISONAM,ISOBIS,HCAL,MAT,VOL,KEYFLX,CDOOR,
+     2 LEAKSW,IMPX,DEN,MIX,IAPT,IPHASE,NGRP,IGRMIN,IGRMAX,NBNRS,IREX,
+     3 TITR,ICORR,MAXST,GOLD,UNGAR,PHGAR,STGAR,SFGAR,SSGAR,S0GAR,SAGAR,
+     4 SDGAR,MASKG,SIGGAR)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Compute the self-shielded cross sections in each energy group using
+* the resonance spectrum expansion 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.
+* NL      number of legendre orders required in the calculation
+*         (NL=1 or higher).
+* NED     number of extra vector edits.
+* NDEL    number of delayed neutron precursor groups.
+* ISONAM  alias name of isotopes in IPLIB.
+* ISOBIS  alias name of isotopes in IPLI0.
+* HCAL    name of the self-shielding calculation.
+* 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.
+* IGRMAX  most thermal group where the self-shielding is applied.
+* NBNRS   number of correlated fuel regions. Note that NBNRS=max(IREX).
+* IREX    fuel region index assigned to each mixture. Equal to zero
+*         in non-resonant mixtures or in mixtures not used.
+* TITR    title.
+* ICORR   mutual resonance shielding flag (=1 to suppress the model
+*         in cases it is required in LIB operator).
+* MAXST   maximum number of fixed point iterations for the ST scattering
+*         source.
+*
+*Parameters: output
+* GOLD    Goldstein-Cohen parameters.
+* UNGAR   averaged flux unknowns.
+* PHGAR   averaged fluxes in correlated fuel regions.
+* STGAR   microscopic self-shielded total x-s.
+* SFGAR   microscopic self-shielded fission x-s.
+* SSGAR   microscopic self-shielded scattering x-s.
+* S0GAR   microscopic transfer scattering xs (isotope,secondary,
+*         primary).
+* SAGAR   microscopic self-shielded additional xs.
+* SDGAR   microscopic self-shielded delayed nu-sigf xs.
+* MASKG   energy group mask pointing on self-shielded groups.
+* SIGGAR  macroscopic x-s of the non-resonant isotopes in each mixture:
+*         (*,*,*,1) total; (*,*,*,2) transport correction; 
+*         (*,*,*,3) P0 scattering; (*,*,*,4) flux times P0 scattering.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPTRK,IPLIB,IPLI0
+      INTEGER IFTRAK,NREG,NUN,NBMIX,NBISO,NIRES,NL,NED,NDEL,
+     1 ISONAM(3,NBISO),ISOBIS(3,NBISO),MAT(NREG),KEYFLX(NREG),IMPX,
+     2 MIX(NBISO),IAPT(NBISO),IPHASE,NGRP,IGRMIN,IGRMAX,NBNRS,
+     3 IREX(NBMIX),ICORR,MAXST
+      REAL VOL(NREG),DEN(NBISO),GOLD(NIRES,NGRP),UNGAR(NUN,NIRES,NGRP),
+     1 PHGAR(NBNRS,NIRES,NGRP),STGAR(NBNRS,NIRES,NGRP),
+     2 SFGAR(NBNRS,NIRES,NGRP),SSGAR(NBNRS,NIRES,NL,NGRP),
+     3 S0GAR(NBNRS,NIRES,NL,NGRP,NGRP),SAGAR(NBNRS,NIRES,NED,NGRP),
+     4 SDGAR(NBNRS,NIRES,NDEL,NGRP),SIGGAR(NBMIX,0:NIRES,NGRP,4)
+      LOGICAL LEAKSW,MASKG(NGRP,NIRES)
+      CHARACTER HCAL*12,CDOOR*12,TITR*72
+*----
+*  LOCAL VARIABLES
+*----
+      TYPE(C_PTR) IPP,KPLIB,LPLIB,MPLIB,JPLI0,KPLI0,IOFSET
+      LOGICAL LLIB
+      PARAMETER (MAXED=50,MAXNOR=20)
+      CHARACTER TEXT12*12,HVECT(MAXED)*8,CBDPNM*12,HSMG*131
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      TYPE(C_PTR), ALLOCATABLE, DIMENSION(:) :: IPPT1,IPISO1,IPISO2
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: IWRK
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) :: NOR,IPPT2,ISM
+      REAL, ALLOCATABLE, DIMENSION(:) :: GAS,GA1,VOLMER,DELTAU,GOLD2
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: GA2,CONR,XFLUX
+      TYPE VECTOR_ARRAY
+        DOUBLE PRECISION, POINTER, DIMENSION(:) :: VECTOR
+      END TYPE VECTOR_ARRAY
+      TYPE(VECTOR_ARRAY), ALLOCATABLE, DIMENSION(:,:) :: GAMMA_V
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(IPPT1(NIRES))
+      ALLOCATE(NOR(NIRES,NGRP),IPPT2(NIRES,5),IWRK(NGRP),ISM(2,NL))
+      ALLOCATE(GAS(NGRP),GA1(NGRP),GA2(NGRP,NGRP),CONR(NBNRS,NIRES),
+     1 VOLMER(0:NBNRS),DELTAU(NGRP),GAMMA_V(NGRP,NIRES))
+      ALLOCATE(IPISO1(NBISO),IPISO2(NBISO))
+*
+      CALL KDRCPU(TK1)
+      PHGAR(:NBNRS,:NIRES,:NGRP)=1.0
+      NOR(:NIRES,1)=-1
+*
+      IF(NED.GT.0) THEN
+        IF(NED.GT.MAXED) CALL XABORT('USSRSE: INVALID VALUE OF MAXED.')
+        CALL LCMGTC(IPLIB,'ADDXSNAME-P0',8,NED,HVECT)
+      ENDIF
+*
+      CALL LIBIPS(IPLIB,NBISO,IPISO1)
+      CALL LIBIPS(IPLI0,NBISO,IPISO2)
+      SIGGAR(:NBMIX,0:NIRES,:NGRP,:4)=0.0
+      DO 190 ISO=1,NBISO
+      IBM=MIX(ISO)
+      DO 30 I=1,NREG
+      IF(MAT(I).EQ.IBM) GO TO 35
+   30 CONTINUE
+      GO TO 190
+   35 IRES=IAPT(ISO)
+      DENN=DEN(ISO)
+      JRES=IRES
+      IF(IRES.EQ.NIRES+1) JRES=0
+*----
+*  RECOVER INFINITE DILUTION OR SELF-SHIELDED CROSS SECTIONS AND
+*  COMPUTE OUT-OF-FUEL MACROSCOPIC CROSS SECTIONS.
+*----
+      KPLI0=IPISO2(ISO) ! set ISO-th 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 USSRSE: ISOTOPE '',3A4,7H'' (ISO=,I8,5H) IS ,
+     1   39HNOT AVAILABLE IN THE ORIGINAL MICROLIB.)') (ISONAM(I0,ISO),
+     2   I0=1,3),ISO
+         CALL XABORT(HSMG)
+      ENDIF
+      IF((.NOT.LLIB).AND.(IMPX.GT.2)) WRITE(6,'(/18H USSRSE: RECOVER I,
+     1 8HSOTOPE '',3A4,23H'' FROM THE NEW LIBRARY.)') (ISOBIS(I0,ISO),
+     2 I0=1,3)
+      IF((DENN.NE.0.0).AND.(IBM.NE.0)) THEN
+         CALL LCMLEN(IPP,'NTOT0',ILENGT,ITYLCM)
+         IF(ILENGT.NE.NGRP) THEN
+           CALL LCMLIB(IPP)
+           CALL XABORT('USSRSE: INVALID X-SECTIONS.')
+         ENDIF
+         CALL LCMGET(IPP,'NTOT0',GA1)
+         CALL XDRLGS(IPP,-1,IMPX,0,0,1,NGRP,GAS,GA2,ITYPRO)
+         DO 40 IGRP=1,NGRP
+         SIGGAR(IBM,JRES,IGRP,1)=SIGGAR(IBM,JRES,IGRP,1)+DENN*GA1(IGRP)
+   40    CONTINUE
+         CALL LCMGET(IPP,'SIGS00',GA1)
+         CALL LCMLEN(IPP,'NWT0',ILENGT,ITYLCM)
+         IF(ILENGT.GT.0) THEN
+            CALL LCMGET(IPP,'NWT0',GAS)
+         ELSE
+            GAS(:NGRP)=1.0
+         ENDIF
+         DO 45 IGRP=1,NGRP
+         SIGGAR(IBM,JRES,IGRP,3)=SIGGAR(IBM,JRES,IGRP,3)+DENN*GA1(IGRP)
+         DO 44 JGRP=1,IGRP
+         SIGGAR(IBM,JRES,IGRP,4)=SIGGAR(IBM,JRES,IGRP,4)+DENN*
+     1   GA2(IGRP,JGRP)*GAS(JGRP)
+   44    CONTINUE
+   45    CONTINUE
+         CALL LCMLEN(IPP,'TRANC',ILENGT,ITYLCM)
+         IF(ILENGT.GT.0) THEN
+            CALL LCMGET(IPP,'TRANC',GA1)
+         ELSE
+            GA1(:NGRP)=0.0
+         ENDIF
+         DO 50 IGRP=1,NGRP
+         SIGGAR(IBM,JRES,IGRP,2)=SIGGAR(IBM,JRES,IGRP,2)+DENN*GA1(IGRP)
+   50    CONTINUE
+      ENDIF
+      CALL LCMGET(IPLI0,'DELTAU',DELTAU)
+*----
+*  RECOVER PROBABILITY TABLE INFORMATION.
+*----
+      IF((IRES.GT.0).AND.(IRES.LE.NIRES)) THEN
+         IF(NOR(IRES,1).EQ.-1) THEN
+            KPLIB=IPISO1(ISO) ! set ISO-th isotope
+*
+*           RECOVER INFINITE DILUTION VALUES.
+            CALL LCMGET(KPLIB,'NTOT0',GAS)
+            DO 55 IGRP=1,NGRP
+            STGAR(:NBNRS,IRES,IGRP)=0.0
+            STGAR(:NBNRS,IRES,IGRP)=GAS(IGRP)
+            SFGAR(:NBNRS,IRES,IGRP)=0.0
+            SAGAR(:NBNRS,IRES,:NED,IGRP)=0.0
+            SDGAR(:NBNRS,IRES,:NDEL,IGRP)=0.0
+   55       CONTINUE
+            CALL LCMLEN(KPLIB,'NUSIGF',ILENGT,ITYLCM)
+            IF(ILENGT.GT.0) THEN
+               CALL LCMGET(KPLIB,'NUSIGF',GAS)
+               DO 60 IGRP=1,NGRP
+               SFGAR(:NBNRS,IRES,IGRP)=GAS(IGRP)
+   60          CONTINUE
+            ENDIF
+            DO 80 IL=1,NL
+            CALL XDRLGS(KPLIB,-1,IMPX,IL-1,IL-1,1,NGRP,GAS,GA2,ITYPRO)
+*           JG IS THE SECONDARY GROUP.
+            DO 72 IGRP=1,NGRP
+            SSGAR(:NBNRS,IRES,IL,IGRP)=GAS(IGRP)
+            DO 70 JGRP=1,NGRP
+            S0GAR(:NBNRS,IRES,IL,JGRP,IGRP)=GA2(JGRP,IGRP)
+   70       CONTINUE
+   72       CONTINUE
+   80       CONTINUE
+            DO 110 IED=1,NED
+            CALL LCMLEN(KPLIB,HVECT(IED),ILENGT,ITYLCM)
+            IF(ILENGT.GT.0) THEN
+               CALL LCMGET(KPLIB,HVECT(IED),GAS)
+               DO 105 IGRP=1,NGRP
+               SAGAR(:NBNRS,IRES,IED,IGRP)=GAS(IGRP)
+  105          CONTINUE
+            ENDIF
+  110       CONTINUE
+            DO 130 IDEL=1,NDEL
+            WRITE(TEXT12,'(6HNUSIGF,I2.2)') IDEL
+            CALL LCMLEN(KPLIB,TEXT12,ILENGT,ITYLCM)
+            IF(ILENGT.GT.0) THEN
+               CALL LCMGET(KPLIB,TEXT12,GAS)
+               DO 125 IGRP=1,NGRP
+               SDGAR(:NBNRS,IRES,IDEL,IGRP)=GAS(IGRP)
+  125          CONTINUE
+            ENDIF
+  130       CONTINUE
+*
+            GOLD(IRES,:NGRP)=1.0
+            NOR(IRES,:NGRP)=0
+            CALL LCMLEN(KPLIB,'NGOLD',ILENGT,ITYLCM)
+            IF(ILENGT.GT.0) THEN
+              ALLOCATE(GOLD2(NGRP))
+              CALL LCMGET(KPLIB,'NGOLD',GOLD2)
+              GOLD(IRES,IGRMIN:IGRMAX)=GOLD2(IGRMIN:IGRMAX)
+              DEALLOCATE(GOLD2)
+            ENDIF
+            CALL LCMLEN(KPLIB,'PT-TABLE',ILENGT,ITYLCM)
+            IF(ILENGT.EQ.-1) THEN
+              CALL LCMSIX(KPLIB,'PT-TABLE',1)
+                CALL LCMGET(KPLIB,'NOR',IWRK)
+                CALL LCMLEN(KPLIB,'GROUP-RSE',ILENGT,ITYLCM)
+                IF(ILENGT.GT.0) THEN
+                  LPLIB=LCMGID(KPLIB,'GROUP-RSE')
+                  DO IGRP=IGRMIN,IGRMAX
+                    MI=NOR(IRES,IGRP)
+                    CALL LCMLEL(LPLIB,IGRP,ILENGT,ITYLCM)
+                    IF(ILENGT.EQ.0) GOLD(IRES,IGRP)=1.0
+                    NOR(IRES,IGRP)=IWRK(IGRP)
+                    CALL LCMLEL(LPLIB,IGRP,ILENGT,ITYLCM)
+                    IF(GOLD(IRES,IGRP).EQ.-1001.0) THEN
+                      MPLIB=LCMGIL(LPLIB,IGRP)
+                      CALL LCMGPD(MPLIB,'GAMMA_V',IOFSET)
+                      CALL C_F_POINTER(IOFSET,GAMMA_V(IGRP,IRES)%VECTOR,
+     1                (/MI/))
+                    ENDIF
+                  ENDDO
+                ENDIF
+              CALL LCMSIX(KPLIB,' ',2)
+            ENDIF
+         ENDIF
+      ENDIF
+  190 CONTINUE
+      CALL KDRCPU(TK2)
+      IF(IMPX.GT.1) WRITE(6,'(/34H USSRSE: CPU TIME SPENT TO RECOVER,
+     1 23H INFINITE-DILUTION XS =,F8.1,8H SECOND./)') TK2-TK1
+*
+      CALL KDRCPU(TK1)
+      TK4=0.0
+      TK5=0.0
+      ICPIJ=0
+*----
+*  COMPUTE THE MERGED VOLUMES AND NUMBER DENSITIES.
+*----
+      VOLMER(0:NBNRS)=0.0
+      DO 210 I=1,NREG
+      IBM=MAT(I)
+      IF(IBM.GT.0) VOLMER(IREX(IBM))=VOLMER(IREX(IBM))+VOL(I)
+  210 CONTINUE
+      CONR(:NBNRS,:NIRES)=0.0
+      DO 240 ISO=1,NBISO
+      JRES=IAPT(ISO)
+      IF((JRES.GT.0).AND.(JRES.LE.NIRES)) THEN
+         DENN=DEN(ISO)
+         DO 230 IREG=1,NREG
+         IBM=MAT(IREG)
+         IF(MIX(ISO).EQ.IBM) THEN
+            IND=IREX(IBM)
+            IF(IND.EQ.0) CALL XABORT('USSRSE: IREX FAILURE.')
+            CONR(IND,JRES)=CONR(IND,JRES)+DENN*VOL(IREG)/VOLMER(IND)
+         ENDIF
+  230    CONTINUE
+      ENDIF
+  240 CONTINUE
+*----
+*  RECOVER POSITION OF PROBABILITY TABLES AND NAME OF RESONANT ISOTOPE.
+*----
+      DO 270 IRES=1,NIRES
+      ISOT=0
+      DO 250 JSOT=1,NBISO
+      IF(IAPT(JSOT).EQ.IRES) THEN
+         ISOT=JSOT
+         GO TO 260
+      ENDIF
+  250 CONTINUE
+      CALL XABORT('USSRSE: UNABLE TO FIND A RESONANT ISOTOPE.')
+  260 KPLIB=IPISO1(ISOT) ! set ISOT-th isotope
+      CALL LCMLEN(KPLIB,'PT-TABLE',ILONG,ITYLCM)
+      IF(ILONG.EQ.0) CALL XABORT('USSRSE: BUG1.')
+      CALL LCMSIX(KPLIB,'PT-TABLE',1)
+      CALL LCMGET(KPLIB,'NDEL',NDEL0)
+      IF(NDEL0.GT.NDEL) CALL XABORT('USSRSE: NDEL OVERFLOW.')
+      IPPT1(IRES)=KPLIB
+      CALL LCMSIX(KPLIB,' ',2)
+      IPPT2(IRES,1)=IREX(MIX(ISOT))
+      IPPT2(IRES,2)=ISONAM(1,ISOT)
+      IPPT2(IRES,3)=ISONAM(2,ISOT)
+      IPPT2(IRES,4)=ISONAM(3,ISOT)
+      IPPT2(IRES,5)=NDEL0
+      IF(IPPT2(IRES,1).LE.0) CALL XABORT('USSRSE: BUG3.')
+  270 CONTINUE
+*----
+*  DETERMINE WHICH GROUPS ARE SELF-SHIELDED.
+*----
+      DO 290 IGRP=1,NGRP
+      DO 280 IRES=1,NIRES
+      MASKG(IGRP,IRES)=((IGRP.GE.IGRMIN).AND.(IGRP.LE.IGRMAX).AND.
+     1 (NOR(IRES,IGRP).GT.0))
+  280 CONTINUE
+  290 CONTINUE
+*----
+*  INITIALIZATION OF THE MULTIBAND FLUXES AND SOURCES.
+*----
+      CALL LCMSIX(IPLI0,'SHIBA_SG',1)
+      CALL LCMSIX(IPLI0,HCAL,1)
+      DO 310 IRES=1,NIRES
+      WRITE(CBDPNM,'(3HCOR,I4.4,1H/,I4.4)') IRES,NIRES
+      CALL LCMSIX(IPLI0,CBDPNM,1)
+      JPLI0=LCMLID(IPLI0,'NWT0-PT',NGRP)
+      DO 300 IGRP=1,NGRP
+      IF(MASKG(IGRP,IRES)) THEN
+         CALL LCMLEL(JPLI0,IGRP,ILENGT,ITYLCM)
+         IF(ILENGT.EQ.0) THEN
+            MI=NOR(IRES,IGRP)
+            ALLOCATE(XFLUX(NBNRS,MI))
+            IF(GOLD(IRES,IGRP).EQ.1.0) THEN
+              XFLUX(:NBNRS,:MI)=1.0
+            ELSE IF(GOLD(IRES,IGRP).EQ.-1001.0) THEN
+              DO IM=1,MI
+                XFLUX(:NBNRS,IM)=REAL(GAMMA_V(IGRP,IRES)%VECTOR(IM))
+              ENDDO
+            ENDIF
+            CALL LCMPDL(JPLI0,IGRP,NBNRS*MI,2,XFLUX)
+            DEALLOCATE(XFLUX)
+         ENDIF
+      ENDIF
+  300 CONTINUE
+      CALL LCMSIX(IPLI0,' ',2)
+  310 CONTINUE
+*
+      DO 340 IRES=1,NIRES
+      DO 330 IGRP=1,NGRP
+      IF(MASKG(IGRP,IRES)) ICPIJ=ICPIJ+NOR(IRES,IGRP)
+  330 CONTINUE
+      CALL KDRCPU(TKA)
+*----
+*  ITERATIVE APPROACH FOR THE HELIOS/WIMS-7 METHOD.
+*----
+      MAX_R=12
+      CALL USSIT1(MAX_R,NGRP,MASKG(1,IRES),IRES,IPLI0,IPTRK,IFTRAK,
+     1 CDOOR,IMPX,NBMIX,NREG,NUN,NL,IPHASE,MAXST,MAT,VOL,KEYFLX,LEAKSW,
+     2 IREX,SIGGAR,TITR,NIRES,NBNRS,NOR,CONR,GOLD,IPPT1,IPPT2,STGAR,
+     3 SSGAR,VOLMER,UNGAR)
+*----
+*  ITERATIVE APPROACH FOR THE RESONANCE SPECTRUM EXPANSION METHOD.
+*----
+      CALL USSIT3(MAXNOR,NGRP,MASKG(1,IRES),IRES,IPLI0,IPTRK,IFTRAK,
+     1 CDOOR,IMPX,NBMIX,NREG,NUN,IPHASE,MAXST,MAT,VOL,KEYFLX,LEAKSW,
+     2 IREX,SIGGAR,TITR,ICORR,NIRES,NBNRS,CONR,GOLD,IPPT1,IPPT2,
+     3 VOLMER,UNGAR)
+      CALL KDRCPU(TKB)
+      TK4=TK4+(TKB-TKA)
+  340 CONTINUE
+*----
+*  COMPUTE THE SELF-SHIELDED REACTION RATES.
+*----
+      CALL USSIT4(MAXNOR,IPLI0,IPPT1,IPPT2,NGRP,NIRES,NBNRS,NL,NED,
+     1 NDEL,PHGAR,STGAR,SFGAR,SSGAR,S0GAR,SAGAR,SDGAR)
+      CALL LCMSIX(IPLI0,' ',2)
+      CALL LCMSIX(IPLI0,' ',2)
+      CALL LCMVAL(IPLI0,' ')
+*----
+*  RESET MASKG FOR SPH CALCULATION IN SMALL LETHARGY WIDTH GROUPS.
+*----
+      IF(NIRES.GT.1) THEN
+        DO 360 IGRP=1,NGRP
+        DO 350 IRES=1,NIRES
+        IF(MASKG(IGRP,IRES)) THEN
+          MASKG(IGRP,IRES)=.NOT.(GOLD(IRES,IGRP).EQ.-1001.)
+          IF(DELTAU(IGRP).GT.0.1) MASKG(IGRP,IRES)=.TRUE.
+        ENDIF
+  350   CONTINUE
+  360   CONTINUE
+      ENDIF
+      CALL KDRCPU(TK2)
+      IF(IMPX.GT.1) WRITE(6,'(/34H USSRSE: CPU TIME SPENT TO COMPUTE,
+     1 31H SELF-SHIELDED REACTION RATES =,F8.1,19H SECOND, INCLUDING:
+     2 /9X,F8.1,46H SECOND TO BUILD/SOLVE SUBGROUP MATRIX SYSTEM;/9X,
+     4 F8.1,38H SECOND TO COMPUTE THE REACTION RATES./9X,9HNUMBER OF,
+     5 23H ASSEMBLY DOORS CALLS =,I5,1H.)') TK2-TK1,TK4,TK5,ICPIJ
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(IPISO2,IPISO1)
+      DEALLOCATE(GAMMA_V,DELTAU,VOLMER,CONR,GA2,GA1,GAS)
+      DEALLOCATE(ISM,IWRK,IPPT2,NOR)
+      DEALLOCATE(IPPT1)
+      RETURN
+      END