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

       INCLUDE '(DBLPRC)'
       INCLUDE '(DIMPAR)'
       INCLUDE '(IOUNIT)'
*
       INCLUDE '(BEAMCM)'
       INCLUDE '(FHEAVY)'
       INCLUDE '(FLKSTK)'
       INCLUDE '(IOIOCM)'
       INCLUDE '(LTCLCM)'
       INCLUDE '(PAPROP)'
       INCLUDE '(SOURCM)'
       INCLUDE '(SUMCOU)'
*
       LOGICAL LFIRST
*
       PARAMETER (NCF = 51)
       DIMENSION CF(0:NCF), ENERCF(0:NCF)
       DIMENSION CUMUCF(0:NCF)
       SAVE LFIRST
       DATA LFIRST / .TRUE. /
       DATA ENERCF / 4.14d-10,1.00d-09,1.00d-08,5.00d-08,1.00d-07,
     & 2.00d-07,4.00d-07,7.00d-07,1.00d-06,3.00d-06,6.00d-06,1.00d-05,
     & 2.00d-05,4.00d-05,6.00d-05,8.00d-05,1.00d-04,1.50d-04,2.00d-04,
     & 2.50d-04,3.00d-04,3.50d-04,4.00d-04,4.50d-04,5.00d-04,5.50d-04,
     & 6.00d-04,7.00d-04,8.00d-04,9.00d-04,1.00d-03,1.20d-03,1.40d-03,
     & 1.60d-03,1.80d-03,2.00d-03,2.30d-03,2.60d-03,3.00d-03,3.50d-03,
     & 4.00d-03,4.50d-03,5.00d-03,6.00d-03,7.00d-03,8.00d-03,9.00d-03,
     & 1.00d-02,1.10d-02,1.20d-02,1.30d-02,1.40d-02 /
* spectra are in n/cm2
       DATA CF / 3.10E-10,1.11E-08,1.27E-07,2.76E-07,7.82E-07,2.21E-06,
     & 4.53E-06,5.68E-06,5.51E-05,1.28E-04,2.30E-04,7.74E-04,2.17E-03,
     & 2.80E-03,3.29E-03,3.68E-03,1.05E-02,1.21E-02,1.33E-02,1.42E-02,
     & 1.49E-02,1.55E-02,1.60E-02,1.63E-02,1.66E-02,1.68E-02,3.38E-02,
     & 3.39E-02,3.37E-02,3.33E-02,6.46E-02,6.12E-02,5.73E-02,5.31E-02,
     & 4.88E-02,6.55E-02,5.67E-02,6.33E-02,6.21E-02,4.68E-02,3.49E-02,
     & 2.58E-02,3.30E-02,1.74E-02,9.01E-03,4.61E-03,2.33E-03,1.17E-03,
     & 5.83E-04,2.88E-04,1.42E-04,6.94E-05 /
       NOMORE = 0
*  +-------------------------------------------------------------------*
*  |  First call initializations:
       IF ( LFIRST ) THEN
*  |  *** The following 3 cards are mandatory ***
          TKESUM = ZERZER
          LFIRST = .FALSE.
          LUSSRC = .TRUE.
*  |  *** User initialization ***
          WRITE(LUNOUT,*)
          WRITE(LUNOUT,'(A,132A)') ("*",I=1,132)
          WRITE(LUNOUT,*)
          WRITE(LUNOUT,*)
          cumucf(0) = 0.d0
*           prepare the cumulative function
             do i = 1, ncf
                cumucf(i) = cumucf(i-1) +
     &                     (cf(i-1)+cf(i))*(enercf(i)-enercf(i-1))*0.5d0
             end do
*               WRITE(88,*) "Lethargies, Cumulative, non normalized"
*               DO I = 0, NCF
*                  WRITE(88,'(I3,1X,1P,2G15.9)') I, ENERCF(I), CUMUCF(I)
*               END DO
*           normalize the cumulative function
             cunorm = cumucf(ncf)
             do i = 1, ncf
                cumucf(i) = cumucf(i)/cunorm
             end do
*               WRITE(88,*) "Lethargies, Cumulative, normalized"
*               DO I = 0, NCF
*                  WRITE(88,'(I3,1X,1P,2G15.9)') I, ENERCF(I), CUMUCF(I)
*               END DO
*         kount = 0
*      WRITE(88,*)
*     & "Sampled: xi, cum(i-1), cum(i), leth(i-1), leth(i), leth, energy"
       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
*      kount = kount + 1
*      if(kount .gt. 100) nomore = 1
       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). Sample from the cumulative
       xi = flrndm(dummy)
       if(nint(whasou(1)) .eq. 1) then
          do i = 1, ncf
             if(xi .gt. cumucf(i-1) .and. xi .le. cumucf(i)) then
                enkin = enercf(i-1)+(enercf(i)-enercf(i-1))*
     &         (sqrt(xi) - sqrt(cumucf(i-1)))/
     &                               (sqrt(cumucf(i))-sqrt(cumucf(i-1)))
*               WRITE(88,'(3G12.4,1P,4G12.4)')
*     & XI,CUMUCF(I-1),CUMUCF(I),ENERCF(I-1),ENERCF(I),ENKIN,EXP(ENKIN)
                go to 1
             end if
          end do
       else if(nint(whasou(1)) .eq. 2) then
          do i = 1, nambe
             if(xi .gt. cumuam(i-1) .and. xi .le. cumuam(i)) then
                enkin = eneram(i-1)+(eneram(i)-eneram(i-1))*
     &        (sqrt(xi) - sqrt(cumuam(i-1)))/
     &                               (sqrt(cumuam(i))-sqrt(cumuam(i-1)))
                go to 1
             end if
          end do
       end if
  1    continue
*     back from lethargy to energy
       tkeflk (npflka) = exp(enkin)
*     TKEFLK (NPFLKA) = SQRT ( PBEAM**2 + AM (IONID)**2 ) - AM (IONID)
* 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
* 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 )
*  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

