*$ CREATE SOURCE.FOR
*COPY SOURCE
*
*=== source ===========================================================*
*
      SUBROUTINE SOURCE ( NOMORE )      
*
*----------------------------------------------------------------------*
*  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 '(DBLPRC)'
      INCLUDE '(DIMPAR)'
      INCLUDE '(IOUNIT)'
      INCLUDE '(BEAMCM)'
      INCLUDE '(FHEAVY)'
      INCLUDE '(FLKSTK)'
      INCLUDE '(IOIOCM)'
      INCLUDE '(LTCLCM)'
      INCLUDE '(PAPROP)'
      INCLUDE '(SOURCM)'
      INCLUDE '(SUMCOU)'
      include '(CASLIM)'
      include '(SOUEVT)'
*
      LOGICAL LFIRST
*
      SAVE LFIRST
      DATA LFIRST / .TRUE. /

c      parameter(nlines=1000000)
c      save line
c      double precision enrgy(nlines),x(nlines)
c      double precision y(nlines),z(nlines),cx(nlines),cy(nlines)
c      double precision cz(nlines)
c      double precision Nstop
*======================================================================*
*                                                                      *
*                 FRA VERSION                                          *
*                                                                      *
*======================================================================*
      NOMORE = 0
*  +-------------------------------------------------------------------*
******************************************************************
************ Francesca declarations ******************
*Luke, you can pass info to the source.f file from the SOURCE card 

************ END declarations ************************
******************************************************************
*
*  |  First call initializations:
      IF ( LFIRST ) THEN
*  |  *** The following 3 cards are mandatory ***
         TKESUM = ZERZER
         LFIRST = .FALSE.
         LUSSRC = .TRUE.
***** Fra initialization: reading particle info from a file.dat *******
c         line = 1.0
c         CALL OAUXFI('xtube2.dat',LUNRDB,'OLD',IERR)
c         DO I=1,nlines
c         READ(LUNRDB,*,END=197)enrgy(I),x(I),y(I),z(I),
c     &           cx(I),cy(I),cz(I)
c 22      format(d14.7,d16.8,d16.8,d16.8,d16.8,d16.8,d16.8)
c*          write(*,*)enrgy(I),
c*     &      x(I),y(I),z(I),cx(I),cy(I),cz(I)
c         ENDDO
c 197     CONTINUE
c         Nstop=i-1
c         NCASES=Nstop
c         write(*,*)'ncases', NCASES

      write(*,*) "reading data from ", SDUSOU
      idata = WHASOU(6)
      call oauxfi(SDUSOU, idata, 'old', ierr)
      if (ierr.gt.0) then
      	call flabrt('source', 'can''t open SDUSOU')
      end if

************* END my initialization ************************************
      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
***********************************************************************
      NPFLKA = NPFLKA + 1

      read(idata, 22) enrgyy,xx,yy,zz,cxx,cyy,czz
22    format(d14.7,d16.8,d16.8,d16.8,d16.8,d16.8,d16.8)

      write(71,*) enrgyy,xx,yy,zz,cxx,cyy,czz

***********************************************************************
***********************************************************************
* Wt is the weight of the particle
      WTFLK  (NPFLKA) = ONEONE
      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  = 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:
      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
***** source particle age (no kaons) ***
      AGESTK (NPFLKA) = +ZERZER
* AGE for Kaons*************************
      AKNSHR (NPFLKA) = -TWOTWO
************************************************************************
* Group number for "low" energy neutrons, set to 0 anyway
      IGROUP (NPFLKA) = 0
************************************************************************
************************************************************************
*.......................................................................
**** Kinetic energy, direction cosines and position of the particles ***
c      TKEFLK (NPFLKA) = enrgy(line)
c      write(*,*)'Ef=',TKEFLK(NPFLKA),' ncase=',ncase
c      TXFLK  (NPFLKA) = cx(line)
c      TYFLK  (NPFLKA) = cy(line)
c      TZFLK  (NPFLKA) = cz(line)      
c*      write(*,*)'txf=',TXFLK(NPFLKA),' tyf=',TYFLK(NPFLKA),
c*     & ' tzf=',TZFLK(NPFLKA)
c      XFLK(NPFLKA)= x(line)
c      YFLK(NPFLKA)= y(line)
c      ZFLK(NPFLKA) = z(line)
c*      write(*,*)'xf=',XFLK(NPFLKA),' yf=',YFLK(NPFLKA),
c*     &     ' zf=',ZFLK(NPFLKA)
c      if(ncase.gt.Ncases)then
c         nomore=1
c      else
c         line = line+1
c      endif
************************************************************************
      TKEFLK (NPFLKA) = enrgyy
      TXFLK  (NPFLKA) = cxx
      TYFLK  (NPFLKA) = cyy
      TZFLK  (NPFLKA) = czz
      XFLK(NPFLKA) = xx
      YFLK(NPFLKA) = yy
      ZFLK(NPFLKA) = zz
*....................................................................... 
************************************************************************   
* Particle momentum ****************************************************
      PMOFLK (NPFLKA) = SQRT ( TKEFLK (NPFLKA) * ( TKEFLK (NPFLKA)
     &                       + TWOTWO * AM (IONID) ) )
************************************************************************
* Polarization cosines *************************************************
      TXPOL  (NPFLKA) = -TWOTWO
      TYPOL  (NPFLKA) = +ZERZER
      TZPOL  (NPFLKA) = +ZERZER
******* ***************************************************************
*********************************************************************
*  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
 555  continue
      RETURN
*=== End of subroutine Source =========================================*
      END

