TFluka/fluscw_1mevn.f

   1 *$ CREATE FLUSCW.FOR
   2 *COPY FLUSCW
   3 *                                                                      *
   4 *=== fluscw ===========================================================*
   5 *                                                                      *
   6       DOUBLE PRECISION FUNCTION FLUSCW
   7      #(IJ,PLA,TXX,TYY,TZZ,WEE,XX,YY,ZZ,NREG,IOLREG,LLO,NSURF)
   8
   9       INCLUDE '(DBLPRC)'
  10       INCLUDE '(DIMPAR)'
  11       INCLUDE '(IOUNIT)'
  12       SAVE
  13
  14       INCLUDE '(PAPROP)'
  15       INCLUDE '(USRBIN)'
  16       INCLUDE '(USRBDX)'
  17       INCLUDE '(USRTRC)'
  18       INCLUDE '(SCOHLP)'
  19 *
  20 *----------------------------------------------------------------------*
  21 *                                                                      *
  22 * This functions returns neutron, proton and pion  displacement damage *
  23 * weight factors.                                                      *
  24 *                                                                      *
  25 *----------------------------------------------------------------------*
  26 *                                                                      *
  27 *     Input variables:                                                 *
  28 *                                                                      *
  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!)                              *
  38 *                                                                      *
  39 *     Output variables:                                                *
  40 *                                                                      *
  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                                  *
  44 *                                                                      *
  45 *     Useful variables (common SCOHLP):                                *
  46 *                                                                      *
  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                         *
  54 *                                                                      *
  55 *----------------------------------------------------------------------*
  56 *
  57       LOGICAL LFIRST
  58
  59
  60       DATA LFIRST /.TRUE./
  61 *
  62 * the default is unit weight
  63       FLUSCW = ONEONE
  64       LSCZER = .FALSE.
  65 *
  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
  70 *
  71 * read displacement damage factors
  72       IF (LFIRST) THEN
  73          WRITE(LUNOUT,1000)
  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')
  78          CALL SKIP(19,7)
  79          CALL RDXSC(19,LUNOUT,IDNDAM,-3)
  80          CLOSE(19)
  81          WRITE(LUNOUT,*) ' protons '
  82          OPEN(19,FILE='disdam_p.dat',STATUS='UNKNOWN')
  83          CALL SKIP(19,7)
  84          CALL RDXSC(19,LUNOUT,IDPDAM,-2)
  85          CLOSE(19)
  86          WRITE(LUNOUT,*) ' pions   '
  87          OPEN(19,FILE='disdam_pi.dat',STATUS='UNKNOWN')
  88          CALL SKIP(19,7)
  89          CALL RDXSC(19,LUNOUT,IDODAM,-2)
  90          CLOSE(19)
  91          LFIRST = .FALSE.
  92       ENDIF
  93 *
  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
  97          EKIN = ABS(PLA)
  98       ELSEIF (PLA.GT.0.0D0) THEN
  99          EKIN = SQRT(PLA**2+AM(IJ)**2)-AM(IJ)
 100       ELSE
 101          RETURN
 102       ENDIF
 103 *
 104 * calculate the weight
 105       IF (IJ.EQ.1) THEN
 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)
 111       ELSE
 112          STOP ' FLUSCW: inconsistent IJ '
 113       ENDIF
 114
 115       RETURN
 116       END
 117 *
 118 *===rdxsc==============================================================*
 119 *
 120       SUBROUTINE RDXSC(LINP,LCHK,IDXSEC,MODE)
 121
 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.      *
 126 *                   interpolation                                      *
 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)                            *
 131 *                         E_hi(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 ************************************************************************
 139
 140       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 141       SAVE
 142
 143       PARAMETER (MAXSDT = 10,
 144      &           MAXSBI = 1400)
 145       COMMON /CSECT/ XSEAV(MAXSDT,MAXSBI),XS(MAXSDT,MAXSBI),
 146      &               NXSBIN(MAXSDT),NXSDT
 147       DIMENSION TMPXS(2,MAXSBI)
 148
 149       LOGICAL LFIRST
 150       DATA LFIRST /.TRUE./
 151
 152       IF (LFIRST) THEN
 153          NXSDT = 0
 154          DO 1 IDT=1,MAXSDT
 155             NXSBIN(IDT) = 0
 156             DO 2 IBIN=1,MAXSBI
 157                XSEAV(IDT,IBIN) = 0.0D0
 158                XS(IDT,IBIN)    = 0.0D0
 159     2       CONTINUE
 160     1    CONTINUE
 161          LFIRST = .FALSE.
 162       ENDIF
 163
 164       NXSDT = NXSDT+1
 165
 166       NEBIN = 0
 167    10 CONTINUE
 168          NEBIN = NEBIN+1
 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)
 181          ELSE
 182             STOP ' RDXSC: unsupported mode ! '
 183          ENDIF
 184          GOTO 10
 185  9000 CONTINUE
 186
 187       IF (IABS(MODE).EQ.3) THEN
 188          DO 3 I=1,NEBIN-1
 189             IDX = NEBIN-I
 190             XSEAV(NXSDT,IDX) = TMPXS(1,I)
 191             XS(NXSDT,IDX) = TMPXS(2,I)
 192     3    CONTINUE
 193       ENDIF
 194
 195       NXSBIN(NXSDT) = NEBIN-1
 196       IDXSEC = NXSDT
 197
 198       IF (LCHK.GT.0) THEN
 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
 204             E  = 10.0D0**AE
 205             WRITE(LCHK,'(1X,2E15.5)') E,XSECT(IDXSEC,E,0)
 206     4    CONTINUE
 207       ENDIF
 208
 209       RETURN
 210       END
 211 *
 212 *===xsect==============================================================*
 213 *
 214       DOUBLE PRECISION FUNCTION XSECT(IDXSEC,E,MODE)
 215
 216 ************************************************************************
 217 *   IDXSEC         index of stored cross section in common             *
 218 *   E              energy                                              *
 219 ************************************************************************
 220
 221       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 222       SAVE
 223
 224       PARAMETER (MAXSDT = 10,
 225      &           MAXSBI = 1400)
 226       COMMON /CSECT/ XSEAV(MAXSDT,MAXSBI),XS(MAXSDT,MAXSBI),
 227      &               NXSBIN(MAXSDT),NXSDT
 228
 229       IF (E.LE.XSEAV(IDXSEC,1)) THEN
 230          CS = 0.0D0
 231       ELSEIF (E.GT.XSEAV(IDXSEC,NXSBIN(IDXSEC))) THEN
 232          N  = NXSBIN(IDXSEC)-1
 233          CS = XSCINT(IDXSEC,E,N)
 234       ELSE
 235          DO 1 J=1,NXSBIN(IDXSEC)-1
 236             IF ((E.GT.XSEAV(IDXSEC,J)).AND.(E.LE.XSEAV(IDXSEC,J+1)))
 237      &                                                          THEN
 238                CS = XSCINT(IDXSEC,E,J)
 239                GOTO 2
 240             ENDIF
 241    1     CONTINUE
 242          STOP ' xsection value not found '
 243    2     CONTINUE
 244       ENDIF
 245       XSECT = CS
 246
 247       RETURN
 248       END
 249 *
 250 *===xscint=============================================================*
 251 *
 252       DOUBLE PRECISION FUNCTION XSCINT(IDXSEC,E,J)
 253
 254       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 255       SAVE
 256
 257       PARAMETER (MAXSDT = 10,
 258      &           MAXSBI = 1400)
 259       COMMON /CSECT/ XSEAV(MAXSDT,MAXSBI),XS(MAXSDT,MAXSBI),
 260      &               NXSBIN(MAXSDT),NXSDT
 261
 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)
 267
 268       RETURN
 269       END
 270 *
 271 *===skip===============================================================*
 272 *
 273       SUBROUTINE SKIP(LUNIT,NSKIP)
 274
 275       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
 276       SAVE
 277       PARAMETER (LOUT=6,LLOOK=9)
 278
 279       CHARACTER*132 ACARD
 280
 281       IF (NSKIP.EQ.0) RETURN
 282
 283       DO 1 K=1,NSKIP
 284          READ(LUNIT,'(A132)') ACARD
 285     1 CONTINUE
 286
 287       RETURN
 288       END