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
