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

      INCLUDE '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
*
*----------------------------------------------------------------------*
*                                                                      *
*     Copyright (C) 1990-2010      by    Alfredo Ferrari & Paola Sala  *
*     All Rights Reserved.                                             *
*                                                                      *
*                                                                      *                               *
*     Source reading input from file made by usrrnc.f                  *
*     Argument list:						       *
*         IZ            : atomic number of the residual nucleus        *
*         IA            : mass number of the residual nucleus	       *
*         IS            : isomeric state of the residual nucleus       *
*         X, Y, Z       : particle position			       *
*         MREG          : number of the current region 		       *
*         WEE           : particle weight			       *
*         ICALL         : internal code calling flag		       *
*				(not for general use)                  * 
*                                                                      *
*     Created on 07 January 1990   by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*     Last change on  15-Sep-19    by    Enrico Mazza	               *
*                                                                      *
*  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)'
*
      LOGICAL LFIRST

*
      SAVE LFIRST
      DATA LFIRST / .TRUE. /

      INTEGER, DIMENSION(:), ALLOCATABLE :: IZ, IA, IS
      INTEGER :: A, B, C
      REAL, DIMENSION(:), ALLOCATABLE :: X, Y, Z
      REAL :: D, E, F
      INTEGER, DIMENSION(:), ALLOCATABLE :: MREG
      INTEGER :: G
      REAL, DIMENSION(:), ALLOCATABLE :: WEE
      REAL :: H
      INTEGER, DIMENSION(:), ALLOCATABLE :: ICALL
      INTEGER :: I
      INTEGER :: NNN, NMAX
      REAL :: SUM_WEE
      LOGICAL :: file_exists
      CHARACTER*250 LINE ! the line buffer
      CHARACTER*8 REGNAM

      SAVE IZ, IA, IS
      SAVE X, Y, Z
      SAVE MREG, WEE, ICALL
      SAVE NMAX

*======================================================================*
*                                                                      *
*                 BASIC VERSION                                        *
*                                                                      *
*======================================================================*
      NOMORE = 0
*  +-------------------------------------------------------------------*
*  |  First call initializations:
      IF ( LFIRST ) THEN
         WRITE (LUNOUT,*)
         WRITE(LUNOUT,*) '--- USER source called SourceISOTOPE v1.4 ---'
*  |  *** The following 3 cards are mandatory ***
         TKESUM = ZERZER
         LFIRST = .FALSE.
         LUSSRC = .TRUE.
*  |  *** User initialization ***
*
*  |  WHASOU(1) = UNIT number (File name is linked with OPEN)
         LUNRD = NINT(WHASOU(1))
         INQUIRE (UNIT = 99, EXIST=file_exists)
         WRITE (LUNOUT,*) '*** file inquire existence:', file_exists
	 WRITE (LUNOUT,*) '*** Source file OPEN: UNIT',INT(WHASOU(1))

* count # of lines in the input file         
         NMAX = 0
 11      CONTINUE
         READ(LUNRD, '(A)', ERR=9999, END=102 )
         NMAX = NMAX + 1
         GOTO 11
 102  CONTINUE

* use this number to allocate memory for arrays
         ALLOCATE (IZ(NMAX), IA(NMAX), IS(NMAX)) 
         ALLOCATE (X(NMAX), Y(NMAX), Z(NMAX)) 
         ALLOCATE (MREG(NMAX), WEE(NMAX), ICALL(NMAX))

* go on reading input file from the beginning
         REWIND(UNIT=99)      
         NNN = 0
         SUM_WEE = 0
 10      CONTINUE
                IF (NNN.EQ.0) THEN
                READ(LUNRD,*)
                ELSE
		READ(LUNRD, '(A)', ERR=9999, END=101 ) LINE
            	READ (LINE,*) A, B, C, D, E, F, G, H, I
		IF (NNN.GT.NMAX) CALL FLABRT('SOURCE','Error in NMAX')
		IZ(NNN) = A
		IA(NNN) = B
		IS(NNN) = C
		X(NNN) = D
		Y(NNN) = E
		Z(NNN) = F
		MREG(NNN) = G
		WEE(NNN) = H
		ICALL(NNN) = I
                END IF
                SUM_WEE = SUM_WEE + H
                NNN = NNN + 1
         GOTO 10
 101      CONTINUE
         IF (NNN.EQ.0) CALL FLABRT('SOURCE','Error reading file')
	 WRITE (LUNOUT,*) '*** Source file READ ***'
         WRITE (LUNOUT,*)
         WRITE (LUNOUT,*) '*** ISOTOPE source: ',NNN,' particles loaded'
         WRITE (LUNOUT,*) '*** SUM OF WEIGHTS: ',SUM_WEE
         WRITE (LUNOUT,*)      
       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
 300  CONTINUE
      RNDSIG = FLRNDM () 
      ISAMPL = INT(NMAX*RNDSIG) + 1 
      IF (IZ(ISAMPL) .EQ. ZERZER .OR. IA(ISAMPL) .EQ. ZERZER) GO TO 300
      
      NPFLKA = NPFLKA + 1


*  Wt is the weight of the particle (FROM FILE)
      WTFLK  (NPFLKA) = WEE(ISAMPL)
      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  = IA(ISAMPL)
         IZRES  = IZ(ISAMPL)
         IISRES = IS(ISAMPL)
         CALL STISBM ( IARES, IZRES, IISRES )
         IJHION = IZ(ISAMPL)  * 1000 + IA(ISAMPL)
         IJHION = IJHION * 100 + KXHEAV
         IONID  = IJHION
         CALL DCDION ( IONID )
         CALL SETION ( IONID )
*  |
*  +-------------------------------------------------------------------*
*  |  Heavy ion:
      ELSE IF ( IJBEAM .EQ. -2 ) THEN
         IJHION = IZ(ISAMPL)  * 1000 + IA(ISAMPL)
         IJHION = IJHION * 100 + KXHEAV
         IONID  = IJHION
         CALL DCDION ( IONID )
         CALL SETION ( 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
*  |
*  +-------------------------------------------------------------------*
*  |  Normal hadron:
      ELSE
         IONID = 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
      END IF
*  |
*  +-------------------------------------------------------------------*
*  From this point .....
*  Particle generation (1 for primaries)
      LOFLK  (NPFLKA) = 1
*  User dependent flag:
      LOUSE  (NPFLKA) = 0
*  No channeling:
      LCHFLK (NPFLKA) = .FALSE.
      DCHFLK (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
*  Particle age (s)
      AGESTK (NPFLKA) = +ZERZER
      AKNSHR (NPFLKA) = -TWOTWO
*  Particle momentum
      PMOFLK (NPFLKA) = +ZERZER
*  Kinetic energy of the particle (GeV)
      TKEFLK (NPFLKA) = SQRT ( PMOFLK(NPFLKA)**2 + AM (IONID)**2 )
     &                         - AM (IONID)

*  Cosines (tx,ty,tz)
      TXFLK  (NPFLKA) = UBEAM
      TYFLK  (NPFLKA) = VBEAM
      TZFLK  (NPFLKA) = WBEAM
*     TZFLK  (NPFLKA) = SQRT ( ONEONE - TXFLK (NPFLKA)**2
*    &                       - TYFLK (NPFLKA)**2 )
*  Polarization cosines:
      TXPOL  (NPFLKA) = -TWOTWO
      TYPOL  (NPFLKA) = +ZERZER
      TZPOL  (NPFLKA) = +ZERZER
*  Particle coordinates
      XFLK   (NPFLKA) = X(ISAMPL)
      YFLK   (NPFLKA) = Y(ISAMPL)
      ZFLK   (NPFLKA) = Z(ISAMPL)
*  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 =========================================*
 9999 CONTINUE
      CALL FLABRT('SOURCE','Error reading source file')
 8888 CONTINUE
      WRITE (LUNOUT,*) '--- WARNING: region name conversion FAILED ---'
      END
