*$ 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.                                             *
*                                                                      *
*                                                                      *
*     New source for FLUKA9x-FLUKA20xy:                                *
*                                                                      *
*     Created on 07 January 1990   by    Alfredo Ferrari & Paola Sala  *
*                                                   Infn - Milan       *
*                                                                      *
*     Last change on  17-Oct-10    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.            *
*                                                                      *
*       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)'
*
      PARAMETER (NMAX=1000000)
      LOGICAL LFIRST
*       
      CHARACTER*250 LINE
      INTEGER NNN,IJ(NMAX)
      DIMENSION XXX(NMAX),YYY(NMAX),ZZZ(NMAX)
      DIMENSION UUU(NMAX),VVV(NMAX),WWW(NMAX)
      DIMENSION ERG(NMAX),WGT(NMAX)

      SAVE LFIRST
      SAVE XXX,YYY,ZZZ
      SAVE UUU,VVV,WWW
      SAVE ERG,WGT
      
      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 ***
         
         LUNRD = NINT(WHASOU(1))
         NNN = 0
 10      CONTINUE
         READ( LUNRD, '(A)', ERR=10, END=20 ) LINE
         READ (LINE,*,ERR=10) I, X, Y, Z, U, V, W, E, WG
         NNN = NNN + 1
         IF (NNN.GT.NMAX) CALL FLABRT('SOURCE','Increase NMAX')
         IJ(NNN)  = I
         XXX(NNN) = X
         YYY(NNN) = Y
         ZZZ(NNN) = Z
*     Normalize direction to 1.0
         UVW = SQRT(U**2 + V**2 + W**2)
         UUU(NNN) = U / UVW
         VVV(NNN) = V / UVW
         WWW(NNN) = W / UVW
         ERG(NNN) = E
         WGT(NNN) = WG
         PRINT*,"I,X,Y,Z,U,V,W,E,WG",I, X, Y, Z, U, V, W, E, WG
         GOTO 10
 20      CONTINUE
         IF (NNN.EQ.0) CALL FLABRT('SOURCE','Error reading file')
         WRITE (LUNOUT,*)
         WRITE (LUNOUT,*) '*** rdsource: ',NNN,' particles loaded'
         WRITE (LUNOUT,*)
      END IF
      
      NPFLKA = NPFLKA + 1
      RNDSIG = FLRNDM (RNDSIG)
      N = INT(NNN*RNDSIG)+1
      WTFLK (NPFLKA) = WGT(N)
      WEIPRI = WEIPRI + WTFLK (NPFLKA)
      PRINT*,"...DEBUG...DEBUG...WEIPRI:",WEIPRI,WGT(N)
      ILOFLK (NPFLKA) = IJ(N)

*     From this point ....
      LOFLK (NPFLKA) = 1           !     Generation is 1 for source particles
      LOUSE (NPFLKA) = 0           !     User variables: the user can set
      DO 100 ISPR = 1, MKBMX1      !     different values in the STUPRF or
         SPAREK (1,NPFLKA) = ZERZER !    STUPRE routine, but it is better
 100  CONTINUE                      !   not to do it here
      DO 200 ISPR = 1, MKBMX2       ! 
         ISPARK (ISPR,NPFLKA) = 0   !    More user variables (integer)
 200  CONTINUE
      NPARMA = NPARMA + 1           !    Updating the maximum particle number
      NUMPAR (NPFLKA) = NPARMA      !    Setting the current particle number
      NEVENT (NPFLKA) = 0           !    Resetting the current event number
      DFNEAR (NPFLKA) = +ZERZER     !    Resetting the distance to the
                                    !       nearest boundary
*     ... to this point: don't change anything
      AGESTK (NPFLKA) = +ZERZER       !  Particle age is zero by default
      AKNSHR (NPFLKA) = -TWOTWO       !  Resets the Kshort component of
                                      ! K0/K0bar. Usually not to be changed.
      IGROUP (NPFLKA) = 0             !  Group number for low energy
                                      ! neutrons: if set to 0, the program
                                      ! derives it from the kinetic energy

      IJBEAM = IJ(N)
      
*     +-------------------------------------------------------------------*
*     |  (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.
*     |  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
      
      TKEFLK (NPFLKA) = ERG(N)
      PMOFLK (NPFLKA) = SQRT(ERG(N)*(ERG(N)+TWOTWO*AM(IONID)))
      
      PRINT*,"...DEBUG...DEBUG...",RNDSIG,"N",N,'NPFLKA',NPFLKA
      PRINT*,"UUU(N)",UUU(N),"VVV(N)",VVV(N),"WWW(N)",WWW(N)
      PRINT*,"XXX(N)",XXX(N),"YYY(N)",YYY(N),"ZZZ(N)",ZZZ(N)
*     Cosines (tx,ty,tz)
      TXFLK (NPFLKA) = UUU(N)
      TYFLK (NPFLKA) = VVV(N)
      TZFLK (NPFLKA) = WWW(N)
*     Particle coordinates
*     
      XFLK (NPFLKA) = XXX(N)    !+ XBEAM
      YFLK (NPFLKA) = YYY(N) + TWOTWO*TWOTWO    !+ YBEAM
      ZFLK (NPFLKA) = ZZZ(N)    !+ ZBEAM

      TXPOL  (NPFLKA) = -TWOTWO
      TYPOL  (NPFLKA) = +ZERZER
      TZPOL  (NPFLKA) = +ZERZER

      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

