       SUBROUTINE SOURCE (NOMORE)
       INCLUDE '(DBLPRC)' !CONTAINING PHYSICAL CONSTANTS ,THEIR VALUES AND DEFINITION FOR DOUBLE AND SINGLE PRECISION.
       INCLUDE '(IOUNIT)'
       INCLUDE '(DIMPAR)'
       INCLUDE '(BEAMCM)' !CONTAINING BEAM INFORMATION LIKE MOMENTUM, ENERGY, SPECTRUM SHAPE, RADIUS ETC
       INCLUDE '(FHEAVY)'
       INCLUDE '(FLKSTK)'
       INCLUDE '(SOURCM)'
       INCLUDE '(LTCLCM)'
       INCLUDE '(PAPROP)'
       INCLUDE '(SUMCOU)'
       INCLUDE '(IOIOCM)'
       LOGICAL LFIRST
       PARAMETER (NDP=25) !NUMBER OF DATA POINTS=24+1
       DIMENSION ENERN(0:NDP),DATA0(0:NDP), DATA5(0:NDP)
       DIMENSION DATA10(0:NDP), DATA20(0:NDP-9), ANG(0:3),CUMUDAT0(0:NDP)
       DIMENSION CUMUDAT5(0:NDP),CUMUDAT10(0:NDP),CUMUDAT20(0:NDP)
*WHERE ENERN=ENERGY SPECTRUM OF NEUTRONS, CUMUN=COMMULATIVE DISTRIBUTION OF NEUTRONS(NORMALIZED)
       SAVE LFIRST
       DATA LFIRST / .TRUE. /
       DATA ENERN /0.0D0, 3.9700E-3,
      & 5.0235D-3, 5.9617D-3, 7.1081E-3, 8.0485D-3,
      & 9.1061D-3, 1.0169D-2, 1.1016D-2, 1.2070D-2,
      & 1.3140D-2, 1.4133D-2, 1.5118D-2, 1.6082D-2,
      & 1.7154D-2, 1.8130D-2, 1.9105D-2, 2.0069D-2,
      & 2.1136D-2, 2.1989D-2, 2.3057D-2, 2.4128D-2,
      & 2.5094D-2, 2.6056D-2, 2.7017D-2, 2.7976D-2 /
*NUMBER OF NEUTRONS PER STERADIAN (FOR 28 MEV INCIDENT DEUTRONS)
* DATA FOR 0 DEGREE
       DATA DATA0 /0.0D0,
      & 1.07401E+11, 1.22756E+11, 1.51138E+11, 1.86422E+11,
      & 2.11738E+11, 2.20962E+11, 2.22523E+11, 2.28676E+11,
      & 2.42498E+11, 2.33329E+11, 1.81237E+11, 1.39875E+11,
      & 1.29936E+11, 1.19234E+11, 9.09011E+10, 6.56341E+10,
      & 5.49288E+10, 5.11240E+10, 4.80800E+10, 4.19759E+10,
      & 3.12734E+10, 1.90353E+10, 1.13956E+10, 6.82159E+9,
      & 3.78042E+9 /
*DATA FOR 5 DEGREE
       DATA DATA5 /0.0D0,
      & 1.41573D+11, 1.10896D+11, 1.22469D+11, 1.56327D+11,
      & 1.86338D+11, 2.01752D+11, 1.94892D+11, 1.84959D+11,
      & 1.83467D+11, 1.81984D+11, 1.65905D+11, 1.38301D+11,
      & 1.07620D+11, 8.46151D+10, 7.08407D+10, 5.93656D+10,
      & 4.86588D+10, 3.94938D+10, 3.26282D+10, 2.65311D+10,
      & 2.12074D+10, 1.58784D+10, 1.05548D+10, 6.76237D+9,
      & 2.96991D+9 /
*DATA FOR  10 DEGREE
       DATA DATA10 /0.0D0,
      & 5.42155D+10, 7.76617D+10, 1.00339D+11, 1.17255D+11,
      & 1.24183D+11, 1.25347D+11, 1.24593D+11, 1.18842D+11,
      & 1.07332D+11, 1.00045D+11, 1.00442D+11, 8.85469D+10,
      & 5.89804D+10, 4.09367D+10, 3.86468D+10, 4.05785D+10,
      & 3.63666D+10, 2.94656D+10, 2.37144D+10, 1.79646D+10,
      & 1.18334D+10, 7.23462D+9,  4.55928D+9,  3.03497D+9,
      & 3.04724D+9 /
*DATA FOR 20 DEGREE
       DATA DATA20 /0.0D0,
      & 6.05087D+10, 6.09061D+10, 6.32303D+10, 6.55543D+10,
      & 6.67219D+10, 6.32699D+10, 5.55805D+10, 4.86616D+10,
      & 4.55945D+10, 4.36816D+10, 4.02281D+10, 3.44643D+10,
      & 2.79316D+10, 2.33244D+10,2.02561D+10,  1.79580D+10 /
*DATA FOR DIRECTION COSINES
       DATA ANG / 0.939692621D0,
      & 0.984807753D0, 0.996194698D0, 1.0D0 /
*CUMULATIVE ANGULAR DISTRIBUTION
*      DATA CUMUANG / 0.0D0,
*    & 0.3266372518D0, 0.6148960479D0, 1.0D0



*---------------------------------------------------------------------*
       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,*)
       END IF
       NPFLKA=NPFLKA + 1 !INCREAMENTING COUNTER
       WTFLK (NPFLKA)=ONEONE
       WEIPRI=WEIPRI + WTFLK (NPFLKA)
*--------------------CUMULATIVE DATA DISTRIBUTIONS---------------------------------*
*      COMMULATIVE DISTRIBUTION CONSTRUCTION
       CUMUDAT0(0)=0.0D0
       CUMUDAT5(0)=0.0D0
       CUMUDAT10(0)=0.0D0
       CUMUDAT20(0)=0.0D0
       DO I=1,NDP
       CUMUDAT0(I)=CUMUDAT0(I-1)+(DATA0(I-1)+DATA0(I))*
      & (ENERN(I)-ENERN(I-1))*0.5
*COMMULATIVE DISTRIBUTION FOR 5 DEGREE
      CUMUDAT5(I)=CUMUDAT5(I-1)+(DATA5(I-1)+DATA5(I))*
      & (ENERN(I)-ENERN(I-1))*0.5
*COMMULATIVE DISTRIBUTION FOR 10 DEGREE
      CUMUDAT10(I)=CUMUDAT0(I-1)+(DATA0(I-1)+DATA0(I))*
      & (ENERN(I)-ENERN(I-1))*0.5
*COMMULATIVE DISTRIBUTION FOR 20 DEGREE
       CUMUDAT20=CUMUDAT0(I-1)+(DATA0(I-1)+DATA0(I))*
      & (ENERN(I)-ENERN(I-1))*0.5
       END DO
*-----------NORMALIZATION OF CUMULATIVE DISTRIBUTIONS------------------*
       NORM0=CUMUDAT0(NDP)
       NORM5=CUMUDAT5(NDP)
       NORM10=CUMUDAT10(NDP)
       NORM20=CUMUDAT20(NDP)
       DO I=1,NDP
       CUMUDAT0(I)=CUMUDAT0(I)/NORM0
       CUMUDAT5(I)=CUMUDAT5(I)/NORM5
       CUMUDAT10(I)=CUMUDAT10(I)/NORM10
       CUMUDAT20(I)=CUMUDAT20(I)/NORM20
       END DO
       WRITE(90,*) "INTEGERS, ENERGIES CALCULATED, NORMALIZED CULATIVE DISTRIBUTION"
       DO I=1, NDP
       WRITE (90,'(13,1X,1P,2G15.9)') I, ENERN(I), CUMUDAT0(I)
      &                     CUMUDAT5(I), CUMUDAT10(I), CUMUDAT20 (I)
       END DO
*----------------------K.E SAMPLING------------------------------*
* KINETIC ENERGY OF THE PARTICLE (GEV). SAMPLE FROM THE CUMULATIVE
*----------------------------------------------------------------*
       PHI=TWOPIP*FLRNDM(DUMMY)
       COSTHE=ANG(3)
       SINTHE=SQRT(ONEONE-COSTHE**2)
* COSINES (TX,TY,TZ)
       TXFLK (NPFLKA) = COSTHE*COS(PHI)
       TYFLK (NPFLKA) = SINTHE*SIN(PHI)
       TZFLK (NPFLKA) = SQRT(ONEONE-TXFLK (NPFLKA)**2-TYFLK (NPFLKA)**2)
       AI = FLRNDM(DUMMY)
          DO 100 I = 1, NDP
             IF(AI .GT. CUMUNDAT0(I-1) .AND. AI .LE. CUMUNDAT0(I)) THEN
       KIENER = ENERN(I-1)+(ENERN(I)-ENERN(I-1))*
      & (SQRT(AI) - SQRT(CUMUN(I-1)))/
      & (SQRT(CUMUN(I))-SQRT(CUMUN(I-1)))
       GO TO 110
       END IF
  100  CONTINUE
       STOP 'FAILED TO SAMPLE ENERGY'
  110  CONTINUE
       KOUNT=KOUNT + 1
*---------------------
*      BI = FLRNDM(DUMMY)
       COSTHE=ANG(2)
       SINTHE=SQRT(ONEONE-COSTHE**2)
       DO 200 I = 1, NDP
       IF(AI .GT. CUMUNDAT5(I-1) .AND. AI .LE. CUMUNDAT5(I)) THEN
       KIENER = ENERN(I-1)+(ENERN(I)-ENERN(I-1))*
      & (SQRT(AI) - SQRT(CUMUN(I-1)))/
      & (SQRT(CUMUN(I))-SQRT(CUMUN(I-1)))
       GO TO 210
       END IF
  200 CONTINUE
       STOP 'FAILED TO SAMPLE ENERGY'
  210 CONTINUE
     KOUNT=KOUNT + 1
       TXFLK (NPFLKA) = COSTHE*COS(PHI)
       TYFLK (NPFLKA) = SINTHE*SIN(PHI)
       TZFLK (NPFLKA) = SQRT(ONEONE-TXFLK (NPFLKA)**2-TYFLK (NPFLKA)**2)
*--------------------------
*      CI = FLRNDM(DUMMY)
        COSTHE=ANG(1)
        SINTHE=SQRT(ONEONE-COSTHE**2)
       DO 300 I = 1, NDP
       IF(AI .GT. CUMUNDAT10(I-1) .AND. AI .LE. CUMUNDAT10(I)) THEN
       KIENER = ENERN(I-1)+(ENERN(I)-ENERN(I-1))*
      & (SQRT(AI) - SQRT(CUMUN(I-1)))/
      & (SQRT(CUMUN(I))-SQRT(CUMUN(I-1)))
       GO TO 310
       END IF
  300 CONTINUE
       STOP 'FAILED TO SAMPLE ENERGY'
  310 CONTINUE
       KOUNT=KOUNT + 1
       TXFLK (NPFLKA) = COSTHE*COS(PHI)
       TYFLK (NPFLKA) = SINTHE*SIN(PHI)
       TZFLK (NPFLKA) = SQRT(ONEONE-TXFLK (NPFLKA)**2-TYFLK (NPFLKA)**2)
*-----    ------------------
*         DI = FLRNDM(DUMMY)
        COSTHE=ANG(0)
        SINTHE=SQRT(ONEONE-COSTHE**2)
       DO 400 I = 1, NDP
       IF(AI .GT. CUMUNDAT20(I-1) .AND. DI .LE. CUMUNDAT20(I)) THEN
       KIENER = ENERN(I-1)+(ENERN(I)-ENERN(I-1))*
      & (SQRT(AI) - SQRT(CUMUN(I-1)))/
      & (SQRT(CUMUN(I))-SQRT(CUMUN(I-1)))
       GO TO 410
       END IF
  400 CONTINUE
       STOP 'FAILED TO SAMPLE ENERGY'
  410 CONTINUE
       KOUNT=KOUNT + 1
       TXFLK (NPFLKA) = COSTHE*COS(PHI)
       TYFLK (NPFLKA) = SINTHE*SIN(PHI)
       TZFLK (NPFLKA) = SQRT(ONEONE-TXFLK (NPFLKA)**2-TYFLK (NPFLKA)**2)
*--------------------------------------------------------------------*
*--------------------------------------------------------------------*
* 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:
  100  DO ISPR = 1, MKBMX1
       SPAREK (ISPR,NPFLKA) = ZERZER
  100  CONTINUE
* USER DEPENDENT SPARE FLAGS:
  200  DO 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 PARTICLES
       TKEFLK (NPFLKA)=KIENER
       PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
                                 *  + TWOTWO * AM (IONID) ) )
*--------------------------------------------------------------

* POLARIZATION COSINES:
       TXPOL (NPFLKA) = -TWOTWO
       TYPOL (NPFLKA) = +ZERZER
       TZPOL (NPFLKA) = +ZERZER
* PARTICLE COORDINATES, COORDINATES SPECIFIED IN BEAMPOS
       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

