TAmpt/AMPT/zpc.f

   1 c.................... zpc.f
   2 c       PROGRAM ZPC
   3       SUBROUTINE ZPCMN
   4 c       Version: 1.0.1
   5 c       Author: Bin Zhang
   6 c       (suggestions, problems -> bzhang@nt1.phys.columbia.edu)
   7 cms
   8 cms     dlw & gsfs Comments our writing of output files
   9 cms
  10         implicit double precision (a-h, o-z)
  11 clin-4/20/01        PARAMETER (NMAXGL = 16000)
  12         parameter (MAXPTN=400001)
  13         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
  14 cc      SAVE /para3/
  15         SAVE
  16 c
  17 c       loop over events
  18         do 1000 i = 1, nevnt
  19            ievt = i
  20 c       generation of the initial condition for one event
  21            call inievt
  22 c      loop over many runs of the same event
  23            do 2000 j = 1, nsbrun
  24               isbrun = j
  25 c       initialization for one run of an event
  26               call inirun
  27 clin-4/2008 not used:
  28 c             CALL HJAN1A
  29  3000         continue
  30 c       do one collision
  31               call zpcrun(*4000)
  32               call zpca1
  33               goto 3000
  34  4000         continue
  35               call zpca2
  36  2000      continue
  37  1000   continue
  38         call zpcou
  39 clin-5/2009 ctest off
  40 c     5/17/01 calculate v2 for parton already frozen out:
  41 c        call flowp(3)
  42 c.....to get average values for different strings
  43         CALL zpstrg
  44         RETURN
  45         end
  46
  47 ******************************************************************************
  48 ******************************************************************************
  49
  50         block data zpcbdt
  51 c       set initial values in block data
  52
  53         implicit double precision (a-h, o-z)
  54         parameter (MAXPTN=400001)
  55         PARAMETER (MAXSTR=150001)
  56         common /para1/ mul
  57 cc      SAVE /para1/
  58         common /para2/ xmp, xmu, alpha, rscut2, cutof2
  59 cc      SAVE /para2/
  60         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
  61 cc      SAVE /para3/
  62         common /para4/ iftflg, ireflg, igeflg, ibstfg
  63 cc      SAVE /para4/
  64         common /para5/ iconfg, iordsc
  65 cc      SAVE /para5/
  66         common /para6/ centy
  67 cc      SAVE /para6/
  68 clin-6/2009 nsmbbbar and nsmmeson respectively give the total number of
  69 c     baryons/anti-baryons and mesons for each event:
  70 c        common /para7/ ioscar
  71         common /para7/ ioscar,nsmbbbar,nsmmeson
  72 cc      SAVE /para7/
  73         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
  74      &       PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
  75      &       XMASS0(MAXPTN), ITYP0(MAXPTN)
  76 cc      SAVE /prec1/
  77         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
  78      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
  79      &       xmass(MAXPTN), ityp(MAXPTN)
  80 cc      SAVE /prec2/
  81         common /prec3/gxs(MAXPTN),gys(MAXPTN),gzs(MAXPTN),fts(MAXPTN),
  82      &     pxs(MAXPTN), pys(MAXPTN), pzs(MAXPTN), es(MAXPTN),
  83      &     xmasss(MAXPTN), ityps(MAXPTN)
  84 cc      SAVE /prec3/
  85         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
  86 cc      SAVE /prec4/
  87         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
  88 cc      SAVE /prec5/
  89         common /prec6/ etas(MAXPTN), raps(MAXPTN), taus(MAXPTN)
  90 cc      SAVE /prec6/
  91         common /aurec1/ jxa, jya, jza
  92 cc      SAVE /aurec1/
  93         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
  94 cc      SAVE /aurec2/
  95         common /ilist1/
  96      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
  97      &     ictype, icsta(MAXPTN),
  98      &     nic(MAXPTN), icels(MAXPTN)
  99 cc      SAVE /ilist1/
 100         common /ilist2/ icell, icel(10,10,10)
 101 cc      SAVE /ilist2/
 102         common /ilist3/ size1, size2, size3, v1, v2, v3, size
 103 cc      SAVE /ilist3/
 104         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
 105 cc      SAVE /ilist4/
 106 c     6/07/02 initialize in ftime to expedite compiling:
 107 c        common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
 108 cc      SAVE /ilist5/
 109         common /ilist6/ t, iopern, icolln
 110 cc      SAVE /ilist6/
 111         COMMON /ilist7/ LSTRG0(MAXPTN), LPART0(MAXPTN)
 112 cc      SAVE /ilist7/
 113         COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
 114 cc      SAVE /ilist8/
 115         common /rndm1/ number
 116 cc      SAVE /rndm1/
 117         common /rndm2/ iff
 118 cc      SAVE /rndm2/
 119         common /rndm3/ iseedp
 120 cc      SAVE /rndm3/
 121         common /ana1/ ts(12)
 122 cc      SAVE /ana1/
 123         common /ana2/
 124      &     det(12), dn(12), detdy(12), detdn(12), dndy(12),
 125      &     det1(12), dn1(12), detdy1(12), detdn1(12), dndy1(12),
 126      &     det2(12), dn2(12), detdy2(12), detdn2(12), dndy2(12)
 127 cc      SAVE /ana2/
 128         common /ana3/ em(4, 4, 12)
 129 cc      SAVE /ana3/
 130         common /ana4/ fdetdy(24), fdndy(24), fdndpt(12)
 131 cc      SAVE /ana4/
 132         SAVE
 133         data centy/0d0/
 134 c     6/07/02 initialize in ftime to expedite compiling:
 135 c        data (ct(i), i = 1, MAXPTN)/MAXPTN*0d0/
 136 c        data (ot(i), i = 1, MAXPTN)/MAXPTN*0d0/
 137 c        data tlarge/1000000.d0/
 138         data number/0/
 139         data ts/0.11d0, 0.12d0, 0.15d0, 0.2d0, 0.3d0, 0.4d0, 0.6d0,
 140      &     0.8d0, 1d0, 2d0, 4d0, 6d0/
 141 c
 142         end
 143
 144 ******************************************************************************
 145 ******************************************************************************
 146
 147         subroutine inizpc
 148
 149         implicit double precision (a-h, o-z)
 150         SAVE
 151
 152         call readpa
 153
 154         call inipar
 155
 156         call inian1
 157
 158         return
 159         end
 160
 161         subroutine readpa
 162
 163         implicit double precision (a-h, o-z)
 164
 165 cc        external ran1
 166
 167         character*50 str
 168
 169         common /para2/ xmp, xmu, alpha, rscut2, cutof2
 170 cc      SAVE /para2/
 171         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
 172 cc      SAVE /para3/
 173         common /para4/ iftflg, ireflg, igeflg, ibstfg
 174 cc      SAVE /para4/
 175         common /para5/ iconfg, iordsc
 176 cc      SAVE /para5/
 177         common /para7/ ioscar,nsmbbbar,nsmmeson
 178 cc      SAVE /para7/
 179         common /ilist3/ size1, size2, size3, v1, v2, v3, size
 180 cc      SAVE /ilist3/
 181         common /rndm1/ number
 182 cc      SAVE /rndm1/
 183         common /rndm2/ iff
 184 cc      SAVE /rndm2/
 185         common /rndm3/ iseedp
 186 cc      SAVE /rndm3/
 187         SAVE
 188
 189         str=str
 190         iseed=iseedp
 191 c       this is the initialization file containing the initial values of
 192 c          the parameters
 193 cbz1/31/99
 194 c        open (5, file = 'zpc.ini', status = 'unknown')
 195 cbz1/31/99end
 196
 197 c       this is the final data file containing general info about the cascade
 198 cbz1/31/99
 199 c        open (6, file = 'zpc.res', status = 'unknown')
 200 cms     open (25, file = 'ana/zpc.res', status = 'unknown')
 201 cbz1/31/99end
 202
 203 c       this is the input file containing initial particle records
 204 cbz1/25/99
 205 c        open (7, file = 'zpc.inp', status = 'unknown')
 206 cbz1/25/99end
 207
 208 c       this gives the optional OSCAR standard output
 209 cbz1/31/99
 210 c        open (8, file = 'zpc.oscar', status = 'unknown')
 211         if(ioscar.eq.1) then
 212 cms        open (26, file = 'ana/parton.oscar', status = 'unknown')
 213 cms        open (19, file = 'ana/hadron.oscar', status = 'unknown')
 214         endif
 215 cbz1/31/99end
 216
 217 c     2/11/03 combine zpc initialization into ampt.ini:
 218 c        open (29, file = 'zpc.ini', status = 'unknown')
 219 c        read (29, *) str, xmp
 220         xmp=0d0
 221 c        read (29, *) str, xmu
 222 c        read (29, *) str, alpha
 223         cutof2 = 4.5d0 * (alpha / xmu) ** 2
 224 c        read (29, *) str, rscut2
 225         rscut2=0.01d0
 226 c        read (29, *) str, nsevt
 227         nsevt=1
 228 c        read (29, *) str, nevnt
 229         nevnt=1
 230 c        read (29, *) str, nsbrun
 231         nsbrun=1
 232 c        read (29, *) str, iftflg
 233         iftflg=0
 234 c        read (29, *) str, ireflg
 235         ireflg=1
 236 cbz1/31/99
 237         IF (ireflg .EQ. 0) THEN
 238 cms        OPEN (27, FILE = 'zpc.inp', STATUS = 'UNKNOWN')
 239         END IF
 240 cbz1/31/99end
 241 c        read (29, *) str, igeflg
 242         igeflg=0
 243 c        read (29, *) str, ibstfg
 244         ibstfg=0
 245 c        read (29, *) str, iconfg
 246         iconfg=1
 247 c        read (29, *) str, iordsc
 248         iordsc=11
 249 c        read (29, *) str, ioscar
 250 c        read (29, *) str, v1, v2, v3
 251         v1=0.2d0
 252         v2=0.2d0
 253         v3=0.2d0
 254 c        read (29, *) str, size1, size2, size3
 255         size1=1.5d0
 256         size2=1.5d0
 257         size3=0.7d0
 258         if (size1 .eq. 0d0 .or. size2 .eq. 0d0 .or.
 259      &     size3 .eq. 0d0) then
 260            if (size1 .ne. 0d0 .or. size2 .ne. 0d0 .or. size3 .ne. 0d0
 261      &        .or. v1 .ne. 0d0 .or. v2 .ne. 0d0 .or. v3 .ne. 0d0) then
 262               print *, 'to get rid of space division:'
 263               print *, 'set all sizes and vs to 0'
 264               stop 'chker'
 265            end if
 266         end if
 267         size = min(size1, size2, size3)
 268 c        read (29, *) str, iff
 269         iff=-1
 270 c        read (29, *) str, iseed
 271
 272 c     10/24/02 get rid of argument usage mismatch in ran1():
 273         isedng=-iseed
 274 c        a = ran1(-iseed)
 275         a = ran1(isedng)
 276 c        read (29, *) str, irused
 277         irused=2
 278         do 1001 i = 1, irused - 1
 279 c           a = ran1(2)
 280            iseed2=2
 281            a = ran1(iseed2)
 282  1001   continue
 283 c     10/24/02-end
 284
 285         if (iconfg .eq. 2 .or. iconfg .eq. 3) then
 286            v1 = 0d0
 287            v2 = 0d0
 288         end if
 289
 290         if (iconfg .eq. 4 .or. iconfg .eq. 5) then
 291            v1 = 0d0
 292            v2 = 0d0
 293            v3 = 0d0
 294         end if
 295
 296 c        close(5)
 297
 298         return
 299         end
 300
 301         subroutine inipar
 302
 303         implicit double precision (a-h,o-z)
 304
 305         common /para4/ iftflg, ireflg, igeflg, ibstfg
 306 cc      SAVE /para4/
 307         common /para6/ centy
 308 cc      SAVE /para6/
 309         SAVE
 310
 311         if (ibstfg .ne. 0) then
 312            centy = -6d0
 313         end if
 314
 315         return
 316         end
 317
 318         subroutine inian1
 319
 320         implicit double precision (a-h,o-z)
 321
 322         common /para4/ iftflg, ireflg, igeflg, ibstfg
 323 cc      SAVE /para4/
 324         common /ana1/ ts(12)
 325 cc      SAVE /ana1/
 326         SAVE
 327         if (ibstfg .ne. 0) then
 328            a = cosh(6d0)
 329            do 1001 i = 1, 12
 330               ts(i) = ts(i) * a
 331  1001      continue
 332         end if
 333
 334         return
 335         end
 336
 337 ******************************************************************************
 338
 339         subroutine inievt
 340
 341         implicit double precision (a-h, o-z)
 342
 343         COMMON /para1/ mul
 344 cc      SAVE /para1/
 345         common /para4/ iftflg, ireflg, igeflg, ibstfg
 346 cc      SAVE /para4/
 347         SAVE
 348
 349 cbz1/25/99
 350 c        mul = 0
 351 cbz1/25/99
 352         if (ireflg .eq. 0) call readi
 353         if (igeflg .ne. 0) call genei
 354         if (ibstfg .ne. 0) call boosti
 355
 356         return
 357         end
 358
 359         subroutine readi
 360
 361         implicit double precision (a-h, o-z)
 362         parameter (MAXPTN=400001)
 363         double precision field(9)
 364         common /para1/ mul
 365 cc      SAVE /para1/
 366         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
 367 cc      SAVE /para3/
 368         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
 369      &       PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
 370      &       XMASS0(MAXPTN), ITYP0(MAXPTN)
 371 cc      SAVE /prec1/
 372         SAVE
 373         do 1001 i = 1, MAXPTN
 374            if (ievt .ne. 1 .and. i .eq. 1) then
 375               ityp0(i) = ntyp
 376               gx0(1) = field(1)
 377               gy0(1) = field(2)
 378               gz0(1) = field(3)
 379               ft0(1) = field(4)
 380               px0(1) = field(5)
 381               py0(1) = field(6)
 382               pz0(1) = field(7)
 383               e0(1) = field(8)
 384               xmass0(i) = field(9)
 385               mul = 1
 386            else
 387  900              read (27, *, end = 1000) neve, ntyp, field
 388               if (neve .lt. nsevt) goto 900
 389               if (neve .gt.
 390      &           nsevt + ievt - 1) goto 1000
 391               ityp0(i) = ntyp
 392               gx0(i) = field(1)
 393               gy0(i) = field(2)
 394               gz0(i) = field(3)
 395               ft0(i) = field(4)
 396               px0(i) = field(5)
 397               py0(i) = field(6)
 398               pz0(i) = field(7)
 399               e0(i) = field(8)
 400               xmass0(i) = field(9)
 401               mul = mul + 1
 402            end if
 403  1001   continue
 404
 405  1000        continue
 406
 407         return
 408         end
 409
 410         subroutine genei
 411
 412         implicit double precision (a-h, o-z)
 413         parameter (MAXPTN=400001)
 414         common /para1/ mul
 415 cc      SAVE /para1/
 416         common /para2/ xmp, xmu, alpha, rscut2, cutof2
 417 cc      SAVE /para2/
 418         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
 419 cc      SAVE /para3/
 420         common /para5/ iconfg, iordsc
 421 cc      SAVE /para5/
 422         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
 423      &       PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
 424      &       XMASS0(MAXPTN), ITYP0(MAXPTN)
 425 cc      SAVE /prec1/
 426         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
 427 cc      SAVE /prec5/
 428         common /lora/ enenew, pxnew, pynew, pznew
 429 cc      SAVE /lora/
 430         common /rndm3/ iseedp
 431 cc      SAVE /rndm3/
 432         SAVE
 433
 434 cc        external ran1
 435
 436         iseed=iseedp
 437         incmul = 4000
 438         temp = 0.5d0
 439         etamin = -5d0
 440         etamax = 5d0
 441         r0 = 5d0
 442         tau0 = 0.1d0
 443         deta = etamax - etamin
 444
 445         do 1001 i = mul + 1, mul + incmul
 446            ityp0(i) = 21
 447            xmass0(i) = xmp
 448            call energya(e, temp)
 449            call momntm(px, py, pz, e)
 450 c     7/20/01:
 451 c           e = sqrt(e ** 2 + xmp ** 2)
 452            e = dsqrt(e ** 2 + xmp ** 2)
 453            if (iconfg .le. 3) then
 454               eta(i) = etamin + deta * ran1(iseed)
 455               bex = 0d0
 456               bey = 0d0
 457               bez = -tanh(eta(i))
 458               call lorenza(e, px, py, pz, bex, bey, bez)
 459               px0(i) = pxnew
 460               py0(i) = pynew
 461               pz0(i) = pznew
 462               e0(i) = enenew
 463            else
 464               px0(i) = px
 465               py0(i) = py
 466               pz0(i) = pz
 467               e0(i) = e
 468            end if
 469  1001   continue
 470
 471         do 1002 i = mul + 1, mul + incmul
 472            if (iconfg .le. 3) then
 473               gz0(i) = tau0 * sinh(eta(i))
 474               ft0(i) = tau0 * cosh(eta(i))
 475               if (iconfg .eq. 1) then
 476                  call posit1(x, y, r0)
 477                  gx0(i) = x + px0(i) * ft0(i)/e0(i)
 478                  gy0(i) = y + py0(i) * ft0(i)/e0(i)
 479               else if (iconfg .eq. 2 .or. iconfg .eq. 3) then
 480                  call posit2(x, y)
 481                  gx0(i) = x
 482                  gy0(i) = y
 483               end if
 484            else
 485               ft0(i) = 0d0
 486               call posit3(x, y, z)
 487               gx0(i) = x
 488               gy0(i) = y
 489               gz0(i) = z
 490            end if
 491  1002   continue
 492
 493         mul = mul + incmul
 494
 495 c       check if it's necessary to adjust array size 'adarr'
 496             if (mul .ge. MAXPTN .or. mul .eq. 0) then
 497            print *, 'event',ievt,'has',mul,'number of gluon',
 498      &          'adjusting counting is necessary'
 499            stop 'adarr'
 500         end if
 501
 502         return
 503         end
 504
 505         subroutine posit1(x, y, r0)
 506
 507         implicit double precision (a-h, o-z)
 508
 509 cc        external ran1
 510         common /rndm3/ iseedp
 511 cc      SAVE /rndm3/
 512         SAVE
 513
 514         iseed=iseedp
 515  10        x = 2d0 * ran1(iseed) - 1d0
 516         y = 2d0 * ran1(iseed) - 1d0
 517         if (x ** 2 + y ** 2 .gt. 1d0) goto 10
 518         x = x * r0
 519         y = y * r0
 520
 521         return
 522         end
 523
 524         subroutine posit2(x, y)
 525
 526         implicit double precision (a-h, o-z)
 527
 528 c        external ran1
 529
 530         common /ilist3/ size1, size2, size3, v1, v2, v3, size
 531 cc      SAVE /ilist3/
 532         common /rndm3/ iseedp
 533 cc      SAVE /rndm3/
 534         SAVE
 535         iseed=iseedp
 536          x = 2d0 * ran1(iseed) - 1d0
 537         y = 2d0 * ran1(iseed) - 1d0
 538         x = x * 5d0 * size1
 539         y = y * 5d0 * size2
 540
 541         return
 542         end
 543
 544         subroutine posit3(x, y, z)
 545
 546         implicit double precision (a-h, o-z)
 547
 548 cc        external ran1
 549
 550         common /ilist3/ size1, size2, size3, v1, v2, v3, size
 551 cc      SAVE /ilist3/
 552         common /rndm3/ iseedp
 553 cc      SAVE /rndm3/
 554         SAVE
 555
 556         iseed=iseedp
 557          x = 2d0 * ran1(iseed) - 1d0
 558         y = 2d0 * ran1(iseed) - 1d0
 559         z = 2d0 * ran1(iseed) - 1d0
 560         x = x * 5d0 * size1
 561         y = y * 5d0 * size2
 562         z = z * 5d0 * size3
 563
 564         return
 565         end
 566
 567         subroutine energya(e, temp)
 568
 569 c       to generate the magnitude of the momentum e,
 570 c       knowing the temperature of the local thermal distribution temp
 571
 572         implicit double precision (a-h, o-z)
 573
 574 cc        external ran1
 575
 576         common /para2/ xmp, xmu, alpha, rscut2, cutof2
 577 cc      SAVE /para2/
 578         common /rndm3/ iseedp
 579 cc      SAVE /rndm3/
 580         SAVE
 581
 582         iseed=iseedp
 583  1000        continue
 584
 585         e = ran1(iseed)
 586         e = e * ran1(iseed)
 587         e = e * ran1(iseed)
 588
 589         if (e .le. 0d0) goto 1000
 590         e = - temp * log(e)
 591         if (ran1(iseed) .gt.
 592      &     exp((e - dsqrt(e ** 2 + xmp ** 2))/temp)) then
 593            goto 1000
 594         end if
 595
 596         return
 597         end
 598
 599         subroutine momntm(px, py, pz, e)
 600
 601 c       to generate the 3 components of the momentum px, py, pz,
 602 c       from the magnitude of the momentum e
 603
 604         implicit double precision (a-h,o-z)
 605
 606 cc        external ran1
 607
 608         parameter (pi = 3.14159265358979d0)
 609         common /rndm3/ iseedp
 610 cc      SAVE /rndm3/
 611         SAVE
 612
 613         iseed=iseedp
 614         cost = 2d0 * ran1(iseed) - 1d0
 615 c     7/20/01:
 616 c        sint = sqrt(1d0 - cost ** 2)
 617         sint = dsqrt(1d0 - cost ** 2)
 618         phi = 2d0 * pi * ran1(iseed)
 619
 620         px = e * sint * cos(phi)
 621         py = e * sint * sin(phi)
 622         pz = e * cost
 623
 624         return
 625         end
 626
 627         subroutine boosti
 628
 629         implicit double precision (a-h, o-z)
 630         parameter (MAXPTN=400001)
 631         common /para1/ mul
 632 cc      SAVE /para1/
 633         common /para6/ centy
 634 cc      SAVE /para6/
 635         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
 636      &       PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
 637      &       XMASS0(MAXPTN), ITYP0(MAXPTN)
 638 cc      SAVE /prec1/
 639         common /lora/ enenew, pxnew, pynew, pznew
 640 cc      SAVE /lora/
 641         SAVE
 642
 643         external lorenza
 644
 645         bex = 0d0
 646         bey = 0d0
 647         bez = - tanh(centy)
 648
 649 c       save data for many runs of the same initial condition
 650         do 1001 i = 1, mul
 651            px1 = gx0(i)
 652            py1 = gy0(i)
 653            pz1 = gz0(i)
 654            e1 = ft0(i)
 655            call lorenza(e1, px1, py1, pz1, bex, bey, bez)
 656            gx0(i) = pxnew
 657            gy0(i) = pynew
 658            gz0(i) = pznew
 659            ft0(i) = enenew
 660            px1 = px0(i)
 661            py1 = py0(i)
 662            pz1 = pz0(i)
 663            e1 = e0(i)
 664            call lorenza(e1, px1, py1, pz1, bex, bey, bez)
 665            px0(i) = pxnew
 666            py0(i) = pynew
 667            pz0(i) = pznew
 668            e0(i) = enenew
 669  1001   continue
 670
 671         return
 672         end
 673
 674 ******************************************************************************
 675
 676         subroutine inirun
 677         SAVE
 678
 679 c       sort prec2 according to increasing formation time
 680         call ftime
 681         call inirec
 682         call iilist
 683         call inian2
 684
 685         return
 686         end
 687
 688         subroutine ftime
 689 c       this subroutine generates formation time for the particles
 690 c       indexing ft(i)
 691 c       input e(i)
 692 c       output ft(i), indx(i)
 693
 694         implicit double precision (a-h, o-z)
 695
 696         external ftime1
 697         parameter (MAXPTN=400001)
 698         common /para1/ mul
 699 cc      SAVE /para1/
 700         common /para2/ xmp, xmu, alpha, rscut2, cutof2
 701 cc      SAVE /para2/
 702         common /para4/ iftflg, ireflg, igeflg, ibstfg
 703 cc      SAVE /para4/
 704         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
 705      &       PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
 706      &       XMASS0(MAXPTN), ITYP0(MAXPTN)
 707 cc      SAVE /prec1/
 708         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
 709 cc      SAVE /prec4/
 710         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
 711 cc      SAVE /ilist4/
 712         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
 713 cc      SAVE /ilist5/
 714         common /par1/ formt
 715 cc      SAVE /par1/
 716         common/anim/nevent,isoft,isflag,izpc
 717 cc      SAVE /anim/
 718         common /rndm3/ iseedp
 719 cc      SAVE /rndm3/
 720         SAVE
 721
 722         iseed=iseedp
 723 clin-6/07/02 initialize here to expedite compiling, instead in zpcbdt:
 724         do 1001 i = 1, MAXPTN
 725            ct(i)=0d0
 726            ot(i)=0d0
 727  1001   continue
 728         tlarge=1000000.d0
 729 clin-6/07/02-end
 730
 731         if (iftflg .eq. 0) then
 732 c     5/01/01 different prescription for parton initial formation time:
 733            if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
 734               do 1002 i = 1, mul
 735                  if (ft0(i) .gt. tlarge) ft0(i) = tlarge
 736  1002         continue
 737               goto 150
 738            else
 739 c     5/01/01-end
 740
 741            do 1003 i = 1, MAXPTN
 742               ft0(i) = tlarge
 743  1003      continue
 744            do 1004 i = 1, mul
 745               xmt2 = px0(i) ** 2 + py0(i) ** 2 + xmp ** 2
 746               formt = xmt2 / e0(i)
 747               ft0(i) = ftime1(iseed)
 748               if (ft0(i) .gt. tlarge) ft0(i) = tlarge
 749  1004      continue
 750 c     5/01/01:
 751         endif
 752
 753         end if
 754
 755 c     5/01/01:
 756  150        continue
 757
 758 c        call index1(MAXPTN, mul, ft0, indx)
 759         if (mul .gt. 1) then
 760            call index1(MAXPTN, mul, ft0, indx)
 761         else
 762 clin-7/09/03: need to set value for mul=1:
 763            indx(1)=1
 764         end if
 765 c
 766         return
 767         end
 768
 769         subroutine inirec
 770
 771         implicit double precision (a-h, o-z)
 772 cc        external ran1
 773         parameter (MAXPTN=400001)
 774         common /para1/ mul
 775 cc      SAVE /para1/
 776         common /para4/ iftflg, ireflg, igeflg, ibstfg
 777 cc      SAVE /para4/
 778         common /para5/ iconfg, iordsc
 779 cc      SAVE /para5/
 780         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
 781      &       PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
 782      &       XMASS0(MAXPTN), ITYP0(MAXPTN)
 783 cc      SAVE /prec1/
 784         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
 785      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
 786      &       xmass(MAXPTN), ityp(MAXPTN)
 787 cc      SAVE /prec2/
 788         common /prec3/gxs(MAXPTN),gys(MAXPTN),gzs(MAXPTN),fts(MAXPTN),
 789      &     pxs(MAXPTN), pys(MAXPTN), pzs(MAXPTN), es(MAXPTN),
 790      &     xmasss(MAXPTN), ityps(MAXPTN)
 791 cc      SAVE /prec3/
 792         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
 793 cc      SAVE /prec4/
 794         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
 795 cc      SAVE /prec5/
 796         common /prec6/ etas(MAXPTN), raps(MAXPTN), taus(MAXPTN)
 797 cc      SAVE /prec6/
 798         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
 799 cc      SAVE /ilist4/
 800 cbz1/25/99
 801         COMMON /ilist7/ LSTRG0(MAXPTN), LPART0(MAXPTN)
 802 cc      SAVE /ilist7/
 803         COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
 804 cc      SAVE /ilist8/
 805 cbz1/25/99end
 806         COMMON /smearz/smearp,smearh
 807 cc      SAVE /smearz/
 808         dimension vxp(MAXPTN), vyp(MAXPTN), vzp(MAXPTN)
 809         common /precpa/ vxp0(MAXPTN), vyp0(MAXPTN), vzp0(MAXPTN)
 810 cc      SAVE /precpa/
 811         common/anim/nevent,isoft,isflag,izpc
 812 cc      SAVE /anim/
 813 clin-6/06/02 local parton freezeout:
 814         common /frzprc/
 815      &       gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
 816      &       pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
 817      &       xmfrz(MAXPTN),
 818      &       tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
 819 cc      SAVE /frzprc/
 820         common /rndm3/ iseedp
 821 cc      SAVE /rndm3/
 822         common /para7/ ioscar,nsmbbbar,nsmmeson
 823         COMMON /AREVT/ IAEVT, IARUN, MISS
 824         SAVE
 825         iseed=iseedp
 826 clin-6/06/02 local freezeout initialization:
 827         if(isoft.eq.5) then
 828            itlast=0
 829            call inifrz
 830         endif
 831
 832         do 1001 i = 1, mul
 833 clin-7/09/01 define indx(i) to save time:
 834 c           ityp(i) = ityp0(indx(i))
 835 c           gx(i) = gx0(indx(i))
 836 c           gy(i) = gy0(indx(i))
 837 c           gz(i) = gz0(indx(i))
 838 c           ft(i) = ft0(indx(i))
 839 c           px(i) = px0(indx(i))
 840 c           py(i) = py0(indx(i))
 841 c           pz(i) = pz0(indx(i))
 842 c           e(i) = e0(indx(i))
 843 c           xmass(i) = xmass0(indx(i))
 844 ccbz1/25/99
 845 c           LSTRG1(I) = LSTRG0(INDX(I))
 846 c           LPART1(I) = LPART0(INDX(I))
 847 ccbz1/25/99end
 848            indxi=indx(i)
 849            ityp(i) = ityp0(indxi)
 850            gx(i) = gx0(indxi)
 851            gy(i) = gy0(indxi)
 852            gz(i) = gz0(indxi)
 853            ft(i) = ft0(indxi)
 854            px(i) = px0(indxi)
 855            py(i) = py0(indxi)
 856            pz(i) = pz0(indxi)
 857            e(i) = e0(indxi)
 858            xmass(i) = xmass0(indxi)
 859            LSTRG1(I) = LSTRG0(INDXI)
 860            LPART1(I) = LPART0(INDXI)
 861            vxp(I) = vxp0(INDXI)
 862            vyp(I) = vyp0(INDXI)
 863            vzp(I) = vzp0(INDXI)
 864 clin-7/09/01-end
 865 c
 866 clin-6/06/02 local freezeout initialization:
 867          if(isoft.eq.5) then
 868             idfrz(i)=ityp(i)
 869             gxfrz(i)=gx(i)
 870             gyfrz(i)=gy(i)
 871             gzfrz(i)=gz(i)
 872             ftfrz(i)=ft(i)
 873             pxfrz(i)=px(i)
 874             pyfrz(i)=py(i)
 875             pzfrz(i)=pz(i)
 876             efrz(i)=e(i)
 877             xmfrz(i)=xmass(i)
 878             ifrz(i)=0
 879          endif
 880 clin-6/06/02-end
 881  1001 continue
 882
 883 c       save particle info for fixed time analysis
 884         do 1002 i = 1, mul
 885            ityps(i) = ityp(i)
 886            gxs(i) = gx(i)
 887            gys(i) = gy(i)
 888            gzs(i) = gz(i)
 889            fts(i) = ft(i)
 890            pxs(i) = px(i)
 891            pys(i) = py(i)
 892            pzs(i) = pz(i)
 893            es(i) = e(i)
 894            xmasss(i) = xmass(i)
 895  1002   continue
 896
 897 clin-6/2009
 898 cms     if(isoft.eq.1.and.(ioscar.eq.2.or.ioscar.eq.3))
 899 cms  1       write(92,*) iaevt,miss,mul
 900
 901         do 1003 i = 1, mul
 902            energy = e(i)
 903            vx(i) = px(i) / energy
 904            vy(i) = py(i) / energy
 905            vz(i) = pz(i) / energy
 906            if (iftflg .eq. 0) then
 907               formt = ft(i)
 908 c     7/09/01 propagate partons with parent velocity till formation
 909 c     so that partons in same hadron have 0 distance:
 910 c            gx(i) = gx(i) + vx(i) * formt
 911 c            gy(i) = gy(i) + vy(i) * formt
 912 c            gz(i) = gz(i) + vz(i) * formt
 913             if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
 914                gx(i) = gx(i) + vxp(i) * formt
 915                gy(i) = gy(i) + vyp(i) * formt
 916                gz(i) = gz(i) + vzp(i) * formt
 917             else
 918                gx(i) = gx(i) + vx(i) * formt
 919                gy(i) = gy(i) + vy(i) * formt
 920                gz(i) = gz(i) + vz(i) * formt
 921             endif
 922 c     7/09/01-end
 923 c
 924 c     3/27/00-ctest off no smear z on partons to avoid eta overflow:
 925 c              gz(i) = gz(i)+smearp*(2d0 * ran1(iseed) - 1d0)
 926 c     to give eta=y +- smearp*random:
 927 c              smeary=smearp*(2d0 * ran1(iseed) - 1d0)
 928 c              smearf=dexp(2*smeary)*(1+vz(i))/(1-vz(i)+1.d-8)
 929 c              gz(i) = gz(i)+formt*(smearf-1)/(smearf+1)
 930 c     3/27/00-end
 931            end if
 932
 933 clin-6/2009 write out initial parton information after string melting
 934 c     and after propagating to its format time:
 935            if(ioscar.eq.2.or.ioscar.eq.3) then
 936               if(dmax1(abs(gx(i)),abs(gy(i)),
 937      1             abs(gz(i)),abs(ft(i))).lt.9999) then
 938 cms              write(92,200) ityp(i),px(i),py(i),pz(i),xmass(i),
 939 cms  1                gx(i),gy(i),gz(i),ft(i)
 940               else
 941 cms              write(92,201) ityp(i),px(i),py(i),pz(i),xmass(i),
 942 cms  1                gx(i),gy(i),gz(i),ft(i)
 943               endif
 944            endif
 945 cyy 200      format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2))
 946 cyy 201      format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2))
 947 c
 948  1003   continue
 949
 950         if (iconfg .le. 3) then
 951            do 1004 i = 1, mul
 952               if (ft(i) .le. abs(gz(i))) then
 953                  eta(i) = 1000000.d0
 954               else
 955                  eta(i) = 0.5d0 * log((ft(i) + gz(i)) / (ft(i) - gz(i)))
 956               end if
 957               if (e(i) .le. abs(pz(i))) then
 958                  rap(i) = 1000000.d0
 959               else
 960                  rap(i) = 0.5d0 * log((e(i) + pz(i)) / (e(i) - pz(i)))
 961               end if
 962               tau(i) = ft(i) / cosh(eta(i))
 963  1004      continue
 964
 965            do 1005 i = 1, mul
 966               etas(i) = eta(i)
 967               raps(i) = rap(i)
 968               taus(i) = tau(i)
 969  1005      continue
 970         end if
 971
 972         return
 973         end
 974
 975         subroutine iilist
 976
 977         implicit double precision (a-h, o-z)
 978         parameter (MAXPTN=400001)
 979         common /para1/ mul
 980 cc      SAVE /para1/
 981         common /ilist1/
 982      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
 983      &     ictype, icsta(MAXPTN),
 984      &     nic(MAXPTN), icels(MAXPTN)
 985 cc      SAVE /ilist1/
 986         common /ilist2/ icell, icel(10,10,10)
 987 cc      SAVE /ilist2/
 988         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
 989 cc      SAVE /ilist4/
 990         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
 991 cc      SAVE /ilist5/
 992         common /ilist6/ t, iopern, icolln
 993 cc      SAVE /ilist6/
 994         SAVE
 995
 996         iscat = MAXPTN
 997         jscat = MAXPTN
 998
 999         do 1001 i = 1, mul
1000            next(i) = 0
1001            last(i) = 0
1002            icsta(i) = 0
1003            nic(i) = 0
1004            icels(i) = 0
1005  1001   continue
1006
1007         icell = 0
1008         do 1004 i1 = 1, 10
1009            do 1003 i2 = 1, 10
1010               do 1002 i3 = 1, 10
1011                  icel(i1, i2, i3) = 0
1012  1002         continue
1013  1003      continue
1014  1004   continue
1015
1016         ichkpt = 0
1017         ifmpt = 1
1018
1019         do 1005 i = 1, mul
1020            ct(i) = tlarge
1021            ot(i) = tlarge
1022  1005   continue
1023
1024         iopern = 0
1025         icolln = 0
1026         t = 0.d0
1027
1028         return
1029         end
1030
1031         subroutine inian2
1032
1033         implicit double precision (a-h, o-z)
1034
1035         common /para5/ iconfg, iordsc
1036 cc      SAVE /para5/
1037         common /ana2/
1038      &     det(12), dn(12), detdy(12), detdn(12), dndy(12),
1039      &     det1(12), dn1(12), detdy1(12), detdn1(12), dndy1(12),
1040      &     det2(12), dn2(12), detdy2(12), detdn2(12), dndy2(12)
1041 cc      SAVE /ana2/
1042         SAVE
1043
1044         if (iconfg .le. 3) then
1045            do 1001 i = 1, 12
1046               det(i) = 0d0
1047               dn(i) = 0d0
1048               det1(i) = 0d0
1049               dn1(i) = 0d0
1050               det2(i) = 0d0
1051               dn2(i) = 0d0
1052  1001      continue
1053         end if
1054
1055         return
1056         end
1057
1058 ******************************************************************************
1059 ******************************************************************************
1060
1061         subroutine zpcrun(*)
1062
1063         implicit double precision (a-h, o-z)
1064         parameter (MAXPTN=400001)
1065         parameter (tend1 = 250d0)
1066         parameter (tend2 = 6.1d0)
1067         common /para1/ mul
1068 cc      SAVE /para1/
1069         common /para5/ iconfg, iordsc
1070         common /para7/ ioscar,nsmbbbar,nsmmeson
1071 cc      SAVE /para5/
1072         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1073      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1074      &       xmass(MAXPTN), ityp(MAXPTN)
1075 cc      SAVE /prec2/
1076         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1077 cc      SAVE /prec4/
1078         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1079 cc      SAVE /prec5/
1080         common /ilist1/
1081      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1082      &     ictype, icsta(MAXPTN),
1083      &     nic(MAXPTN), icels(MAXPTN)
1084 cc      SAVE /ilist1/
1085         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
1086 cc      SAVE /ilist4/
1087         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1088 cc      SAVE /ilist5/
1089         common /ilist6/ t, iopern, icolln
1090 cc      SAVE /ilist6/
1091         common /ana1/ ts(12)
1092 cc      SAVE /ana1/
1093         common/anim/nevent,isoft,isflag,izpc
1094 cc      SAVE /anim/
1095         COMMON /AREVT/ IAEVT, IARUN, MISS
1096         SAVE
1097
1098 c       save last collision info
1099         if (mod(ictype, 2) .eq. 0) then
1100            call savrec(iscat)
1101            call savrec(jscat)
1102         end if
1103
1104 c1      get operation type
1105         call getict(t1)
1106 c2      check freezeout condition
1107         if (iconfg .eq. 1 .and. t1 .gt. tlarge / 2d0) return 1
1108         if (iconfg .eq. 2 .or. iconfg .eq. 3) then
1109            if (t1 .gt. tend1) return 1
1110 c           if (ichkpt .eq. mul) then
1111 c              ii = 0
1112 c              do i = 1, mul
1113 c                 gztemp = gz(i) + vz(i) * (t1 - ft(i))
1114 c                 if (sqrt(t1 ** 2 - gztemp ** 2) .lt. tend) then
1115 c                    ii = 1
1116 c                    goto 1000
1117 c                 end if
1118 c              end do
1119 c 1000              continue
1120 c              if (ii .eq. 0) return 1
1121 c           end if
1122         end if
1123         if (iconfg .eq. 4 .or. iconfg .eq. 5) then
1124            if (t1 .gt. tend2) return 1
1125         end if
1126
1127 clin-6/06/02 local freezeout for string melting,
1128 c     decide what partons have frozen out at time t1:
1129       if(isoft.eq.5) then
1130          call local(t1)
1131       endif
1132
1133 c3      update iopern, t
1134
1135         iopern = iopern + 1
1136         t = t1
1137         if (mod(ictype, 2) .eq. 0) then
1138            icolln = icolln + 1
1139
1140 c     4/18/01-ctest off
1141 c           write (2006, 1233) 'iscat=', iscat, 'jscat=', jscat,
1142 c           write (2006, *) 'iscat=', iscat, ' jscat=', jscat,
1143 c     1 ityp(iscat), ityp(jscat)
1144 c           write (2006, 1233) 'iscat=', max(indx(iscat), indx(jscat)),
1145 c     &        'jscat=', min(indx(iscat), indx(jscat))
1146
1147 c           write (2006, 1234) ' icolln=', icolln, 't=', t
1148
1149 c 1233           format (a10, i10, a10, i10)
1150 c 1234           format (a15, i10, a5, f23.17, a5, f23.17)
1151         end if
1152
1153 c4.1    deal with formation
1154         if (iconfg .eq. 1
1155      &     .or. iconfg .eq. 2
1156      &     .or. iconfg .eq. 4) then
1157            if (ictype .eq. 1 .or. ictype .eq. 2 .or.
1158      &        ictype .eq. 5 .or. ictype .eq. 6) then
1159               call celasn
1160            end if
1161         end if
1162
1163 c4.2    deal with collisions
1164
1165         if (ictype .ne. 1) then
1166
1167            iscat0 = iscat
1168            jscat0 = jscat
1169
1170 c        iscat is the larger one so that if it's a wall collision,
1171 c       it's still ok
1172            iscat = max0(iscat0, jscat0)
1173            jscat = min0(iscat0, jscat0)
1174
1175 ctest off check icsta(i): 0 with f77 compiler
1176 c        write(9,*) 'BB:ictype,t1,iscat,jscat,icsta(i)=',
1177 c     1 ictype,t1,iscat,jscat,icsta(iscat)
1178
1179 c       check collision time table error 'tterr'
1180 clin-4/2008 to avoid out-of-bound error in next():
1181 c           if (jscat .ne. 0 .and. next(jscat) .ne. iscat)
1182 c     &        then
1183 c              print *, 'iscat=', iscat, 'jscat=', jscat,
1184 c     &             'next(', jscat, ')=', next(jscat)
1185 c
1186 c              if (ct(iscat) .lt. tlarge / 2d0) stop 'tterr'
1187 c              if (ct(jscat) .lt. tlarge / 2d0) stop 'tterr'
1188 c           end if
1189            if (jscat .ne. 0) then
1190               if(next(jscat) .ne. iscat) then
1191                  print *, 'iscat=', iscat, 'jscat=', jscat,
1192      &                'next(', jscat, ')=', next(jscat)
1193                  if (ct(iscat) .lt. tlarge / 2d0) stop 'tterr'
1194                  if (ct(jscat) .lt. tlarge / 2d0) stop 'tterr'
1195               endif
1196            end if
1197 clin-4/2008-end
1198
1199 c4.2.1     collisions with wall
1200
1201 c     8/19/02 avoid actual argument in common blocks of cellre:
1202          niscat=iscat
1203          njscat=jscat
1204 c           if (icsta(iscat) .ne. 0) call cellre(iscat, t)
1205 c           if (jscat .ne. 0) then
1206 c              if (icsta(jscat) .ne. 0) call cellre(jscat, t)
1207 c           end if
1208            if (icsta(iscat) .ne. 0) call cellre(niscat, t)
1209            if (jscat .ne. 0) then
1210               if (icsta(jscat) .ne. 0) call cellre(njscat, t)
1211            end if
1212
1213 c4.2.2     collision between particles
1214
1215 clin-6/2009 write out info for each collision:
1216 c           if (mod(ictype, 2) .eq. 0) call scat(t, iscat, jscat)
1217            if (mod(ictype, 2) .eq. 0) then
1218               if(ioscar.eq.3) then
1219 cms              write(95,*) 'event,miss,iscat,jscat=',iaevt,miss,iscat,jscat
1220                  if(dmax1(abs(gx(iscat)),abs(gy(iscat)),
1221      1                abs(gz(iscat)),abs(ft(iscat)),abs(gx(jscat)),
1222      2                abs(gy(jscat)),abs(gz(jscat)),abs(ft(jscat)))
1223      3                .lt.9999) then
1224 cms                 write(95,200) ityp(iscat),px(iscat),py(iscat),
1225 cms  1                   pz(iscat),xmass(iscat),gx(iscat),gy(iscat),
1226 cms  2                   gz(iscat),ft(iscat)
1227 cms                 write(95,200) ityp(jscat),px(jscat),py(jscat),
1228 cms  1                   pz(jscat),xmass(jscat),gx(jscat),gy(jscat),
1229 cms  2                   gz(jscat),ft(jscat)
1230                  else
1231 cms                 write(95,201) ityp(iscat),px(iscat),py(iscat),
1232 cms  1                   pz(iscat),xmass(iscat),gx(iscat),gy(iscat),
1233 cms  2                   gz(iscat),ft(iscat)
1234 cms                 write(95,201) ityp(jscat),px(jscat),py(jscat),
1235 cms  1                   pz(jscat),xmass(jscat),gx(jscat),gy(jscat),
1236 cms  2                   gz(jscat),ft(jscat)
1237                  endif
1238               endif
1239 c
1240               call scat(t, iscat, jscat)
1241 c
1242               if(ioscar.eq.3) then
1243                  if(dmax1(abs(gx(iscat)),abs(gy(iscat)),
1244      1                abs(gz(iscat)),abs(ft(iscat)),abs(gx(jscat)),
1245      2                abs(gy(jscat)),abs(gz(jscat)),abs(ft(jscat)))
1246      3                .lt.9999) then
1247 cms                 write(95,200) ityp(iscat),px(iscat),py(iscat),
1248 cms  1                   pz(iscat),xmass(iscat),gx(iscat),gy(iscat),
1249 cms  2                   gz(iscat),ft(iscat)
1250 cms                 write(95,200) ityp(jscat),px(jscat),py(jscat),
1251 cms  1                   pz(jscat),xmass(jscat),gx(jscat),gy(jscat),
1252 cms  2                   gz(jscat),ft(jscat)
1253                  else
1254 cms                 write(95,201) ityp(iscat),px(iscat),py(iscat),
1255 cms  1                   pz(iscat),xmass(iscat),gx(iscat),gy(iscat),
1256 cms  2                   gz(iscat),ft(iscat)
1257 cms                 write(95,201) ityp(jscat),px(jscat),py(jscat),
1258 cms  1                   pz(jscat),xmass(jscat),gx(jscat),gy(jscat),
1259 cms  2                   gz(jscat),ft(jscat)
1260                  endif
1261               endif
1262            endif
1263 cyy           format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2))
1264 cyy           format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2))
1265
1266         end if
1267
1268 c5      update the interaction list
1269         call ulist(t)
1270
1271 c6      update ifmpt. ichkpt
1272 c       old ichkpt and ifmpt are more conveniently used in ulist
1273         if (ifmpt .le. mul) then
1274            if (ictype .ne. 0 .and. ictype .ne. 3
1275      &        .and. ictype .ne. 4) then
1276               ichkpt = ichkpt + 1
1277               ifmpt = ifmpt + 1
1278            end if
1279         end if
1280
1281         return
1282         end
1283
1284         subroutine savrec(i)
1285
1286         implicit double precision (a-h, o-z)
1287         parameter (MAXPTN=400001)
1288         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1289      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1290      &       xmass(MAXPTN), ityp(MAXPTN)
1291 cc      SAVE /prec2/
1292         common /prec3/gxs(MAXPTN),gys(MAXPTN),gzs(MAXPTN),fts(MAXPTN),
1293      &     pxs(MAXPTN), pys(MAXPTN), pzs(MAXPTN), es(MAXPTN),
1294      &     xmasss(MAXPTN), ityps(MAXPTN)
1295 cc      SAVE /prec3/
1296         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1297 cc      SAVE /prec5/
1298         common /prec6/ etas(MAXPTN), raps(MAXPTN), taus(MAXPTN)
1299 cc      SAVE /prec6/
1300         SAVE
1301
1302         ityps(i) = ityp(i)
1303         gxs(i) = gx(i)
1304         gys(i) = gy(i)
1305         gzs(i) = gz(i)
1306         fts(i) = ft(i)
1307         pxs(i) = px(i)
1308         pys(i) = py(i)
1309         pzs(i) = pz(i)
1310         es(i) = e(i)
1311         xmasss(i) = xmass(i)
1312         etas(i) = eta(i)
1313         raps(i) = rap(i)
1314         taus(i) = tau(i)
1315
1316         return
1317         end
1318
1319         subroutine getict(t1)
1320         implicit double precision (a-h, o-z)
1321         parameter (MAXPTN=400001)
1322         common /para1/ mul
1323 cc      SAVE /para1/
1324         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1325      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1326      &       xmass(MAXPTN), ityp(MAXPTN)
1327 cc      SAVE /prec2/
1328         common /ilist1/
1329      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1330      &     ictype, icsta(MAXPTN),
1331      &     nic(MAXPTN), icels(MAXPTN)
1332 cc      SAVE /ilist1/
1333         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
1334 cc      SAVE /ilist4/
1335         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1336 cc      SAVE /ilist5/
1337         SAVE
1338
1339 c       neglect possibility of 2 collisions at the same time
1340 c0       set initial conditions
1341
1342         t1 = tlarge
1343         iscat = 0
1344         jscat = 0
1345
1346 c1      get next collision between particles
1347         do 1001 i = 1, ichkpt
1348            if (ot(i) .lt. t1) then
1349               t1 = ot(i)
1350               iscat = i
1351            end if
1352  1001   continue
1353         if (iscat .ne. 0) jscat = next(iscat)
1354
1355 c2      get ictype
1356 c     10/30/02 ictype=0:collision; 1:parton formation
1357         if (iscat .ne. 0 .and. jscat .ne. 0) then
1358            if (icsta(iscat) .eq. 0 .and. icsta(jscat) .eq. 0) then
1359               ictype = 0
1360            else
1361               ictype = 4
1362            end if
1363         else if (iscat .ne. 0 .or. jscat .ne. 0) then
1364            ictype = 3
1365         end if
1366 c
1367         if (ifmpt .le. mul) then
1368            if (ft(ifmpt) .lt. t1) then
1369               ictype = 1
1370               t1 = ft(ifmpt)
1371            else if (ft(ifmpt) .eq. t1) then
1372               if (ictype .eq. 0) ictype = 2
1373               if (ictype .eq. 3) ictype = 5
1374               if (ictype .eq. 4) ictype = 6
1375            end if
1376         end if
1377
1378         return
1379         end
1380
1381         subroutine celasn
1382 c       this subroutine is used to assign a cell for a newly formed particle
1383 c       output: nic(MAXPTN) icels(MAXPTN) in the common /ilist1/
1384 c       icell, and icel(10,10,10) in the common /ilist2/
1385
1386         implicit double precision (a-h, o-z)
1387         parameter (MAXPTN=400001)
1388         common /para1/ mul
1389 cc      SAVE /para1/
1390         common /para5/ iconfg, iordsc
1391 cc      SAVE /para5/
1392         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1393      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1394      &       xmass(MAXPTN), ityp(MAXPTN)
1395 cc      SAVE /prec2/
1396         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1397 cc      SAVE /prec4/
1398         common /ilist1/
1399      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1400      &     ictype, icsta(MAXPTN),
1401      &     nic(MAXPTN), icels(MAXPTN)
1402 cc      SAVE /ilist1/
1403         common /ilist2/ icell, icel(10,10,10)
1404 cc      SAVE /ilist2/
1405         common /ilist3/ size1, size2, size3, v1, v2, v3, size
1406 cc      SAVE /ilist3/
1407         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
1408 cc      SAVE /ilist4/
1409         SAVE
1410
1411         external integ
1412
1413         i = ifmpt
1414         tt = ft(i)
1415         td = tt - size
1416         if (iconfg .eq. 1 .and. (size1 .eq. 0d0 .or.
1417      &     size2 .eq. 0d0 .or. size3 .eq. 0d0)) then
1418            i1 = 11
1419            i2 = 11
1420            i3 = 11
1421         else if (iconfg .eq. 4 .or. td .le. 0d0) then
1422            i1 = integ(gx(i) / size1) + 6
1423            i2 = integ(gy(i) / size2) + 6
1424            i3 = integ(gz(i) / size3) + 6
1425            if (integ(gx(i) / size1) .eq. gx(i) / size1 .and.
1426      &        vx(i) .lt. 0d0)
1427      &        i1 = i1 - 1
1428            if (integ(gy(i) / size2) .eq. gy(i) / size2 .and.
1429      &        vy(i) .lt. 0d0)
1430      &        i2 = i2 - 1
1431            if (integ(gz(i) / size3) .eq. gz(i) / size3 .and.
1432      &        vz(i) .lt. 0d0)
1433      &        i3 = i3 - 1
1434         else
1435            i1 = integ(gx(i) / (size1 + v1 * td)) + 6
1436            i2 = integ(gy(i) / (size2 + v2 * td)) + 6
1437            i3 = integ(gz(i) / (size3 + v3 * td)) + 6
1438            if (integ(gx(i) / (size1 + v1 * td)) .eq. gx(i) /
1439      &        (size1 + v1 * td) .and. vx(i) .lt. (i1 - 6) * v1)
1440      &        i1 = i1 - 1
1441            if (integ(gy(i) / (size2 + v2 * td)) .eq. gy(i)/
1442      &        (size2 + v2 * td) .and. vy(i) .lt. (i2 - 6) * v2)
1443      &        i2 = i2 - 1
1444            if (integ(gz(i) / (size3 + v3 * td)) .eq. gz(i)/
1445      &        (size3 + v3 * td) .and. vz(i) .lt. (i3 - 6) * v3)
1446      &        i3 = i3 - 1
1447         end if
1448
1449         if (i1 .le. 0 .or. i1 .ge. 11 .or. i2 .le. 0 .or.
1450      &     i2 .ge. 11 .or. i3 .le. 0 .or. i3 .ge. 11) then
1451            i1 = 11
1452            i2 = 11
1453            i3 = 11
1454         end if
1455
1456         if (i1 .eq. 11) then
1457            j = icell
1458            call newcre(i, j)
1459            icell = j
1460            icels(i) = 111111
1461         else
1462            j = icel(i1, i2, i3)
1463            call newcre(i, j)
1464            icel(i1, i2, i3) = j
1465            icels(i) = i1 * 10000 + i2 * 100 + i3
1466         end if
1467
1468         return
1469         end
1470
1471         integer function integ(x)
1472 c       this function is used to get the largest integer that is smaller than
1473 c       x
1474
1475         implicit double precision (a-h, o-z)
1476         SAVE
1477
1478         if (x .lt. 0d0) then
1479            integ = int(x - 1d0)
1480         else
1481            integ = int( x )
1482         end if
1483
1484         return
1485         end
1486
1487         subroutine cellre(i, t)
1488 c       this subroutine is used for changing the cell of a particle that
1489 c       collide with the wall
1490
1491         implicit double precision (a-h, o-z)
1492         parameter (MAXPTN=400001)
1493         common /para5/ iconfg, iordsc
1494 cc      SAVE /para5/
1495         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1496      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1497      &       xmass(MAXPTN), ityp(MAXPTN)
1498 cc      SAVE /prec2/
1499         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1500 cc      SAVE /prec4/
1501         common /aurec1/ jxa, jya, jza
1502 cc      SAVE /aurec1/
1503         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
1504 cc      SAVE /aurec2/
1505         common /ilist1/
1506      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1507      &     ictype, icsta(MAXPTN),
1508      &     nic(MAXPTN), icels(MAXPTN)
1509 cc      SAVE /ilist1/
1510         common /ilist2/ icell, icel(10,10,10)
1511 cc      SAVE /ilist2/
1512         common /ilist3/ size1, size2, size3, v1, v2, v3, size
1513 cc      SAVE /ilist3/
1514         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
1515 cc      SAVE /ilist4/
1516         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1517 cc      SAVE /ilist5/
1518         SAVE
1519
1520         logical good
1521
1522         external integ
1523
1524 c       this happens before update the /prec2/ common; in contrast with
1525 c       scat which happens after updating the glue common
1526
1527         t0 = t
1528
1529  1000        continue
1530
1531         if (iconfg .eq. 3 .or. iconfg .eq. 5) then
1532            k = mod(icsta(i), 10)
1533
1534            if (k .eq. 1) then
1535               gx(i) = gx(i) - 10d0 * size1
1536               dgxa(i) = dgxa(i) + 10d0 * size1
1537               do 1001 ii = 1, ichkpt
1538                  if (next(ii) .eq. i) then
1539                     dgxa(ii) = dgxa(ii) - 10d0 * size1
1540                  end if
1541  1001         continue
1542            end if
1543            if (k .eq. 2) then
1544               gx(i) = gx(i) + 10d0 * size1
1545               dgxa(i) = dgxa(i) - 10d0 * size1
1546               do 1002 ii = 1, ichkpt
1547                  if (next(ii) .eq. i) then
1548                     dgxa(ii) = dgxa(ii) + 10d0 * size1
1549                  end if
1550  1002         continue
1551            end if
1552            if (k .eq. 3) then
1553               gy(i) = gy(i) - 10d0 * size2
1554               dgya(i) = dgya(i) + 10d0 * size2
1555               do 1003 ii = 1, ichkpt
1556                  if (next(ii) .eq. i) then
1557                     dgya(ii) = dgya(ii) - 10d0 * size2
1558                  end if
1559  1003         continue
1560            end if
1561            if (k .eq. 4) then
1562               gy(i) = gy(i) + 10d0 * size2
1563               dgya(i) = dgya(i) - 10d0 * size2
1564               do 1004 ii = 1, ichkpt
1565                  if (next(ii) .eq. i) then
1566                     dgya(ii) = dgya(ii) + 10d0 * size2
1567                  end if
1568  1004         continue
1569            end if
1570            if (iconfg .eq. 5) then
1571               if (k .eq. 5) then
1572                  gz(i) = gz(i) - 10d0 * size3
1573                  dgza(i) = dgza(i) + 10d0 * size3
1574                  do 1005 ii = 1, ichkpt
1575                     if (next(ii) .eq. i) then
1576                        dgza(ii) = dgza(ii) - 10d0 * size3
1577                     end if
1578  1005            continue
1579               end if
1580               if (k .eq. 6) then
1581                  gz(i) = gz(i) + 10d0 * size3
1582                  dgza(i) = dgza(i) - 10d0 * size3
1583                  do 1006 ii = 1, ichkpt
1584                     if (next(ii) .eq. i) then
1585                        dgza(ii) = dgza(ii) + 10d0 * size3
1586                     end if
1587  1006               continue
1588               end if
1589            end if
1590         else
1591            icels0 = icels(i)
1592
1593            i1 = icels0 / 10000
1594            i2 = (icels0 - i1 * 10000) / 100
1595            i3 = icels0 - i1 * 10000 - i2 * 100
1596
1597 cc       for particle inside the cube
1598            if (i1 .ge. 1 .and. i1 .le. 10
1599      &        .and. i2 .ge. 1 .and. i2 .le. 10
1600      &        .and. i3 .ge. 1 .and. i3 .le. 10) then
1601
1602 c       this assignment takes care of nic(i)=0 automatically
1603               if (icel(i1, i2, i3) .eq. i) icel(i1, i2, i3) = nic(i)
1604
1605 c1      rearrange the old cell
1606
1607               call oldcre(i)
1608
1609 c2      rearrange the new cell
1610
1611               k = mod(icsta(i), 10)
1612
1613 c2.1    particle goes out of the cube
1614               if (iconfg .eq. 1) then
1615                  good = (i1 .eq. 1 .and. k .eq. 2)
1616      &              .or. (i1 .eq. 10 .and. k .eq. 1)
1617      &              .or. (i2 .eq. 1 .and. k .eq. 4)
1618      &              .or. (i2 .eq. 10 .and. k .eq. 3)
1619      &              .or. (i3 .eq. 1 .and. k .eq. 6)
1620      &              .or. (i3 .eq. 10 .and. k .eq. 5)
1621               end if
1622               if (iconfg .eq. 2) then
1623                  good = (i3 .eq. 1 .and. k .eq. 6)
1624      &              .or. (i3 .eq. 10 .and. k .eq. 5)
1625               end if
1626               if (good) then
1627
1628 c                j = icell
1629                  call newcre(i, icell)
1630 c                 icell = j
1631
1632                  icels(i) = 111111
1633
1634 c2.2    particle moves inside the cube
1635               else
1636
1637                  if (k .eq. 1) i1 = i1 + 1
1638                  if (k .eq. 2) i1 = i1 - 1
1639                  if (k .eq. 3) i2 = i2 + 1
1640                  if (k .eq. 4) i2 = i2 - 1
1641                  if (k .eq. 5) i3 = i3 + 1
1642                  if (k .eq. 6) i3 = i3 - 1
1643
1644                  if (iconfg .eq. 2 .or. iconfg .eq. 4) then
1645                     if (i1 .eq. 0) then
1646                        i1 = 10
1647                        gx(i) = gx(i) + 10d0 * size1
1648                     end if
1649                     if (i1 .eq. 11) then
1650                        i1 = 1
1651                        gx(i) = gx(i) - 10d0 * size1
1652                     end if
1653                     if (i2 .eq. 0) then
1654                        i2 = 10
1655                        gy(i) = gy(i) + 10d0 * size2
1656                     end if
1657                     if (i2 .eq. 11) then
1658                        i2 = 1
1659                        gy(i) = gy(i) - 10d0 * size2
1660                     end if
1661                     if (iconfg .eq. 4) then
1662                        if (i3 .eq. 0) then
1663                           i3 = 10
1664                           gz(i) = gz(i) + 10d0 * size3
1665                        end if
1666                        if (i3 .eq. 11) then
1667                           i3 = 1
1668                           gz(i) = gz(i) - 10d0 * size3
1669                        end if
1670                     end if
1671                  end if
1672
1673                  j = icel(i1, i2, i3)
1674
1675                  call newcre(i, j)
1676 c       in case icel changes
1677
1678                  icel(i1 ,i2, i3) = j
1679
1680                  icels(i) = i1 * 10000 + i2 * 100 + i3
1681
1682               end if
1683
1684 cc       for particles outside the cube
1685            else
1686
1687               if (icell .eq. i) icell = nic(i)
1688
1689               call oldcre(i)
1690
1691               k = mod(icsta(i), 10)
1692
1693               ddt = t - ft(i)
1694               dtt = t - size
1695               if (dtt .le. 0d0) then
1696                  i1 = integ((gx(i) + vx(i) * ddt) / size1) + 6
1697                  i2 = integ((gy(i) + vy(i) * ddt) / size2) + 6
1698                  i3 = integ((gz(i) + vz(i) * ddt) / size3) + 6
1699               else
1700                  i1 = integ((gx(i) + vx(i) * ddt) /
1701      &               (size1 + v1 * dtt)) + 6
1702                  i2 = integ((gy(i) + vy(i) * ddt) /
1703      &               (size2 + v2 * dtt)) + 6
1704                  i3 = integ((gz(i) + vz(i) * ddt) /
1705      &               (size3 + v3 * dtt)) + 6
1706               end if
1707
1708
1709               if (k .eq. 1) i1 = 1
1710               if (k .eq. 2) i1 = 10
1711               if (k .eq. 3) i2 = 1
1712               if (k .eq. 4) i2 = 10
1713               if (k .eq. 5) i3 = 1
1714               if (k .eq. 6) i3 = 10
1715
1716               j = icel(i1, i2, i3)
1717               call newcre(i, j)
1718               icel(i1, i2, i3) = j
1719
1720               icels(i) = i1 * 10000 + i2 * 100 + i3
1721
1722            end if
1723         end if
1724
1725         if (next(i) .ne. 0) then
1726            otmp = ot(next(i))
1727            ctmp = ct(next(i))
1728         end if
1729
1730         if (i1 .eq. 11 .and. i2 .eq. 11 .and. i3 .eq. 11) then
1731            call dchout(i, k, t)
1732         else
1733            if (iconfg .eq. 1) then
1734               call dchin1(i, k, i1, i2, i3, t)
1735            else if (iconfg .eq. 2) then
1736               call dchin2(i, k, i1, i2, i3, t)
1737            else if (iconfg .eq. 4) then
1738               call dchin3(i, k, i1, i2, i3, t)
1739            end if
1740         end if
1741
1742         if (icsta(i) / 10 .eq. 11) then
1743            ot(next(i)) = otmp
1744            ct(next(i)) = ctmp
1745            next(next(i)) = i
1746            call wallc(i, i1, i2, i3, t0, tmin1)
1747            if (tmin1 .lt. ct(i)) then
1748               icsta(i) = icsta(i) + 10
1749               t0 = tmin1
1750               goto 1000
1751            end if
1752         end if
1753
1754         return
1755         end
1756
1757         subroutine oldcre(i)
1758 c       this subroutine is used to rearrange the old cell nic when a particle
1759 c       goes out of the cell
1760
1761         implicit double precision (a-h, o-z)
1762         parameter (MAXPTN=400001)
1763         common /ilist1/
1764      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1765      &     ictype, icsta(MAXPTN),
1766      &     nic(MAXPTN), icels(MAXPTN)
1767 cc      SAVE /ilist1/
1768         SAVE
1769
1770         if (nic(i) .eq. 0) return
1771
1772         j = nic(i)
1773
1774         if (nic(j) .eq. i) then
1775            nic(j) = 0
1776            return
1777         end if
1778
1779         do 10 while (nic(j) .ne. i)
1780            j = nic(j)
1781  10        continue
1782
1783         nic(j) = nic(i)
1784
1785         return
1786         end
1787
1788
1789         subroutine newcre(i, k)
1790 c       this subroutine is used to mk rearrange of the new cell a particle
1791 c       enters,
1792 c       input i
1793 c       output nic(i)
1794
1795         implicit double precision (a-h, o-z)
1796         parameter (MAXPTN=400001)
1797         common /ilist1/
1798      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1799      &     ictype, icsta(MAXPTN),
1800      &     nic(MAXPTN), icels(MAXPTN)
1801 cc      SAVE /ilist1/
1802         SAVE
1803
1804         if (k .eq. 0) then
1805            k = i
1806            nic(i) = 0
1807         else if (nic(k) .eq. 0) then
1808            nic(k) = i
1809            nic(i) = k
1810         else
1811            j = k
1812            do 10 while (nic(j) .ne. k)
1813               j = nic(j)
1814  10           continue
1815
1816            nic(j) = i
1817            nic(i) = k
1818
1819         end if
1820
1821         return
1822         end
1823
1824         subroutine scat(t, iscat, jscat)
1825
1826 c       this subroutine is used to calculate the 2 particle scattering
1827
1828         implicit double precision (a-h, o-z)
1829         SAVE
1830
1831         call newpos(t, iscat)
1832         call newpos(t, jscat)
1833         call newmom(t)
1834
1835         return
1836         end
1837
1838         subroutine newpos(t, i)
1839
1840 c       this subroutine is used to calculate the 2 particle scattering
1841 c       get new position
1842
1843         implicit double precision (a-h, o-z)
1844         parameter (MAXPTN=400001)
1845         common /para5/ iconfg, iordsc
1846 cc      SAVE /para5/
1847         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1848      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1849      &       xmass(MAXPTN), ityp(MAXPTN)
1850 cc      SAVE /prec2/
1851         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1852 cc      SAVE /prec4/
1853         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1854 cc      SAVE /prec5/
1855         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1856 cc      SAVE /ilist5/
1857         SAVE
1858
1859         t=t
1860         dt1 = ct(i) - ft(i)
1861
1862         gx(i) = gx(i) + vx(i) * dt1
1863         gy(i) = gy(i) + vy(i) * dt1
1864         gz(i) = gz(i) + vz(i) * dt1
1865         ft(i) = ct(i)
1866
1867         if (iconfg .le. 3) then
1868            if (ft(i) .le. abs(gz(i))) then
1869               eta(i) = 1000000.d0
1870            else
1871               eta(i) = 0.5d0 * log((ft(i) + gz(i)) / (ft(i) - gz(i)))
1872            end if
1873            tau(i) = ft(i) / cosh(eta(i))
1874         end if
1875
1876         return
1877         end
1878
1879         subroutine newmom(t)
1880
1881 c       this subroutine is used to calculate the 2 particle scattering
1882
1883         implicit double precision (a-h, o-z)
1884
1885         parameter (hbarc = 0.197327054d0)
1886         parameter (MAXPTN=400001)
1887         parameter (pi = 3.14159265358979d0)
1888         COMMON /para1/ mul
1889 cc      SAVE /para1/
1890         common /para2/ xmp, xmu, alpha, rscut2, cutof2
1891 cc      SAVE /para2/
1892         common /para5/ iconfg, iordsc
1893 cc      SAVE /para5/
1894 ctrans
1895         common /para6/ centy
1896 cc      SAVE /para6/
1897 ctransend
1898         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1899      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1900      &       xmass(MAXPTN), ityp(MAXPTN)
1901 cc      SAVE /prec2/
1902         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1903 cc      SAVE /prec4/
1904         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1905 cc      SAVE /prec5/
1906         common /aurec1/ jxa, jya, jza
1907 cc      SAVE /aurec1/
1908         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
1909 cc      SAVE /aurec2/
1910         common /ilist1/
1911      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1912      &     ictype, icsta(MAXPTN),
1913      &     nic(MAXPTN), icels(MAXPTN)
1914 cc      SAVE /ilist1/
1915         common /ilist3/ size1, size2, size3, v1, v2, v3, size
1916 cc      SAVE /ilist3/
1917         common /lora/ enenew, pxnew, pynew, pznew
1918 cc      SAVE /lora/
1919         common /cprod/ xn1, xn2, xn3
1920 cc      SAVE /cprod/
1921         common /rndm2/ iff
1922 cc      SAVE /rndm2/
1923         common/anim/nevent,isoft,isflag,izpc
1924 cc      SAVE /anim/
1925         common /frzprc/
1926      &       gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
1927      &       pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
1928      &       xmfrz(MAXPTN),
1929      &       tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
1930 cc      SAVE /frzprc/
1931         SAVE
1932
1933         t=t
1934 clin-6/06/02 no momentum change for partons already frozen out,
1935 c     however, spatial upgrade is needed to ensure overall system freezeout:
1936       if(isoft.eq.5) then
1937          if(ifrz(iscat).eq.1.or.ifrz(jscat).eq.1) then
1938             last(iscat) = jscat
1939             last(jscat) = iscat
1940             return
1941          endif
1942       endif
1943 clin-6/06/02-end
1944
1945 c       iff is used to randomize the interaction to have both attractive and
1946 c        repulsive
1947
1948         iff = - iff
1949
1950         if (iconfg .eq. 2 .or. iconfg .eq. 4) then
1951            icels1 = icels(iscat)
1952            i1 = icels1 / 10000
1953            j1 = (icels1 - i1 * 10000) / 100
1954            icels2 = icels(jscat)
1955            i2 = icels2 / 10000
1956            j2 = (icels2 - i2 * 10000) / 100
1957            if (iconfg .eq. 4) then
1958               k1 = icels1 - i1 * 10000 - j1 * 100
1959               k2 = icels2 - i2 * 10000 - j2 * 100
1960            end if
1961         end if
1962
1963         px1 = px(iscat)
1964         py1 = py(iscat)
1965         pz1 = pz(iscat)
1966         e1 = e(iscat)
1967         x1 = gx(iscat)
1968         y1 = gy(iscat)
1969         z1 = gz(iscat)
1970         t1 = ft(iscat)
1971         px2 = px(jscat)
1972         py2 = py(jscat)
1973         pz2 = pz(jscat)
1974         e2 = e(jscat)
1975
1976         if (iconfg .eq. 1) then
1977            x2 = gx(jscat)
1978            y2 = gy(jscat)
1979            z2 = gz(jscat)
1980         else if (iconfg .eq. 2 .or. iconfg .eq. 4) then
1981            if (i1 - i2 .gt. 5) then
1982               x2 = gx(jscat) + 10d0 * size1
1983            else if (i1 - i2 .lt. -5) then
1984               x2 = gx(jscat) - 10d0 * size1
1985            else
1986               x2 = gx(jscat)
1987            end if
1988            if (j1 - j2 .gt. 5) then
1989               y2 = gy(jscat) + 10d0 * size2
1990            else if (j1 - j2 .lt. -5) then
1991               y2 = gy(jscat) - 10d0 * size2
1992            else
1993               y2 = gy(jscat)
1994            end if
1995            if (iconfg .eq. 4) then
1996               if (k1 - k2 .gt. 5) then
1997                  z2 = gz(jscat) + 10d0 * size3
1998               else if (k1 - k2 .lt. -5) then
1999                  z2 = gz(jscat) - 10d0 * size3
2000               else
2001                  z2 = gz(jscat)
2002               end if
2003            else
2004               z2 = gz(jscat)
2005            end if
2006         else if (iconfg .eq. 3 .or. iconfg .eq. 5) then
2007            x2 = gx(jscat) + dgxa(jscat)
2008            y2 = gy(jscat) + dgya(jscat)
2009            if (iconfg .eq. 5) then
2010               z2 = gz(jscat) + dgza(jscat)
2011            else
2012               z2 = gz(jscat)
2013            end if
2014         end if
2015         t2 = ft(jscat)
2016 ctrans
2017         rts2 = (e1 + e2) ** 2 - (px1 + px2) ** 2 -
2018      &     (py1 + py2) ** 2 - (pz1 + pz2) ** 2
2019 ctransend
2020         bex = (px1 + px2) / (e1 + e2)
2021         bey = (py1 + py2) / (e1 + e2)
2022         bez = (pz1 + pz2) / (e1 + e2)
2023
2024         call lorenza(e1, px1, py1, pz1, bex, bey, bez)
2025 cc      SAVE pxnew, ..., values for later use.
2026         px1 = pxnew
2027         py1 = pynew
2028         pz1 = pznew
2029         e1 = enenew
2030
2031         pp2 = pxnew ** 2 + pynew ** 2 + pznew ** 2
2032         call getht(iscat, jscat, pp2, that)
2033         theta = dacos(that / (2d0 * pp2) + 1d0)
2034         theta = dble(iff) * theta
2035
2036 c       we boost to the cm frame, get rotation axis, and rotate 1 particle
2037 c       momentum
2038
2039         call lorenza(t1, x1, y1, z1, bex, bey, bez)
2040
2041         x1 = pxnew
2042         y1 = pynew
2043         z1 = pznew
2044
2045         call lorenza(t2, x2, y2, z2, bex, bey, bez)
2046
2047         x2 = pxnew
2048         y2 = pynew
2049         z2 = pznew
2050
2051 c       notice now pxnew, ..., are new positions
2052         call cropro(x1-x2, y1-y2, z1-z2, px1, py1, pz1)
2053
2054         call xnormv(xn1, xn2, xn3)
2055
2056 cbz1/29/99
2057 c        call rotate(xn1, xn2, xn3, theta, px1, py1, pz1)
2058         call zprota(xn1, xn2, xn3, theta, px1, py1, pz1)
2059 cbz1/29/99end
2060
2061 c       we invert the momentum to get the other particle's momentum
2062         px2 = -px1
2063         py2 = -py1
2064         pz2 = -pz1
2065 clin-4/13/01: modify in case m1, m2 are different:
2066 c        e2 = e1
2067         e2 = dsqrt(px2**2+py2**2+pz2**2+xmass(jscat)**2)
2068
2069 c       boost the 2 particle 4 momentum back to lab frame
2070         call lorenza(e1, px1, py1, pz1, -bex, -bey, -bez)
2071         px(iscat) = pxnew
2072         py(iscat) = pynew
2073         pz(iscat) = pznew
2074         e(iscat) = enenew
2075         call lorenza(e2, px2, py2, pz2, -bex, -bey, -bez)
2076         px(jscat) = pxnew
2077         py(jscat) = pynew
2078         pz(jscat) = pznew
2079         e(jscat) = enenew
2080
2081         vx(iscat) = px(iscat) / e(iscat)
2082         vy(iscat) = py(iscat) / e(iscat)
2083         vz(iscat) = pz(iscat) / e(iscat)
2084         vx(jscat) = px(jscat) / e(jscat)
2085         vy(jscat) = py(jscat) / e(jscat)
2086         vz(jscat) = pz(jscat) / e(jscat)
2087
2088         last(iscat) = jscat
2089         last(jscat) = iscat
2090
2091         if (iconfg .le. 3) then
2092            if (e(iscat) .le. abs(pz(iscat))) then
2093               rap(iscat) = 1000000.d0
2094            else
2095               rap(iscat) = 0.5d0 * log((e(iscat) + pz(iscat)) /
2096      &           (e(iscat) - pz(iscat)))
2097            end if
2098
2099            if (e(jscat) .le. abs(pz(jscat))) then
2100               rap(jscat) = 1000000.d0
2101            else
2102               rap(jscat) = 0.5d0 * log((e(jscat) + pz(jscat)) /
2103      &           (e(jscat) - pz(jscat)))
2104            end if
2105
2106 ctrans
2107            rap1 = rap(iscat)
2108            rap2 = rap(jscat)
2109
2110            if ((rap1 .lt. centy + 0.5d0 .and.
2111      &        rap1 .gt. centy - 0.5d0)) then
2112 c              write (9, *) sqrt(ft(iscat) ** 2 - gz(iscat) ** 2), rts2
2113            end if
2114            if ((rap2 .lt. centy + 0.5d0 .and.
2115      &        rap2 .gt. centy - 0.5d0)) then
2116 c              write (9, *) sqrt(ft(jscat) ** 2 - gz(jscat) ** 2), rts2
2117            end if
2118 ctransend
2119         end if
2120
2121         return
2122         end
2123
2124         subroutine getht(iscat, jscat, pp2, that)
2125
2126 c       this subroutine is used to get \hat{t} for a particular processes
2127
2128         implicit double precision (a-h, o-z)
2129
2130         parameter (hbarc = 0.197327054d0)
2131         parameter (MAXPTN=400001)
2132         common /para2/ xmp, xmu, alpha, rscut2, cutof2
2133 cc      SAVE /para2/
2134         common/anim/nevent,isoft,isflag,izpc
2135 cc      SAVE /anim/
2136 cc        external ran1
2137         common /rndm3/ iseedp
2138 cc      SAVE /rndm3/
2139         SAVE
2140
2141         iscat=iscat
2142         jscat=jscat
2143         iseed=iseedp
2144         xmu2 = (hbarc * xmu) ** 2
2145         xmp2 = xmp ** 2
2146         xm2 = xmu2 + xmp2
2147         rx=ran1(iseed)
2148         that = xm2*(1d0+1d0/((1d0-xm2/(4d0*pp2+xm2))*rx-1d0))
2149 ctest off isotropic scattering:
2150 c     &     + 1d0/((1d0 - xm2 / (4d0 * pp2 + xm2)) * ran1(2) - 1d0))
2151 c        if(izpc.eq.100) that=-4d0*pp2*ran1(2)
2152         if(izpc.eq.100) that=-4d0*pp2*rx
2153
2154         return
2155         end
2156
2157       subroutine ulist(t)
2158 c     this subroutine is used to update a new collision time list
2159 c       notice this t has been updated
2160
2161         implicit double precision (a-h, o-z)
2162         parameter (MAXPTN=400001)
2163         common /ilist1/
2164      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
2165      &     ictype, icsta(MAXPTN),
2166      &     nic(MAXPTN), icels(MAXPTN)
2167 cc      SAVE /ilist1/
2168         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
2169 cc      SAVE /ilist4/
2170         SAVE
2171
2172         if (ictype .eq. 1 .or. ictype .eq. 2 .or. ictype .eq. 5
2173      &     .or. ictype .eq. 6) then
2174            l = ifmpt
2175            call ulist1(l, t)
2176         end if
2177         if (ictype .ne. 1) then
2178            l = iscat
2179            call ulist1(l, t)
2180            if (jscat .ne. 0) then
2181               l = jscat
2182               call ulist1(l, t)
2183            end if
2184         end if
2185
2186         return
2187         end
2188
2189         subroutine ulist1(l, t)
2190 c       this subroutine is used to update the interaction list when particle
2191 c       l is disturbed.
2192
2193         implicit double precision (a-h, o-z)
2194         parameter (MAXPTN=400001)
2195         common /para5/ iconfg, iordsc
2196 cc      SAVE /para5/
2197         common /ilist1/
2198      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
2199      &     ictype, icsta(MAXPTN),
2200      &     nic(MAXPTN), icels(MAXPTN)
2201 cc      SAVE /ilist1/
2202         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2203 cc      SAVE /ilist5/
2204         SAVE
2205
2206         icels0 = icels(l)
2207         i1 = icels0 / 10000
2208         i2 = (icels0 - i1 * 10000) / 100
2209         i3 = icels0 - i1 * 10000 - i2 * 100
2210 c       save collision info for use when the collision is a collision with wall
2211 c       otherwise wallc will change icsta
2212         k = mod(icsta(l), 10)
2213
2214         call wallc(l, i1, i2, i3, t, tmin1)
2215         tmin = tmin1
2216         nc = 0
2217
2218         if (i1 .eq. 11 .and. i2 .eq. 11 .and. i3 .eq. 11) then
2219            call chkout(l, t, tmin, nc)
2220         else
2221            if (iconfg .eq. 1) then
2222               call chkin1(l, i1, i2, i3, t, tmin, nc)
2223            else if (iconfg .eq. 2) then
2224               call chkin2(l, i1, i2, i3, t, tmin, nc)
2225            else if (iconfg .eq. 4) then
2226               call chkin3(l, i1, i2, i3, t, tmin, nc)
2227            else if (iconfg .eq. 3 .or. iconfg .eq. 5) then
2228               call chkcel(l, i1, i2, i3, t, tmin, nc)
2229            end if
2230         end if
2231
2232         call fixtim(l, t, tmin1, tmin, nc)
2233
2234         return
2235         end
2236
2237         subroutine wallc(i, i1, i2, i3, t, tmin)
2238 c       this subroutine calculates the next time for collision with wall
2239 c       for particle i
2240 c       input particle label i,t
2241 c       output tmin collision time with wall, icsta(i) wall collision
2242 c       information
2243
2244         implicit double precision (a-h, o-z)
2245         parameter (MAXPTN=400001)
2246         common /para5/ iconfg, iordsc
2247 cc      SAVE /para5/
2248         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2249 cc      SAVE /ilist5/
2250         SAVE
2251
2252         tmin = tlarge
2253
2254         if (iconfg .le. 2 .or. iconfg .eq. 4) then
2255 c       if particle is inside the cube
2256            if ((i1 .ge. 1 .and. i1 .le. 10)
2257      &          .or. (i2 .ge. 1 .and. i2 .le. 10)
2258      &          .or. (i3 .ge. 1 .and. i3 .le. 10)) then
2259               call wallc1(i, i1, i2, i3, t, tmin)
2260 c       if particle is outside the cube
2261            else
2262               call wallcb(i, t, tmin)
2263            end if
2264         else if (iconfg .eq. 3 .or. iconfg .eq. 5) then
2265            call wallc2(i, i1, i2, i3, t, tmin)
2266         end if
2267
2268         return
2269         end
2270
2271         subroutine wallc1(i, i1, i2, i3, t, tmin)
2272 c       this subroutine is used to get wall collision time
2273 c       when particle is inside the cube, it sets the icsta at the same time
2274 c       input i,i1,i2,i3,t
2275 c       output tmin, icsta(i)
2276 c       note the icsta is not finally set. we need further judgement in
2277 c       fixtim
2278
2279         implicit double precision (a-h, o-z)
2280         parameter (MAXPTN=400001)
2281         common /para5/ iconfg, iordsc
2282 cc      SAVE /para5/
2283         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
2284      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
2285      &       xmass(MAXPTN), ityp(MAXPTN)
2286 cc      SAVE /prec2/
2287         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
2288 cc      SAVE /prec4/
2289         common /ilist1/
2290      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
2291      &     ictype, icsta(MAXPTN),
2292      &     nic(MAXPTN), icels(MAXPTN)
2293 cc      SAVE /ilist1/
2294         common /ilist3/ size1, size2, size3, v1, v2, v3, size
2295 cc      SAVE /ilist3/
2296         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2297 cc      SAVE /ilist5/
2298         SAVE
2299
2300         x1p = gx(i)
2301         x2p = gy(i)
2302         x3p = gz(i)
2303         tf = ft(i)
2304         v1p = vx(i)
2305         v2p = vy(i)
2306         v3p = vz(i)
2307
2308         if (t .lt. size .and. tf .lt. size) then
2309
2310            if (v1p .gt. 0d0) then
2311               t1 = ((dble(i1) - 5d0) * size1 - x1p) / v1p + tf
2312            else if (v1p .lt. 0d0) then
2313               t1 = ((dble(i1) - 6d0) * size1 - x1p) / v1p + tf
2314            else
2315               t1 = tlarge
2316            end if
2317
2318            if (v2p .gt. 0d0) then
2319               t2 = ((dble(i2) - 5d0) * size2 - x2p) / v2p + tf
2320            else if (v2p .lt. 0d0) then
2321               t2 = ((dble(i2) - 6d0) * size2 - x2p) / v2p + tf
2322            else
2323               t2 = tlarge
2324            end if
2325
2326            if (v3p .gt. 0d0) then
2327               t3 = ((dble(i3) - 5d0) * size3 - x3p) / v3p + tf
2328            else if (v3p .lt. 0d0) then
2329               t3 = ((dble(i3) - 6d0) * size3 - x3p) / v3p + tf
2330            else
2331               t3 = tlarge
2332            end if
2333
2334 c       if a particle is on the wall, we don't collide it on the same wall
2335
2336 c        if (t1 .eq. 0d0) t1 = tlarge
2337 c        if (t2 .eq. 0d0) t2 = tlarge
2338 c        if (t3 .eq. 0d0) t3 = tlarge
2339
2340            tmin = min(t1, t2, t3)
2341
2342 c       set icsta,
2343 c       after checking this is not an earlier collision comparing with
2344 c       a collision with another particle, we need to set icsta=0
2345 c       after checking whether there is also a particle collision
2346 c       at the same time, we need to reset the second bit of icsta
2347
2348            if (tmin .eq. t1) then
2349               if (v1p .gt. 0d0) then
2350                  icsta(i) = 101
2351               else
2352                  icsta(i) = 102
2353               end if
2354            end if
2355
2356            if (tmin .eq. t2) then
2357               if (v2p .gt. 0d0) then
2358                  icsta(i) = 103
2359               else
2360                  icsta(i) = 104
2361               end if
2362            end if
2363
2364            if (tmin .eq. t3) then
2365               if (v3p .gt. 0d0) then
2366                  icsta(i) = 105
2367               else
2368                  icsta(i) = 106
2369               end if
2370            end if
2371
2372         if (tmin .le. size) return
2373
2374         end if
2375
2376         if (v1p .gt. (i1 - 5) * v1) then
2377            t1 = ((i1 - 5) * (size1 - v1 * size) +
2378      &          v1p * tf - x1p) / (v1p - (i1 - 5) * v1)
2379         else if (v1p .lt. (i1 - 6) * v1) then
2380            t1 = ((i1 - 6) * (size1 - v1 * size) +
2381      &          v1p * tf - x1p) / (v1p - (i1 - 6) * v1)
2382         else
2383            t1 = tlarge
2384         end if
2385
2386         if (v2p .gt. (i2 - 5) * v2) then
2387            t2 = ((i2 - 5) * (size2 - v2 * size) +
2388      &          v2p * tf - x2p) / (v2p - (i2 - 5) * v2)
2389         else if (v2p .lt. (i2 - 6) * v2) then
2390            t2 = ((i2 - 6) * (size2 - v2 * size) +
2391      &          v2p * tf - x2p) / (v2p - (i2 - 6) * v2)
2392         else
2393            t2 = tlarge
2394         end if
2395
2396         if (v3p .gt. (i3 - 5) * v3) then
2397            t3 = ((i3 - 5) * (size3 - v3 * size) +
2398      &          v3p * tf - x3p) / (v3p - (i3 - 5) * v3)
2399         else if (v3p .lt. (i3 - 6) * v3) then
2400            t3 = ((i3 - 6) * (size3 - v3 * size) +
2401      &          v3p * tf - x3p) / (v3p - (i3 - 6) * v3)
2402         else
2403            t3 = tlarge
2404         end if
2405
2406 c       if a particle is on the wall, we don't collide it on the same wall
2407
2408 c        if (t1 .eq. 0d0) t1 = tlarge
2409 c        if (t2 .eq. 0d0) t2 = tlarge
2410 c        if (t3 .eq. 0d0) t3 = tlarge
2411
2412         tmin = min(t1, t2, t3)
2413
2414 c       set icsta,
2415 c       after checking this is not an earlier collision comparing with
2416 c       a collision with another particle, we need to set icsta=0
2417 c       after checking whether there is also a particle collision
2418 c       at the same time, we need to reset the second bit of icsta
2419
2420         if (tmin .eq. t1) then
2421            if (v1p .gt. (i1 - 5) * v1) then
2422               icsta(i) = 101
2423            else
2424               icsta(i) = 102
2425            end if
2426         end if
2427
2428         if (tmin .eq. t2) then
2429            if (v2p .gt. (i2 - 5) * v2) then
2430               icsta(i) = 103
2431            else
2432               icsta(i) = 104
2433            end if
2434         end if
2435
2436         if (tmin .eq. t3) then
2437            if (v3p .gt. (i3 - 5) * v3) then
2438               icsta(i) = 105
2439            else
2440               icsta(i) = 106
2441            end if
2442         end if
2443
2444         return
2445         end
2446
2447         subroutine wallc2(i, i1, i2, i3, t, tmin)
2448 c       this subroutine is used to get wall collision time
2449 c       when particle is inside the cube, it sets the icsta at the same time
2450 c       input i,i1,i2,i3,t
2451 c       output tmin, icsta(i)
2452 c       note the icsta is not finally set. we need further judgement in
2453 c       fixtim
2454
2455         implicit double precision (a-h, o-z)
2456         parameter (MAXPTN=400001)
2457
2458         common /para5/ iconfg, iordsc
2459 cc      SAVE /para5/
2460         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
2461      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
2462      &       xmass(MAXPTN), ityp(MAXPTN)
2463 cc      SAVE /prec2/
2464         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
2465 cc      SAVE /prec4/
2466         common /ilist1/
2467      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
2468      &     ictype, icsta(MAXPTN),
2469      &     nic(MAXPTN), icels(MAXPTN)
2470 cc      SAVE /ilist1/
2471         common /ilist3/ size1, size2, size3, v1, v2, v3, size
2472 cc      SAVE /ilist3/
2473         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2474 cc      SAVE /ilist5/
2475         SAVE
2476
2477         i1=i1
2478         i2=i2
2479         i3=i3
2480         t=t
2481         x1p = gx(i)
2482         x2p = gy(i)
2483         x3p = gz(i)
2484         tf = ft(i)
2485         v1p = vx(i)
2486         v2p = vy(i)
2487         v3p = vz(i)
2488
2489         if (v1p .gt. 0d0) then
2490            t1 = (5d0 * size1 - x1p) / v1p + tf
2491         else if (v1p .lt. 0d0) then
2492            t1 = (-5d0 * size1 - x1p) / v1p + tf
2493         else
2494            t1 = tlarge
2495         end if
2496
2497         if (v2p .gt. 0d0) then
2498            t2 = (5d0 * size2 - x2p) / v2p + tf
2499         else if (v2p .lt. 0d0) then
2500            t2 = (- 5d0 * size2 - x2p) / v2p +tf
2501         else
2502            t2 = tlarge
2503         end if
2504
2505         if (iconfg .eq. 5) then
2506            if (v3p .gt. 0d0) then
2507               t3 = (5d0 * size3 - x3p) / v3p + tf
2508            else if (v3p .lt. 0d0) then
2509               t3 = (- 5d0 * size3 - x3p) / v3p +tf
2510            else
2511               t3 = tlarge
2512            end if
2513         else
2514            t3 = tlarge
2515         end if
2516
2517         tmin = min(t1, t2, t3)
2518
2519 c       set icsta,
2520 c       after checking this is not an earlier collision comparing with
2521 c       a collision with another particle, we need to set icsta=0
2522 c       after checking whether there is also a particle collision
2523 c       at the same time, we need to reset the second bit of icsta
2524
2525         if (tmin .eq. t1) then
2526            if (v1p .gt. 0d0) then
2527               icsta(i) = 101
2528            else
2529               icsta(i) = 102
2530            end if
2531         end if
2532
2533         if (tmin .eq. t2) then
2534            if (v2p .gt. 0d0) then
2535               icsta(i) = 103
2536            else
2537               icsta(i) = 104
2538            end if
2539         end if
2540
2541         if (tmin .eq. t3) then
2542            if (v3p .gt. 0d0) then
2543               icsta(i) = 105
2544            else
2545               icsta(i) = 106
2546            end if
2547         end if
2548
2549         return
2550         end
2551
2552         subroutine wallcb(i, t, tmin)
2553 c       this subroutine is used to calculate the wall collision time
2554 c       when the particle is outside the cube
2555 c       input i,t
2556 c       output tmin,icsta(i)
2557 c       note the icsta is not finally set. we need further judgement in
2558 c       fixtim
2559
2560         implicit double precision (a-h, o-z)
2561         parameter (MAXPTN=400001)
2562
2563         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
2564      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
2565      &       xmass(MAXPTN), ityp(MAXPTN)
2566 cc      SAVE /prec2/
2567         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
2568 cc      SAVE /prec4/
2569         common /ilist1/
2570      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
2571      &     ictype, icsta(MAXPTN),
2572      &     nic(MAXPTN), icels(MAXPTN)
2573 cc      SAVE /ilist1/
2574         common /ilist3/ size1, size2, size3, v1, v2, v3, size
2575 cc      SAVE /ilist3/
2576         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2577 cc      SAVE /ilist5/
2578         SAVE
2579
2580 c       check if there is a collision by looking at the closest approach point
2581 c       and see if it's inside the cube
2582
2583         if (size1 .eq. 0d0 .or. size2 .eq. 0d0 .or.
2584      &     size3 .eq. 0d0) return
2585
2586         x1p = gx(i)
2587         x2p = gy(i)
2588         x3p = gz(i)
2589         v1p = vx(i)
2590         v2p = vy(i)
2591         v3p = vz(i)
2592         tf = ft(i)
2593
2594         if (t .lt. size .and. tf .lt. size) then
2595            if (x1p .lt. - 5d0 * size1 .and. v1p .gt. 0d0) then
2596               t1 = (- 5d0 * size1 - x1p) / v1p + tf
2597            else if(x1p .gt. 5d0 * size1 .and. v1p .lt. 0d0) then
2598               t1 = - (x1p - 5d0 * size1) / v1p + tf
2599            else
2600               t1 = tlarge
2601            end if
2602
2603            if (t1 .ne. tlarge) then
2604               x2pp = x2p + v2p * (t1 - tf)
2605               x3pp = x3p + v3p * (t1 - tf)
2606               if (x2pp .le. - 5d0 * size2 .or. x2pp .ge. 5d0 * size2
2607      &             .or. x3pp .le. - 5d0 * size3
2608      &             .or. x3pp .ge. 5d0 * size3)
2609      &             t1 = tlarge
2610            end if
2611
2612            if (x2p .lt. - 5d0 * size2 .and. v2p .gt. 0d0) then
2613               t2 = (- 5d0 * size2 - x2p) / v2p + tf
2614            else if(x2p .gt. 5d0 * size2 .and. v2p .lt. 0d0) then
2615               t2 = - (x2p - 5d0 * size2) / v2p + tf
2616            else
2617               t2 = tlarge
2618            end if
2619
2620            if (t2 .ne. tlarge) then
2621               x1pp = x1p + v1p * (t2 - tf)
2622               x3pp = x3p + v3p * (t2 - tf)
2623               if (x1pp .le. - 5d0 * size1 .or. x1pp .ge. 5d0 * size1
2624      &          .or. x3pp .le. - 5d0 * size3 .or. x3pp .ge. 5d0 * size3)
2625      &             t2 = tlarge
2626            end if
2627
2628            if (x3p .lt. - 5d0 * size3 .and. v3p .gt. 0d0) then
2629               t3 = (- 5d0 * size3 - x3p) / v3p + tf
2630            else if(x3p .gt. 5d0 * size3 .and. v3p .lt. 0d0) then
2631               t3 = - (x3p - 5d0 * size3) / v3p + tf
2632            else
2633               t3 = tlarge
2634            end if
2635
2636            if (t3 .ne. tlarge) then
2637               x1pp = x1p + v1p * (t3 - tf)
2638               x2pp = x2p + v2p * (t3 - tf)
2639               if (x1pp .le. - 5d0 * size1 .or. x1pp .ge. 5d0 * size1
2640      &          .or. x2pp .le. - 5d0 * size2 .or. x2pp .ge. 5d0 * size2)
2641      &             t3 = tlarge
2642            end if
2643
2644            tmin = min(t1, t2, t3)
2645
2646 c       set icsta,
2647 c       after checking this is not an earlier collision comparing with
2648 c       a collision with another particle, we need to set icsta=0
2649 c       after checking whether there is also a particle collision
2650 c       at the same time, we need to reset the second bit of icsta
2651
2652            if (tmin .eq. t1) then
2653               if (v1p .gt. 0d0) then
2654                  icsta(i) = 101
2655               else
2656                  icsta(i) = 102
2657               end if
2658            end if
2659
2660            if (tmin .eq. t2) then
2661               if (v2p .gt. 0d0) then
2662                  icsta(i) = 103
2663               else
2664                  icsta(i) = 104
2665               end if
2666            end if
2667
2668            if (tmin .eq. t3) then
2669               if (v3p .gt. 0d0) then
2670                  icsta(i) = 105
2671               else
2672                  icsta(i) = 106
2673               end if
2674            end if
2675
2676         if (tmin .le. size) return
2677
2678         end if
2679
2680 c       notice now x1q, x2q, x3q are coordinates at time t
2681         x1q = x1p + v1p * (t - tf)
2682         x2q = x2p + v2p * (t - tf)
2683         x3q = x3p + v3p * (t - tf)
2684
2685         if (x1q .lt. - 5d0 * (size1 + v1 * (t - size)) .and.
2686      &      v1p .gt. - 5d0 * v1) then
2687            t1 = (- 5d0 * (size1 - v1 * size) + v1p * tf - x1p) /
2688      &          (v1p - (- 5d0) * v1)
2689            icsta1 = 101
2690         else if (x1q .gt. 5d0 * (size1 + v1 * (t-size)) .and.
2691      &     v1p .lt. 5d0 * v1) then
2692            t1 = (5d0 * (size1 - v1 * size) + v1p * tf - x1p) /
2693      &          (v1p - 5d0 * v1)
2694            icsta1 = 102
2695         else
2696            t1 = tlarge
2697         end if
2698
2699         if (t1 .ne. tlarge) then
2700            x2pp = x2p + v2p * (t1 - tf)
2701            x3pp = x3p + v3p * (t1 - tf)
2702            if (x2pp .le. - 5d0 * (size2 + v2 * (t1 - size))
2703      &        .or. x2pp .ge. 5d0 * (size2 + v2 * (t1 - size))
2704      &        .or. x3pp .le. - 5d0 * (size3 + v3 * (t1 - size))
2705      &        .or. x3pp .ge. 5d0 * (size3 + v3 * (t1 - size)))
2706      &        t1 = tlarge
2707         end if
2708
2709         if (x2q .lt. - 5d0 * (size2 + v2 * (t - size)) .and.
2710      &     v2p .gt. - 5d0 * v2) then
2711            t2 = (- 5d0 * (size2 - v2 * size) + v2p * tf - x2p) /
2712      &          (v2p - (- 5d0) * v2)
2713            icsta2 = 103
2714         else if (x2q .gt. 5d0 * (size2 + v2 * (t - size)) .and.
2715      &     v2p .lt. 5d0 * v2) then
2716            t2 = (5d0 * (size2 - v2 * size) + v2p * tf - x2p) /
2717      &          (v2p - 5d0 * v2)
2718            icsta2 = 104
2719         else
2720            t2 = tlarge
2721         end if
2722
2723         if (t2 .ne. tlarge) then
2724            x1pp = x1p + v1p * (t2 - tf)
2725            x3pp = x3p + v3p * (t2 - tf)
2726            if (x1pp .le. - 5d0 * (size1 + v1 * (t2 - size))
2727      &        .or. x1pp .ge. 5d0 * (size1 + v1 * (t2 - size))
2728      &        .or. x3pp .le. - 5d0 * (size3 + v3 * (t2 - size))
2729      &        .or. x3pp .ge. 5d0 * (size3 + v3 * (t2 - size)))
2730      &        t2 = tlarge
2731         end if
2732
2733         if (x3q .lt. - 5d0 * (size3 + v3 * (t - size)) .and.
2734      &     v3p .gt. - 5d0 * v3) then
2735            t3 = (- 5d0 * (size3 - v3 * size) + v3p * tf - x3p) /
2736      &          (v3p - (- 5d0) * v3)
2737            icsta3 = 105
2738         else if (x3q .gt. 5d0 * (size3 + v3 * (t - size)) .and.
2739      &     v3p .lt. 5d0 * v3) then
2740            t3 = (5d0 * (size3 - v3 * size) + v3p * tf - x3p) /
2741      &          (v3p - 5d0 * v3)
2742            icsta3 = 106
2743         else
2744            t3 = tlarge
2745         end if
2746
2747         if (t3 .ne. tlarge) then
2748            x2pp = x2p + v2p * (t3 - tf)
2749            x1pp = x1p + v1p * (t3 - tf)
2750            if (x2pp .le. - 5d0 * (size2 + v2 * (t3 - size))
2751      &        .or. x2pp .ge. 5d0 * (size2 + v2 * (t3 - size))
2752      &        .or. x1pp .le. - 5d0 * (size1 + v1 * (t3 - size))
2753      &        .or. x1pp .ge. 5d0 * (size1 + v1 * (t3 - size)))
2754      &        t3 = tlarge
2755         end if
2756
2757         tmin = min(t1, t2, t3)
2758
2759 c       set icsta,
2760 c       after checking this is not an earlier collision comparing with
2761 c       a collision with another particle, we need to set icsta=0
2762 c       after checking whether there is also a particle collision
2763 c       at the same time, we need to reset the second bit of icsta
2764
2765         if (tmin .eq. t1) then
2766            icsta(i) = icsta1
2767         else if (tmin .eq. t2) then
2768            icsta(i) = icsta2
2769         else if (tmin .eq. t3) then
2770            icsta(i) = icsta3
2771         end if
2772
2773         return
2774         end
2775
2776         subroutine chkout(l, t, tmin, nc)
2777 c       this subroutine is used to check the collisions with particles in
2778 c       surface cells to see if we can get a smaller collision time than tmin
2779 c       with particle nc, when the colliding particle is outside the cube
2780 c       input l,t,tmin,nc
2781 c       output tmin, nc
2782
2783         implicit double precision (a-h, o-z)
2784         parameter (MAXPTN=400001)
2785
2786         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
2787      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
2788      &       xmass(MAXPTN), ityp(MAXPTN)
2789 cc      SAVE /prec2/
2790         SAVE
2791
2792         m1 = 11
2793         m2 = 11
2794         m3 = 11
2795         call chkcel(l, m1, m2, m3, t, tmin, nc)
2796
2797         do 1003 i = 1, 10
2798            do 1002 j = 1, 10
2799               do 1001 k = 1, 10
2800                  if (i .eq. 1 .or. i .eq. 10 .or. j .eq. 1
2801      &              .or. j .eq. 10 .or. k .eq. 1 .or. k .eq. 10)
2802      &                    call chkcel(l, i, j, k, t, tmin, nc)
2803  1001         continue
2804  1002      continue
2805  1003   continue
2806
2807         return
2808         end
2809
2810         subroutine chkin1(l, i1, i2, i3, t, tmin, nc)
2811 c       this subroutine is used to check collisions for particle inside
2812 c       the cube
2813 c       and update the afftected particles through chkcel
2814
2815         implicit double precision (a-h, o-z)
2816         SAVE
2817
2818 c       itest is a flag to make sure the 111111 cell is checked only once
2819         itest = 0
2820
2821         do 1003 i = i1 - 1, i1 + 1
2822            do 1002 j = i2 - 1, i2 + 1
2823               do 1001 k =  i3 - 1, i3 + 1
2824                  if (i .ge. 1 .and. i .le. 10 .and. j .ge. 1 .and.
2825      &               j .le. 10 .and. k .ge. 1 .and. k .le. 10) then
2826                     call chkcel(l, i, j, k, t, tmin, nc)
2827                  else if (itest .eq. 0) then
2828                     m1 = 11
2829                     m2 = 11
2830                     m3 = 11
2831                     call chkcel(l, m1, m2, m3, t, tmin, nc)
2832                     itest = 1
2833                  end if
2834  1001         continue
2835  1002      continue
2836  1003   continue
2837
2838         return
2839         end
2840
2841         subroutine chkin2(l, i1, i2, i3, t, tmin, nc)
2842 c       this subroutine is used to check collisions for particle inside
2843 c       the cube
2844 c       and update the afftected particles through chkcel
2845
2846         implicit double precision (a-h, o-z)
2847         SAVE
2848
2849 c       itest is a flag to make sure the 111111 cell is checked only once
2850         itest = 0
2851
2852         do 1003 i = i1 - 1, i1 + 1
2853            do 1002 j = i2 - 1, i2 + 1
2854               do 1001 k =  i3 - 1, i3 + 1
2855                  ia = i
2856                  ib = j
2857                  ic = k
2858                  if (k .ge. 1 .and. k .le. 10) then
2859                     if (i .eq. 0) ia = 10
2860                     if (i .eq. 11) ia = 1
2861                     if (j .eq. 0) ib = 10
2862                     if (j .eq. 11) ib = 1
2863                     call chkcel(l, ia, ib, ic, t, tmin, nc)
2864                  end if
2865  1001         continue
2866  1002      continue
2867  1003   continue
2868
2869         return
2870         end
2871
2872         subroutine chkin3(l, i1, i2, i3, t, tmin, nc)
2873 c       this subroutine is used to check collisions for particle inside
2874 c       the cube
2875 c       and update the afftected particles through chkcel
2876
2877         implicit double precision (a-h, o-z)
2878         SAVE
2879
2880 c       itest is a flag to make sure the 111111 cell is checked only once
2881         itest = 0
2882
2883         do 1003 i = i1 - 1, i1 + 1
2884            do 1002 j = i2 - 1, i2 + 1
2885               do 1001 k =  i3 - 1, i3 + 1
2886                  if (i .eq. 0) then
2887                     ia = 10
2888                  else if (i .eq. 11) then
2889                     ia = 1
2890                  else
2891                     ia = i
2892                  end if
2893                  if (j .eq. 0) then
2894                     ib = 10
2895                  else if (j .eq. 11) then
2896                     ib = 1
2897                  else
2898                     ib = j
2899                  end if
2900                  if (k .eq. 0) then
2901                     ic = 10
2902                  else if (k .eq. 11) then
2903                     ic = 1
2904                  else
2905                     ic = k
2906                  end if
2907                  call chkcel(l, ia, ib, ic, t, tmin, nc)
2908  1001         continue
2909  1002      continue
2910  1003   continue
2911
2912         return
2913         end
2914
2915         subroutine chkcel(il, i1, i2, i3, t, tmin, nc)
2916 c       this program is used to check through all the particles
2917 c       in the cell (i1,i2,i3) and see if we can get a particle collision
2918 c       with time less than the original input tmin ( the collision time of
2919 c       il with the wall
2920 c       and update the affected particles
2921
2922         implicit double precision (a-h, o-z)
2923         parameter (MAXPTN=400001)
2924
2925         common /para5/ iconfg, iordsc
2926 cc      SAVE /para5/
2927         common /ilist1/
2928      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
2929      &     ictype, icsta(MAXPTN),
2930      &     nic(MAXPTN), icels(MAXPTN)
2931 cc      SAVE /ilist1/
2932         common /ilist2/ icell, icel(10, 10, 10)
2933 cc      SAVE /ilist2/
2934         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
2935 cc      SAVE /ilist4/
2936         SAVE
2937
2938         if (iconfg .eq. 3 .or. iconfg .eq. 5) then
2939            jj = ichkpt
2940            do 1001 j = 1, jj
2941               call ck(j, ick)
2942 c     10/24/02 get rid of argument usage mismatch in ud2():
2943                             jud2=j
2944 c              if (ick .eq. 1) call ud2(j, il, t, tmin, nc)
2945               if (ick .eq. 1) call ud2(jud2, il, t, tmin, nc)
2946  1001      continue
2947            return
2948         end if
2949
2950         if (i1 .eq. 11 .and. i2 .eq. 11 .and. i3 .eq. 11) then
2951            l = icell
2952         else
2953            l = icel(i1, i2, i3)
2954         end if
2955
2956 c       if there is no particle
2957         if (l .eq. 0) then
2958            return
2959         end if
2960         j = nic(l)
2961 c       if there is only one particle
2962         if (j .eq. 0) then
2963            call ck(l, ick)
2964            if (ick .eq. 1) call ud2(l, il, t, tmin, nc)
2965
2966 c       if there are many particles
2967         else
2968
2969 c       we don't worry about the other colliding particle because it's
2970 c       set in last(), and will be checked in ud2
2971
2972            call ck(l, ick)
2973            if (ick .eq. 1) call ud2(l, il, t, tmin, nc)
2974
2975            do 10 while(j .ne. l)
2976               call ck(j, ick)
2977               if (ick .eq. 1) call ud2(j, il, t, tmin, nc)
2978               j = nic(j)
2979  10           continue
2980         end if
2981
2982         return
2983         end
2984
2985         subroutine ck(l, ick)
2986 c       this subroutine is used for chcell to check whether l should be
2987 c       checked or not for updating tmin, nc
2988 c       input l
2989 c       output ick
2990 c       if ick=1, l should be checked, otherwise it should not be.
2991
2992         implicit double precision (a-h, o-z)
2993         parameter (MAXPTN=400001)
2994
2995         common /ilist1/
2996      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
2997      &     ictype, icsta(MAXPTN),
2998      &     nic(MAXPTN), icels(MAXPTN)
2999 cc      SAVE /ilist1/
3000         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
3001 cc      SAVE /ilist4/
3002         SAVE
3003
3004         ick = 1
3005         if (ictype .eq. 1) then
3006            if (l .eq. ifmpt) ick = 0
3007         else if (ictype .eq. 0 .or. ictype .eq. 3 .or.
3008      &     ictype .eq. 4) then
3009            if (l .eq. iscat .or. l .eq. jscat) ick = 0
3010         else
3011            if (l .eq. iscat .or. l .eq. jscat .or.
3012      &         l .eq. ifmpt) ick = 0
3013         end if
3014 c       notice il is either iscat or jscat, or ifmpt, we deal with them
3015 c       seperately according to ictype
3016
3017         return
3018         end
3019
3020         subroutine dchout(l, ii, t)
3021 c       this subroutine is used to check collisions of l with particles when
3022 c       l is outside the cube and the collision just happened is a collision
3023 c       including a collision with wall (hence we need to use dcheck to throw
3024 c       away old collisions that are not in the new neighboring cells.
3025
3026 c       input l,t
3027
3028         implicit double precision (a-h, o-z)
3029         parameter (MAXPTN=400001)
3030
3031         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
3032      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
3033      &       xmass(MAXPTN), ityp(MAXPTN)
3034 cc      SAVE /prec2/
3035         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
3036 cc      SAVE /prec4/
3037         common /ilist3/ size1, size2, size3, v1, v2, v3, size
3038 cc      SAVE /ilist3/
3039         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3040 cc      SAVE /ilist5/
3041         SAVE
3042
3043         external integ
3044
3045         tt = ft(l)
3046         td = t - size
3047         x1 = gx(l) + vx(l) * (t - tt)
3048         x2 = gy(l) + vy(l) * (t - tt)
3049         x3 = gz(l) + vz(l) * (t - tt)
3050         if (td .le. 0d0) then
3051            i1 = integ(x1 / size1) + 6
3052            i2 = integ(x2 / size2 ) + 6
3053            i3 = integ(x3 / size3 ) + 6
3054            if (integ(x1 / size1) .eq. x1 / size1 .and. vx(l) .lt. 0d0)
3055      &        i1 = i1 - 1
3056            if (integ(x2 / size2) .eq. x2 / size2 .and. vy(l) .lt. 0d0)
3057      &        i2 = i2 - 1
3058            if (integ(x3 / size3) .eq. x3 / size3 .and. vz(l) .lt. 0d0)
3059      &        i3 = i3 - 1
3060         else
3061            i1 = integ(x1 / (size1 + v1 * td)) + 6
3062            i2 = integ(x2 / (size2 + v2 * td)) + 6
3063            i3 = integ(x3 / (size3 + v3 * td)) + 6
3064 c     10/24/02 (i) below should be (l):
3065            if (integ(x1 / (size1 + v1 * td)) .eq.
3066      &        x1 / (size1 +v1 * td) .and.
3067      &        vx(l) .lt. (i1 - 6) * v1) i1 = i1 - 1
3068 c     &        vx(i) .lt. (i1 - 6) * v1) i1 = i1 - 1
3069            if (integ(x2 / (size2 + v2 * td)) .eq.
3070      &        x2 / (size2 + v2 * td) .and.
3071      &        vy(l) .lt. (i2 - 6) * v2) i2 = i2 - 1
3072 c     &        vy(i) .lt. (i2 - 6) * v2) i2 = i2 - 1
3073            if (integ(x3 / (size3 + v3 * td)) .eq.
3074      &        x3 / (size3 + v3 * td) .and.
3075      &        vz(l) .lt. (i3 - 6) * v3) i3 = i3 - 1
3076 c     &        vz(i) .lt. (i3 - 6) * v3) i3 = i3 - 1
3077         end if
3078
3079         if (ii .eq. 1) then
3080            i = 9
3081            do 1002 j = i2 - 1, i2 + 1
3082               do 1001 k = i3 - 1, i3 + 1
3083                  if (j .ge. 1 .and. j .le. 10 .and. k .ge. 1 .and.
3084      &              k .le. 10) then
3085                     call dchcel(l, i, j, k, t)
3086                  end if
3087  1001         continue
3088  1002      continue
3089         end if
3090
3091         if (ii .eq. 2) then
3092            i = 2
3093            do 1004 j = i2 - 1, i2 + 1
3094               do 1003 k = i3 - 1, i3 + 1
3095                  if (j .ge. 1 .and. j .le. 10 .and. k .ge. 1 .and.
3096      &              k .le. 10) then
3097                     call dchcel(l, i, j, k, t)
3098                  end if
3099  1003         continue
3100  1004      continue
3101         end if
3102
3103         if (ii .eq. 3) then
3104            j = 9
3105            do 1006 i = i1 - 1, i1 + 1
3106               do 1005 k = i3 - 1, i3 + 1
3107                  if (i .ge. 1 .and. i .le. 10 .and. k .ge. 1 .and.
3108      &              k .le. 10) then
3109                     call dchcel(l, i, j, k, t)
3110                  end if
3111  1005         continue
3112  1006      continue
3113         end if
3114
3115         if (ii .eq. 4) then
3116            j = 2
3117            do 1008 i = i1 - 1, i1 + 1
3118               do 1007 k = i3 - 1, i3 + 1
3119                  if (i .ge. 1 .and. i .le. 10 .and. k .ge. 1 .and.
3120      &              k .le. 10) then
3121                     call dchcel(l, i, j, k, t)
3122                  end if
3123  1007         continue
3124  1008      continue
3125         end if
3126
3127         if (ii .eq. 5) then
3128            k = 9
3129            do 1010 i = i1 - 1, i1 + 1
3130               do 1009 j = i2 - 1, i2 + 1
3131                  if (i .ge. 1 .and. i .le. 10 .and. j .ge. 1 .and.
3132      &              j .le. 10) then
3133                     call dchcel(l, i, j, k, t)
3134                  end if
3135  1009         continue
3136  1010      continue
3137         end if
3138
3139         if (ii .eq. 6) then
3140            k = 2
3141            do 1012 i = i1 - 1, i1 + 1
3142               do 1011 j = i2 - 1, i2 + 1
3143                  if (i .ge. 1 .and. i .le. 10 .and. j .ge. 1 .and.
3144      &              j .le. 10) then
3145                     call dchcel(l, i, j, k, t)
3146                  end if
3147  1011         continue
3148  1012      continue
3149         end if
3150
3151         return
3152         end
3153
3154         subroutine dchin1(l, ii, i1, i2, i3, t)
3155 c       this subroutine is used to check collisions for particle inside
3156 c       the cube when the collision just happened is a collision including
3157 c       collision with wall
3158 c       and update the afftected particles through chkcel
3159
3160 c       input l,ii(specifying the direction of the wall collision),
3161 c          i1,i2,i3, (specifying the position of the cell
3162 c                    we are going to check)
3163 c          t
3164
3165         implicit double precision (a-h, o-z)
3166         SAVE
3167
3168 c       itest is a flag to make sure the 111111 cell is checked only once
3169         itest = 0
3170
3171         if (ii .eq. 1) then
3172            if (i1 .eq. 1) goto 100
3173            if (i1 .eq. 2) then
3174               if (i2 .ge. 2 .and. i2 .le. 9 .and. i3 .ge. 2 .and.
3175      &           i3 .le. 9) then
3176                  i = 11
3177                  j = 11
3178                  k = 11
3179                  call dchcel(l, i, j, k, t)
3180               end if
3181               goto 100
3182            end if
3183            i = i1 - 2
3184            do 1002 j = i2 - 1, i2 + 1
3185               do 1001 k = i3 - 1, i3 + 1
3186                  if (j .ge. 1 .and. j .le. 10 .and. k .ge. 1 .and.
3187      &              k .le. 10)
3188      &                    call dchcel(l, i, j, k, t)
3189  1001         continue
3190  1002      continue
3191         end if
3192
3193         if (ii .eq. 2) then
3194            if (i1 .eq. 10) goto 100
3195            if (i1 .eq. 9) then
3196               if (i2 .ge. 2 .and. i2 .le. 9 .and. i3 .ge. 2 .and.
3197      &           i3 .le. 9) then
3198                  i = 11
3199                  j = 11
3200                  k = 11
3201                  call dchcel(l, i, j, k, t)
3202               end if
3203               goto 100
3204            end if
3205            i = i1 + 2
3206            do 1004 j = i2 - 1, i2 + 1
3207               do 1003 k = i3 - 1, i3 + 1
3208                  if (j .ge. 1 .and. j .le. 10 .and. k .ge. 1 .and.
3209      &              k .le. 10)
3210      &                    call dchcel(l, i, j, k, t)
3211  1003         continue
3212  1004      continue
3213         end if
3214
3215         if (ii .eq. 3) then
3216            if (i2 .eq. 1) goto 100
3217            if (i2 .eq. 2) then
3218               if (i1 .ge. 2 .and. i1 .le. 9 .and. i3 .ge. 2 .and.
3219      &           i3 .le. 9) then
3220                  i = 11
3221                  j = 11
3222                  k = 11
3223                  call dchcel(l, i, j, k, t)
3224               end if
3225               goto 100
3226            end if
3227            j = i2 - 2
3228            do 1006 i = i1 - 1, i1 + 1
3229               do 1005 k = i3 - 1, i3 + 1
3230                  if (i .ge. 1 .and. i .le. 10 .and. k .ge. 1 .and.
3231      &              k .le. 10)
3232      &              call dchcel(l, i, j, k, t)
3233  1005         continue
3234  1006      continue
3235         end if
3236
3237         if (ii .eq. 4) then
3238            if (i2 .eq. 10) goto 100
3239            if (i2 .eq. 9) then
3240               if (i1 .ge. 2 .and. i1 .le. 9 .and. i3 .ge. 2 .and.
3241      &           i3 .le. 9) then
3242                  i = 11
3243                  j = 11
3244                  k = 11
3245                  call dchcel(l, i, j, k, t)
3246               end if
3247               goto 100
3248            end if
3249            j = i2 + 2
3250            do 1008 i = i1 - 1, i1 + 1
3251               do 1007 k = i3 - 1, i3 + 1
3252                  if (i .ge. 1 .and. i .le. 10 .and. k .ge. 1 .and.
3253      &           k .le. 10)
3254      &                 call dchcel(l, i, j, k, t)
3255  1007         continue
3256  1008      continue
3257         end if
3258
3259         if (ii .eq. 5) then
3260            if (i3 .eq. 1) goto 100
3261            if (i3 .eq. 2) then
3262               if (i1 .ge. 2 .and. i1 .le. 9 .and. i2 .ge. 2 .and.
3263      &           i2 .le. 9) then
3264                  i = 11
3265                  j = 11
3266                  k = 11
3267                  call dchcel(l, i, j, k, t)
3268               end if
3269               goto 100
3270            end if
3271            k = i3 - 2
3272            do 1010 i = i1 - 1, i1 + 1
3273               do 1009 j = i2 - 1, i2 + 1
3274                  if (i .ge. 1 .and. i .le. 10 .and. j .ge. 1 .and.
3275      &           j .le. 10)
3276      &                 call dchcel(l, i, j, k, t)
3277  1009         continue
3278  1010      continue
3279         end if
3280
3281         if (ii .eq. 6) then
3282            if (i3 .eq. 10) goto 100
3283            if (i3 .eq. 9) then
3284               if (i1 .ge. 2 .and. i1 .le. 9 .and. i2 .ge. 2 .and.
3285      &           i2 .le. 9) then
3286                  i = 11
3287                  j = 11
3288                  k = 11
3289                  call dchcel(l, i, j, k, t)
3290               end if
3291               goto 100
3292            end if
3293            k = i3 + 2
3294            do 1012 i = i1 - 1, i1 + 1
3295               do 1011 j = i2 - 1, i2 + 1
3296                  if (i .ge. 1 .and. i .le. 10 .and. j .ge. 1 .and.
3297      &           j .le. 10)
3298      &                 call dchcel(l, i, j, k, t)
3299  1011         continue
3300  1012      continue
3301         end if
3302
3303  100        continue
3304
3305         return
3306         end
3307
3308         subroutine dchin2(l, ii, i1, i2, i3, t)
3309 c       this subroutine is used to check collisions for particle inside
3310 c       the cube when the collision just happened is a collision including
3311 c       collision with wall
3312 c       and update the afftected particles through chkcel
3313
3314 c       input l,ii(specifying the direction of the wall collision),
3315 c          i1,i2,i3, (specifying the position of the cell
3316 c                    we are going to check)
3317 c          t
3318
3319         implicit double precision (a-h, o-z)
3320         SAVE
3321
3322         if (ii .eq. 1) then
3323            i = i1 - 2
3324            if (i .le. 0) i = i + 10
3325            ia = i
3326            do 1002 j = i2 - 1, i2 + 1
3327               do 1001 k = i3 - 1, i3 + 1
3328                  ib = j
3329                  ic = k
3330                  if (j .eq. 0) ib = 10
3331                  if (j .eq. 11) ib = 1
3332                  if (k .ge. 1 .and. k .le. 10) then
3333                     call dchcel(l, ia, ib, ic, t)
3334                  end if
3335  1001         continue
3336  1002      continue
3337         end if
3338
3339         if (ii .eq. 2) then
3340            i = i1 + 2
3341            if (i .ge. 11) i = i - 10
3342            ia = i
3343            do 1004 j = i2 - 1, i2 + 1
3344               do 1003 k = i3 - 1, i3 + 1
3345                  ib = j
3346                  ic = k
3347                  if (j .eq. 0) ib = 10
3348                  if (j .eq. 11) ib = 1
3349                  if (k .ge. 1 .and. k .le. 10) then
3350                     call dchcel(l, ia, ib, ic, t)
3351                  end if
3352  1003         continue
3353  1004      continue
3354         end if
3355
3356         if (ii .eq. 3) then
3357            j = i2 - 2
3358            if (j .le. 0) j = j + 10
3359            ib = j
3360            do 1006 i = i1 - 1, i1 + 1
3361               do 1005 k = i3 - 1, i3 + 1
3362                  ia = i
3363                  ic = k
3364                  if (i .eq. 0) ia = 10
3365                  if (i .eq. 11) ia = 1
3366                  if (k .ge. 1 .and. k .le. 10) then
3367                     call dchcel(l, ia, ib, ic, t)
3368                  end if
3369  1005         continue
3370  1006      continue
3371         end if
3372
3373         if (ii .eq. 4) then
3374            j = i2 + 2
3375            if (j .ge. 11) j = j - 10
3376            ib = j
3377            do 1008 i = i1 - 1, i1 + 1
3378               do 1007 k = i3 - 1, i3 + 1
3379                  ia = i
3380                  ic = k
3381                  if (i .eq. 0) ia = 10
3382                  if (i .eq. 11) ia = 1
3383                  if (k .ge. 1 .and. k .le. 10) then
3384                     call dchcel(l, ia, ib, ic, t)
3385                  end if
3386  1007         continue
3387  1008      continue
3388         end if
3389
3390         if (ii .eq. 5) then
3391            if (i3 .eq. 2) goto 100
3392            k = i3 - 2
3393            ic = k
3394            do 1010 i = i1 - 1, i1 + 1
3395               do 1009 j = i2 - 1, i2 + 1
3396                  ia = i
3397                  ib = j
3398                  if (i .eq. 0) ia = 10
3399                  if (i .eq. 11) ia = 1
3400                  if (j .eq. 0) ib = 10
3401                  if (j .eq. 11) ib = 1
3402                      call dchcel(l, ia, ib, ic, t)
3403  1009         continue
3404  1010      continue
3405         end if
3406
3407         if (ii .eq. 6) then
3408            if (i3 .eq. 9) goto 100
3409            k = i3 + 2
3410            ic = k
3411            do 1012 i = i1 - 1, i1 + 1
3412               do 1011 j = i2 - 1, i2 + 1
3413                  ia = i
3414                  ib = j
3415                  if (i .eq. 0) ia = 10
3416                  if (i .eq. 11) ia = 1
3417                  if (j .eq. 0) ib = 10
3418                  if (j .eq. 11) ib = 1
3419                      call dchcel(l, ia, ib, ic, t)
3420  1011         continue
3421  1012      continue
3422         end if
3423
3424  100        continue
3425
3426         return
3427         end
3428
3429         subroutine dchin3(l, ii, i1, i2, i3, t)
3430 c       this subroutine is used to check collisions for particle inside
3431 c       the cube when the collision just happened is a collision including
3432 c       collision with wall
3433 c       and update the afftected particles through chkcel
3434
3435 c       input l,ii(specifying the direction of the wall collision),
3436 c          i1,i2,i3, (specifying the position of the cell
3437 c                    we are going to check)
3438 c          t
3439
3440         implicit double precision (a-h, o-z)
3441         SAVE
3442
3443         if (ii .eq. 1) then
3444            i = i1 - 2
3445            if (i .le. 0) i = i + 10
3446            ia = i
3447            do 1002 j = i2 - 1, i2 + 1
3448               do 1001 k = i3 - 1, i3 + 1
3449                  ib = j
3450                  ic = k
3451                  if (j .eq. 0) ib = 10
3452                  if (j .eq. 11) ib = 1
3453                  if (k .eq. 0) ic = 10
3454                  if (k .eq. 11) ic = 1
3455                  call dchcel(l, ia, ib, ic, t)
3456  1001         continue
3457  1002      continue
3458         end if
3459
3460         if (ii .eq. 2) then
3461            i = i1 + 2
3462            if (i .ge. 11) i = i - 10
3463            ia = i
3464            do 1004 j = i2 - 1, i2 + 1
3465               do 1003 k = i3 - 1, i3 + 1
3466                  ib = j
3467                  ic = k
3468                  if (j .eq. 0) ib = 10
3469                  if (j .eq. 11) ib = 1
3470                  if (k .eq. 0) ic = 10
3471                  if (k .eq. 11) ic = 1
3472                  call dchcel(l, ia, ib, ic, t)
3473  1003         continue
3474  1004      continue
3475         end if
3476
3477         if (ii .eq. 3) then
3478            j = i2 - 2
3479            if (j .le. 0) j = j + 10
3480            ib = j
3481            do 1006 i = i1 - 1, i1 + 1
3482               do 1005 k = i3 - 1, i3 + 1
3483                  ia = i
3484                  ic = k
3485                  if (i .eq. 0) ia = 10
3486                  if (i .eq. 11) ia = 1
3487                  if (k .eq. 0) ic = 10
3488                  if (k .eq. 11) ic = 1
3489                  call dchcel(l, ia, ib, ic, t)
3490  1005         continue
3491  1006      continue
3492         end if
3493
3494         if (ii .eq. 4) then
3495            j = i2 + 2
3496            if (j .ge. 11) j = j - 10
3497            ib = j
3498            do 1008 i = i1 - 1, i1 + 1
3499               do 1007 k = i3 - 1, i3 + 1
3500                  ia = i
3501                  ic = k
3502                  if (i .eq. 0) ia = 10
3503                  if (i .eq. 11) ia = 1
3504                  if (k .eq. 0) ic = 10
3505                  if (k .eq. 11) ic = 1
3506                  call dchcel(l, ia, ib, ic, t)
3507  1007         continue
3508  1008      continue
3509         end if
3510
3511         if (ii .eq. 5) then
3512            k = i3 - 2
3513            if (k .le. 0) k = k + 10
3514            ic = k
3515            do 1010 i = i1 - 1, i1 + 1
3516               do 1009 j = i2 - 1, i2 + 1
3517                  ia = i
3518                  ib = j
3519                  if (i .eq. 0) ia = 10
3520                  if (i .eq. 11) ia = 1
3521                  if (j .eq. 0) ib = 10
3522                  if (j .eq. 11) ib = 1
3523                      call dchcel(l, ia, ib, ic, t)
3524  1009         continue
3525  1010      continue
3526         end if
3527
3528         if (ii .eq. 6) then
3529            k = i3 + 2
3530            if (k .ge. 11) k = k - 10
3531            ic = k
3532            do 1012 i = i1 - 1, i1 + 1
3533               do 1011 j = i2 - 1, i2 + 1
3534                  ia = i
3535                  ib = j
3536                  if (i .eq. 0) ia = 10
3537                  if (i .eq. 11) ia = 1
3538                  if (j .eq. 0) ib = 10
3539                  if (j .eq. 11) ib = 1
3540                      call dchcel(l, ia, ib, ic, t)
3541  1011         continue
3542  1012      continue
3543         end if
3544 c
3545         return
3546         end
3547
3548         subroutine dchcel(l, i, j, k, t)
3549 c       this subroutine is used to recalculate next collision time for
3550 c       particles in the cell i,j,k if the next collision partener is l
3551
3552         implicit double precision (a-h, o-z)
3553         parameter (MAXPTN=400001)
3554
3555         common /ilist1/
3556      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
3557      &     ictype, icsta(MAXPTN),
3558      &     nic(MAXPTN), icels(MAXPTN)
3559 cc      SAVE /ilist1/
3560         common /ilist2/ icell, icel(10, 10, 10)
3561 cc      SAVE /ilist2/
3562         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3563 cc      SAVE /ilist5/
3564         SAVE
3565
3566         if (i .eq. 11 .or. j .eq. 11 .or. k .eq. 11) then
3567            if ( .not. (i .eq. 11 .and. j .eq. 11 .and.
3568      &     k .eq. 11)) stop 'cerr'
3569            m = icell
3570         else
3571            m = icel(i, j, k)
3572         end if
3573
3574         if (m .eq. 0) return
3575         if (next(m) .eq. l) then
3576            tm = tlarge
3577            last0 = 0
3578            call reor(t, tm, m, last0)
3579         end if
3580         n = nic(m)
3581         if (n .eq. 0) return
3582         do 10 while(n .ne. m)
3583            if (next(n) .eq. l) then
3584               tm = tlarge
3585               last0 = 0
3586               call reor(t, tm, n, last0)
3587            end if
3588            n = nic(n)
3589  10        continue
3590
3591         return
3592         end
3593
3594         subroutine fixtim(l, t, tmin1, tmin, nc)
3595 c       this subroutine is used to compare the collision time with wall tmin1
3596 c       and new collision time with particles for particle l
3597 c       when used in ulist, input nc may be 0, which indicates no particle
3598 c       collisions happen before wall collision, of course, then tmin=tmin1
3599
3600         implicit double precision (a-h, o-z)
3601         parameter (MAXPTN=400001)
3602
3603         common /ilist1/
3604      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
3605      &     ictype, icsta(MAXPTN),
3606      &     nic(MAXPTN), icels(MAXPTN)
3607 cc      SAVE /ilist1/
3608         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3609 cc      SAVE /ilist5/
3610         SAVE
3611
3612         t=t
3613         k = nc
3614         if (tmin .lt. tmin1) then
3615            ot(l) = tmin
3616            if (ct(l) .lt. tmin1) then
3617               icsta(l) = 0
3618            else
3619               icsta(l) = icsta(l) + 10
3620            end if
3621            next(l) = k
3622         else if (tmin .eq. tmin1) then
3623            ot(l) = tmin
3624            if (nc .eq. 0) then
3625               next(l) = 0
3626            else
3627               icsta(l) = icsta(l) + 10
3628               next(l) = k
3629            end if
3630         else
3631            ot(l) = tmin1
3632            next(l) = 0
3633         end if
3634
3635         return
3636         end
3637
3638         subroutine ud2(i, j, t, tmin, nc)
3639 c       this subroutine is used to update next(i), ct(i), ot(i),
3640 c        and get tmin, nc for j
3641
3642         implicit double precision (a-h, o-z)
3643         parameter (MAXPTN=400001)
3644
3645         common /para5/ iconfg, iordsc
3646 cc      SAVE /para5/
3647         common /aurec1/ jxa, jya, jza
3648 cc      SAVE /aurec1/
3649         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
3650 cc      SAVE /aurec2/
3651         common /ilist1/
3652      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
3653      &     ictype, icsta(MAXPTN),
3654      &     nic(MAXPTN), icels(MAXPTN)
3655 cc      SAVE /ilist1/
3656         common /ilist3/ size1, size2, size3, v1, v2, v3, size
3657 cc      SAVE /ilist3/
3658         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3659 cc      SAVE /ilist5/
3660         SAVE
3661
3662         logical allok
3663
3664         call isco(i, j, allok, tm, t1, t2)
3665
3666         if (allok) then
3667 c       tm eq tmin, change nc to make sure fixtime get the collision with both
3668 c       wall and particle
3669
3670              if (tm .lt. tmin) then
3671                 tmin = tm
3672                 ct(j) = t2
3673                 nc = i
3674                 if (iconfg .eq. 3 .or. iconfg .eq. 5) then
3675                    dgxa(j) = jxa * 10d0 * size1
3676                    dgya(j) = jya * 10d0 * size2
3677                    if (iconfg .eq. 5) then
3678                       dgza(j) = jza * 10d0 * size3
3679                    end if
3680                 end if
3681              end if
3682
3683              if (tm .le. ot(i)) then
3684                 ct(i) = t1
3685                 icels0 = icels(i)
3686                 i1 = icels0 / 10000
3687                 i2 = (icels0 - i1 * 10000) / 100
3688                 i3 = icels0 - i1 * 10000 - i2 * 100
3689                 call wallc(i, i1, i2, i3, t, tmin1)
3690                 call fixtim(i, t, tmin1, tm, j)
3691                 if (iconfg .eq. 3 .or. iconfg .eq. 5) then
3692                    dgxa(i) = - jxa * 10d0 * size1
3693                    dgya(i) = - jya * 10d0 * size2
3694                    if (iconfg .eq. 5) then
3695                       dgza(i) = - jza * 10d0 * size3
3696                    end if
3697                 end if
3698              end if
3699
3700              if (tm .gt. ot(i) .and. next(i) .eq. j) then
3701                 ct(i) = t1
3702                 call reor(t, tm, i, j)
3703              end if
3704
3705            else if (next(i) .eq. j) then
3706
3707              tm = tlarge
3708
3709              call reor(t, tm, i, j)
3710
3711           end if
3712
3713         return
3714         end
3715
3716         subroutine isco(i, j, allok, tm, t1, t2)
3717
3718         implicit double precision (a-h, o-z)
3719
3720         common /para5/ iconfg, iordsc
3721 cc      SAVE /para5/
3722         SAVE
3723
3724         logical allok
3725
3726         iorder = iordsc / 10
3727         if (iconfg .eq. 1) then
3728            if (iorder .eq. 1) then
3729               call isco1(i, j, allok, tm, t1, t2)
3730            else if (iorder .eq. 2) then
3731               call isco2(i, j, allok, tm, t1, t2)
3732            else if (iorder .eq. 3) then
3733               call isco3(i, j, allok, tm, t1, t2)
3734            end if
3735         else if (iconfg .eq. 2 .or. iconfg .eq. 4) then
3736            if (iorder .eq. 1) then
3737               call isco4(i, j, allok, tm, t1, t2)
3738            else if (iorder .eq. 2) then
3739               call isco5(i, j, allok, tm, t1, t2)
3740            else if (iorder .eq. 3) then
3741               call isco6(i, j, allok, tm, t1, t2)
3742            end if
3743         else if (iconfg .eq. 3) then
3744            if (iorder .eq. 1) then
3745               call isco7(i, j, allok, tm, t1, t2)
3746            else if (iorder .eq. 2) then
3747               call isco8(i, j, allok, tm, t1, t2)
3748            else if (iorder .eq. 3) then
3749               call isco9(i, j, allok, tm, t1, t2)
3750            end if
3751         else if (iconfg .eq. 5) then
3752            if (iorder .eq. 1) then
3753               call isco10(i, j, allok, tm, t1, t2)
3754            else if (iorder .eq. 2) then
3755               call isco11(i, j, allok, tm, t1, t2)
3756            else if (iorder .eq. 3) then
3757               call isco12(i, j, allok, tm, t1, t2)
3758            end if
3759         end if
3760
3761         return
3762         end
3763
3764         subroutine isco1(i, j, allok, tm, t1, t2)
3765 c       this subroutine is used to decide whether there is a collision between
3766 c       particle i and j, if there is one allok=1, and tm gives the
3767 c       collision time, t1 the collision time for i,
3768 c       t2 the collision time for j
3769
3770         implicit double precision (a-h, o-z)
3771         parameter (MAXPTN=400001)
3772
3773         common /para2/ xmp, xmu, alpha, rscut2, cutof2
3774 cc      SAVE /para2/
3775         common /para5/ iconfg, iordsc
3776 cc      SAVE /para5/
3777         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
3778      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
3779      &       xmass(MAXPTN), ityp(MAXPTN)
3780 cc      SAVE /prec2/
3781         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
3782 cc      SAVE /prec4/
3783         common /ilist1/
3784      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
3785      &     ictype, icsta(MAXPTN),
3786      &     nic(MAXPTN), icels(MAXPTN)
3787 cc      SAVE /ilist1/
3788         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3789 cc      SAVE /ilist5/
3790         SAVE
3791
3792         logical allok
3793
3794 c       preventing consecutive collisions
3795         allok = last(i) .ne. j .or. last(j) .ne. i
3796
3797 c       set up numbers for later calculations
3798         i1 = i
3799         i2 = j
3800
3801         p4 = ft(i2) - ft(i1)
3802         p1 = gx(i2) - gx(i1)
3803         p2 = gy(i2) - gy(i1)
3804         p3 = gz(i2) - gz(i1)
3805
3806         q4 = e(i1)
3807         q1 = px(i1)
3808         q2 = py(i1)
3809         q3 = pz(i1)
3810
3811         r4 = e(i2)
3812         r1 = px(i2)
3813         r2 = py(i2)
3814         r3 = pz(i2)
3815
3816         a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
3817         b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
3818         c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
3819         d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
3820         ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
3821         f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
3822
3823 c       make sure particle 2 formed early
3824         h = a + b
3825         if (h .gt. 0d0) then
3826            g = a
3827            a = -b
3828            b = -g
3829
3830            g = c
3831            c = d
3832            d = g
3833
3834            i1 = j
3835            i2 = i
3836         end if
3837
3838 c       check the approaching criteria
3839         if (allok) then
3840
3841            vp = a * d - b * ee
3842
3843            allok = allok .and. vp .lt. 0d0
3844
3845         end if
3846
3847 c       check the closest approach distance criteria
3848          if (allok) then
3849
3850            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
3851      &           (ee ** 2 - c * d)
3852
3853            allok = allok .and. dm2 .lt. cutof2
3854
3855         end if
3856
3857 c       check the time criteria
3858         if (allok) then
3859
3860            tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
3861            tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
3862            tm = 0.5d0 * (tc1 + tc2)
3863
3864            allok = allok .and. tm .gt. ft(i) .and. tm .gt. ft(j)
3865
3866         end if
3867
3868 c        check rts cut
3869         if (allok) then
3870
3871            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
3872      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
3873
3874            allok = allok .and. rts2 .gt. rscut2
3875         end if
3876
3877         if (.not. allok) then
3878            tm = tlarge
3879            t1 = tlarge
3880            t2 = tlarge
3881         else if (h .gt. 0d0) then
3882            t1 = tm
3883            t2 = tm
3884         else
3885            t1 = tm
3886            t2 = tm
3887         end if
3888
3889         return
3890         end
3891
3892         subroutine isco2(i, j, allok, tm, t1, t2)
3893 c       this subroutine is used to decide whether there is a collision between
3894 c       particle i and j, if there is one allok=1, and tm gives the
3895 c       collision time, t1 the collision time for i,
3896 c       t2 the collision time for j
3897
3898         implicit double precision (a-h, o-z)
3899         parameter (MAXPTN=400001)
3900
3901         common /para2/ xmp, xmu, alpha, rscut2, cutof2
3902 cc      SAVE /para2/
3903         common /para5/ iconfg, iordsc
3904 cc      SAVE /para5/
3905         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
3906      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
3907      &       xmass(MAXPTN), ityp(MAXPTN)
3908 cc      SAVE /prec2/
3909         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
3910 cc      SAVE /prec4/
3911         common /ilist1/
3912      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
3913      &     ictype, icsta(MAXPTN),
3914      &     nic(MAXPTN), icels(MAXPTN)
3915 cc      SAVE /ilist1/
3916         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3917 cc      SAVE /ilist5/
3918         SAVE
3919
3920         logical allok
3921
3922 c       preventing consecutive collisions
3923         allok = last(i) .ne. j .or. last(j) .ne. i
3924
3925 c       set up numbers for later calculations
3926         i1 = i
3927         i2 = j
3928
3929         p4 = ft(i2) - ft(i1)
3930         p1 = gx(i2) - gx(i1)
3931         p2 = gy(i2) - gy(i1)
3932         p3 = gz(i2) - gz(i1)
3933
3934         q4 = e(i1)
3935         q1 = px(i1)
3936         q2 = py(i1)
3937         q3 = pz(i1)
3938
3939         r4 = e(i2)
3940         r1 = px(i2)
3941         r2 = py(i2)
3942         r3 = pz(i2)
3943
3944         a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
3945         b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
3946         c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
3947         d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
3948         ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
3949         f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
3950
3951 c       make sure particle 2 formed early
3952         h = a + b
3953         if (h .gt. 0d0) then
3954            g = a
3955            a = -b
3956            b = -g
3957
3958            g = c
3959            c = d
3960            d = g
3961
3962            i1 = j
3963            i2 = i
3964         end if
3965
3966 c       check the approaching criteria
3967         if (allok) then
3968
3969            vp = a * d - b * ee
3970
3971            allok = allok .and. vp .lt. 0d0
3972
3973         end if
3974
3975 c       check the closest approach distance criteria
3976          if (allok) then
3977
3978            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
3979      &          (ee ** 2 - c * d)
3980
3981            allok = allok .and. dm2 .lt. cutof2
3982
3983         end if
3984
3985 c       check the time criteria
3986         if (allok) then
3987
3988            tc1 = ft(i1) - e(i1) * (a * d - b * ee)/(ee ** 2 - c * d)
3989            tc2 = ft(i2) + e(i2) * (b * c - a * ee)/(ee ** 2 - c * d)
3990            if (iordsc .eq. 20) then
3991               tm = min(tc1, tc2)
3992            else if (iordsc .eq. 21) then
3993               tm = 0.5d0 * (tc1 + tc2)
3994            else
3995               tm = max(tc1, tc2)
3996            end if
3997
3998            allok = allok .and. tm .gt. ft(i) .and. tm .gt. ft(j)
3999
4000         end if
4001
4002 c        check rts cut
4003         if (allok) then
4004
4005            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
4006      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
4007
4008            allok = allok .and. rts2 .gt. rscut2
4009         end if
4010
4011         if (.not. allok) then
4012            tm = tlarge
4013            t1 = tlarge
4014            t2 = tlarge
4015         else if (h .gt. 0d0) then
4016            t1 = tc2
4017            t2 = tc1
4018         else
4019            t1 = tc1
4020            t2 = tc2
4021         end if
4022
4023         return
4024         end
4025
4026         subroutine isco3(i, j, allok, tm, t1, t2)
4027 c       this subroutine is used to decide whether there is a collision between
4028 c       particle i and j, if there is one allok=1, and tm gives the
4029 c       collision time, t1 the collision time for i,
4030 c       t2 the collision time for j
4031
4032         implicit double precision (a-h, o-z)
4033         parameter (MAXPTN=400001)
4034
4035         common /para2/ xmp, xmu, alpha, rscut2, cutof2
4036 cc      SAVE /para2/
4037         common /para5/ iconfg, iordsc
4038 cc      SAVE /para5/
4039         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4040      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4041      &       xmass(MAXPTN), ityp(MAXPTN)
4042 cc      SAVE /prec2/
4043         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4044 cc      SAVE /prec4/
4045         common /ilist1/
4046      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
4047      &     ictype, icsta(MAXPTN),
4048      &     nic(MAXPTN), icels(MAXPTN)
4049 cc      SAVE /ilist1/
4050         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4051 cc      SAVE /ilist5/
4052         SAVE
4053
4054         logical allok
4055
4056 c       preventing consecutive collisions
4057         allok = last(i) .ne. j .or. last(j) .ne. i
4058
4059         if (ft(i) .ge. ft(j)) then
4060            i1 = j
4061            i2 = i
4062         else
4063            i1 = i
4064            i2 = j
4065         end if
4066
4067         if (allok) then
4068
4069            t1 = ft(i1)
4070            vx1 = vx(i1)
4071            vy1 = vy(i1)
4072            vz1 = vz(i1)
4073
4074            t2 = ft(i2)
4075
4076            dvx = vx(i2) - vx1
4077            dvy = vy(i2) - vy1
4078            dvz = vz(i2) - vz1
4079
4080            dt = t2 - t1
4081
4082            dx = gx(i2) - gx(i1) - vx1 * dt
4083            dy = gy(i2) - gy(i1) - vy1 * dt
4084            dz = gz(i2) - gz(i1) - vz1 * dt
4085
4086            vp = dvx * dx + dvy * dy + dvz * dz
4087
4088            allok = allok .and. vp .lt. 0d0
4089
4090         end if
4091
4092         if (allok) then
4093
4094            v2= dvx * dvx + dvy * dvy + dvz * dvz
4095
4096            if (v2 .eq. 0d0) then
4097               tm = tlarge
4098            else
4099               tm = t2 - vp / v2
4100            end if
4101
4102 c       note now tm is the absolute time
4103
4104            allok = allok .and. tm .gt. t1 .and. tm .gt. t2
4105
4106         end if
4107
4108         if (allok) then
4109
4110            dgx = dx - dvx * t2
4111            dgy = dy - dvy * t2
4112            dgz = dz - dvz * t2
4113
4114            dm2 = - v2 * tm ** 2  + dgx * dgx + dgy * dgy + dgz * dgz
4115
4116            allok = allok .and. dm2 .lt. cutof2
4117
4118         end if
4119
4120         if (allok) then
4121
4122            e1 = e(i1)
4123            px1 = px(i1)
4124            py1 = py(i1)
4125            pz1 = pz(i1)
4126            e2 = e(i2)
4127            px2 = px(i2)
4128            py2 = py(i2)
4129            pz2 = pz(i2)
4130
4131            rts2 = (e1 + e2) ** 2 - (px1 + px2) ** 2
4132      &          - (py1 + py2) ** 2 - (pz1 + pz2) ** 2
4133
4134            allok = allok .and. rts2 .gt. rscut2
4135         end if
4136
4137         if (.not. allok) then
4138            tm = tlarge
4139            t1 = tlarge
4140            t2 = tlarge
4141         else
4142            t1 = tm
4143            t2 = tm
4144         end if
4145
4146         return
4147         end
4148
4149         subroutine isco4(i, j, allok, tm, t1, t2)
4150 c       this subroutine is used to decide whether there is a collision between
4151 c       particle i and j, if there is one allok=1, and tm gives the
4152 c       collision time, t1 the collision time for i,
4153 c       t2 the collision time for j
4154
4155         implicit double precision (a-h, o-z)
4156         parameter (MAXPTN=400001)
4157
4158         common /para2/ xmp, xmu, alpha, rscut2, cutof2
4159 cc      SAVE /para2/
4160         common /para5/ iconfg, iordsc
4161 cc      SAVE /para5/
4162         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4163      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4164      &       xmass(MAXPTN), ityp(MAXPTN)
4165 cc      SAVE /prec2/
4166         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4167 cc      SAVE /prec4/
4168         common /ilist1/
4169      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
4170      &     ictype, icsta(MAXPTN),
4171      &     nic(MAXPTN), icels(MAXPTN)
4172 cc      SAVE /ilist1/
4173         common /ilist3/ size1, size2, size3, v1, v2, v3, size
4174 cc      SAVE /ilist3/
4175         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4176 cc      SAVE /ilist5/
4177         SAVE
4178
4179         logical allok
4180
4181 c       preventing consecutive collisions
4182         allok = last(i) .ne. j .or. last(j) .ne. i
4183
4184 c       set up numbers for later calculations
4185
4186         icels1 = icels(i)
4187         ii1 = icels1 / 10000
4188         jj1 = (icels1 - ii1 * 10000) / 100
4189         kk1 = icels1 - ii1 * 10000 - jj1 * 100
4190         icels2 = icels(j)
4191         ii2 = icels2 / 10000
4192         jj2 = (icels2 - ii2 * 10000) / 100
4193         kk2 = icels2 - ii2 * 10000 - jj2 * 100
4194
4195         i1 = i
4196         i2 = j
4197
4198         p4 = ft(i2) - ft(i1)
4199         p1 = gx(i2) - gx(i1)
4200         p2 = gy(i2) - gy(i1)
4201         p3 = gz(i2) - gz(i1)
4202
4203         if (ii1 - ii2 .gt. 5) then
4204            p1 = p1 + 10d0 * size1
4205         else if (ii1 - ii2 .lt. -5) then
4206            p1 = p1 - 10d0 * size1
4207         end if
4208         if (jj1 - jj2 .gt. 5) then
4209            p2 = p2 + 10d0 * size2
4210         else if (jj1 - jj2 .lt. -5) then
4211            p2 = p2 - 10d0 * size2
4212         end if
4213         if (kk1 - kk2 .gt. 5) then
4214            p3 = p3 + 10d0 * size3
4215         else if (kk1 - kk2 .lt. -5) then
4216            p3 = p3 - 10d0 * size3
4217         end if
4218
4219         q4 = e(i1)
4220         q1 = px(i1)
4221         q2 = py(i1)
4222         q3 = pz(i1)
4223
4224         r4 = e(i2)
4225         r1 = px(i2)
4226         r2 = py(i2)
4227         r3 = pz(i2)
4228
4229         a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
4230         b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
4231         c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
4232         d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
4233         ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
4234         f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
4235
4236 c       make sure particle 2 formed early
4237         h = a + b
4238         if (h .gt. 0d0) then
4239            g = a
4240            a = -b
4241            b = -g
4242
4243            g = c
4244            c = d
4245            d = g
4246
4247            i1 = j
4248            i2 = i
4249         end if
4250
4251 c       check the approaching criteria
4252         if (allok) then
4253
4254            vp = a * d - b * ee
4255
4256            allok = allok .and. vp .lt. 0d0
4257
4258         end if
4259
4260 c       check the closest approach distance criteria
4261          if (allok) then
4262
4263            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
4264      &           (ee ** 2 - c * d)
4265
4266            allok = allok .and. dm2 .lt. cutof2
4267
4268         end if
4269
4270 c       check the time criteria
4271         if (allok) then
4272
4273            tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
4274            tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
4275            tm = 0.5d0 * (tc1 + tc2)
4276
4277            allok = allok .and. tm .gt. ft(i) .and. tm .gt. ft(j)
4278
4279         end if
4280
4281 c        check rts cut
4282         if (allok) then
4283
4284            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
4285      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
4286
4287            allok = allok .and. rts2 .gt. rscut2
4288         end if
4289
4290         if (.not. allok) then
4291            tm = tlarge
4292            t1 = tlarge
4293            t2 = tlarge
4294         else if (h .gt. 0d0) then
4295            t1 = tm
4296            t2 = tm
4297         else
4298            t1 = tm
4299            t2 = tm
4300         end if
4301
4302         return
4303         end
4304
4305         subroutine isco5(i, j, allok, tm, t1, t2)
4306 c       this subroutine is used to decide whether there is a collision between
4307 c       particle i and j, if there is one allok=1, and tm gives the
4308 c       collision time, t1 the collision time for i,
4309 c       t2 the collision time for j
4310
4311         implicit double precision (a-h, o-z)
4312         parameter (MAXPTN=400001)
4313
4314         common /para2/ xmp, xmu, alpha, rscut2, cutof2
4315 cc      SAVE /para2/
4316         common /para5/ iconfg, iordsc
4317 cc      SAVE /para5/
4318         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4319      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4320      &       xmass(MAXPTN), ityp(MAXPTN)
4321 cc      SAVE /prec2/
4322         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4323 cc      SAVE /prec4/
4324         common /ilist1/
4325      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
4326      &     ictype, icsta(MAXPTN),
4327      &     nic(MAXPTN), icels(MAXPTN)
4328 cc      SAVE /ilist1/
4329         common /ilist3/ size1, size2, size3, v1, v2, v3, size
4330 cc      SAVE /ilist3/
4331         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4332 cc      SAVE /ilist5/
4333         SAVE
4334
4335         logical allok
4336
4337 c       preventing consecutive collisions
4338         allok = last(i) .ne. j .or. last(j) .ne. i
4339
4340 c       set up numbers for later calculations
4341
4342         icels1 = icels(i)
4343         ii1 = icels1 / 10000
4344         jj1 = (icels1 - ii1 * 10000) / 100
4345         kk1 = icels1 - ii1 * 10000 - jj1 * 100
4346         icels2 = icels(j)
4347         ii2 = icels2 / 10000
4348         jj2 = (icels2 - ii2 * 10000) / 100
4349         kk2 = icels2 - ii2 * 10000 - jj2 * 100
4350
4351         i1 = i
4352         i2 = j
4353
4354         p4 = ft(i2) - ft(i1)
4355         p1 = gx(i2) - gx(i1)
4356         p2 = gy(i2) - gy(i1)
4357         p3 = gz(i2) - gz(i1)
4358
4359         if (ii1 - ii2 .gt. 5) then
4360            p1 = p1 + 10d0 * size1
4361         else if (ii1 - ii2 .lt. -5) then
4362            p1 = p1 - 10d0 * size1
4363         end if
4364         if (jj1 - jj2 .gt. 5) then
4365            p2 = p2 + 10d0 * size2
4366         else if (jj1 - jj2 .lt. -5) then
4367            p2 = p2 - 10d0 * size2
4368         end if
4369         if (kk1 - kk2 .gt. 5) then
4370            p3 = p3 + 10d0 * size3
4371         else if (kk1 - kk2 .lt. -5) then
4372            p3 = p3 - 10d0 * size3
4373         end if
4374
4375         q4 = e(i1)
4376         q1 = px(i1)
4377         q2 = py(i1)
4378         q3 = pz(i1)
4379
4380         r4 = e(i2)
4381         r1 = px(i2)
4382         r2 = py(i2)
4383         r3 = pz(i2)
4384
4385         a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
4386         b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
4387         c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
4388         d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
4389         ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
4390         f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
4391
4392 c       make sure particle 2 formed early
4393         h = a + b
4394         if (h .gt. 0d0) then
4395            g = a
4396            a = -b
4397            b = -g
4398
4399            g = c
4400            c = d
4401            d = g
4402
4403            i1 = j
4404            i2 = i
4405         end if
4406
4407 c       check the approaching criteria
4408         if (allok) then
4409
4410            vp = a * d - b * ee
4411
4412            allok = allok .and. vp .lt. 0d0
4413
4414         end if
4415
4416 c       check the closest approach distance criteria
4417          if (allok) then
4418
4419            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
4420      &           (ee ** 2 - c * d)
4421
4422            allok = allok .and. dm2 .lt. cutof2
4423
4424         end if
4425
4426 c       check the time criteria
4427         if (allok) then
4428
4429            tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
4430            tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
4431            if (iordsc .eq. 20) then
4432               tm = min(tc1, tc2)
4433            else if (iordsc .eq. 21) then
4434               tm = 0.5d0 * (tc1 + tc2)
4435            else
4436               tm = max(tc1, tc2)
4437            end if
4438
4439            allok = allok .and. tm .gt. ft(i) .and. tm .gt. ft(j)
4440
4441         end if
4442
4443 c        check rts cut
4444         if (allok) then
4445
4446            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
4447      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
4448
4449            allok = allok .and. rts2 .gt. rscut2
4450         end if
4451
4452         if (.not. allok) then
4453            tm = tlarge
4454            t1 = tlarge
4455            t2 = tlarge
4456         else if (h .gt. 0d0) then
4457            t1 = tc2
4458            t2 = tc1
4459         else
4460            t1 = tc1
4461            t2 = tc2
4462         end if
4463
4464         return
4465         end
4466
4467         subroutine isco6(i, j, allok, tm, t1, t2)
4468 c       this subroutine is used to decide whether there is a collision between
4469 c       particle i and j, if there is one allok=1, and tm gives the
4470 c       collision time, t1 the collision time for i,
4471 c       t2 the collision time for j
4472
4473         implicit double precision (a-h, o-z)
4474         parameter (MAXPTN=400001)
4475
4476         common /para2/ xmp, xmu, alpha, rscut2, cutof2
4477 cc      SAVE /para2/
4478         common /para5/ iconfg, iordsc
4479 cc      SAVE /para5/
4480         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4481      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4482      &       xmass(MAXPTN), ityp(MAXPTN)
4483 cc      SAVE /prec2/
4484         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4485 cc      SAVE /prec4/
4486         common /ilist1/
4487      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
4488      &     ictype, icsta(MAXPTN),
4489      &     nic(MAXPTN), icels(MAXPTN)
4490 cc      SAVE /ilist1/
4491         common /ilist3/ size1, size2, size3, v1, v2, v3, size
4492 cc      SAVE /ilist3/
4493         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4494 cc      SAVE /ilist5/
4495         SAVE
4496
4497         logical allok
4498
4499 c       preventing consecutive collisions
4500         allok = last(i) .ne. j .or. last(j) .ne. i
4501
4502         if (ft(i) .ge. ft(j)) then
4503            i1 = j
4504            i2 = i
4505         else
4506            i1 = i
4507            i2 = j
4508         end if
4509
4510         icels1 = icels(i1)
4511         ii1 = icels1 / 10000
4512         jj1 = (icels1 - ii1 * 10000) / 100
4513         kk1 = icels1 - ii1 * 10000 - jj1 * 100
4514         icels2 = icels(i2)
4515         ii2 = icels2 / 10000
4516         jj2 = (icels2 - ii2 * 10000) / 100
4517         kk2 = icels2 - ii2 * 10000 - jj2 * 100
4518
4519         if (allok) then
4520
4521            t1 = ft(i1)
4522            vx1 = vx(i1)
4523            vy1 = vy(i1)
4524            vz1 = vz(i1)
4525
4526            t2 = ft(i2)
4527
4528            dvx = vx(i2) - vx1
4529            dvy = vy(i2) - vy1
4530            dvz = vz(i2) - vz1
4531
4532            dt = t2 - t1
4533
4534            dx = gx(i2) - gx(i1) - vx1 * dt
4535            dy = gy(i2) - gy(i1) - vy1 * dt
4536            dz = gz(i2) - gz(i1) - vz1 * dt
4537
4538            if (ii1 - ii2 .gt. 5) then
4539               dx = dx + 10d0 * size1
4540            else if (ii1 - ii2 .lt. - 5) then
4541               dx = dx - 10d0 * size1
4542            end if
4543
4544            if (jj1 - jj2 .gt. 5) then
4545               dy = dy + 10d0 * size2
4546            else if (jj1 - jj2 .lt. - 5) then
4547               dy = dy - 10d0 * size2
4548            end if
4549
4550            if (kk1 - kk2 .gt. 5) then
4551               dz = dz + 10d0 * size3
4552            else if (kk1 - kk2 .lt. -5) then
4553               dz = dz - 10d0 * size3
4554            end if
4555
4556            vp = dvx * dx + dvy * dy + dvz * dz
4557
4558            allok = allok .and. vp .lt. 0d0
4559
4560         end if
4561
4562         if (allok) then
4563
4564            v2p = dvx * dvx + dvy * dvy + dvz * dvz
4565
4566            if (v2p .eq. 0d0) then
4567               tm = tlarge
4568            else
4569               tm = t2 - vp / v2p
4570            end if
4571
4572 c       note now tm is the absolute time
4573
4574            allok = allok .and. tm .gt. t1 .and. tm .gt. t2
4575
4576         end if
4577
4578         if (allok) then
4579
4580            dgx = dx - dvx * t2
4581            dgy = dy - dvy * t2
4582            dgz = dz - dvz * t2
4583
4584            dm2 = - v2p * tm ** 2  + dgx * dgx + dgy * dgy + dgz * dgz
4585
4586            allok = allok .and. dm2 .lt. cutof2
4587
4588         end if
4589
4590         if (allok) then
4591
4592            e1 = e(i1)
4593            px1 = px(i1)
4594            py1 = py(i1)
4595            pz1 = pz(i1)
4596            e2 = e(i2)
4597            px2 = px(i2)
4598            py2 = py(i2)
4599            pz2 = pz(i2)
4600
4601            rts2 = (e1 + e2) ** 2 - (px1 + px2) ** 2
4602      &          - (py1 + py2) ** 2 - (pz1 + pz2) ** 2
4603
4604            allok = allok .and. rts2 .gt. rscut2
4605         end if
4606
4607         if (.not. allok) then
4608            tm = tlarge
4609            t1 = tlarge
4610            t2 = tlarge
4611         else
4612            t1 = tm
4613            t2 = tm
4614         end if
4615
4616         return
4617         end
4618
4619         subroutine isco7(i, j, allok, tm, t1, t2)
4620 c       this subroutine is used to decide whether there is a collision between
4621 c       particle i and j, if there is one allok=1, and tm gives the
4622 c       collision time, t1 the collision time for i,
4623 c       t2 the collision time for j
4624
4625         implicit double precision (a-h, o-z)
4626         parameter (MAXPTN=400001)
4627
4628         common /para2/ xmp, xmu, alpha, rscut2, cutof2
4629 cc      SAVE /para2/
4630         common /para5/ iconfg, iordsc
4631 cc      SAVE /para5/
4632         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4633      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4634      &       xmass(MAXPTN), ityp(MAXPTN)
4635 cc      SAVE /prec2/
4636         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4637 cc      SAVE /prec4/
4638         common /aurec1/ jxa, jya, jza
4639 cc      SAVE /aurec1/
4640         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
4641 cc      SAVE /aurec2/
4642         common /ilist1/
4643      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
4644      &     ictype, icsta(MAXPTN),
4645      &     nic(MAXPTN), icels(MAXPTN)
4646 cc      SAVE /ilist1/
4647         common /ilist3/ size1, size2, size3, v1, v2, v3, size
4648 cc      SAVE /ilist3/
4649         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4650 cc      SAVE /ilist5/
4651         SAVE
4652
4653         logical allok, allokp
4654
4655 c       preventing consecutive collisions
4656         allok = last(i) .ne. j .or. last(j) .ne. i
4657
4658 c       set up numbers for later calculations
4659
4660         tm = tlarge
4661
4662         if (allok) then
4663            do 1000 ii = - 1, 1
4664               do 2000 jj = - 1, 1
4665
4666                  allokp = .true.
4667
4668                  i1 = i
4669                  i2 = j
4670
4671                  p4 = ft(j) - ft(i)
4672                  p1 = gx(j) - gx(i)
4673                  p2 = gy(j) - gy(i)
4674                  p3 = gz(j) - gz(i)
4675
4676                  p1 = p1 + ii * 10d0 * size1
4677                  p2 = p2 + jj * 10d0 * size2
4678
4679                  q4 = e(i)
4680                  q1 = px(i)
4681                  q2 = py(i)
4682                  q3 = pz(i)
4683
4684                  r4 = e(j)
4685                  r1 = px(j)
4686                  r2 = py(j)
4687                  r3 = pz(j)
4688
4689                  a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
4690                  b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
4691                  c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
4692                  d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
4693                  ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
4694                  f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
4695
4696 c       make sure particle 2 formed early
4697                  h = a + b
4698                  if (h .gt. 0d0) then
4699                     g = a
4700                     a = -b
4701                     b = -g
4702                     g = c
4703                     c = d
4704                     d = g
4705                     i1 = j
4706                     i2 = i
4707                  end if
4708
4709 c       check the approaching criteria
4710                  if (allokp) then
4711                     vp = a * d - b * ee
4712                     allokp = allokp .and. vp .lt. 0d0
4713                  end if
4714
4715 c       check the closest approach distance criteria
4716                  if (allokp) then
4717            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
4718      &            (ee ** 2 - c * d)
4719                     allokp = allokp .and. dm2 .lt. cutof2
4720                  end if
4721
4722 c       check the time criteria
4723                  if (allokp) then
4724            tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
4725            tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
4726            tmp = 0.5d0 * (tc1 + tc2)
4727            allokp = allokp .and. tmp .gt. ft(i) .and. tmp .gt. ft(j)
4728                  end if
4729
4730                  if (allokp .and. tmp .lt. tm) then
4731                     tm = tmp
4732                     jxa = ii
4733                     jya = jj
4734 cd                    dgxa(j) = ii * 10d0 * size1
4735 cd                    dgya(j) = jj * 10d0 * size2
4736 cd                    dgxa(i) = - dgxa(j)
4737 cd                    dgya(i) = - dgya(j)
4738                  end if
4739
4740  2000              continue
4741  1000           continue
4742
4743            if (tm .eq. tlarge) then
4744               allok = .false.
4745            end if
4746
4747         end if
4748
4749 c        check rts cut
4750         if (allok) then
4751
4752            q4 = e(i1)
4753            q1 = px(i1)
4754            q2 = py(i1)
4755            q3 = pz(i1)
4756
4757            r4 = e(i2)
4758            r1 = px(i2)
4759            r2 = py(i2)
4760            r3 = pz(i2)
4761
4762            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
4763      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
4764
4765            allok = allok .and. rts2 .gt. rscut2
4766         end if
4767
4768         if (.not. allok) then
4769            tm = tlarge
4770            t1 = tlarge
4771            t2 = tlarge
4772         else if (h .gt. 0d0) then
4773            t1 = tm
4774            t2 = tm
4775         else
4776            t1 = tm
4777            t2 = tm
4778         end if
4779
4780         return
4781         end
4782
4783         subroutine isco8(i, j, allok, tm, t1, t2)
4784 c       this subroutine is used to decide whether there is a collision between
4785 c       particle i and j, if there is one allok=1, and tm gives the
4786 c       collision time, t1 the collision time for i,
4787 c       t2 the collision time for j
4788
4789         implicit double precision (a-h, o-z)
4790         parameter (MAXPTN=400001)
4791
4792         common /para2/ xmp, xmu, alpha, rscut2, cutof2
4793 cc      SAVE /para2/
4794         common /para5/ iconfg, iordsc
4795 cc      SAVE /para5/
4796         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4797      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4798      &       xmass(MAXPTN), ityp(MAXPTN)
4799 cc      SAVE /prec2/
4800         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4801 cc      SAVE /prec4/
4802         common /aurec1/ jxa, jya, jza
4803 cc      SAVE /aurec1/
4804         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
4805 cc      SAVE /aurec2/
4806         common /ilist1/
4807      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
4808      &     ictype, icsta(MAXPTN),
4809      &     nic(MAXPTN), icels(MAXPTN)
4810 cc      SAVE /ilist1/
4811         common /ilist3/ size1, size2, size3, v1, v2, v3, size
4812 cc      SAVE /ilist3/
4813         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4814 cc      SAVE /ilist5/
4815         SAVE
4816
4817         logical allok, allokp
4818
4819 c       preventing consecutive collisions
4820         allok = last(i) .ne. j .or. last(j) .ne. i
4821
4822 c       set up numbers for later calculations
4823
4824         tm = tlarge
4825
4826         if (allok) then
4827            do 1000 ii = - 1, 1
4828               do 2000 jj = - 1, 1
4829
4830                  allokp = .true.
4831
4832                  i1 = i
4833                  i2 = j
4834
4835                  p4 = ft(j) - ft(i)
4836                  p1 = gx(j) - gx(i)
4837                  p2 = gy(j) - gy(i)
4838                  p3 = gz(j) - gz(i)
4839
4840                  p1 = p1 + ii * 10d0 * size1
4841                  p2 = p2 + jj * 10d0 * size2
4842
4843                  q4 = e(i)
4844                  q1 = px(i)
4845                  q2 = py(i)
4846                  q3 = pz(i)
4847
4848                  r4 = e(j)
4849                  r1 = px(j)
4850                  r2 = py(j)
4851                  r3 = pz(j)
4852
4853                  a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
4854                  b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
4855                  c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
4856                  d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
4857                  ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
4858                  f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
4859
4860 c       make sure particle 2 formed early
4861                  h = a + b
4862                  if (h .gt. 0d0) then
4863                     g = a
4864                     a = -b
4865                     b = -g
4866                     g = c
4867                     c = d
4868                     d = g
4869                     i1 = j
4870                     i2 = i
4871                  end if
4872
4873 c       check the approaching criteria
4874                  if (allokp) then
4875                     vp = a * d - b * ee
4876                     allokp = allokp .and. vp .lt. 0d0
4877                  end if
4878
4879 c       check the closest approach distance criteria
4880                  if (allokp) then
4881            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
4882      &            (ee ** 2 - c * d)
4883                     allokp = allokp .and. dm2 .lt. cutof2
4884                  end if
4885
4886 c       check the time criteria
4887                  if (allokp) then
4888            tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
4889            tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
4890                     if (iordsc .eq. 20) then
4891                        tmp = min(tc1, tc2)
4892                     else if (iordsc .eq. 21) then
4893                        tmp = 0.5d0 * (tc1 + tc2)
4894                     else
4895                        tmp = max(tc1, tc2)
4896                     end if
4897            allokp = allokp .and. tmp .gt. ft(i) .and. tmp .gt. ft(j)
4898                  end if
4899
4900                  if (allokp .and. tmp .lt. tm) then
4901                     tm = tmp
4902                     jxa = ii
4903                     jya = jj
4904                     ha = h
4905                     tc1a = tc1
4906                     tc2a = tc2
4907 cd                    dgxa(j) = ii * 10d0 * size1
4908 cd                    dgya(j) = jj * 10d0 * size2
4909 cd                    dgxa(i) = - dgxa(j)
4910 cd                    dgya(i) = - dgya(j)
4911                  end if
4912
4913  2000              continue
4914  1000           continue
4915
4916            if (tm .eq. tlarge) then
4917               allok = .false.
4918            end if
4919
4920         end if
4921
4922 c        check rts cut
4923         if (allok) then
4924
4925            q4 = e(i1)
4926            q1 = px(i1)
4927            q2 = py(i1)
4928            q3 = pz(i1)
4929
4930            r4 = e(i2)
4931            r1 = px(i2)
4932            r2 = py(i2)
4933            r3 = pz(i2)
4934
4935            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
4936      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
4937
4938            allok = allok .and. rts2 .gt. rscut2
4939         end if
4940
4941         if (.not. allok) then
4942            tm = tlarge
4943            t1 = tlarge
4944            t2 = tlarge
4945         else if (ha .gt. 0d0) then
4946            t1 = tc2a
4947            t2 = tc1a
4948         else
4949            t1 = tc1a
4950            t2 = tc2a
4951         end if
4952
4953         return
4954         end
4955
4956         subroutine isco9(i, j, allok, tm, t1, t2)
4957 c       this subroutine is used to decide whether there is a collision between
4958 c       particle i and j, if there is one allok=1, and tm gives the
4959 c       collision time, t1 the collision time for i,
4960 c       t2 the collision time for j
4961
4962         implicit double precision (a-h, o-z)
4963         parameter (MAXPTN=400001)
4964
4965         common /para2/ xmp, xmu, alpha, rscut2, cutof2
4966 cc      SAVE /para2/
4967         common /para5/ iconfg, iordsc
4968 cc      SAVE /para5/
4969         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4970      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4971      &       xmass(MAXPTN), ityp(MAXPTN)
4972 cc      SAVE /prec2/
4973         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4974 cc      SAVE /prec4/
4975         common /aurec1/ jxa, jya, jza
4976 cc      SAVE /aurec1/
4977         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
4978 cc      SAVE /aurec2/
4979         common /ilist1/
4980      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
4981      &     ictype, icsta(MAXPTN),
4982      &     nic(MAXPTN), icels(MAXPTN)
4983 cc      SAVE /ilist1/
4984         common /ilist3/ size1, size2, size3, v1, v2, v3, size
4985 cc      SAVE /ilist3/
4986         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4987 cc      SAVE /ilist5/
4988         SAVE
4989
4990         logical allok, allokp
4991
4992 c       preventing consecutive collisions
4993         allok = last(i) .ne. j .or. last(j) .ne. i
4994
4995         if (ft(i) .ge. ft(j)) then
4996            i1 = j
4997            i2 = i
4998            isign = -1
4999         else
5000            i1 = i
5001            i2 = j
5002            isign = 1
5003         end if
5004
5005         if (allok) then
5006            tm = tlarge
5007
5008            t1 = ft(i1)
5009            vx1 = vx(i1)
5010            vy1 = vy(i1)
5011            vz1 = vz(i1)
5012
5013            t2 = ft(i2)
5014
5015            dvx = vx(i2) - vx1
5016            dvy = vy(i2) - vy1
5017            dvz = vz(i2) - vz1
5018
5019            dt = t2 - t1
5020
5021            do 1000 ii = - 1, 1
5022               do 2000 jj = - 1, 1
5023
5024                  allokp = .true.
5025
5026                  dx = gx(i2) - gx(i1) - vx1 * dt
5027                  dy = gy(i2) - gy(i1) - vy1 * dt
5028                  dz = gz(i2) - gz(i1) - vz1 * dt
5029
5030                  dx = dx + ii * 10d0 * size1
5031                  dy = dy + jj * 10d0 * size2
5032
5033                  vp = dvx * dx + dvy * dy + dvz * dz
5034
5035                  allokp = allokp .and. vp .lt. 0d0
5036
5037                  if (allokp) then
5038
5039                     v2 = dvx * dvx + dvy * dvy + dvz * dvz
5040
5041                     if (v2 .eq. 0d0) then
5042                        tmp = tlarge
5043                     else
5044                        tmp = t2 - vp / v2
5045                     end if
5046
5047 c       note now tm is the absolute time
5048
5049                     allokp = allokp .and. tmp .gt. t1 .and.
5050      &                         tmp .gt. t2
5051
5052                  end if
5053
5054                  if (allokp) then
5055
5056                     dgx = dx - dvx * t2
5057                     dgy = dy - dvy * t2
5058                     dgz = dz - dvz * t2
5059
5060                     dm2 = - v2 * tmp ** 2  + dgx * dgx +
5061      &                    dgy * dgy + dgz * dgz
5062
5063                     allokp = allokp .and. dm2 .lt. cutof2
5064
5065                  end if
5066
5067                  if (allokp .and. tmp .lt. tm) then
5068                     tm = tmp
5069                     jxa = isign * ii
5070                     jya = isign * jj
5071                  end if
5072
5073  2000              continue
5074  1000           continue
5075
5076            if (tm .eq. tlarge) then
5077               allok = .false.
5078            end if
5079         end if
5080
5081         if (allok) then
5082
5083            e1 = e(i1)
5084            px1 = px(i1)
5085            py1 = py(i1)
5086            pz1 = pz(i1)
5087            e2 = e(i2)
5088            px2 = px(i2)
5089            py2 = py(i2)
5090            pz2 = pz(i2)
5091
5092            rts2 = (e1 + e2) ** 2 - (px1 + px2) ** 2
5093      &          - (py1 + py2) ** 2 - (pz1 + pz2) ** 2
5094
5095            allok = allok .and. rts2 .gt. rscut2
5096         end if
5097
5098         if (.not. allok) then
5099            tm = tlarge
5100            t1 = tlarge
5101            t2 = tlarge
5102         else
5103            t1 = tm
5104            t2 = tm
5105         end if
5106
5107         return
5108         end
5109
5110         subroutine isco10(i, j, allok, tm, t1, t2)
5111 c       this subroutine is used to decide whether there is a collision between
5112 c       particle i and j, if there is one allok=1, and tm gives the
5113 c       collision time, t1 the collision time for i,
5114 c       t2 the collision time for j
5115
5116         implicit double precision (a-h, o-z)
5117         parameter (MAXPTN=400001)
5118
5119         common /para2/ xmp, xmu, alpha, rscut2, cutof2
5120 cc      SAVE /para2/
5121         common /para5/ iconfg, iordsc
5122 cc      SAVE /para5/
5123         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
5124      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
5125      &       xmass(MAXPTN), ityp(MAXPTN)
5126 cc      SAVE /prec2/
5127         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
5128 cc      SAVE /prec4/
5129         common /aurec1/ jxa, jya, jza
5130 cc      SAVE /aurec1/
5131         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
5132 cc      SAVE /aurec2/
5133         common /ilist1/
5134      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
5135      &     ictype, icsta(MAXPTN),
5136      &     nic(MAXPTN), icels(MAXPTN)
5137 cc      SAVE /ilist1/
5138         common /ilist3/ size1, size2, size3, v1, v2, v3, size
5139 cc      SAVE /ilist3/
5140         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
5141 cc      SAVE /ilist5/
5142         SAVE
5143
5144         logical allok, allokp
5145
5146 c       preventing consecutive collisions
5147         allok = last(i) .ne. j .or. last(j) .ne. i
5148
5149 c       set up numbers for later calculations
5150
5151         tm = tlarge
5152
5153         if (allok) then
5154            do 1000 ii = - 1, 1
5155               do 2000 jj = - 1, 1
5156                  do 3000 kk = -1, 1
5157                  allokp = .true.
5158
5159                  i1 = i
5160                  i2 = j
5161
5162                  p4 = ft(j) - ft(i)
5163                  p1 = gx(j) - gx(i)
5164                  p2 = gy(j) - gy(i)
5165                  p3 = gz(j) - gz(i)
5166
5167                  p1 = p1 + ii * 10d0 * size1
5168                  p2 = p2 + jj * 10d0 * size2
5169                  p3 = p3 + kk * 10d0 * size3
5170
5171                  q4 = e(i)
5172                  q1 = px(i)
5173                  q2 = py(i)
5174                  q3 = pz(i)
5175
5176                  r4 = e(j)
5177                  r1 = px(j)
5178                  r2 = py(j)
5179                  r3 = pz(j)
5180
5181                  a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
5182                  b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
5183                  c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
5184                  d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
5185                  ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
5186                  f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
5187
5188 c       make sure particle 2 formed early
5189                  h = a + b
5190                  if (h .gt. 0d0) then
5191                     g = a
5192                     a = -b
5193                     b = -g
5194                     g = c
5195                     c = d
5196                     d = g
5197                     i1 = j
5198                     i2 = i
5199                  end if
5200
5201 c       check the approaching criteria
5202                  if (allokp) then
5203                     vp = a * d - b * ee
5204                     allokp = allokp .and. vp .lt. 0d0
5205                  end if
5206
5207 c       check the closest approach distance criteria
5208                  if (allokp) then
5209            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
5210      &            (ee ** 2 - c * d)
5211                     allokp = allokp .and. dm2 .lt. cutof2
5212                  end if
5213
5214 c       check the time criteria
5215                  if (allokp) then
5216            tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
5217            tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
5218            tmp = 0.5d0 * (tc1 + tc2)
5219            allokp = allokp .and. tmp .gt. ft(i) .and. tmp .gt. ft(j)
5220                  end if
5221
5222                  if (allokp .and. tmp .lt. tm) then
5223                     tm = tmp
5224                     jxa = ii
5225                     jya = jj
5226                     jza = kk
5227 cd                    dgxa(j) = ii * 10d0 * size1
5228 cd                    dgya(j) = jj * 10d0 * size2
5229 cd                    dgxa(i) = - dgxa(j)
5230 cd                    dgya(i) = - dgya(j)
5231                  end if
5232
5233  3000                 continue
5234  2000              continue
5235  1000           continue
5236
5237            if (tm .eq. tlarge) then
5238               allok = .false.
5239            end if
5240
5241         end if
5242
5243 c        check rts cut
5244         if (allok) then
5245
5246            q4 = e(i1)
5247            q1 = px(i1)
5248            q2 = py(i1)
5249            q3 = pz(i1)
5250
5251            r4 = e(i2)
5252            r1 = px(i2)
5253            r2 = py(i2)
5254            r3 = pz(i2)
5255
5256            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
5257      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
5258
5259            allok = allok .and. rts2 .gt. rscut2
5260         end if
5261
5262         if (.not. allok) then
5263            tm = tlarge
5264            t1 = tlarge
5265            t2 = tlarge
5266         else if (h .gt. 0d0) then
5267            t1 = tm
5268            t2 = tm
5269         else
5270            t1 = tm
5271            t2 = tm
5272         end if
5273
5274         return
5275         end
5276
5277         subroutine isco11(i, j, allok, tm, t1, t2)
5278 c       this subroutine is used to decide whether there is a collision between
5279 c       particle i and j, if there is one allok=1, and tm gives the
5280 c       collision time, t1 the collision time for i,
5281 c       t2 the collision time for j
5282
5283         implicit double precision (a-h, o-z)
5284         parameter (MAXPTN=400001)
5285
5286         common /para2/ xmp, xmu, alpha, rscut2, cutof2
5287 cc      SAVE /para2/
5288         common /para5/ iconfg, iordsc
5289 cc      SAVE /para5/
5290         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
5291      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
5292      &       xmass(MAXPTN), ityp(MAXPTN)
5293 cc      SAVE /prec2/
5294         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
5295 cc      SAVE /prec4/
5296         common /aurec1/ jxa, jya, jza
5297 cc      SAVE /aurec1/
5298         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
5299 cc      SAVE /aurec2/
5300         common /ilist1/
5301      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
5302      &     ictype, icsta(MAXPTN),
5303      &     nic(MAXPTN), icels(MAXPTN)
5304 cc      SAVE /ilist1/
5305         common /ilist3/ size1, size2, size3, v1, v2, v3, size
5306 cc      SAVE /ilist3/
5307         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
5308 cc      SAVE /ilist5/
5309         SAVE
5310
5311         logical allok, allokp
5312
5313 c       preventing consecutive collisions
5314         allok = last(i) .ne. j .or. last(j) .ne. i
5315
5316 c       set up numbers for later calculations
5317
5318         tm = tlarge
5319
5320         if (allok) then
5321            do 1000 ii = - 1, 1
5322               do 2000 jj = - 1, 1
5323                  do 3000 kk = - 1, 1
5324
5325                  allokp = .true.
5326
5327                  i1 = i
5328                  i2 = j
5329
5330                  p4 = ft(j) - ft(i)
5331                  p1 = gx(j) - gx(i)
5332                  p2 = gy(j) - gy(i)
5333                  p3 = gz(j) - gz(i)
5334
5335                  p1 = p1 + ii * 10d0 * size1
5336                  p2 = p2 + jj * 10d0 * size2
5337                  p3 = p3 + kk * 10d0 * size3
5338
5339                  q4 = e(i)
5340                  q1 = px(i)
5341                  q2 = py(i)
5342                  q3 = pz(i)
5343
5344                  r4 = e(j)
5345                  r1 = px(j)
5346                  r2 = py(j)
5347                  r3 = pz(j)
5348
5349                  a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
5350                  b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
5351                  c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
5352                  d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
5353                  ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
5354                  f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
5355
5356 c       make sure particle 2 formed early
5357                  h = a + b
5358                  if (h .gt. 0d0) then
5359                     g = a
5360                     a = -b
5361                     b = -g
5362                     g = c
5363                     c = d
5364                     d = g
5365                     i1 = j
5366                     i2 = i
5367                  end if
5368
5369 c       check the approaching criteria
5370                  if (allokp) then
5371                     vp = a * d - b * ee
5372                     allokp = allokp .and. vp .lt. 0d0
5373                  end if
5374
5375 c       check the closest approach distance criteria
5376                  if (allokp) then
5377            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
5378      &            (ee ** 2 - c * d)
5379                     allokp = allokp .and. dm2 .lt. cutof2
5380                  end if
5381
5382 c       check the time criteria
5383                  if (allokp) then
5384            tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
5385            tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
5386                     if (iordsc .eq. 20) then
5387                        tmp = min(tc1, tc2)
5388                     else if (iordsc .eq. 21) then
5389                        tmp = 0.5d0 * (tc1 + tc2)
5390                     else
5391                        tmp = max(tc1, tc2)
5392                     end if
5393            allokp = allokp .and. tmp .gt. ft(i) .and. tmp .gt. ft(j)
5394                  end if
5395
5396                  if (allokp .and. tmp .lt. tm) then
5397                     tm = tmp
5398                     jxa = ii
5399                     jya = jj
5400                     jza = kk
5401                     ha = h
5402                     tc1a = tc1
5403                     tc2a = tc2
5404 cd                    dgxa(j) = ii * 10d0 * size1
5405 cd                    dgya(j) = jj * 10d0 * size2
5406 cd                    dgxa(i) = - dgxa(j)
5407 cd                    dgya(i) = - dgya(j)
5408                  end if
5409
5410  3000                 continue
5411  2000              continue
5412  1000           continue
5413
5414            if (tm .eq. tlarge) then
5415               allok = .false.
5416            end if
5417
5418         end if
5419
5420 c        check rts cut
5421         if (allok) then
5422
5423            q4 = e(i1)
5424            q1 = px(i1)
5425            q2 = py(i1)
5426            q3 = pz(i1)
5427
5428            r4 = e(i2)
5429            r1 = px(i2)
5430            r2 = py(i2)
5431            r3 = pz(i2)
5432
5433            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
5434      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
5435
5436            allok = allok .and. rts2 .gt. rscut2
5437         end if
5438
5439         if (.not. allok) then
5440            tm = tlarge
5441            t1 = tlarge
5442            t2 = tlarge
5443         else if (ha .gt. 0d0) then
5444            t1 = tc2a
5445            t2 = tc1a
5446         else
5447            t1 = tc1a
5448            t2 = tc2a
5449         end if
5450
5451         return
5452         end
5453
5454         subroutine isco12(i, j, allok, tm, t1, t2)
5455 c       this subroutine is used to decide whether there is a collision between
5456 c       particle i and j, if there is one allok=1, and tm gives the
5457 c       collision time, t1 the collision time for i,
5458 c       t2 the collision time for j
5459
5460         implicit double precision (a-h, o-z)
5461         parameter (MAXPTN=400001)
5462
5463         common /para2/ xmp, xmu, alpha, rscut2, cutof2
5464 cc      SAVE /para2/
5465         common /para5/ iconfg, iordsc
5466 cc      SAVE /para5/
5467         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
5468      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
5469      &       xmass(MAXPTN), ityp(MAXPTN)
5470 cc      SAVE /prec2/
5471         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
5472 cc      SAVE /prec4/
5473         common /aurec1/ jxa, jya, jza
5474 cc      SAVE /aurec1/
5475         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
5476 cc      SAVE /aurec2/
5477         common /ilist1/
5478      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
5479      &     ictype, icsta(MAXPTN),
5480      &     nic(MAXPTN), icels(MAXPTN)
5481 cc      SAVE /ilist1/
5482         common /ilist3/ size1, size2, size3, v1, v2, v3, size
5483 cc      SAVE /ilist3/
5484         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
5485 cc      SAVE /ilist5/
5486         SAVE
5487
5488         logical allok, allokp
5489
5490 c       preventing consecutive collisions
5491         allok = last(i) .ne. j .or. last(j) .ne. i
5492
5493         if (ft(i) .ge. ft(j)) then
5494            i1 = j
5495            i2 = i
5496            isign = -1
5497         else
5498            i1 = i
5499            i2 = j
5500            isign = 1
5501         end if
5502
5503         if (allok) then
5504            tm = tlarge
5505
5506            t1 = ft(i1)
5507            vx1 = vx(i1)
5508            vy1 = vy(i1)
5509            vz1 = vz(i1)
5510
5511            t2 = ft(i2)
5512
5513            dvx = vx(i2) - vx1
5514            dvy = vy(i2) - vy1
5515            dvz = vz(i2) - vz1
5516
5517            dt = t2 - t1
5518
5519            do 1000 ii = - 1, 1
5520               do 2000 jj = - 1, 1
5521                  do 3000 kk = -1, 1
5522
5523                  allokp = .true.
5524
5525                  dx = gx(i2) - gx(i1) - vx1 * dt
5526                  dy = gy(i2) - gy(i1) - vy1 * dt
5527                  dz = gz(i2) - gz(i1) - vz1 * dt
5528
5529                  dx = dx + ii * 10d0 * size1
5530                  dy = dy + jj * 10d0 * size2
5531                  dz = dz + kk * 10d0 * size3
5532
5533                  vp = dvx * dx + dvy * dy + dvz * dz
5534
5535                  allokp = allokp .and. vp .lt. 0d0
5536
5537                  if (allokp) then
5538
5539                     v2 = dvx * dvx + dvy * dvy + dvz * dvz
5540
5541                     if (v2 .eq. 0d0) then
5542                        tmp = tlarge
5543                     else
5544                        tmp = t2 - vp / v2
5545                     end if
5546
5547 c       note now tm is the absolute time
5548
5549                     allokp = allokp .and. tmp .gt. t1 .and.
5550      &                         tmp .gt. t2
5551
5552                  end if
5553
5554                  if (allokp) then
5555
5556                     dgx = dx - dvx * t2
5557                     dgy = dy - dvy * t2
5558                     dgz = dz - dvz * t2
5559
5560                     dm2 = - v2 * tmp ** 2  + dgx * dgx +
5561      &                    dgy * dgy + dgz * dgz
5562
5563                     allokp = allokp .and. dm2 .lt. cutof2
5564
5565                  end if
5566
5567                  if (allokp .and. tmp .lt. tm) then
5568                     tm = tmp
5569                     jxa = isign * ii
5570                     jya = isign * jj
5571                     jza = isign * kk
5572                  end if
5573
5574  3000                 continue
5575  2000              continue
5576  1000           continue
5577
5578            if (tm .eq. tlarge) then
5579               allok = .false.
5580            end if
5581         end if
5582
5583         if (allok) then
5584
5585            e1 = e(i1)
5586            px1 = px(i1)
5587            py1 = py(i1)
5588            pz1 = pz(i1)
5589            e2 = e(i2)
5590            px2 = px(i2)
5591            py2 = py(i2)
5592            pz2 = pz(i2)
5593
5594            rts2 = (e1 + e2) ** 2 - (px1 + px2) ** 2
5595      &          - (py1 + py2) ** 2 - (pz1 + pz2) ** 2
5596
5597            allok = allok .and. rts2 .gt. rscut2
5598         end if
5599
5600         if (.not. allok) then
5601            tm = tlarge
5602            t1 = tlarge
5603            t2 = tlarge
5604         else
5605            t1 = tm
5606            t2 = tm
5607         end if
5608
5609         return
5610         end
5611
5612         subroutine reor(t, tmin, j, last0)
5613 c       this subroutine is used to fix ct(i) when tm is greater than ct(i)
5614 c       next(i) is last1 or last2
5615
5616         implicit double precision (a-h, o-z)
5617         parameter (MAXPTN=400001)
5618
5619         common /para5/ iconfg, iordsc
5620 cc      SAVE /para5/
5621         common /ilist1/
5622      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
5623      &     ictype, icsta(MAXPTN),
5624      &     nic(MAXPTN), icels(MAXPTN)
5625 cc      SAVE /ilist1/
5626 cd        common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
5627 cc      SAVE /ilist5/
5628         SAVE
5629
5630         icels0 = icels(j)
5631
5632         i1 = icels0 / 10000
5633         i2 = (icels0 - i1 * 10000) / 100
5634         i3 = icels0 - i1 * 10000 - i2 * 100
5635
5636         call wallc(j, i1, i2, i3, t, tmin1)
5637
5638         if (tmin .le. tmin1) then
5639            nc = last0
5640         else
5641            tmin = tmin1
5642            nc = 0
5643         end if
5644
5645         if (iconfg .eq. 3 .or. iconfg .eq. 5) then
5646            call chcell(j, i1, i2, i3, last0, t, tmin, nc)
5647         else
5648            if (i1 .eq. 11 .and. i2 .eq. 11 .and. i3 .eq. 11) then
5649               call chout(j, last0, t, tmin, nc)
5650            else
5651               if (iconfg .eq. 1) then
5652                  call chin1(j, i1, i2, i3, last0, t, tmin, nc)
5653               else if (iconfg .eq. 2) then
5654                  call chin2(j, i1, i2, i3, last0, t, tmin, nc)
5655               else if (iconfg .eq. 4) then
5656                  call chin3(j, i1, i2, i3, last0, t, tmin, nc)
5657               end if
5658            end if
5659         end if
5660
5661         call fixtim(j, t, tmin1, tmin, nc)
5662
5663         return
5664         end
5665
5666         subroutine chout(l, last0, t, tmin, nc)
5667 c       this subroutine is used to check the surface when the colliding
5668 c       particle is outside the cube
5669
5670         implicit double precision (a-h, o-z)
5671         parameter (MAXPTN=400001)
5672
5673         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
5674      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
5675      &       xmass(MAXPTN), ityp(MAXPTN)
5676 cc      SAVE /prec2/
5677         SAVE
5678
5679         m1 = 11
5680         m2 = 11
5681         m3 = 11
5682         call chcell(l, m1, m2, m3, last0, t, tmin, nc)
5683
5684         do 1003 i = 1, 10
5685            do 1002 j = 1, 10
5686               do 1001 k = 1, 10
5687                  if (i .eq. 1 .or. i .eq. 10 .or. j .eq. 1 .or.
5688      &              j .eq. 10 .or. k .eq. 1 .or. k. eq. 10)
5689      &               call chcell(l, i, j, k, last0, t, tmin, nc)
5690  1001         continue
5691  1002      continue
5692  1003   continue
5693
5694         return
5695         end
5696
5697         subroutine chin1(l, i1, i2, i3, last0, t, tmin, nc)
5698 c       this subroutine is used to check collisions for particle inside
5699 c       the cube
5700 c       and update the afftected particles through chcell
5701
5702         implicit double precision (a-h, o-z)
5703         SAVE
5704
5705 c       itest is a flag to make sure the 111111 cell is checked only once
5706         itest = 0
5707
5708         do 1003 i = i1 - 1, i1 + 1
5709            do 1002 j = i2 - 1, i2 + 1
5710               do 1001 k = i3 - 1, i3 + 1
5711                  if (i .ge. 1 .and. i .le. 10
5712      &              .and. j .ge. 1 .and. j .le. 10
5713      &              .and. k .ge. 1 .and. k .le. 10) then
5714                     call chcell(l, i, j, k, last0, t, tmin, nc)
5715                  else if (itest .eq. 0) then
5716                     m1 = 11
5717                     m2 = 11
5718                     m3 = 11
5719                     call chcell(l, m1, m2, m3, last0, t, tmin, nc)
5720                     itest = 1
5721                  end if
5722  1001         continue
5723  1002      continue
5724  1003   continue
5725
5726         return
5727         end
5728
5729         subroutine chin2(l, i1, i2, i3, last0, t, tmin, nc)
5730 c       this subroutine is used to check collisions for particle inside
5731 c       the cube
5732 c       and update the afftected particles through chcell
5733
5734         implicit double precision (a-h, o-z)
5735         SAVE
5736
5737 c       itest is a flag to make sure the 111111 cell is checked only once
5738         itest = 0
5739
5740         do 1003 i = i1 - 1, i1 + 1
5741            do 1002 j = i2 - 1, i2 + 1
5742               do 1001 k = i3 - 1, i3 + 1
5743                  ia = i
5744                  ib = j
5745                  ic = k
5746                  if (k .ge. 1 .and. k .le. 10) then
5747                     if (i .eq. 0) ia = 10
5748                     if (i .eq. 11) ia = 1
5749                     if (j .eq. 0) ib = 10
5750                     if (j .eq. 11) ib = 1
5751                     call chcell(l, ia, ib, ic, last0, t, tmin, nc)
5752                  end if
5753  1001         continue
5754  1002      continue
5755  1003   continue
5756
5757         return
5758         end
5759
5760         subroutine chin3(l, i1, i2, i3, last0, t, tmin, nc)
5761 c       this subroutine is used to check collisions for particle inside
5762 c       the cube
5763 c       and update the afftected particles through chcell
5764
5765         implicit double precision (a-h, o-z)
5766         SAVE
5767
5768 c       itest is a flag to make sure the 111111 cell is checked only once
5769         itest = 0
5770
5771         do 1003 i = i1 - 1, i1 + 1
5772            do 1002 j = i2 - 1, i2 + 1
5773               do 1001 k = i3 - 1, i3 + 1
5774                  if (i .eq. 0) then
5775                     ia = 10
5776                  else if (i .eq. 11) then
5777                     ia = 1
5778                  else
5779                     ia = i
5780                  end if
5781                  if (j .eq. 0) then
5782                     ib = 10
5783                  else if (j .eq. 11) then
5784                     ib = 1
5785                  else
5786                     ib = j
5787                  end if
5788                  if (k .eq. 0) then
5789                     ic = 10
5790                  else if (k .eq. 11) then
5791                     ic = 1
5792                  else
5793                     ic = k
5794                  end if
5795                  call chcell(l, ia, ib, ic, last0, t, tmin, nc)
5796  1001         continue
5797  1002      continue
5798  1003   continue
5799
5800         return
5801         end
5802
5803         subroutine chcell(il, i1, i2, i3, last0, t, tmin, nc)
5804 c       this program is used to check through all the particles, except last0
5805 c       in the cell (i1,i2,i3) and see if we can get a particle collision
5806 c       with time less than the original input tmin ( the collision time of
5807 c       il with the wall
5808 c       last0 cas be set to 0 if we don't want to exclude last0
5809
5810         implicit double precision (a-h, o-z)
5811         parameter (MAXPTN=400001)
5812
5813         common /para5/ iconfg, iordsc
5814 cc      SAVE /para5/
5815         common /ilist1/
5816      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
5817      &     ictype, icsta(MAXPTN),
5818      &     nic(MAXPTN), icels(MAXPTN)
5819 cc      SAVE /ilist1/
5820         common /ilist2/ icell, icel(10, 10, 10)
5821 cc      SAVE /ilist2/
5822         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
5823 cc      SAVE /ilist4/
5824         SAVE
5825
5826         t=t
5827         if (iconfg .eq. 3 .or. iconfg .eq. 5) then
5828            jj = ichkpt
5829            do 1001 j = 1, jj
5830 c     10/24/02 get rid of argument usage mismatch in mintm():
5831               jmintm=j
5832               if (j .ne. il .and. j .ne. last0)
5833      &          call mintm(il, jmintm, tmin, nc)
5834 c     &          call mintm(il, j, tmin, nc)
5835
5836  1001         continue
5837            return
5838         end if
5839
5840 c       set l
5841         if (i1 .eq. 11 .and. i2 .eq. 11 .and. i3 .eq. 11) then
5842            l = icell
5843         else
5844            l = icel(i1 ,i2, i3)
5845         end if
5846
5847         if (l .eq. 0) return
5848
5849         j = nic(l)
5850
5851 c       if there is only one particle
5852         if (j .eq. 0) then
5853
5854 c       if it's not il or last0,when last is not wall
5855            if (l .eq. il .or. l .eq. last0) return
5856            call mintm(il, l, tmin, nc)
5857
5858 c       if there are many particles
5859         else
5860            if (l .ne. il .and. l .ne. last0)
5861      &        call mintm(il, l, tmin, nc)
5862            do 10 while(j .ne. l)
5863               if (j .ne. il .and. j .ne. last0)
5864      &             call mintm(il, j, tmin, nc)
5865               j = nic(j)
5866  10           continue
5867         end if
5868
5869         return
5870         end
5871
5872         subroutine mintm(i, j, tmin, nc)
5873 c       this subroutine is used to check whether particle j has smaller
5874 c       collision time with particle i than other particles
5875 c       or in other words, update next(i)
5876
5877 c       input i,j,tmin,nc
5878 c       output tmin,nc
5879
5880         implicit double precision (a-h, o-z)
5881         parameter (MAXPTN=400001)
5882
5883         common /para5/ iconfg, iordsc
5884 cc      SAVE /para5/
5885         common /aurec1/ jxa, jya, jza
5886 cc      SAVE /aurec1/
5887         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
5888 cc      SAVE /aurec2/
5889         common /ilist1/
5890      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
5891      &     ictype, icsta(MAXPTN),
5892      &     nic(MAXPTN), icels(MAXPTN)
5893 cc      SAVE /ilist1/
5894         common /ilist3/ size1, size2, size3, v1, v2, v3, size
5895 cc      SAVE /ilist3/
5896         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
5897 cc      SAVE /ilist5/
5898         SAVE
5899
5900         logical allok
5901
5902         call isco(i, j, allok, tm, t1, t2)
5903
5904         if (allok .and. tm .lt. tmin) then
5905            tmin = tm
5906            ct(i) = t1
5907            nc = j
5908            if (iconfg .eq. 3 .or. iconfg .eq. 5) then
5909               dgxa(i) = - jxa * 10d0 * size1
5910               dgya(i) = - jya * 10d0 * size2
5911               if (iconfg .eq. 5) then
5912                  dgza(i) = - jza * 10d0 * size3
5913               end if
5914            end if
5915         end if
5916
5917          return
5918         end
5919
5920 ******************************************************************************
5921 ******************************************************************************
5922
5923         subroutine zpca1
5924
5925         implicit double precision (a-h, o-z)
5926         parameter (MAXPTN=400001)
5927
5928         common /ilist1/
5929      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
5930      &     ictype, icsta(MAXPTN),
5931      &     nic(MAXPTN), icels(MAXPTN)
5932 cc      SAVE /ilist1/
5933         SAVE
5934
5935         if (mod(ictype,2) .eq. 0) then
5936            call zpca1a(iscat)
5937            call zpca1a(jscat)
5938 clin-5/2009 ctest off v2 for parton:
5939 c           call flowp(1)
5940         end if
5941
5942         return
5943         end
5944
5945         subroutine zpca1a(i)
5946
5947         implicit double precision (a-h, o-z)
5948         parameter (MAXPTN=400001)
5949
5950         common /para2/ xmp, xmu, alpha, rscut2, cutof2
5951 cc      SAVE /para2/
5952         common /para5/ iconfg, iordsc
5953 cc      SAVE /para5/
5954         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
5955      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
5956      &       xmass(MAXPTN), ityp(MAXPTN)
5957 cc      SAVE /prec2/
5958         common /prec3/gxs(MAXPTN),gys(MAXPTN),gzs(MAXPTN),fts(MAXPTN),
5959      &     pxs(MAXPTN), pys(MAXPTN), pzs(MAXPTN), es(MAXPTN),
5960      &     xmasss(MAXPTN), ityps(MAXPTN)
5961 cc      SAVE /prec3/
5962         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
5963 cc      SAVE /prec5/
5964         common /prec6/ etas(MAXPTN), raps(MAXPTN), taus(MAXPTN)
5965 cc      SAVE /prec6/
5966         common /ana1/ ts(12)
5967 cc      SAVE /ana1/
5968         SAVE
5969
5970         if (iconfg .eq. 1) then
5971            t1 = fts(i)
5972            t2 = ft(i)
5973            ipic = 11
5974         else if (iconfg .eq. 2 .or.
5975      &     iconfg .eq. 3) then
5976 cd           t1 = fts(i)
5977 cd           t2 = ft(i)
5978            t1 = taus(i)
5979            t2 = tau(i)
5980            ipic = 12
5981         else if (iconfg .eq. 4 .or.
5982      &     iconfg .eq. 5) then
5983            t1 = fts(i)
5984            t2 = ft(i)
5985            ipic = 12
5986         end if
5987
5988         if (iconfg .le. 3) then
5989            do 1002 ian = 1, ipic
5990               if (t1 .le. ts(ian) .and.
5991      &           t2 .gt. ts(ian)) then
5992                  rapi = raps(i)
5993 c     7/20/01:
5994 c                 et = sqrt(pxs(i) ** 2 + pys(i) ** 2 + xmp ** 2)
5995                  et = dsqrt(pxs(i) ** 2 + pys(i) ** 2 + xmp ** 2)
5996                  call zpca1b(rapi, et, ian)
5997               end if
5998  1002      continue
5999         else
6000            do 1003 ian = 1, ipic
6001               if (t1 .le. ts(ian) .and.
6002      &           t2 .gt. ts(ian)) then
6003                  p0 = es(i)
6004                  p1 = pxs(i)
6005                  p2 = pys(i)
6006                  p3 = pzs(i)
6007                  call zpca1c(p0, p1, p2, p3, ian)
6008               end if
6009  1003      continue
6010         end if
6011
6012         return
6013         end
6014
6015         subroutine zpca1b(rapi, et, ian)
6016
6017         implicit double precision (a-h, o-z)
6018         parameter (MAXPTN=400001)
6019
6020         common /para6/ centy
6021 cc      SAVE /para6/
6022         common /ilist6/ t, iopern, icolln
6023 cc      SAVE /ilist6/
6024         common /ana2/
6025      &     det(12), dn(12), detdy(12), detdn(12), dndy(12),
6026      &     det1(12), dn1(12), detdy1(12), detdn1(12), dndy1(12),
6027      &     det2(12), dn2(12), detdy2(12), detdn2(12), dndy2(12)
6028 cc      SAVE /ana2/
6029         SAVE
6030
6031         if (rapi .gt. centy - 0.5d0 .and.
6032      &     rapi .lt. centy + 0.5d0) then
6033            det2(ian) = det2(ian) + et
6034            dn2(ian) = dn2(ian) + 1d0
6035 cdtrans
6036            if (ian .eq. 10) then
6037 cd              write (10, *) t, det2(ian)
6038            end if
6039            if (ian .eq. 11) then
6040 cd              write (11, *) t, det2(ian)
6041            end if
6042            if (ian .eq. 12) then
6043 cd              write (12, *) t, det2(ian)
6044            end if
6045 cdtransend
6046            if (rapi .gt. centy - 0.25d0 .and.
6047      &        rapi .lt. centy + 0.25d0) then
6048               det1(ian) = det1(ian) + et
6049               dn1(ian) = dn1(ian) + 1d0
6050               if (rapi .gt. centy - 0.1d0 .and.
6051      &           rapi .lt. centy + 0.1d0) then
6052                  det(ian) = det(ian) + et
6053                  dn(ian) = dn(ian) + 1d0
6054               end if
6055            end if
6056         end if
6057
6058         return
6059         end
6060
6061         subroutine zpca1c(p0, p1, p2, p3, ian)
6062
6063         implicit double precision (a-h, o-z)
6064
6065         common /ana3/ em(4, 4, 12)
6066 cc      SAVE /ana3/
6067
6068         dimension en(4)
6069         SAVE
6070
6071         en(1) = p0
6072         en(2) = p1
6073         en(3) = p2
6074         en(4) = p3
6075
6076         do 1002 i = 1, 4
6077            do 1001 j = 1, 4
6078               em(i, j, ian) = em(i, j, ian) + en(i) * en(j) / p0
6079  1001      continue
6080  1002   continue
6081
6082         return
6083         end
6084
6085 ******************************************************************************
6086
6087         subroutine zpca2
6088
6089         implicit double precision (a-h, o-z)
6090         parameter (MAXPTN=400001)
6091
6092         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
6093 cc      SAVE /para3/
6094         common /para5/ iconfg, iordsc
6095 cc      SAVE /para5/
6096         common /para7/ ioscar,nsmbbbar,nsmmeson
6097 cc      SAVE /para7/
6098         common /ilist6/ t, iopern, icolln
6099 cc      SAVE /ilist6/
6100         common /rndm1/ number
6101 cc      SAVE /rndm1/
6102         common /rndm2/ iff
6103 cc      SAVE /rndm2/
6104         common /rndm3/ iseedp
6105 cc      SAVE /rndm3/
6106         COMMON /AREVT/ IAEVT, IARUN, MISS
6107 cc      SAVE /AREVT/
6108         SAVE
6109
6110         if (iconfg .le. 3) then
6111            call zpca2a
6112         else
6113            call zpca2b
6114         end if
6115
6116         if (ioscar .eq. 1) then
6117            call zpca2c
6118         end if
6119
6120 cbzdbg2/17/99
6121 c        write (25, *) 'Event', nsevt - 1 + ievt,
6122 c    &         ', run', isbrun,
6123 c        WRITE (25, *) ' Event ', IAEVT, ', run ', IARUN,
6124 c     &     ',\n\t number of operations = ', iopern,
6125 c     &     ',\n\t number of collisions between particles = ',
6126 c     &         icolln,
6127 c     &     ',\n\t freezeout time=', t,
6128 c     &     ',\n\t ending at the ', number, 'th random number',
6129 c     &     ',\n\t ending collision iff=', iff
6130 cms     WRITE (25, *) ' Event ', IAEVT, ', run ', IARUN
6131 cms     WRITE (25, *) '    number of operations = ', iopern
6132 cms     WRITE (25, *) '    number of collisions between particles = ',
6133 cms  &       icolln
6134 cms     WRITE (25, *) '    freezeout time=', t
6135 cms     WRITE (25, *) '    ending at the ', number, 'th random number'
6136 cms     WRITE (25, *) '    ending collision iff=', iff
6137
6138         return
6139         end
6140
6141         subroutine zpca2a
6142
6143         implicit double precision (a-h, o-z)
6144         parameter (MAXPTN=400001)
6145
6146         common /para1/ mul
6147 cc      SAVE /para1/
6148         common /para2/ xmp, xmu, alpha, rscut2, cutof2
6149 cc      SAVE /para2/
6150         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
6151 cc      SAVE /para3/
6152         common /para5/ iconfg, iordsc
6153 cc      SAVE /para5/
6154         common /para6/ centy
6155 cc      SAVE /para6/
6156         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
6157      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
6158      &       xmass(MAXPTN), ityp(MAXPTN)
6159 cc      SAVE /prec2/
6160         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
6161 cc      SAVE /prec5/
6162         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
6163 cc      SAVE /ilist4/
6164         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
6165 cc      SAVE /ilist5/
6166         common /ilist6/ t, iopern, icolln
6167 cc      SAVE /ilist6/
6168         common /rndm1/ number
6169 cc      SAVE /rndm1/
6170         common /rndm2/ iff
6171 cc      SAVE /rndm2/
6172         common /rndm3/ iseedp
6173 cc      SAVE /rndm3/
6174         common /ana1/ ts(12)
6175 cc      SAVE /ana1/
6176         common /ana2/
6177      &     det(12), dn(12), detdy(12), detdn(12), dndy(12),
6178      &     det1(12), dn1(12), detdy1(12), detdn1(12), dndy1(12),
6179      &     det2(12), dn2(12), detdy2(12), detdn2(12), dndy2(12)
6180 cc      SAVE /ana2/
6181         common /ana4/ fdetdy(24), fdndy(24), fdndpt(12)
6182 cc      SAVE /ana4/
6183         SAVE
6184
6185         do 1004 i = 1, ichkpt
6186            rapi = rap(i)
6187 c     7/20/01:
6188 c           et = sqrt(px(i) ** 2 + py(i) ** 2 + xmp ** 2)
6189            et = dsqrt(px(i) ** 2 + py(i) ** 2 + xmp ** 2)
6190
6191            do 1001 j = 1, 24
6192               if (rapi .gt. j + centy - 13d0
6193      &           .and. rapi .lt. j  + centy - 12d0) then
6194                  fdetdy(j) = fdetdy(j) + et
6195                  fdndy(j) = fdndy(j) + 1d0
6196               end if
6197  1001      continue
6198
6199            do 1002 j = 1, 12
6200               if (et .gt. 0.5d0 * (j - 1) .and.
6201      &           et .lt. 0.5d0 * j ) then
6202                  fdndpt(j) = fdndpt(j) + 1d0
6203               end if
6204  1002      continue
6205
6206            if (iconfg .eq. 1) then
6207               t1 = ft(i)
6208               t2 = tlarge
6209               ipic = 11
6210            else
6211               t1 = tau(i)
6212               t2 = tlarge
6213               ipic = 12
6214            end if
6215
6216            do 1003 ian = 1, ipic
6217               if (t1 .le. ts(ian) .and.
6218      &           t2 .gt. ts(ian)) then
6219                  call zpca1b(rapi, et, ian)
6220               end if
6221  1003      continue
6222
6223            if (iconfg .eq. 1) then
6224               call zpca1b(rapi, et, 12)
6225            end if
6226  1004   continue
6227
6228         do 1005 ian = 1, 12
6229 c           if (dn(ian) .eq. 0d0 .or. dn1(ian) .eq. 0d0 .or.
6230 c     &        dn2(ian) .eq. 0d0) then
6231 c              print *, 'event=', ievt
6232 c              print *, 'dn(', ian, ')=', dn(ian), 'dn1(', ian,
6233 c     &           ')=', dn1(ian), 'dn2(', ian, ')=', dn2(ian)
6234 c           end if
6235            detdy(ian) = detdy(ian) + det(ian)
6236            if (dn(ian) .ne. 0) then
6237               detdn(ian) = detdn(ian) + det(ian) / dn(ian)
6238            end if
6239            dndy(ian) = dndy(ian) + dn(ian)
6240            detdy1(ian) = detdy1(ian) + det1(ian)
6241            if (dn1(ian) .ne. 0) then
6242               detdn1(ian) = detdn1(ian) + det1(ian) / dn1(ian)
6243            end if
6244            dndy1(ian) = dndy1(ian) + dn1(ian)
6245            detdy2(ian) = detdy2(ian) + det2(ian)
6246            if (dn2(ian) .ne. 0) then
6247               detdn2(ian) = detdn2(ian) + det2(ian) / dn2(ian)
6248            end if
6249            dndy2(ian) = dndy2(ian) + dn2(ian)
6250  1005   continue
6251
6252         return
6253         end
6254
6255         subroutine zpca2b
6256
6257         implicit double precision (a-h, o-z)
6258         parameter (MAXPTN=400001)
6259
6260         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
6261      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
6262      &       xmass(MAXPTN), ityp(MAXPTN)
6263 cc      SAVE /prec2/
6264         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
6265 cc      SAVE /ilist4/
6266         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
6267 cc      SAVE /ilist5/
6268         common /ana1/ ts(12)
6269 cc      SAVE /ana1/
6270         SAVE
6271
6272         do 1002 i = 1, ichkpt
6273            t1 = ft(i)
6274            t2 = tlarge
6275            ipic = 12
6276
6277            do 1001 ian = 1, ipic
6278               if (t1 .le. ts(ian) .and.
6279      &           t2 .gt. ts(ian)) then
6280                  p0 = e(i)
6281                  p1 = px(i)
6282                  p2 = py(i)
6283                  p3 = pz(i)
6284                  call zpca1c(p0, p1, p2, p3, ian)
6285               end if
6286  1001      continue
6287  1002   continue
6288
6289         return
6290         end
6291
6292         subroutine zpca2c
6293
6294         implicit double precision (a-h, o-z)
6295
6296         character*8 code, versn
6297         character*4 reffra
6298         integer aproj, zproj, atarg, ztarg, event
6299         parameter (MAXPTN=400001)
6300
6301         common /para1/ mul
6302 cc      SAVE /para1/
6303         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
6304      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
6305      &       xmass(MAXPTN), ityp(MAXPTN)
6306 cc      SAVE /prec2/
6307         SAVE
6308         data nff/0/
6309
6310 c       file header
6311         if (nff .eq. 0) then
6312 cms        write (26, 101) 'OSCAR1997A'
6313 cms        write (26, 101) 'final_id_p_x'
6314            code = 'ZPC'
6315            versn = '1.0.1'
6316            aproj = -1
6317            zproj = -1
6318            atarg = -1
6319            ztarg = -1
6320            reffra = 'cm'
6321            ebeam = 0d0
6322            ntestp = 1
6323 cms        write (26, 102) code, versn, aproj, zproj, atarg, ztarg,
6324 cms  &        reffra, ebeam, ntestp
6325            nff = 1
6326            event = 1
6327            bimp = 0d0
6328            phi = 0d0
6329         end if
6330
6331 c       comment
6332
6333 c       event header
6334 cms     write (26, 103) event, mul, bimp, phi
6335
6336 c       particles
6337         do 99 i = 1, mul
6338 cms        write (26, 104) i, ityp(i),
6339 cms  &        px(i), py(i), pz(i), e(i), xmass(i),
6340 cms  &        gx(i), gy(i), gz(i), ft(i)
6341  99         continue
6342
6343          event = event + 1
6344
6345 cyy 101        format (a12)
6346 cyy 102        format (2(a8, 2x), '(',i3, ',',i6, ')+(',i3, ',', i6, ')',
6347 cyy     &     2x, a4, 2x, e10.4, 2x, i8)
6348 cyy 103        format (i10, 2x, i10, 2x, f8.3, 2x, f8.3)
6349 cyy 104        format (i10, 2x, i10, 2x, 9(e12.6, 2x))
6350
6351         return
6352         end
6353
6354 ******************************************************************************
6355
6356         subroutine zpcou
6357
6358         implicit double precision (a-h, o-z)
6359
6360         common /para5/ iconfg, iordsc
6361 cc      SAVE /para5/
6362         SAVE
6363
6364         if (iconfg .le. 3) then
6365            call zpcou1
6366         else
6367            call zpcou2
6368         end if
6369
6370         return
6371         end
6372
6373         subroutine zpcou1
6374
6375         implicit double precision (a-h, o-z)
6376
6377         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
6378 cc      SAVE /para3/
6379         common /ana1/ ts(12)
6380 cc      SAVE /ana1/
6381         common /ana2/
6382      &     det(12), dn(12), detdy(12), detdn(12), dndy(12),
6383      &     det1(12), dn1(12), detdy1(12), detdn1(12), dndy1(12),
6384      &     det2(12), dn2(12), detdy2(12), detdn2(12), dndy2(12)
6385 cc      SAVE /ana2/
6386         common /ana4/ fdetdy(24), fdndy(24), fdndpt(12)
6387 cc      SAVE /ana4/
6388         SAVE
6389 c
6390         dpt = 0.5d0
6391         dy2 = 1d0
6392         dy1 = 0.5d0
6393         dy = 0.2d0
6394         ntotal = nevnt * nsbrun
6395 c
6396         return
6397         end
6398
6399         subroutine zpcou2
6400
6401         implicit double precision (a-h, o-z)
6402
6403         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
6404 cc      SAVE /para3/
6405         common /ilist3/ size1, size2, size3, v1, v2, v3, size
6406 cc      SAVE /ilist3/
6407         common /ana1/ ts(12)
6408 cc      SAVE /ana1/
6409         common /ana3/ em(4, 4, 12)
6410 cc      SAVE /ana3/
6411         SAVE
6412 c
6413 cms     open (28, file = 'ana4/em.dat', status = 'unknown')
6414         vol = 1000.d0 * size1 * size2 * size3
6415         ntotal = nevnt * nsbrun
6416
6417         do 1002 ian = 1, 12
6418 cms        write (28, *) '*** for time ', ts(ian), 'fm(s)'
6419            do 1001 i = 1, 4
6420 cms           write (28, *) em(i, 1, ian) / vol / ntotal,
6421 cms  &                        em(i, 2, ian) / vol / ntotal,
6422 cms  &                        em(i, 3, ian) / vol / ntotal,
6423 cms  &                        em(i, 4, ian) / vol / ntotal
6424  1001      continue
6425  1002   continue
6426
6427         return
6428         end
6429
6430 ******************************************************************************
6431
6432       subroutine lorenza(energy, px, py, pz, bex, bey, bez)
6433
6434 c     add in a cut for beta2 to prevent gam to be nan (infinity)
6435
6436       implicit double precision (a-h, o-z)
6437
6438       common /lora/ enenew, pxnew, pynew, pznew
6439 cc      SAVE /lora/
6440       SAVE
6441
6442       beta2 = bex ** 2 + bey ** 2 + bez ** 2
6443       if (beta2 .eq. 0d0) then
6444          enenew = energy
6445          pxnew = px
6446          pynew = py
6447          pznew = pz
6448       else
6449          if (beta2 .gt. 0.999999999999999d0) then
6450             beta2 = 0.999999999999999d0
6451             print *,'beta2=0.999999999999999'
6452          end if
6453 clin-7/20/01:
6454 c         gam = 1.d0 / sqrt(1.d0 - beta2)
6455          gam = 1.d0 / dsqrt(1.d0 - beta2)
6456          enenew = gam * (energy - bex * px - bey * py - bez * pz)
6457          pxnew = - gam * bex * energy + (1.d0
6458      &        + (gam - 1.d0) * bex ** 2 / beta2) * px
6459      &        + (gam - 1.d0) * bex * bey/beta2 * py
6460      &        + (gam - 1.d0) * bex * bez/beta2 * pz
6461          pynew = - gam * bey * energy
6462      &        + (gam - 1.d0) * bex * bey / beta2 * px
6463      &        + (1.d0 + (gam - 1.d0) * bey ** 2 / beta2) * py
6464      &        + (gam - 1.d0) * bey * bez / beta2 * pz
6465          pznew = - gam * bez * energy
6466      &        +  (gam - 1.d0) * bex * bez / beta2 * px
6467      &        + (gam - 1.d0) * bey * bez / beta2 * py
6468      &        + (1.d0 + (gam - 1.d0) * bez ** 2 / beta2) * pz
6469       endif
6470
6471       return
6472       end
6473
6474       subroutine index1(n, m, arrin, indx)
6475 c     indexes the first m elements of ARRIN of length n, i.e., outputs INDX
6476 c     such that ARRIN(INDEX(J)) is in ascending order for J=1,...,m
6477
6478       implicit double precision (a-h, o-z)
6479
6480       dimension arrin(n), indx(n)
6481       SAVE
6482       do 1001 j = 1, m
6483          indx(j) = j
6484  1001   continue
6485       l = m / 2 + 1
6486       ir = m
6487  10   continue
6488       if (l .gt. 1) then
6489          l = l - 1
6490          indxt = indx(l)
6491          q = arrin(indxt)
6492       else
6493          indxt = indx(ir)
6494          q = arrin(indxt)
6495          indx(ir) = indx(1)
6496          ir = ir - 1
6497          if (ir .eq. 1) then
6498             indx(1) = indxt
6499             return
6500          end if
6501       end if
6502       i = l
6503       j = l + l
6504  20   if (j .le. ir) then
6505          if (j .lt. ir) then
6506             if (arrin(indx(j)) .lt. arrin(indx(j + 1))) j = j + 1
6507          end if
6508          if (q .lt. arrin(indx(j))) then
6509             indx(i) = indx(j)
6510             i = j
6511             j = j + j
6512          else
6513             j = ir + 1
6514          end if
6515       goto 20
6516       end if
6517       indx(i) = indxt
6518       goto 10
6519
6520       end
6521
6522
6523         double precision function ftime1(iseed)
6524
6525 c       this program is used to generate formation time
6526 c       the calling program needs a common /par1/
6527 c       and declare external ftime1
6528
6529 clin-8/19/02
6530         implicit double precision (a-h, o-z)
6531
6532 cc        external ran1
6533
6534         parameter (hbarc = 0.197327054d0)
6535
6536         common /par1/ formt
6537 cc      SAVE /par1/
6538         SAVE
6539
6540         aa = hbarc / formt
6541
6542 clin7/20/01:
6543 c        ftime1 = aa * sqrt(1d0 / ran1(iseed) - 1d0)
6544         ftime1 = aa * dsqrt(1d0 / ran1(iseed) - 1d0)
6545         return
6546         end
6547
6548
6549       subroutine cropro(vx1, vy1, vz1, vx2, vy2, vz2)
6550
6551 c     this subroutine is used to calculate the cross product of
6552 c     (vx1,vy1,vz1) and (vx2,vy2,vz2) and get the result (vx3,vy3,vz3)
6553 c     and put the vector into common /cprod/
6554
6555       implicit double precision (a-h, o-z)
6556
6557       common/cprod/ vx3, vy3, vz3
6558 cc      SAVE /cprod/
6559       SAVE
6560
6561       vx3 = vy1 * vz2 - vz1 * vy2
6562       vy3 = vz1 * vx2 - vx1 * vz2
6563       vz3 = vx1 * vy2 - vy1 * vx2
6564
6565       return
6566       end
6567
6568       subroutine xnormv(vx, vy, vz)
6569
6570 c      this subroutine is used to get a normalized vector
6571
6572       implicit double precision (a-h, o-z)
6573       SAVE
6574
6575 clin-7/20/01:
6576 c      vv = sqrt(vx ** 2 + vy ** 2 + vz ** 2)
6577       vv = dsqrt(vx ** 2 + vy ** 2 + vz ** 2)
6578       vx = vx / vv
6579       vy = vy / vv
6580       vz = vz / vv
6581
6582       return
6583       end
6584
6585 cbz1/29/99
6586 c      subroutine rotate(xn1, xn2, xn3, theta, v1, v2, v3)
6587       subroutine zprota(xn1, xn2, xn3, theta, v1, v2, v3)
6588 cbz1/29/99end
6589
6590 c     this subroutine is used to rotate the vector (v1,v2,v3) by an angle theta
6591 c     around the unit vector (xn1, xn2, xn3)
6592
6593       implicit double precision (a-h, o-z)
6594       SAVE
6595
6596       vx = v1
6597       vy = v2
6598       vz = v3
6599       c = cos(theta)
6600       omc = 1d0 - c
6601       s = sin(theta)
6602       a11 = xn1 ** 2 * omc + c
6603       a12 = xn1 * xn2 * omc - s * xn3
6604       a13 = xn1 * xn3 * omc + s * xn2
6605       a21 = xn1 * xn2 * omc + s * xn3
6606       a22 = xn2 **2 * omc + c
6607       a23 = xn2 * xn3 * omc - s * xn1
6608       a31 = xn1 * xn3 * omc - s * xn2
6609       a32 = xn3 * xn2 * omc + s * xn1
6610       a33 = xn3 ** 2 * omc + c
6611       v1 = vx * a11 + vy * a12 + vz * a13
6612       v2 = vx * a21 + vy * a22 + vz * a23
6613       v3 = vx * a31 + vy * a32 + vz * a33
6614
6615       return
6616       end
6617
6618 c      double precision function ran1(idum)
6619
6620 c*     return a uniform random deviate between 0.0 and 1.0. set idum to
6621 c*     any negative value to initialize or reinitialize the sequence.
6622
6623 c      implicit double precision (a-h, o-z)
6624
6625 c      dimension r(97)
6626
6627 c      common /rndm1/ number
6628 cc      SAVE /rndm1/
6629 c      parameter (m1 = 259200, ia1 = 7141, ic1 = 54773, rm1 = 1d0 / m1)
6630 c      parameter (m2 = 134456, ia2 = 8121, ic2 = 28411, rm2 = 1d0 / m2)
6631 c      parameter (m3 = 243000, ia3 = 4561, ic3 = 51349)
6632 clin-6/23/00 save ix1-3:
6633 clin-10/30/02 r unsaved, causing wrong values for ran1 when compiled with f77:
6634 cc      SAVE ix1,ix2,ix3,r
6635 c      SAVE
6636 c      data iff/0/
6637
6638 c      if (idum .lt. 0 .or. iff .eq. 0) then
6639 c         iff = 1
6640 c         ix1 = mod(ic1 - idum, m1)
6641 c         ix1 = mod(ia1 * ix1 + ic1, m1)
6642 c         ix2 = mod(ix1, m2)
6643 c         ix1 = mod(ia1 * ix1 + ic1, m1)
6644 c         ix3 = mod(ix1, m3)
6645 c         do 11 j = 1, 97
6646 c            ix1 = mod(ia1 * ix1 + ic1, m1)
6647 c            ix2 = mod(ia2 * ix2 + ic2, m2)
6648 c            r(j) = (dble(ix1) + dble(ix2) * rm2) * rm1
6649 c 11         continue
6650 c         idum = 1
6651 c      end if
6652 c      ix1 = mod(ia1 * ix1 + ic1, m1)
6653 c      ix2 = mod(ia2 * ix2 + ic2, m2)
6654 c      ix3 = mod(ia3 * ix3 + ic3, m3)
6655 clin-7/01/02       j = 1 + (97 * i x 3) / m3
6656 c      j=1+(97*ix3)/m3
6657 clin-4/2008:
6658 c      if (j .gt. 97 .or. j .lt. 1) pause
6659 c      if (j .gt. 97 .or. j .lt. 1) print *, 'In zpc ran1, j<1 or j>97',j
6660 c      ran1 = r(j)
6661 c      r(j) = (dble(ix1) + dble(ix2) * rm2) * rm1
6662
6663 clin-6/23/00 check random number generator:
6664 c      number = number + 1
6665 c      if(number.le.100000) write(99,*) 'number, ran1=', number,ran1
6666
6667 c      return
6668 c      end