1 files changed, 112 insertions, 0 deletions
diff --git a/Dragon/src/MCGFCS.f b/Dragon/src/MCGFCS.f
new file mode 100644
index 0000000..759ff03
--- /dev/null
+++ b/Dragon/src/MCGFCS.f
@@ -0,0 +1,112 @@
+*DECK MCGFCS
+      SUBROUTINE MCGFCS(N,NDIM,NZON,QN,FI,M,NANI,NLIN,NFUNL,SC,S,KPN,
+     1                  NREG,IPRINT,KEYFLX,KEYCUR,IBC,SIGAL,STIS)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Calculation of source for collision at iteration iter.
+*
+*Copyright:
+* Copyright (C) 2002 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): I. Suslov and R. Le Tellier
+*
+*Parameters: input
+* N       number of spatial unknowns.
+* NDIM    number of dimensions for the geometry.
+* NZON    index-number of the mixture type assigned to each volume.
+* QN      input source (fission-other groups) vector.
+* FI      unknown vector.
+* M       number of material mixtures.
+* NANI    scattering anisotropy (=1 for isotropic scattering).
+* NLIN    linear discontinuous flag (=1 SC/DD0; =3 LDC/DD1).
+* NFUNL   number of spherical harmonics components.
+* SC      macroscopic scattering cross section.
+* KPN     total number of unknowns in vectors QN and FI.
+* NREG    number of volumes.
+* IPRINT  print parameter (equal to zero for no print).
+* KEYFLX  position of flux elements in FI vector.
+* KEYCUR  position of current elements in FI.
+* IBC     index for boundary condition to connect a surface to another.
+* STIS    integration strategy flag.
+* SIGAL   total cross-section and albedo array.
+*
+*Parameters: output
+* S       source elements vector.
+*
+*-----------------------------------------------------------------------
+*
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER N,NDIM,NZON(N),M,NANI,NLIN,NFUNL,KPN,NREG,IPRINT,
+     1 KEYFLX(NREG,NLIN,NFUNL),KEYCUR(N-NREG),IBC(N-NREG),STIS
+      REAL QN(KPN),FI(KPN),SC(0:M,NANI),SIGAL(-6:M)
+      DOUBLE PRECISION S(KPN)
+*
+      IF(NDIM.EQ.2) THEN
+*     2D geometry
+      DO IR=1,N
+         IBM=NZON(IR)
+         IF(IBM.LT.0) THEN
+*        Boundary condition
+            ISUR=IR-NREG
+            ISUR2=IBC(ISUR)
+            IND=KEYCUR(ISUR)
+            IND2=KEYCUR(ISUR2)
+            IF(IND.GT.0) S(IND)=SIGAL(IBM)*FI(IND2)
+         ELSEIF(IBM.GE.0) THEN
+*        Volume cell
+            DO IL=0,NANI-1
+               XSC=REAL(2*IL+1)*SC(IBM,IL+1)
+               DO IM=0,IL
+                  DO IE=1,NLIN
+                     IND=KEYFLX(IR,IE,1+IL*(IL+1)/2+IM)
+                     IF(IND.GT.0) THEN
+                       S(IND)=QN(IND)+XSC*FI(IND)
+                       IF(STIS.EQ.-1) S(IND)=S(IND)/SIGAL(IBM)
+                     ENDIF
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+      ENDDO
+      ELSE ! NDIM.EQ.3   
+*     3D geometry     
+      DO IR=1,N
+         IBM=NZON(IR)
+         IF(IBM.LT.0) THEN
+*        Boundary condition
+            ISUR=IR-NREG
+            ISUR2=IBC(ISUR)
+            IND=KEYCUR(ISUR)
+            IND2=KEYCUR(ISUR2)
+            IF(IND.GT.0) S(IND)=SIGAL(IBM)*FI(IND2)
+         ELSEIF(IBM.GE.0) THEN
+*        Volume cell
+            INDA=0
+            DO IL=0,NANI-1
+               XSC=REAL(2*IL+1)*SC(IBM,IL+1)
+               DO IM=-IL,IL
+                  DO IE=1,NLIN
+                     INDA=INDA+1
+                     IND=KEYFLX(IR,IE,INDA)
+                     IF(IND.GT.0) THEN
+                       S(IND)=QN(IND)+XSC*FI(IND)
+                       IF(STIS.EQ.-1) S(IND)=S(IND)/SIGAL(IBM)
+                     ENDIF
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+      ENDDO
+      ENDIF
+*
+      IF(IPRINT.GT.6) CALL PRINDM ('S     ',S,KPN)
+      RETURN
+      END