source: ZHANGProjects/ICOSIM/CPP/trunk/flukaIntegration/fluka/source.f @ 2

Last change on this file since 2 was 2, checked in by zhangj, 10 years ago

Initial import
File size: 13.8 KB
Line 
1*$ CREATE SOURCE.FOR
2*COPY SOURCE
3*
4*=== source ===========================================================*
5*
6      SUBROUTINE SOURCE ( NOMORE )
7
8      INCLUDE '(DBLPRC)'
9      INCLUDE '(DIMPAR)'
10      INCLUDE '(IOUNIT)'
11*
12*----------------------------------------------------------------------*
13*                                                                      *
14*     Copyright (C) 1990-2009      by    Alfredo Ferrari & Paola Sala  *
15*     All Rights Reserved.                                             *
16*                                                                      *
17*                                                                      *
18*     New source for FLUKA9x-FLUKA20xy:                                *
19*                                                                      *
20*     Created on 07 january 1990   by    Alfredo Ferrari & Paola Sala  *
21*                                                   Infn - Milan       *
22*                                                                      *
23*     Last change on 08-feb-09     by    Alfredo Ferrari               *
24*                                                                      *
25*  This is just an example of a possible user written source routine.  *
26*  note that the beam card still has some meaning - in the scoring the *
27*  maximum momentum used in deciding the binning is taken from the     *
28*  beam momentum.  Other beam card parameters are obsolete.            *
29*                                                                      *
30*       Output variables:                                              *
31*                                                                      *
32*              Nomore = if > 0 the run will be terminated              *
33*                                                                      *
34*----------------------------------------------------------------------*
35*
36
37      INCLUDE '(BEAMCM)'
38      INCLUDE '(FHEAVY)'
39      INCLUDE '(FLKSTK)'
40      INCLUDE '(IOIOCM)'
41      INCLUDE '(LTCLCM)'
42      INCLUDE '(PAPROP)'
43      INCLUDE '(SOURCM)'
44      INCLUDE '(SUMCOU)'
45      INCLUDE 'flukaio.inc'
46
47      LOGICAL LFIRST
48
49      SAVE LFIRST
50      DATA LFIRST / .TRUE. /
51
52      CHARACTER*13 FILNAM
53      SAVE FILNAM
54      DATA FILNAM / 'coupling.stop' /
55
56      LOGICAL LEXIST
57      SAVE LEXIST
58      DATA LEXIST / .FALSE. /
59     
60      INTEGER ICOUNT, ITOTAL
61      COMMON /TRACKING/ ICOUNT, ITOTAL
62     
63      INTEGER IPRINT
64      LOGICAL LPRINT
65     
66
67*======================================================================*
68*                                                                      *
69*                 BASIC VERSION                                        *
70*                                                                      *
71*======================================================================*
72      NOMORE = 0
73*  +-------------------------------------------------------------------*
74*  |  First call initializations:
75      IF ( LFIRST ) THEN
76*  |  *** The following 3 cards are mandatory ***
77         TKESUM = ZERZER
78         LFIRST = .FALSE.
79         LUSSRC = .TRUE.
80*  |  *** User initialization ***
81
82*        Server Initialization
83
84         IPORT  = 0  ! Random assigned
85         ITOTAL = 0
86         ICOUNT = 0
87         ITURN  = 0
88
89         IPORT = 12345
90         IPORT = NTSTART(IPORT)
91         if (IPORT .EQ. -1 ) THEN
92           WRITE(LUNOUT,*) ' Error initializing server'
93           CALL FLUSH(LUNOUT)
94           GOTO 3900
95         END IF
96
97         WRITE(LUNOUT,*) ' '
98         WRITE(LUNOUT,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
99         WRITE(LUNOUT,*) ' !!       FlukaIO server initialized     !!'
100         WRITE(LUNOUT,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
101         WRITE(LUNOUT,*) ' '
102         WRITE(LUNOUT,*) ' Listening on port ', IPORT
103         WRITE(LUNOUT,*) ' '
104         WRITE(LUNOUT,*) ' To stop the simulation, touch file: '//
105     &                   FILNAM
106         WRITE(LUNOUT,*) ' '
107
108*        Write server parameters
109         CALL OAUXFI( '../network.nfo', 95, 'APPEND', IERR )
110         WRITE(95,*) IPORT                 ! Port number where the server is listening
111         CALL FFLUSH
112         CALL FLUSH(95)
113         CLOSE(95)
114
115*        parse configuration file:
116c         CALL PARSECONF()
117
118*        Initialize scrapers
119c         CALL SCR_INIT()
120
121*        Accept network connection
122         WRITE(LUNOUT,*) ' Waiting for connection'
123         CALL FLUSH(LUNOUT)
124
125         ICID = NTACCEPT()
126         IF (ICID .lt. 0) THEN
127           WRITE(LUNOUT,*) ' Error accepting connection'
128           CALL FLUSH(LUNOUT)
129           GOTO 3800
130         END IF
131
132c         NTTIMEOUT(ICID, 10000)
133         WRITE(LUNOUT,*) ' Connection accepted'
134         CALL FLUSH(LUNOUT)
135
136      END IF ! End of initialization
137
138 1972 CONTINUE ! Next message
139
140* A.Mereghetti, 04/05/2011
141* since BEAMPOS card is used, X/Y/Z/U/V/WBEAM are used
142*       as expected:
143        N = NTWAIT(
144     &        ICID,
145     &        MTYPE, 
146     &        IDP, IDGEN, WGT, 
147     &        XCBEAM, YCBEAM, ZCBEAM, 
148     &        UCBEAM, VCBEAM, WCBEAM,
149     &        IAA, IZZ, AMBEAM, PCBEAM,
150     &        T)
151        if(n.lt.0) then
152          write(LUNOUT,*) 'Client timeout'
153          CALL FLUSH(LUNOUT)
154          goto 3900
155        end if
156
157*A enlever...
158*       WRITE(LUNOUT,*)'PARAM: ', ICID, MTYPE, ITURN,
159*     &  IDP,IDGEN, WGT, XCBEAM, YCBEAM, ZCBEAM,
160*     &  UCBEAM, VCBEAM, WCBEAM, IAA, IZZ,
161*     &  AMBEAM, PCBEAM
162
163* +--   Parse message type and act consecuently
164* |
165        IF(MTYPE.EQ.N_PART) THEN
166c         Particle read, continue processing
167          goto 3000
168        end if
169
170        IF(MTYPE.EQ.N_TURN) THEN
171          goto 3100
172        end if
173       
174        IF(MTYPE.EQ.N_END) THEN
175          write(LUNOUT,*) '  Read end of computation'
176          goto 3700
177        end if
178* |
179* +--   End of message parsing
180     
181
1823000  continue ! A particle was read
183*  |
184*  +-------------------------------------------------------------------*
185*  Push one source particle to the stack. Note that you could as well
186*  push many but this way we reserve a maximum amount of space in the
187*  stack for the secondaries to be generated
188
189*     First particle of the turn, start of turn
190      IF(ICOUNT.EQ.0) THEN
191
192c       User requested run stop by creating coupling.stop file
193        INQUIRE ( FILE = FILNAM, EXIST = LEXIST )
194        IF (LEXIST) goto 4200 ! Inform error and close
195       
196*       Scraper start of turn
197c        CALL SCR_SOT()
198      END IF
199
200* Npflka is the stack counter: of course any time source is called it
201* must be =0
202      NPFLKA = NPFLKA + 1
203      ICOUNT = ICOUNT + 1
204
205* Wt is the weight of the particle
206      WTFLK  (NPFLKA) = ONEONE
207      WEIPRI = WEIPRI + WTFLK (NPFLKA)
208
209* Particle type (1=proton.....). Ijbeam is the type set by the BEAM
210* card
211*  +-------------------------------------------------------------------*
212*  |  (Radioactive) isotope:
213      IF ( IJBEAM .EQ. -2 .AND. LRDBEA ) THEN
214         IARES  = IPROA
215         IZRES  = IPROZ
216         IISRES = IPROM
217         CALL STISBM ( IARES, IZRES, IISRES )
218         IJHION = IPROZ  * 1000 + IPROA
219         IJHION = IJHION * 100 + KXHEAV
220         IONID  = IJHION
221         CALL DCDION ( IONID )
222         CALL SETION ( IONID )
223*  |
224*  +-------------------------------------------------------------------*
225*  |  Heavy ion:
226      ELSE IF ( IJBEAM .EQ. -2 ) THEN
227         IJHION = IPROZ  * 1000 + IPROA
228         IJHION = IJHION * 100 + KXHEAV
229         IONID  = IJHION
230         CALL DCDION ( IONID )
231         CALL SETION ( IONID )
232         ILOFLK (NPFLKA) = IJHION
233*  |  Flag this is prompt radiation
234         LRADDC (NPFLKA) = .FALSE.
235*  |  Group number for "low" energy neutrons, set to 0 anyway
236         IGROUP (NPFLKA) = 0
237*  |
238*  +-------------------------------------------------------------------*
239*  |  Normal hadron:
240      ELSE
241         IONID = IJBEAM
242         ILOFLK (NPFLKA) = IJBEAM
243*  |  Flag this is prompt radiation
244         LRADDC (NPFLKA) = .FALSE.
245*  |  Group number for "low" energy neutrons, set to 0 anyway
246         IGROUP (NPFLKA) = 0
247      END IF
248*  |
249*  +-------------------------------------------------------------------*
250
251
252
253* From this point .....
254* Particle generation (1 for primaries)
255      LOFLK  (NPFLKA) = 1
256* User dependent flag:
257      LOUSE  (NPFLKA) = 0
258*  No channeling:
259      LCHFLK (NPFLKA) = .FALSE.
260      DCHFLK (NPFLKA) = ZERZER
261* User dependent spare variables:
262      DO 100 ISPR = 1, MKBMX1
263         SPAREK (ISPR,NPFLKA) = ZERZER
264 100  CONTINUE
265* User dependent spare flags:
266      DO 200 ISPR = 1, MKBMX2
267         ISPARK (ISPR,NPFLKA) = 0
268 200  CONTINUE
269* Save the track number of the stack particle:
270      ISPARK (MKBMX2,NPFLKA) = NPFLKA
271      NPARMA = NPARMA + 1
272      NUMPAR (NPFLKA) = NPARMA
273      NEVENT (NPFLKA) = 0
274      DFNEAR (NPFLKA) = +ZERZER
275* ... to this point: don't change anything
276* Particle age (s)
277      AGESTK (NPFLKA) = +ZERZER
278      AKNSHR (NPFLKA) = -TWOTWO
279* Kinetic energy of the particle (GeV)
280      TKEFLK (NPFLKA) = SQRT ( PCBEAM**2 + AMBEAM**2 ) - AMBEAM
281* Particle momentum
282      PMOFLK (NPFLKA) = PCBEAM
283* Cosines (tx,ty,tz)
284      TXFLK  (NPFLKA) = UCBEAM
285      TYFLK  (NPFLKA) = VCBEAM
286      TZFLK  (NPFLKA) = WCBEAM
287* Polarization cosines:
288      TXPOL  (NPFLKA) = -TWOTWO
289      TYPOL  (NPFLKA) = +ZERZER
290      TZPOL  (NPFLKA) = +ZERZER
291* Particle coordinates
292      XFLK   (NPFLKA) = XBEAM + XCBEAM
293      YFLK   (NPFLKA) = YBEAM + YCBEAM
294      ZFLK   (NPFLKA) = ZBEAM + ZCBEAM
295
296*Test, a enlever apres!!!
297
298*      WRITE(LUNOUT,*)'Coordonnees:',XFLK (NPFLKA), YFLK (NPFLKA),
299*     &  ZFLK (NPFLKA), ZBEAM, ZCBEAM
300
301*
302* A.Mereghetti, 19/02/2010
303* the fort.99 file is written every IPRINT revolutions, as
304*     requested from the config.txt file:
305      IF ( LPRINT ) THEN
306*        first turn: print anyway
307         IF ( ITURN .EQ. 1 ) THEN
308            WRITE(99,'(2i7,1p,6(1x,e16.9))') IDP,ITURN,
309     &          XFLK(NPFLKA),YFLK(NPFLKA),
310     &          TXFLK(NPFLKA),TYFLK(NPFLKA),TZFLK(NPFLKA),TKEFLK(NPFLKA)
311            CALL FLUSH(99)
312           
313*        moving scraper?
314c         ELSEIF ( LMVSCR ) THEN
315*           print if:
316c            IF ( (  MOD( ITURN, IPRINT ).EQ.1 .AND.
317c     &          ( ITURN .GE. ISCINS ) ) .OR.
318c     &          ( IPRINT .EQ. 1 ) ) THEN
319c               WRITE(99,'(2i7,1p,6(1x,e16.9))') IDP,ITURN,
320c     &              XFLK(NPFLKA),YFLK(NPFLKA),
321c     &              TXFLK(NPFLKA),TYFLK(NPFLKA),TZFLK(NPFLKA),
322c     &              TKEFLK(NPFLKA)
323c               CALL FLUSH(99)
324c            ENDIF
325         
326*        not moving scraper?
327         ELSEIF ( ( MOD( ITURN, IPRINT ).EQ.1 ) .OR.
328     &          ( IPRINT .EQ. 1 ) ) THEN
329            WRITE(99,'(2i7,1p,6(1x,e16.9))') IDP,ITURN,
330     &          XFLK(NPFLKA),YFLK(NPFLKA),
331     &          TXFLK(NPFLKA),TYFLK(NPFLKA),TZFLK(NPFLKA),TKEFLK(NPFLKA)
332            CALL FLUSH(99)
333         ENDIF
334      ENDIF
335*
336*  Calculate the total kinetic energy of the primaries: don't change
337      IF ( ILOFLK (NPFLKA) .EQ. -2 .OR. ILOFLK (NPFLKA) .GT. 100000 )
338     &   THEN
339         TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
340      ELSE IF ( ILOFLK (NPFLKA) .NE. 0 ) THEN
341         TKESUM = TKESUM + ( TKEFLK (NPFLKA) + AMDISC (ILOFLK(NPFLKA)) )
342     &          * WTFLK (NPFLKA)
343      ELSE
344         TKESUM = TKESUM + TKEFLK (NPFLKA) * WTFLK (NPFLKA)
345      END IF
346      RADDLY (NPFLKA) = ZERZER
347
348
349*  Here we ask for the region number of the hitting point.
350*     NREG (NPFLKA) = ...
351*  The following line makes the starting region search much more
352*  robust if particles are starting very close to a boundary:
353      CALL GEOCRS ( TXFLK (NPFLKA), TYFLK (NPFLKA), TZFLK (NPFLKA) )
354      CALL GEOREG ( XFLK  (NPFLKA), YFLK  (NPFLKA), ZFLK  (NPFLKA),
355     &              NRGFLK(NPFLKA), IDISC )
356
357
358*  Do not change these cards:
359      CALL GEOHSM ( NHSPNT (NPFLKA), 1, -11, MLATTC )
360      NLATTC (NPFLKA) = MLATTC
361      CMPATH (NPFLKA) = ZERZER
362      CALL SOEVSV
363
364!     Scraper particle
365!     statistics on received particles
366c      CALL SCR_PART(XCBEAM, YCBEAM, ZCBEAM, UCBEAM, VCBEAM, WCBEAM)
367
368      ! Skip the rest, return with the particle
369      RETURN
370
3713100  continue ! End of turn
372      ITURN = IDP
373c      CALL SCR_EOT()
374      N = NTSENDEOT(ICID, ITURN)
375
376      write(LUNOUT,*) '  End of turn ', ITURN, ' received ', ICOUNT, 
377     &                ' particles'
378
379      ITOTAL = ITOTAL + ICOUNT
380      ICOUNT = 0
381
382      ! iterate again
383      goto 1972
384
385
3863700  continue ! the tracking is over, FLUKA is going to end
387
388      write(LUNOUT,*) '** Sending end of connection'
389      N = NTSENDEOC(ICID)
390
391c      CALL SCR_EOC()
392
3933800  continue ! Close connection
394
395      write(LUNOUT,*) '** Closing connection'
396      N = NTEND(ICID)
397      ICID = -1
398
3993900  continue ! Close server
400
401      write(LUNOUT,*) '** Shuting server down'
402      N = NTSHDWN()
403      NOMORE = 1
404
405      IF (ITOTAL.eq.0) THEN
406        write(LUNOUT, *) '** 0 Particles received STOPPING Fluka'
407        STOP
408      END IF
409
4104000  continue
411
412      RETURN
413
414*  +-------------------------------------------------------------------*
415*  Exceptional exit, send and wait for end of computation
4164200  continue
417
418      write(LUNOUT,*) '** coupling.stop detected, ',
419     &    'ending simulation'
420      write(LUNOUT,*) "** Sending End Of Computation"
421
422      N = NTSENDEOC(ICID)
423      if(N.lt.0) then
424        write(LUNOUT,*) 'Error sending end of turn'
425        CALL FLUSH(LUNOUT)
426        goto 3900
427      end if
428
4294500  continue
430
431*     Wait for an end of computation message
432      N = NTWAIT(
433     &        ICID,
434     &        MTYPE, 
435     &        IDP, IDGEN, WGT, 
436     &        XCBEAM, YCBEAM, ZCBEAM, 
437     &        UCBEAM, VCBEAM, WCBEAM,
438     &        IAA, IZZ, AMBEAM, PCBEAM,
439     &        T)
440      if(N.lt.0) then
441        write(LUNOUT,*) 'Client timeout'
442        CALL FLUSH(LUNOUT)
443        goto 3900
444      end if
445
446      IF(MTYPE.NE.N_END) THEN
447        goto 4500 ! Read next message
448      end if
449
450      goto 3800
451*  +-------------------------------------------------------------------*
452
453*=== End of subroutine Source =========================================*
454      END
Note: See TracBrowser for help on using the repository browser.