summaryrefslogtreecommitdiff
path: root/Dragon/src/SYB001.f
blob: dfd99039abda9b76bcba5cfda4e3cd34b236586f (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
*DECK SYB001
      SUBROUTINE SYB001 (NREG,NSUPCE,NPIJ,SIGT,SIGW,IMPX,IQUAD,NMC,
     1 RAYRE,PIJW,PISW,PSJW,PSSW)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute the cellwise scattering-reduced collision, escape and
* transmission probabilities for the 'do-it-yourself' approach.
*
*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): A. Hebert
*
*Parameters: input
* NREG    total number of regions (NREG=NMC(NSUPCE+1)).
* NSUPCE  total number of cells.
* NPIJ    length of cellwise scattering-reduced collision probability
*         matrices.
* SIGT    total macroscopic cross sections.
* SIGW    P0 within-group scattering macroscopic cross sections.
* IMPX    print flag (equal to 0 for no print).
* IQUAD   quadrature parameter.
* NMC     offset of the first volume in each cell.
* RAYRE   radius of the tubes in each cell.
*
*Parameters: output
* PIJW    cellwise scattering-reduced collision probability matrices.
* PISW    cellwise scattering-reduced escape probability matrices.
* PSJW    cellwise scattering-reduced collision probability matrices
*         for incoming neutrons.
* PSSW    cellwise scattering-reduced transmission probability
*         matrices.
*
*-----------------------------------------------------------------------
*
*----
*  SUBROUTINE ARGUMENTS
*----
      INTEGER NSUPCE,NREG,NPIJ,IMPX,IQUAD,NMC(NSUPCE+1)
      REAL SIGT(NREG),SIGW(NREG),RAYRE(NREG+NSUPCE),PIJW(NPIJ),
     1 PISW(NREG),PSJW(NREG),PSSW(NSUPCE)
*----
*  LOCAL VARIABLES
*----
      PARAMETER (PI=3.141592654)
      LOGICAL LSKIP
      REAL, ALLOCATABLE, DIMENSION(:) :: PIS,PSJ,ZTR,WORK
      REAL, ALLOCATABLE, DIMENSION(:,:) :: PP
*----
*  SCRATCH STORAGE ALLOCATION
*----
      ALLOCATE(PIS(NREG),PSJ(NREG))
*
      IPIJ=0
      DO 160 IKK=1,NSUPCE
      J1=NMC(IKK)
      J2=NMC(IKK+1)-J1
*----
*  COMPUTE THE REDUCED COLLISION PROBABILITY MATRIX
*----
      ALLOCATE(PP(J2,J2),ZTR(1+IQUAD*((J2*(5+J2))/2)))
      CALL SYBT1D(J2,RAYRE(J1+IKK),.FALSE.,IQUAD,ZTR)
      CALL SYBALC(J2,J2,RAYRE(J1+IKK),SIGT(J1+1),IQUAD,0.0,ZTR,PP)
      DEALLOCATE(ZTR)
      SURFA=2.0*PI*RAYRE(J1+J2+IKK)
      PSS=0.0
      RJN=0.0
      DO 20 I=1,J2
      PIS(I)=0.0
      DO 10 J=1,J2
      PIS(I)=PIS(I)+PP(I,J)*SIGT(J+J1)
   10 CONTINUE
      PIS(I)=1.0-PIS(I)
      RJN1=RAYRE(I+J1+IKK)**2
      PSJ(I)=4.0*PI*(RJN1-RJN)*PIS(I)/SURFA
      PSS=PSS+PSJ(I)*SIGT(I+J1)
      RJN=RJN1
   20 CONTINUE
      PSS=1.0-PSS
      IF(IMPX.GE.8) THEN 
         CALL SYBPRX(1,1,J2,IKK,SIGT(J1+1),SIGW(J1+1),PP(1,1),PIS(1),
     1   PSJ(1),PSS)
      ENDIF
*----
*  CHECK IF SCATTERING REDUCTION IS REQUIRED
*----
      LSKIP=.TRUE.
      DO 30 I=1,J2
      LSKIP=LSKIP.AND.(SIGW(J1+I).EQ.0.0)
   30 CONTINUE
*----
*  SCATTERING REDUCTION IF LSKIP=.FALSE.
*----
      IF(LSKIP) THEN
*        DO NOT PERFORM SCATTERING REDUCTION.
         DO 45 I=1,J2
         DO 40 J=1,J2
         PIJW(IPIJ+(J-1)*J2+I)=PP(I,J)
   40    CONTINUE
   45    CONTINUE
         DO 50 I=1,J2
         PISW(J1+I)=PIS(I)
         PSJW(J1+I)=PSJ(I)
   50    CONTINUE
         PSSW(IKK)=PSS
      ELSE
*        COMPUTE THE SCATTERING-REDUCED COLLISION AND ESCAPE MATRICES.
         DO 70 I=1,J2
         DO 60 J=1,J2
         PIJW(IPIJ+(J-1)*J2+I)=-PP(I,J)*SIGW(J1+J)
   60    CONTINUE
         PIJW(IPIJ+(I-1)*J2+I)=1.0+PIJW(IPIJ+(I-1)*J2+I)
   70    CONTINUE
         CALL ALINV(J2,PIJW(IPIJ+1),J2,IER)
         IF(IER.NE.0) CALL XABORT('SYB001: SINGULAR MATRIX.')
         ALLOCATE(WORK(J2))
         DO 120 I=1,J2
         DO 80 K=1,J2
         WORK(K)=PIJW(IPIJ+(K-1)*J2+I)
   80    CONTINUE
         DO 100 J=1,J2
         WGAR=0.0
         DO 90 K=1,J2
         WGAR=WGAR+WORK(K)*PP(K,J)
   90    CONTINUE
         PIJW(IPIJ+(J-1)*J2+I)=WGAR
  100    CONTINUE
         WGAR=0.0
         DO 110 K=1,J2
         WGAR=WGAR+WORK(K)*PIS(K)
  110    CONTINUE
         PISW(J1+I)=WGAR
  120    CONTINUE
         DEALLOCATE(WORK)
*
*        COMPUTE THE SCATTERING-REDUCED COLLISION PROBABILITY MATRIX
*        FOR INCOMING NEUTRONS.
         DO 140 J=1,J2
         WGAR=PSJ(J)
         DO 130 K=1,J2
         WGAR=WGAR+PSJ(K)*SIGW(J1+K)*PIJW(IPIJ+(J-1)*J2+K)
  130    CONTINUE
         PSJW(J1+J)=WGAR
  140    CONTINUE
*
*        COMPUTE THE SCATTERING-REDUCED TRANSMISSION PROBABILITY MATRIX.
         WGAR=PSS
         DO 150 K=1,J2
         WGAR=WGAR+PSJ(K)*SIGW(J1+K)*PISW(J1+K)
  150    CONTINUE
         PSSW(IKK)=WGAR
      ENDIF
      DEALLOCATE(PP)
      IF(IMPX.GE.10) THEN 
         CALL SYBPRX(2,1,J2,IKK,SIGT(J1+1),SIGW(J1+1),PIJW(IPIJ+1),
     1   PISW(J1+1),PSJW(J1+1),PSSW(J1+1))
      ENDIF
      IPIJ=IPIJ+J2*J2
  160 CONTINUE
*----
*  SCRATCH STORAGE DEALLOCATION
*----
      DEALLOCATE(PSJ,PIS)
      RETURN
      END