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

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Created on 07 january 1990   by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*     Last change on 21-jun-98     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.            *
*                                                                      *
*----------------------------------------------------------------------*
*
      INCLUDE '(AACOLL)'
      INCLUDE '(BEAM)'
      INCLUDE '(EPISOR)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(PAPROP)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(STACK)'
      INCLUDE '(STARS)'
*
      LOGICAL LFIRST
*
      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 ***
         WRITE(LUNOUT,*) 'Version 1 of Routine mysourcetest.f called'
      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
* 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) = ONEONE
      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) = SQRT ( PBEAM**2 + AM (IJBEAM)**2 ) - AM (IJBEAM)
* Particle momentum
      PMOM   (LSTACK) = PBEAM
*     PMOM (LSTACK) = SQRT ( TKE (LSTACK) * ( TKE (LSTACK) + TWOTWO
*    &                     * AM (ILO(LSTACK)) ) )
* Cosines (tx,ty,tz)
      TX     (LSTACK) = 0.0
      TY     (LSTACK) = -1.0
      TZ     (LSTACK) = 0.0
*     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) = 0.0
      YA (LSTACK) = 55.0
      ZA (LSTACK) = 0.0
*  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
*  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
      CALL SOEVSV
      RETURN
*=== End of subroutine Source =========================================*
      END