Initial commit from Polytechnique Montreal

1 files changed, 136 insertions, 0 deletions

diff --git a/Dragon/src/MCGFFIS.f b/Dragon/src/MCGFFIS.f
new file mode 100644
index 0000000..9f68a5e
--- /dev/null
+++ b/Dragon/src/MCGFFIS.f
@@ -0,0 +1,136 @@
+*DECK MCGFFIS
+      SUBROUTINE MCGFFIS(SUBSCH,K,KPN,M,N,H,NOM,NZON,XST,S,NREG,KEYFLX,
+     1           KEYCUR,F,B,W,OMEGA2,IDIR,NSOUT,XSI)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Solution of transport equation on a track
+* ray-tracing (isotropic scattering case,'source term isolation' off).
+*
+*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
+* SUBSCH  Track coefficients calculation subroutine.
+* K       total number of volumes for which specific values
+*         of the neutron flux and reactions rates are required.
+* KPN     total number of unknowns in vectors F.
+* M       number of material mixtures.
+* N       number of elements in the current track.
+* H       vector containing the lenght of the different segments of this
+*         track.
+* NOM     vector containing the region number of the different segments
+*         of this track.
+* NZON    index-number of the mixture type assigned to each volume.
+* XST     macroscopic total cross section.
+* S       total source vector.
+* NREG    number of volumes.
+* KEYFLX  position of flux elements in PHI vector.
+* KEYCUR  position of current elements in PHI vector.
+* W       weight associated with this track.
+* OMEGA2  square x, y and z-component of the direction
+*         Omega for 2D geometry.
+* IDIR    direction of fundamental current for TIBERE with MoC
+*         (=0,1,2,3).
+* NSOUT   number of outer surfaces.
+* XSI     x,y and z component of the shape parameter for TIBERE. 
+*
+*Parameters: input/output
+* F       vector containing the zonal scalar flux.
+*
+*Parameters: scratch
+* B       undefined.
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER K,KPN,M,N,NOM(N),NZON(K),NREG,KEYFLX(NREG,1),
+     1 KEYCUR(K-NREG),IDIR
+      REAL XST(0:M)
+      DOUBLE PRECISION W,H(N),S(KPN),F(KPN),B(2,N),OMEGA2(3)
+      INTEGER NSOUT
+      DOUBLE PRECISION XSI(NSOUT)
+      EXTERNAL SUBSCH
+*---
+* LOCAL VARIABLES
+*---
+      DOUBLE PRECISION F0,SI,SJ,RM,RP,WW
+      INTEGER I,J,NOMI,IND,IND1,INDN,NOMJ,INDC
+*
+      WW=DBLE(W)
+*----
+*     Calculation coefficients for this track.
+*----
+*          MCGSCAS: Step-Characteristics Scheme with Tabulated Exponentials
+*          MCGDDFS: Diamond-Differencing Scheme
+*          MCGSCES: Step-Characteristics Scheme with Exact Exponentials
+      CALL SUBSCH(N,K,M,NOM,NZON,H,XST,B)
+*----
+*     Summation along the track in both directions
+*----
+*     Calculation of the incoming current on surface with the correction XSI
+*     incoming flux in + direction
+      IND1=KEYCUR(-NOM(1))
+      RP=S(IND1)
+*     incoming flux in - direction
+      INDN=KEYCUR(-NOM(N))
+      RM=S(INDN)
+      IF(IDIR.GT.0) THEN
+*     incoming flux in + direction
+        RP=RP*OMEGA2(IDIR)/XSI(IND1-NREG)
+*     incoming flux in - direction
+        RM=RM*OMEGA2(IDIR)/XSI(INDN-NREG)
+      ENDIF
+*                        
+      DO I=2,N-1
+         NOMI=NOM(I)
+         SI=S(NOMI)
+         J=N+1-I
+         NOMJ=NOM(J)
+         SJ=S(NOMJ)
+         IF(IDIR.EQ.0) THEN
+*        + direction
+           F0=B(1,I)*RP+B(2,I)*SI
+           RP=RP+B(1,I)*(SI-XST(NZON(NOMI))*RP)
+           IND=KEYFLX(NOMI,1)
+           F(IND)=F(IND)+F0*WW
+*        - direction
+           F0=B(1,J)*RM+B(2,J)*SJ
+           RM=RM+B(1,J)*(SJ-XST(NZON(NOMJ))*RM)
+           IND=KEYFLX(NOMJ,1)
+           F(IND)=F(IND)+F0*WW
+         ELSE
+*        Solve current equation for TIBERE
+*        and calculate Xi, Yi and Zi 
+*        + direction
+           F0=B(1,I)*RP+B(2,I)*SI*OMEGA2(IDIR)
+           RP=RP+B(1,I)*(SI*OMEGA2(IDIR)-XST(NZON(NOMI))*RP)
+           IND=KEYFLX(NOMI,1)
+           F(IND)=F(IND)+F0*WW
+           INDC=KPN/2+IND
+           F(INDC)=F(INDC)+F0*WW*OMEGA2(IDIR)
+*        - direction
+           F0=B(1,J)*RM+B(2,J)*SJ*OMEGA2(IDIR)
+           RM=RM+B(1,J)*(SJ*OMEGA2(IDIR)-XST(NZON(NOMJ))*RM)
+           IND=KEYFLX(NOMJ,1)
+           F(IND)=F(IND)+F0*WW
+           INDC=KPN/2+IND
+           F(INDC)=F(INDC)+F0*WW*OMEGA2(IDIR)
+         ENDIF
+      ENDDO
+*     outgoing flux in + direction
+       F(INDN)=F(INDN)+RP*WW
+*     outgoing flux in - direction
+       F(IND1)=F(IND1)+RM*WW    
+      RETURN
+      END