*$ CREATE SOURCE.FOR
*COPY SOURCE
*
*=== source ===========================================================*
*
      SUBROUTINE SOURCE ( NOMORE )

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1990-2010      by    Alfredo Ferrari & Paola Sala  *
*     All Rights Reserved.                                             *
*                                                                      *
*                                                                      *
*     New source for FLUKA9x-FLUKA20xy:                                *
*                                                                      *
*     Created on 07 January 1990   by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*     Last change on  17-Oct-10    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.            *
*                                                                      *
*       Output variables:                                              *
*                                                                      *
*              Nomore = if > 0 the run will be terminated              *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(BEAMCM)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(IOIOCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(SOURCM)'
      INCLUDE '(SUMCOU)'

      LOGICAL LFIRST

*======================================================================*
* NSIZE is a number of energy intervals for generating neutron spectra
* BOUND is an array which contains energy boundaries
* FLUX is an array which contains cumulative fluxes 
      PARAMETER (NSIZE = 254)
      DOUBLE PRECISION BOUND
      DOUBLE PRECISION FLUX
      DIMENSION BOUND(1:NSIZE)
      DIMENSION FLUX(1:NSIZE)
*======================================================================*

      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 ***
      END IF

      DATA BOUND /1.000d-14, 2.000d-13, 2.930d-13, 4.292d-13,
     & 6.287d-13, 9.209d-13, 1.349d-12, 1.976d-12, 2.895d-12,
     & 4.240d-12, 6.212d-12, 9.099d-12, 1.333d-11, 1.953d-11,
     & 2.860d-11, 4.190d-11, 6.138d-11, 8.991d-11, 1.317d-10,
     & 1.929d-10, 2.826d-10, 4.140d-10, 5.316d-10, 6.251d-10,
     & 6.826d-10, 8.337d-10, 8.764d-10, 1.125d-09, 1.445d-09,
     & 1.855d-09, 2.382d-09, 3.059d-09, 3.928d-09, 5.043d-09,
     & 6.476d-09, 8.315d-09, 1.068d-08, 1.371d-08, 1.760d-08,
     & 2.260d-08, 2.902d-08, 3.727d-08, 4.785d-08, 6.144d-08,
     & 7.889d-08, 1.013d-07, 1.301d-07, 1.670d-07, 2.145d-07,
     & 2.754d-07, 3.536d-07, 4.540d-07, 5.829d-07, 7.485d-07,
     & 9.611d-07, 1.234d-06, 1.364d-06, 1.507d-06, 1.585d-06,
     & 1.666d-06, 1.841d-06, 2.035d-06, 2.249d-06, 2.485d-06,
     & 2.613d-06, 2.747d-06, 2.863d-06, 3.035d-06, 3.355d-06,
     & 3.707d-06, 4.097d-06, 4.307d-06, 4.528d-06, 5.005d-06,
     & 5.531d-06, 6.267d-06, 7.102d-06, 8.047d-06, 9.119d-06,
     & 1.033d-05, 1.171d-05, 1.327d-05, 1.503d-05, 1.704d-05,
     & 1.930d-05, 2.133d-05, 2.187d-05, 2.358d-05, 2.418d-05,
     & 2.479d-05, 2.606d-05, 2.809d-05, 3.183d-05, 3.431d-05,
     & 3.518d-05, 3.607d-05, 4.087d-05, 4.631d-05, 5.248d-05,
     & 5.656d-05, 5.946d-05, 6.251d-05, 6.738d-05, 7.635d-05,
     & 8.652d-05, 9.804d-05, 1.111d-04, 1.168d-04, 1.228d-04,
     & 1.291d-04, 1.357d-04, 1.426d-04, 1.500d-04, 1.576d-04,
     & 1.616d-04, 1.657d-04, 1.699d-04, 1.742d-04, 1.832d-04,
     & 1.925d-04, 2.024d-04, 2.128d-04, 2.237d-04, 2.352d-04,
     & 2.472d-04, 2.599d-04, 2.732d-04, 2.802d-04, 2.872d-04,
     & 2.945d-04, 3.020d-04, 3.096d-04, 3.175d-04, 3.337d-04,
     & 3.508d-04, 3.688d-04, 3.877d-04, 4.076d-04, 4.285d-04,
     & 4.505d-04, 4.736d-04, 4.979d-04, 5.105d-04, 5.234d-04,
     & 5.366d-04, 5.502d-04, 5.784d-04, 6.081d-04, 6.393d-04,
     & 6.721d-04, 7.065d-04, 7.427d-04, 7.808d-04, 8.209d-04,
     & 8.629d-04, 9.072d-04, 9.537d-04, 9.616d-04, 9.778d-04,
     & 1.003d-03, 1.054d-03, 1.108d-03, 1.165d-03, 1.194d-03,
     & 1.225d-03, 1.287d-03, 1.353d-03, 1.423d-03, 1.496d-03,
     & 1.534d-03, 1.572d-03, 1.612d-03, 1.653d-03, 1.738d-03,
     & 1.827d-03, 1.873d-03, 1.920d-03, 1.969d-03, 2.019d-03,
     & 2.122d-03, 2.231d-03, 2.269d-03, 2.307d-03, 2.346d-03,
     & 2.365d-03, 2.385d-03, 2.425d-03, 2.466d-03, 2.592d-03,
     & 2.725d-03, 2.865d-03, 3.012d-03, 3.088d-03, 3.166d-03,
     & 3.247d-03, 3.329d-03, 3.499d-03, 3.679d-03, 3.867d-03,
     & 4.066d-03, 4.274d-03, 4.493d-03, 4.607d-03, 4.724d-03,
     & 4.843d-03, 4.966d-03, 5.220d-03, 5.488d-03, 5.770d-03,
     & 5.916d-03, 6.065d-03, 6.219d-03, 6.376d-03, 6.538d-03,
     & 6.592d-03, 6.648d-03, 6.703d-03, 6.873d-03, 7.047d-03,
     & 7.225d-03, 7.408d-03, 7.596d-03, 7.788d-03, 7.985d-03,
     & 8.187d-03, 8.395d-03, 8.607d-03, 8.825d-03, 9.048d-03,
     & 9.277d-03, 9.512d-03, 9.753d-03, 1.000d-02, 1.025d-02,
     & 1.051d-02, 1.078d-02, 1.105d-02, 1.133d-02, 1.162d-02,
     & 1.191d-02, 1.221d-02, 1.252d-02, 1.284d-02, 1.317d-02,
     & 1.350d-02, 1.384d-02, 1.419d-02, 1.455d-02, 1.492d-02,
     & 1.530d-02, 1.568d-02, 1.608d-02, 1.649d-02, 1.690d-02/

      DATA FLUX /0.000000d+00, 2.763974d-10, 3.024944d-08,
     & 1.274075d-07, 7.721650d-07, 2.896046d-06, 1.164409d-05,
     & 3.609419d-05, 1.133372d-04, 3.291885d-04, 8.622172d-04,
     & 2.084253d-03, 4.881060d-03, 1.063404d-02, 2.143630d-02,
     & 3.936589d-02, 6.361617d-02, 8.808168d-02, 1.057764d-01,
     & 1.159331d-01, 1.227520d-01, 1.289157d-01, 1.330855d-01,
     & 1.358113d-01, 1.372625d-01, 1.406727d-01, 1.415266d-01,
     & 1.457864d-01, 1.500435d-01, 1.543078d-01, 1.586179d-01,
     & 1.630235d-01, 1.674072d-01, 1.717942d-01, 1.762092d-01,
     & 1.806820d-01, 1.852370d-01, 1.898078d-01, 1.944402d-01,
     & 1.991852d-01, 2.038825d-01, 2.087344d-01, 2.136129d-01,
     & 2.185466d-01, 2.234902d-01, 2.285560d-01, 2.336089d-01,
     & 2.386977d-01, 2.439657d-01, 2.492125d-01, 2.544960d-01,
     & 2.598240d-01, 2.652464d-01, 2.707554d-01, 2.762206d-01,
     & 2.817859d-01, 2.840920d-01, 2.864237d-01, 2.875669d-01,
     & 2.886813d-01, 2.909963d-01, 2.933245d-01, 2.957008d-01,
     & 2.980939d-01, 2.992674d-01, 3.004606d-01, 3.014589d-01,
     & 3.028680d-01, 3.053674d-01, 3.078542d-01, 3.103743d-01,
     & 3.116212d-01, 3.129099d-01, 3.155237d-01, 3.181473d-01,
     & 3.215044d-01, 3.248988d-01, 3.284565d-01, 3.320075d-01,
     & 3.356997d-01, 3.394483d-01, 3.433703d-01, 3.473944d-01,
     & 3.516203d-01, 3.559950d-01, 3.596463d-01, 3.606101d-01,
     & 3.634580d-01, 3.644533d-01, 3.654133d-01, 3.673756d-01,
     & 3.703792d-01, 3.755912d-01, 3.788391d-01, 3.799490d-01,
     & 3.810540d-01, 3.868167d-01, 3.929802d-01, 3.994367d-01,
     & 4.036644d-01, 4.065120d-01, 4.093689d-01, 4.138443d-01,
     & 4.215444d-01, 4.299451d-01, 4.389332d-01, 4.487399d-01,
     & 4.528877d-01, 4.570903d-01, 4.614329d-01, 4.659518d-01,
     & 4.705113d-01, 4.753742d-01, 4.804055d-01, 4.829820d-01,
     & 4.855182d-01, 4.881542d-01, 4.908271d-01, 4.962377d-01,
     & 5.018885d-01, 5.078612d-01, 5.139489d-01, 5.202551d-01,
     & 5.267474d-01, 5.335440d-01, 5.405140d-01, 5.478104d-01,
     & 5.515167d-01, 5.552277d-01, 5.590226d-01, 5.628922d-01,
     & 5.668279d-01, 5.708795d-01, 5.790029d-01, 5.872481d-01,
     & 5.954402d-01, 6.034829d-01, 6.111587d-01, 6.167375d-01,
     & 6.205586d-01, 6.278544d-01, 6.377849d-01, 6.433680d-01,
     & 6.492103d-01, 6.551967d-01, 6.614707d-01, 6.743845d-01,
     & 6.882024d-01, 7.021523d-01, 7.153714d-01, 7.288379d-01,
     & 7.428803d-01, 7.571754d-01, 7.716760d-01, 7.860657d-01,
     & 7.994571d-01, 8.116445d-01, 8.133637d-01, 8.165945d-01,
     & 8.203192d-01, 8.278732d-01, 8.384109d-01, 8.496698d-01,
     & 8.557990d-01, 8.621782d-01, 8.743563d-01, 8.823092d-01,
     & 8.921957d-01, 9.021439d-01, 9.069365d-01, 9.115568d-01,
     & 9.157733d-01, 9.196909d-01, 9.280575d-01, 9.356690d-01,
     & 9.389622d-01, 9.420476d-01, 9.452690d-01, 9.485372d-01,
     & 9.546892d-01, 9.602348d-01, 9.620149d-01, 9.637930d-01,
     & 9.656544d-01, 9.665662d-01, 9.674717d-01, 9.691090d-01,
     & 9.706249d-01, 9.746630d-01, 9.781464d-01, 9.811428d-01,
     & 9.837588d-01, 9.848672d-01, 9.859324d-01, 9.868093d-01,
     & 9.875702d-01, 9.889074d-01, 9.900179d-01, 9.910270d-01,
     & 9.919972d-01, 9.928967d-01, 9.936570d-01, 9.940383d-01,
     & 9.944312d-01, 9.948062d-01, 9.951554d-01, 9.957701d-01,
     & 9.963755d-01, 9.968748d-01, 9.971175d-01, 9.973375d-01,
     & 9.975541d-01, 9.977368d-01, 9.979431d-01, 9.979889d-01,
     & 9.980582d-01, 9.981159d-01, 9.982940d-01, 9.984768d-01,
     & 9.986176d-01, 9.987685d-01, 9.988765d-01, 9.989870d-01,
     & 9.990861d-01, 9.991979d-01, 9.992767d-01, 9.993642d-01,
     & 9.994369d-01, 9.995084d-01, 9.995671d-01, 9.996314d-01,
     & 9.996844d-01, 9.997344d-01, 9.997668d-01, 9.998030d-01,
     & 9.998370d-01, 9.998669d-01, 9.998863d-01, 9.999011d-01,
     & 9.999147d-01, 9.999285d-01, 9.999437d-01, 9.999542d-01,
     & 9.999609d-01, 9.999667d-01, 9.999739d-01, 9.999798d-01,
     & 9.999843d-01, 9.999915d-01, 9.999929d-01, 9.999954d-01,
     & 9.999981d-01, 9.999987d-01, 1.000000d+00/

*  |
*  +-------------------------------------------------------------------*
*  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.
*  |  Group number for "low" energy neutrons, set to 0 anyway
         IGROUP (NPFLKA) = 0
*  |
*  +-------------------------------------------------------------------*
*  |  Normal hadron:
      ELSE
         IONID = IJBEAM
         ILOFLK (NPFLKA) = IJBEAM
*  |  Flag this is prompt radiation
         LRADDC (NPFLKA) = .FALSE.
*  |  Group number for "low" energy neutrons, set to 0 anyway
         IGROUP (NPFLKA) = 0
      END IF
*  |
*  +-------------------------------------------------------------------*
*  From this point .....
*  Particle generation (1 for primaries)
      LOFLK  (NPFLKA) = 1
*  User dependent flag:
      LOUSE  (NPFLKA) = 0
*  No channeling:
      LCHFLK (NPFLKA) = .FALSE.
      DCHFLK (NPFLKA) = ZERZER
*  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
*  Kinetic energy of the particle (GeV)

*======================================================================*
* kinetic energy sampling 
      PART=FLRNDM(X)
      DO i=1, NSIZE-1
         IF (PART.GE.FLUX(i).AND.PART.LT.FLUX(i+1)) THEN
            TKEFLK (NPFLKA) = BOUND(i)+(BOUND(i+1)-BOUND(i))*
     &(PART-FLUX(i))/(FLUX(i+1)-FLUX(i))
            EXIT
         END IF
      END DO
*======================================================================*

*  Particle momentum
      PMOFLK (NPFLKA) = PBEAM
*     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

*======================================================================*
* cylindrical beam sampling
      RRR = (2.49D+01*FLRNDM(R))**0.5D+00
      AAA = FLRNDM(A) * TWOPIP
      XBEAM = RRR * DCOS(AAA)
      YBEAM = RRR * DSIN(AAA)
      ZBEAM = -1.0D+00
      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 )
*  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

