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

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1990-2020      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  03-Jan-20    by             Paola Sala           *
*                                                Infn - Milan          *
*                                                                      *
*  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)'
*     
      real*8 ParRndm                                  !My Code*********
      integer,parameter::neutron_dimension=1012,gamma_dimension=1001
      real*8 energy_low,
     & energy_high,spec_count,
     & neutron_spec_count(neutron_dimension),
     & gamma_spec_count(gamma_dimension),
     & increment_count
      character(len=1000)::temp,temp1
      integer i
      real*8 neu_spec_pro(1012),gamma_spec_pro(1001)
      real*8 energy
      real*8 DirRndm
      real*8 EneRndm 
      real*8 thi     
*
      LOGICAL LFIRST, LISNUT
      SAVE LFIRST
      DATA LFIRST / .TRUE. /
*  Statement function:
      LISNUT (IJ) = INDEX ( PRNAME (IJ), 'NEUTRI' ) .GT. 0
*======================================================================*
*                                                                      *
*                 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
*  |
*  +-------------------------------------------------------------------*
*  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)  
      neutron_spec_count(1)=0
      gamma_spec_count(1)=0    
      ParRndm=FLRNDM(XDUMMY)                        !My Code*********
      IF(ParRndm.LT.1.076E-3/(1.076E-3+6.91E-6))THEN    !My Code*********
      IJBEAM=8                                      !My Code*********
      ELSE                                          !My Code*********
      IJBEAM=7                                      !My Code*********
      END IF                                        !My Code*********
*  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 = IPROM  * 100000 + MOD ( IPROZ, 100 ) * 1000 + IPROA
         IJHION = IJHION * 100    + KXHEAV
         IONID  = IJHION
         CALL DCDION ( IONID )
         CALL SETION ( IONID )
         EEXSTK (NPFLKA) = EXENRG (IONID)
         TMNSTK (NPFLKA) = TMNLF  (IONID)
         ILVSTK (NPFLKA) = IEXLVL (IONID)
         LFRPHN (NPFLKA) = .FALSE.
*  |
*  +-------------------------------------------------------------------*
*  |  Heavy ion:
      ELSE IF ( IJBEAM .EQ. -2 ) THEN
         IJHION = IPROM  * 100000 + MOD ( IPROZ, 100 ) * 1000 + IPROA
         IJHION = IJHION * 100    + KXHEAV
         IONID  = IJHION
         CALL DCDION ( IONID )
         CALL SETION ( IONID )
         EEXSTK (NPFLKA) = EXENRG (IONID)
         TMNSTK (NPFLKA) = TMNLF  (IONID)
         ILVSTK (NPFLKA) = IEXLVL (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
*  |  Parent radioactive isotope:
         IRDAZM (NPFLKA) = 0
*  |  Particle age (s)
         AGESTK (NPFLKA) = +ZERZER
*  |  Kinetic energy of the particle (GeV)
         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) ) )
         LFRPHN (NPFLKA) = .FALSE.
*  |
*  +-------------------------------------------------------------------*
*  |  Normal hadron:
      ELSE
         IONID = IJBEAM
         EEXSTK (NPFLKA) = EXENRG (IJBEAM)
         TMNSTK (NPFLKA) = TMNLF  (IJBEAM)
         ILVSTK (NPFLKA) = IEXLVL (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
*  |  Parent radioactive isotope:
         IRDAZM (NPFLKA) = 0
*  |  Particle age (s)
         AGESTK (NPFLKA) = +ZERZER
*  |  Kinetic energy of the particle (GeV)
      OPEN(UNIT=11,FILE="/home/llh/flukaworks/np_spectrum/80M/n31")     
      OPEN(UNIT=12,FILE="/home/llh/flukaworks/np_spectrum/80M/n32")     
      OPEN(UNIT=13,FILE="/home/llh/flukaworks/np_spectrum/80M/n33")     
      OPEN(UNIT=14,FILE="/home/llh/flukaworks/np_spectrum/80M/n34")     
      OPEN(UNIT=15,FILE="/home/llh/flukaworks/np_spectrum/80M/n35")     
      OPEN(UNIT=16,FILE="/home/llh/flukaworks/np_spectrum/80M/n36")   
      OPEN(UNIT=17,FILE="/home/llh/flukaworks/np_spectrum/80M/n37")     
      OPEN(UNIT=21,FILE="/home/llh/flukaworks/np_spectrum/80M/p51")     
      OPEN(UNIT=22,FILE="/home/llh/flukaworks/np_spectrum/80M/p52")    
*  Cosines (tx,ty,tz)
******************************************************************
* when the particle is neutron
******************************************************************
      IF(IJBEAM.eq.8)THEN
      DirRndm=FLRNDM(XDUMMY)     
*****************************************************************
*neutron moment angle is 0-2.5 
*****************************************************************
      IF(DirRndm.le.3.093E-5/1.076E-3)THEN
      WBEAM=cos(FLRNDM(XDUMMY)*2.5/180*3.1415926)
      read(11,"(A1000)")temp1
      read(11,"(A1000)")temp1
      do i=1,1011
      read(11,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      read(temp(25:33),*)spec_count
      increment_count=increment_count+spec_count*(energy_high
     & -energy_low)
      neutron_spec_count(i+1)=increment_count
      end do
      do i=1,1012
      neu_spec_pro(i)=neutron_spec_count(i)/neutron_spec_count(1012)
      end do
      rewind(11)
      read(11,"(A1000)")temp1
      read(11,"(A1000)")temp1
      EneRndm=FLRNDM(XDUMMY)
      do i=1,1011     
      read(11,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      if(EneRndm.gt.neu_spec_pro(i).and.EneRndm.le.neu_spec_pro(i+1))
     & then
      energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY)
      end if
      end do
************************************************************************
*the neutron moment angle is 0-2.5
************************************************************************
*************************************************************************
*the neutron moment angle is 2.5-5
************************************************************************
      else if(DirRndm.gt.3.093E-5/1.076E-3.and.
     & DirRndm.le.1.146e-4/1.076E-3)then
      increment_count=0
      WBEAM=cos((FLRNDM(XDUMMY)*2.5+2.5)/180*3.1415926)
      read(12,"(A1000)")temp1
      read(12,"(A1000)")temp1
      do i=1,1011
      read(12,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      read(temp(25:33),*)spec_count
      increment_count=increment_count+spec_count*(energy_high
     & -energy_low)
      neutron_spec_count(i+1)=increment_count
      end do
      do i=1,1012
      neu_spec_pro(i)=neutron_spec_count(i)/neutron_spec_count(1012)
      end do
      rewind(12)
      read(12,"(A1000)")temp1
      read(12,"(A1000)")temp1
      EneRndm=FLRNDM(XDUMMY)
      do i=1,1011     
      read(12,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      if(EneRndm.gt.neu_spec_pro(i).and.EneRndm.le.neu_spec_pro(i+1))
     & then
      energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY)
      end if
      end do
*************************************************************************
*the neutron moment angle is 2.5-5
************************************************************************
*************************************************************************
*the neutron moment angle is 5-7.5
************************************************************************
      else IF(DirRndm.gt.1.146e-4/1.076E-3.and.
     & DirRndm.le.2.645e-4/1.076E-3)THEN
      WBEAM=cos((FLRNDM(XDUMMY)*2.5+5)/180*3.1415926)
      increment_count=0
      read(13,"(A1000)")temp1
      read(13,"(A1000)")temp1
      do i=1,1011
      read(13,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      read(temp(25:33),*)spec_count
      increment_count=increment_count+spec_count*(energy_high
     & -energy_low)
      neutron_spec_count(i+1)=increment_count
      end do
      do i=1,1012
      neu_spec_pro(i)=neutron_spec_count(i)/neutron_spec_count(1012)
      end do
      rewind(13)
      read(13,"(A1000)")temp1
      read(13,"(A1000)")temp1
      EneRndm=FLRNDM(XDUMMY)
      do i=1,1011     
      read(13,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      if(EneRndm.gt.neu_spec_pro(i).and.EneRndm.le.neu_spec_pro(i+1))
     & then
      energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY)
      end if
      end do
*************************************************************************
*the neutron moment angle is 5-7.5
************************************************************************
*************************************************************************
*the neutron moment angle is 7.5-10
************************************************************************
      else IF(DirRndm.gt.2.645e-4/1.076E-3.and.
     & DirRndm.le.4.17e-4/1.076E-3)THEN
      WBEAM=cos((FLRNDM(XDUMMY)*2.5+7.5)/180*3.1415926)
      increment_count=0
      read(14,"(A1000)")temp1
      read(14,"(A1000)")temp1
      do i=1,1011
      read(14,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      read(temp(25:33),*)spec_count
      increment_count=increment_count+spec_count*(energy_high
     & -energy_low)
      neutron_spec_count(i+1)=increment_count
      end do
      do i=1,1012
      neu_spec_pro(i)=neutron_spec_count(i)/neutron_spec_count(1012)
      end do
      rewind(14)
      read(14,"(A1000)")temp1
      read(14,"(A1000)")temp1
      EneRndm=FLRNDM(XDUMMY)
      do i=1,1011     
      read(14,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      if(EneRndm.gt.neu_spec_pro(i).and.EneRndm.le.neu_spec_pro(i+1))
     & then
      energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY)
      end if
      end do
*************************************************************************
*the neutron moment angle is 7.5-10
************************************************************************
*************************************************************************
*the neutron moment angle is 10-12.5
************************************************************************
      else IF(DirRndm.gt.4.17e-4/1.076E-3.and.
     & DirRndm.le.6.173e-4/1.076E-3)THEN
      WBEAM=cos((FLRNDM(XDUMMY)*2.5+10)/180*3.1415926)
      increment_count=0
      read(15,"(A1000)")temp1
      read(15,"(A1000)")temp1
      do i=1,1011
      read(15,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      read(temp(25:33),*)spec_count
      increment_count=increment_count+spec_count*(energy_high
     & -energy_low)
      neutron_spec_count(i+1)=increment_count
      end do
      do i=1,1012
      neu_spec_pro(i)=neutron_spec_count(i)/neutron_spec_count(1012)
      end do
      rewind(15)
      read(15,"(A1000)")temp1
      read(15,"(A1000)")temp1
      EneRndm=FLRNDM(XDUMMY)
      do i=1,1011     
      read(15,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      if(EneRndm.gt.neu_spec_pro(i).and.EneRndm.le.neu_spec_pro(i+1))
     & then
      energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY)
      end if
      end do
*************************************************************************
*the neutron moment angle is 10-12.5
************************************************************************
*************************************************************************
*the neutron moment angle is 12.5-15
************************************************************************
      else IF(DirRndm.gt.6.173e-4/1.076E-3.and.
     & DirRndm.le.8.385e-4/1.076E-3)THEN
      WBEAM=cos((FLRNDM(XDUMMY)*2.5+12.5)/180*3.1415926)
      increment_count=0
      read(16,"(A1000)")temp1
      read(16,"(A1000)")temp1
      do i=1,1011
      read(16,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      read(temp(25:33),*)spec_count
      increment_count=increment_count+spec_count*(energy_high
     & -energy_low)
      neutron_spec_count(i+1)=increment_count
      end do
      do i=1,1012
      neu_spec_pro(i)=neutron_spec_count(i)/neutron_spec_count(1012)
      end do
      rewind(16)
      read(16,"(A1000)")temp1
      read(16,"(A1000)")temp1
      EneRndm=FLRNDM(XDUMMY)
      do i=1,1011     
      read(16,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      if(EneRndm.gt.neu_spec_pro(i).and.EneRndm.le.neu_spec_pro(i+1))
     & then
      energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY)
      end if
      end do
*************************************************************************
*the neutron moment angle is 12.5-15
************************************************************************
*************************************************************************
*the neutron moment angle is 15-17.5
************************************************************************
      else IF(DirRndm.gt.8.385e-4/1.076E-3)THEN
      WBEAM=cos((FLRNDM(XDUMMY)*2.5+15)/180*3.1415926)
      increment_count=0
      read(17,"(A1000)")temp1
      read(17,"(A1000)")temp1
      do i=1,1011
      read(17,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      read(temp(25:33),*)spec_count
      increment_count=increment_count+spec_count*(energy_high
     & -energy_low)
      neutron_spec_count(i+1)=increment_count
      end do
      do i=1,1012
      neu_spec_pro(i)=neutron_spec_count(i)/neutron_spec_count(1012)
      end do
      rewind(17)
      read(17,"(A1000)")temp1
      read(17,"(A1000)")temp1
      EneRndm=FLRNDM(XDUMMY)
      do i=1,1011     
      read(17,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      if(EneRndm.gt.neu_spec_pro(i).and.EneRndm.le.neu_spec_pro(i+1))
     & then
      energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY)
      end if
      end do
      end if
*************************************************************************
*the neutron moment angle is 15-17.5
************************************************************************
************************************************************************
* when the particle is neutron
************************************************************************
************************************************************************
* when the particle is photon
************************************************************************
************************************************************************
*the photon moment angle is 0-15
************************************************************************
      else if(ijbeam.eq.7)then
      DirRndm=FLRNDM(XDUMMY)
      if(DirRndm.le.1.762/6.906)then
      wbeam=cos(Flrndm(xdummy)*15/180*3.1415926)
      increment_count=0
      read(21,"(A1000)")temp1
      read(21,"(A1000)")temp1
      do i=1,1000
      read(21,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      read(temp(25:33),*)spec_count
      increment_count=increment_count+spec_count*(energy_high
     & -energy_low)
      gamma_spec_count(i+1)=increment_count
      end do
      do i=1,1001
      gamma_spec_pro(i)=gamma_spec_count(i)/gamma_spec_count(1001)
      end do
      rewind(21)
      read(21,"(A1000)")temp1
      read(21,"(A1000)")temp1
      EneRndm=FLRNDM(XDUMMY)
      do i=1,1000     
      read(21,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      if(EneRndm.gt.gamma_spec_pro(i).and.EneRndm.le.
     & gamma_spec_pro(i+1))then
      energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY)
      end if
      end do
************************************************************************
*the photon moment angle is 0-15
************************************************************************
************************************************************************
*the photon moment angle is 15-30
************************************************************************
      else if(DirRndm.gt.1.762/6.906)then
      wbeam=cos((Flrndm(xdummy)*15+15)/180*3.1415926)
      increment_count=0
      read(22,"(A1000)")temp1
      read(22,"(A1000)")temp1
      do i=1,1000
      read(22,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      read(temp(25:33),*)spec_count
      increment_count=increment_count+spec_count*(energy_high
     & -energy_low)
      gamma_spec_count(i+1)=increment_count
      end do
      do i=1,1001
      gamma_spec_pro(i)=gamma_spec_count(i)/gamma_spec_count(1001)
      end do
      rewind(22)
      read(22,"(A1000)")temp1
      read(22,"(A1000)")temp1
      EneRndm=FLRNDM(XDUMMY)
      do i=1,1000     
      read(22,"(A1000)")temp
      read(temp(3:11),*)energy_low
      read(temp(14:22),*)energy_high
      if(EneRndm.gt.gamma_spec_pro(i).and.EneRndm.le.
     & gamma_spec_pro(i+1))then
      energy=energy_low+(energy_high-energy_low)*FLRNDM(XDUMMY)
      end if
      end do
      end if
      end if
************************************************************************
*the photon moment angle is 15-30
************************************************************************      
************************************************************************
* when the particle is photon
************************************************************************
      PBEAM=SQRT(energy*(energy+2*am(ionid)))
*        PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
*    &                          + TWOTWO * AM (IONID) ) )      
      thi=flrndm(xdummy)*360
      ubeam=sin(thi/180.*3.1415926)*sqrt(1-wbeam*wbeam)
      vbeam=cos(thi/180.*3.1415926)*sqrt(1-wbeam*wbeam)
      TXFLK  (NPFLKA) = UBEAM
      TYFLK  (NPFLKA) = VBEAM
      TZFLK  (NPFLKA) = WBEAM
*     TZFLK  (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
*    &                       - TYFLK (NPFLKA)**2 )                             
         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) ) )
*  |  +----------------------------------------------------------------*
*  |  |  Check if it is a neutrino, if so force the interaction
*  |  |  (unless the relevant flag has been disabled)
         IF ( LISNUT (IJBEAM) .AND. LNUFIN ) THEN
            LFRPHN (NPFLKA) = .TRUE.
*  |  |
*  |  +----------------------------------------------------------------*
*  |  |  Not a neutrino
         ELSE
            LFRPHN (NPFLKA) = .FALSE.
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
      END IF
*  |
*  +-------------------------------------------------------------------*
*  From this point .....
*  Particle generation (1 for primaries)
      LOFLK  (NPFLKA) = 1
*  User dependent flag:
      LOUSE  (NPFLKA) = 0
*  No channeling:
      KCHFLK (NPFLKA) = 0
      ECRFLK (NPFLKA) = ZERZER
*  Extra infos:
      INFSTK (NPFLKA) = 0
      LNFSTK (NPFLKA) = 0
      ANFSTK (NPFLKA) = ZERZER
*  Parent variables:
      IPRSTK (NPFLKA) = 0
      EKPSTK (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
      AKNSHR (NPFLKA) = -TWOTWO
*  ahead particle Moment here
*  Polarization cosines:
      TXPOL  (NPFLKA) = -TWOTWO
      TYPOL  (NPFLKA) = +ZERZER
      TZPOL  (NPFLKA) = +ZERZER
*  Particle coordinates
      xbeam=0
      ybeam=0
      zbeam=0
      XFLK   (NPFLKA) = XBEAM
      YFLK   (NPFLKA) = YBEAM
      ZFLK   (NPFLKA) = ZBEAM
*  Calculate the total kinetic energy of the primaries: don't change
*  +-------------------------------------------------------------------*
*  |  (Radioactive) isotope:
      IF ( IJBEAM .EQ. -2 .AND. LRDBEA ) THEN
*  |
*  +-------------------------------------------------------------------*
*  |  Heavy ion:
      ELSE IF ( ILOFLK (NPFLKA) .EQ. -2 .OR.
     &          ILOFLK (NPFLKA) .GT. 100000 ) THEN
         TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
*  |
*  +-------------------------------------------------------------------*
*  |  Standard particle:
      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

