Re: source.f

Dear Francesca
the write statements are correct, but your source routine is not well
written and produces a lot of errors.
In order to help you in the right way we need also your input file.
In any case, in attachment I put your source modified in a few points.
If you a make a diff you will see immediately the relevant points.
Using a generic input I succeeded in running the
job and having unit 80 filled (at the end of a run
it will be called xxxx00n_data.dat, where xxxx is the input name and n
is the cycle numnber).
In my opinion the most serious problem was that you were sampling the
energy but you did not fill the
momentum variable PMOFLK.
I am not sure that you can understand completely all the points, and I
am not sure to understand exactly which was you goal, therefore I suggest
that you contact me in Milano so I can explain to you all details.
    Best regards
       Giuseppe

P.S.: remember to subscribe to fluka-discuss and send mails to the list
without html

Botta Francesca wrote:
> dear fluka users,
>
> I'm just a beginner with the use of user routines, in particular source.f.
> for debug purposes I was trying to put in my source-1cm.f file
> (attached) a "write" command to print the coordinates and the energy
> of each primary particle on a file, but it doesn't work
>
> could you please tell me if it's possible to do it and how?
>
> thanks so much in advance
>
> all the best
>
> francesca

*$ 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)'
      DIMENSION F(21),K(21),ES(21),PTS(21)
*
      LOGICAL LFIRST
*
      SAVE LFIRST
      DATA LFIRST / .TRUE. /
*======================================================================*
* *
* BASIC VERSION *
* *
*======================================================================*
      NOMORE = 0
* +-------------------------------------------------------------------*
* | First call initializations:
      IF ( LFIRST ) THEN
* | *** The following 3 cards are mandatory ***
         TKESUM = ZERZER
         LFIRST = .FALSE.
         LUSSRC = .TRUE.
* | *** User initialization ***
* pigreco=3.141593
*
* sphere radius
*
         R=1.0D+00 ! in cm
*
* open file
*
      OPEN(UNIT=80,FILE='data.dat')
      WRITE(80,*)'coordinates(cm),energy(GeV)'
*
      END IF
* |
* +-------------------------------------------------------------------*
* 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) - uniform in (0-1] MeV
      RN=ONEONE-FLRNDM()
      ESAMP = RN*0.001D+00
      TKEFLK (NPFLKA) = ESAMP
* Particle momentum
      PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
     & + TWOTWO * AM (IONID) ) )
* Cosines (tx,ty,tz)
      TXFLK (NPFLKA) = UBEAM
      TYFLK (NPFLKA) = VBEAM
      TZFLK (NPFLKA) = WBEAM
* TZFLK (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
* & - TYFLK (NPFLKA)**2 )
* Polarization cosines:
      TXPOL (NPFLKA) = -TWOTWO
      TYPOL (NPFLKA) = +ZERZER
      TZPOL (NPFLKA) = +ZERZER
* Particle coordinates - SAMPLED INSIDE SPHERE RADIUS R
      CALL SFLOOD(XSPH,YSPH,ZSPH,UUUX,VVVY,WWWZ)
      XFLK (NPFLKA) = XSPH
      YFLK (NPFLKA) = YSPH
      ZFLK (NPFLKA) = ZSPH
*
* print energy and coordinates on file
*
      write(80,*)ESAMP,XSPH,YSPH,ZSPH
*
* 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
      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
      RETURN
*=== End of subroutine Source =========================================*
      END
Received on Mon Oct 26 2009 - 10:21:20 CET