Initial commit from Polytechnique Montreal

1 files changed, 138 insertions, 0 deletions

diff --git a/Dragon/src/MCGFFAL.f b/Dragon/src/MCGFFAL.f
new file mode 100644
index 0000000..bda8917
--- /dev/null
+++ b/Dragon/src/MCGFFAL.f
@@ -0,0 +1,138 @@
+*DECK MCGFFAL
+      SUBROUTINE MCGFFAL(SUBSCH,K,KPN,M,N,H,NOM,NZON,WEIGHT,XST,SP,SM,
+     1           DSP,DSM,NREG,NMU,NANI,NFUNLX,TRHAR,KEYCUR,IMU,B,F,
+     2           PHIV,DPHIV)
+*
+*-----------------------------------------------------------------------
+*
+*Purpose:
+* Solution of transport equation on a finite track.
+* Linear-discontinuous-characteristics approximation.
+* Ray-tracing (anisotropic scattering case,'source term isolation' off).
+*
+*Copyright:
+* Copyright (C) 2015 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
+* 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.
+* WEIGHT  track weight.
+* XST     macroscopic total cross section.
+* SP      total source vector for + direction.
+* SM      total source vector for - direction.
+* DSP     linear component of the total source vector for + direction.
+* DSM     linear component of the total source vector for - direction.
+* NREG    number of volumes.
+* NMU     order of the polar quadrature set.
+* NANI    scattering anisotropy (=1 for isotropic scattering).
+* NFUNLX  number of moments of the spherical harmonics.
+* NMOD    third dimension of TRHAR.
+* TRHAR   spherical harmonics components for this angle in the plane.
+* KEYCUR  position of current elements in PHI vector.
+* IMU     azimuthal angle corresponding to this track.
+*
+*Parameters: input/output
+* F       vector containing the zonal scalar flux (surface components).
+* PHIV    vector containing the zonal scalar flux (component 1).
+* DPHIV   vector containing the zonal scalar flux (components 2 and 3).
+*
+*Parameters: scratch
+* B       undefined.
+*
+*-----------------------------------------------------------------------
+*
+      IMPLICIT NONE
+*----
+*  SUBROUTINE ARGUMENTS
+*----
+      INTEGER K,KPN,M,N,NOM(N),NZON(K),NMU,NFUNLX,NREG,
+     1 KEYCUR(K-NREG),IMU,NANI
+      REAL XST(0:M),TRHAR(NMU,NFUNLX,2)
+      DOUBLE PRECISION WEIGHT,H(N),SP(N),SM(N),DSP(N),DSM(N),B(0:5,N),
+     1 F(KPN),PHIV(NFUNLX,NREG),DPHIV(2*NFUNLX,NREG)
+      EXTERNAL SUBSCH
+*----
+* LOCAL VARIABLES
+*---
+      REAL ETA,XI
+      DOUBLE PRECISION F0,DF0,SI,SJ,DSI,DSJ,RM,RP,DSIG,DH
+      INTEGER I,NOMI,IND1,INDN,JF,J,NOMJ
+*----
+*     Calculation of the coefficients for this track.
+*----
+*       MCGSCAL: Linear discontinuous-Characteristics Scheme with
+*                Tabulated Exponentials
+*       MCGDDFL: Diamond-Differencing DD1 Scheme
+*       MCGSCEL: Linear discontinuous-Characteristics Scheme with
+*                Exact Exponentials
+      IF(NANI.LE.0) CALL XABORT('MCGFFAL: INVALID VALUE OF NANI.')
+      CALL SUBSCH(N,K,M,NOM,NZON,H,XST,B)
+*----
+*     Summation along the track in both directions
+*----
+*     incoming flux in + direction
+      IND1=KEYCUR(-NOM(1))
+      RP=SP(1)
+*     incoming flux in - direction
+      INDN=KEYCUR(-NOM(N))
+      RM=SM(N)
+*     track angles in 3D
+      ETA=TRHAR(IMU,3,1)
+      XI=TRHAR(IMU,2,1)
+      DO I=2,N-1
+*        + direction
+         NOMI=NOM(I)
+         DSIG=XST(NZON(NOMI))
+         DH=H(I)
+         SI=SP(I)
+         DSI=DSP(I)
+         F0=B(1,I)*RP+B(2,I)*SI+B(3,I)*DSI
+         DF0=B(4,I)*RP-B(3,I)*SI/(DH*DH)+B(5,I)*DSIG*DSI
+         RP=B(0,I)*RP+B(1,I)*SI-B(3,I)*DSIG*DSI
+         DO JF=1,NFUNLX
+            PHIV(JF,NOMI)=PHIV(JF,NOMI)+WEIGHT*F0*TRHAR(IMU,JF,1)
+            DPHIV(JF,NOMI)=DPHIV(JF,NOMI)+WEIGHT*DF0*ETA*
+     >                     TRHAR(IMU,JF,1)
+            DPHIV(NFUNLX+JF,NOMI)=DPHIV(NFUNLX+JF,NOMI)+WEIGHT*DF0*XI*
+     >                            TRHAR(IMU,JF,1)
+         ENDDO
+*        - direction
+         J=N+1-I
+         NOMJ=NOM(J)
+         DSIG=XST(NZON(NOMJ))
+         DH=H(J)
+         SJ=SM(J)
+         DSJ=DSM(J)
+         F0=B(1,J)*RM+B(2,J)*SJ+B(3,J)*DSJ
+         DF0=B(4,J)*RM-B(3,J)*SJ/(DH*DH)+B(5,J)*DSIG*DSJ
+         RM=B(0,J)*RM+B(1,J)*SJ-B(3,J)*DSIG*DSJ
+         DO JF=1,NFUNLX
+            PHIV(JF,NOMJ)=PHIV(JF,NOMJ)+WEIGHT*F0*TRHAR(IMU,JF,2)
+            DPHIV(JF,NOMJ)=DPHIV(JF,NOMJ)-WEIGHT*DF0*ETA*
+     >                     TRHAR(IMU,JF,2)
+            DPHIV(NFUNLX+JF,NOMJ)=DPHIV(NFUNLX+JF,NOMJ)-WEIGHT*DF0*XI*
+     >                            TRHAR(IMU,JF,2)
+         ENDDO
+      ENDDO
+*     outgoing flux in + direction
+      F(INDN)=F(INDN)+WEIGHT*RP
+*     outgoing flux in - direction
+      F(IND1)=F(IND1)+WEIGHT*RM
+*
+      RETURN
+      END