TFluka/comscw_activity.f

   1 *$ CREATE COMSCW.FOR
   2 *COPY COMSCW
   3 *
   4 *=== comscw  ===========================================================*
   5 *
   6       DOUBLE PRECISION FUNCTION COMSCW ( IJ    , XA    , YA    , ZA    ,
   7      &                                   MREG  , RULL  , LLO   , ICALL )
   8
   9       INCLUDE '(DBLPRC)'
  10       INCLUDE '(DIMPAR)'
  11       INCLUDE '(IOUNIT)'
  12 *                                                                      *
  13 *----------------------------------------------------------------------*
  14 *                                                                      *
  15 * Special multiplication factors for activity scoring.                 *
  16 *                                                                      *
  17 * Note: Data files 'LE-CH.dat' and 'LE-EC.dat' read on unit 20.        *
  18 *                                                                      *
  19 *                                                                      *
  20 * SDUM = 'abcdefgh'                                              Lesco *
  21 *                                                                      *
  22 *   'ab' = 'LE'  division by exemption limits activated                *
  23 *                                                                      *
  24 *        'c' = 'C' : division by CERN zonage limits                    *
  25 *            'def' = '10 '    : division by 1/10 of the value      10  *
  26 *            'def' = '100'    : division by 1/100 of the value    100  *
  27 *                    otherwise: division by the value itself        1  *
  28 *        'c' = 'S' : division by Swiss LE values                    2  *
  29 *        'c' = 'E' : division by European concentration limits      3  *
  30 *                                                                      *
  31 *                Comscw = 1. for selected isotope                      *
  32 *                       = 0. otherwise                                 *
  33 *                                                                      *
  34 *   'ab' = 'IS'  selection of individual isotops activated       9999  *
  35 *                                                                      *
  36 *        'cd'  = symbol of isotope                                     *
  37 *        'efg' = mass number of isotope                                *
  38 *        'h'   = 'm' : metastable state                                *
  39 *                                                                      *
  40 *                Comscw = 1. for selected isotope                      *
  41 *                       = 0. otherwise                                 *
  42 *                                                                      *
  43 *   otherwise    Comscw = 1.                                       -1  *
  44 *                                                                      *
  45 *   Example(s):                                                        *
  46 *                                                                      *
  47 *    SDUM = 'LEC10xxx' : division by 1/10 of LE(EU) or by Swiss LE     *
  48 *         = 'LESxxxxx' : division by Swiss LE                          *
  49 *         = 'ISV  48x' : scoring of activity for V48                   *
  50 *         = 'ISMn 52m' : scoring of activity for mMn52                 *
  51 *                                                                      *
  52 *    Note: 'x' refers to characters without specific meaning           *
  53 *                                                                      *
  54 *     Input variables:                                                 *
  55 *                                                                      *
  56 *           Ij = (generalized) particle code                           *
  57 *     Xa,Ya,Za = position                                              *
  58 *         Mreg = region number                                         *
  59 *         Rull = amount to be deposited                                *
  60 *          Llo = particle generation                                   *
  61 *        Icall = call id                                               *
  62 *                                                                      *
  63 *     Output variables:                                                *
  64 *                                                                      *
  65 *       Comscw = factor the scored amount will be multiplied by        *
  66 *       Lsczer = logical flag, if true no amount will be scored        *
  67 *                regardless of Comscw                                  *
  68 *                                                                      *
  69 *     Useful variables (common SCOHLP):                                *
  70 *                                                                      *
  71 *     Energy/Star binnings/scorings (Comscw):                          *
  72 *          ISCRNG = 1 --> Energy density  binning                      *
  73 *          ISCRNG = 2 --> Star   density  binning                      *
  74 *          ISCRNG = 3 --> Residual nuclei scoring                      *
  75 *          ISCRNG = 4 --> Momentum transfer density binning            *
  76 *          ISCRNG = 5 --> Activity density binning                     *
  77 *          JSCRNG = # of the binning                                   *
  78 *                                                                      *
  79 *     Useful variables (common SOUEVT):                                *
  80 *                                                                      *
  81 *          X,Y,Zsoevt(i) = position    of the i_th source particle     *
  82 *          TX,Y,Zsoev(i) = direction   of the i_th source particle     *
  83 *              Wtsoev(i) = weight      of the i_th source particle     *
  84 *              Pmsoev(i) = momentum    of the i_th source particle     *
  85 *              Tksoev(i) = kin. energy of the i_th source particle     *
  86 *              Agsoev(i) = age         of the i_th source particle     *
  87 *              Aksoev(i) = Kaon ampl.  of the i_th source particle     *
  88 *              Ussoev(i) = user var.   of the i_th source particle     *
  89 *              Ijsoev(i) = identity    of the i_th source particle     *
  90 *              Nrsoev(i) = region      of the i_th source particle     *
  91 *              Nlsoev(i) = lattice     of the i_th source particle     *
  92 *                Npsoev  = number of the source particles              *
  93 *----------------------------------------------------------------------*
  94 *
  95       INCLUDE '(FLKMAT)'
  96       INCLUDE '(SCOHLP)'
  97       INCLUDE '(SOUEVT)'
  98 *
  99       INCLUDE '(RSNCCM)'
 100       INCLUDE '(USRBIN)'
 101       DIMENSION IZSCO(MXUSBN),IASCO(MXUSBN),ISSCO(MXUSBN),
 102      &          LESCO(MXUSBN)
 103 *
 104       CHARACTER CISOIN*2,CISO*2,CSET*10,CDUM*4,CA*3
 105       PARAMETER (IZMAX  = 109,
 106      &           IAZMIN = 2,
 107      &           IAZMAX = 160)
 108       DIMENSION CISO(IZMAX)
 109       DIMENSION XLESWS(IZMAX,IAZMIN:IAZMAX,2),
 110      &          XLEECO(IZMAX,IAZMIN:IAZMAX,2),
 111      &          XLEO10(IZMAX,IAZMIN:IAZMAX,2)
 112       DATA CISO /
 113      &    'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne','Na',
 114      &    'Mg','Al','Si','P ','S ','Cl','Ar','K ','Ca','Sc','Ti',
 115      &    'V ','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As',
 116      &    'Se','Br','Kr','Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru',
 117      &    'Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I ','Xe','Cs',
 118      &    'Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy',
 119      &    'Ho','Er','Tm','Yb','Lu','Hf','Ta','W ','Re','Os','Ir',
 120      &    'Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra',
 121      &    'Ac','Th','Pa','U ','Np','Pu','Am','Cm','Bk','Cf','Es',
 122      &    'Fm','Md','No','Lr','Rf','Ha','Sg','Ns','Hs','Mt'/
 123 *
 124       LOGICAL LFIRST,LFOUND
 125       DATA LFIRST /.TRUE./
 126 *
 127       LSCZER = .FALSE.
 128       COMSCW = ONEONE
 129
 130       IF (LFIRST) THEN
 131
 132          WRITE(LUNOUT,'(A)') ' COMSCW: activity weighting activated'
 133          LFIRST = .FALSE.
 134 *
 135 *----------------------------------------------------------------------*
 136 * Initialization
 137 *
 138 * LE data taken from Appendix 3 column 9,
 139 * Ordonnance sur la radioprotection (ORaP) du 22 juin 1994
 140 * (etat au 4 avril 2000)
 141 *
 142 * (data file in Bq/kg!)
 143 *
 144          DO 104 IZ=1,IZMAX
 145             DO 105 IAZ=IAZMIN,IAZMAX
 146                XLESWS(IZ,IAZ,1) = ZERZER
 147                XLESWS(IZ,IAZ,2) = ZERZER
 148   105       CONTINUE
 149   104    CONTINUE
 150 *
 151          OPEN(20,FILE='LE-CH.dat',STATUS='UNKNOWN')
 152          NISO   = 0
 153          IZ0    = 1
 154   100    CONTINUE
 155          READ(20,*,END=101) CISOIN,IA,XLIMIT
 156 * convert from Bq/kg into Bq/g
 157          XLIMIT = 1.0D-3*XLIMIT
 158          LFOUND = .FALSE.
 159          DO 102 IZ=IZ0,IZMAX
 160             IF (CISOIN.EQ.CISO(IZ)) THEN
 161                IF (IA.LT.0) THEN
 162                   IS = 2
 163                   IA = ABS(IA)
 164                ELSE
 165                   IS = 1
 166                ENDIF
 167                IAZ = IA-IZ
 168                IF ((IAZ.LT.IAZMIN).OR.(IAZ.GT.IAZMAX)) THEN
 169                   WRITE(LUNOUT,*)
 170      &               ' COMSCW: warning! Iaz out of allowed range: ',
 171      &               IAZ, IAZMIN,IAZMAX
 172                   STOP
 173                ENDIF
 174                IF (XLESWS(IZ,IAZ,IS).GT.ZERZER) THEN
 175                   WRITE(LUNOUT,*)
 176      &               ' COMSCW: warning! two entries for this isotope: ',
 177      &               CISOIN,IA,XLIMIT
 178                   STOP
 179                ELSE
 180                   XLESWS(IZ,IAZ,IS) = XLIMIT
 181                ENDIF
 182                NISO = NISO+1
 183                IZ0  = IZ
 184                LFOUND = .TRUE.
 185                GOTO 103
 186             ENDIF
 187   102    CONTINUE
 188          WRITE(LUNOUT,*)
 189      &      ' COMSCW: isotope not recognized: ',CISOIN,IA,XLIMIT
 190          STOP
 191   103    CONTINUE
 192          GOTO 100
 193   101    CONTINUE
 194          CLOSE(20)
 195 *
 196 * LE data taken from
 197 * Official journal of the European Communities L159, 29 June 1996
 198 * Council Directive 96/29/Euratom
 199 *
 200 * (data file in Bq/g!)
 201 *
 202          DO 204 IZ=1,IZMAX
 203             DO 205 IAZ=IAZMIN,IAZMAX
 204                XLEECO(IZ,IAZ,1) = ZERZER
 205                XLEECO(IZ,IAZ,2) = ZERZER
 206   205       CONTINUE
 207   204    CONTINUE
 208 *
 209          OPEN(20,FILE='LE-EC.dat',STATUS='UNKNOWN')
 210          NISO   = 0
 211          IZ0    = 1
 212   200    CONTINUE
 213          READ(20,*,END=201) CISOIN,IA,XLIMIT,IFLAG
 214          LFOUND = .FALSE.
 215          DO 202 IZ=IZ0,IZMAX
 216             IF (CISOIN.EQ.CISO(IZ)) THEN
 217                IF (IA.LT.0) THEN
 218                   IS = 2
 219                   IA = ABS(IA)
 220                ELSE
 221                   IS = 1
 222                ENDIF
 223                IAZ = IA-IZ
 224                IF ((IAZ.LT.IAZMIN).OR.(IAZ.GT.IAZMAX)) THEN
 225                   WRITE(LUNOUT,*)
 226      &               ' COMSCW: warning! Iaz out of allowed range: ',
 227      &               IAZ, IAZMIN,IAZMAX
 228                   STOP
 229                ENDIF
 230                IF (XLEECO(IZ,IAZ,IS).GT.ZERZER) THEN
 231                   WRITE(LUNOUT,*)
 232      &               ' COMSCW: warning! two entries for this isotope: ',
 233      &               CISOIN,IA,XLIMIT
 234                   STOP
 235                ELSE
 236                   XLEECO(IZ,IAZ,IS) = XLIMIT
 237 * zero entries with Swiss values
 238                   IF (IFLAG.EQ.1) XLEECO(IZ,IAZ,IS) = ZERZER
 239                ENDIF
 240                NISO = NISO+1
 241                IZ0  = IZ
 242                LFOUND = .TRUE.
 243                GOTO 203
 244             ENDIF
 245   202    CONTINUE
 246          WRITE(LUNOUT,*)
 247      &      ' COMSCW: isotope not recognized: ',CISOIN,IA,XLIMIT
 248          STOP
 249   203    CONTINUE
 250          GOTO 200
 251   201    CONTINUE
 252          CLOSE(20)
 253          DO 404 IZ=1,IZMAX
 254             DO 405 IAZ=IAZMIN,IAZMAX
 255                DO 406 IS=1,2
 256                   XLEO10(IZ,IAZ,IS) = ZERZER
 257                   IF (XLEECO(IZ,IAZ,IS).GT.ZERZER) THEN
 258                      XLEO10(IZ,IAZ,IS) = XLEECO(IZ,IAZ,IS)/10.0D0
 259                   ELSE
 260                      XLEO10(IZ,IAZ,IS) = XLESWS(IZ,IAZ,IS)
 261                   ENDIF
 262   406          CONTINUE
 263   405       CONTINUE
 264   404    CONTINUE
 265 *
 266          DO 500 I=1,MXUSBN
 267             IZSCO(I) = 0
 268             IASCO(I) = 0
 269             ISSCO(I) = 0
 270             LESCO(I) = 0
 271   500    CONTINUE
 272
 273       ENDIF
 274 *
 275 *----------------------------------------------------------------------*
 276 * Online weighting
 277 *
 278       IF ( ISCRNG .EQ. 5 ) THEN
 279 *
 280 * determine type of weighting
 281          IF (LESCO(JSCRNG).EQ.0) THEN
 282             CSET = TITUSB(JSCRNG)
 283             IF (CSET(1:2).EQ.'IS') THEN
 284                LESCO(JSCRNG) = 9999
 285                DO 700 IZ=1,IZMAX
 286                   IF (CSET(3:4).EQ.CISO(IZ)) THEN
 287                      IZSCO(JSCRNG) = IZ
 288                      READ(CSET,'(A4,A3)') CDUM,CA
 289                      READ(CA,'(I3)') IASCO(JSCRNG)
 290                      IF (CSET(8:8).EQ.'m') THEN
 291                         ISSCO(JSCRNG) = 2
 292                      ELSE
 293                         ISSCO(JSCRNG) = 1
 294                      ENDIF
 295                   ENDIF
 296   700          CONTINUE
 297                IF ((IZSCO(JSCRNG).LE.0).OR.(IASCO(JSCRNG).LE.0)
 298      &                                 .OR.(ISSCO(JSCRNG).LE.0)) THEN
 299                   WRITE(LUNOUT,*) ' COMSCW: unknown isotope, Z,A,S = ',
 300      &               IZSCO(JSCRNG),IASCO(JSCRNG),ISSCO(JSCRNG)
 301                   STOP
 302                ENDIF
 303             ELSEIF(CSET(1:2).EQ.'LE') THEN
 304                IF (CSET(3:3).EQ.'C') THEN
 305                   IF (CSET(4:6).EQ.'100') THEN
 306                      LESCO(JSCRNG) = 100
 307                   ELSEIF (CSET(4:6).EQ.'10') THEN
 308                      LESCO(JSCRNG) = 10
 309                   ELSE
 310                      LESCO(JSCRNG) = 1
 311                   ENDIF
 312                ELSEIF (CSET(3:3).EQ.'S') THEN
 313                   LESCO(JSCRNG) = 2
 314                ELSEIF (CSET(3:3).EQ.'E') THEN
 315                   LESCO(JSCRNG) = 3
 316                ELSE
 317                   WRITE(LUNOUT,*) ' COMSCW: unknown LE set ',CSET(3:3)
 318                   STOP
 319                ENDIF
 320             ELSE
 321                LESCO(JSCRNG) = -1
 322             ENDIF
 323             WRITE(LUNOUT,1000) ' COMSCW: scoring ',JSCRNG,CSET,
 324      &         ' weighted with properties ',
 325      &         LESCO(JSCRNG),IZSCO(JSCRNG),IASCO(JSCRNG),ISSCO(JSCRNG)
 326  1000       FORMAT(A,I3,2A,I5,3I4)
 327          ENDIF
 328 *
 329 * obtain present isotope from common block
 330          JA  = IARSDL(1)
 331          JZ  = IZRSDL(1)
 332          JS  = ISRSDL(1)+1
 333          IF (JS.GT.2) JS = 2
 334          JAZ = JA-JZ
 335 *
 336          COMSCW = ZERZER
 337          IF (LESCO(JSCRNG).EQ.-1) THEN
 338             COMSCW = ONEONE
 339          ELSEIF (LESCO(JSCRNG).EQ.9999) THEN
 340             IF ((JA.EQ.IASCO(JSCRNG)).AND.(JZ.EQ.IZSCO(JSCRNG)).AND.
 341      &                                    (JS.EQ.ISSCO(JSCRNG))) THEN
 342                COMSCW = ONEONE
 343             ELSE
 344                COMSCW = ZERZER
 345             ENDIF
 346          ELSEIF ((LESCO(JSCRNG).EQ.1).OR.(LESCO(JSCRNG).EQ.10).OR.
 347      &           (LESCO(JSCRNG).EQ.100)) THEN
 348             FACT = 10.0D0/DBLE(LESCO(JSCRNG))
 349             IF (XLEO10(JZ,JAZ,JS).GT.ZERZER)
 350      &         COMSCW = ONEONE/(FACT*XLEO10(JZ,JAZ,JS))
 351          ELSEIF (LESCO(JSCRNG).EQ.2) THEN
 352             IF (XLESWS(JZ,JAZ,JS).GT.ZERZER)
 353      &         COMSCW = ONEONE/XLESWS(JZ,JAZ,JS)
 354          ELSEIF (LESCO(JSCRNG).EQ.3) THEN
 355             IF (XLEECO(JZ,JAZ,JS).GT.ZERZER)
 356      &         COMSCW = ONEONE/XLEECO(JZ,JAZ,JS)
 357          ELSE
 358             WRITE(LUNOUT,*) ' COMSCW: invalid option ',LESCO(JSCRNG)
 359             STOP
 360          ENDIF
 361
 362       ENDIF
 363
 364       RETURN
 365 *=== End of function Comscw ===========================================*
 366       END