OK, in response to the replies, I am posting my source file. I
wouldn't understand however, how it would run on my associate's
computer but not on mine. To check, I had her email me the source
file and input file, and it still gave me the same error. Thanks for
looking,
Ariel
*$ CREATE SOURCE.FOR
*COPY SOURCE
*
*=== source ===========================================================*
*
SUBROUTINE SOURCE ( NOMORE )
INCLUDE '(DBLPRC)'
INCLUDE '(DIMPAR)'
INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
* *
* Copyright (C) 1990-2006 by Alfredo Ferrari & Paola Sala *
* All Rights Reserved. *
* *
* *
* New source for FLUKA9x-FLUKA200x: *
* *
* Created on 07 january 1990 by Alfredo Ferrari & Paola Sala *
* Infn - Milan *
* *
* Last change on 03-mar-06 by Alfredo Ferrari *
* *
* This is just an example of a possible user written source routine. *
* note that the beam card still has some meaning - in the scoring the *
* maximum momentum used in deciding the binning is taken from the *
* beam momentum. Other beam card parameters are obsolete. *
* *
*----------------------------------------------------------------------*
*
INCLUDE '(BEAMCM)'
INCLUDE '(FHEAVY)'
INCLUDE '(FLKSTK)'
INCLUDE '(IOIOCM)'
INCLUDE '(LTCLCM)'
INCLUDE '(PAPROP)'
INCLUDE '(SOURCM)'
INCLUDE '(SUMCOU)'
*
LOGICAL LFIRST
*
SAVE LFIRST
DATA LFIRST / .TRUE. /
PARAMETER (COEFC = 1.230E+12)
PARAMETER (COEFK = 2.115E+00)
PARAMETER (ALPHA = 0.2815E+00)
PARAMETER (IBNUM = 999) ! number of energy bins
DIMENSION ENERGY(0:IBNUM) ! Initial Energy Of Each Bin
DIMENSION SPECTR(0:IBNUM) ! differential fluence function
DIMENSION CUMUL(0:IBNUM) ! Cumulative function
*======================================================================*
* *
* BASIC VERSION *
* *
*======================================================================*
NOMORE = 0
* +-------------------------------------------------------------------*
* | First call initializations:
IF ( LFIRST ) THEN
* | *** The following 3 cards are mandatory ***
TKESUM = ZERZER
LFIRST = .FALSE.
LUSSRC = .TRUE.
* | *** User initialization ***
WRITE(*,*)
WRITE(*,*) "source foror Al slab/ monodirectional beam"
WRITE(*,*) "Version 1 called"
WRITE(*,*)
DO I = 0, IBNUM
! WRITE(*,*) I
ENERGY(I) = 1.0E-3 + 1.0E-3 * I
! WRITE(*,*) I, ENERGY(I)
SPECTR(I) = COEFC * COEFK * ALPHA *
& EXP(-COEFK * (ENERGY(I) ** ALPHA))*
& ENERGY(I) ** (ALPHA - 1.D0)
WRITE(*,*) I, ENERGY(I), SPECTR(I)
END DO
! WRITE(*,*) "enegry, spectr initialized"
CUMUL(0) = ZERZER
* calculate the cumulative function
DO I = 1, IBNUM
CUMUL(I) = CUMUL(I-1) + HLFHLF * (SPECTR(I) + SPECTR(I-1))
& * (ENERGY(I) - ENERGY(I-1))
WRITE(*,*) I, ENERGY(I), CUMUL(I)
END DO
! WRITE(*,*) "cumulative function initialized", CUMUL(IBNUM)
* normalize the cumulative function to 1
SUM = CUMUL(IBNUM)
DO I = 1, IBNUM
CUMUL(I) = CUMUL(I)/SUM
WRITE(*,*) I, ENERGY(I), CUMUL(I)
END DO
! WRITE(*,*) "cumulative function norm", CUMUL(IBNUM)
END IF
* IJBEAM = 1
* |
* +-------------------------------------------------------------------*
* Push one source particle to the stack. Note that you could as well
* push many but this way we reserve a maximum amount of space in the
* stack for the secondaries to be generated
* Npflka is the stack counter: of course any time source is called it
* must be =0
NPFLKA = NPFLKA + 1
* Wt is the weight of the particle
WTFLK (NPFLKA) = ONEONE
WEIPRI = WEIPRI + WTFLK (NPFLKA)
* Particle type (1=proton.....). Ijbeam is the type set by the BEAM
* card
* +-------------------------------------------------------------------*
* | (Radioactive) isotope:
IF ( IJBEAM .EQ. -2 .AND. LRDBEA ) THEN
IARES = IPROA
IZRES = IPROZ
IISRES = IPROM
CALL STISBM ( IARES, IZRES, IISRES )
IJHION = IPROZ * 1000 + IPROA
IJHION = IJHION * 100 + KXHEAV
IONID = IJHION
CALL DCDION ( IONID )
CALL SETION ( IONID )
* |
* +-------------------------------------------------------------------*
* | Heavy ion:
ELSE IF ( IJBEAM .EQ. -2 ) THEN
IJHION = IPROZ * 1000 + IPROA
IJHION = IJHION * 100 + KXHEAV
IONID = IJHION
CALL DCDION ( IONID )
CALL SETION ( IONID )
ILOFLK (NPFLKA) = IJHION
* | Flag this is prompt radiation
LRADDC (NPFLKA) = .FALSE.
* |
* +-------------------------------------------------------------------*
* | Normal hadron:
ELSE
IONID = IJBEAM
ILOFLK (NPFLKA) = IJBEAM
* | Flag this is prompt radiation
LRADDC (NPFLKA) = .FALSE.
END IF
* |
* +-------------------------------------------------------------------*
* From this point .....
* Particle generation (1 for primaries)
LOFLK (NPFLKA) = 1
* User dependent flag:
LOUSE (NPFLKA) = 0
* User dependent spare variables:
DO 100 ISPR = 1, MKBMX1
SPAREK (ISPR,NPFLKA) = ZERZER
100 CONTINUE
* User dependent spare flags:
DO 200 ISPR = 1, MKBMX2
ISPARK (ISPR,NPFLKA) = 0
200 CONTINUE
* Save the track number of the stack particle:
ISPARK (MKBMX2,NPFLKA) = NPFLKA
NPARMA = NPARMA + 1
NUMPAR (NPFLKA) = NPARMA
NEVENT (NPFLKA) = 0
DFNEAR (NPFLKA) = +ZERZER
* ... to this point: don't change anything
* Particle age (s)
AGESTK (NPFLKA) = +ZERZER
AKNSHR (NPFLKA) = -TWOTWO
* Group number for "low" energy neutrons, set to 0 anyway
IGROUP (NPFLKA) = 0
* Kinetic energy of the particle (GeV)
! TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID)
!
! Sampling from the distribution
!
XI = FLRNDM(DUMMY)
WRITE(*,*) " XI", XI
DO 111 I = 0, IBNUM - 1
IF ( XI .LT. CUMUL(I+1) ) THEN
* come here when cumul(i) =< xi < cumul(i+1). Interpolate
write(*,*) ENERGY(I), ENERGY(I+1)
SAMPLE = ENERGY(I) + (ENERGY(I+1) - ENERGY(I))
& * (XI - CUMUL(I)) / (CUMUL(I+1) - CUMUL(I))
WRITE(*, *) "Sample", SAMPLE
GO TO 222
END IF
111 CONTINUE
222 CONTINUE
TKEFLK (NPFLKA) = SAMPLE
!
!
* Particle momentum
* PMOFLK (NPFLKA) = PBEAM
!
PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
& + TWOTWO * AM (IONID) ) )
WRITE(*, *) " momentum", PMOFLK (NPFLKA)
!
* Cosines (tx,ty,tz)
! TXFLK (NPFLKA) = UBEAM
! TYFLK (NPFLKA) = VBEAM
! TZFLK (NPFLKA) = WBEAM
TXFLK (NPFLKA) = 0.D0
TYFLK (NPFLKA) = 0.D0
TZFLK (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
& - TYFLK (NPFLKA)**2 )
WRITE(*, *) " cos z", TZFLK (NPFLKA)
!
* Polarization cosines:
TXPOL (NPFLKA) = -TWOTWO
TYPOL (NPFLKA) = +ZERZER
TZPOL (NPFLKA) = +ZERZER
* Particle coordinates
! XFLK (NPFLKA) = XBEAM
! YFLK (NPFLKA) = YBEAM
! ZFLK (NPFLKA) = ZBEAM
XFLK (NPFLKA) = -40.D0 + 80.D0 * FLRNDM(DUMMY)
YFLK (NPFLKA) = -40.D0 + 80.D0 * FLRNDM(DUMMY)
ZFLK (NPFLKA) = -10.D0
WRITE(*, *) " cordinates",XFLK(NPFLKA),YFLK(NPFLKA),ZFLK (NPFLKA)
WRITE(*, *) "***************************************************"
* Calculate the total kinetic energy of the primaries: don't change
IF ( ILOFLK (NPFLKA) .EQ. -2 .OR. ILOFLK (NPFLKA) .GT. 100000 )
& THEN
TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
ELSE IF ( ILOFLK (NPFLKA) .NE. 0 ) THEN
TKESUM = TKESUM + ( TKEFLK (NPFLKA) + AMDISC (ILOFLK(NPFLKA)) )
& * WTFLK (NPFLKA)
ELSE
TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
END IF
WRITE(*, *) " total kinetic energy", TKESUM
RADDLY (NPFLKA) = ZERZER
* Here we ask for the region number of the hitting point.
* NREG (NPFLKA) = ...
* The following line makes the starting region search much more
* robust if particles are starting very close to a boundary:
CALL GEOCRS ( TXFLK (NPFLKA), TYFLK (NPFLKA), TZFLK (NPFLKA) )
CALL GEOREG ( XFLK (NPFLKA), YFLK (NPFLKA), ZFLK (NPFLKA),
& NRGFLK(NPFLKA), IDISC )
* Do not change these cards:
CALL GEOHSM ( NHSPNT (NPFLKA), 1, -11, MLATTC )
NLATTC (NPFLKA) = MLATTC
CMPATH (NPFLKA) = ZERZER
CALL SOEVSV
* WRITE(*, *) " just before exiting tracker"
RETURN
*=== End of subroutine Source =========================================*
END
Received on Fri Jul 11 17:40:35 2008
This archive was generated by hypermail 2.1.8 : Fri Jul 11 2008 - 17:40:35 CEST