*$ CREATE SOURCE.FOR
*COPY SOURCE
*
* M. Kirk: copied from source_cpt.f
*
*=== 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.            *
*                                                                      *
*     version 2005 modified 24.11.2005 by FRANK GUTERMUTH              * 
*     transferred to version 2006 by TORSTEN RADON                     *
*     modified 11.06.2007 by MARKUS KIRK                               *
*----------------------------------------------------------------------*
*
      INCLUDE '(BEAMCM)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(IOIOCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(SOURCM)'
      INCLUDE '(SUMCOU)'
*
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* FG: Use option card SOURCE in the input file with the following WHATs
*
* WHAT(1) number of I/O unit to which the formatted data file with the
* source spectrum has been connected (use option card OPEN with
* positive WHAT(1) (formatted file ) and SDUM READONLY to do this)
*
* WHAT(2) number of 'spectral particles' per primary particle
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*
* FG: Structure of data file has to be the following:
*
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* FG: Variable description
*
* ENNUM      = number of energy bins of spectrum
*
* ENDIS(1,J) = lower energy of bin number J of spectrum
* ENDIS(2,J) = number of particles in a bin
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*
*
*
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* FG: Declaration of additional variables
*
      LOGICAL LENERG
*                                You may edit this to suit your needs.
*     was all REAL before...
**      DOUBLE PRECISION ENPRB, ENERG, ENDIS(2,200)
**      INTEGER ENNUM
**      DOUBLE PRECISION WEIGHT
**      DOUBLE PRECISION MAXEN, MINEN, MAXEPB
*
      DOUBLE PRECISION ENPRB, ENERG, ENDIS(2,200)
      INTEGER ENNUM
      DOUBLE PRECISION WEIGHT
      DOUBLE PRECISION MAXEN, MINEN, MAXEPB
*
      SAVE ENPRB, ENERG, ENDIS
      SAVE ENNUM
      SAVE WEIGHT
      SAVE MAXEN, MINEN, MAXEPB
*
      PARAMETER ( TWOPI = TWOTWO * 3.141592653589793238D+00 )
*
* Copied from fluka-discussion forum...
**       DOUBLE PRECISION tkeffff(MXKK),cxtrkkk(MXKK),cytrkkk(MXKK)
**       DATA tkeffff / MXKK*ZERZER /
**       DATA cxtrkkk / MXKK*ZERZER /
**       DATA cytrkkk / MXKK*ZERZER /
**       SAVE tkeffff, cxtrkkk, cytrkkk 
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*
*
      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 ***
* FG: Notify that external SOURCE is used during first particle run 
      WRITE(UNIT=LUNOUT,FMT=*)
     $ 'USER-ROUTINE SOURCE last modified Markus Kirk 11.06.2007'
      WRITE(UNIT=LUNOUT,FMT=*)
     $ 'IJBEAM = ', IJBEAM, ' AND LRDBEA = ', LRDBEA
*
* FG: read data from file SDSOU from directory
*
       WRITE(UNIT=LUNOUT,FMT=*)
     $ 'read data from file SDSOU from directory...'
**       REWIND(WHASOU(1))
       READ(UNIT=WHASOU(1), FMT=*) ENNUM
       DO 50 J=1, ENNUM
*                                    E           N(E)
         READ(UNIT=WHASOU(1), FMT=*) ENDIS(1,J), ENDIS(2,J)
         WRITE(UNIT=LUNOUT,FMT=*) ENDIS(1,J), ENDIS(2,J)
  50   CONTINUE
       WRITE(UNIT=LUNOUT,FMT=*) '...done.'
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      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 I.e. it IS set to zero the first time START is called.
*
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*
* Reading the spectrum from a file was done here. See .bak1 version.
*
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* FG: Find max and min of energetic prob. distribut.
*
        MAXEN=ENDIS(1,ENNUM)
        MINEN=ENDIS(1,1)
        MAXEPB=ENDIS(2,1)
*
        DO 90 J=1, ENNUM
          MAXEPB=MAX(ENDIS(2,J), MAXEPB)
  90    CONTINUE
        WRITE(UNIT=LUNOUT,FMT=*) 'MAXEPB = ', MAXEPB
*
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* FG: Finding particle angle and energy using the 'rejection method'
*
      WRITE(UNIT=LUNOUT,FMT=*) 'WHASOU(2) = ', WHASOU(2)
      DO 300 K=1, WHASOU(2)
      LENERG = .FALSE.
 130  ENERG=MINEN+((MAXEN-MINEN)*FLRNDM(DUMMY))
      ENPRB=MAXEPB*FLRNDM(DUMMY)
      DO 150 J=1, ENNUM-1
        IF ((ENDIS(1,J) .LE. ENERG) .AND.
     $  (ENERG .LT. ENDIS(1,J+1))) THEN
        IF(ENPRB .LE. ENDIS(2,J)) THEN
        LENERG = .TRUE.
        ENDIF
        ENDIF
 150    CONTINUE
        IF (.NOT. LENERG) GO TO 130
        WRITE(UNIT=LUNOUT,FMT=*) 'ENERG = ', ENERG
*
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* FG: Updating stack counter
      NPFLKA = NPFLKA + 1
*
* FG: Set particle weight to the integral of the given distribution
      WEIGHT=0.
      DO 170 I=1, ENNUM
        WEIGHT=WEIGHT+ENDIS(2,I)
 170  CONTINUE
      WRITE(UNIT=LUNOUT,FMT=*) 'WEIGHT = ', WEIGHT
      WTFLK (NPFLKA) = ONEONE
      WEIPRI = WEIPRI + WTFLK (NPFLKA)
      WRITE(UNIT=LUNOUT,FMT=*) 'WTFLK (NPFLKA) = ', WTFLK (NPFLKA)
      WRITE(UNIT=LUNOUT,FMT=*) 'WEIPRI = ', WEIPRI
*
*  +-------------------------------------------------------------------*
*  |  (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
         WRITE(UNIT=LUNOUT,FMT=*) 'IONID  = IJHION = ', IONID
         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 180 ISPR = 1, MKBMX1
         SPAREK (ISPR,NPFLKA) = ZERZER
 180  CONTINUE
* User dependent spare flags:
      DO 190 ISPR = 1, MKBMX2
         ISPARK (ISPR,NPFLKA) = 0
 190  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
*
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* FG: Kinetic energy / momentum / angle of the particle
*
      TKEFLK(NPFLKA) = ENERG
      PMOFLK(NPFLKA) = SQRT(ENERG*(ENERG+TWOTWO*AM(ILOFLK(NPFLKA))))
      TXFLK(NPFLKA)  = UBEAM
      TYFLK(NPFLKA)  = VBEAM
      TZFLK(NPFLKA)  = WBEAM
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*
* Polarization cosines:
      TXPOL  (NPFLKA) = -TWOTWO
      TYPOL  (NPFLKA) = +ZERZER
      TZPOL  (NPFLKA) = +ZERZER
* Particle coordinates
      XFLK   (NPFLKA) = XBEAM
      YFLK   (NPFLKA) = YBEAM
      ZFLK   (NPFLKA) = ZBEAM
*  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 )
*
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
* FG: Load next particle onto the particle stack
*
 300  CONTINUE
* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
*
*  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

