Initial commit from Polytechnique Montreal

1 files changed, 383 insertions, 0 deletions

diff --git a/Dragon/src/USSEXC.f b/Dragon/src/USSEXC.f
new file mode 100644
index 0000000..92b5db4
--- /dev/null
+++ b/Dragon/src/USSEXC.f
@@ -0,0 +1,383 @@
+*DECK USSEXC
+      SUBROUTINE USSEXC(MAXNOR,CDOOR,IPLI0,IPTRK,IFTRAK,IMPX,NGRP,IGRP,
+     1 IASM,NBMIX,NREG,NUN,NL,IPHASE,MAT,VOL,KEYFLX,IREX,SIGGAR,TITR,
+     2 NIRES,IRES,NBNRS,NOR,CONR,IPPT1,IPPT2,STGAR,SSGAR,VOLMER,XFLUX,
+     3 UNGAR)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Solution of the flux for the subgroup projection method (SPM) using
+* the response matrix method. This is a non-iterative approach which is
+* useful in exceptional cases where the fixed-point approach fails.
+*
+*Copyright:
+* Copyright (C) 2010 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
+* MAXNOR  maximum order of the probability tables (PT).
+* CDOOR   name of the geometry/solution operator.
+* IPLI0   pointer to the internal microscopic cross section library
+*         builded by the self-shielding module.
+* IPTRK   pointer to the tracking (L_TRACK signature).
+* IFTRAK  file unit number used to store the tracks.
+* IMPX    print flag (equal to zero for no print).
+* NGRP    number of energy groups.
+* IGRP    index of energy group being processed.
+* IASM    offset of information computed by DOORAV or DOORPV.
+* NBMIX   number of mixtures in the internal library.
+* NREG    number of regions.
+* NUN     number of unknowns in the flux or source vector in one
+*         energy group and one band.
+* NL      number of Legendre orders required in the calculation
+*         (NL=1 or higher).
+* IPHASE  type of flux solution (=1 use a native flux solution door;
+*         =2 use collision probabilities).
+* MAT     index-number of the mixture type assigned to each volume.
+* VOL     volumes.
+* KEYFLX  pointers of fluxes in unknown vector.
+* IREX    fuel region index assigned to each mixture. Equal to zero
+*         in non-resonant mixtures or in mixtures not used.
+* 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.
+* TITR    title.
+* NIRES   exact number of correlated resonant isotopes.
+* IRES    index of the resonant isotope being processed.
+* NBNRS   number of correlated fuel regions.
+* NOR     exact order of the probability table.
+* CONR    number density of the resonant isotopes.
+* IPPT1   pointer to LCM directory of each resonant isotope.
+* IPPT2   information related to each resonant isotope:
+*         IPPT2(:,1) index of a resonant region (used with infinite
+*         dilution case);
+*         IPPT2(:,2:4) alias name of resonant isotope.
+* STGAR   averaged microscopic total xs in resonant region.
+* SSGAR   averaged microscopic scattering xs in resonant region.
+* VOLMER  volumes of the resonant and non-resonant regions.
+*
+*Parameters: input/output
+* XFLUX   subgroup flux.
+*
+*Parameters: output
+* UNGAR   averaged flux unknowns.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+      USE DOORS_MOD
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) IPLI0,IPTRK,IPPT1(NIRES)
+      INTEGER MAXNOR,IFTRAK,IMPX,NGRP,IGRP,IASM,NBMIX,NREG,NUN,NL,
+     1 IPHASE,MAT(NREG),KEYFLX(NREG),IREX(NBMIX),NIRES,IRES,NBNRS,
+     2 NOR(NIRES,NGRP),IPPT2(NIRES,4)
+      REAL VOL(NREG),SIGGAR(NBMIX,0:NIRES,NGRP,4),CONR(NBNRS,NIRES),
+     1 STGAR(NBNRS,NIRES,NGRP),SSGAR(NBNRS,NIRES,NL,NGRP),
+     2 VOLMER(0:NBNRS),XFLUX(NBNRS,MAXNOR),UNGAR(NUN,NIRES,NGRP)
+      CHARACTER CDOOR*12,TITR*72
+*----
+*  LOCAL VARIABLES
+*----
+      TYPE(C_PTR) IPSYS,KPSYS,JPLIB,KPLIB,IPMACR,IPSOU
+      LOGICAL EMPTY,LCM,LEXAC,REBFLG
+      CHARACTER CBDPNM*12,TEXT12*12,TEXX12*12
+*----
+*  ALLOCATABLE ARRAYS
+*----
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: NPSYS
+      REAL, ALLOCATABLE, DIMENSION(:) :: AWPHI,FUN,SUN,SIGG
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: WEIGH,TOTPT,SIGWS,PAV
+      DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: MATRIX
+      TYPE(C_PTR) SIGP_PTR
+      REAL, POINTER, DIMENSION(:) :: SIGP
+*----
+*  STATEMENT FUNCTIONS
+*----
+      INM(IND,INOR,NBNRS)=(INOR-1)*NBNRS+IND
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      NORI=NOR(IRES,IGRP)
+      ALLOCATE(WEIGH(MAXNOR,NIRES),TOTPT(MAXNOR,NIRES),
+     1 SIGWS(MAXNOR,NIRES),PAV(0:NBNRS,0:NBNRS),AWPHI(0:NBNRS),
+     2 SIGG(0:NBMIX))
+      ALLOCATE(MATRIX(NBNRS*NORI,NBNRS*NORI+1))
+*----
+*  RECOVER THE SPECIFIC DIRECTORY FOR IRES-TH RESONANT ISOTOPE.
+*----
+      WRITE(CBDPNM,'(3HCOR,I4.4,1H/,I4.4)') IRES,NIRES
+      CALL LCMSIX(IPLI0,CBDPNM,1)
+      IPSYS=LCMGID(IPLI0,'ASSEMB-RIBON')
+      CALL LCMSIX(IPLI0,' ',2)
+*----
+*  COMPUTE THE AVERAGED COLLISION PROBABILITY MATRIX.
+*----
+      ALLOCATE(NPSYS(NORI*(NBNRS+1)))
+      ALLOCATE(FUN(NUN*NORI*(NBNRS+1)),SUN(NUN*NORI*(NBNRS+1)))
+      FUN(:NUN*NORI*(NBNRS+1))=0.0
+      SUN(:NUN*NORI*(NBNRS+1))=0.0
+      DO 50 INOR=1,NORI
+      DO 40 JNBN=0,NBNRS
+      NPSYS((INOR-1)*(NBNRS+1)+JNBN+1)=IASM+INOR
+      T1=0.0
+      DO 10 I=1,NREG
+      IBM=MAT(I)
+      IF(IBM.EQ.0) GO TO 10
+      IND=IREX(IBM)
+      IF((JNBN.EQ.0).AND.(IND.EQ.0)) THEN
+         T1=T1+SIGGAR(IBM,0,IGRP,3)*VOL(I)
+      ELSE IF(IND.EQ.JNBN) THEN
+         T1=T1+VOL(I)
+      ENDIF
+   10 CONTINUE
+      SIGG(0:NBMIX)=0.0
+      DO 20 IBM=1,NBMIX
+      IND=IREX(IBM)
+      IF((JNBN.EQ.0).AND.(IND.EQ.0)) THEN
+         SIGG(IBM)=SIGG(IBM)+SIGGAR(IBM,0,IGRP,3)
+      ELSE IF(IND.EQ.JNBN) THEN
+         SIGG(IBM)=SIGG(IBM)+1.0
+      ENDIF
+   20 ENDDO
+      IOF=(INOR-1)*NUN*(NBNRS+1)+JNBN*NUN
+      CALL DOORS(CDOOR,IPTRK,NBMIX,0,NUN,SIGG,SUN(IOF+1))
+      DO 30 I=1,NUN
+      IF(T1.NE.0.0) SUN(IOF+I)=SUN(IOF+I)/T1
+   30 CONTINUE
+   40 CONTINUE
+   50 CONTINUE
+*----
+*  SOLVE FOR THE MULTIBAND FLUX.
+*----
+      IDIR=0
+      NABS=NORI*(NBNRS+1)
+      LEXAC=.FALSE.
+      IPMACR=C_NULL_PTR
+      IPSOU=C_NULL_PTR
+      REBFLG=.FALSE.
+      CALL DOORFV(CDOOR,IPSYS,NPSYS,IPTRK,IFTRAK,IMPX,NABS,NBMIX,
+     1 IDIR,NREG,NUN,IPHASE,LEXAC,MAT,VOL,KEYFLX,TITR,SUN,FUN,IPMACR,
+     2 IPSOU,REBFLG)
+*----
+*  HOMOGENIZE THE MULTIBAND FLUX.
+*----
+      DO 170 INOR=1,NORI
+      PAV(0:NBNRS,0:NBNRS)=0.0
+      DO 155 JNBN=0,NBNRS
+      DO 150 I=1,NREG
+      IBM=MAT(I)
+      IF(IBM.EQ.0) GO TO 150
+      IOF=(INOR-1)*NUN*(NBNRS+1)+JNBN*NUN+KEYFLX(I)-1
+      PAV(IREX(IBM),JNBN)=PAV(IREX(IBM),JNBN)+FUN(IOF+1)*VOL(I)
+  150 CONTINUE
+  155 CONTINUE
+      DO 165 I=0,NBNRS
+      DO 160 J=0,NBNRS
+      IF(VOLMER(I).NE.0.0) PAV(I,J)=PAV(I,J)*VOLMER(J)/VOLMER(I)
+  160 CONTINUE
+  165 CONTINUE
+      KPSYS=LCMGIL(IPSYS,IASM+INOR)
+      CALL LCMPUT(KPSYS,'DRAGON-PAV',(NBNRS+1)**2,2,PAV(0,0))
+  170 CONTINUE
+      DEALLOCATE(SUN,FUN,NPSYS)
+*----
+*  COLLECT THE BASE POINTS IN TOTAL AND PARTIAL CROSS SECTION.
+*----
+      DO 200 JRES=1,NIRES
+      JPLIB=LCMGID(IPPT1(JRES),'GROUP-PT')
+      CALL LCMLEL(JPLIB,IGRP,ILONG,ITYLCM)
+      IF(ILONG.NE.0) THEN
+         KPLIB=LCMGIL(JPLIB,IGRP)
+         CALL LCMINF(KPLIB,TEXT12,TEXX12,EMPTY,ILONG,LCM)
+         CALL LCMLEN(KPLIB,'PROB-TABLE',LENG,ITYLCM)
+         NPART=LENG/MAXNOR
+         IF(LCM) THEN
+            CALL LCMGPD(KPLIB,'PROB-TABLE',SIGP_PTR)
+            CALL C_F_POINTER(SIGP_PTR,SIGP,(/ MAXNOR*NPART /))
+         ELSE
+            ALLOCATE(SIGP(MAXNOR*NPART))
+            CALL LCMGET(KPLIB,'PROB-TABLE',SIGP)
+         ENDIF
+         DO 190 INOR=1,NOR(JRES,IGRP)
+         WEIGH(INOR,JRES)=SIGP(INOR)
+         TOTPT(INOR,JRES)=SIGP(MAXNOR+INOR)
+         SIGWS(INOR,JRES)=SIGP(3*MAXNOR+INOR)
+  190    CONTINUE
+         IF(.NOT.LCM) DEALLOCATE(SIGP)
+      ELSE
+         WEIGH(1,JRES)=1.0
+         TOTPT(1,JRES)=STGAR(IPPT2(JRES,1),JRES,IGRP)
+         SIGWS(1,JRES)=SSGAR(IPPT2(JRES,1),JRES,1,IGRP)
+      ENDIF
+  200 CONTINUE
+*----
+*  RESPONSE MATRIX APPROACH. LOOP OVER THE SECONDARY SUBGROUPS.
+*----
+      DO 272 INOR=1,NORI
+      KPSYS=LCMGIL(IPSYS,IASM+INOR)
+      CALL LCMGET(KPSYS,'DRAGON-PAV',PAV(0,0))
+*----
+*  LOOP OVER THE PRIMARY SUBGROUPS. NORI+1 IS THE SOURCE.
+*----
+      DO 271 JNOR=1,NORI+1
+      IF(JNOR.LE.NORI) THEN
+         JNBMAX=NBNRS
+      ELSE
+         JNBMAX=1
+      ENDIF
+      DO 270 JNBN=1,JNBMAX
+      AWPHI(1:NBNRS)=0.0
+      DO 250 I=1,NREG
+      IBM=MAT(I)
+      IF(IBM.EQ.0) GO TO 250
+      JND=IREX(IBM)
+      QQQ=0.0
+      IF(JNOR.EQ.NORI+1) THEN
+         DO 230 JRES=0,NIRES
+         IF(JRES.EQ.0) THEN
+            QQQ=QQQ+SIGGAR(IBM,0,IGRP,3)
+         ELSE IF((JRES.NE.IRES).AND.(JND.GT.0)) THEN
+            QQQ=QQQ+SIGGAR(IBM,JRES,IGRP,4)
+         ENDIF
+  230    CONTINUE
+      ELSE IF((JND.EQ.JNBN).AND.(JNOR.NE.INOR)) THEN
+         QQQ=QQQ-WEIGH(JNOR,IRES)*CONR(JND,IRES)*SIGWS(JNOR,IRES)
+      ENDIF
+      DO 240 IND=1,NBNRS
+      AWPHI(IND)=AWPHI(IND)+PAV(IND,JND)*QQQ*VOL(I)/VOLMER(JND)
+  240 CONTINUE
+  250 CONTINUE
+      DO 260 IND=1,NBNRS
+      MATRIX(INM(IND,INOR,NBNRS),INM(JNBN,JNOR,NBNRS))=AWPHI(IND)
+  260 CONTINUE
+  270 CONTINUE
+  271 CONTINUE
+  272 CONTINUE
+*
+      DO 280 I=1,NBNRS*NORI
+      MATRIX(I,I)=MATRIX(I,I)+1.0D0
+  280 CONTINUE
+      CALL ALSBD(NBNRS*NORI,1,MATRIX,IER,NBNRS*NORI)
+      IF(IER.NE.0) CALL XABORT('USSEXC: SINGULAR MATRIX.')
+      XFLUX(:NBNRS,:NORI)=0.0
+      DO 295 IND=1,NBNRS
+      DO 290 INOR=1,NORI
+      I1=INM(IND,INOR,NBNRS)
+      XFLUX(IND,INOR)=REAL(MATRIX(I1,NBNRS*NORI+1))
+  290 CONTINUE
+  295 CONTINUE
+* END OF RESPONSE MATRIX APPROACH.
+*
+*----
+*  COLLECT THE BASE POINTS IN PARTIAL CROSS SECTION.
+*----
+      DO 340 JRES=1,NIRES
+      JPLIB=LCMGID(IPPT1(JRES),'GROUP-PT')
+      CALL LCMLEL(JPLIB,IGRP,ILONG,ITYLCM)
+      IF(ILONG.NE.0) THEN
+         KPLIB=LCMGIL(JPLIB,IGRP)
+         CALL LCMINF(KPLIB,TEXT12,TEXX12,EMPTY,ILONG,LCM)
+         CALL LCMLEN(KPLIB,'PROB-TABLE',LENG,ITYLCM)
+         NPART=LENG/MAXNOR
+         IF(LCM) THEN
+            CALL LCMGPD(KPLIB,'PROB-TABLE',SIGP_PTR)
+            CALL C_F_POINTER(SIGP_PTR,SIGP,(/ MAXNOR*NPART /))
+         ELSE
+            ALLOCATE(SIGP(MAXNOR*NPART))
+            CALL LCMGET(KPLIB,'PROB-TABLE',SIGP)
+         ENDIF
+         DO 330 INOR=1,NOR(JRES,IGRP)
+         WEIGH(INOR,JRES)=SIGP(INOR)
+         SIGWS(INOR,JRES)=SIGP(3*MAXNOR+INOR)
+  330    CONTINUE
+         IF(.NOT.LCM) DEALLOCATE(SIGP)
+      ELSE
+         WEIGH(1,JRES)=1.0
+         SIGWS(1,JRES)=SSGAR(IPPT2(JRES,1),JRES,1,IGRP)
+      ENDIF
+  340 CONTINUE
+*----
+*  COMPUTE THE AVERAGED SOURCE.
+*----
+      ALLOCATE(FUN(NUN*NORI),SUN(NUN*NORI))
+      SUN(:NUN*NORI)=0.0
+      ALLOCATE(NPSYS(NORI))
+      DO 385 INOR=1,NORI
+      NPSYS(INOR)=IASM+INOR
+      KPSYS=LCMGIL(IPSYS,IASM+INOR)
+      CALL LCMLEN(KPSYS,'FUNKNO$USS',ILENG,ITYLCM)
+      IF(ILENG.EQ.NUN) THEN
+        CALL LCMGET(KPSYS,'FUNKNO$USS',FUN((INOR-1)*NUN+1))
+      ELSE
+        FUN((INOR-1)*NUN+1:INOR*NUN)=0.0
+      ENDIF
+      SIGG(0)=0.0
+      DO 380 IBM=1,NBMIX
+      QQQ=SIGGAR(IBM,0,IGRP,3)
+      IND=IREX(IBM)
+      DO 360 JRES=1,NIRES
+      IF((JRES.NE.IRES).AND.(IND.GT.0)) THEN
+         QQQ=QQQ+SIGGAR(IBM,JRES,IGRP,4)
+      ENDIF
+  360 CONTINUE
+      IF(IND.GT.0) THEN
+         DO 370 JNOR=1,NORI
+         IF(JNOR.NE.INOR) THEN
+            QQQ=QQQ+WEIGH(JNOR,IRES)*CONR(IND,IRES)*SIGWS(JNOR,IRES)*
+     1      XFLUX(IND,JNOR)
+         ENDIF
+  370    CONTINUE
+      ENDIF
+      SIGG(IBM)=QQQ*WEIGH(INOR,IRES)
+  380 CONTINUE
+      IOF=(INOR-1)*NUN
+      CALL DOORS(CDOOR,IPTRK,NBMIX,0,NUN,SIGG,SUN(IOF+1))
+  385 CONTINUE
+*
+      IF(IMPX.GT.0) THEN
+         WRITE(TEXT12,'(3A4)') (IPPT2(IRES,I),I=2,4)
+         WRITE(6,'(15H USSEXC: GROUP=,I5,24H. SUBGROUP CALCULATION B,
+     1   37HASED ON RESPONSE MATRICES.  ISOTOPE='',A12,2H''.)') IGRP,
+     2   TEXT12
+      ENDIF
+*----
+*  SOLVE FOR THE MULTIBAND FLUX (VECTOR OF LENGTH NREG).
+*----
+      IPMACR=C_NULL_PTR
+      IPSOU=C_NULL_PTR
+      REBFLG=.FALSE.
+      CALL DOORFV(CDOOR,IPSYS,NPSYS,IPTRK,IFTRAK,IMPX,NORI,NBMIX,IDIR,
+     1 NREG,NUN,IPHASE,LEXAC,MAT,VOL,KEYFLX,TITR,SUN,FUN,IPMACR,IPSOU,
+     2 REBFLG)
+      DEALLOCATE(NPSYS)
+*----
+*  INTEGRATE THE REGION-ORDERED FLUX OVER SUBGROUPS AND COMPUTE UNGAR,
+*  THE REGION-ORDERED FLUX.
+*----
+      UNGAR(:NREG,IRES,IGRP)=0.0
+      DO 420 INOR=1,NORI
+      KPSYS=LCMGIL(IPSYS,IASM+INOR)
+      IOF=(INOR-1)*NUN
+      CALL LCMPUT(KPSYS,'FUNKNO$USS',NUN,2,FUN(IOF+1))
+*
+      DO 410 I=1,NREG
+      IOF=(INOR-1)*NUN+KEYFLX(I)-1
+      UNGAR(I,IRES,IGRP)=UNGAR(I,IRES,IGRP)+FUN(IOF+1)
+  410 CONTINUE
+  420 CONTINUE
+      DEALLOCATE(SUN,FUN)
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(MATRIX)
+      DEALLOCATE(SIGG,AWPHI,PAV,SIGWS,TOTPT,WEIGH)
+      RETURN
+      END