*$ CREATE FLUSCW.FOR *COPY FLUSCW * * *=== fluscw ===========================================================* * * DOUBLE PRECISION FUNCTION FLUSCW #(IJ,PLA,TXX,TYY,TZZ,WEE,XX,YY,ZZ,NREG,IOLREG,LLO,NSURF) INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' SAVE INCLUDE '(PAPROP)' INCLUDE '(USRBIN)' INCLUDE '(USRBDX)' INCLUDE '(USRTRC)' INCLUDE '(SCOHLP)' * *----------------------------------------------------------------------* * * * This functions returns neutron, proton and pion displacement damage * * weight factors. * * * *----------------------------------------------------------------------* * * * Input variables: * * * * Ij = (generalized) particle code * * Pla = particle momentum (if > 0), or kinetic energy (if <0 )* * Txx,yy,zz = particle direction cosines * * Wee = particle weight * * Xx,Yy,Zz = position * * Nreg = (new) region number * * Iolreg = (old) region number * * Llo = particle generation * * Nsurf = transport flag (ignore!) * * * * Output variables: * * * * Fluscw = factor the scored amount will be multiplied by * * Lsczer = logical flag, if true no amount will be scored * * regardless of Fluscw * * * * Useful variables (common SCOHLP): * * * * Flux like binnings/estimators (Fluscw): * * ISCRNG = 1 --> Boundary crossing estimator * * ISCRNG = 2 --> Track length binning * * ISCRNG = 3 --> Track length estimator * * ISCRNG = 4 --> Collision density estimator * * ISCRNG = 5 --> Yield estimator * * JSCRNG = # of the binning/estimator * * * *----------------------------------------------------------------------* * LOGICAL LFIRST DATA LFIRST /.TRUE./ * * the default is unit weight FLUSCW = ONEONE LSCZER = .FALSE. * * skip all scorings other than tracklength IF (ISCRNG.NE.2) RETURN * skip all particles other than protons, pions and neutrons IF ((IJ.NE.1).AND.(IJ.NE.13).AND.(IJ.NE.14).AND.(IJ.NE.8)) RETURN * * read displacement damage factors IF (LFIRST) THEN WRITE(LUNOUT,1000) 1000 FORMAT(1X,'FLUSCW: direct conversion of neutron fluence to', & ' 1 MeV neutron equivalent requested') WRITE(LUNOUT,*) ' neutrons' OPEN(19,FILE='disdam_n.dat',STATUS='UNKNOWN') CALL SKIP(19,7) CALL RDXSC(19,LUNOUT,IDNDAM,-3) CLOSE(19) WRITE(LUNOUT,*) ' protons ' OPEN(19,FILE='disdam_p.dat',STATUS='UNKNOWN') CALL SKIP(19,7) CALL RDXSC(19,LUNOUT,IDPDAM,-2) CLOSE(19) WRITE(LUNOUT,*) ' pions ' OPEN(19,FILE='disdam_pi.dat',STATUS='UNKNOWN') CALL SKIP(19,7) CALL RDXSC(19,LUNOUT,IDODAM,-2) CLOSE(19) LFIRST = .FALSE. ENDIF * * should be always called with ekin ( pla < 0 ), * but we leave the check for the moment.. IF (PLA.LT.0.0D0) THEN EKIN = ABS(PLA) ELSEIF (PLA.GT.0.0D0) THEN EKIN = SQRT(PLA**2+AM(IJ)**2)-AM(IJ) ELSE RETURN ENDIF * * calculate the weight IF (IJ.EQ.1) THEN FLUSCW = XSECT(IDPDAM,EKIN,0) ELSEIF ((IJ.EQ.13).OR.(IJ.EQ.14)) THEN FLUSCW = XSECT(IDODAM,EKIN,0) ELSEIF (IJ.EQ.8) THEN FLUSCW = XSECT(IDNDAM,EKIN,0) ELSE STOP ' FLUSCW: inconsistent IJ ' ENDIF RETURN END * *===rdxsc==============================================================* * SUBROUTINE RDXSC(LINP,LCHK,IDXSEC,MODE) ************************************************************************ * LINP logical input unit * * LCHK > 0 cross section data are plotted to unit LCHK * * in equidistant log. binning using double-log. * * interpolation * * (otherwise they are not plotted) * * IDXSEC index of stored cross section in common * * |MODE| = 1 cross section data given as histogram * * i.e. E_lo(i) sigma(i) * * E_hi(i) sigma(i) * * with increasing energy * * = 2 cross section data given for bin averages with * * increasing energy * * = 3 cross section data given for bin averages with * * decreasing energy * * if MODE < 0 the energy unit is assumed to be MeV (otherwise GeV) * ************************************************************************ IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE PARAMETER (MAXSDT = 10, & MAXSBI = 1400) COMMON /CSECT/ XSEAV(MAXSDT,MAXSBI),XS(MAXSDT,MAXSBI), & NXSBIN(MAXSDT),NXSDT DIMENSION TMPXS(2,MAXSBI) LOGICAL LFIRST DATA LFIRST /.TRUE./ IF (LFIRST) THEN NXSDT = 0 DO 1 IDT=1,MAXSDT NXSBIN(IDT) = 0 DO 2 IBIN=1,MAXSBI XSEAV(IDT,IBIN) = 0.0D0 XS(IDT,IBIN) = 0.0D0 2 CONTINUE 1 CONTINUE LFIRST = .FALSE. ENDIF NXSDT = NXSDT+1 NEBIN = 0 10 CONTINUE NEBIN = NEBIN+1 IF (NEBIN.GT.MAXSBI) STOP ' nebin > maxsbi !!!' IF (IABS(MODE).EQ.1) THEN READ(LINP,*,END=9000) ELOW,XS(NXSDT,NEBIN) READ(LINP,*) ELHI,XS(NXSDT,NEBIN) XSEAV(NXSDT,NEBIN) = SQRT(ELOW*ELHI) IF (MODE.LT.0) XSEAV(NXSDT,NEBIN) = 1.D-3*XSEAV(NXSDT,NEBIN) ELSEIF (IABS(MODE).EQ.2) THEN READ(LINP,*,END=9000) XSEAV(NXSDT,NEBIN),XS(NXSDT,NEBIN) IF (MODE.LT.0) XSEAV(NXSDT,NEBIN) = 1.D-3*XSEAV(NXSDT,NEBIN) ELSEIF (IABS(MODE).EQ.3) THEN READ(LINP,*,END=9000) TMPXS(1,NEBIN),TMPXS(2,NEBIN) IF (MODE.LT.0) TMPXS(1,NEBIN) = 1.D-3*TMPXS(1,NEBIN) ELSE STOP ' RDXSC: unsupported mode ! ' ENDIF GOTO 10 9000 CONTINUE IF (IABS(MODE).EQ.3) THEN DO 3 I=1,NEBIN-1 IDX = NEBIN-I XSEAV(NXSDT,IDX) = TMPXS(1,I) XS(NXSDT,IDX) = TMPXS(2,I) 3 CONTINUE ENDIF NXSBIN(NXSDT) = NEBIN-1 IDXSEC = NXSDT IF (LCHK.GT.0) THEN AELO = LOG10(XSEAV(IDXSEC,1)) AEHI = LOG10(XSEAV(IDXSEC,NXSBIN(IDXSEC))) DAE = (AEHI-AELO)/DBLE(NXSBIN(IDXSEC)) DO 4 I=1,NXSBIN(IDXSEC)+1 AE = AELO+DBLE(I-1)*DAE E = 10.0D0**AE WRITE(LCHK,'(1X,2E15.5)') E,XSECT(IDXSEC,E,0) 4 CONTINUE ENDIF RETURN END * *===xsect==============================================================* * DOUBLE PRECISION FUNCTION XSECT(IDXSEC,E,MODE) ************************************************************************ * IDXSEC index of stored cross section in common * * E energy * ************************************************************************ IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE PARAMETER (MAXSDT = 10, & MAXSBI = 1400) COMMON /CSECT/ XSEAV(MAXSDT,MAXSBI),XS(MAXSDT,MAXSBI), & NXSBIN(MAXSDT),NXSDT IF (E.LE.XSEAV(IDXSEC,1)) THEN CS = 0.0D0 ELSEIF (E.GT.XSEAV(IDXSEC,NXSBIN(IDXSEC))) THEN N = NXSBIN(IDXSEC)-1 CS = XSCINT(IDXSEC,E,N) ELSE DO 1 J=1,NXSBIN(IDXSEC)-1 IF ((E.GT.XSEAV(IDXSEC,J)).AND.(E.LE.XSEAV(IDXSEC,J+1))) & THEN CS = XSCINT(IDXSEC,E,J) GOTO 2 ENDIF 1 CONTINUE STOP ' xsection value not found ' 2 CONTINUE ENDIF XSECT = CS RETURN END * *===xscint=============================================================* * DOUBLE PRECISION FUNCTION XSCINT(IDXSEC,E,J) IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE PARAMETER (MAXSDT = 10, & MAXSBI = 1400) COMMON /CSECT/ XSEAV(MAXSDT,MAXSBI),XS(MAXSDT,MAXSBI), & NXSBIN(MAXSDT),NXSDT FAC = (LOG10(E) -LOG10(XSEAV(IDXSEC,J)))/ & (LOG10(XSEAV(IDXSEC,J+1))-LOG10(XSEAV(IDXSEC,J))) XSCINT = LOG10(XS(IDXSEC,J)) & +(LOG10(XS(IDXSEC,J+1))-LOG10(XS(IDXSEC,J)))*FAC XSCINT = 10**(XSCINT) RETURN END * *===skip===============================================================* * SUBROUTINE SKIP(LUNIT,NSKIP) IMPLICIT DOUBLE PRECISION (A-H,O-Z) SAVE PARAMETER (LOUT=6,LLOOK=9) CHARACTER*132 ACARD IF (NSKIP.EQ.0) RETURN DO 1 K=1,NSKIP READ(LUNIT,'(A132)') ACARD 1 CONTINUE RETURN END