summaryrefslogtreecommitdiff
path: root/Dragon/src/MCGFST.f
blob: c1e25fe1ee351f4590889baa5cc41370509e1b35 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
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