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
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
|
*DECK HEADRV
SUBROUTINE HEADRV(IPDEP,NPART,IPMAC,NBMIX,NGRP,ZNORM,IMPX,ESUM,
1 CSUM,IBC,RHO)
*
*-----------------------------------------------------------------------
*
*Purpose:
* Compute energy and charge deposition from many particles.
*
*Copyright:
* Copyright (C) 2020 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
* IPDEP L_DEPOSITION pointer to the deposition information object.
* NPART number of particles contributing to energy and charge
* deposition.
* IPMAC L_MACROLIB pointers to the extended macrolibs.
* NBMIX number of material mixtures.
* NGRP total number of energy groups.
* ZNORM flux normalization factor.
* IMPX print parameter.
*
*Parameters: output
* ESUM total energy deposition (MeV/cc).
* CSUM total charge deposition (electron/cc).
*
*-----------------------------------------------------------------------
*
USE GANLIB
*----
* SUBROUTINE ARGUMENTS
*----
INTEGER NPART,NBMIX,NGRP,IMPX
TYPE(C_PTR) IPDEP,IPMAC(NPART)
REAL RHO(NBMIX)
DOUBLE PRECISION ZNORM,ESUM,CSUM
*----
* LOCAL VARIABLES
*----
TYPE(C_PTR) JPMAC,KPMAC
PARAMETER(NSTATE=40,IOUT=6)
INTEGER ISTATE(NSTATE),MAT(NBMIX)
CHARACTER HSMG*131,TEXT1*1
DOUBLE PRECISION VTOT
LOGICAL LCHARG
REAL FLUXC(NBMIX)
*----
* ALLOCATABLE ARRAYS
*----
REAL, ALLOCATABLE, DIMENSION(:) :: VOL,SGD,FLIN,ESTOPW
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:) :: EDEPOT,CDEPOT
DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:) :: EDEPO,CDEPO
CHARACTER(LEN=1), ALLOCATABLE, DIMENSION(:) :: TEXT1V
CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:) :: SNAME
CHARACTER(LEN=8), ALLOCATABLE, DIMENSION(:) :: FUNA8
*----
* MEMORY ALLOCATION
*----
ALLOCATE(EDEPO(NBMIX,NPART),CDEPO(NBMIX,NPART),EDEPOT(NBMIX),
1 CDEPOT(NBMIX),VOL(NBMIX),SGD(NBMIX),FLIN(NBMIX))
ALLOCATE(FUNA8(2*NPART+2),SNAME(2*NPART+2))
*----
* RECOVER ENERGY AND CHARGE DEPOSITION
*----
CALL LCMLEN(IPDEP,'EDEPOS',LENGT,ITYLCM)
IF(LENGT.NE.0) THEN
IF(LENGT.NE.NBMIX*NPART) CALL XABORT('HEADRV: INVALID EDEPOS R'
1 //'ECORD LENGTH.')
CALL LCMGET(IPDEP,'EDEPOS',EDEPO)
CALL LCMGET(IPDEP,'FLUX-NORM',ZNORM)
ELSE
EDEPO(:NBMIX,:NPART)=0.0D0
ENDIF
CALL LCMLEN(IPDEP,'CDEPOS',LENGT,ITYLCM)
IF(LENGT.NE.0) THEN
IF(LENGT.NE.NBMIX*NPART) CALL XABORT('HEADRV: INVALID CDEPOS R'
1 //'ECORD LENGTH.')
CALL LCMGET(IPDEP,'CDEPOS',CDEPO)
ELSE
CDEPO(:NBMIX,:NPART)=0.0D0
ENDIF
EDEPOT(:NBMIX)=0.0D0
CDEPOT(:NBMIX)=0.0D0
CHARGE=0.0
DO I=1,NPART
CALL LCMLEN(IPMAC(I),'FLUXC',IBC2,ITYLCM)
IF(IBC.EQ.1.AND.IBC2.NE.0) THEN
CALL LCMGET(IPMAC(I),'FLUXC',FLUXC)
CALL LCMGET(IPMAC(I),'ECUTOFF',ECUTOFF)
CALL LCMLEN(IPMAC(I),'ESTOPW',LENGT,ITYLCM)
ALLOCATE(ESTOPW(LENGT))
CALL LCMGET(IPMAC(I),'ESTOPW',ESTOPW)
CALL LCMGET(IPMAC(I),'MATCOD',MAT)
ENDIF
CALL LCMGET(IPMAC(I),'VOLUME',VOL)
CALL LCMGTC(IPMAC(I),'PARTICLE',1,TEXT1)
SNAME(I)=TEXT1
IF(TEXT1.EQ.'N'.OR.TEXT1.EQ.'G') THEN
CHARGE=0.0
ELSEIF(TEXT1.EQ.'C'.OR.TEXT1.EQ.'P') THEN
CHARGE=1.0
ELSEIF(TEXT1.EQ.'B') THEN
CHARGE=-1.0
ELSE
CALL XABORT('HEADRV: UNKNOWN PARTICLE.')
ENDIF
FUNA8(I)='ENERGDEP'
FUNA8(NPART+I+1)='CHARGDEP'
SNAME(NPART+I+1)=SNAME(I)
JPMAC=LCMGID(IPMAC(I),'GROUP')
DO IGR=1,NGRP
KPMAC=LCMGIL(JPMAC,IGR)
CALL LCMLEN(KPMAC,'H-FACTOR',LENGT,ITYLCM)
IF(LENGT.EQ.0) THEN
WRITE(HSMG,'(42HHEADRV: NO H-FACTOR FOUND IN MACROLIB NUMB,
1 2HER,I3,1H.)') I
CALL XABORT(HSMG)
ENDIF
CALL LCMGET(KPMAC,'FLUX-INTG',FLIN)
CALL LCMGET(KPMAC,'H-FACTOR',SGD)
DO IBM=1,NBMIX
EDEPO(IBM,I)=EDEPO(IBM,I)+FLIN(IBM)*SGD(IBM)*ZNORM/
1 (VOL(IBM)*RHO(IBM))
IF(IBC.EQ.1.AND.IBC2.NE.0) THEN
EDEPO(IBM,I)=EDEPO(IBM,I)+ECUTOFF*ESTOPW(MAT(IBM))
1 *FLUXC(IBM)*ZNORM/RHO(IBM)
ENDIF
ENDDO
CALL LCMLEN(KPMAC,'C-FACTOR',LENGT,ITYLCM)
IF(LENGT.GT.0) THEN
LCHARG=.TRUE.
CALL LCMGET(KPMAC,'C-FACTOR',SGD)
DO IBM=1,NBMIX
CDEPO(IBM,I)=CDEPO(IBM,I)+FLIN(IBM)*SGD(IBM)*ZNORM/
1 (VOL(IBM)*RHO(IBM))
IF(IBC.EQ.1.AND.IBC2.NE.0) THEN
CDEPO(IBM,I)=CDEPO(IBM,I)+ESTOPW(MAT(IBM))*FLUXC(IBM)
1 *ZNORM/RHO(IBM)*(-CHARGE)
ENDIF
ENDDO
ENDIF
ENDDO
IF(IBC.EQ.1.AND.IBC2.NE.0) DEALLOCATE(ESTOPW)
ENDDO
FUNA8(NPART+1)='TOTENDEP'
FUNA8(2*NPART+2)='TOTCHDEP'
SNAME(NPART+1)='EDEPO'
SNAME(2*NPART+2)='CDEPO'
VTOT=0.0D0
ESUM=0.0D0
CSUM=0.0D0
DO IBM=1,NBMIX
DO I=1,NPART
EDEPOT(IBM)=EDEPOT(IBM)+EDEPO(IBM,I)
CDEPOT(IBM)=CDEPOT(IBM)+CDEPO(IBM,I)
ENDDO
VTOT=VTOT+VOL(IBM)
ESUM=ESUM+EDEPOT(IBM)*VOL(IBM)
CSUM=CSUM+CDEPOT(IBM)*VOL(IBM)
ENDDO
ESUM=ESUM/VTOT
CSUM=CSUM/VTOT
*----
* PRINT ENERGY AND CHARGE DEPOSITION
*----
IF(IMPX.GT.0) THEN
WRITE(IOUT,1001) ' VOLUME ',(FUNA8(J),'_',SNAME(J),J=1,
1 2*NPART+2)
DO IBM=1,NBMIX
WRITE(IOUT,1002) VOL(IBM),(EDEPO(IBM,I),I=1,NPART),
1 EDEPOT(IBM),(CDEPO(IBM,I),I=1,NPART),CDEPOT(IBM)
ENDDO
WRITE(IOUT,'(/14H TOTAL VOLUME:,14X,1P,E12.4)') VTOT
WRITE(IOUT,'(25H TOTAL ENERGY DEPOSITION:,3X,1P,E12.4)') ESUM
WRITE(IOUT,'(25H TOTAL CHARGE DEPOSITION:,3X,1P,E12.4)') CSUM
ENDIF
*----
* SAVE ENERGY AND CHARGE DEPOSITION
*----
CALL LCMPUT(IPDEP,'VOLUME',NBMIX,2,VOL)
CALL LCMPUT(IPDEP,'EDEPOS',NBMIX*NPART,4,EDEPO)
CALL LCMPUT(IPDEP,'EDEPOS_TOT',NBMIX,4,EDEPOT)
IF(LCHARG) THEN
CALL LCMPUT(IPDEP,'CDEPOS',NBMIX*NPART,4,CDEPO)
CALL LCMPUT(IPDEP,'CDEPOS_TOT',NBMIX,4,CDEPOT)
ENDIF
CALL LCMPUT(IPDEP,'FLUX-NORM',1,4,ZNORM)
*----
* PROCESS STATE-VECTOR
*----
CALL LCMLEN(IPDEP,'STATE-VECTOR',LENGT,ITYLCM)
IF(LENGT.NE.0) THEN
CALL LCMGET(IPDEP,'STATE-VECTOR',ISTATE)
ALLOCATE(TEXT1V(NPART))
CALL LCMGTC(IPDEP,'PARTICLE-NAM',1,NPART,TEXT1V)
DO I=1,NPART
IF(TEXT1V(I).NE.SNAME(I)(:1)) THEN
WRITE(HSMG,'(22HHEADRV: PARTICLE NAMES,2A2,
1 16H ARE INCOHERENT.)') TEXT1V(I),SNAME(I)(:1)
CALL XABORT(HSMG)
ENDIF
ENDDO
DEALLOCATE(TEXT1V)
ELSE
ISTATE(:NSTATE)=0
ISTATE(1)=NBMIX
ISTATE(2)=NPART
IF(LCHARG) ISTATE(3)=1
ALLOCATE(TEXT1V(NPART))
DO I=1,NPART
TEXT1V(I)=SNAME(I)(:1)
ENDDO
CALL LCMPTC(IPDEP,'PARTICLE-NAM',1,NPART,TEXT1V)
DEALLOCATE(TEXT1V)
ENDIF
ISTATE(4)=ISTATE(4)+1
CALL LCMPUT(IPDEP,'STATE-VECTOR',NSTATE,1,ISTATE)
*----
* MEMORY DEALLOCATION
*----
DEALLOCATE(SNAME,FUNA8)
DEALLOCATE(FLIN,SGD,VOL,CDEPOT,EDEPOT,CDEPO,EDEPO)
RETURN
*
1001 FORMAT(/1X,A11,21(1X,A8,A1,A6))
1002 FORMAT(1X,1P,E11.4,21E16.4)
END
|
