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

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
      INCLUDE '(AACOLL)'
      INCLUDE '(BEAM)'
      INCLUDE '(EPISOR)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(PAPROP)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(STACK)'
      INCLUDE '(STARS)'
*
      LOGICAL LFIRST
      character*80 line
*
      DIMENSION en(0:556), spec(556)
      SAVE LFIRST, KOUNT
      DATA LFIRST / .TRUE. /
      DATA KOUNT / 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 ***
         WRITE(LUNOUT,*)
         WRITE(LUNOUT,'(A,132A)') ("*",i=1,132)
         WRITE(LUNOUT,*)
         WRITE(LUNOUT,*)
         IF(NINT(WHASOU(1)) .EQ. 1) THEN
            OPEN(UNIT = 81, FILE = "../BL1.spectrum", 
     &         STATUS = "OLD")
            WRITE(LUNOUT,'(A,A,A)') "    BL1 white spectrum ",
     &      "(file BL1.spectrum)" 
         END IF
         WRITE(LUNOUT,*)
         WRITE(LUNOUT,*)
         WRITE(LUNOUT,'(A,132A)') ("*",i=1,132)
         WRITE(LUNOUT,*)
         DO 2 I = 1, 231
            READ(81,*,END=999)  EN(I), SPEC(I)
            EN(I) = EN(I)*1.D-3
 2       CONTINUE
 999     CONTINUe
       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
      RAN01 = FLRNDM(TKESUM)
      ENERGY = 1.D-6 + ( 1.99526D-4 - 1.D-6) * RAN01
      DO 3 I = 1, 231
         IF(ENERGY .GT. EN(I)) THEN
            WEIGHT = HLFHLF*(SPEC(I)+SPEC(I+1))
         END IF
 3    CONTINUE
      KOUNT = KOUNT + 1
* Lstack is the stack counter: of course any time source is called it
* must be =0
      LSTACK = LSTACK + 1
* Wt is the weight of the particle
      WT     (LSTACK) = WEIGHT
      WEIPRI = WEIPRI + WT (LSTACK)
* Particle type (1=proton.....). Ijbeam is the type set by the BEAM
* card
*  +-------------------------------------------------------------------*
*  |  Heavy ion:
      IF ( IJBEAM .EQ. -2 ) THEN
         IJHION = IPROZ * 1000 + IPROA
         IJHION = IJHION * 100 + KXHEAV
         IONID  = IJHION
         CALL DCDION ( IONID )
         CALL SETION ( IONID )
         ILO (LSTACK) = IJHION
*  |
*  +-------------------------------------------------------------------*
*  |  Normal hadron:
      ELSE
         ILO (LSTACK) = IJBEAM
      END IF
*  |
*  +-------------------------------------------------------------------*
*
* From this point .....
* Particle generation (1 for primaries)
      LO     (LSTACK) = 1
* User dependent flag:
      LOUSE  (LSTACK) = 0
* User dependent spare variables:
      DO 100 ISPR = 1, MKBMX1
         SPAREK (ISPR,LSTACK) = ZERZER
 100  CONTINUE
* User dependent spare flags:
      DO 200 ISPR = 1, MKBMX2
         ISPARK (ISPR,LSTACK) = 0
 200  CONTINUE
* Save the track number of the stack particle:
      ISPARK (MKBMX2,LSTACK) = LSTACK
      NPARMA = NPARMA + 1
      NUMPAR (LSTACK) = NPARMA
      NEVENT (LSTACK) = 0
      DFNEAR (LSTACK) = +ZERZER
* ... to this point: don't change anything
* Particle age (s)
      AGESTK (LSTACK) = +ZERZER
      AKNSHR (LSTACK) = -TWOTWO
* Group number for "low" energy neutrons, set to 0 anyway
      IGROUP (LSTACK) = 0
* Kinetic energy of the particle (GeV)
      TKE(LSTACK) = ENERGY
* Particle momentum
      PMOM   (LSTACK) = ENERGY
*     PMOM (LSTACK) = SQRT ( TKE (LSTACK) * ( TKE (LSTACK) + TWOTWO
*    &                     * AM (ILO(LSTACK)) ) )
* Cosines (tx,ty,tz)
      TX     (LSTACK) = TINX
      TY     (LSTACK) = TINY
      TZ     (LSTACK) = TINZ
*     TZ     (LSTACK) = SQRT ( ONEONE - TX(LSTACK)**2 - TY(LSTACK)**2 )
* Polarization cosines:
      TXPOL (LSTACK) = -TWOTWO
      TYPOL (LSTACK) = +ZERZER
      TZPOL (LSTACK) = +ZERZER
* Particle coordinates
      XA (LSTACK) = XINA
      YA (LSTACK) = YINA
      ZA (LSTACK) = ZINA
*  Calculate the total kinetic energy of the primaries: don't change
      IF ( ILO(LSTACK) .EQ. -2 .OR. ILO(LSTACK) .GT. 100000 ) THEN
         TKESUM = TKESUM + TKE (LSTACK) * WT (LSTACK)
      ELSE IF ( ILO(LSTACK) .NE. 0 ) THEN
         TKESUM = TKESUM + ( TKE (LSTACK) + AMDISC (ILO(LSTACK)) )
     &          * WT (LSTACK)
      ELSE
         TKESUM = TKESUM + TKE (LSTACK) * WT (LSTACK)
      END IF
*      IF (KOUNT .LT. 100) WRITE (*,*) PMOM(LSTACK),WT(LSTACK)
*  Here we ask for the region number of the hitting point.
*     NREG (LSTACK) = ...
*  The following line makes the starting region search much more
*  robust if particles are starting very close to a boundary:
      CALL GEOCRS ( TX (LSTACK), TY (LSTACK), TZ (LSTACK) )
      CALL GEOREG ( XA (LSTACK), YA (LSTACK), ZA (LSTACK),
     &              NREG (LSTACK), IDISC )
*  Do not change these cards:
      CALL GEOHSM ( NHSPNT (LSTACK), 1, -11, MLATTC )
      NLATTC (LSTACK) = MLATTC
*     CMPATH (LSTACK) = ZERZER
      CALL SOEVSV
      RETURN
*=== End of subroutine Source =========================================*
      END
