c $Id: urqmd.f,v 1.22 2001/04/06 22:38:54 weber Exp $
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
cdh program UrQMD
subroutine UrQMD(iflbmax)
c
c Authors : The UrQMD collaboration
c S.A. Bass, M. Belkacem, M. Bleicher, M. Brandstetter,
c L. Bravina, C. Ernst, L. Gerland, M. Hofmann,
c S. Hofmann, J. Konopka, G. Mao, L. Neise, S. Soff,
c C. Spieles, H. Weber, L.A. Winckelmann, H. Stoecker
c and W. Greiner
c
c Revision: 1.2
c
cc contact address:
cc
cc urqmd@th.physik.uni-frankfurt.de
cc
c
c This is the main module of {\tt urqmd}. It servers as a connection between
c the initialization, the propagation (including the real part of the
c optical potential) and the collision term (imaginary part of the optical
c potential).
c
c iflbmax: flag for retrying interaction after non-interaction
c 0 = do retry until interaction happens
c 1 = do not ( propagate particle, retry then)
c
C modifications for use in connection with CORSIKA by D. Heck
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
implicit none
include 'coms.f'
include 'comres.f'
include 'options.f'
include 'colltab.f'
include 'inputs.f'
include 'newpart.f'
include 'boxinc.f'
c
integer i,j,k,steps,ii,ocharge,ncharge, mc, mp, noc, it1,it2
real*8 sqrts,otime,xdummy,st
logical isstable
integer stidx,iflbmax
real*8 Ekinbar, Ekinmes, ESky2, ESky3, EYuk, ECb, EPau
common /energies/ Ekinbar, Ekinmes, ESky2, ESky3, EYuk, ECb, EPau
integer cti1sav,cti2sav
common/cncc/ncc
integer*8 ncc
c
c numerical/technical initialisation
c
cdh call uinit(0)
cdh call osc_header
cdh call osc99_header
c
c Main program
c
noc=0
1 mc=0
mp=0
c
c loop over all events
c
cdh do 10 event=1,nevents
c start event here
c
c time is the system time at the BEGINNING of every timestep
time = 0.0
c initialize random number generator
c call auto-seed generator only for first event and if no seed was fixed
cdh if(.not.firstseed.and.(.not.fixedseed)) then
cdh ranseed=-(1*abs(ranseed))
cdh call sseed(ranseed)
cdh else
cdh firstseed=.false.
cdh endif
cdh
cdh write(6,*)'event# ',event,ranseed
c
c initialisation of physics quantities
c
call init
c old time if an old fort.14 is used
if(CTOption(40).eq.1)time=acttime
c output preparation
c write headers to file
call output(13)
call output(14)
call output(15)
cdh call output(16)
cdh if(event.eq.1)call output(17)
cdh call osc99_event(-1)
c for CTOption(4)=1 : output of initialization configuration
if(CTOption(4).eq.1)call file14out(0)
c participant/spectator model:
cdh if(CTOption(28).ne.0) call rmspec(0.5d0*bimp,-(0.5d0*bimp))
c compute time of output
otime = outsteps*dtimestep
c reset time step counter
steps = 0
c loop over all timesteps
do 20 steps=1,nsteps
c store coordinates in arrays with *_t
c this is needed for MD type propagation
if (eos.ne.0) then
do 23 j=1,npart
r0_t(j) = r0(j)
rx_t(j) = rx(j)
ry_t(j) = ry(j)
rz_t(j) = rz(j)
23 continue
end if
c we are at the beginning of the timestep, set current time (acttime)
acttime = time
c option for MD without collision term
if(CTOption(16).ne.0) goto 103
c Load collision table with next collisions in current timestep
call colload
c check for collisions in time-step, nct = # of collisions in table
if (nct.gt.0) then
c entry-point for collision loop in case of full colload after every coll.
101 continue
k = 0
c normal entry-point for collision loop
100 continue
c get next collision
call getnext(k)
c exit collision loop if no collisions are left
if (k.eq.0) goto 102
c propagate all particles to next collision time
c store actual time in acttime, propagation time st=cttime(k)-acttime
st=cttime(k)-acttime
call cascstep(acttime,st)
c new actual time (for upcoming collision)
acttime = cttime(k)
c perform collision
if(cti2(k).gt.0.)then
if(abs(sqrts(cti1(k),cti2(k))-ctsqrts(k)).gt.1d-3)then
write(6,*)' ***(E) wrong collision update (col) ***'
write(6,*)cti1(k),cti2(k),
& ctsqrts(k),sqrts(cti1(k),cti2(k))
endif
else if(cti2(k).eq.0.and.
& abs(fmass(cti1(k))-ctsqrts(k)).gt.1d-3) then
write(6,*)' *** main(W) wrong collision update (decay)'
write(6,*)ctag,cti1(k),ityp(cti1(k)),dectime(cti1(k)),
& fmass(cti1(k)),ctsqrts(k)
endif
ocharge=charge(cti1(k))
if(cti2(k).gt.0) ocharge=ocharge+charge(cti2(k))
c store quantities in local variables for charge conservation check
it1= ityp(cti1(k))
if(cti2(k).gt.0)it2= ityp(cti2(k))
c increment "dirty" collision counter
if(cti2(k).gt.0)then !scatter !scatter
mc=mc+1
endif
c perform scattering/decay
ncc = ncc+1
cti1sav = cti1(k) ! hjd
cti2sav = cti2(k) ! hjd
call scatter(cti1(k),cti2(k),ctsigtot(k),ctsqrts(k),
& ctcolfluc(k))
c
c update collision table
c
c normal update mode
if(CTOption(17).eq.0) then
if(nexit.eq.0) then
c new collision partners for pauli-blocked states (nexit=0)
c hjd1 and cdh
if (cti1(k).ne.cti1sav.or.cti2(k).ne.cti2sav) then
goto 1
c cti1(k) = cti1sav !!! hjd1
c cti2(k) = cti2sav !!! hjd1
endif
c hjd1 and cdh
call collupd(cti1(k),1)
if(cti2(k).gt.0) call collupd(cti2(k),1)
else
ncharge=0
c new collision partners for scattered/produced particles (nexit><0)
do 30 i=1,nexit
c ncharge is used for charge conservation check
ncharge=ncharge+charge(inew(i))
call collupd(inew(i),1)
30 continue
c charge conservation check
if(ocharge.ne.ncharge) then
write(6,*)'ch-conservation error coll/dec ',ctag
write(6,*)' it1:',it1,' it2:',it2
write(6,*)' ch:',ocharge,ncharge
write(6,*)'cti1(k),cti2(k),ctsigtot(k),ctsqrts(k)'
write(6,*)cti1(k),cti2(k),ctsigtot(k),ctsqrts(k)
endif
endif
c update collisions for partners of annihilated particles
do 55 ii=1,nsav
call collupd(ctsav(ii),1)
55 continue
nsav=0
else ! (CTOption(17).ne.0)
c full collision load
call colload
endif
if (CTOption(17).eq.0) goto 100
goto 101
c this is the point to jump to after all collisions in the timestep
c have been taken care of
102 continue
endif ! (nct.gt.0)
c After all collisions in the timestep are done, propagate to end of
c the timestep.
c point to jump to in case of MD without collision term
103 continue
c increment timestep
time = time+dtimestep
c After all collisions in the timestep are done, propagate to end of
c the timestep.
call cascstep(acttime,time-acttime)
c in case of potential interaction, do MD propagation step
if (eos.ne.0) then
c set initial conditions for MD propagation-step
do 24 j=1,npart
r0(j) = r0_t(j)
rx(j) = rx_t(j)
ry(j) = ry_t(j)
rz(j) = rz_t(j)
24 continue
c now molecular dynamics trajectories
call proprk(time,dtimestep)
end if ! (eos.ne.0)
c perform output if desired
cdh if(mod(steps,outsteps).eq.0.and.steps.lt.nsteps)then
cdh if(CTOption(28).eq.2)call spectrans(otime)
cdh call file14out(steps)
cdh endif ! output handling
20 continue ! time step loop
c
acttime=time
c optional decay of all unstable particles before final output
c DANGER: pauli-blocked decays are not performed !!!
if(CTOption(18).eq.0) then
c no do-loop is used because npart changes in loop-structure
i=0
nct=0
actcol=0
c disable Pauli-Blocker for final decays
CTOption(10)=1
c decay loop structure starts here
40 continue
i=i+1
c is particle unstable
if(dectime(i).lt.1.d30) then
41 continue
isstable = .false.
do 44 stidx=1,nstable
if (ityp(i).eq.stabvec(stidx)) then
cdh write (6,*) 'no decay of particle ',ityp(i)
isstable = .true.
endif
44 enddo
if (.not.isstable) then
c perform decay
call scatter(i,0,0.d0,fmass(i),xdummy)
c backtracing if decay-product is unstable itself
if(dectime(i).lt.1.d30) goto 41
endif
endif
c check next particle
if(i.lt.npart) goto 40
endif ! final decay
c final output
cdh if(CTOption(28).eq.2)call spectrans(otime)
call file13out(nsteps)
call file14out(nsteps)
cdh call file16out
cdh call osc_event
cdh call osc99_event(1)
cdh call osc99_eoe
c print *,"noc",noc,bdist,ctag,bimp
mp=mp+npart
if(iflbmax.eq.0.and.ctag.eq.0)then
cdh write(*,*)'(W) No collision in event ',event
noc=noc+1
if(noc.ge.1000 .and. mod(noc,1000) .eq. 0)then
c RU Mi 23. Jun 08:47:32 CEST 2021
c switch off printout of number of collision counter warning:
c print *,'no collision problem in UrQMD'
c~ stop
endif
if (noc.ge.50000) then
print *,'UrQMD terminating without collision !? ...'
print *,' ebeam=', ebeam
print *,' projectile mass=', Ap
print *,' projectile charge=', Zp
print *,' target mass=', At
print *, 'target charge=', Zt
write(*,*) ' iterations =', noc
c RU Mi 23. Jun 08:53:43 CEST 2021
c MAYBE do not quit: just assume this is very rate and resembles an elastic FS
call exit(333) ! think of a better way to hand over the error
! to C++
c return ! RU !!!!! back to C8
c RU Mi 23. Jun 08:54:25 CEST 2021
endif
goto 1
endif
c RU Mi 23. Jun 08:37:22 CEST 2021
c switch off printout of number of collision counter:
c write(*,*) 'iterations =', noc
c end of event loop
cdh 10 continue
cdh write(6,*)'no. of collisions = ',mc/dble(nevents), ' (per event)'
cdh write(6,*)'final particles = ',mp/dble(nevents), ' (per event)'
cdh write(6,*)'empty events : ', noc, ' = ',
cdh + noc*1d2/dble(nevents), '%'
return
end