From 7dfcc480ba1e19bd3232349fc733caef94034292 Mon Sep 17 00:00:00 2001
From: stainer_t <thomas.stainer@oecd-nea.org>
Date: Mon, 8 Sep 2025 13:48:49 +0200
Subject: Initial commit from Polytechnique Montreal

---
 Dragon/src/MCGFST.f | 135 ++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 135 insertions(+)
 create mode 100644 Dragon/src/MCGFST.f

(limited to 'Dragon/src/MCGFST.f')

diff --git a/Dragon/src/MCGFST.f b/Dragon/src/MCGFST.f
new file mode 100644
index 0000000..c1e25fe
--- /dev/null
+++ b/Dragon/src/MCGFST.f
@@ -0,0 +1,135 @@
+*DECK MCGFST
+      SUBROUTINE MCGFST(NGEFF,KPSYS,NCONV,KPN,K,NREG,NANI,NFUNL,NPJJM,
+     1           KEYFLX,KEYCUR,PJJIND,NZON,V,S,PHIOUT,IDIR)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Addition of the contribution to the flux of the regional source
+* when the 'source term isolation' option is turned on for the 
+* method of characteristics integration.
+*
+*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): R. Le Tellier
+*
+*Parameters: input
+* NGEFF   number of groups to process.
+* KPSYS   pointer array for each group properties.
+* NCONV   logical array of convergence status for each group (.TRUE.
+*         not converged).
+* KPN     total number of unknowns per group in flux vector.
+* K       total number of volumes for which specific values
+*         of the neutron flux and reactions rates are required.
+* NREG    number of volumes.
+* NANI    scattering anisotropy (=1 for isotropic scattering).
+* NFUNL   number of moments of the flux (in 2D : NFUNL=NANI*(NANI+1)/2).
+* NPJJM   number of pjj modes to store for STIS option.
+* KEYFLX  position of flux elements in flux vector.
+* KEYCUR  position of current elements in flux vector.
+* PJJIND  index for pjj(nu <- nu') modes.
+* NZON    index-number of the mixture type assigned to each volume.
+* V       volumes.
+* S       source vector.
+* IDIR    direction of fundamental current for TIBERE with MoC
+*         (=0,1,2,3).
+*
+*Parameters: input/output
+* PHIOUT  flux vector.
+*
+*-----------------------------------------------------------------------
+*
+      USE GANLIB
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      TYPE(C_PTR) KPSYS(NGEFF)
+      INTEGER NGEFF,KPN,K,NREG,NANI,NFUNL,NPJJM,KEYFLX(NREG,NFUNL),
+     1 KEYCUR(K-NREG),PJJIND(NPJJM,2),NZON(K),IDIR
+      REAL V(K)
+      DOUBLE PRECISION S(KPN,NGEFF),PHIOUT(KPN,NGEFF)
+      LOGICAL NCONV(NGEFF)
+*----
+*  LOCAL VARIABLES
+*----
+      TYPE(C_PTR) JPSYS
+      INTEGER II,I,IND,IMOD,INU,INUP,INDP
+      REAL, ALLOCATABLE, DIMENSION(:,:) :: PJJ,PJJI
+      CHARACTER*12 NPJJT,NPJJIT
+*----
+*  SCRATCH STORAGE ALLOCATION
+*----
+      ALLOCATE(PJJ(NREG,NPJJM),PJJI(NREG,NPJJM))
+*
+      IF(IDIR .EQ.0) THEN
+        NPJJT='PJJ$MCCG'
+        NPJJIT='            '
+      ELSEIF(IDIR .EQ. 1) THEN
+        NPJJT='PJJX$MCCG'
+        NPJJIT='PJJXI$MCCG'
+      ELSEIF(IDIR .EQ. 2) THEN
+        NPJJT='PJJY$MCCG'
+        NPJJIT='PJJYI$MCCG'
+      ELSE
+        NPJJT='PJJZ$MCCG'
+        NPJJIT='PJJZI$MCCG'
+      ENDIF
+      IF(NANI.LE.0) CALL XABORT('MCGFST: INVALID VALUE OF NANI.')
+      DO II=1,NGEFF
+      IF(NCONV(II)) THEN
+         JPSYS=KPSYS(II)
+         CALL LCMGET(JPSYS,NPJJT,PJJ)
+         IF(IDIR.GT.0) CALL LCMGET(JPSYS,NPJJIT,PJJI)
+         DO I=1,K
+         IF(V(I).GT.0.) THEN 
+            IF(NZON(I).LT.0) THEN
+               IND=KEYCUR(I-NREG)
+               PHIOUT(IND,II)=PHIOUT(IND,II)/V(I)
+            ELSE
+               DO INU=1,NFUNL
+                  IND=KEYFLX(I,INU)
+                  PHIOUT(IND,II)=PHIOUT(IND,II)/V(I)
+               ENDDO
+*              DIVIDE THE EXTRA TERMS XI, YI, AND ZI BY THE VOLUME
+               IF(IDIR.NE.0) THEN
+                 IND=KEYFLX(I,NFUNL)
+                 PHIOUT(IND+KPN/2,II)=PHIOUT(IND+KPN/2,II)/V(I)
+               ENDIF
+               DO IMOD=1,NPJJM
+                  INU=PJJIND(IMOD,1)
+                  INUP=PJJIND(IMOD,2)
+                  IND=KEYFLX(I,INU)
+                  INDP=KEYFLX(I,INUP)
+                  PHIOUT(IND,II)=PHIOUT(IND,II)+
+     1                           PJJ(I,IMOD)*S(INDP,II)
+                  IF(IDIR .GT. 0) THEN
+                     PHIOUT(IND+KPN/2,II)=PHIOUT(IND+KPN/2,II)+
+     1                           PJJI(I,IMOD)*S(INDP,II)
+                  ENDIF
+                  IF(INU.NE.INUP) THEN
+                     PHIOUT(INDP,II)=PHIOUT(INDP,II)+
+     1                           PJJ(I,IMOD)*S(IND,II)
+                     IF(IDIR .GT. 0) THEN
+                       PHIOUT(INDP+KPN/2,II)=PHIOUT(INDP+KPN/2,II)+
+     1                              PJJI(I,IMOD)*S(IND,II)
+                     ENDIF
+                  ENDIF
+               ENDDO
+            ENDIF
+         ENDIF
+         ENDDO
+      ENDIF
+      ENDDO
+*----
+*  SCRATCH STORAGE DEALLOCATION
+*----
+      DEALLOCATE(PJJI,PJJ)
+*
+      RETURN
+      END
-- 
cgit v1.2.3