summaryrefslogtreecommitdiff
path: root/Trivac/src/TRIRCA.f
blob: 8202f297f2e2c6d4a93c92b7d6f4a268d0c86188 (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
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
*DECK TRIRCA
      SUBROUTINE TRIRCA(IPMACR,IPMACP,NGRP,NBMIX,NANI,NW,LDIFF,IL,IPR,
     1 RCAT)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute the RCAT removal matrix in SPN cases.
*
*Copyright:
* Copyright (C) 2012 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
* IPMACR  L_MACROLIB pointer to the unperturbed cross sections.
* IPMACP  L_MACROLIB pointer to the perturbed cross sections if
*         IPR.gt.0. Equal to IPMACR if IPR=0.
* NGRP    number of energy groups.
* NBMIX   total number of material mixtures in the macrolib.
* NANI    maximum scattering order recovered from tracking and macrolib.
* NW      maximum Legendre order (0 or 1) for the total cross sections.
* LDIFF   flag set to .true. to use 1/3D as 'NTOT1' cross sections.
* IL      scattering Legendre order.
* IPR     type of assembly:
*         =0: calculation of the system matrices;
*         =1: calculation of the derivative of these matrices;
*         =2: calculation of the first variation of these matrices;
*         =3: identical to IPR=2, but these variation are added to
*         unperturbed system matrices.
*
*Parameters: output
* RCAT    removal matrix.
*
*-----------------------------------------------------------------------
*
      USE GANLIB
*----
*  SUBROUTINE ARGUMENTS
*----
      TYPE(C_PTR) IPMACR,IPMACP
      INTEGER NGRP,NBMIX,NANI,NW,IL,IPR
      LOGICAL LDIFF
      DOUBLE PRECISION RCAT(NGRP,NGRP,NBMIX)
*----
*  LOCAL VARIABLES
*----
      TYPE(C_PTR) JPMACR,JPMACP,KPMACR,KPMACP
      CHARACTER TEXT12*12,CM*2,HSMG*131
      DOUBLE PRECISION OTH
      INTEGER, ALLOCATABLE, DIMENSION(:) :: IJJ,NJJ,IPOS
      REAL, ALLOCATABLE, DIMENSION(:) :: WORK
      REAL, ALLOCATABLE, DIMENSION(:,:) :: SGD
      PARAMETER(OTH=1.0D0/3.0D0)
*----
*  SCRATCH STORAGE ALLOCATION
*----
      ALLOCATE(IJJ(NBMIX),NJJ(NBMIX),IPOS(NBMIX))
      ALLOCATE(SGD(NBMIX,3),WORK(NBMIX*NGRP))
*
      JPMACR=LCMGID(IPMACR,'GROUP')
      JPMACP=LCMGID(IPMACP,'GROUP')
      WRITE(CM,'(I2.2)') IL-1
      RCAT(:NGRP,:NGRP,:NBMIX)=0.0D0
      DO 100 IGR=1,NGRP
*     PROCESS SECONDARY GROUP IGR.
      KPMACP=LCMGIL(JPMACP,IGR)
      SGD(:NBMIX,1)=0.0
      CALL LCMLEN(KPMACP,'SIGW'//CM,LENGT,ITYLCM)
      IF((LENGT.GT.0).AND.(IL.LE.NANI)) THEN
         IF(LENGT.GT.NBMIX) CALL XABORT('TRIRCA: INVALID LENGTH FOR'
     1   //' SIGW'//CM//' CROSS SECTIONS.')
         CALL LCMGET(KPMACP,'SIGW'//CM,SGD(1,1))
      ENDIF
      WRITE(TEXT12,'(4HNTOT,I1)') MIN(IL-1,9)
      CALL LCMLEN(KPMACP,TEXT12,LENGT,ITYLCM)
      CALL LCMLEN(KPMACP,'NTOT1',LENGT1,ITYLCM)
      IF((IL.EQ.1).AND.(LENGT.NE.NBMIX)) CALL XABORT('TRIRCA: NO NTOT0'
     1 //' CROSS SECTIONS.')
      IF(MOD(IL-1,2).EQ.0) THEN
*        macroscopic total cross section in even-parity equations.
         IF(LENGT.EQ.NBMIX) THEN
            CALL LCMGET(KPMACP,TEXT12,SGD(1,2))
         ELSE
            CALL LCMGET(KPMACP,'NTOT0',SGD(1,2))
         ENDIF
         DO 10 IBM=1,NBMIX
         IF((SGD(IBM,2)-SGD(IBM,1).LT.0.0).AND.(IPR.EQ.0)) THEN
            WRITE(HSMG,'(28HTRIRCA: NEGATIVE XS IN GROUP,I5)') IGR
            CALL XABORT(HSMG)
         ENDIF
         RCAT(IGR,IGR,IBM)=SGD(IBM,2)-SGD(IBM,1)
   10    CONTINUE
      ELSE
*        macroscopic total cross section in odd-parity equations.
         IF(LDIFF) THEN
            CALL LCMLEN(KPMACP,'DIFF',LENGT,ITYLCM)
            IF(LENGT.EQ.0) CALL XABORT('TRIRCA: DIFFUSION COEFFICIENTS'
     1      //' EXPECTED IN THE MACROLIB.')
            IF(LENGT.GT.NBMIX) CALL XABORT('TRIRCA: INVALID LENGTH FOR'
     1      //' DIFFUSION COEFFICIENTS.')
            CALL LCMGET(KPMACP,'DIFF',SGD(1,2))
            IF(IPR.EQ.0) THEN
              DO 20 IBM=1,NBMIX
              RCAT(IGR,IGR,IBM)=OTH/SGD(IBM,2)
   20         CONTINUE
            ELSE IF(IPR.EQ.1) THEN
              KPMACR=LCMGIL(JPMACR,IGR)
              CALL LCMGET(KPMACR,'DIFF',SGD(1,3))
              DO 30 IBM=1,NBMIX
              RCAT(IGR,IGR,IBM)=-OTH*SGD(IBM,2)/SGD(IBM,3)**2
   30         CONTINUE
            ELSE IF(IPR.EQ.2) THEN
              KPMACR=LCMGIL(JPMACR,IGR)
              CALL LCMGET(KPMACR,'DIFF',SGD(1,3))
              DO 40 IBM=1,NBMIX
              RCAT(IGR,IGR,IBM)=OTH/(SGD(IBM,2)+SGD(IBM,3))
     1                         -OTH/SGD(IBM,3)
   40         CONTINUE
            ELSE IF(IPR.EQ.3) THEN
              KPMACR=LCMGIL(JPMACR,IGR)
              CALL LCMGET(KPMACR,'DIFF',SGD(1,3))
              DO 50 IBM=1,NBMIX
              RCAT(IGR,IGR,IBM)=OTH/(SGD(IBM,2)+SGD(IBM,3))
   50         CONTINUE
            ENDIF
            GO TO 100
         ELSE
            IF(NW.EQ.0) THEN
               CALL LCMGET(KPMACP,'NTOT0',SGD(1,2))
            ELSE IF(LENGT.EQ.NBMIX) THEN
               CALL LCMGET(KPMACP,TEXT12,SGD(1,2))
            ELSE IF((NW.GE.1).AND.(LENGT1.EQ.NBMIX)) THEN
               CALL LCMGET(KPMACP,'NTOT1',SGD(1,2))
            ELSE
               CALL LCMGET(KPMACP,'NTOT0',SGD(1,2))
            ENDIF
            DO 60 IBM=1,NBMIX
            RCAT(IGR,IGR,IBM)=SGD(IBM,2)-SGD(IBM,1)
   60       CONTINUE
         ENDIF
         IF(IPR.EQ.0) THEN
            DO 65 IBM=1,NBMIX
            IF(RCAT(IGR,IGR,IBM).LT.0.0) THEN
               WRITE(HSMG,'(39HTRIRCA: INVALID CROSS-SECTION DATA (IL=,
     1         I3,2H).)') IL
               CALL XABORT(HSMG)
            ENDIF
   65       CONTINUE
         ENDIF
      ENDIF
      CALL LCMLEN(KPMACP,'NJJS'//CM,LENGT,ITYLCM)
      IF(LENGT.GT.NBMIX) CALL XABORT('TRIRCA: INVALID LENGTH FOR NJJS'
     1 //CM//' INFORMATION.')
      IF((LENGT.GT.0).AND.(IL.LE.NANI)) THEN
         CALL LCMGET(KPMACP,'NJJS'//CM,NJJ)
         CALL LCMGET(KPMACP,'IJJS'//CM,IJJ)
         IGMIN=IGR
         IGMAX=IGR
         DO 70 IBM=1,NBMIX
         IGMIN=MIN(IGMIN,IJJ(IBM)-NJJ(IBM)+1)
         IGMAX=MAX(IGMAX,IJJ(IBM))
   70    CONTINUE
         CALL LCMGET(KPMACP,'IPOS'//CM,IPOS)
         CALL LCMGET(KPMACP,'SCAT'//CM,WORK)
         DO 90 JGR=IGMAX,IGMIN,-1
         IF(JGR.EQ.IGR) GO TO 90
         DO 80 IBM=1,NBMIX
         IF((JGR.GT.IJJ(IBM)-NJJ(IBM)).AND.(JGR.LE.IJJ(IBM))) THEN
            RCAT(IGR,JGR,IBM)=-WORK(IPOS(IBM)+IJJ(IBM)-JGR)
         ENDIF
   80    CONTINUE
   90    CONTINUE
      ENDIF
  100 CONTINUE
*----
*  SCRATCH STORAGE DEALLOCATION
*----
      DEALLOCATE(WORK,SGD)
      DEALLOCATE(IPOS,NJJ,IJJ)
      RETURN
      END