4 *=== fluscw ===========================================================*
6 DOUBLE PRECISION FUNCTION FLUSCW
7 #(IJ,PLA,TXX,TYY,TZZ,WEE,XX,YY,ZZ,NREG,IOLREG,LLO,NSURF)
20 *----------------------------------------------------------------------*
22 * This functions returns neutron, proton and pion displacement damage *
25 *----------------------------------------------------------------------*
29 * Ij = (generalized) particle code *
30 * Pla = particle momentum (if > 0), or kinetic energy (if <0 )*
31 * Txx,yy,zz = particle direction cosines *
32 * Wee = particle weight *
33 * Xx,Yy,Zz = position *
34 * Nreg = (new) region number *
35 * Iolreg = (old) region number *
36 * Llo = particle generation *
37 * Nsurf = transport flag (ignore!) *
41 * Fluscw = factor the scored amount will be multiplied by *
42 * Lsczer = logical flag, if true no amount will be scored *
43 * regardless of Fluscw *
45 * Useful variables (common SCOHLP): *
47 * Flux like binnings/estimators (Fluscw): *
48 * ISCRNG = 1 --> Boundary crossing estimator *
49 * ISCRNG = 2 --> Track length binning *
50 * ISCRNG = 3 --> Track length estimator *
51 * ISCRNG = 4 --> Collision density estimator *
52 * ISCRNG = 5 --> Yield estimator *
53 * JSCRNG = # of the binning/estimator *
55 *----------------------------------------------------------------------*
62 * the default is unit weight
66 * skip all scorings other than tracklength
67 IF (ISCRNG.NE.2) RETURN
68 * skip all particles other than protons, pions and neutrons
69 IF ((IJ.NE.1).AND.(IJ.NE.13).AND.(IJ.NE.14).AND.(IJ.NE.8)) RETURN
71 * read displacement damage factors
74 1000 FORMAT(1X,'FLUSCW: direct conversion of neutron fluence to',
75 & ' 1 MeV neutron equivalent requested')
76 WRITE(LUNOUT,*) ' neutrons'
77 OPEN(19,FILE='disdam_n.dat',STATUS='UNKNOWN')
79 CALL RDXSC(19,LUNOUT,IDNDAM,-3)
81 WRITE(LUNOUT,*) ' protons '
82 OPEN(19,FILE='disdam_p.dat',STATUS='UNKNOWN')
84 CALL RDXSC(19,LUNOUT,IDPDAM,-2)
86 WRITE(LUNOUT,*) ' pions '
87 OPEN(19,FILE='disdam_pi.dat',STATUS='UNKNOWN')
89 CALL RDXSC(19,LUNOUT,IDODAM,-2)
94 * should be always called with ekin ( pla < 0 ),
95 * but we leave the check for the moment..
96 IF (PLA.LT.0.0D0) THEN
98 ELSEIF (PLA.GT.0.0D0) THEN
99 EKIN = SQRT(PLA**2+AM(IJ)**2)-AM(IJ)
104 * calculate the weight
106 FLUSCW = XSECT(IDPDAM,EKIN,0)
107 ELSEIF ((IJ.EQ.13).OR.(IJ.EQ.14)) THEN
108 FLUSCW = XSECT(IDODAM,EKIN,0)
109 ELSEIF (IJ.EQ.8) THEN
110 FLUSCW = XSECT(IDNDAM,EKIN,0)
112 STOP ' FLUSCW: inconsistent IJ '
118 *===rdxsc==============================================================*
120 SUBROUTINE RDXSC(LINP,LCHK,IDXSEC,MODE)
122 ************************************************************************
123 * LINP logical input unit *
124 * LCHK > 0 cross section data are plotted to unit LCHK *
125 * in equidistant log. binning using double-log. *
127 * (otherwise they are not plotted) *
128 * IDXSEC index of stored cross section in common *
129 * |MODE| = 1 cross section data given as histogram *
130 * i.e. E_lo(i) sigma(i) *
132 * with increasing energy *
133 * = 2 cross section data given for bin averages with *
134 * increasing energy *
135 * = 3 cross section data given for bin averages with *
136 * decreasing energy *
137 * if MODE < 0 the energy unit is assumed to be MeV (otherwise GeV) *
138 ************************************************************************
140 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
143 PARAMETER (MAXSDT = 10,
145 COMMON /CSECT/ XSEAV(MAXSDT,MAXSBI),XS(MAXSDT,MAXSBI),
146 & NXSBIN(MAXSDT),NXSDT
147 DIMENSION TMPXS(2,MAXSBI)
157 XSEAV(IDT,IBIN) = 0.0D0
169 IF (NEBIN.GT.MAXSBI) STOP ' nebin > maxsbi !!!'
170 IF (IABS(MODE).EQ.1) THEN
171 READ(LINP,*,END=9000) ELOW,XS(NXSDT,NEBIN)
172 READ(LINP,*) ELHI,XS(NXSDT,NEBIN)
173 XSEAV(NXSDT,NEBIN) = SQRT(ELOW*ELHI)
174 IF (MODE.LT.0) XSEAV(NXSDT,NEBIN) = 1.D-3*XSEAV(NXSDT,NEBIN)
175 ELSEIF (IABS(MODE).EQ.2) THEN
176 READ(LINP,*,END=9000) XSEAV(NXSDT,NEBIN),XS(NXSDT,NEBIN)
177 IF (MODE.LT.0) XSEAV(NXSDT,NEBIN) = 1.D-3*XSEAV(NXSDT,NEBIN)
178 ELSEIF (IABS(MODE).EQ.3) THEN
179 READ(LINP,*,END=9000) TMPXS(1,NEBIN),TMPXS(2,NEBIN)
180 IF (MODE.LT.0) TMPXS(1,NEBIN) = 1.D-3*TMPXS(1,NEBIN)
182 STOP ' RDXSC: unsupported mode ! '
187 IF (IABS(MODE).EQ.3) THEN
190 XSEAV(NXSDT,IDX) = TMPXS(1,I)
191 XS(NXSDT,IDX) = TMPXS(2,I)
195 NXSBIN(NXSDT) = NEBIN-1
199 AELO = LOG10(XSEAV(IDXSEC,1))
200 AEHI = LOG10(XSEAV(IDXSEC,NXSBIN(IDXSEC)))
201 DAE = (AEHI-AELO)/DBLE(NXSBIN(IDXSEC))
202 DO 4 I=1,NXSBIN(IDXSEC)+1
203 AE = AELO+DBLE(I-1)*DAE
205 WRITE(LCHK,'(1X,2E15.5)') E,XSECT(IDXSEC,E,0)
212 *===xsect==============================================================*
214 DOUBLE PRECISION FUNCTION XSECT(IDXSEC,E,MODE)
216 ************************************************************************
217 * IDXSEC index of stored cross section in common *
219 ************************************************************************
221 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
224 PARAMETER (MAXSDT = 10,
226 COMMON /CSECT/ XSEAV(MAXSDT,MAXSBI),XS(MAXSDT,MAXSBI),
227 & NXSBIN(MAXSDT),NXSDT
229 IF (E.LE.XSEAV(IDXSEC,1)) THEN
231 ELSEIF (E.GT.XSEAV(IDXSEC,NXSBIN(IDXSEC))) THEN
233 CS = XSCINT(IDXSEC,E,N)
235 DO 1 J=1,NXSBIN(IDXSEC)-1
236 IF ((E.GT.XSEAV(IDXSEC,J)).AND.(E.LE.XSEAV(IDXSEC,J+1)))
238 CS = XSCINT(IDXSEC,E,J)
242 STOP ' xsection value not found '
250 *===xscint=============================================================*
252 DOUBLE PRECISION FUNCTION XSCINT(IDXSEC,E,J)
254 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
257 PARAMETER (MAXSDT = 10,
259 COMMON /CSECT/ XSEAV(MAXSDT,MAXSBI),XS(MAXSDT,MAXSBI),
260 & NXSBIN(MAXSDT),NXSDT
262 FAC = (LOG10(E) -LOG10(XSEAV(IDXSEC,J)))/
263 & (LOG10(XSEAV(IDXSEC,J+1))-LOG10(XSEAV(IDXSEC,J)))
264 XSCINT = LOG10(XS(IDXSEC,J))
265 & +(LOG10(XS(IDXSEC,J+1))-LOG10(XS(IDXSEC,J)))*FAC
266 XSCINT = 10**(XSCINT)
271 *===skip===============================================================*
273 SUBROUTINE SKIP(LUNIT,NSKIP)
275 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
277 PARAMETER (LOUT=6,LLOOK=9)
281 IF (NSKIP.EQ.0) RETURN
284 READ(LUNIT,'(A132)') ACARD