Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:13:23 

0001 c.................... zpc.f
0002 c       PROGRAM ZPC
0003       SUBROUTINE ZPCMN
0004 c       Version: 1.0.1
0005 c       Author: Bin Zhang 
0006 c       (suggestions, problems -> bzhang@nt1.phys.columbia.edu)
0007         implicit double precision (a-h, o-z)
0008 clin-4/20/01        PARAMETER (NMAXGL = 16000)
0009         parameter (MAXPTN=400001)
0010         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
0011 cc      SAVE /para3/
0012         SAVE   
0013 c
0014 c       loop over events
0015         do 1000 i = 1, nevnt
0016            ievt = i
0017 c       generation of the initial condition for one event
0018            call inievt
0019 c      loop over many runs of the same event
0020            do 2000 j = 1, nsbrun
0021               isbrun = j
0022 c       initialization for one run of an event
0023               call inirun
0024 clin-4/2008 not used:
0025 c             CALL HJAN1A
0026  3000         continue
0027 c       do one collision
0028               call zpcrun(*4000)
0029               call zpca1
0030               goto 3000
0031  4000         continue
0032               call zpca2
0033  2000      continue
0034  1000   continue
0035         call zpcou
0036 clin-5/2009 ctest off
0037 c     5/17/01 calculate v2 for parton already frozen out:
0038 c        call flowp(3)
0039 c.....to get average values for different strings
0040         CALL zpstrg
0041         RETURN
0042         end
0043 
0044 ******************************************************************************
0045 ******************************************************************************
0046 
0047         block data zpcbdt
0048 c       set initial values in block data
0049 
0050         implicit double precision (a-h, o-z)
0051         parameter (MAXPTN=400001)
0052         PARAMETER (MAXSTR=150001)
0053         common /para1/ mul
0054 cc      SAVE /para1/
0055         common /para2/ xmp, xmu, alpha, rscut2, cutof2
0056 cc      SAVE /para2/
0057         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
0058 cc      SAVE /para3/
0059         common /para4/ iftflg, ireflg, igeflg, ibstfg
0060 cc      SAVE /para4/
0061         common /para5/ iconfg, iordsc
0062 cc      SAVE /para5/
0063         common /para6/ centy
0064 cc      SAVE /para6/
0065 clin-6/2009 nsmbbbar and nsmmeson respectively give the total number of 
0066 c     baryons/anti-baryons and mesons for each event:
0067 c        common /para7/ ioscar
0068         common /para7/ ioscar,nsmbbbar,nsmmeson
0069 cc      SAVE /para7/
0070         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
0071      &       PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
0072      &       XMASS0(MAXPTN), ITYP0(MAXPTN)
0073 cc      SAVE /prec1/
0074         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
0075      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
0076      &       xmass(MAXPTN), ityp(MAXPTN)
0077 cc      SAVE /prec2/
0078         common /prec3/gxs(MAXPTN),gys(MAXPTN),gzs(MAXPTN),fts(MAXPTN),
0079      &     pxs(MAXPTN), pys(MAXPTN), pzs(MAXPTN), es(MAXPTN),
0080      &     xmasss(MAXPTN), ityps(MAXPTN)
0081 cc      SAVE /prec3/
0082         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
0083 cc      SAVE /prec4/
0084         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
0085 cc      SAVE /prec5/
0086         common /prec6/ etas(MAXPTN), raps(MAXPTN), taus(MAXPTN)
0087 cc      SAVE /prec6/
0088         common /aurec1/ jxa, jya, jza
0089 cc      SAVE /aurec1/
0090         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
0091 cc      SAVE /aurec2/
0092         common /ilist1/
0093      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
0094      &     ictype, icsta(MAXPTN),
0095      &     nic(MAXPTN), icels(MAXPTN)
0096 cc      SAVE /ilist1/
0097         common /ilist2/ icell, icel(10,10,10)
0098 cc      SAVE /ilist2/
0099         common /ilist3/ size1, size2, size3, v1, v2, v3, size
0100 cc      SAVE /ilist3/
0101         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
0102 cc      SAVE /ilist4/
0103 c     6/07/02 initialize in ftime to expedite compiling:
0104 c        common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
0105 cc      SAVE /ilist5/
0106         common /ilist6/ t, iopern, icolln
0107 cc      SAVE /ilist6/
0108         COMMON /ilist7/ LSTRG0(MAXPTN), LPART0(MAXPTN)
0109 cc      SAVE /ilist7/
0110         COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
0111 cc      SAVE /ilist8/
0112         common /rndm1/ number
0113 cc      SAVE /rndm1/
0114         common /rndm2/ iff
0115 cc      SAVE /rndm2/
0116         common /rndm3/ iseedp
0117 cc      SAVE /rndm3/
0118         common /ana1/ ts(12)
0119 cc      SAVE /ana1/
0120         common /ana2/
0121      &     det(12), dn(12), detdy(12), detdn(12), dndy(12),
0122      &     det1(12), dn1(12), detdy1(12), detdn1(12), dndy1(12),
0123      &     det2(12), dn2(12), detdy2(12), detdn2(12), dndy2(12)
0124 cc      SAVE /ana2/
0125         common /ana3/ em(4, 4, 12)
0126 cc      SAVE /ana3/
0127         common /ana4/ fdetdy(24), fdndy(24), fdndpt(12)
0128 cc      SAVE /ana4/
0129         SAVE   
0130         data centy/0d0/
0131 c     6/07/02 initialize in ftime to expedite compiling:
0132 c        data (ct(i), i = 1, MAXPTN)/MAXPTN*0d0/
0133 c        data (ot(i), i = 1, MAXPTN)/MAXPTN*0d0/
0134 c        data tlarge/1000000.d0/
0135         data number/0/
0136         data ts/0.11d0, 0.12d0, 0.15d0, 0.2d0, 0.3d0, 0.4d0, 0.6d0,
0137      &     0.8d0, 1d0, 2d0, 4d0, 6d0/
0138 c
0139         end
0140 
0141 ******************************************************************************
0142 ******************************************************************************
0143 
0144         subroutine inizpc
0145 
0146         implicit double precision (a-h, o-z)
0147         SAVE   
0148 
0149         call readpa
0150 
0151         call inipar
0152 
0153         call inian1
0154 
0155         return
0156         end
0157 
0158         subroutine readpa
0159 
0160         implicit double precision (a-h, o-z)
0161 
0162 cc        external ran1
0163 
0164 cms        character*50 str
0165 
0166         common /para2/ xmp, xmu, alpha, rscut2, cutof2
0167 cc      SAVE /para2/
0168         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
0169 cc      SAVE /para3/
0170         common /para4/ iftflg, ireflg, igeflg, ibstfg
0171 cc      SAVE /para4/
0172         common /para5/ iconfg, iordsc
0173 cc      SAVE /para5/
0174         common /para7/ ioscar,nsmbbbar,nsmmeson
0175 cc      SAVE /para7/
0176         common /ilist3/ size1, size2, size3, v1, v2, v3, size
0177 cc      SAVE /ilist3/
0178         common /rndm1/ number
0179 cc      SAVE /rndm1/
0180         common /rndm2/ iff
0181 cc      SAVE /rndm2/
0182         common /rndm3/ iseedp
0183 cc      SAVE /rndm3/
0184         SAVE   
0185 
0186         iseed=iseedp
0187 c       this is the initialization file containing the initial values of 
0188 c          the parameters
0189 cbz1/31/99
0190 c        open (5, file = 'zpc.ini', status = 'unknown')
0191 cbz1/31/99end
0192 
0193 c       this is the final data file containing general info about the cascade
0194 cbz1/31/99
0195 c        open (6, file = 'zpc.res', status = 'unknown')
0196 cms        open (25, file = 'ana/zpc.res', status = 'unknown')
0197 cbz1/31/99end
0198 
0199 c       this is the input file containing initial particle records
0200 cbz1/25/99
0201 c        open (7, file = 'zpc.inp', status = 'unknown')
0202 cbz1/25/99end
0203 
0204 c       this gives the optional OSCAR standard output
0205 cbz1/31/99
0206 c        open (8, file = 'zpc.oscar', status = 'unknown')
0207         if(ioscar.eq.1) then
0208 cms           open (26, file = 'ana/parton.oscar', status = 'unknown')
0209 cms           open (19, file = 'ana/hadron.oscar', status = 'unknown')
0210         endif
0211 cbz1/31/99end
0212 
0213 c     2/11/03 combine zpc initialization into ampt.ini:
0214 c        open (29, file = 'zpc.ini', status = 'unknown')
0215 c        read (29, *) str, xmp
0216         xmp=0d0
0217 c        read (29, *) str, xmu
0218 c        read (29, *) str, alpha
0219         cutof2 = 4.5d0 * (alpha / xmu) ** 2
0220 c        read (29, *) str, rscut2
0221         rscut2=0.01d0
0222 c        read (29, *) str, nsevt
0223         nsevt=1
0224 c        read (29, *) str, nevnt
0225         nevnt=1
0226 c        read (29, *) str, nsbrun
0227         nsbrun=1
0228 c        read (29, *) str, iftflg
0229         iftflg=0
0230 c        read (29, *) str, ireflg
0231         ireflg=1
0232 cbz1/31/99
0233         IF (ireflg .EQ. 0) THEN
0234 cms           OPEN (27, FILE = 'zpc.inp', STATUS = 'UNKNOWN')
0235         END IF
0236 cbz1/31/99end
0237 c        read (29, *) str, igeflg
0238         igeflg=0
0239 c        read (29, *) str, ibstfg
0240         ibstfg=0
0241 c        read (29, *) str, iconfg
0242         iconfg=1
0243 c        read (29, *) str, iordsc
0244         iordsc=11
0245 c        read (29, *) str, ioscar
0246 c        read (29, *) str, v1, v2, v3
0247         v1=0.2d0
0248         v2=0.2d0
0249         v3=0.2d0
0250 c        read (29, *) str, size1, size2, size3
0251         size1=1.5d0
0252         size2=1.5d0
0253         size3=0.7d0
0254         if (size1 .eq. 0d0 .or. size2 .eq. 0d0 .or. 
0255      &     size3 .eq. 0d0) then
0256            if (size1 .ne. 0d0 .or. size2 .ne. 0d0 .or. size3 .ne. 0d0
0257      &        .or. v1 .ne. 0d0 .or. v2 .ne. 0d0 .or. v3 .ne. 0d0) then
0258               print *, 'to get rid of space division:'
0259               print *, 'set all sizes and vs to 0'
0260               stop 'chker'
0261            end if
0262         end if
0263         size = min(size1, size2, size3)
0264 c        read (29, *) str, iff
0265         iff=-1
0266 c        read (29, *) str, iseed
0267 
0268 c     10/24/02 get rid of argument usage mismatch in ran1():
0269         isedng=-iseed
0270 c        a = ran1(-iseed)
0271         a = ran1(isedng)
0272 c        read (29, *) str, irused
0273         irused=2
0274         do 1001 i = 1, irused - 1
0275 c           a = ran1(2)
0276            iseed2=2
0277            a = ran1(iseed2)
0278  1001   continue
0279 c     10/24/02-end
0280 
0281         if (iconfg .eq. 2 .or. iconfg .eq. 3) then
0282            v1 = 0d0
0283            v2 = 0d0
0284         end if
0285 
0286         if (iconfg .eq. 4 .or. iconfg .eq. 5) then
0287            v1 = 0d0
0288            v2 = 0d0
0289            v3 = 0d0
0290         end if
0291 
0292         close(5)
0293 
0294         return
0295         end
0296 
0297         subroutine inipar
0298 
0299         implicit double precision (a-h,o-z)
0300 
0301         common /para4/ iftflg, ireflg, igeflg, ibstfg
0302 cc      SAVE /para4/
0303         common /para6/ centy
0304 cc      SAVE /para6/
0305         SAVE   
0306 
0307         if (ibstfg .ne. 0) then
0308            centy = -6d0
0309         end if
0310 
0311         return
0312         end
0313 
0314         subroutine inian1
0315 
0316         implicit double precision (a-h,o-z)
0317 
0318         common /para4/ iftflg, ireflg, igeflg, ibstfg
0319 cc      SAVE /para4/
0320         common /ana1/ ts(12)
0321 cc      SAVE /ana1/   
0322         SAVE   
0323         if (ibstfg .ne. 0) then
0324            a = cosh(6d0)
0325            do 1001 i = 1, 12
0326               ts(i) = ts(i) * a
0327  1001      continue
0328         end if
0329 
0330         return
0331         end
0332 
0333 ******************************************************************************
0334 
0335         subroutine inievt
0336 
0337         implicit double precision (a-h, o-z)
0338 
0339         COMMON /para1/ mul
0340 cc      SAVE /para1/
0341         common /para4/ iftflg, ireflg, igeflg, ibstfg
0342 cc      SAVE /para4/
0343         SAVE   
0344 
0345 cbz1/25/99
0346 c        mul = 0
0347 cbz1/25/99
0348         if (ireflg .eq. 0) call readi
0349         if (igeflg .ne. 0) call genei
0350         if (ibstfg .ne. 0) call boosti
0351 
0352         return
0353         end
0354 
0355         subroutine readi
0356 
0357         implicit double precision (a-h, o-z)
0358         parameter (MAXPTN=400001)
0359         double precision field(9)
0360         common /para1/ mul
0361 cc      SAVE /para1/
0362         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
0363 cc      SAVE /para3/
0364         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
0365      &       PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
0366      &       XMASS0(MAXPTN), ITYP0(MAXPTN)
0367 cc      SAVE /prec1/
0368         SAVE   
0369         do 1001 i = 1, MAXPTN
0370            if (ievt .ne. 1 .and. i .eq. 1) then
0371               ityp0(i) = ntyp
0372               gx0(1) = field(1)
0373               gy0(1) = field(2)
0374               gz0(1) = field(3)
0375               ft0(1) = field(4)
0376               px0(1) = field(5)
0377               py0(1) = field(6)
0378               pz0(1) = field(7)
0379               e0(1) = field(8)
0380               xmass0(i) = field(9)
0381               mul = 1
0382            else
0383  900              read (27, *, end = 1000) neve, ntyp, field
0384               if (neve .lt. nsevt) goto 900
0385               if (neve .gt.
0386      &           nsevt + ievt - 1) goto 1000
0387               ityp0(i) = ntyp
0388               gx0(i) = field(1)
0389               gy0(i) = field(2)
0390               gz0(i) = field(3)
0391               ft0(i) = field(4)
0392               px0(i) = field(5)
0393               py0(i) = field(6)
0394               pz0(i) = field(7)
0395               e0(i) = field(8)
0396               xmass0(i) = field(9)
0397               mul = mul + 1
0398            end if
0399  1001   continue
0400         
0401  1000        continue
0402         
0403         return
0404         end
0405 
0406         subroutine genei
0407 
0408         implicit double precision (a-h, o-z)
0409         parameter (MAXPTN=400001)
0410         common /para1/ mul
0411 cc      SAVE /para1/
0412         common /para2/ xmp, xmu, alpha, rscut2, cutof2
0413 cc      SAVE /para2/
0414         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
0415 cc      SAVE /para3/
0416         common /para5/ iconfg, iordsc
0417 cc      SAVE /para5/
0418         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
0419      &       PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
0420      &       XMASS0(MAXPTN), ITYP0(MAXPTN)
0421 cc      SAVE /prec1/
0422         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
0423 cc      SAVE /prec5/
0424         common /lor/ enenew, pxnew, pynew, pznew
0425 cc      SAVE /lor/
0426         common /rndm3/ iseedp
0427 cc      SAVE /rndm3/
0428         SAVE   
0429 
0430 cc        external ran1
0431 
0432         iseed=iseedp
0433         incmul = 4000
0434         temp = 0.5d0
0435         etamin = -5d0        
0436         etamax = 5d0
0437         r0 = 5d0
0438         tau0 = 0.1d0
0439         deta = etamax - etamin
0440 
0441         do 1001 i = mul + 1, mul + incmul
0442            ityp0(i) = 21
0443            xmass0(i) = xmp
0444            call energy(e, temp)
0445            call momntm(px, py, pz, e)
0446 c     7/20/01:
0447 c           e = sqrt(e ** 2 + xmp ** 2)
0448            e = dsqrt(e ** 2 + xmp ** 2)
0449            if (iconfg .le. 3) then
0450               eta(i) = etamin + deta * ran1(iseed)
0451               bex = 0d0
0452               bey = 0d0
0453               bez = -tanh(eta(i))
0454               call lorenz(e, px, py, pz, bex, bey, bez)
0455               px0(i) = pxnew
0456               py0(i) = pynew
0457               pz0(i) = pznew
0458               e0(i) = enenew
0459            else
0460               px0(i) = px
0461               py0(i) = py
0462               pz0(i) = pz
0463               e0(i) = e
0464            end if
0465  1001   continue
0466 
0467         do 1002 i = mul + 1, mul + incmul
0468            if (iconfg .le. 3) then
0469               gz0(i) = tau0 * sinh(eta(i))
0470               ft0(i) = tau0 * cosh(eta(i))
0471               if (iconfg .eq. 1) then
0472                  call posit1(x, y, r0)
0473                  gx0(i) = x + px0(i) * ft0(i)/e0(i)
0474                  gy0(i) = y + py0(i) * ft0(i)/e0(i)
0475               else if (iconfg .eq. 2 .or. iconfg .eq. 3) then
0476                  call posit2(x, y)
0477                  gx0(i) = x
0478                  gy0(i) = y
0479               end if
0480            else
0481               ft0(i) = 0d0
0482               call posit3(x, y, z)
0483               gx0(i) = x
0484               gy0(i) = y
0485               gz0(i) = z
0486            end if
0487  1002   continue
0488 
0489         mul = mul + incmul
0490             
0491 c       check if it's necessary to adjust array size 'adarr'
0492             if (mul .ge. MAXPTN .or. mul .eq. 0) then
0493 cms           print *, 'event',ievt,'has',mul,'number of gluon',
0494 cms     &          'adjusting counting is necessary'
0495            stop 'adarr'
0496         end if
0497         
0498         return
0499         end
0500 
0501         subroutine posit1(x, y, r0)
0502         
0503         implicit double precision (a-h, o-z)
0504 
0505 cc        external ran1
0506         common /rndm3/ iseedp
0507 cc      SAVE /rndm3/
0508         SAVE   
0509 
0510         iseed=iseedp
0511  10        x = 2d0 * ran1(iseed) - 1d0
0512         y = 2d0 * ran1(iseed) - 1d0
0513         if (x ** 2 + y ** 2 .gt. 1d0) goto 10
0514         x = x * r0
0515         y = y * r0
0516         
0517         return
0518         end
0519 
0520         subroutine posit2(x, y)
0521         
0522         implicit double precision (a-h, o-z)
0523 
0524 cc        external ran1
0525 
0526         common /ilist3/ size1, size2, size3, v1, v2, v3, size
0527 cc      SAVE /ilist3/
0528         common /rndm3/ iseedp
0529 cc      SAVE /rndm3/
0530         SAVE   
0531         iseed=iseedp
0532          x = 2d0 * ran1(iseed) - 1d0
0533         y = 2d0 * ran1(iseed) - 1d0
0534         x = x * 5d0 * size1
0535         y = y * 5d0 * size2
0536         
0537         return
0538         end
0539 
0540         subroutine posit3(x, y, z)
0541         
0542         implicit double precision (a-h, o-z)
0543 
0544 cc        external ran1
0545 
0546         common /ilist3/ size1, size2, size3, v1, v2, v3, size
0547 cc      SAVE /ilist3/
0548         common /rndm3/ iseedp
0549 cc      SAVE /rndm3/
0550         SAVE   
0551 
0552         iseed=iseedp
0553          x = 2d0 * ran1(iseed) - 1d0
0554         y = 2d0 * ran1(iseed) - 1d0
0555         z = 2d0 * ran1(iseed) - 1d0
0556         x = x * 5d0 * size1
0557         y = y * 5d0 * size2
0558         z = z * 5d0 * size3
0559         
0560         return
0561         end
0562         
0563         subroutine energy(e, temp)
0564 
0565 c       to generate the magnitude of the momentum e,
0566 c       knowing the temperature of the local thermal distribution temp
0567         
0568         implicit double precision (a-h, o-z)
0569         
0570 cc        external ran1
0571 
0572         common /para2/ xmp, xmu, alpha, rscut2, cutof2
0573 cc      SAVE /para2/
0574         common /rndm3/ iseedp
0575 cc      SAVE /rndm3/
0576         SAVE   
0577 
0578         iseed=iseedp
0579  1000        continue
0580         
0581         e = ran1(iseed)
0582         e = e * ran1(iseed)
0583         e = e * ran1(iseed)
0584 
0585         if (e .le. 0d0) goto 1000
0586         e = - temp * log(e)
0587         if (ran1(iseed) .gt. 
0588      &     exp((e - dsqrt(e ** 2 + xmp ** 2))/temp)) then
0589            goto 1000
0590         end if
0591 
0592         return
0593         end
0594         
0595         subroutine momntm(px, py, pz, e)
0596 
0597 c       to generate the 3 components of the momentum px, py, pz,
0598 c       from the magnitude of the momentum e
0599         
0600         implicit double precision (a-h,o-z)
0601         
0602 cc        external ran1
0603         
0604         parameter (pi = 3.14159265358979d0)
0605         common /rndm3/ iseedp
0606 cc      SAVE /rndm3/
0607         SAVE   
0608 
0609         iseed=iseedp
0610         cost = 2d0 * ran1(iseed) - 1d0
0611 c     7/20/01:
0612 c        sint = sqrt(1d0 - cost ** 2)
0613         sint = dsqrt(1d0 - cost ** 2)
0614         phi = 2d0 * pi * ran1(iseed)
0615       
0616         px = e * sint * cos(phi)
0617         py = e * sint * sin(phi)
0618         pz = e * cost
0619         
0620         return
0621         end
0622 
0623         subroutine boosti
0624 
0625         implicit double precision (a-h, o-z)
0626         parameter (MAXPTN=400001)
0627         common /para1/ mul
0628 cc      SAVE /para1/
0629         common /para6/ centy
0630 cc      SAVE /para6/
0631         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
0632      &       PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
0633      &       XMASS0(MAXPTN), ITYP0(MAXPTN)
0634 cc      SAVE /prec1/
0635         common /lor/ enenew, pxnew, pynew, pznew
0636 cc      SAVE /lor/
0637         SAVE   
0638 
0639         external lorenz
0640 
0641         bex = 0d0 
0642         bey = 0d0
0643         bez = - tanh(centy)
0644         
0645 c       save data for many runs of the same initial condition           
0646         do 1001 i = 1, mul
0647            px1 = gx0(i)
0648            py1 = gy0(i)
0649            pz1 = gz0(i)
0650            e1 = ft0(i)
0651            call lorenz(e1, px1, py1, pz1, bex, bey, bez)
0652            gx0(i) = pxnew
0653            gy0(i) = pynew
0654            gz0(i) = pznew
0655            ft0(i) = enenew
0656            px1 = px0(i)
0657            py1 = py0(i)
0658            pz1 = pz0(i)
0659            e1 = e0(i)
0660            call lorenz(e1, px1, py1, pz1, bex, bey, bez)
0661            px0(i) = pxnew
0662            py0(i) = pynew
0663            pz0(i) = pznew
0664            e0(i) = enenew
0665  1001   continue
0666         
0667         return
0668         end
0669 
0670 ******************************************************************************
0671 
0672         subroutine inirun
0673         SAVE   
0674 
0675 c       sort prec2 according to increasing formation time
0676         call ftime
0677         call inirec
0678         call iilist
0679         call inian2
0680 
0681         return
0682         end
0683 
0684         subroutine ftime
0685 c       this subroutine generates formation time for the particles
0686 c       indexing ft(i)
0687 c       input e(i)
0688 c       output ft(i), indx(i)
0689 
0690         implicit double precision (a-h, o-z)
0691 
0692         external ftime1
0693         parameter (MAXPTN=400001)
0694         common /para1/ mul
0695 cc      SAVE /para1/
0696         common /para2/ xmp, xmu, alpha, rscut2, cutof2
0697 cc      SAVE /para2/
0698         common /para4/ iftflg, ireflg, igeflg, ibstfg
0699 cc      SAVE /para4/
0700         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
0701      &       PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
0702      &       XMASS0(MAXPTN), ITYP0(MAXPTN)
0703 cc      SAVE /prec1/
0704         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
0705 cc      SAVE /prec4/
0706         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
0707 cc      SAVE /ilist4/
0708         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
0709 cc      SAVE /ilist5/
0710         common /par1/ formt
0711 cc      SAVE /par1/
0712         common/anim/nevent,isoft,isflag,izpc
0713 cc      SAVE /anim/
0714         common /rndm3/ iseedp
0715 cc      SAVE /rndm3/
0716         SAVE   
0717 
0718         iseed=iseedp
0719 clin-6/07/02 initialize here to expedite compiling, instead in zpcbdt:
0720         do 1001 i = 1, MAXPTN
0721            ct(i)=0d0
0722            ot(i)=0d0
0723  1001   continue
0724         tlarge=1000000.d0
0725 clin-6/07/02-end
0726 
0727         if (iftflg .eq. 0) then
0728 c     5/01/01 different prescription for parton initial formation time:
0729            if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
0730               do 1002 i = 1, mul
0731                  if (ft0(i) .gt. tlarge) ft0(i) = tlarge
0732  1002         continue
0733               goto 150
0734            else
0735 c     5/01/01-end
0736 
0737            do 1003 i = 1, MAXPTN
0738               ft0(i) = tlarge
0739  1003      continue
0740            do 1004 i = 1, mul
0741               xmt2 = px0(i) ** 2 + py0(i) ** 2 + xmp ** 2
0742               formt = xmt2 / e0(i)           
0743               ft0(i) = ftime1(iseed)
0744               if (ft0(i) .gt. tlarge) ft0(i) = tlarge
0745  1004      continue
0746 c     5/01/01:
0747         endif
0748 
0749         end if
0750 
0751 c     5/01/01:
0752  150        continue
0753 
0754 c        call index1(MAXPTN, mul, ft0, indx)
0755         if (mul .gt. 1) then
0756            call index1(MAXPTN, mul, ft0, indx)
0757         else
0758 clin-7/09/03: need to set value for mul=1:
0759            indx(1)=1
0760         end if
0761 c
0762         return
0763         end
0764 
0765         subroutine inirec
0766 
0767         implicit double precision (a-h, o-z)
0768 cc        external ran1
0769         parameter (MAXPTN=400001)
0770         common /para1/ mul
0771 cc      SAVE /para1/
0772         common /para4/ iftflg, ireflg, igeflg, ibstfg
0773 cc      SAVE /para4/
0774         common /para5/ iconfg, iordsc
0775 cc      SAVE /para5/
0776         COMMON /prec1/GX0(MAXPTN),GY0(MAXPTN),GZ0(MAXPTN),FT0(MAXPTN),
0777      &       PX0(MAXPTN), PY0(MAXPTN), PZ0(MAXPTN), E0(MAXPTN),
0778      &       XMASS0(MAXPTN), ITYP0(MAXPTN)
0779 cc      SAVE /prec1/
0780         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
0781      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
0782      &       xmass(MAXPTN), ityp(MAXPTN)
0783 cc      SAVE /prec2/
0784         common /prec3/gxs(MAXPTN),gys(MAXPTN),gzs(MAXPTN),fts(MAXPTN),
0785      &     pxs(MAXPTN), pys(MAXPTN), pzs(MAXPTN), es(MAXPTN),
0786      &     xmasss(MAXPTN), ityps(MAXPTN)
0787 cc      SAVE /prec3/
0788         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
0789 cc      SAVE /prec4/
0790         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
0791 cc      SAVE /prec5/
0792         common /prec6/ etas(MAXPTN), raps(MAXPTN), taus(MAXPTN)
0793 cc      SAVE /prec6/
0794         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
0795 cc      SAVE /ilist4/
0796 cbz1/25/99
0797         COMMON /ilist7/ LSTRG0(MAXPTN), LPART0(MAXPTN)
0798 cc      SAVE /ilist7/
0799         COMMON /ilist8/ LSTRG1(MAXPTN), LPART1(MAXPTN)
0800 cc      SAVE /ilist8/
0801 cbz1/25/99end
0802         COMMON /smearz/smearp,smearh
0803 cc      SAVE /smearz/
0804         dimension vxp(MAXPTN), vyp(MAXPTN), vzp(MAXPTN)
0805         common /precpa/ vxp0(MAXPTN), vyp0(MAXPTN), vzp0(MAXPTN)
0806 cc      SAVE /precpa/
0807         common/anim/nevent,isoft,isflag,izpc
0808 cc      SAVE /anim/
0809 clin-6/06/02 local parton freezeout:
0810         common /frzprc/ 
0811      &       gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
0812      &       pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
0813      &       xmfrz(MAXPTN), 
0814      &       tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
0815 cc      SAVE /frzprc/
0816         common /rndm3/ iseedp
0817 cc      SAVE /rndm3/
0818         common /para7/ ioscar,nsmbbbar,nsmmeson
0819         COMMON /AREVT/ IAEVT, IARUN, MISS
0820         SAVE   
0821         iseed=iseedp
0822 clin-6/06/02 local freezeout initialization:
0823         if(isoft.eq.5) then
0824            itlast=0
0825            call inifrz
0826         endif
0827 
0828         do 1001 i = 1, mul
0829 clin-7/09/01 define indx(i) to save time:
0830 c           ityp(i) = ityp0(indx(i))
0831 c           gx(i) = gx0(indx(i))
0832 c           gy(i) = gy0(indx(i))
0833 c           gz(i) = gz0(indx(i))
0834 c           ft(i) = ft0(indx(i))
0835 c           px(i) = px0(indx(i))
0836 c           py(i) = py0(indx(i))
0837 c           pz(i) = pz0(indx(i))
0838 c           e(i) = e0(indx(i))
0839 c           xmass(i) = xmass0(indx(i))
0840 ccbz1/25/99
0841 c           LSTRG1(I) = LSTRG0(INDX(I))
0842 c           LPART1(I) = LPART0(INDX(I))
0843 ccbz1/25/99end
0844            indxi=indx(i)
0845            ityp(i) = ityp0(indxi)
0846            gx(i) = gx0(indxi)
0847            gy(i) = gy0(indxi)
0848            gz(i) = gz0(indxi)
0849            ft(i) = ft0(indxi)
0850            px(i) = px0(indxi)
0851            py(i) = py0(indxi)
0852            pz(i) = pz0(indxi)
0853            e(i) = e0(indxi)
0854            xmass(i) = xmass0(indxi)
0855            LSTRG1(I) = LSTRG0(INDXI)
0856            LPART1(I) = LPART0(INDXI)
0857            vxp(I) = vxp0(INDXI)
0858            vyp(I) = vyp0(INDXI)
0859            vzp(I) = vzp0(INDXI)
0860 clin-7/09/01-end
0861 c
0862 clin-6/06/02 local freezeout initialization:
0863          if(isoft.eq.5) then
0864             idfrz(i)=ityp(i)
0865             gxfrz(i)=gx(i)
0866             gyfrz(i)=gy(i)
0867             gzfrz(i)=gz(i)
0868             ftfrz(i)=ft(i)
0869             pxfrz(i)=px(i)
0870             pyfrz(i)=py(i)
0871             pzfrz(i)=pz(i)
0872             efrz(i)=e(i)
0873             xmfrz(i)=xmass(i)
0874             ifrz(i)=0
0875          endif
0876 clin-6/06/02-end
0877  1001 continue
0878 
0879 c       save particle info for fixed time analysis
0880         do 1002 i = 1, mul
0881            ityps(i) = ityp(i)
0882            gxs(i) = gx(i)
0883            gys(i) = gy(i)
0884            gzs(i) = gz(i)
0885            fts(i) = ft(i)
0886            pxs(i) = px(i)
0887            pys(i) = py(i)
0888            pzs(i) = pz(i)
0889            es(i) = e(i)
0890            xmasss(i) = xmass(i)
0891  1002   continue
0892 
0893 clin-6/2009
0894         if(isoft.eq.1.and.(ioscar.eq.2.or.ioscar.eq.3))
0895      1       write(92,*) iaevt,miss,mul
0896 
0897         do 1003 i = 1, mul
0898            energy = e(i)
0899            vx(i) = px(i) / energy
0900            vy(i) = py(i) / energy
0901            vz(i) = pz(i) / energy
0902            if (iftflg .eq. 0) then
0903               formt = ft(i)
0904 c     7/09/01 propagate partons with parent velocity till formation
0905 c     so that partons in same hadron have 0 distance:
0906 c            gx(i) = gx(i) + vx(i) * formt
0907 c            gy(i) = gy(i) + vy(i) * formt
0908 c            gz(i) = gz(i) + vz(i) * formt
0909             if(isoft.eq.3.or.isoft.eq.4.or.isoft.eq.5) then
0910                gx(i) = gx(i) + vxp(i) * formt
0911                gy(i) = gy(i) + vyp(i) * formt
0912                gz(i) = gz(i) + vzp(i) * formt
0913             else
0914                gx(i) = gx(i) + vx(i) * formt
0915                gy(i) = gy(i) + vy(i) * formt
0916                gz(i) = gz(i) + vz(i) * formt
0917             endif
0918 c     7/09/01-end
0919 c
0920 c     3/27/00-ctest off no smear z on partons to avoid eta overflow:
0921 c              gz(i) = gz(i)+smearp*(2d0 * ran1(iseed) - 1d0)
0922 c     to give eta=y +- smearp*random:
0923 c              smeary=smearp*(2d0 * ran1(iseed) - 1d0)
0924 c              smearf=dexp(2*smeary)*(1+vz(i))/(1-vz(i)+1.d-8)
0925 c              gz(i) = gz(i)+formt*(smearf-1)/(smearf+1)
0926 c     3/27/00-end
0927            end if
0928 
0929 clin-6/2009 write out initial parton information after string melting
0930 c     and after propagating to its format time:
0931            if(ioscar.eq.2.or.ioscar.eq.3) then
0932               if(dmax1(abs(gx(i)),abs(gy(i)),
0933      1             abs(gz(i)),abs(ft(i))).lt.9999) then
0934                  write(92,200) ityp(i),px(i),py(i),pz(i),xmass(i),
0935      1                gx(i),gy(i),gz(i),ft(i)
0936               else
0937                  write(92,201) ityp(i),px(i),py(i),pz(i),xmass(i),
0938      1                gx(i),gy(i),gz(i),ft(i)
0939               endif
0940            endif
0941  200       format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2))
0942  201       format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2))
0943 c
0944  1003   continue
0945 
0946         if (iconfg .le. 3) then
0947            do 1004 i = 1, mul
0948               if (ft(i) .le. abs(gz(i))) then
0949                  eta(i) = 1000000.d0
0950               else
0951                  eta(i) = 0.5d0 * log((ft(i) + gz(i)) / (ft(i) - gz(i)))
0952               end if
0953               if (e(i) .le. abs(pz(i))) then
0954                  rap(i) = 1000000.d0
0955               else
0956                  rap(i) = 0.5d0 * log((e(i) + pz(i)) / (e(i) - pz(i)))
0957               end if
0958               tau(i) = ft(i) / cosh(eta(i))
0959  1004      continue
0960            
0961            do 1005 i = 1, mul
0962               etas(i) = eta(i)
0963               raps(i) = rap(i)
0964               taus(i) = tau(i)
0965  1005      continue
0966         end if
0967 
0968         return
0969         end
0970 
0971         subroutine iilist
0972 
0973         implicit double precision (a-h, o-z)
0974         parameter (MAXPTN=400001)
0975         common /para1/ mul
0976 cc      SAVE /para1/
0977         common /ilist1/
0978      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
0979      &     ictype, icsta(MAXPTN),
0980      &     nic(MAXPTN), icels(MAXPTN)
0981 cc      SAVE /ilist1/
0982         common /ilist2/ icell, icel(10,10,10)
0983 cc      SAVE /ilist2/
0984         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
0985 cc      SAVE /ilist4/
0986         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
0987 cc      SAVE /ilist5/
0988         common /ilist6/ t, iopern, icolln
0989 cc      SAVE /ilist6/
0990         SAVE   
0991 
0992         iscat = MAXPTN
0993         jscat = MAXPTN
0994 
0995         do 1001 i = 1, mul
0996            next(i) = 0
0997            last(i) = 0
0998            icsta(i) = 0
0999            nic(i) = 0
1000            icels(i) = 0
1001  1001   continue
1002 
1003         icell = 0
1004         do 1004 i1 = 1, 10
1005            do 1003 i2 = 1, 10
1006               do 1002 i3 = 1, 10
1007                  icel(i1, i2, i3) = 0
1008  1002         continue
1009  1003      continue
1010  1004   continue
1011 
1012         ichkpt = 0
1013         ifmpt = 1
1014 
1015         do 1005 i = 1, mul
1016            ct(i) = tlarge
1017            ot(i) = tlarge
1018  1005   continue
1019 
1020         iopern = 0
1021         icolln = 0
1022         t = 0.d0
1023 
1024         return
1025         end
1026 
1027         subroutine inian2
1028 
1029         implicit double precision (a-h, o-z)
1030 
1031         common /para5/ iconfg, iordsc
1032 cc      SAVE /para5/
1033         common /ana2/
1034      &     det(12), dn(12), detdy(12), detdn(12), dndy(12),
1035      &     det1(12), dn1(12), detdy1(12), detdn1(12), dndy1(12),
1036      &     det2(12), dn2(12), detdy2(12), detdn2(12), dndy2(12)
1037 cc      SAVE /ana2/
1038         SAVE   
1039 
1040         if (iconfg .le. 3) then
1041            do 1001 i = 1, 12
1042               det(i) = 0d0
1043               dn(i) = 0d0
1044               det1(i) = 0d0
1045               dn1(i) = 0d0
1046               det2(i) = 0d0
1047               dn2(i) = 0d0
1048  1001      continue
1049         end if
1050 
1051         return
1052         end
1053 
1054 ******************************************************************************
1055 ******************************************************************************
1056 
1057         subroutine zpcrun(*)
1058 
1059         implicit double precision (a-h, o-z)
1060         parameter (MAXPTN=400001)
1061         parameter (tend1 = 250d0)
1062         parameter (tend2 = 6.1d0)
1063         common /para1/ mul
1064 cc      SAVE /para1/
1065         common /para5/ iconfg, iordsc
1066         common /para7/ ioscar,nsmbbbar,nsmmeson
1067 cc      SAVE /para5/
1068         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1069      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1070      &       xmass(MAXPTN), ityp(MAXPTN)
1071 cc      SAVE /prec2/
1072         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1073 cc      SAVE /prec4/
1074         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1075 cc      SAVE /prec5/
1076         common /ilist1/
1077      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1078      &     ictype, icsta(MAXPTN),
1079      &     nic(MAXPTN), icels(MAXPTN)
1080 cc      SAVE /ilist1/
1081         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
1082 cc      SAVE /ilist4/
1083         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1084 cc      SAVE /ilist5/
1085         common /ilist6/ t, iopern, icolln
1086 cc      SAVE /ilist6/
1087         common /ana1/ ts(12)
1088 cc      SAVE /ana1/
1089         common/anim/nevent,isoft,isflag,izpc
1090 cc      SAVE /anim/
1091         COMMON /AREVT/ IAEVT, IARUN, MISS
1092         SAVE   
1093 
1094 c       save last collision info
1095         if (mod(ictype, 2) .eq. 0) then
1096            call savrec(iscat)
1097            call savrec(jscat)
1098         end if
1099 
1100 c1      get operation type
1101         call getict(t1)
1102 c2      check freezeout condition
1103         if (iconfg .eq. 1 .and. t1 .gt. tlarge / 2d0) return 1
1104         if (iconfg .eq. 2 .or. iconfg .eq. 3) then
1105            if (t1 .gt. tend1) return 1
1106 c           if (ichkpt .eq. mul) then
1107 c              ii = 0
1108 c              do i = 1, mul
1109 c                 gztemp = gz(i) + vz(i) * (t1 - ft(i))
1110 c                 if (sqrt(t1 ** 2 - gztemp ** 2) .lt. tend) then
1111 c                    ii = 1
1112 c                    goto 1000
1113 c                 end if
1114 c              end do
1115 c 1000              continue
1116 c              if (ii .eq. 0) return 1
1117 c           end if
1118         end if
1119         if (iconfg .eq. 4 .or. iconfg .eq. 5) then
1120            if (t1 .gt. tend2) return 1
1121         end if
1122 
1123 clin-6/06/02 local freezeout for string melting,
1124 c     decide what partons have frozen out at time t1:
1125       if(isoft.eq.5) then
1126          call local(t1)
1127       endif
1128 
1129 c3      update iopern, t
1130 
1131         iopern = iopern + 1
1132         t = t1
1133         if (mod(ictype, 2) .eq. 0) then
1134            icolln = icolln + 1
1135 
1136 c     4/18/01-ctest off
1137 c           write (2006, 1233) 'iscat=', iscat, 'jscat=', jscat,
1138 c           write (2006, *) 'iscat=', iscat, ' jscat=', jscat,
1139 c     1 ityp(iscat), ityp(jscat)
1140 c           write (2006, 1233) 'iscat=', max(indx(iscat), indx(jscat)),
1141 c     &        'jscat=', min(indx(iscat), indx(jscat))
1142 
1143 c           write (2006, 1234) ' icolln=', icolln, 't=', t
1144 
1145 c 1233           format (a10, i10, a10, i10)
1146 c 1234           format (a15, i10, a5, f23.17, a5, f23.17)
1147         end if
1148 
1149 c4.1    deal with formation
1150         if (iconfg .eq. 1
1151      &     .or. iconfg .eq. 2
1152      &     .or. iconfg .eq. 4) then
1153            if (ictype .eq. 1 .or. ictype .eq. 2 .or. 
1154      &        ictype .eq. 5 .or. ictype .eq. 6) then
1155               call celasn
1156            end if
1157         end if
1158 
1159 c4.2    deal with collisions
1160 
1161         if (ictype .ne. 1) then
1162 
1163            iscat0 = iscat
1164            jscat0 = jscat
1165            
1166 c        iscat is the larger one so that if it's a wall collision,
1167 c       it's still ok
1168            iscat = max0(iscat0, jscat0)
1169            jscat = min0(iscat0, jscat0)
1170 
1171 ctest off check icsta(i): 0 with f77 compiler
1172 c        write(9,*) 'BB:ictype,t1,iscat,jscat,icsta(i)=',
1173 c     1 ictype,t1,iscat,jscat,icsta(iscat)
1174            
1175 c       check collision time table error 'tterr'
1176 clin-4/2008 to avoid out-of-bound error in next():
1177 c           if (jscat .ne. 0 .and. next(jscat) .ne. iscat)
1178 c     &        then
1179 c              print *, 'iscat=', iscat, 'jscat=', jscat,
1180 c     &             'next(', jscat, ')=', next(jscat)
1181 c
1182 c              if (ct(iscat) .lt. tlarge / 2d0) stop 'tterr'
1183 c              if (ct(jscat) .lt. tlarge / 2d0) stop 'tterr'
1184 c           end if 
1185            if (jscat .ne. 0) then
1186               if(next(jscat) .ne. iscat) then
1187                  print *, 'iscat=', iscat, 'jscat=', jscat,
1188      &                'next(', jscat, ')=', next(jscat)
1189                  if (ct(iscat) .lt. tlarge / 2d0) stop 'tterr'
1190                  if (ct(jscat) .lt. tlarge / 2d0) stop 'tterr'
1191               endif
1192            end if 
1193 clin-4/2008-end
1194            
1195 c4.2.1     collisions with wall
1196 
1197 c     8/19/02 avoid actual argument in common blocks of cellre:
1198          niscat=iscat
1199          njscat=jscat
1200 c           if (icsta(iscat) .ne. 0) call cellre(iscat, t)
1201 c           if (jscat .ne. 0) then
1202 c              if (icsta(jscat) .ne. 0) call cellre(jscat, t)
1203 c           end if
1204            if (icsta(iscat) .ne. 0) call cellre(niscat, t)
1205            if (jscat .ne. 0) then
1206               if (icsta(jscat) .ne. 0) call cellre(njscat, t)
1207            end if
1208 
1209 c4.2.2     collision between particles     
1210 
1211 clin-6/2009 write out info for each collision:
1212 c           if (mod(ictype, 2) .eq. 0) call scat(t, iscat, jscat)
1213            if (mod(ictype, 2) .eq. 0) then
1214               if(ioscar.eq.3) then
1215             write(95,*) 'event,miss,iscat,jscat=',iaevt,miss,iscat,jscat
1216                  if(dmax1(abs(gx(iscat)),abs(gy(iscat)),
1217      1                abs(gz(iscat)),abs(ft(iscat)),abs(gx(jscat)),
1218      2                abs(gy(jscat)),abs(gz(jscat)),abs(ft(jscat)))
1219      3                .lt.9999) then
1220                     write(95,200) ityp(iscat),px(iscat),py(iscat),
1221      1                   pz(iscat),xmass(iscat),gx(iscat),gy(iscat),
1222      2                   gz(iscat),ft(iscat)
1223                     write(95,200) ityp(jscat),px(jscat),py(jscat),
1224      1                   pz(jscat),xmass(jscat),gx(jscat),gy(jscat),
1225      2                   gz(jscat),ft(jscat)
1226                  else
1227                     write(95,201) ityp(iscat),px(iscat),py(iscat),
1228      1                   pz(iscat),xmass(iscat),gx(iscat),gy(iscat),
1229      2                   gz(iscat),ft(iscat)
1230                     write(95,201) ityp(jscat),px(jscat),py(jscat),
1231      1                   pz(jscat),xmass(jscat),gx(jscat),gy(jscat),
1232      2                   gz(jscat),ft(jscat)
1233                  endif
1234               endif
1235 c     
1236               call scat(t, iscat, jscat)
1237 c     
1238               if(ioscar.eq.3) then
1239                  if(dmax1(abs(gx(iscat)),abs(gy(iscat)),
1240      1                abs(gz(iscat)),abs(ft(iscat)),abs(gx(jscat)),
1241      2                abs(gy(jscat)),abs(gz(jscat)),abs(ft(jscat)))
1242      3                .lt.9999) then
1243                     write(95,200) ityp(iscat),px(iscat),py(iscat),
1244      1                   pz(iscat),xmass(iscat),gx(iscat),gy(iscat),
1245      2                   gz(iscat),ft(iscat)
1246                     write(95,200) ityp(jscat),px(jscat),py(jscat),
1247      1                   pz(jscat),xmass(jscat),gx(jscat),gy(jscat),
1248      2                   gz(jscat),ft(jscat)
1249                  else
1250                     write(95,201) ityp(iscat),px(iscat),py(iscat),
1251      1                   pz(iscat),xmass(iscat),gx(iscat),gy(iscat),
1252      2                   gz(iscat),ft(iscat)
1253                     write(95,201) ityp(jscat),px(jscat),py(jscat),
1254      1                   pz(jscat),xmass(jscat),gx(jscat),gy(jscat),
1255      2                   gz(jscat),ft(jscat)
1256                  endif
1257               endif
1258            endif
1259  200       format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,f8.2))
1260  201       format(I6,2(1x,f8.3),1x,f10.3,1x,f6.3,4(1x,e8.2))
1261            
1262         end if
1263 
1264 c5      update the interaction list
1265         call ulist(t)
1266 
1267 c6      update ifmpt. ichkpt
1268 c       old ichkpt and ifmpt are more conveniently used in ulist
1269         if (ifmpt .le. mul) then
1270            if (ictype .ne. 0 .and. ictype .ne. 3 
1271      &        .and. ictype .ne. 4) then
1272               ichkpt = ichkpt + 1
1273               ifmpt = ifmpt + 1
1274            end if
1275         end if
1276 
1277         return
1278         end
1279 
1280         subroutine savrec(i)
1281 
1282         implicit double precision (a-h, o-z)
1283         parameter (MAXPTN=400001)
1284         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1285      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1286      &       xmass(MAXPTN), ityp(MAXPTN)
1287 cc      SAVE /prec2/
1288         common /prec3/gxs(MAXPTN),gys(MAXPTN),gzs(MAXPTN),fts(MAXPTN),
1289      &     pxs(MAXPTN), pys(MAXPTN), pzs(MAXPTN), es(MAXPTN),
1290      &     xmasss(MAXPTN), ityps(MAXPTN)
1291 cc      SAVE /prec3/
1292         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1293 cc      SAVE /prec5/
1294         common /prec6/ etas(MAXPTN), raps(MAXPTN), taus(MAXPTN)
1295 cc      SAVE /prec6/
1296         SAVE   
1297 
1298         ityps(i) = ityp(i)
1299         gxs(i) = gx(i)
1300         gys(i) = gy(i)
1301         gzs(i) = gz(i)
1302         fts(i) = ft(i)
1303         pxs(i) = px(i)
1304         pys(i) = py(i)
1305         pzs(i) = pz(i)
1306         es(i) = e(i)
1307         xmasss(i) = xmass(i)
1308         etas(i) = eta(i)
1309         raps(i) = rap(i)
1310         taus(i) = tau(i)
1311 
1312         return
1313         end
1314 
1315         subroutine getict(t1)
1316         implicit double precision (a-h, o-z)
1317         parameter (MAXPTN=400001)
1318         common /para1/ mul
1319 cc      SAVE /para1/
1320         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1321      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1322      &       xmass(MAXPTN), ityp(MAXPTN)
1323 cc      SAVE /prec2/
1324         common /ilist1/
1325      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1326      &     ictype, icsta(MAXPTN),
1327      &     nic(MAXPTN), icels(MAXPTN)
1328 cc      SAVE /ilist1/
1329         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
1330 cc      SAVE /ilist4/
1331         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1332 cc      SAVE /ilist5/
1333         SAVE   
1334 
1335 c       neglect possibility of 2 collisions at the same time
1336 c0       set initial conditions
1337 
1338         t1 = tlarge
1339         iscat = 0
1340         jscat = 0
1341 
1342 c1      get next collision between particles
1343         do 1001 i = 1, ichkpt
1344            if (ot(i) .lt. t1) then
1345               t1 = ot(i)
1346               iscat = i
1347            end if
1348  1001   continue
1349         if (iscat .ne. 0) jscat = next(iscat)
1350 
1351 c2      get ictype
1352 c     10/30/02 ictype=0:collision; 1:parton formation
1353         if (iscat .ne. 0 .and. jscat .ne. 0) then
1354            if (icsta(iscat) .eq. 0 .and. icsta(jscat) .eq. 0) then
1355               ictype = 0
1356            else
1357               ictype = 4
1358            end if
1359         else if (iscat .ne. 0 .or. jscat .ne. 0) then
1360            ictype = 3
1361         end if
1362 c
1363         if (ifmpt .le. mul) then
1364            if (ft(ifmpt) .lt. t1) then
1365               ictype = 1
1366               t1 = ft(ifmpt)
1367            else if (ft(ifmpt) .eq. t1) then
1368               if (ictype .eq. 0) ictype = 2
1369               if (ictype .eq. 3) ictype = 5
1370               if (ictype .eq. 4) ictype = 6
1371            end if
1372         end if
1373 
1374         return
1375         end
1376 
1377         subroutine celasn
1378 c       this subroutine is used to assign a cell for a newly formed particle
1379 c       output: nic(MAXPTN) icels(MAXPTN) in the common /ilist1/
1380 c       icell, and icel(10,10,10) in the common /ilist2/
1381 
1382         implicit double precision (a-h, o-z)
1383         parameter (MAXPTN=400001)
1384         common /para1/ mul
1385 cc      SAVE /para1/
1386         common /para5/ iconfg, iordsc
1387 cc      SAVE /para5/
1388         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1389      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1390      &       xmass(MAXPTN), ityp(MAXPTN)
1391 cc      SAVE /prec2/
1392         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1393 cc      SAVE /prec4/
1394         common /ilist1/
1395      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1396      &     ictype, icsta(MAXPTN),
1397      &     nic(MAXPTN), icels(MAXPTN)
1398 cc      SAVE /ilist1/
1399         common /ilist2/ icell, icel(10,10,10)
1400 cc      SAVE /ilist2/
1401         common /ilist3/ size1, size2, size3, v1, v2, v3, size
1402 cc      SAVE /ilist3/
1403         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
1404 cc      SAVE /ilist4/
1405         SAVE   
1406 
1407         external integ
1408 
1409         i = ifmpt
1410         tt = ft(i)
1411         td = tt - size
1412         if (iconfg .eq. 1 .and. (size1 .eq. 0d0 .or.
1413      &     size2 .eq. 0d0 .or. size3 .eq. 0d0)) then
1414            i1 = 11
1415            i2 = 11
1416            i3 = 11
1417         else if (iconfg .eq. 4 .or. td .le. 0d0) then
1418            i1 = integ(gx(i) / size1) + 6
1419            i2 = integ(gy(i) / size2) + 6
1420            i3 = integ(gz(i) / size3) + 6
1421            if (integ(gx(i) / size1) .eq. gx(i) / size1 .and. 
1422      &        vx(i) .lt. 0d0)
1423      &        i1 = i1 - 1
1424            if (integ(gy(i) / size2) .eq. gy(i) / size2 .and. 
1425      &        vy(i) .lt. 0d0)
1426      &        i2 = i2 - 1
1427            if (integ(gz(i) / size3) .eq. gz(i) / size3 .and. 
1428      &        vz(i) .lt. 0d0)
1429      &        i3 = i3 - 1
1430         else
1431            i1 = integ(gx(i) / (size1 + v1 * td)) + 6
1432            i2 = integ(gy(i) / (size2 + v2 * td)) + 6
1433            i3 = integ(gz(i) / (size3 + v3 * td)) + 6
1434            if (integ(gx(i) / (size1 + v1 * td)) .eq. gx(i) / 
1435      &        (size1 + v1 * td) .and. vx(i) .lt. (i1 - 6) * v1)
1436      &        i1 = i1 - 1
1437            if (integ(gy(i) / (size2 + v2 * td)) .eq. gy(i)/
1438      &        (size2 + v2 * td) .and. vy(i) .lt. (i2 - 6) * v2)
1439      &        i2 = i2 - 1
1440            if (integ(gz(i) / (size3 + v3 * td)) .eq. gz(i)/
1441      &        (size3 + v3 * td) .and. vz(i) .lt. (i3 - 6) * v3)
1442      &        i3 = i3 - 1
1443         end if
1444 
1445         if (i1 .le. 0 .or. i1 .ge. 11 .or. i2 .le. 0 .or.
1446      &     i2 .ge. 11 .or. i3 .le. 0 .or. i3 .ge. 11) then
1447            i1 = 11
1448            i2 = 11
1449            i3 = 11
1450         end if
1451 
1452         if (i1 .eq. 11) then
1453            j = icell
1454            call newcre(i, j)
1455            icell = j
1456            icels(i) = 111111
1457         else
1458            j = icel(i1, i2, i3)
1459            call newcre(i, j)
1460            icel(i1, i2, i3) = j
1461            icels(i) = i1 * 10000 + i2 * 100 + i3
1462         end if
1463 
1464         return
1465         end
1466 
1467         integer function integ(x)
1468 c       this function is used to get the largest integer that is smaller than
1469 c       x
1470 
1471         implicit double precision (a-h, o-z)
1472         SAVE   
1473 
1474         if (x .lt. 0d0) then
1475            integ = int(x - 1d0)
1476         else
1477            integ = int( x )
1478         end if
1479 
1480         return
1481         end
1482 
1483         subroutine cellre(i, t)
1484 c       this subroutine is used for changing the cell of a particle that
1485 c       collide with the wall
1486 
1487         implicit double precision (a-h, o-z)
1488         parameter (MAXPTN=400001)
1489         common /para5/ iconfg, iordsc
1490 cc      SAVE /para5/
1491         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1492      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1493      &       xmass(MAXPTN), ityp(MAXPTN)
1494 cc      SAVE /prec2/
1495         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1496 cc      SAVE /prec4/
1497         common /aurec1/ jxa, jya, jza
1498 cc      SAVE /aurec1/
1499         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
1500 cc      SAVE /aurec2/
1501         common /ilist1/
1502      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1503      &     ictype, icsta(MAXPTN),
1504      &     nic(MAXPTN), icels(MAXPTN)
1505 cc      SAVE /ilist1/
1506         common /ilist2/ icell, icel(10,10,10)
1507 cc      SAVE /ilist2/
1508         common /ilist3/ size1, size2, size3, v1, v2, v3, size
1509 cc      SAVE /ilist3/
1510         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
1511 cc      SAVE /ilist4/
1512         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1513 cc      SAVE /ilist5/
1514         SAVE   
1515 
1516         logical good
1517 
1518         external integ
1519 
1520 c       this happens before update the /prec2/ common; in contrast with 
1521 c       scat which happens after updating the glue common
1522 
1523         t0 = t
1524 
1525  1000        continue
1526 
1527         if (iconfg .eq. 3 .or. iconfg .eq. 5) then
1528            k = mod(icsta(i), 10)
1529 
1530            if (k .eq. 1) then
1531               gx(i) = gx(i) - 10d0 * size1
1532               dgxa(i) = dgxa(i) + 10d0 * size1
1533               do 1001 ii = 1, ichkpt
1534                  if (next(ii) .eq. i) then
1535                     dgxa(ii) = dgxa(ii) - 10d0 * size1
1536                  end if
1537  1001         continue
1538            end if
1539            if (k .eq. 2) then
1540               gx(i) = gx(i) + 10d0 * size1
1541               dgxa(i) = dgxa(i) - 10d0 * size1
1542               do 1002 ii = 1, ichkpt
1543                  if (next(ii) .eq. i) then
1544                     dgxa(ii) = dgxa(ii) + 10d0 * size1
1545                  end if
1546  1002         continue
1547            end if
1548            if (k .eq. 3) then
1549               gy(i) = gy(i) - 10d0 * size2
1550               dgya(i) = dgya(i) + 10d0 * size2
1551               do 1003 ii = 1, ichkpt
1552                  if (next(ii) .eq. i) then
1553                     dgya(ii) = dgya(ii) - 10d0 * size2
1554                  end if
1555  1003         continue
1556            end if
1557            if (k .eq. 4) then
1558               gy(i) = gy(i) + 10d0 * size2
1559               dgya(i) = dgya(i) - 10d0 * size2
1560               do 1004 ii = 1, ichkpt
1561                  if (next(ii) .eq. i) then
1562                     dgya(ii) = dgya(ii) + 10d0 * size2
1563                  end if
1564  1004         continue
1565            end if
1566            if (iconfg .eq. 5) then
1567               if (k .eq. 5) then
1568                  gz(i) = gz(i) - 10d0 * size3
1569                  dgza(i) = dgza(i) + 10d0 * size3
1570                  do 1005 ii = 1, ichkpt
1571                     if (next(ii) .eq. i) then
1572                        dgza(ii) = dgza(ii) - 10d0 * size3
1573                     end if
1574  1005            continue
1575               end if
1576               if (k .eq. 6) then
1577                  gz(i) = gz(i) + 10d0 * size3
1578                  dgza(i) = dgza(i) - 10d0 * size3
1579                  do 1006 ii = 1, ichkpt
1580                     if (next(ii) .eq. i) then
1581                        dgza(ii) = dgza(ii) + 10d0 * size3
1582                     end if
1583  1006               continue
1584               end if
1585            end if
1586         else
1587            icels0 = icels(i)
1588 
1589            i1 = icels0 / 10000
1590            i2 = (icels0 - i1 * 10000) / 100
1591            i3 = icels0 - i1 * 10000 - i2 * 100
1592            
1593 cc       for particle inside the cube
1594            if (i1 .ge. 1 .and. i1 .le. 10
1595      &        .and. i2 .ge. 1 .and. i2 .le. 10
1596      &        .and. i3 .ge. 1 .and. i3 .le. 10) then
1597 
1598 c       this assignment takes care of nic(i)=0 automatically
1599               if (icel(i1, i2, i3) .eq. i) icel(i1, i2, i3) = nic(i)
1600 
1601 c1      rearrange the old cell
1602 
1603               call oldcre(i)
1604 
1605 c2      rearrange the new cell
1606 
1607               k = mod(icsta(i), 10)
1608               
1609 c2.1    particle goes out of the cube       
1610               if (iconfg .eq. 1) then
1611                  good = (i1 .eq. 1 .and. k .eq. 2)
1612      &              .or. (i1 .eq. 10 .and. k .eq. 1)
1613      &              .or. (i2 .eq. 1 .and. k .eq. 4)
1614      &              .or. (i2 .eq. 10 .and. k .eq. 3)
1615      &              .or. (i3 .eq. 1 .and. k .eq. 6)
1616      &              .or. (i3 .eq. 10 .and. k .eq. 5)
1617               end if
1618               if (iconfg .eq. 2) then
1619                  good = (i3 .eq. 1 .and. k .eq. 6)
1620      &              .or. (i3 .eq. 10 .and. k .eq. 5)
1621               end if
1622               if (good) then
1623 
1624 c                j = icell
1625                  call newcre(i, icell)
1626 c                 icell = j
1627 
1628                  icels(i) = 111111
1629 
1630 c2.2    particle moves inside the cube
1631               else
1632 
1633                  if (k .eq. 1) i1 = i1 + 1
1634                  if (k .eq. 2) i1 = i1 - 1
1635                  if (k .eq. 3) i2 = i2 + 1
1636                  if (k .eq. 4) i2 = i2 - 1
1637                  if (k .eq. 5) i3 = i3 + 1
1638                  if (k .eq. 6) i3 = i3 - 1
1639                  
1640                  if (iconfg .eq. 2 .or. iconfg .eq. 4) then
1641                     if (i1 .eq. 0) then
1642                        i1 = 10
1643                        gx(i) = gx(i) + 10d0 * size1
1644                     end if
1645                     if (i1 .eq. 11) then
1646                        i1 = 1
1647                        gx(i) = gx(i) - 10d0 * size1
1648                     end if
1649                     if (i2 .eq. 0) then
1650                        i2 = 10
1651                        gy(i) = gy(i) + 10d0 * size2
1652                     end if
1653                     if (i2 .eq. 11) then
1654                        i2 = 1
1655                        gy(i) = gy(i) - 10d0 * size2
1656                     end if
1657                     if (iconfg .eq. 4) then
1658                        if (i3 .eq. 0) then
1659                           i3 = 10
1660                           gz(i) = gz(i) + 10d0 * size3
1661                        end if
1662                        if (i3 .eq. 11) then
1663                           i3 = 1
1664                           gz(i) = gz(i) - 10d0 * size3
1665                        end if
1666                     end if
1667                  end if
1668                  
1669                  j = icel(i1, i2, i3)
1670                  
1671                  call newcre(i, j)
1672 c       in case icel changes
1673                  
1674                  icel(i1 ,i2, i3) = j
1675                  
1676                  icels(i) = i1 * 10000 + i2 * 100 + i3
1677                  
1678               end if
1679               
1680 cc       for particles outside the cube
1681            else
1682               
1683               if (icell .eq. i) icell = nic(i)
1684               
1685               call oldcre(i)
1686               
1687               k = mod(icsta(i), 10)
1688               
1689               ddt = t - ft(i)
1690               dtt = t - size
1691               if (dtt .le. 0d0) then
1692                  i1 = integ((gx(i) + vx(i) * ddt) / size1) + 6
1693                  i2 = integ((gy(i) + vy(i) * ddt) / size2) + 6
1694                  i3 = integ((gz(i) + vz(i) * ddt) / size3) + 6
1695               else
1696                  i1 = integ((gx(i) + vx(i) * ddt) / 
1697      &               (size1 + v1 * dtt)) + 6
1698                  i2 = integ((gy(i) + vy(i) * ddt) /
1699      &               (size2 + v2 * dtt)) + 6
1700                  i3 = integ((gz(i) + vz(i) * ddt) /
1701      &               (size3 + v3 * dtt)) + 6
1702               end if 
1703 
1704 
1705               if (k .eq. 1) i1 = 1
1706               if (k .eq. 2) i1 = 10
1707               if (k .eq. 3) i2 = 1
1708               if (k .eq. 4) i2 = 10
1709               if (k .eq. 5) i3 = 1
1710               if (k .eq. 6) i3 = 10
1711 
1712               j = icel(i1, i2, i3)
1713               call newcre(i, j)
1714               icel(i1, i2, i3) = j
1715               
1716               icels(i) = i1 * 10000 + i2 * 100 + i3
1717               
1718            end if
1719         end if
1720 
1721         if (next(i) .ne. 0) then
1722            otmp = ot(next(i))
1723            ctmp = ct(next(i))
1724         end if
1725 
1726         if (i1 .eq. 11 .and. i2 .eq. 11 .and. i3 .eq. 11) then
1727            call dchout(i, k, t)
1728         else
1729            if (iconfg .eq. 1) then
1730               call dchin1(i, k, i1, i2, i3, t)
1731            else if (iconfg .eq. 2) then
1732               call dchin2(i, k, i1, i2, i3, t)
1733            else if (iconfg .eq. 4) then
1734               call dchin3(i, k, i1, i2, i3, t)              
1735            end if
1736         end if
1737 
1738         if (icsta(i) / 10 .eq. 11) then
1739            ot(next(i)) = otmp
1740            ct(next(i)) = ctmp
1741            next(next(i)) = i
1742            call wallc(i, i1, i2, i3, t0, tmin1)
1743            if (tmin1 .lt. ct(i)) then
1744               icsta(i) = icsta(i) + 10
1745               t0 = tmin1
1746               goto 1000
1747            end if
1748         end if
1749 
1750         return
1751         end
1752            
1753         subroutine oldcre(i) 
1754 c       this subroutine is used to rearrange the old cell nic when a particle 
1755 c       goes out of the cell
1756 
1757         implicit double precision (a-h, o-z)
1758         parameter (MAXPTN=400001)
1759         common /ilist1/
1760      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1761      &     ictype, icsta(MAXPTN),
1762      &     nic(MAXPTN), icels(MAXPTN)
1763 cc      SAVE /ilist1/
1764         SAVE   
1765 
1766         if (nic(i) .eq. 0) return
1767 
1768         j = nic(i)
1769 
1770         if (nic(j) .eq. i) then
1771            nic(j) = 0
1772            return
1773         end if
1774 
1775         do 10 while (nic(j) .ne. i)
1776            j = nic(j)
1777  10        continue
1778         
1779         nic(j) = nic(i)
1780 
1781         return
1782         end
1783 
1784 
1785         subroutine newcre(i, k)
1786 c       this subroutine is used to mk rearrange of the new cell a particle
1787 c       enters,
1788 c       input i
1789 c       output nic(i)
1790 
1791         implicit double precision (a-h, o-z)
1792         parameter (MAXPTN=400001)
1793         common /ilist1/
1794      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1795      &     ictype, icsta(MAXPTN),
1796      &     nic(MAXPTN), icels(MAXPTN)
1797 cc      SAVE /ilist1/
1798         SAVE   
1799 
1800         if (k .eq. 0) then
1801            k = i
1802            nic(i) = 0
1803         else if (nic(k) .eq. 0) then
1804            nic(k) = i
1805            nic(i) = k
1806         else
1807            j = k
1808            do 10 while (nic(j) .ne. k)
1809               j = nic(j)
1810  10           continue
1811 
1812            nic(j) = i
1813            nic(i) = k
1814 
1815         end if
1816         
1817         return
1818         end
1819 
1820         subroutine scat(t, iscat, jscat)
1821 
1822 c       this subroutine is used to calculate the 2 particle scattering
1823 
1824         implicit double precision (a-h, o-z)
1825         SAVE   
1826 
1827         call newpos(t, iscat)
1828         call newpos(t, jscat)
1829         call newmom(t)
1830 
1831         return
1832         end
1833 
1834         subroutine newpos(t, i)
1835 
1836 c       this subroutine is used to calculate the 2 particle scattering
1837 c       get new position
1838 
1839         implicit double precision (a-h, o-z)
1840         parameter (MAXPTN=400001)
1841         common /para5/ iconfg, iordsc
1842 cc      SAVE /para5/
1843         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1844      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1845      &       xmass(MAXPTN), ityp(MAXPTN)
1846 cc      SAVE /prec2/
1847         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1848 cc      SAVE /prec4/
1849         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1850 cc      SAVE /prec5/
1851         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
1852 cc      SAVE /ilist5/
1853         SAVE   
1854 
1855         dt1 = ct(i) - ft(i)
1856         
1857         gx(i) = gx(i) + vx(i) * dt1
1858         gy(i) = gy(i) + vy(i) * dt1
1859         gz(i) = gz(i) + vz(i) * dt1
1860         ft(i) = ct(i)
1861            
1862         if (iconfg .le. 3) then
1863            if (ft(i) .le. abs(gz(i))) then
1864               eta(i) = 1000000.d0
1865            else
1866               eta(i) = 0.5d0 * log((ft(i) + gz(i)) / (ft(i) - gz(i)))
1867            end if
1868            tau(i) = ft(i) / cosh(eta(i))
1869         end if
1870 
1871         return
1872         end
1873 
1874         subroutine newmom(t)
1875 
1876 c       this subroutine is used to calculate the 2 particle scattering
1877 
1878         implicit double precision (a-h, o-z)
1879 
1880         parameter (hbarc = 0.197327054d0)
1881         parameter (MAXPTN=400001)
1882         parameter (pi = 3.14159265358979d0)
1883         COMMON /para1/ mul
1884 cc      SAVE /para1/
1885         common /para2/ xmp, xmu, alpha, rscut2, cutof2
1886 cc      SAVE /para2/
1887         common /para5/ iconfg, iordsc
1888 cc      SAVE /para5/
1889 ctrans
1890         common /para6/ centy
1891 cc      SAVE /para6/
1892 ctransend
1893         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
1894      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
1895      &       xmass(MAXPTN), ityp(MAXPTN)
1896 cc      SAVE /prec2/
1897         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
1898 cc      SAVE /prec4/
1899         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
1900 cc      SAVE /prec5/
1901         common /aurec1/ jxa, jya, jza
1902 cc      SAVE /aurec1/
1903         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
1904 cc      SAVE /aurec2/
1905         common /ilist1/
1906      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
1907      &     ictype, icsta(MAXPTN),
1908      &     nic(MAXPTN), icels(MAXPTN)
1909 cc      SAVE /ilist1/
1910         common /ilist3/ size1, size2, size3, v1, v2, v3, size
1911 cc      SAVE /ilist3/
1912         common /lor/ enenew, pxnew, pynew, pznew
1913 cc      SAVE /lor/
1914         common /cprod/ xn1, xn2, xn3
1915 cc      SAVE /cprod/
1916         common /rndm2/ iff
1917 cc      SAVE /rndm2/
1918         common/anim/nevent,isoft,isflag,izpc
1919 cc      SAVE /anim/
1920         common /frzprc/ 
1921      &       gxfrz(MAXPTN), gyfrz(MAXPTN), gzfrz(MAXPTN), ftfrz(MAXPTN),
1922      &       pxfrz(MAXPTN), pyfrz(MAXPTN), pzfrz(MAXPTN), efrz(MAXPTN),
1923      &       xmfrz(MAXPTN), 
1924      &       tfrz(302), ifrz(MAXPTN), idfrz(MAXPTN), itlast
1925 cc      SAVE /frzprc/
1926         SAVE   
1927 
1928 clin-6/06/02 no momentum change for partons already frozen out,
1929 c     however, spatial upgrade is needed to ensure overall system freezeout:
1930       if(isoft.eq.5) then
1931          if(ifrz(iscat).eq.1.or.ifrz(jscat).eq.1) then
1932             last(iscat) = jscat
1933             last(jscat) = iscat
1934             return
1935          endif
1936       endif
1937 clin-6/06/02-end
1938 
1939 c       iff is used to randomize the interaction to have both attractive and
1940 c        repulsive
1941 
1942         iff = - iff
1943 
1944         if (iconfg .eq. 2 .or. iconfg .eq. 4) then
1945            icels1 = icels(iscat)
1946            i1 = icels1 / 10000
1947            j1 = (icels1 - i1 * 10000) / 100
1948            icels2 = icels(jscat)
1949            i2 = icels2 / 10000
1950            j2 = (icels2 - i2 * 10000) / 100
1951            if (iconfg .eq. 4) then
1952               k1 = icels1 - i1 * 10000 - j1 * 100
1953               k2 = icels2 - i2 * 10000 - j2 * 100
1954            end if
1955         end if
1956 
1957         px1 = px(iscat)
1958         py1 = py(iscat)
1959         pz1 = pz(iscat)
1960         e1 = e(iscat)
1961         x1 = gx(iscat)
1962         y1 = gy(iscat)
1963         z1 = gz(iscat)
1964         t1 = ft(iscat)
1965         px2 = px(jscat)
1966         py2 = py(jscat)
1967         pz2 = pz(jscat)
1968         e2 = e(jscat)
1969 
1970         if (iconfg .eq. 1) then
1971            x2 = gx(jscat)
1972            y2 = gy(jscat)
1973            z2 = gz(jscat)
1974         else if (iconfg .eq. 2 .or. iconfg .eq. 4) then
1975            if (i1 - i2 .gt. 5) then
1976               x2 = gx(jscat) + 10d0 * size1
1977            else if (i1 - i2 .lt. -5) then
1978               x2 = gx(jscat) - 10d0 * size1
1979            else
1980               x2 = gx(jscat)
1981            end if
1982            if (j1 - j2 .gt. 5) then
1983               y2 = gy(jscat) + 10d0 * size2
1984            else if (j1 - j2 .lt. -5) then
1985               y2 = gy(jscat) - 10d0 * size2
1986            else
1987               y2 = gy(jscat)
1988            end if
1989            if (iconfg .eq. 4) then
1990               if (k1 - k2 .gt. 5) then
1991                  z2 = gz(jscat) + 10d0 * size3
1992               else if (k1 - k2 .lt. -5) then
1993                  z2 = gz(jscat) - 10d0 * size3
1994               else
1995                  z2 = gz(jscat)
1996               end if
1997            else
1998               z2 = gz(jscat)
1999            end if
2000         else if (iconfg .eq. 3 .or. iconfg .eq. 5) then
2001            x2 = gx(jscat) + dgxa(jscat)
2002            y2 = gy(jscat) + dgya(jscat)
2003            if (iconfg .eq. 5) then
2004               z2 = gz(jscat) + dgza(jscat)
2005            else
2006               z2 = gz(jscat)
2007            end if
2008         end if
2009         t2 = ft(jscat)
2010 ctrans
2011         rts2 = (e1 + e2) ** 2 - (px1 + px2) ** 2 -
2012      &     (py1 + py2) ** 2 - (pz1 + pz2) ** 2
2013 ctransend
2014         bex = (px1 + px2) / (e1 + e2)
2015         bey = (py1 + py2) / (e1 + e2)
2016         bez = (pz1 + pz2) / (e1 + e2)
2017 
2018         call lorenz(e1, px1, py1, pz1, bex, bey, bez)
2019 cc      SAVE pxnew, ..., values for later use.
2020         px1 = pxnew
2021         py1 = pynew
2022         pz1 = pznew
2023         e1 = enenew
2024 
2025         pp2 = pxnew ** 2 + pynew ** 2 + pznew ** 2
2026         call getht(iscat, jscat, pp2, that)
2027         theta = dacos(that / (2d0 * pp2) + 1d0)
2028         theta = dble(iff) * theta
2029 
2030 c       we boost to the cm frame, get rotation axis, and rotate 1 particle 
2031 c       momentum
2032 
2033         call lorenz(t1, x1, y1, z1, bex, bey, bez)
2034 
2035         x1 = pxnew
2036         y1 = pynew
2037         z1 = pznew
2038 
2039         call lorenz(t2, x2, y2, z2, bex, bey, bez)
2040 
2041         x2 = pxnew
2042         y2 = pynew
2043         z2 = pznew
2044 
2045 c       notice now pxnew, ..., are new positions
2046         call cropro(x1-x2, y1-y2, z1-z2, px1, py1, pz1)
2047 
2048         call xnormv(xn1, xn2, xn3)
2049 
2050 cbz1/29/99
2051 c        call rotate(xn1, xn2, xn3, theta, px1, py1, pz1)
2052         call zprota(xn1, xn2, xn3, theta, px1, py1, pz1)
2053 cbz1/29/99end
2054 
2055 c       we invert the momentum to get the other particle's momentum
2056         px2 = -px1
2057         py2 = -py1
2058         pz2 = -pz1
2059 clin-4/13/01: modify in case m1, m2 are different:
2060 c        e2 = e1
2061         e2 = dsqrt(px2**2+py2**2+pz2**2+xmass(jscat)**2)
2062 
2063 c       boost the 2 particle 4 momentum back to lab frame
2064         call lorenz(e1, px1, py1, pz1, -bex, -bey, -bez)
2065         px(iscat) = pxnew
2066         py(iscat) = pynew
2067         pz(iscat) = pznew
2068         e(iscat) = enenew
2069         call lorenz(e2, px2, py2, pz2, -bex, -bey, -bez)        
2070         px(jscat) = pxnew
2071         py(jscat) = pynew
2072         pz(jscat) = pznew
2073         e(jscat) = enenew
2074 
2075         vx(iscat) = px(iscat) / e(iscat)
2076         vy(iscat) = py(iscat) / e(iscat)
2077         vz(iscat) = pz(iscat) / e(iscat)
2078         vx(jscat) = px(jscat) / e(jscat)
2079         vy(jscat) = py(jscat) / e(jscat)
2080         vz(jscat) = pz(jscat) / e(jscat)
2081         
2082         last(iscat) = jscat
2083         last(jscat) = iscat
2084 
2085         if (iconfg .le. 3) then
2086            if (e(iscat) .le. abs(pz(iscat))) then
2087               rap(iscat) = 1000000.d0
2088            else
2089               rap(iscat) = 0.5d0 * log((e(iscat) + pz(iscat)) /
2090      &           (e(iscat) - pz(iscat)))
2091            end if
2092 
2093            if (e(jscat) .le. abs(pz(jscat))) then
2094               rap(jscat) = 1000000.d0
2095            else
2096               rap(jscat) = 0.5d0 * log((e(jscat) + pz(jscat)) /
2097      &           (e(jscat) - pz(jscat)))
2098            end if
2099 
2100 ctrans
2101            rap1 = rap(iscat)
2102            rap2 = rap(jscat)
2103 
2104            if ((rap1 .lt. centy + 0.5d0 .and.
2105      &        rap1 .gt. centy - 0.5d0)) then
2106 c              write (9, *) sqrt(ft(iscat) ** 2 - gz(iscat) ** 2), rts2
2107            end if
2108            if ((rap2 .lt. centy + 0.5d0 .and.
2109      &        rap2 .gt. centy - 0.5d0)) then
2110 c              write (9, *) sqrt(ft(jscat) ** 2 - gz(jscat) ** 2), rts2
2111            end if
2112 ctransend
2113         end if
2114 
2115         return
2116         end
2117 
2118         subroutine getht(iscat, jscat, pp2, that)
2119 
2120 c       this subroutine is used to get \hat{t} for a particular processes
2121 
2122         implicit double precision (a-h, o-z)
2123 
2124         parameter (hbarc = 0.197327054d0)
2125         parameter (MAXPTN=400001)
2126         common /para2/ xmp, xmu, alpha, rscut2, cutof2
2127 cc      SAVE /para2/
2128         common/anim/nevent,isoft,isflag,izpc
2129 cc      SAVE /anim/
2130 cc        external ran1
2131         common /rndm3/ iseedp
2132 cc      SAVE /rndm3/
2133         SAVE   
2134 
2135         iseed=iseedp
2136         xmu2 = (hbarc * xmu) ** 2
2137         xmp2 = xmp ** 2
2138         xm2 = xmu2 + xmp2
2139         rx=ran1(iseed)
2140         that = xm2*(1d0+1d0/((1d0-xm2/(4d0*pp2+xm2))*rx-1d0))
2141 ctest off isotropic scattering:
2142 c     &     + 1d0/((1d0 - xm2 / (4d0 * pp2 + xm2)) * ran1(2) - 1d0))
2143 c        if(izpc.eq.100) that=-4d0*pp2*ran1(2)
2144         if(izpc.eq.100) that=-4d0*pp2*rx
2145 
2146         return
2147         end
2148 
2149       subroutine ulist(t)
2150 c     this subroutine is used to update a new collision time list
2151 c       notice this t has been updated
2152 
2153         implicit double precision (a-h, o-z)
2154         parameter (MAXPTN=400001)
2155         common /ilist1/
2156      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
2157      &     ictype, icsta(MAXPTN),
2158      &     nic(MAXPTN), icels(MAXPTN)
2159 cc      SAVE /ilist1/
2160         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
2161 cc      SAVE /ilist4/
2162         SAVE   
2163 
2164         if (ictype .eq. 1 .or. ictype .eq. 2 .or. ictype .eq. 5
2165      &     .or. ictype .eq. 6) then
2166            l = ifmpt
2167            call ulist1(l, t)
2168         end if
2169         if (ictype .ne. 1) then
2170            l = iscat
2171            call ulist1(l, t)
2172            if (jscat .ne. 0) then
2173               l = jscat
2174               call ulist1(l, t)
2175            end if
2176         end if
2177 
2178         return
2179         end
2180 
2181         subroutine ulist1(l, t)
2182 c       this subroutine is used to update the interaction list when particle
2183 c       l is disturbed.
2184 
2185         implicit double precision (a-h, o-z)
2186         parameter (MAXPTN=400001)
2187         common /para5/ iconfg, iordsc
2188 cc      SAVE /para5/
2189         common /ilist1/
2190      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
2191      &     ictype, icsta(MAXPTN),
2192      &     nic(MAXPTN), icels(MAXPTN)
2193 cc      SAVE /ilist1/
2194         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2195 cc      SAVE /ilist5/
2196         SAVE   
2197 
2198         icels0 = icels(l)
2199         i1 = icels0 / 10000
2200         i2 = (icels0 - i1 * 10000) / 100
2201         i3 = icels0 - i1 * 10000 - i2 * 100
2202 c       save collision info for use when the collision is a collision with wall
2203 c       otherwise wallc will change icsta
2204         k = mod(icsta(l), 10)
2205 
2206         call wallc(l, i1, i2, i3, t, tmin1)
2207         tmin = tmin1
2208         nc = 0
2209 
2210         if (i1 .eq. 11 .and. i2 .eq. 11 .and. i3 .eq. 11) then
2211            call chkout(l, t, tmin, nc)
2212         else
2213            if (iconfg .eq. 1) then
2214               call chkin1(l, i1, i2, i3, t, tmin, nc)
2215            else if (iconfg .eq. 2) then
2216               call chkin2(l, i1, i2, i3, t, tmin, nc)
2217            else if (iconfg .eq. 4) then
2218               call chkin3(l, i1, i2, i3, t, tmin, nc)
2219            else if (iconfg .eq. 3 .or. iconfg .eq. 5) then
2220               call chkcel(l, i1, i2, i3, t, tmin, nc)
2221            end if
2222         end if
2223         
2224         call fixtim(l, t, tmin1, tmin, nc)
2225 
2226         return
2227         end
2228         
2229         subroutine wallc(i, i1, i2, i3, t, tmin)
2230 c       this subroutine calculates the next time for collision with wall 
2231 c       for particle i
2232 c       input particle label i,t
2233 c       output tmin collision time with wall, icsta(i) wall collision
2234 c       information
2235 
2236         implicit double precision (a-h, o-z)
2237         parameter (MAXPTN=400001)
2238         common /para5/ iconfg, iordsc
2239 cc      SAVE /para5/
2240         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2241 cc      SAVE /ilist5/
2242         SAVE   
2243 
2244         tmin = tlarge
2245 
2246         if (iconfg .le. 2 .or. iconfg .eq. 4) then
2247 c       if particle is inside the cube
2248            if ((i1 .ge. 1 .and. i1 .le. 10)
2249      &          .or. (i2 .ge. 1 .and. i2 .le. 10)
2250      &          .or. (i3 .ge. 1 .and. i3 .le. 10)) then
2251               call wallc1(i, i1, i2, i3, t, tmin)
2252 c       if particle is outside the cube
2253            else
2254               call wallcb(i, t, tmin)              
2255            end if
2256         else if (iconfg .eq. 3 .or. iconfg .eq. 5) then
2257            call wallc2(i, i1, i2, i3, t, tmin)
2258         end if
2259 
2260         return
2261         end
2262 
2263         subroutine wallc1(i, i1, i2, i3, t, tmin)
2264 c       this subroutine is used to get wall collision time
2265 c       when particle is inside the cube, it sets the icsta at the same time
2266 c       input i,i1,i2,i3,t
2267 c       output tmin, icsta(i)
2268 c       note the icsta is not finally set. we need further judgement in 
2269 c       fixtim
2270 
2271         implicit double precision (a-h, o-z)
2272         parameter (MAXPTN=400001)
2273         common /para5/ iconfg, iordsc
2274 cc      SAVE /para5/
2275         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
2276      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
2277      &       xmass(MAXPTN), ityp(MAXPTN)
2278 cc      SAVE /prec2/
2279         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
2280 cc      SAVE /prec4/
2281         common /ilist1/
2282      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
2283      &     ictype, icsta(MAXPTN),
2284      &     nic(MAXPTN), icels(MAXPTN)
2285 cc      SAVE /ilist1/
2286         common /ilist3/ size1, size2, size3, v1, v2, v3, size
2287 cc      SAVE /ilist3/
2288         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2289 cc      SAVE /ilist5/
2290         SAVE   
2291 
2292         x1p = gx(i)
2293         x2p = gy(i)
2294         x3p = gz(i)
2295         tf = ft(i)
2296         v1p = vx(i)
2297         v2p = vy(i)
2298         v3p = vz(i)
2299 
2300         if (t .lt. size .and. tf .lt. size) then
2301 
2302            if (v1p .gt. 0d0) then
2303               t1 = ((dble(i1) - 5d0) * size1 - x1p) / v1p + tf
2304            else if (v1p .lt. 0d0) then
2305               t1 = ((dble(i1) - 6d0) * size1 - x1p) / v1p + tf
2306            else
2307               t1 = tlarge
2308            end if
2309            
2310            if (v2p .gt. 0d0) then
2311               t2 = ((dble(i2) - 5d0) * size2 - x2p) / v2p + tf
2312            else if (v2p .lt. 0d0) then
2313               t2 = ((dble(i2) - 6d0) * size2 - x2p) / v2p + tf
2314            else
2315               t2 = tlarge
2316            end if
2317            
2318            if (v3p .gt. 0d0) then
2319               t3 = ((dble(i3) - 5d0) * size3 - x3p) / v3p + tf
2320            else if (v3p .lt. 0d0) then
2321               t3 = ((dble(i3) - 6d0) * size3 - x3p) / v3p + tf
2322            else
2323               t3 = tlarge
2324            end if
2325            
2326 c       if a particle is on the wall, we don't collide it on the same wall
2327            
2328 c        if (t1 .eq. 0d0) t1 = tlarge
2329 c        if (t2 .eq. 0d0) t2 = tlarge
2330 c        if (t3 .eq. 0d0) t3 = tlarge
2331            
2332            tmin = min(t1, t2, t3)
2333            
2334 c       set icsta,
2335 c       after checking this is not an earlier collision comparing with 
2336 c       a collision with another particle, we need to set icsta=0
2337 c       after checking whether there is also a particle collision 
2338 c       at the same time, we need to reset the second bit of icsta
2339            
2340            if (tmin .eq. t1) then
2341               if (v1p .gt. 0d0) then
2342                  icsta(i) = 101
2343               else
2344                  icsta(i) = 102
2345               end if
2346            end if
2347            
2348            if (tmin .eq. t2) then
2349               if (v2p .gt. 0d0) then
2350                  icsta(i) = 103
2351               else
2352                  icsta(i) = 104
2353               end if
2354            end if
2355            
2356            if (tmin .eq. t3) then
2357               if (v3p .gt. 0d0) then
2358                  icsta(i) = 105
2359               else
2360                  icsta(i) = 106
2361               end if
2362            end if
2363            
2364         if (tmin .le. size) return
2365 
2366         end if
2367 
2368         if (v1p .gt. (i1 - 5) * v1) then
2369            t1 = ((i1 - 5) * (size1 - v1 * size) +
2370      &          v1p * tf - x1p) / (v1p - (i1 - 5) * v1)
2371         else if (v1p .lt. (i1 - 6) * v1) then
2372            t1 = ((i1 - 6) * (size1 - v1 * size) +
2373      &          v1p * tf - x1p) / (v1p - (i1 - 6) * v1)
2374         else
2375            t1 = tlarge
2376         end if
2377         
2378         if (v2p .gt. (i2 - 5) * v2) then
2379            t2 = ((i2 - 5) * (size2 - v2 * size) +
2380      &          v2p * tf - x2p) / (v2p - (i2 - 5) * v2)
2381         else if (v2p .lt. (i2 - 6) * v2) then
2382            t2 = ((i2 - 6) * (size2 - v2 * size) +
2383      &          v2p * tf - x2p) / (v2p - (i2 - 6) * v2)
2384         else
2385            t2 = tlarge
2386         end if
2387         
2388         if (v3p .gt. (i3 - 5) * v3) then
2389            t3 = ((i3 - 5) * (size3 - v3 * size) +
2390      &          v3p * tf - x3p) / (v3p - (i3 - 5) * v3)
2391         else if (v3p .lt. (i3 - 6) * v3) then
2392            t3 = ((i3 - 6) * (size3 - v3 * size) +
2393      &          v3p * tf - x3p) / (v3p - (i3 - 6) * v3)
2394         else
2395            t3 = tlarge
2396         end if
2397         
2398 c       if a particle is on the wall, we don't collide it on the same wall
2399         
2400 c        if (t1 .eq. 0d0) t1 = tlarge
2401 c        if (t2 .eq. 0d0) t2 = tlarge
2402 c        if (t3 .eq. 0d0) t3 = tlarge
2403         
2404         tmin = min(t1, t2, t3)
2405         
2406 c       set icsta,
2407 c       after checking this is not an earlier collision comparing with 
2408 c       a collision with another particle, we need to set icsta=0
2409 c       after checking whether there is also a particle collision 
2410 c       at the same time, we need to reset the second bit of icsta
2411         
2412         if (tmin .eq. t1) then
2413            if (v1p .gt. (i1 - 5) * v1) then
2414               icsta(i) = 101
2415            else
2416               icsta(i) = 102
2417            end if
2418         end if
2419         
2420         if (tmin .eq. t2) then
2421            if (v2p .gt. (i2 - 5) * v2) then
2422               icsta(i) = 103
2423            else
2424               icsta(i) = 104
2425            end if
2426         end if
2427         
2428         if (tmin .eq. t3) then
2429            if (v3p .gt. (i3 - 5) * v3) then
2430               icsta(i) = 105
2431            else
2432               icsta(i) = 106
2433            end if
2434         end if
2435         
2436         return
2437         end
2438 
2439         subroutine wallc2(i, i1, i2, i3, t, tmin)
2440 c       this subroutine is used to get wall collision time
2441 c       when particle is inside the cube, it sets the icsta at the same time
2442 c       input i,i1,i2,i3,t
2443 c       output tmin, icsta(i)
2444 c       note the icsta is not finally set. we need further judgement in 
2445 c       fixtim
2446 
2447         implicit double precision (a-h, o-z)
2448         parameter (MAXPTN=400001)
2449 
2450         common /para5/ iconfg, iordsc
2451 cc      SAVE /para5/
2452         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
2453      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
2454      &       xmass(MAXPTN), ityp(MAXPTN)
2455 cc      SAVE /prec2/
2456         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
2457 cc      SAVE /prec4/
2458         common /ilist1/
2459      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
2460      &     ictype, icsta(MAXPTN),
2461      &     nic(MAXPTN), icels(MAXPTN)
2462 cc      SAVE /ilist1/
2463         common /ilist3/ size1, size2, size3, v1, v2, v3, size
2464 cc      SAVE /ilist3/
2465         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2466 cc      SAVE /ilist5/
2467         SAVE   
2468 
2469         x1p = gx(i)
2470         x2p = gy(i)
2471         x3p = gz(i)
2472         tf = ft(i)
2473         v1p = vx(i)
2474         v2p = vy(i)
2475         v3p = vz(i)
2476 
2477         if (v1p .gt. 0d0) then
2478            t1 = (5d0 * size1 - x1p) / v1p + tf
2479         else if (v1p .lt. 0d0) then
2480            t1 = (-5d0 * size1 - x1p) / v1p + tf
2481         else
2482            t1 = tlarge
2483         end if
2484         
2485         if (v2p .gt. 0d0) then
2486            t2 = (5d0 * size2 - x2p) / v2p + tf
2487         else if (v2p .lt. 0d0) then
2488            t2 = (- 5d0 * size2 - x2p) / v2p +tf
2489         else
2490            t2 = tlarge
2491         end if
2492 
2493         if (iconfg .eq. 5) then
2494            if (v3p .gt. 0d0) then
2495               t3 = (5d0 * size3 - x3p) / v3p + tf
2496            else if (v3p .lt. 0d0) then
2497               t3 = (- 5d0 * size3 - x3p) / v3p +tf
2498            else
2499               t3 = tlarge
2500            end if
2501         else
2502            t3 = tlarge
2503         end if
2504            
2505         tmin = min(t1, t2, t3)
2506         
2507 c       set icsta,
2508 c       after checking this is not an earlier collision comparing with 
2509 c       a collision with another particle, we need to set icsta=0
2510 c       after checking whether there is also a particle collision 
2511 c       at the same time, we need to reset the second bit of icsta
2512            
2513         if (tmin .eq. t1) then
2514            if (v1p .gt. 0d0) then
2515               icsta(i) = 101
2516            else
2517               icsta(i) = 102
2518            end if
2519         end if
2520         
2521         if (tmin .eq. t2) then
2522            if (v2p .gt. 0d0) then
2523               icsta(i) = 103
2524            else
2525               icsta(i) = 104
2526            end if
2527         end if
2528 
2529         if (tmin .eq. t3) then
2530            if (v3p .gt. 0d0) then
2531               icsta(i) = 105
2532            else
2533               icsta(i) = 106
2534            end if
2535         end if
2536            
2537         return
2538         end
2539 
2540         subroutine wallcb(i, t, tmin)
2541 c       this subroutine is used to calculate the wall collision time 
2542 c       when the particle is outside the cube
2543 c       input i,t
2544 c       output tmin,icsta(i)
2545 c       note the icsta is not finally set. we need further judgement in 
2546 c       fixtim
2547 
2548         implicit double precision (a-h, o-z)
2549         parameter (MAXPTN=400001)
2550 
2551         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
2552      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
2553      &       xmass(MAXPTN), ityp(MAXPTN)
2554 cc      SAVE /prec2/
2555         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
2556 cc      SAVE /prec4/
2557         common /ilist1/
2558      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
2559      &     ictype, icsta(MAXPTN),
2560      &     nic(MAXPTN), icels(MAXPTN)
2561 cc      SAVE /ilist1/
2562         common /ilist3/ size1, size2, size3, v1, v2, v3, size
2563 cc      SAVE /ilist3/
2564         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
2565 cc      SAVE /ilist5/
2566         SAVE   
2567 
2568 c       check if there is a collision by looking at the closest approach point
2569 c       and see if it's inside the cube
2570 
2571         if (size1 .eq. 0d0 .or. size2 .eq. 0d0 .or. 
2572      &     size3 .eq. 0d0) return
2573 
2574         x1p = gx(i)
2575         x2p = gy(i)
2576         x3p = gz(i)
2577         v1p = vx(i)
2578         v2p = vy(i)
2579         v3p = vz(i)
2580         tf = ft(i)
2581 
2582         if (t .lt. size .and. tf .lt. size) then
2583            if (x1p .lt. - 5d0 * size1 .and. v1p .gt. 0d0) then
2584               t1 = (- 5d0 * size1 - x1p) / v1p + tf
2585            else if(x1p .gt. 5d0 * size1 .and. v1p .lt. 0d0) then
2586               t1 = - (x1p - 5d0 * size1) / v1p + tf
2587            else
2588               t1 = tlarge 
2589            end if
2590 
2591            if (t1 .ne. tlarge) then
2592               x2pp = x2p + v2p * (t1 - tf)
2593               x3pp = x3p + v3p * (t1 - tf)
2594               if (x2pp .le. - 5d0 * size2 .or. x2pp .ge. 5d0 * size2
2595      &             .or. x3pp .le. - 5d0 * size3 
2596      &             .or. x3pp .ge. 5d0 * size3)
2597      &             t1 = tlarge
2598            end if
2599            
2600            if (x2p .lt. - 5d0 * size2 .and. v2p .gt. 0d0) then
2601               t2 = (- 5d0 * size2 - x2p) / v2p + tf
2602            else if(x2p .gt. 5d0 * size2 .and. v2p .lt. 0d0) then
2603               t2 = - (x2p - 5d0 * size2) / v2p + tf
2604            else
2605               t2 = tlarge 
2606            end if
2607            
2608            if (t2 .ne. tlarge) then
2609               x1pp = x1p + v1p * (t2 - tf)
2610               x3pp = x3p + v3p * (t2 - tf)
2611               if (x1pp .le. - 5d0 * size1 .or. x1pp .ge. 5d0 * size1
2612      &          .or. x3pp .le. - 5d0 * size3 .or. x3pp .ge. 5d0 * size3)
2613      &             t2 = tlarge
2614            end if
2615            
2616            if (x3p .lt. - 5d0 * size3 .and. v3p .gt. 0d0) then
2617               t3 = (- 5d0 * size3 - x3p) / v3p + tf
2618            else if(x3p .gt. 5d0 * size3 .and. v3p .lt. 0d0) then
2619               t3 = - (x3p - 5d0 * size3) / v3p + tf
2620            else
2621               t3 = tlarge 
2622            end if
2623            
2624            if (t3 .ne. tlarge) then
2625               x1pp = x1p + v1p * (t3 - tf)
2626               x2pp = x2p + v2p * (t3 - tf)
2627               if (x1pp .le. - 5d0 * size1 .or. x1pp .ge. 5d0 * size1
2628      &          .or. x2pp .le. - 5d0 * size2 .or. x2pp .ge. 5d0 * size2)
2629      &             t3 = tlarge
2630            end if
2631            
2632            tmin = min(t1, t2, t3)
2633 
2634 c       set icsta,
2635 c       after checking this is not an earlier collision comparing with 
2636 c       a collision with another particle, we need to set icsta=0
2637 c       after checking whether there is also a particle collision 
2638 c       at the same time, we need to reset the second bit of icsta
2639 
2640            if (tmin .eq. t1) then
2641               if (v1p .gt. 0d0) then
2642                  icsta(i) = 101
2643               else
2644                  icsta(i) = 102
2645               end if
2646            end if
2647            
2648            if (tmin .eq. t2) then
2649               if (v2p .gt. 0d0) then
2650                  icsta(i) = 103
2651               else
2652                  icsta(i) = 104
2653               end if
2654            end if
2655         
2656            if (tmin .eq. t3) then
2657               if (v3p .gt. 0d0) then
2658                  icsta(i) = 105
2659               else
2660                  icsta(i) = 106
2661               end if
2662            end if
2663            
2664         if (tmin .le. size) return
2665 
2666         end if
2667 
2668 c       notice now x1q, x2q, x3q are coordinates at time t
2669         x1q = x1p + v1p * (t - tf)
2670         x2q = x2p + v2p * (t - tf)
2671         x3q = x3p + v3p * (t - tf)
2672 
2673         if (x1q .lt. - 5d0 * (size1 + v1 * (t - size)) .and. 
2674      &      v1p .gt. - 5d0 * v1) then
2675            t1 = (- 5d0 * (size1 - v1 * size) + v1p * tf - x1p) /
2676      &          (v1p - (- 5d0) * v1)
2677            icsta1 = 101
2678         else if (x1q .gt. 5d0 * (size1 + v1 * (t-size)) .and. 
2679      &     v1p .lt. 5d0 * v1) then
2680            t1 = (5d0 * (size1 - v1 * size) + v1p * tf - x1p) /
2681      &          (v1p - 5d0 * v1)
2682            icsta1 = 102
2683         else
2684            t1 = tlarge 
2685         end if
2686         
2687         if (t1 .ne. tlarge) then
2688            x2pp = x2p + v2p * (t1 - tf)
2689            x3pp = x3p + v3p * (t1 - tf)
2690            if (x2pp .le. - 5d0 * (size2 + v2 * (t1 - size))
2691      &        .or. x2pp .ge. 5d0 * (size2 + v2 * (t1 - size))
2692      &        .or. x3pp .le. - 5d0 * (size3 + v3 * (t1 - size))
2693      &        .or. x3pp .ge. 5d0 * (size3 + v3 * (t1 - size)))
2694      &        t1 = tlarge
2695         end if
2696 
2697         if (x2q .lt. - 5d0 * (size2 + v2 * (t - size)) .and.
2698      &     v2p .gt. - 5d0 * v2) then
2699            t2 = (- 5d0 * (size2 - v2 * size) + v2p * tf - x2p) /
2700      &          (v2p - (- 5d0) * v2)
2701            icsta2 = 103
2702         else if (x2q .gt. 5d0 * (size2 + v2 * (t - size)) .and.
2703      &     v2p .lt. 5d0 * v2) then
2704            t2 = (5d0 * (size2 - v2 * size) + v2p * tf - x2p) / 
2705      &          (v2p - 5d0 * v2)
2706            icsta2 = 104
2707         else
2708            t2 = tlarge 
2709         end if
2710         
2711         if (t2 .ne. tlarge) then
2712            x1pp = x1p + v1p * (t2 - tf)
2713            x3pp = x3p + v3p * (t2 - tf)
2714            if (x1pp .le. - 5d0 * (size1 + v1 * (t2 - size))
2715      &        .or. x1pp .ge. 5d0 * (size1 + v1 * (t2 - size))
2716      &        .or. x3pp .le. - 5d0 * (size3 + v3 * (t2 - size))
2717      &        .or. x3pp .ge. 5d0 * (size3 + v3 * (t2 - size)))
2718      &        t2 = tlarge
2719         end if
2720 
2721         if (x3q .lt. - 5d0 * (size3 + v3 * (t - size)) .and. 
2722      &     v3p .gt. - 5d0 * v3) then
2723            t3 = (- 5d0 * (size3 - v3 * size) + v3p * tf - x3p) /
2724      &          (v3p - (- 5d0) * v3)
2725            icsta3 = 105
2726         else if (x3q .gt. 5d0 * (size3 + v3 * (t - size)) .and.
2727      &     v3p .lt. 5d0 * v3) then
2728            t3 = (5d0 * (size3 - v3 * size) + v3p * tf - x3p) /
2729      &          (v3p - 5d0 * v3)
2730            icsta3 = 106
2731         else
2732            t3 = tlarge 
2733         end if
2734         
2735         if (t3 .ne. tlarge) then
2736            x2pp = x2p + v2p * (t3 - tf)
2737            x1pp = x1p + v1p * (t3 - tf)
2738            if (x2pp .le. - 5d0 * (size2 + v2 * (t3 - size))
2739      &        .or. x2pp .ge. 5d0 * (size2 + v2 * (t3 - size))
2740      &        .or. x1pp .le. - 5d0 * (size1 + v1 * (t3 - size))
2741      &        .or. x1pp .ge. 5d0 * (size1 + v1 * (t3 - size)))
2742      &        t3 = tlarge
2743         end if
2744         
2745         tmin = min(t1, t2, t3)
2746         
2747 c       set icsta,
2748 c       after checking this is not an earlier collision comparing with 
2749 c       a collision with another particle, we need to set icsta=0
2750 c       after checking whether there is also a particle collision 
2751 c       at the same time, we need to reset the second bit of icsta
2752         
2753         if (tmin .eq. t1) then
2754            icsta(i) = icsta1
2755         else if (tmin .eq. t2) then
2756            icsta(i) = icsta2
2757         else if (tmin .eq. t3) then
2758            icsta(i) = icsta3
2759         end if
2760         
2761         return
2762         end
2763            
2764         subroutine chkout(l, t, tmin, nc)
2765 c       this subroutine is used to check the collisions with particles in 
2766 c       surface cells to see if we can get a smaller collision time than tmin
2767 c       with particle nc, when the colliding particle is outside the cube
2768 c       input l,t,tmin,nc
2769 c       output tmin, nc
2770 
2771         implicit double precision (a-h, o-z)
2772         parameter (MAXPTN=400001)
2773 
2774         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
2775      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
2776      &       xmass(MAXPTN), ityp(MAXPTN)
2777 cc      SAVE /prec2/
2778         SAVE   
2779 
2780         m1 = 11
2781         m2 = 11
2782         m3 = 11
2783         call chkcel(l, m1, m2, m3, t, tmin, nc)
2784 
2785         do 1003 i = 1, 10
2786            do 1002 j = 1, 10
2787               do 1001 k = 1, 10
2788                  if (i .eq. 1 .or. i .eq. 10 .or. j .eq. 1
2789      &              .or. j .eq. 10 .or. k .eq. 1 .or. k .eq. 10) 
2790      &                    call chkcel(l, i, j, k, t, tmin, nc)
2791  1001         continue
2792  1002      continue
2793  1003   continue
2794 
2795         return
2796         end
2797 
2798         subroutine chkin1(l, i1, i2, i3, t, tmin, nc)
2799 c       this subroutine is used to check collisions for particle inside
2800 c       the cube
2801 c       and update the afftected particles through chkcel
2802 
2803         implicit double precision (a-h, o-z)
2804         SAVE   
2805 
2806 c       itest is a flag to make sure the 111111 cell is checked only once
2807         itest = 0
2808         
2809         do 1003 i = i1 - 1, i1 + 1
2810            do 1002 j = i2 - 1, i2 + 1
2811               do 1001 k =  i3 - 1, i3 + 1
2812                  if (i .ge. 1 .and. i .le. 10 .and. j .ge. 1 .and.
2813      &               j .le. 10 .and. k .ge. 1 .and. k .le. 10) then
2814                     call chkcel(l, i, j, k, t, tmin, nc)
2815                  else if (itest .eq. 0) then
2816                     m1 = 11
2817                     m2 = 11
2818                     m3 = 11
2819                     call chkcel(l, m1, m2, m3, t, tmin, nc)
2820                     itest = 1
2821                  end if   
2822  1001         continue
2823  1002      continue
2824  1003   continue
2825 
2826         return
2827         end
2828 
2829         subroutine chkin2(l, i1, i2, i3, t, tmin, nc)
2830 c       this subroutine is used to check collisions for particle inside
2831 c       the cube
2832 c       and update the afftected particles through chkcel
2833 
2834         implicit double precision (a-h, o-z)
2835         SAVE   
2836 
2837 c       itest is a flag to make sure the 111111 cell is checked only once
2838         itest = 0
2839         
2840         do 1003 i = i1 - 1, i1 + 1
2841            do 1002 j = i2 - 1, i2 + 1
2842               do 1001 k =  i3 - 1, i3 + 1
2843                  ia = i
2844                  ib = j
2845                  ic = k
2846                  if (k .ge. 1 .and. k .le. 10) then
2847                     if (i .eq. 0) ia = 10
2848                     if (i .eq. 11) ia = 1
2849                     if (j .eq. 0) ib = 10
2850                     if (j .eq. 11) ib = 1
2851                     call chkcel(l, ia, ib, ic, t, tmin, nc)
2852                  end if
2853  1001         continue
2854  1002      continue
2855  1003   continue
2856 
2857         return
2858         end
2859 
2860         subroutine chkin3(l, i1, i2, i3, t, tmin, nc)
2861 c       this subroutine is used to check collisions for particle inside
2862 c       the cube
2863 c       and update the afftected particles through chkcel
2864 
2865         implicit double precision (a-h, o-z)
2866         SAVE   
2867 
2868 c       itest is a flag to make sure the 111111 cell is checked only once
2869         itest = 0
2870         
2871         do 1003 i = i1 - 1, i1 + 1
2872            do 1002 j = i2 - 1, i2 + 1
2873               do 1001 k =  i3 - 1, i3 + 1
2874                  if (i .eq. 0) then
2875                     ia = 10
2876                  else if (i .eq. 11) then
2877                     ia = 1
2878                  else
2879                     ia = i
2880                  end if
2881                  if (j .eq. 0) then
2882                     ib = 10
2883                  else if (j .eq. 11) then
2884                     ib = 1
2885                  else
2886                     ib = j
2887                  end if
2888                  if (k .eq. 0) then
2889                     ic = 10
2890                  else if (k .eq. 11) then
2891                     ic = 1
2892                  else
2893                     ic = k
2894                  end if
2895                  call chkcel(l, ia, ib, ic, t, tmin, nc)
2896  1001         continue
2897  1002      continue
2898  1003   continue
2899 
2900         return
2901         end
2902 
2903         subroutine chkcel(il, i1, i2, i3, t, tmin, nc)
2904 c       this program is used to check through all the particles
2905 c       in the cell (i1,i2,i3) and see if we can get a particle collision 
2906 c       with time less than the original input tmin ( the collision time of 
2907 c       il with the wall
2908 c       and update the affected particles
2909 
2910         implicit double precision (a-h, o-z)
2911         parameter (MAXPTN=400001)
2912 
2913         common /para5/ iconfg, iordsc
2914 cc      SAVE /para5/
2915         common /ilist1/
2916      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
2917      &     ictype, icsta(MAXPTN),
2918      &     nic(MAXPTN), icels(MAXPTN)
2919 cc      SAVE /ilist1/
2920         common /ilist2/ icell, icel(10, 10, 10)
2921 cc      SAVE /ilist2/
2922         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
2923 cc      SAVE /ilist4/
2924         SAVE   
2925 
2926         if (iconfg .eq. 3 .or. iconfg .eq. 5) then
2927            jj = ichkpt
2928            do 1001 j = 1, jj
2929               call ck(j, ick)
2930 c     10/24/02 get rid of argument usage mismatch in ud2():
2931                             jud2=j
2932 c              if (ick .eq. 1) call ud2(j, il, t, tmin, nc)
2933               if (ick .eq. 1) call ud2(jud2, il, t, tmin, nc)
2934  1001      continue
2935            return
2936         end if
2937 
2938         if (i1 .eq. 11 .and. i2 .eq. 11 .and. i3 .eq. 11) then
2939            l = icell
2940         else
2941            l = icel(i1, i2, i3)
2942         end if
2943 
2944 c       if there is no particle
2945         if (l .eq. 0) then
2946            return
2947         end if
2948         j = nic(l)
2949 c       if there is only one particle
2950         if (j .eq. 0) then
2951            call ck(l, ick)
2952            if (ick .eq. 1) call ud2(l, il, t, tmin, nc)
2953 
2954 c       if there are many particles
2955         else
2956 
2957 c       we don't worry about the other colliding particle because it's
2958 c       set in last(), and will be checked in ud2
2959 
2960            call ck(l, ick)
2961            if (ick .eq. 1) call ud2(l, il, t, tmin, nc)
2962 
2963            do 10 while(j .ne. l)
2964               call ck(j, ick)
2965               if (ick .eq. 1) call ud2(j, il, t, tmin, nc)
2966               j = nic(j)
2967  10           continue
2968         end if
2969 
2970         return
2971         end
2972 
2973         subroutine ck(l, ick)
2974 c       this subroutine is used for chcell to check whether l should be
2975 c       checked or not for updating tmin, nc
2976 c       input l
2977 c       output ick
2978 c       if ick=1, l should be checked, otherwise it should not be.
2979 
2980         implicit double precision (a-h, o-z)
2981         parameter (MAXPTN=400001)
2982 
2983         common /ilist1/
2984      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
2985      &     ictype, icsta(MAXPTN),
2986      &     nic(MAXPTN), icels(MAXPTN)
2987 cc      SAVE /ilist1/
2988         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
2989 cc      SAVE /ilist4/
2990         SAVE   
2991 
2992         ick = 1
2993         if (ictype .eq. 1) then
2994            if (l .eq. ifmpt) ick = 0
2995         else if (ictype .eq. 0 .or. ictype .eq. 3 .or. 
2996      &     ictype .eq. 4) then
2997            if (l .eq. iscat .or. l .eq. jscat) ick = 0
2998         else
2999            if (l .eq. iscat .or. l .eq. jscat .or.
3000      &         l .eq. ifmpt) ick = 0
3001         end if
3002 c       notice il is either iscat or jscat, or ifmpt, we deal with them
3003 c       seperately according to ictype
3004 
3005         return
3006         end
3007            
3008         subroutine dchout(l, ii, t)
3009 c       this subroutine is used to check collisions of l with particles when 
3010 c       l is outside the cube and the collision just happened is a collision
3011 c       including a collision with wall (hence we need to use dcheck to throw
3012 c       away old collisions that are not in the new neighboring cells.
3013 
3014 c       input l,t
3015 
3016         implicit double precision (a-h, o-z)
3017         parameter (MAXPTN=400001)
3018 
3019         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
3020      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
3021      &       xmass(MAXPTN), ityp(MAXPTN)
3022 cc      SAVE /prec2/
3023         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
3024 cc      SAVE /prec4/
3025         common /ilist3/ size1, size2, size3, v1, v2, v3, size
3026 cc      SAVE /ilist3/
3027         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3028 cc      SAVE /ilist5/
3029         SAVE   
3030 
3031         external integ
3032 
3033         tt = ft(l)
3034         td = t - size
3035         x1 = gx(l) + vx(l) * (t - tt)
3036         x2 = gy(l) + vy(l) * (t - tt)
3037         x3 = gz(l) + vz(l) * (t - tt)
3038         if (td .le. 0d0) then
3039            i1 = integ(x1 / size1) + 6
3040            i2 = integ(x2 / size2 ) + 6
3041            i3 = integ(x3 / size3 ) + 6
3042            if (integ(x1 / size1) .eq. x1 / size1 .and. vx(l) .lt. 0d0)
3043      &        i1 = i1 - 1
3044            if (integ(x2 / size2) .eq. x2 / size2 .and. vy(l) .lt. 0d0)
3045      &        i2 = i2 - 1
3046            if (integ(x3 / size3) .eq. x3 / size3 .and. vz(l) .lt. 0d0)
3047      &        i3 = i3 - 1
3048         else
3049            i1 = integ(x1 / (size1 + v1 * td)) + 6
3050            i2 = integ(x2 / (size2 + v2 * td)) + 6
3051            i3 = integ(x3 / (size3 + v3 * td)) + 6
3052 c     10/24/02 (i) below should be (l):
3053            if (integ(x1 / (size1 + v1 * td)) .eq. 
3054      &        x1 / (size1 +v1 * td) .and. 
3055      &        vx(l) .lt. (i1 - 6) * v1) i1 = i1 - 1
3056 c     &        vx(i) .lt. (i1 - 6) * v1) i1 = i1 - 1
3057            if (integ(x2 / (size2 + v2 * td)) .eq.
3058      &        x2 / (size2 + v2 * td) .and.
3059      &        vy(l) .lt. (i2 - 6) * v2) i2 = i2 - 1
3060 c     &        vy(i) .lt. (i2 - 6) * v2) i2 = i2 - 1
3061            if (integ(x3 / (size3 + v3 * td)) .eq. 
3062      &        x3 / (size3 + v3 * td) .and.
3063      &        vz(l) .lt. (i3 - 6) * v3) i3 = i3 - 1
3064 c     &        vz(i) .lt. (i3 - 6) * v3) i3 = i3 - 1
3065         end if
3066 
3067         if (ii .eq. 1) then
3068            i = 9
3069            do 1002 j = i2 - 1, i2 + 1
3070               do 1001 k = i3 - 1, i3 + 1
3071                  if (j .ge. 1 .and. j .le. 10 .and. k .ge. 1 .and.
3072      &              k .le. 10) then
3073                     call dchcel(l, i, j, k, t)
3074                  end if
3075  1001         continue
3076  1002      continue
3077         end if
3078 
3079         if (ii .eq. 2) then
3080            i = 2
3081            do 1004 j = i2 - 1, i2 + 1
3082               do 1003 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  1003         continue
3088  1004      continue
3089         end if
3090 
3091         if (ii .eq. 3) then
3092            j = 9
3093            do 1006 i = i1 - 1, i1 + 1
3094               do 1005 k = i3 - 1, i3 + 1
3095                  if (i .ge. 1 .and. i .le. 10 .and. k .ge. 1 .and.
3096      &              k .le. 10) then
3097                     call dchcel(l, i, j, k, t)
3098                  end if
3099  1005         continue
3100  1006      continue
3101         end if
3102 
3103         if (ii .eq. 4) then
3104            j = 2
3105            do 1008 i = i1 - 1, i1 + 1
3106               do 1007 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  1007         continue
3112  1008      continue
3113         end if
3114 
3115         if (ii .eq. 5) then
3116            k = 9
3117            do 1010 i = i1 - 1, i1 + 1
3118               do 1009 j = i2 - 1, i2 + 1
3119                  if (i .ge. 1 .and. i .le. 10 .and. j .ge. 1 .and.
3120      &              j .le. 10) then
3121                     call dchcel(l, i, j, k, t)
3122                  end if
3123  1009         continue
3124  1010      continue
3125         end if
3126 
3127         if (ii .eq. 6) then
3128            k = 2
3129            do 1012 i = i1 - 1, i1 + 1
3130               do 1011 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  1011         continue
3136  1012      continue
3137         end if
3138 
3139         return
3140         end
3141 
3142         subroutine dchin1(l, ii, i1, i2, i3, t)
3143 c       this subroutine is used to check collisions for particle inside
3144 c       the cube when the collision just happened is a collision including 
3145 c       collision with wall
3146 c       and update the afftected particles through chkcel
3147 
3148 c       input l,ii(specifying the direction of the wall collision),
3149 c          i1,i2,i3, (specifying the position of the cell 
3150 c                    we are going to check)
3151 c          t
3152 
3153         implicit double precision (a-h, o-z)
3154         SAVE   
3155 
3156 c       itest is a flag to make sure the 111111 cell is checked only once
3157         itest = 0
3158         
3159         if (ii .eq. 1) then
3160            if (i1 .eq. 1) goto 100
3161            if (i1 .eq. 2) then
3162               if (i2 .ge. 2 .and. i2 .le. 9 .and. i3 .ge. 2 .and.
3163      &           i3 .le. 9) then
3164                  i = 11
3165                  j = 11
3166                  k = 11
3167                  call dchcel(l, i, j, k, t)
3168               end if
3169               goto 100
3170            end if
3171            i = i1 - 2
3172            do 1002 j = i2 - 1, i2 + 1
3173               do 1001 k = i3 - 1, i3 + 1
3174                  if (j .ge. 1 .and. j .le. 10 .and. k .ge. 1 .and.
3175      &              k .le. 10)
3176      &                    call dchcel(l, i, j, k, t)
3177  1001         continue
3178  1002      continue
3179         end if
3180 
3181         if (ii .eq. 2) then
3182            if (i1 .eq. 10) goto 100
3183            if (i1 .eq. 9) then
3184               if (i2 .ge. 2 .and. i2 .le. 9 .and. i3 .ge. 2 .and.
3185      &           i3 .le. 9) then
3186                  i = 11
3187                  j = 11
3188                  k = 11
3189                  call dchcel(l, i, j, k, t)
3190               end if
3191               goto 100
3192            end if
3193            i = i1 + 2
3194            do 1004 j = i2 - 1, i2 + 1
3195               do 1003 k = i3 - 1, i3 + 1
3196                  if (j .ge. 1 .and. j .le. 10 .and. k .ge. 1 .and.
3197      &              k .le. 10)
3198      &                    call dchcel(l, i, j, k, t)
3199  1003         continue
3200  1004      continue
3201         end if
3202 
3203         if (ii .eq. 3) then
3204            if (i2 .eq. 1) goto 100
3205            if (i2 .eq. 2) then
3206               if (i1 .ge. 2 .and. i1 .le. 9 .and. i3 .ge. 2 .and.
3207      &           i3 .le. 9) then
3208                  i = 11
3209                  j = 11
3210                  k = 11
3211                  call dchcel(l, i, j, k, t)
3212               end if
3213               goto 100
3214            end if
3215            j = i2 - 2
3216            do 1006 i = i1 - 1, i1 + 1
3217               do 1005 k = i3 - 1, i3 + 1
3218                  if (i .ge. 1 .and. i .le. 10 .and. k .ge. 1 .and.
3219      &              k .le. 10)
3220      &              call dchcel(l, i, j, k, t)
3221  1005         continue
3222  1006      continue
3223         end if
3224 
3225         if (ii .eq. 4) then
3226            if (i2 .eq. 10) goto 100
3227            if (i2 .eq. 9) then
3228               if (i1 .ge. 2 .and. i1 .le. 9 .and. i3 .ge. 2 .and.
3229      &           i3 .le. 9) then
3230                  i = 11
3231                  j = 11
3232                  k = 11
3233                  call dchcel(l, i, j, k, t)
3234               end if
3235               goto 100
3236            end if
3237            j = i2 + 2
3238            do 1008 i = i1 - 1, i1 + 1
3239               do 1007 k = i3 - 1, i3 + 1
3240                  if (i .ge. 1 .and. i .le. 10 .and. k .ge. 1 .and.
3241      &           k .le. 10)
3242      &                 call dchcel(l, i, j, k, t)
3243  1007         continue
3244  1008      continue
3245         end if
3246 
3247         if (ii .eq. 5) then
3248            if (i3 .eq. 1) goto 100
3249            if (i3 .eq. 2) then
3250               if (i1 .ge. 2 .and. i1 .le. 9 .and. i2 .ge. 2 .and.
3251      &           i2 .le. 9) then
3252                  i = 11
3253                  j = 11
3254                  k = 11
3255                  call dchcel(l, i, j, k, t)
3256               end if
3257               goto 100
3258            end if
3259            k = i3 - 2
3260            do 1010 i = i1 - 1, i1 + 1
3261               do 1009 j = i2 - 1, i2 + 1
3262                  if (i .ge. 1 .and. i .le. 10 .and. j .ge. 1 .and.
3263      &           j .le. 10)
3264      &                 call dchcel(l, i, j, k, t)
3265  1009         continue
3266  1010      continue
3267         end if
3268 
3269         if (ii .eq. 6) then
3270            if (i3 .eq. 10) goto 100
3271            if (i3 .eq. 9) then
3272               if (i1 .ge. 2 .and. i1 .le. 9 .and. i2 .ge. 2 .and.
3273      &           i2 .le. 9) then
3274                  i = 11
3275                  j = 11
3276                  k = 11
3277                  call dchcel(l, i, j, k, t)
3278               end if
3279               goto 100
3280            end if
3281            k = i3 + 2
3282            do 1012 i = i1 - 1, i1 + 1
3283               do 1011 j = i2 - 1, i2 + 1
3284                  if (i .ge. 1 .and. i .le. 10 .and. j .ge. 1 .and.
3285      &           j .le. 10)
3286      &                 call dchcel(l, i, j, k, t)
3287  1011         continue
3288  1012      continue
3289         end if
3290 
3291  100        continue
3292 
3293         return
3294         end
3295 
3296         subroutine dchin2(l, ii, i1, i2, i3, t)
3297 c       this subroutine is used to check collisions for particle inside
3298 c       the cube when the collision just happened is a collision including 
3299 c       collision with wall
3300 c       and update the afftected particles through chkcel
3301 
3302 c       input l,ii(specifying the direction of the wall collision),
3303 c          i1,i2,i3, (specifying the position of the cell 
3304 c                    we are going to check)
3305 c          t
3306 
3307         implicit double precision (a-h, o-z)
3308         SAVE   
3309 
3310         if (ii .eq. 1) then
3311            i = i1 - 2
3312            if (i .le. 0) i = i + 10
3313            ia = i
3314            do 1002 j = i2 - 1, i2 + 1
3315               do 1001 k = i3 - 1, i3 + 1
3316                  ib = j
3317                  ic = k
3318                  if (j .eq. 0) ib = 10
3319                  if (j .eq. 11) ib = 1
3320                  if (k .ge. 1 .and. k .le. 10) then
3321                     call dchcel(l, ia, ib, ic, t)
3322                  end if
3323  1001         continue
3324  1002      continue
3325         end if
3326 
3327         if (ii .eq. 2) then
3328            i = i1 + 2
3329            if (i .ge. 11) i = i - 10
3330            ia = i
3331            do 1004 j = i2 - 1, i2 + 1
3332               do 1003 k = i3 - 1, i3 + 1
3333                  ib = j
3334                  ic = k
3335                  if (j .eq. 0) ib = 10
3336                  if (j .eq. 11) ib = 1
3337                  if (k .ge. 1 .and. k .le. 10) then
3338                     call dchcel(l, ia, ib, ic, t)
3339                  end if
3340  1003         continue
3341  1004      continue
3342         end if
3343 
3344         if (ii .eq. 3) then
3345            j = i2 - 2
3346            if (j .le. 0) j = j + 10
3347            ib = j
3348            do 1006 i = i1 - 1, i1 + 1
3349               do 1005 k = i3 - 1, i3 + 1
3350                  ia = i
3351                  ic = k
3352                  if (i .eq. 0) ia = 10
3353                  if (i .eq. 11) ia = 1
3354                  if (k .ge. 1 .and. k .le. 10) then
3355                     call dchcel(l, ia, ib, ic, t)
3356                  end if
3357  1005         continue
3358  1006      continue
3359         end if
3360 
3361         if (ii .eq. 4) then
3362            j = i2 + 2
3363            if (j .ge. 11) j = j - 10
3364            ib = j
3365            do 1008 i = i1 - 1, i1 + 1
3366               do 1007 k = i3 - 1, i3 + 1
3367                  ia = i
3368                  ic = k
3369                  if (i .eq. 0) ia = 10
3370                  if (i .eq. 11) ia = 1
3371                  if (k .ge. 1 .and. k .le. 10) then
3372                     call dchcel(l, ia, ib, ic, t)
3373                  end if
3374  1007         continue
3375  1008      continue
3376         end if
3377 
3378         if (ii .eq. 5) then
3379            if (i3 .eq. 2) goto 100
3380            k = i3 - 2
3381            ic = k
3382            do 1010 i = i1 - 1, i1 + 1
3383               do 1009 j = i2 - 1, i2 + 1
3384                  ia = i
3385                  ib = j
3386                  if (i .eq. 0) ia = 10
3387                  if (i .eq. 11) ia = 1
3388                  if (j .eq. 0) ib = 10
3389                  if (j .eq. 11) ib = 1
3390                      call dchcel(l, ia, ib, ic, t)
3391  1009         continue
3392  1010      continue
3393         end if
3394 
3395         if (ii .eq. 6) then
3396            if (i3 .eq. 9) goto 100
3397            k = i3 + 2
3398            ic = k
3399            do 1012 i = i1 - 1, i1 + 1
3400               do 1011 j = i2 - 1, i2 + 1
3401                  ia = i
3402                  ib = j
3403                  if (i .eq. 0) ia = 10
3404                  if (i .eq. 11) ia = 1
3405                  if (j .eq. 0) ib = 10
3406                  if (j .eq. 11) ib = 1
3407                      call dchcel(l, ia, ib, ic, t)
3408  1011         continue
3409  1012      continue
3410         end if
3411 
3412  100        continue
3413 
3414         return
3415         end
3416 
3417         subroutine dchin3(l, ii, i1, i2, i3, t)
3418 c       this subroutine is used to check collisions for particle inside
3419 c       the cube when the collision just happened is a collision including 
3420 c       collision with wall
3421 c       and update the afftected particles through chkcel
3422 
3423 c       input l,ii(specifying the direction of the wall collision),
3424 c          i1,i2,i3, (specifying the position of the cell 
3425 c                    we are going to check)
3426 c          t
3427 
3428         implicit double precision (a-h, o-z)
3429         SAVE   
3430 
3431         if (ii .eq. 1) then
3432            i = i1 - 2
3433            if (i .le. 0) i = i + 10
3434            ia = i
3435            do 1002 j = i2 - 1, i2 + 1
3436               do 1001 k = i3 - 1, i3 + 1
3437                  ib = j
3438                  ic = k
3439                  if (j .eq. 0) ib = 10
3440                  if (j .eq. 11) ib = 1
3441                  if (k .eq. 0) ic = 10
3442                  if (k .eq. 11) ic = 1
3443                  call dchcel(l, ia, ib, ic, t)
3444  1001         continue
3445  1002      continue
3446         end if
3447 
3448         if (ii .eq. 2) then
3449            i = i1 + 2
3450            if (i .ge. 11) i = i - 10
3451            ia = i
3452            do 1004 j = i2 - 1, i2 + 1
3453               do 1003 k = i3 - 1, i3 + 1
3454                  ib = j
3455                  ic = k
3456                  if (j .eq. 0) ib = 10
3457                  if (j .eq. 11) ib = 1
3458                  if (k .eq. 0) ic = 10
3459                  if (k .eq. 11) ic = 1
3460                  call dchcel(l, ia, ib, ic, t)
3461  1003         continue
3462  1004      continue
3463         end if
3464 
3465         if (ii .eq. 3) then
3466            j = i2 - 2
3467            if (j .le. 0) j = j + 10
3468            ib = j
3469            do 1006 i = i1 - 1, i1 + 1
3470               do 1005 k = i3 - 1, i3 + 1
3471                  ia = i
3472                  ic = k
3473                  if (i .eq. 0) ia = 10
3474                  if (i .eq. 11) ia = 1
3475                  if (k .eq. 0) ic = 10
3476                  if (k .eq. 11) ic = 1
3477                  call dchcel(l, ia, ib, ic, t)
3478  1005         continue
3479  1006      continue
3480         end if
3481 
3482         if (ii .eq. 4) then
3483            j = i2 + 2
3484            if (j .ge. 11) j = j - 10
3485            ib = j
3486            do 1008 i = i1 - 1, i1 + 1
3487               do 1007 k = i3 - 1, i3 + 1
3488                  ia = i
3489                  ic = k
3490                  if (i .eq. 0) ia = 10
3491                  if (i .eq. 11) ia = 1
3492                  if (k .eq. 0) ic = 10
3493                  if (k .eq. 11) ic = 1
3494                  call dchcel(l, ia, ib, ic, t)
3495  1007         continue
3496  1008      continue
3497         end if
3498 
3499         if (ii .eq. 5) then
3500            k = i3 - 2
3501            if (k .le. 0) k = k + 10
3502            ic = k
3503            do 1010 i = i1 - 1, i1 + 1
3504               do 1009 j = i2 - 1, i2 + 1
3505                  ia = i
3506                  ib = j
3507                  if (i .eq. 0) ia = 10
3508                  if (i .eq. 11) ia = 1
3509                  if (j .eq. 0) ib = 10
3510                  if (j .eq. 11) ib = 1
3511                      call dchcel(l, ia, ib, ic, t)
3512  1009         continue
3513  1010      continue
3514         end if
3515 
3516         if (ii .eq. 6) then
3517            k = i3 + 2
3518            if (k .ge. 11) k = k - 10
3519            ic = k
3520            do 1012 i = i1 - 1, i1 + 1
3521               do 1011 j = i2 - 1, i2 + 1
3522                  ia = i
3523                  ib = j
3524                  if (i .eq. 0) ia = 10
3525                  if (i .eq. 11) ia = 1
3526                  if (j .eq. 0) ib = 10
3527                  if (j .eq. 11) ib = 1
3528                      call dchcel(l, ia, ib, ic, t)
3529  1011         continue
3530  1012      continue
3531         end if
3532 c
3533         return
3534         end
3535 
3536         subroutine dchcel(l, i, j, k, t)
3537 c       this subroutine is used to recalculate next collision time for 
3538 c       particles in the cell i,j,k if the next collision partener is l
3539 
3540         implicit double precision (a-h, o-z)
3541         parameter (MAXPTN=400001)
3542 
3543         common /ilist1/
3544      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
3545      &     ictype, icsta(MAXPTN),
3546      &     nic(MAXPTN), icels(MAXPTN)
3547 cc      SAVE /ilist1/
3548         common /ilist2/ icell, icel(10, 10, 10)
3549 cc      SAVE /ilist2/
3550         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3551 cc      SAVE /ilist5/
3552         SAVE   
3553 
3554         if (i .eq. 11 .or. j .eq. 11 .or. k .eq. 11) then
3555            if ( .not. (i .eq. 11 .and. j .eq. 11 .and.
3556      &     k .eq. 11)) stop 'cerr'
3557            m = icell
3558         else
3559            m = icel(i, j, k)
3560         end if
3561 
3562         if (m .eq. 0) return
3563         if (next(m) .eq. l) then
3564            tm = tlarge
3565            last0 = 0
3566            call reor(t, tm, m, last0)
3567         end if
3568         n = nic(m)
3569         if (n .eq. 0) return
3570         do 10 while(n .ne. m)
3571            if (next(n) .eq. l) then
3572               tm = tlarge
3573               last0 = 0
3574               call reor(t, tm, n, last0)
3575            end if
3576            n = nic(n)
3577  10        continue
3578 
3579         return
3580         end
3581 
3582         subroutine fixtim(l, t, tmin1, tmin, nc)
3583 c       this subroutine is used to compare the collision time with wall tmin1
3584 c       and new collision time with particles for particle l
3585 c       when used in ulist, input nc may be 0, which indicates no particle
3586 c       collisions happen before wall collision, of course, then tmin=tmin1
3587 
3588         implicit double precision (a-h, o-z)
3589         parameter (MAXPTN=400001)
3590 
3591         common /ilist1/
3592      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
3593      &     ictype, icsta(MAXPTN),
3594      &     nic(MAXPTN), icels(MAXPTN)
3595 cc      SAVE /ilist1/
3596         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3597 cc      SAVE /ilist5/
3598         SAVE   
3599 
3600         k = nc
3601         if (tmin .lt. tmin1) then
3602            ot(l) = tmin
3603            if (ct(l) .lt. tmin1) then
3604               icsta(l) = 0
3605            else
3606               icsta(l) = icsta(l) + 10
3607            end if
3608            next(l) = k
3609         else if (tmin .eq. tmin1) then
3610            ot(l) = tmin
3611            if (nc .eq. 0) then
3612               next(l) = 0
3613            else
3614               icsta(l) = icsta(l) + 10
3615               next(l) = k
3616            end if
3617         else
3618            ot(l) = tmin1
3619            next(l) = 0
3620         end if
3621         
3622         return
3623         end
3624 
3625         subroutine ud2(i, j, t, tmin, nc)
3626 c       this subroutine is used to update next(i), ct(i), ot(i),
3627 c        and get tmin, nc for j
3628 
3629         implicit double precision (a-h, o-z)
3630         parameter (MAXPTN=400001)
3631 
3632         common /para5/ iconfg, iordsc
3633 cc      SAVE /para5/
3634         common /aurec1/ jxa, jya, jza
3635 cc      SAVE /aurec1/
3636         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
3637 cc      SAVE /aurec2/
3638         common /ilist1/
3639      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
3640      &     ictype, icsta(MAXPTN),
3641      &     nic(MAXPTN), icels(MAXPTN)
3642 cc      SAVE /ilist1/
3643         common /ilist3/ size1, size2, size3, v1, v2, v3, size
3644 cc      SAVE /ilist3/
3645         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3646 cc      SAVE /ilist5/
3647         SAVE   
3648 
3649         logical allok
3650 
3651         call isco(i, j, allok, tm, t1, t2)
3652 
3653         if (allok) then
3654 c       tm eq tmin, change nc to make sure fixtime get the collision with both 
3655 c       wall and particle
3656 
3657              if (tm .lt. tmin) then
3658                 tmin = tm
3659                 ct(j) = t2
3660                 nc = i
3661                 if (iconfg .eq. 3 .or. iconfg .eq. 5) then
3662                    dgxa(j) = jxa * 10d0 * size1
3663                    dgya(j) = jya * 10d0 * size2
3664                    if (iconfg .eq. 5) then
3665                       dgza(j) = jza * 10d0 * size3
3666                    end if
3667                 end if
3668              end if
3669 
3670              if (tm .le. ot(i)) then
3671                 ct(i) = t1
3672                 icels0 = icels(i)
3673                 i1 = icels0 / 10000
3674                 i2 = (icels0 - i1 * 10000) / 100
3675                 i3 = icels0 - i1 * 10000 - i2 * 100
3676                 call wallc(i, i1, i2, i3, t, tmin1)
3677                 call fixtim(i, t, tmin1, tm, j)
3678                 if (iconfg .eq. 3 .or. iconfg .eq. 5) then
3679                    dgxa(i) = - jxa * 10d0 * size1
3680                    dgya(i) = - jya * 10d0 * size2
3681                    if (iconfg .eq. 5) then
3682                       dgza(i) = - jza * 10d0 * size3
3683                    end if
3684                 end if
3685              end if
3686 
3687              if (tm .gt. ot(i) .and. next(i) .eq. j) then
3688                 ct(i) = t1
3689                 call reor(t, tm, i, j)
3690              end if
3691 
3692            else if (next(i) .eq. j) then
3693 
3694              tm = tlarge
3695                 
3696              call reor(t, tm, i, j)
3697 
3698           end if
3699 
3700         return
3701         end
3702 
3703         subroutine isco(i, j, allok, tm, t1, t2)
3704 
3705         implicit double precision (a-h, o-z)
3706 
3707         common /para5/ iconfg, iordsc
3708 cc      SAVE /para5/
3709         SAVE   
3710 
3711         logical allok
3712 
3713         iorder = iordsc / 10
3714         if (iconfg .eq. 1) then
3715            if (iorder .eq. 1) then
3716               call isco1(i, j, allok, tm, t1, t2)
3717            else if (iorder .eq. 2) then
3718               call isco2(i, j, allok, tm, t1, t2)
3719            else if (iorder .eq. 3) then
3720               call isco3(i, j, allok, tm, t1, t2)
3721            end if
3722         else if (iconfg .eq. 2 .or. iconfg .eq. 4) then
3723            if (iorder .eq. 1) then
3724               call isco4(i, j, allok, tm, t1, t2)
3725            else if (iorder .eq. 2) then
3726               call isco5(i, j, allok, tm, t1, t2)
3727            else if (iorder .eq. 3) then
3728               call isco6(i, j, allok, tm, t1, t2)
3729            end if
3730         else if (iconfg .eq. 3) then
3731            if (iorder .eq. 1) then
3732               call isco7(i, j, allok, tm, t1, t2)
3733            else if (iorder .eq. 2) then
3734               call isco8(i, j, allok, tm, t1, t2)
3735            else if (iorder .eq. 3) then
3736               call isco9(i, j, allok, tm, t1, t2)
3737            end if
3738         else if (iconfg .eq. 5) then
3739            if (iorder .eq. 1) then
3740               call isco10(i, j, allok, tm, t1, t2)
3741            else if (iorder .eq. 2) then
3742               call isco11(i, j, allok, tm, t1, t2)
3743            else if (iorder .eq. 3) then
3744               call isco12(i, j, allok, tm, t1, t2)
3745            end if
3746         end if
3747 
3748         return
3749         end
3750 
3751         subroutine isco1(i, j, allok, tm, t1, t2)
3752 c       this subroutine is used to decide whether there is a collision between
3753 c       particle i and j, if there is one allok=1, and tm gives the 
3754 c       collision time, t1 the collision time for i,
3755 c       t2 the collision time for j
3756 
3757         implicit double precision (a-h, o-z)
3758         parameter (MAXPTN=400001)
3759 
3760         common /para2/ xmp, xmu, alpha, rscut2, cutof2
3761 cc      SAVE /para2/
3762         common /para5/ iconfg, iordsc
3763 cc      SAVE /para5/
3764         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
3765      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
3766      &       xmass(MAXPTN), ityp(MAXPTN)
3767 cc      SAVE /prec2/
3768         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
3769 cc      SAVE /prec4/
3770         common /ilist1/
3771      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
3772      &     ictype, icsta(MAXPTN),
3773      &     nic(MAXPTN), icels(MAXPTN)
3774 cc      SAVE /ilist1/
3775         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3776 cc      SAVE /ilist5/
3777         SAVE   
3778 
3779         logical allok
3780 
3781 c       preventing consecutive collisions
3782         allok = last(i) .ne. j .or. last(j) .ne. i
3783 
3784 c       set up numbers for later calculations
3785         i1 = i
3786         i2 = j
3787 
3788         p4 = ft(i2) - ft(i1)
3789         p1 = gx(i2) - gx(i1)
3790         p2 = gy(i2) - gy(i1)
3791         p3 = gz(i2) - gz(i1)
3792 
3793         q4 = e(i1)
3794         q1 = px(i1)
3795         q2 = py(i1)
3796         q3 = pz(i1)
3797 
3798         r4 = e(i2)
3799         r1 = px(i2)
3800         r2 = py(i2)
3801         r3 = pz(i2)
3802 
3803         a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
3804         b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
3805         c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
3806         d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
3807         ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
3808         f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
3809 
3810 c       make sure particle 2 formed early
3811         h = a + b
3812         if (h .gt. 0d0) then
3813            g = a
3814            a = -b
3815            b = -g
3816 
3817            g = c
3818            c = d
3819            d = g
3820 
3821            i1 = j
3822            i2 = i
3823         end if
3824 
3825 c       check the approaching criteria
3826         if (allok) then
3827 
3828            vp = a * d - b * ee
3829 
3830            allok = allok .and. vp .lt. 0d0
3831 
3832         end if
3833 
3834 c       check the closest approach distance criteria
3835          if (allok) then
3836 
3837            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
3838      &           (ee ** 2 - c * d)
3839 
3840            allok = allok .and. dm2 .lt. cutof2
3841 
3842         end if
3843 
3844 c       check the time criteria
3845         if (allok) then
3846 
3847            tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
3848            tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
3849            tm = 0.5d0 * (tc1 + tc2)
3850 
3851            allok = allok .and. tm .gt. ft(i) .and. tm .gt. ft(j)
3852 
3853         end if
3854 
3855 c        check rts cut
3856         if (allok) then
3857 
3858            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
3859      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
3860 
3861            allok = allok .and. rts2 .gt. rscut2
3862         end if
3863           
3864         if (.not. allok) then
3865            tm = tlarge
3866            t1 = tlarge
3867            t2 = tlarge
3868         else if (h .gt. 0d0) then
3869            t1 = tm
3870            t2 = tm
3871         else
3872            t1 = tm
3873            t2 = tm
3874         end if
3875 
3876         return
3877         end
3878 
3879         subroutine isco2(i, j, allok, tm, t1, t2)
3880 c       this subroutine is used to decide whether there is a collision between
3881 c       particle i and j, if there is one allok=1, and tm gives the 
3882 c       collision time, t1 the collision time for i,
3883 c       t2 the collision time for j
3884 
3885         implicit double precision (a-h, o-z)
3886         parameter (MAXPTN=400001)
3887 
3888         common /para2/ xmp, xmu, alpha, rscut2, cutof2
3889 cc      SAVE /para2/
3890         common /para5/ iconfg, iordsc
3891 cc      SAVE /para5/
3892         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
3893      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
3894      &       xmass(MAXPTN), ityp(MAXPTN)
3895 cc      SAVE /prec2/
3896         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
3897 cc      SAVE /prec4/
3898         common /ilist1/
3899      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
3900      &     ictype, icsta(MAXPTN),
3901      &     nic(MAXPTN), icels(MAXPTN)
3902 cc      SAVE /ilist1/
3903         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
3904 cc      SAVE /ilist5/
3905         SAVE   
3906 
3907         logical allok
3908 
3909 c       preventing consecutive collisions
3910         allok = last(i) .ne. j .or. last(j) .ne. i
3911 
3912 c       set up numbers for later calculations
3913         i1 = i
3914         i2 = j
3915 
3916         p4 = ft(i2) - ft(i1)
3917         p1 = gx(i2) - gx(i1)
3918         p2 = gy(i2) - gy(i1)
3919         p3 = gz(i2) - gz(i1)
3920 
3921         q4 = e(i1)
3922         q1 = px(i1)
3923         q2 = py(i1)
3924         q3 = pz(i1)
3925 
3926         r4 = e(i2)
3927         r1 = px(i2)
3928         r2 = py(i2)
3929         r3 = pz(i2)
3930 
3931         a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
3932         b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
3933         c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
3934         d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
3935         ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
3936         f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
3937 
3938 c       make sure particle 2 formed early
3939         h = a + b
3940         if (h .gt. 0d0) then
3941            g = a
3942            a = -b
3943            b = -g
3944 
3945            g = c
3946            c = d
3947            d = g
3948 
3949            i1 = j
3950            i2 = i
3951         end if
3952 
3953 c       check the approaching criteria
3954         if (allok) then
3955 
3956            vp = a * d - b * ee
3957 
3958            allok = allok .and. vp .lt. 0d0
3959 
3960         end if
3961 
3962 c       check the closest approach distance criteria
3963          if (allok) then
3964 
3965            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
3966      &          (ee ** 2 - c * d)
3967 
3968            allok = allok .and. dm2 .lt. cutof2
3969 
3970         end if
3971 
3972 c       check the time criteria
3973         if (allok) then
3974 
3975            tc1 = ft(i1) - e(i1) * (a * d - b * ee)/(ee ** 2 - c * d)
3976            tc2 = ft(i2) + e(i2) * (b * c - a * ee)/(ee ** 2 - c * d)
3977            if (iordsc .eq. 20) then
3978               tm = min(tc1, tc2)
3979            else if (iordsc .eq. 21) then
3980               tm = 0.5d0 * (tc1 + tc2)
3981            else
3982               tm = max(tc1, tc2)
3983            end if
3984 
3985            allok = allok .and. tm .gt. ft(i) .and. tm .gt. ft(j)
3986 
3987         end if
3988 
3989 c        check rts cut
3990         if (allok) then
3991 
3992            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
3993      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
3994 
3995            allok = allok .and. rts2 .gt. rscut2
3996         end if
3997           
3998         if (.not. allok) then
3999            tm = tlarge
4000            t1 = tlarge
4001            t2 = tlarge
4002         else if (h .gt. 0d0) then
4003            t1 = tc2
4004            t2 = tc1
4005         else
4006            t1 = tc1
4007            t2 = tc2
4008         end if
4009 
4010         return
4011         end
4012 
4013         subroutine isco3(i, j, allok, tm, t1, t2)
4014 c       this subroutine is used to decide whether there is a collision between
4015 c       particle i and j, if there is one allok=1, and tm gives the 
4016 c       collision time, t1 the collision time for i,
4017 c       t2 the collision time for j
4018 
4019         implicit double precision (a-h, o-z)
4020         parameter (MAXPTN=400001)
4021 
4022         common /para2/ xmp, xmu, alpha, rscut2, cutof2
4023 cc      SAVE /para2/
4024         common /para5/ iconfg, iordsc
4025 cc      SAVE /para5/
4026         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4027      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4028      &       xmass(MAXPTN), ityp(MAXPTN)
4029 cc      SAVE /prec2/
4030         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4031 cc      SAVE /prec4/
4032         common /ilist1/
4033      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
4034      &     ictype, icsta(MAXPTN),
4035      &     nic(MAXPTN), icels(MAXPTN)
4036 cc      SAVE /ilist1/
4037         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4038 cc      SAVE /ilist5/  
4039         SAVE   
4040 
4041         logical allok
4042 
4043 c       preventing consecutive collisions
4044         allok = last(i) .ne. j .or. last(j) .ne. i
4045 
4046         if (ft(i) .ge. ft(j)) then
4047            i1 = j
4048            i2 = i
4049         else 
4050            i1 = i
4051            i2 = j
4052         end if
4053         
4054         if (allok) then
4055 
4056            t1 = ft(i1)
4057            vx1 = vx(i1)
4058            vy1 = vy(i1)
4059            vz1 = vz(i1)
4060 
4061            t2 = ft(i2)
4062 
4063            dvx = vx(i2) - vx1
4064            dvy = vy(i2) - vy1
4065            dvz = vz(i2) - vz1
4066 
4067            dt = t2 - t1
4068 
4069            dx = gx(i2) - gx(i1) - vx1 * dt
4070            dy = gy(i2) - gy(i1) - vy1 * dt
4071            dz = gz(i2) - gz(i1) - vz1 * dt
4072 
4073            vp = dvx * dx + dvy * dy + dvz * dz
4074 
4075            allok = allok .and. vp .lt. 0d0
4076 
4077         end if
4078 
4079         if (allok) then
4080 
4081            v2= dvx * dvx + dvy * dvy + dvz * dvz
4082 
4083            if (v2 .eq. 0d0) then
4084               tm = tlarge
4085            else
4086               tm = t2 - vp / v2
4087            end if
4088 
4089 c       note now tm is the absolute time
4090 
4091            allok = allok .and. tm .gt. t1 .and. tm .gt. t2
4092 
4093         end if
4094 
4095         if (allok) then
4096 
4097            dgx = dx - dvx * t2
4098            dgy = dy - dvy * t2
4099            dgz = dz - dvz * t2
4100 
4101            dm2 = - v2 * tm ** 2  + dgx * dgx + dgy * dgy + dgz * dgz
4102 
4103            allok = allok .and. dm2 .lt. cutof2
4104 
4105         end if
4106         
4107         if (allok) then
4108 
4109            e1 = e(i1)
4110            px1 = px(i1)
4111            py1 = py(i1)
4112            pz1 = pz(i1)
4113            e2 = e(i2)
4114            px2 = px(i2)
4115            py2 = py(i2)
4116            pz2 = pz(i2)
4117 
4118            rts2 = (e1 + e2) ** 2 - (px1 + px2) ** 2
4119      &          - (py1 + py2) ** 2 - (pz1 + pz2) ** 2
4120 
4121            allok = allok .and. rts2 .gt. rscut2
4122         end if
4123 
4124         if (.not. allok) then
4125            tm = tlarge
4126            t1 = tlarge
4127            t2 = tlarge
4128         else
4129            t1 = tm
4130            t2 = tm
4131         end if
4132 
4133         return
4134         end
4135 
4136         subroutine isco4(i, j, allok, tm, t1, t2)
4137 c       this subroutine is used to decide whether there is a collision between
4138 c       particle i and j, if there is one allok=1, and tm gives the 
4139 c       collision time, t1 the collision time for i,
4140 c       t2 the collision time for j
4141 
4142         implicit double precision (a-h, o-z)
4143         parameter (MAXPTN=400001)
4144 
4145         common /para2/ xmp, xmu, alpha, rscut2, cutof2
4146 cc      SAVE /para2/
4147         common /para5/ iconfg, iordsc
4148 cc      SAVE /para5/
4149         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4150      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4151      &       xmass(MAXPTN), ityp(MAXPTN)
4152 cc      SAVE /prec2/
4153         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4154 cc      SAVE /prec4/
4155         common /ilist1/
4156      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
4157      &     ictype, icsta(MAXPTN),
4158      &     nic(MAXPTN), icels(MAXPTN)
4159 cc      SAVE /ilist1/
4160         common /ilist3/ size1, size2, size3, v1, v2, v3, size
4161 cc      SAVE /ilist3/
4162         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4163 cc      SAVE /ilist5/
4164         SAVE   
4165 
4166         logical allok
4167 
4168 c       preventing consecutive collisions
4169         allok = last(i) .ne. j .or. last(j) .ne. i
4170 
4171 c       set up numbers for later calculations
4172 
4173         icels1 = icels(i)
4174         ii1 = icels1 / 10000
4175         jj1 = (icels1 - ii1 * 10000) / 100
4176         kk1 = icels1 - ii1 * 10000 - jj1 * 100
4177         icels2 = icels(j)
4178         ii2 = icels2 / 10000
4179         jj2 = (icels2 - ii2 * 10000) / 100
4180         kk2 = icels2 - ii2 * 10000 - jj2 * 100
4181 
4182         i1 = i
4183         i2 = j
4184 
4185         p4 = ft(i2) - ft(i1)
4186         p1 = gx(i2) - gx(i1)
4187         p2 = gy(i2) - gy(i1)
4188         p3 = gz(i2) - gz(i1)
4189 
4190         if (ii1 - ii2 .gt. 5) then
4191            p1 = p1 + 10d0 * size1
4192         else if (ii1 - ii2 .lt. -5) then
4193            p1 = p1 - 10d0 * size1
4194         end if
4195         if (jj1 - jj2 .gt. 5) then
4196            p2 = p2 + 10d0 * size2
4197         else if (jj1 - jj2 .lt. -5) then
4198            p2 = p2 - 10d0 * size2
4199         end if
4200         if (kk1 - kk2 .gt. 5) then
4201            p3 = p3 + 10d0 * size3
4202         else if (kk1 - kk2 .lt. -5) then
4203            p3 = p3 - 10d0 * size3
4204         end if
4205 
4206         q4 = e(i1)
4207         q1 = px(i1)
4208         q2 = py(i1)
4209         q3 = pz(i1)
4210 
4211         r4 = e(i2)
4212         r1 = px(i2)
4213         r2 = py(i2)
4214         r3 = pz(i2)
4215 
4216         a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
4217         b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
4218         c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
4219         d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
4220         ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
4221         f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
4222 
4223 c       make sure particle 2 formed early
4224         h = a + b
4225         if (h .gt. 0d0) then
4226            g = a
4227            a = -b
4228            b = -g
4229 
4230            g = c
4231            c = d
4232            d = g
4233 
4234            i1 = j
4235            i2 = i
4236         end if
4237 
4238 c       check the approaching criteria
4239         if (allok) then
4240 
4241            vp = a * d - b * ee
4242 
4243            allok = allok .and. vp .lt. 0d0
4244 
4245         end if
4246 
4247 c       check the closest approach distance criteria
4248          if (allok) then
4249 
4250            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
4251      &           (ee ** 2 - c * d)
4252 
4253            allok = allok .and. dm2 .lt. cutof2
4254 
4255         end if
4256 
4257 c       check the time criteria
4258         if (allok) then
4259 
4260            tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
4261            tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
4262            tm = 0.5d0 * (tc1 + tc2)
4263 
4264            allok = allok .and. tm .gt. ft(i) .and. tm .gt. ft(j)
4265 
4266         end if
4267 
4268 c        check rts cut
4269         if (allok) then
4270 
4271            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
4272      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
4273 
4274            allok = allok .and. rts2 .gt. rscut2
4275         end if
4276           
4277         if (.not. allok) then
4278            tm = tlarge
4279            t1 = tlarge
4280            t2 = tlarge
4281         else if (h .gt. 0d0) then
4282            t1 = tm
4283            t2 = tm
4284         else
4285            t1 = tm
4286            t2 = tm
4287         end if
4288 
4289         return
4290         end
4291 
4292         subroutine isco5(i, j, allok, tm, t1, t2)
4293 c       this subroutine is used to decide whether there is a collision between
4294 c       particle i and j, if there is one allok=1, and tm gives the 
4295 c       collision time, t1 the collision time for i,
4296 c       t2 the collision time for j
4297 
4298         implicit double precision (a-h, o-z)
4299         parameter (MAXPTN=400001)
4300 
4301         common /para2/ xmp, xmu, alpha, rscut2, cutof2
4302 cc      SAVE /para2/
4303         common /para5/ iconfg, iordsc
4304 cc      SAVE /para5/
4305         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4306      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4307      &       xmass(MAXPTN), ityp(MAXPTN)
4308 cc      SAVE /prec2/
4309         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4310 cc      SAVE /prec4/
4311         common /ilist1/
4312      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
4313      &     ictype, icsta(MAXPTN),
4314      &     nic(MAXPTN), icels(MAXPTN)
4315 cc      SAVE /ilist1/
4316         common /ilist3/ size1, size2, size3, v1, v2, v3, size
4317 cc      SAVE /ilist3/
4318         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4319 cc      SAVE /ilist5/
4320         SAVE   
4321 
4322         logical allok
4323 
4324 c       preventing consecutive collisions
4325         allok = last(i) .ne. j .or. last(j) .ne. i
4326 
4327 c       set up numbers for later calculations
4328 
4329         icels1 = icels(i)
4330         ii1 = icels1 / 10000
4331         jj1 = (icels1 - ii1 * 10000) / 100
4332         kk1 = icels1 - ii1 * 10000 - jj1 * 100
4333         icels2 = icels(j)
4334         ii2 = icels2 / 10000
4335         jj2 = (icels2 - ii2 * 10000) / 100
4336         kk2 = icels2 - ii2 * 10000 - jj2 * 100
4337 
4338         i1 = i
4339         i2 = j
4340 
4341         p4 = ft(i2) - ft(i1)
4342         p1 = gx(i2) - gx(i1)
4343         p2 = gy(i2) - gy(i1)
4344         p3 = gz(i2) - gz(i1)
4345 
4346         if (ii1 - ii2 .gt. 5) then
4347            p1 = p1 + 10d0 * size1
4348         else if (ii1 - ii2 .lt. -5) then
4349            p1 = p1 - 10d0 * size1
4350         end if
4351         if (jj1 - jj2 .gt. 5) then
4352            p2 = p2 + 10d0 * size2
4353         else if (jj1 - jj2 .lt. -5) then
4354            p2 = p2 - 10d0 * size2
4355         end if
4356         if (kk1 - kk2 .gt. 5) then
4357            p3 = p3 + 10d0 * size3
4358         else if (kk1 - kk2 .lt. -5) then
4359            p3 = p3 - 10d0 * size3
4360         end if
4361 
4362         q4 = e(i1)
4363         q1 = px(i1)
4364         q2 = py(i1)
4365         q3 = pz(i1)
4366 
4367         r4 = e(i2)
4368         r1 = px(i2)
4369         r2 = py(i2)
4370         r3 = pz(i2)
4371 
4372         a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
4373         b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
4374         c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
4375         d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
4376         ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
4377         f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
4378 
4379 c       make sure particle 2 formed early
4380         h = a + b
4381         if (h .gt. 0d0) then
4382            g = a
4383            a = -b
4384            b = -g
4385 
4386            g = c
4387            c = d
4388            d = g
4389 
4390            i1 = j
4391            i2 = i
4392         end if
4393 
4394 c       check the approaching criteria
4395         if (allok) then
4396 
4397            vp = a * d - b * ee
4398 
4399            allok = allok .and. vp .lt. 0d0
4400 
4401         end if
4402 
4403 c       check the closest approach distance criteria
4404          if (allok) then
4405 
4406            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
4407      &           (ee ** 2 - c * d)
4408 
4409            allok = allok .and. dm2 .lt. cutof2
4410 
4411         end if
4412 
4413 c       check the time criteria
4414         if (allok) then
4415 
4416            tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
4417            tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
4418            if (iordsc .eq. 20) then
4419               tm = min(tc1, tc2)
4420            else if (iordsc .eq. 21) then
4421               tm = 0.5d0 * (tc1 + tc2)
4422            else
4423               tm = max(tc1, tc2)
4424            end if
4425 
4426            allok = allok .and. tm .gt. ft(i) .and. tm .gt. ft(j)
4427 
4428         end if
4429 
4430 c        check rts cut
4431         if (allok) then
4432 
4433            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
4434      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
4435 
4436            allok = allok .and. rts2 .gt. rscut2
4437         end if
4438           
4439         if (.not. allok) then
4440            tm = tlarge
4441            t1 = tlarge
4442            t2 = tlarge
4443         else if (h .gt. 0d0) then
4444            t1 = tc2
4445            t2 = tc1
4446         else
4447            t1 = tc1
4448            t2 = tc2
4449         end if
4450 
4451         return
4452         end
4453 
4454         subroutine isco6(i, j, allok, tm, t1, t2)
4455 c       this subroutine is used to decide whether there is a collision between
4456 c       particle i and j, if there is one allok=1, and tm gives the 
4457 c       collision time, t1 the collision time for i,
4458 c       t2 the collision time for j
4459 
4460         implicit double precision (a-h, o-z)
4461         parameter (MAXPTN=400001)
4462 
4463         common /para2/ xmp, xmu, alpha, rscut2, cutof2
4464 cc      SAVE /para2/
4465         common /para5/ iconfg, iordsc
4466 cc      SAVE /para5/
4467         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4468      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4469      &       xmass(MAXPTN), ityp(MAXPTN)
4470 cc      SAVE /prec2/
4471         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4472 cc      SAVE /prec4/
4473         common /ilist1/
4474      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
4475      &     ictype, icsta(MAXPTN),
4476      &     nic(MAXPTN), icels(MAXPTN)
4477 cc      SAVE /ilist1/
4478         common /ilist3/ size1, size2, size3, v1, v2, v3, size
4479 cc      SAVE /ilist3/
4480         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4481 cc      SAVE /ilist5/
4482         SAVE   
4483 
4484         logical allok
4485 
4486 c       preventing consecutive collisions
4487         allok = last(i) .ne. j .or. last(j) .ne. i
4488 
4489         if (ft(i) .ge. ft(j)) then
4490            i1 = j
4491            i2 = i
4492         else 
4493            i1 = i
4494            i2 = j
4495         end if
4496 
4497         icels1 = icels(i1)
4498         ii1 = icels1 / 10000
4499         jj1 = (icels1 - ii1 * 10000) / 100
4500         kk1 = icels1 - ii1 * 10000 - jj1 * 100
4501         icels2 = icels(i2)
4502         ii2 = icels2 / 10000
4503         jj2 = (icels2 - ii2 * 10000) / 100
4504         kk2 = icels2 - ii2 * 10000 - jj2 * 100
4505         
4506         if (allok) then
4507 
4508            t1 = ft(i1)
4509            vx1 = vx(i1)
4510            vy1 = vy(i1)
4511            vz1 = vz(i1)
4512 
4513            t2 = ft(i2)
4514 
4515            dvx = vx(i2) - vx1
4516            dvy = vy(i2) - vy1
4517            dvz = vz(i2) - vz1
4518 
4519            dt = t2 - t1
4520 
4521            dx = gx(i2) - gx(i1) - vx1 * dt
4522            dy = gy(i2) - gy(i1) - vy1 * dt
4523            dz = gz(i2) - gz(i1) - vz1 * dt
4524 
4525            if (ii1 - ii2 .gt. 5) then
4526               dx = dx + 10d0 * size1
4527            else if (ii1 - ii2 .lt. - 5) then
4528               dx = dx - 10d0 * size1
4529            end if
4530 
4531            if (jj1 - jj2 .gt. 5) then
4532               dy = dy + 10d0 * size2
4533            else if (jj1 - jj2 .lt. - 5) then
4534               dy = dy - 10d0 * size2
4535            end if
4536 
4537            if (kk1 - kk2 .gt. 5) then
4538               dz = dz + 10d0 * size3
4539            else if (kk1 - kk2 .lt. -5) then
4540               dz = dz - 10d0 * size3
4541            end if
4542 
4543            vp = dvx * dx + dvy * dy + dvz * dz
4544 
4545            allok = allok .and. vp .lt. 0d0
4546 
4547         end if
4548 
4549         if (allok) then
4550 
4551            v2p = dvx * dvx + dvy * dvy + dvz * dvz
4552 
4553            if (v2p .eq. 0d0) then
4554               tm = tlarge
4555            else
4556               tm = t2 - vp / v2p
4557            end if
4558 
4559 c       note now tm is the absolute time
4560 
4561            allok = allok .and. tm .gt. t1 .and. tm .gt. t2
4562 
4563         end if
4564 
4565         if (allok) then
4566 
4567            dgx = dx - dvx * t2
4568            dgy = dy - dvy * t2
4569            dgz = dz - dvz * t2
4570 
4571            dm2 = - v2p * tm ** 2  + dgx * dgx + dgy * dgy + dgz * dgz
4572 
4573            allok = allok .and. dm2 .lt. cutof2
4574 
4575         end if
4576         
4577         if (allok) then
4578 
4579            e1 = e(i1)
4580            px1 = px(i1)
4581            py1 = py(i1)
4582            pz1 = pz(i1)
4583            e2 = e(i2)
4584            px2 = px(i2)
4585            py2 = py(i2)
4586            pz2 = pz(i2)
4587 
4588            rts2 = (e1 + e2) ** 2 - (px1 + px2) ** 2
4589      &          - (py1 + py2) ** 2 - (pz1 + pz2) ** 2
4590 
4591            allok = allok .and. rts2 .gt. rscut2
4592         end if
4593 
4594         if (.not. allok) then
4595            tm = tlarge
4596            t1 = tlarge
4597            t2 = tlarge
4598         else
4599            t1 = tm
4600            t2 = tm
4601         end if
4602 
4603         return
4604         end
4605 
4606         subroutine isco7(i, j, allok, tm, t1, t2)
4607 c       this subroutine is used to decide whether there is a collision between
4608 c       particle i and j, if there is one allok=1, and tm gives the 
4609 c       collision time, t1 the collision time for i,
4610 c       t2 the collision time for j
4611 
4612         implicit double precision (a-h, o-z)
4613         parameter (MAXPTN=400001)
4614 
4615         common /para2/ xmp, xmu, alpha, rscut2, cutof2
4616 cc      SAVE /para2/
4617         common /para5/ iconfg, iordsc
4618 cc      SAVE /para5/
4619         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4620      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4621      &       xmass(MAXPTN), ityp(MAXPTN)
4622 cc      SAVE /prec2/
4623         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4624 cc      SAVE /prec4/
4625         common /aurec1/ jxa, jya, jza
4626 cc      SAVE /aurec1/
4627         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
4628 cc      SAVE /aurec2/
4629         common /ilist1/
4630      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
4631      &     ictype, icsta(MAXPTN),
4632      &     nic(MAXPTN), icels(MAXPTN)
4633 cc      SAVE /ilist1/
4634         common /ilist3/ size1, size2, size3, v1, v2, v3, size
4635 cc      SAVE /ilist3/
4636         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4637 cc      SAVE /ilist5/
4638         SAVE   
4639 
4640         logical allok, allokp
4641 
4642 c       preventing consecutive collisions
4643         allok = last(i) .ne. j .or. last(j) .ne. i
4644 
4645 c       set up numbers for later calculations
4646 
4647         tm = tlarge
4648 
4649         if (allok) then
4650            do 1000 ii = - 1, 1
4651               do 2000 jj = - 1, 1
4652 
4653                  allokp = .true.
4654                  
4655                  i1 = i
4656                  i2 = j
4657 
4658                  p4 = ft(j) - ft(i)
4659                  p1 = gx(j) - gx(i)
4660                  p2 = gy(j) - gy(i)
4661                  p3 = gz(j) - gz(i)
4662 
4663                  p1 = p1 + ii * 10d0 * size1
4664                  p2 = p2 + jj * 10d0 * size2
4665 
4666                  q4 = e(i)
4667                  q1 = px(i)
4668                  q2 = py(i)
4669                  q3 = pz(i)
4670                  
4671                  r4 = e(j)
4672                  r1 = px(j)
4673                  r2 = py(j)
4674                  r3 = pz(j)
4675 
4676                  a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
4677                  b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
4678                  c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
4679                  d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
4680                  ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
4681                  f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
4682 
4683 c       make sure particle 2 formed early
4684                  h = a + b
4685                  if (h .gt. 0d0) then
4686                     g = a
4687                     a = -b
4688                     b = -g
4689                     g = c
4690                     c = d
4691                     d = g
4692                     i1 = j
4693                     i2 = i
4694                  end if
4695                  
4696 c       check the approaching criteria
4697                  if (allokp) then
4698                     vp = a * d - b * ee
4699                     allokp = allokp .and. vp .lt. 0d0
4700                  end if
4701 
4702 c       check the closest approach distance criteria
4703                  if (allokp) then
4704            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
4705      &            (ee ** 2 - c * d)
4706                     allokp = allokp .and. dm2 .lt. cutof2
4707                  end if
4708 
4709 c       check the time criteria
4710                  if (allokp) then
4711            tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
4712            tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
4713            tmp = 0.5d0 * (tc1 + tc2)
4714            allokp = allokp .and. tmp .gt. ft(i) .and. tmp .gt. ft(j)
4715                  end if
4716 
4717                  if (allokp .and. tmp .lt. tm) then
4718                     tm = tmp
4719                     jxa = ii
4720                     jya = jj
4721 cd                    dgxa(j) = ii * 10d0 * size1
4722 cd                    dgya(j) = jj * 10d0 * size2
4723 cd                    dgxa(i) = - dgxa(j)
4724 cd                    dgya(i) = - dgya(j)
4725                  end if
4726 
4727  2000              continue
4728  1000           continue
4729 
4730            if (tm .eq. tlarge) then
4731               allok = .false.
4732            end if
4733            
4734         end if
4735 
4736 c        check rts cut
4737         if (allok) then
4738 
4739            q4 = e(i1)
4740            q1 = px(i1)
4741            q2 = py(i1)
4742            q3 = pz(i1)
4743 
4744            r4 = e(i2)
4745            r1 = px(i2)
4746            r2 = py(i2)
4747            r3 = pz(i2)
4748 
4749            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
4750      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
4751 
4752            allok = allok .and. rts2 .gt. rscut2
4753         end if
4754           
4755         if (.not. allok) then
4756            tm = tlarge
4757            t1 = tlarge
4758            t2 = tlarge
4759         else if (h .gt. 0d0) then
4760            t1 = tm
4761            t2 = tm
4762         else
4763            t1 = tm
4764            t2 = tm
4765         end if
4766 
4767         return
4768         end
4769 
4770         subroutine isco8(i, j, allok, tm, t1, t2)
4771 c       this subroutine is used to decide whether there is a collision between
4772 c       particle i and j, if there is one allok=1, and tm gives the 
4773 c       collision time, t1 the collision time for i,
4774 c       t2 the collision time for j
4775 
4776         implicit double precision (a-h, o-z)
4777         parameter (MAXPTN=400001)
4778 
4779         common /para2/ xmp, xmu, alpha, rscut2, cutof2
4780 cc      SAVE /para2/
4781         common /para5/ iconfg, iordsc
4782 cc      SAVE /para5/
4783         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4784      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4785      &       xmass(MAXPTN), ityp(MAXPTN)
4786 cc      SAVE /prec2/
4787         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4788 cc      SAVE /prec4/
4789         common /aurec1/ jxa, jya, jza
4790 cc      SAVE /aurec1/
4791         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
4792 cc      SAVE /aurec2/
4793         common /ilist1/
4794      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
4795      &     ictype, icsta(MAXPTN),
4796      &     nic(MAXPTN), icels(MAXPTN)
4797 cc      SAVE /ilist1/
4798         common /ilist3/ size1, size2, size3, v1, v2, v3, size
4799 cc      SAVE /ilist3/
4800         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4801 cc      SAVE /ilist5/
4802         SAVE   
4803 
4804         logical allok, allokp
4805 
4806 c       preventing consecutive collisions
4807         allok = last(i) .ne. j .or. last(j) .ne. i
4808 
4809 c       set up numbers for later calculations
4810 
4811         tm = tlarge
4812 
4813         if (allok) then
4814            do 1000 ii = - 1, 1
4815               do 2000 jj = - 1, 1
4816 
4817                  allokp = .true.
4818                  
4819                  i1 = i
4820                  i2 = j
4821 
4822                  p4 = ft(j) - ft(i)
4823                  p1 = gx(j) - gx(i)
4824                  p2 = gy(j) - gy(i)
4825                  p3 = gz(j) - gz(i)
4826 
4827                  p1 = p1 + ii * 10d0 * size1
4828                  p2 = p2 + jj * 10d0 * size2
4829 
4830                  q4 = e(i)
4831                  q1 = px(i)
4832                  q2 = py(i)
4833                  q3 = pz(i)
4834                  
4835                  r4 = e(j)
4836                  r1 = px(j)
4837                  r2 = py(j)
4838                  r3 = pz(j)
4839 
4840                  a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
4841                  b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
4842                  c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
4843                  d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
4844                  ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
4845                  f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
4846 
4847 c       make sure particle 2 formed early
4848                  h = a + b
4849                  if (h .gt. 0d0) then
4850                     g = a
4851                     a = -b
4852                     b = -g
4853                     g = c
4854                     c = d
4855                     d = g
4856                     i1 = j
4857                     i2 = i
4858                  end if
4859                  
4860 c       check the approaching criteria
4861                  if (allokp) then
4862                     vp = a * d - b * ee
4863                     allokp = allokp .and. vp .lt. 0d0
4864                  end if
4865 
4866 c       check the closest approach distance criteria
4867                  if (allokp) then
4868            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
4869      &            (ee ** 2 - c * d)
4870                     allokp = allokp .and. dm2 .lt. cutof2
4871                  end if
4872 
4873 c       check the time criteria
4874                  if (allokp) then
4875            tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
4876            tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
4877                     if (iordsc .eq. 20) then
4878                        tmp = min(tc1, tc2)
4879                     else if (iordsc .eq. 21) then
4880                        tmp = 0.5d0 * (tc1 + tc2)
4881                     else
4882                        tmp = max(tc1, tc2)
4883                     end if
4884            allokp = allokp .and. tmp .gt. ft(i) .and. tmp .gt. ft(j)
4885                  end if
4886 
4887                  if (allokp .and. tmp .lt. tm) then
4888                     tm = tmp
4889                     jxa = ii
4890                     jya = jj
4891                     ha = h
4892                     tc1a = tc1
4893                     tc2a = tc2
4894 cd                    dgxa(j) = ii * 10d0 * size1
4895 cd                    dgya(j) = jj * 10d0 * size2
4896 cd                    dgxa(i) = - dgxa(j)
4897 cd                    dgya(i) = - dgya(j)
4898                  end if
4899 
4900  2000              continue
4901  1000           continue
4902 
4903            if (tm .eq. tlarge) then
4904               allok = .false.
4905            end if
4906            
4907         end if
4908 
4909 c        check rts cut
4910         if (allok) then
4911 
4912            q4 = e(i1)
4913            q1 = px(i1)
4914            q2 = py(i1)
4915            q3 = pz(i1)
4916 
4917            r4 = e(i2)
4918            r1 = px(i2)
4919            r2 = py(i2)
4920            r3 = pz(i2)
4921 
4922            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
4923      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
4924 
4925            allok = allok .and. rts2 .gt. rscut2
4926         end if
4927           
4928         if (.not. allok) then
4929            tm = tlarge
4930            t1 = tlarge
4931            t2 = tlarge
4932         else if (ha .gt. 0d0) then
4933            t1 = tc2a
4934            t2 = tc1a
4935         else
4936            t1 = tc1a
4937            t2 = tc2a
4938         end if
4939 
4940         return
4941         end
4942 
4943         subroutine isco9(i, j, allok, tm, t1, t2)
4944 c       this subroutine is used to decide whether there is a collision between
4945 c       particle i and j, if there is one allok=1, and tm gives the 
4946 c       collision time, t1 the collision time for i,
4947 c       t2 the collision time for j
4948 
4949         implicit double precision (a-h, o-z)
4950         parameter (MAXPTN=400001)
4951 
4952         common /para2/ xmp, xmu, alpha, rscut2, cutof2
4953 cc      SAVE /para2/
4954         common /para5/ iconfg, iordsc
4955 cc      SAVE /para5/
4956         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
4957      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
4958      &       xmass(MAXPTN), ityp(MAXPTN)
4959 cc      SAVE /prec2/
4960         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
4961 cc      SAVE /prec4/
4962         common /aurec1/ jxa, jya, jza
4963 cc      SAVE /aurec1/
4964         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
4965 cc      SAVE /aurec2/
4966         common /ilist1/
4967      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
4968      &     ictype, icsta(MAXPTN),
4969      &     nic(MAXPTN), icels(MAXPTN)
4970 cc      SAVE /ilist1/
4971         common /ilist3/ size1, size2, size3, v1, v2, v3, size
4972 cc      SAVE /ilist3/
4973         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
4974 cc      SAVE /ilist5/
4975         SAVE   
4976 
4977         logical allok, allokp
4978 
4979 c       preventing consecutive collisions
4980         allok = last(i) .ne. j .or. last(j) .ne. i
4981 
4982         if (ft(i) .ge. ft(j)) then
4983            i1 = j
4984            i2 = i
4985            isign = -1
4986         else 
4987            i1 = i
4988            i2 = j
4989            isign = 1
4990         end if
4991 
4992         if (allok) then
4993            tm = tlarge
4994            
4995            t1 = ft(i1)
4996            vx1 = vx(i1)
4997            vy1 = vy(i1)
4998            vz1 = vz(i1)
4999            
5000            t2 = ft(i2)
5001            
5002            dvx = vx(i2) - vx1
5003            dvy = vy(i2) - vy1
5004            dvz = vz(i2) - vz1
5005            
5006            dt = t2 - t1
5007 
5008            do 1000 ii = - 1, 1
5009               do 2000 jj = - 1, 1
5010 
5011                  allokp = .true.
5012 
5013                  dx = gx(i2) - gx(i1) - vx1 * dt
5014                  dy = gy(i2) - gy(i1) - vy1 * dt
5015                  dz = gz(i2) - gz(i1) - vz1 * dt
5016 
5017                  dx = dx + ii * 10d0 * size1
5018                  dy = dy + jj * 10d0 * size2
5019 
5020                  vp = dvx * dx + dvy * dy + dvz * dz
5021 
5022                  allokp = allokp .and. vp .lt. 0d0
5023                  
5024                  if (allokp) then
5025 
5026                     v2 = dvx * dvx + dvy * dvy + dvz * dvz
5027 
5028                     if (v2 .eq. 0d0) then
5029                        tmp = tlarge
5030                     else
5031                        tmp = t2 - vp / v2
5032                     end if
5033 
5034 c       note now tm is the absolute time
5035 
5036                     allokp = allokp .and. tmp .gt. t1 .and.
5037      &                         tmp .gt. t2
5038 
5039                  end if
5040 
5041                  if (allokp) then
5042 
5043                     dgx = dx - dvx * t2
5044                     dgy = dy - dvy * t2
5045                     dgz = dz - dvz * t2
5046 
5047                     dm2 = - v2 * tmp ** 2  + dgx * dgx +
5048      &                    dgy * dgy + dgz * dgz
5049 
5050                     allokp = allokp .and. dm2 .lt. cutof2
5051 
5052                  end if
5053 
5054                  if (allokp .and. tmp .lt. tm) then
5055                     tm = tmp
5056                     jxa = isign * ii
5057                     jya = isign * jj
5058                  end if
5059 
5060  2000              continue
5061  1000           continue
5062            
5063            if (tm .eq. tlarge) then
5064               allok = .false.
5065            end if
5066         end if
5067         
5068         if (allok) then
5069 
5070            e1 = e(i1)
5071            px1 = px(i1)
5072            py1 = py(i1)
5073            pz1 = pz(i1)
5074            e2 = e(i2)
5075            px2 = px(i2)
5076            py2 = py(i2)
5077            pz2 = pz(i2)
5078 
5079            rts2 = (e1 + e2) ** 2 - (px1 + px2) ** 2
5080      &          - (py1 + py2) ** 2 - (pz1 + pz2) ** 2
5081 
5082            allok = allok .and. rts2 .gt. rscut2
5083         end if
5084 
5085         if (.not. allok) then
5086            tm = tlarge
5087            t1 = tlarge
5088            t2 = tlarge
5089         else
5090            t1 = tm
5091            t2 = tm
5092         end if
5093 
5094         return
5095         end
5096 
5097         subroutine isco10(i, j, allok, tm, t1, t2)
5098 c       this subroutine is used to decide whether there is a collision between
5099 c       particle i and j, if there is one allok=1, and tm gives the 
5100 c       collision time, t1 the collision time for i,
5101 c       t2 the collision time for j
5102 
5103         implicit double precision (a-h, o-z)
5104         parameter (MAXPTN=400001)
5105 
5106         common /para2/ xmp, xmu, alpha, rscut2, cutof2
5107 cc      SAVE /para2/
5108         common /para5/ iconfg, iordsc
5109 cc      SAVE /para5/
5110         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
5111      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
5112      &       xmass(MAXPTN), ityp(MAXPTN)
5113 cc      SAVE /prec2/
5114         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
5115 cc      SAVE /prec4/
5116         common /aurec1/ jxa, jya, jza
5117 cc      SAVE /aurec1/
5118         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
5119 cc      SAVE /aurec2/
5120         common /ilist1/
5121      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
5122      &     ictype, icsta(MAXPTN),
5123      &     nic(MAXPTN), icels(MAXPTN)
5124 cc      SAVE /ilist1/
5125         common /ilist3/ size1, size2, size3, v1, v2, v3, size
5126 cc      SAVE /ilist3/
5127         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
5128 cc      SAVE /ilist5/
5129         SAVE   
5130 
5131         logical allok, allokp
5132 
5133 c       preventing consecutive collisions
5134         allok = last(i) .ne. j .or. last(j) .ne. i
5135 
5136 c       set up numbers for later calculations
5137 
5138         tm = tlarge
5139 
5140         if (allok) then
5141            do 1000 ii = - 1, 1
5142               do 2000 jj = - 1, 1
5143                  do 3000 kk = -1, 1
5144                  allokp = .true.
5145                  
5146                  i1 = i
5147                  i2 = j
5148 
5149                  p4 = ft(j) - ft(i)
5150                  p1 = gx(j) - gx(i)
5151                  p2 = gy(j) - gy(i)
5152                  p3 = gz(j) - gz(i)
5153 
5154                  p1 = p1 + ii * 10d0 * size1
5155                  p2 = p2 + jj * 10d0 * size2
5156                  p3 = p3 + kk * 10d0 * size3
5157 
5158                  q4 = e(i)
5159                  q1 = px(i)
5160                  q2 = py(i)
5161                  q3 = pz(i)
5162                  
5163                  r4 = e(j)
5164                  r1 = px(j)
5165                  r2 = py(j)
5166                  r3 = pz(j)
5167 
5168                  a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
5169                  b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
5170                  c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
5171                  d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
5172                  ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
5173                  f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
5174 
5175 c       make sure particle 2 formed early
5176                  h = a + b
5177                  if (h .gt. 0d0) then
5178                     g = a
5179                     a = -b
5180                     b = -g
5181                     g = c
5182                     c = d
5183                     d = g
5184                     i1 = j
5185                     i2 = i
5186                  end if
5187                  
5188 c       check the approaching criteria
5189                  if (allokp) then
5190                     vp = a * d - b * ee
5191                     allokp = allokp .and. vp .lt. 0d0
5192                  end if
5193 
5194 c       check the closest approach distance criteria
5195                  if (allokp) then
5196            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
5197      &            (ee ** 2 - c * d)
5198                     allokp = allokp .and. dm2 .lt. cutof2
5199                  end if
5200 
5201 c       check the time criteria
5202                  if (allokp) then
5203            tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
5204            tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
5205            tmp = 0.5d0 * (tc1 + tc2)
5206            allokp = allokp .and. tmp .gt. ft(i) .and. tmp .gt. ft(j)
5207                  end if
5208 
5209                  if (allokp .and. tmp .lt. tm) then
5210                     tm = tmp
5211                     jxa = ii
5212                     jya = jj
5213                     jza = kk
5214 cd                    dgxa(j) = ii * 10d0 * size1
5215 cd                    dgya(j) = jj * 10d0 * size2
5216 cd                    dgxa(i) = - dgxa(j)
5217 cd                    dgya(i) = - dgya(j)
5218                  end if
5219 
5220  3000                 continue
5221  2000              continue
5222  1000           continue
5223 
5224            if (tm .eq. tlarge) then
5225               allok = .false.
5226            end if
5227            
5228         end if
5229 
5230 c        check rts cut
5231         if (allok) then
5232 
5233            q4 = e(i1)
5234            q1 = px(i1)
5235            q2 = py(i1)
5236            q3 = pz(i1)
5237 
5238            r4 = e(i2)
5239            r1 = px(i2)
5240            r2 = py(i2)
5241            r3 = pz(i2)
5242 
5243            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
5244      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
5245 
5246            allok = allok .and. rts2 .gt. rscut2
5247         end if
5248           
5249         if (.not. allok) then
5250            tm = tlarge
5251            t1 = tlarge
5252            t2 = tlarge
5253         else if (h .gt. 0d0) then
5254            t1 = tm
5255            t2 = tm
5256         else
5257            t1 = tm
5258            t2 = tm
5259         end if
5260 
5261         return
5262         end
5263 
5264         subroutine isco11(i, j, allok, tm, t1, t2)
5265 c       this subroutine is used to decide whether there is a collision between
5266 c       particle i and j, if there is one allok=1, and tm gives the 
5267 c       collision time, t1 the collision time for i,
5268 c       t2 the collision time for j
5269 
5270         implicit double precision (a-h, o-z)
5271         parameter (MAXPTN=400001)
5272 
5273         common /para2/ xmp, xmu, alpha, rscut2, cutof2
5274 cc      SAVE /para2/
5275         common /para5/ iconfg, iordsc
5276 cc      SAVE /para5/
5277         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
5278      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
5279      &       xmass(MAXPTN), ityp(MAXPTN)
5280 cc      SAVE /prec2/
5281         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
5282 cc      SAVE /prec4/
5283         common /aurec1/ jxa, jya, jza
5284 cc      SAVE /aurec1/
5285         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
5286 cc      SAVE /aurec2/
5287         common /ilist1/
5288      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
5289      &     ictype, icsta(MAXPTN),
5290      &     nic(MAXPTN), icels(MAXPTN)
5291 cc      SAVE /ilist1/
5292         common /ilist3/ size1, size2, size3, v1, v2, v3, size
5293 cc      SAVE /ilist3/
5294         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
5295 cc      SAVE /ilist5/
5296         SAVE   
5297 
5298         logical allok, allokp
5299 
5300 c       preventing consecutive collisions
5301         allok = last(i) .ne. j .or. last(j) .ne. i
5302 
5303 c       set up numbers for later calculations
5304 
5305         tm = tlarge
5306 
5307         if (allok) then
5308            do 1000 ii = - 1, 1
5309               do 2000 jj = - 1, 1
5310                  do 3000 kk = - 1, 1
5311 
5312                  allokp = .true.
5313                  
5314                  i1 = i
5315                  i2 = j
5316 
5317                  p4 = ft(j) - ft(i)
5318                  p1 = gx(j) - gx(i)
5319                  p2 = gy(j) - gy(i)
5320                  p3 = gz(j) - gz(i)
5321 
5322                  p1 = p1 + ii * 10d0 * size1
5323                  p2 = p2 + jj * 10d0 * size2
5324                  p3 = p3 + kk * 10d0 * size3
5325 
5326                  q4 = e(i)
5327                  q1 = px(i)
5328                  q2 = py(i)
5329                  q3 = pz(i)
5330                  
5331                  r4 = e(j)
5332                  r1 = px(j)
5333                  r2 = py(j)
5334                  r3 = pz(j)
5335 
5336                  a = p4 * q4 - p1 * q1 - p2 * q2 - p3 * q3
5337                  b = p4 * r4 - p1 * r1 - p2 * r2 - p3 * r3
5338                  c = q4 * q4 - q1 * q1 - q2 * q2 - q3 * q3
5339                  d = r4 * r4 - r1 * r1 - r2 * r2 - r3 * r3
5340                  ee = q4 * r4 - q1 * r1 - q2 * r2 - q3 * r3
5341                  f = p4 * p4 - p1 * p1 - p2 * p2 - p3 * p3
5342 
5343 c       make sure particle 2 formed early
5344                  h = a + b
5345                  if (h .gt. 0d0) then
5346                     g = a
5347                     a = -b
5348                     b = -g
5349                     g = c
5350                     c = d
5351                     d = g
5352                     i1 = j
5353                     i2 = i
5354                  end if
5355                  
5356 c       check the approaching criteria
5357                  if (allokp) then
5358                     vp = a * d - b * ee
5359                     allokp = allokp .and. vp .lt. 0d0
5360                  end if
5361 
5362 c       check the closest approach distance criteria
5363                  if (allokp) then
5364            dm2 = - f - (a ** 2 * d + b ** 2 * c - 2d0 * a * b * ee) /
5365      &            (ee ** 2 - c * d)
5366                     allokp = allokp .and. dm2 .lt. cutof2
5367                  end if
5368 
5369 c       check the time criteria
5370                  if (allokp) then
5371            tc1 = ft(i1) - e(i1) * (a * d - b * ee) / (ee ** 2 - c * d)
5372            tc2 = ft(i2) + e(i2) * (b * c - a * ee) / (ee ** 2 - c * d)
5373                     if (iordsc .eq. 20) then
5374                        tmp = min(tc1, tc2)
5375                     else if (iordsc .eq. 21) then
5376                        tmp = 0.5d0 * (tc1 + tc2)
5377                     else
5378                        tmp = max(tc1, tc2)
5379                     end if
5380            allokp = allokp .and. tmp .gt. ft(i) .and. tmp .gt. ft(j)
5381                  end if
5382 
5383                  if (allokp .and. tmp .lt. tm) then
5384                     tm = tmp
5385                     jxa = ii
5386                     jya = jj
5387                     jza = kk
5388                     ha = h
5389                     tc1a = tc1
5390                     tc2a = tc2
5391 cd                    dgxa(j) = ii * 10d0 * size1
5392 cd                    dgya(j) = jj * 10d0 * size2
5393 cd                    dgxa(i) = - dgxa(j)
5394 cd                    dgya(i) = - dgya(j)
5395                  end if
5396 
5397  3000                 continue
5398  2000              continue
5399  1000           continue
5400 
5401            if (tm .eq. tlarge) then
5402               allok = .false.
5403            end if
5404            
5405         end if
5406 
5407 c        check rts cut
5408         if (allok) then
5409 
5410            q4 = e(i1)
5411            q1 = px(i1)
5412            q2 = py(i1)
5413            q3 = pz(i1)
5414 
5415            r4 = e(i2)
5416            r1 = px(i2)
5417            r2 = py(i2)
5418            r3 = pz(i2)
5419 
5420            rts2 = (q4 + r4) ** 2 - (q1 + r1) ** 2
5421      &          - (q2 + r2) ** 2 - (q3 + r3) ** 2
5422 
5423            allok = allok .and. rts2 .gt. rscut2
5424         end if
5425 
5426         if (.not. allok) then
5427            tm = tlarge
5428            t1 = tlarge
5429            t2 = tlarge
5430         else if (ha .gt. 0d0) then
5431            t1 = tc2a
5432            t2 = tc1a
5433         else
5434            t1 = tc1a
5435            t2 = tc2a
5436         end if
5437 
5438         return
5439         end
5440 
5441         subroutine isco12(i, j, allok, tm, t1, t2)
5442 c       this subroutine is used to decide whether there is a collision between
5443 c       particle i and j, if there is one allok=1, and tm gives the 
5444 c       collision time, t1 the collision time for i,
5445 c       t2 the collision time for j
5446 
5447         implicit double precision (a-h, o-z)
5448         parameter (MAXPTN=400001)
5449 
5450         common /para2/ xmp, xmu, alpha, rscut2, cutof2
5451 cc      SAVE /para2/
5452         common /para5/ iconfg, iordsc
5453 cc      SAVE /para5/
5454         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
5455      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
5456      &       xmass(MAXPTN), ityp(MAXPTN)
5457 cc      SAVE /prec2/
5458         common /prec4/ vx(MAXPTN), vy(MAXPTN), vz(MAXPTN)
5459 cc      SAVE /prec4/
5460         common /aurec1/ jxa, jya, jza
5461 cc      SAVE /aurec1/
5462         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
5463 cc      SAVE /aurec2/
5464         common /ilist1/
5465      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
5466      &     ictype, icsta(MAXPTN),
5467      &     nic(MAXPTN), icels(MAXPTN)
5468 cc      SAVE /ilist1/
5469         common /ilist3/ size1, size2, size3, v1, v2, v3, size
5470 cc      SAVE /ilist3/
5471         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
5472 cc      SAVE /ilist5/
5473         SAVE   
5474 
5475         logical allok, allokp
5476 
5477 c       preventing consecutive collisions
5478         allok = last(i) .ne. j .or. last(j) .ne. i
5479 
5480         if (ft(i) .ge. ft(j)) then
5481            i1 = j
5482            i2 = i
5483            isign = -1
5484         else 
5485            i1 = i
5486            i2 = j
5487            isign = 1
5488         end if
5489 
5490         if (allok) then
5491            tm = tlarge
5492            
5493            t1 = ft(i1)
5494            vx1 = vx(i1)
5495            vy1 = vy(i1)
5496            vz1 = vz(i1)
5497            
5498            t2 = ft(i2)
5499            
5500            dvx = vx(i2) - vx1
5501            dvy = vy(i2) - vy1
5502            dvz = vz(i2) - vz1
5503            
5504            dt = t2 - t1
5505 
5506            do 1000 ii = - 1, 1
5507               do 2000 jj = - 1, 1
5508                  do 3000 kk = -1, 1
5509 
5510                  allokp = .true.
5511 
5512                  dx = gx(i2) - gx(i1) - vx1 * dt
5513                  dy = gy(i2) - gy(i1) - vy1 * dt
5514                  dz = gz(i2) - gz(i1) - vz1 * dt
5515 
5516                  dx = dx + ii * 10d0 * size1
5517                  dy = dy + jj * 10d0 * size2
5518                  dz = dz + kk * 10d0 * size3
5519 
5520                  vp = dvx * dx + dvy * dy + dvz * dz
5521 
5522                  allokp = allokp .and. vp .lt. 0d0
5523                  
5524                  if (allokp) then
5525 
5526                     v2 = dvx * dvx + dvy * dvy + dvz * dvz
5527 
5528                     if (v2 .eq. 0d0) then
5529                        tmp = tlarge
5530                     else
5531                        tmp = t2 - vp / v2
5532                     end if
5533 
5534 c       note now tm is the absolute time
5535 
5536                     allokp = allokp .and. tmp .gt. t1 .and.
5537      &                         tmp .gt. t2
5538 
5539                  end if
5540 
5541                  if (allokp) then
5542 
5543                     dgx = dx - dvx * t2
5544                     dgy = dy - dvy * t2
5545                     dgz = dz - dvz * t2
5546 
5547                     dm2 = - v2 * tmp ** 2  + dgx * dgx +
5548      &                    dgy * dgy + dgz * dgz
5549 
5550                     allokp = allokp .and. dm2 .lt. cutof2
5551 
5552                  end if
5553 
5554                  if (allokp .and. tmp .lt. tm) then
5555                     tm = tmp
5556                     jxa = isign * ii
5557                     jya = isign * jj
5558                     jza = isign * kk
5559                  end if
5560 
5561  3000                 continue
5562  2000              continue
5563  1000           continue
5564            
5565            if (tm .eq. tlarge) then
5566               allok = .false.
5567            end if
5568         end if
5569         
5570         if (allok) then
5571 
5572            e1 = e(i1)
5573            px1 = px(i1)
5574            py1 = py(i1)
5575            pz1 = pz(i1)
5576            e2 = e(i2)
5577            px2 = px(i2)
5578            py2 = py(i2)
5579            pz2 = pz(i2)
5580 
5581            rts2 = (e1 + e2) ** 2 - (px1 + px2) ** 2
5582      &          - (py1 + py2) ** 2 - (pz1 + pz2) ** 2
5583 
5584            allok = allok .and. rts2 .gt. rscut2
5585         end if
5586 
5587         if (.not. allok) then
5588            tm = tlarge
5589            t1 = tlarge
5590            t2 = tlarge
5591         else
5592            t1 = tm
5593            t2 = tm
5594         end if
5595 
5596         return
5597         end
5598 
5599         subroutine reor(t, tmin, j, last0)
5600 c       this subroutine is used to fix ct(i) when tm is greater than ct(i)
5601 c       next(i) is last1 or last2
5602 
5603         implicit double precision (a-h, o-z)
5604         parameter (MAXPTN=400001)
5605 
5606         common /para5/ iconfg, iordsc
5607 cc      SAVE /para5/
5608         common /ilist1/
5609      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
5610      &     ictype, icsta(MAXPTN),
5611      &     nic(MAXPTN), icels(MAXPTN)
5612 cc      SAVE /ilist1/
5613 cd        common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
5614 cc      SAVE /ilist5/
5615         SAVE   
5616 
5617         icels0 = icels(j)
5618 
5619         i1 = icels0 / 10000
5620         i2 = (icels0 - i1 * 10000) / 100
5621         i3 = icels0 - i1 * 10000 - i2 * 100
5622 
5623         call wallc(j, i1, i2, i3, t, tmin1)
5624 
5625         if (tmin .le. tmin1) then
5626            nc = last0
5627         else
5628            tmin = tmin1
5629            nc = 0
5630         end if
5631 
5632         if (iconfg .eq. 3 .or. iconfg .eq. 5) then
5633            call chcell(j, i1, i2, i3, last0, t, tmin, nc)
5634         else
5635            if (i1 .eq. 11 .and. i2 .eq. 11 .and. i3 .eq. 11) then
5636               call chout(j, last0, t, tmin, nc)
5637            else
5638               if (iconfg .eq. 1) then
5639                  call chin1(j, i1, i2, i3, last0, t, tmin, nc)
5640               else if (iconfg .eq. 2) then
5641                  call chin2(j, i1, i2, i3, last0, t, tmin, nc)
5642               else if (iconfg .eq. 4) then
5643                  call chin3(j, i1, i2, i3, last0, t, tmin, nc)
5644               end if
5645            end if
5646         end if
5647         
5648         call fixtim(j, t, tmin1, tmin, nc)
5649 
5650         return
5651         end
5652 
5653         subroutine chout(l, last0, t, tmin, nc)
5654 c       this subroutine is used to check the surface when the colliding
5655 c       particle is outside the cube
5656 
5657         implicit double precision (a-h, o-z)
5658         parameter (MAXPTN=400001)
5659 
5660         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
5661      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
5662      &       xmass(MAXPTN), ityp(MAXPTN)
5663 cc      SAVE /prec2/
5664         SAVE   
5665 
5666         m1 = 11
5667         m2 = 11
5668         m3 = 11
5669         call chcell(l, m1, m2, m3, last0, t, tmin, nc)
5670 
5671         do 1003 i = 1, 10
5672            do 1002 j = 1, 10
5673               do 1001 k = 1, 10
5674                  if (i .eq. 1 .or. i .eq. 10 .or. j .eq. 1 .or.
5675      &              j .eq. 10 .or. k .eq. 1 .or. k. eq. 10)
5676      &               call chcell(l, i, j, k, last0, t, tmin, nc)
5677  1001         continue
5678  1002      continue
5679  1003   continue
5680 
5681         return
5682         end
5683 
5684         subroutine chin1(l, i1, i2, i3, last0, t, tmin, nc)
5685 c       this subroutine is used to check collisions for particle inside
5686 c       the cube
5687 c       and update the afftected particles through chcell
5688 
5689         implicit double precision (a-h, o-z)
5690         SAVE   
5691         
5692 c       itest is a flag to make sure the 111111 cell is checked only once
5693         itest = 0
5694 
5695         do 1003 i = i1 - 1, i1 + 1
5696            do 1002 j = i2 - 1, i2 + 1
5697               do 1001 k = i3 - 1, i3 + 1
5698                  if (i .ge. 1 .and. i .le. 10
5699      &              .and. j .ge. 1 .and. j .le. 10
5700      &              .and. k .ge. 1 .and. k .le. 10) then
5701                     call chcell(l, i, j, k, last0, t, tmin, nc)
5702                  else if (itest .eq. 0) then
5703                     m1 = 11
5704                     m2 = 11
5705                     m3 = 11
5706                     call chcell(l, m1, m2, m3, last0, t, tmin, nc)
5707                     itest = 1
5708                  end if   
5709  1001         continue
5710  1002      continue
5711  1003   continue
5712 
5713         return
5714         end
5715 
5716         subroutine chin2(l, i1, i2, i3, last0, t, tmin, nc)
5717 c       this subroutine is used to check collisions for particle inside
5718 c       the cube
5719 c       and update the afftected particles through chcell
5720 
5721         implicit double precision (a-h, o-z)
5722         SAVE   
5723         
5724 c       itest is a flag to make sure the 111111 cell is checked only once
5725         itest = 0
5726 
5727         do 1003 i = i1 - 1, i1 + 1
5728            do 1002 j = i2 - 1, i2 + 1
5729               do 1001 k = i3 - 1, i3 + 1
5730                  ia = i
5731                  ib = j
5732                  ic = k
5733                  if (k .ge. 1 .and. k .le. 10) then
5734                     if (i .eq. 0) ia = 10
5735                     if (i .eq. 11) ia = 1
5736                     if (j .eq. 0) ib = 10
5737                     if (j .eq. 11) ib = 1
5738                     call chcell(l, ia, ib, ic, last0, t, tmin, nc)
5739                  end if
5740  1001         continue
5741  1002      continue
5742  1003   continue
5743 
5744         return
5745         end
5746 
5747         subroutine chin3(l, i1, i2, i3, last0, t, tmin, nc)
5748 c       this subroutine is used to check collisions for particle inside
5749 c       the cube
5750 c       and update the afftected particles through chcell
5751 
5752         implicit double precision (a-h, o-z)
5753         SAVE   
5754         
5755 c       itest is a flag to make sure the 111111 cell is checked only once
5756         itest = 0
5757 
5758         do 1003 i = i1 - 1, i1 + 1
5759            do 1002 j = i2 - 1, i2 + 1
5760               do 1001 k = i3 - 1, i3 + 1
5761                  if (i .eq. 0) then
5762                     ia = 10
5763                  else if (i .eq. 11) then
5764                     ia = 1
5765                  else
5766                     ia = i
5767                  end if
5768                  if (j .eq. 0) then
5769                     ib = 10
5770                  else if (j .eq. 11) then
5771                     ib = 1
5772                  else
5773                     ib = j
5774                  end if
5775                  if (k .eq. 0) then
5776                     ic = 10
5777                  else if (k .eq. 11) then
5778                     ic = 1
5779                  else
5780                     ic = k
5781                  end if
5782                  call chcell(l, ia, ib, ic, last0, t, tmin, nc)
5783  1001         continue
5784  1002      continue
5785  1003   continue
5786 
5787         return
5788         end
5789 
5790         subroutine chcell(il, i1, i2, i3, last0, t, tmin, nc)
5791 c       this program is used to check through all the particles, except last0
5792 c       in the cell (i1,i2,i3) and see if we can get a particle collision 
5793 c       with time less than the original input tmin ( the collision time of 
5794 c       il with the wall
5795 c       last0 cas be set to 0 if we don't want to exclude last0
5796 
5797         implicit double precision (a-h, o-z)
5798         parameter (MAXPTN=400001)
5799         
5800         common /para5/ iconfg, iordsc
5801 cc      SAVE /para5/
5802         common /ilist1/
5803      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
5804      &     ictype, icsta(MAXPTN),
5805      &     nic(MAXPTN), icels(MAXPTN)
5806 cc      SAVE /ilist1/
5807         common /ilist2/ icell, icel(10, 10, 10)
5808 cc      SAVE /ilist2/
5809         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
5810 cc      SAVE /ilist4/
5811         SAVE   
5812 
5813         if (iconfg .eq. 3 .or. iconfg .eq. 5) then
5814            jj = ichkpt
5815            do 1001 j = 1, jj
5816 c     10/24/02 get rid of argument usage mismatch in mintm():
5817               jmintm=j
5818               if (j .ne. il .and. j .ne. last0)
5819      &          call mintm(il, jmintm, tmin, nc)
5820 c     &          call mintm(il, j, tmin, nc)
5821 
5822  1001         continue
5823            return
5824         end if
5825 
5826 c       set l
5827         if (i1 .eq. 11 .and. i2 .eq. 11 .and. i3 .eq. 11) then
5828            l = icell
5829         else
5830            l = icel(i1 ,i2, i3)
5831         end if
5832 
5833         if (l .eq. 0) return
5834         
5835         j = nic(l)
5836         
5837 c       if there is only one particle
5838         if (j .eq. 0) then
5839            
5840 c       if it's not il or last0,when last is not wall
5841            if (l .eq. il .or. l .eq. last0) return
5842            call mintm(il, l, tmin, nc)
5843            
5844 c       if there are many particles
5845         else
5846            if (l .ne. il .and. l .ne. last0)
5847      &        call mintm(il, l, tmin, nc)
5848            do 10 while(j .ne. l)
5849               if (j .ne. il .and. j .ne. last0)
5850      &             call mintm(il, j, tmin, nc)
5851               j = nic(j)
5852  10           continue
5853         end if
5854         
5855         return
5856         end
5857 
5858         subroutine mintm(i, j, tmin, nc)
5859 c       this subroutine is used to check whether particle j has smaller
5860 c       collision time with particle i than other particles
5861 c       or in other words, update next(i)
5862 
5863 c       input i,j,tmin,nc
5864 c       output tmin,nc
5865 
5866         implicit double precision (a-h, o-z)
5867         parameter (MAXPTN=400001)
5868 
5869         common /para5/ iconfg, iordsc
5870 cc      SAVE /para5/
5871         common /aurec1/ jxa, jya, jza
5872 cc      SAVE /aurec1/
5873         common /aurec2/ dgxa(MAXPTN), dgya(MAXPTN), dgza(MAXPTN)
5874 cc      SAVE /aurec2/
5875         common /ilist1/
5876      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
5877      &     ictype, icsta(MAXPTN),
5878      &     nic(MAXPTN), icels(MAXPTN)
5879 cc      SAVE /ilist1/
5880         common /ilist3/ size1, size2, size3, v1, v2, v3, size
5881 cc      SAVE /ilist3/
5882         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
5883 cc      SAVE /ilist5/
5884         SAVE   
5885 
5886         logical allok
5887 
5888         call isco(i, j, allok, tm, t1, t2)
5889 
5890         if (allok .and. tm .lt. tmin) then
5891            tmin = tm
5892            ct(i) = t1
5893            nc = j
5894            if (iconfg .eq. 3 .or. iconfg .eq. 5) then
5895               dgxa(i) = - jxa * 10d0 * size1
5896               dgya(i) = - jya * 10d0 * size2
5897               if (iconfg .eq. 5) then
5898                  dgza(i) = - jza * 10d0 * size3
5899               end if
5900            end if
5901         end if
5902 
5903          return
5904         end
5905 
5906 ******************************************************************************
5907 ******************************************************************************
5908 
5909         subroutine zpca1
5910 
5911         implicit double precision (a-h, o-z)
5912         parameter (MAXPTN=400001)
5913 
5914         common /ilist1/
5915      &     iscat, jscat, next(MAXPTN), last(MAXPTN),
5916      &     ictype, icsta(MAXPTN),
5917      &     nic(MAXPTN), icels(MAXPTN)
5918 cc      SAVE /ilist1/
5919         SAVE   
5920 
5921         if (mod(ictype,2) .eq. 0) then
5922            call zpca1a(iscat)
5923            call zpca1a(jscat)
5924 clin-5/2009 ctest off v2 for parton:
5925 c           call flowp(1)
5926         end if
5927 
5928         return
5929         end
5930 
5931         subroutine zpca1a(i)
5932 
5933         implicit double precision (a-h, o-z)
5934         parameter (MAXPTN=400001)
5935 
5936         common /para2/ xmp, xmu, alpha, rscut2, cutof2
5937 cc      SAVE /para2/
5938         common /para5/ iconfg, iordsc
5939 cc      SAVE /para5/
5940         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
5941      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
5942      &       xmass(MAXPTN), ityp(MAXPTN)
5943 cc      SAVE /prec2/
5944         common /prec3/gxs(MAXPTN),gys(MAXPTN),gzs(MAXPTN),fts(MAXPTN),
5945      &     pxs(MAXPTN), pys(MAXPTN), pzs(MAXPTN), es(MAXPTN),
5946      &     xmasss(MAXPTN), ityps(MAXPTN)
5947 cc      SAVE /prec3/
5948         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
5949 cc      SAVE /prec5/
5950         common /prec6/ etas(MAXPTN), raps(MAXPTN), taus(MAXPTN)
5951 cc      SAVE /prec6/
5952         common /ana1/ ts(12)
5953 cc      SAVE /ana1/
5954         SAVE   
5955 
5956         if (iconfg .eq. 1) then
5957            t1 = fts(i)
5958            t2 = ft(i)
5959            ipic = 11
5960         else if (iconfg .eq. 2 .or.
5961      &     iconfg .eq. 3) then
5962 cd           t1 = fts(i)
5963 cd           t2 = ft(i)
5964            t1 = taus(i)
5965            t2 = tau(i)
5966            ipic = 12
5967         else if (iconfg .eq. 4 .or.
5968      &     iconfg .eq. 5) then
5969            t1 = fts(i)
5970            t2 = ft(i)
5971            ipic = 12
5972         end if
5973 
5974         if (iconfg .le. 3) then
5975            do 1002 ian = 1, ipic
5976               if (t1 .le. ts(ian) .and.
5977      &           t2 .gt. ts(ian)) then
5978                  rapi = raps(i)
5979 c     7/20/01:
5980 c                 et = sqrt(pxs(i) ** 2 + pys(i) ** 2 + xmp ** 2)
5981                  et = dsqrt(pxs(i) ** 2 + pys(i) ** 2 + xmp ** 2)
5982                  call zpca1b(rapi, et, ian)
5983               end if
5984  1002      continue
5985         else
5986            do 1003 ian = 1, ipic
5987               if (t1 .le. ts(ian) .and.
5988      &           t2 .gt. ts(ian)) then
5989                  p0 = es(i)
5990                  p1 = pxs(i)
5991                  p2 = pys(i)
5992                  p3 = pzs(i)
5993                  call zpca1c(p0, p1, p2, p3, ian)
5994               end if
5995  1003      continue
5996         end if
5997 
5998         return
5999         end
6000 
6001         subroutine zpca1b(rapi, et, ian)
6002 
6003         implicit double precision (a-h, o-z)
6004         parameter (MAXPTN=400001)
6005 
6006         common /para6/ centy
6007 cc      SAVE /para6/
6008         common /ilist6/ t, iopern, icolln
6009 cc      SAVE /ilist6/
6010         common /ana2/
6011      &     det(12), dn(12), detdy(12), detdn(12), dndy(12),
6012      &     det1(12), dn1(12), detdy1(12), detdn1(12), dndy1(12),
6013      &     det2(12), dn2(12), detdy2(12), detdn2(12), dndy2(12)
6014 cc      SAVE /ana2/
6015         SAVE   
6016 
6017         if (rapi .gt. centy - 0.5d0 .and. 
6018      &     rapi .lt. centy + 0.5d0) then
6019            det2(ian) = det2(ian) + et
6020            dn2(ian) = dn2(ian) + 1d0
6021 cdtrans
6022            if (ian .eq. 10) then
6023 cd              write (10, *) t, det2(ian)
6024            end if
6025            if (ian .eq. 11) then
6026 cd              write (11, *) t, det2(ian)
6027            end if
6028            if (ian .eq. 12) then
6029 cd              write (12, *) t, det2(ian)
6030            end if
6031 cdtransend
6032            if (rapi .gt. centy - 0.25d0 .and. 
6033      &        rapi .lt. centy + 0.25d0) then
6034               det1(ian) = det1(ian) + et
6035               dn1(ian) = dn1(ian) + 1d0
6036               if (rapi .gt. centy - 0.1d0 .and.
6037      &           rapi .lt. centy + 0.1d0) then
6038                  det(ian) = det(ian) + et
6039                  dn(ian) = dn(ian) + 1d0
6040               end if
6041            end if
6042         end if
6043 
6044         return
6045         end
6046 
6047         subroutine zpca1c(p0, p1, p2, p3, ian)
6048 
6049         implicit double precision (a-h, o-z)
6050 
6051         common /ana3/ em(4, 4, 12)
6052 cc      SAVE /ana3/
6053 
6054         dimension en(4)
6055         SAVE   
6056 
6057         en(1) = p0
6058         en(2) = p1
6059         en(3) = p2
6060         en(4) = p3
6061 
6062         do 1002 i = 1, 4
6063            do 1001 j = 1, 4
6064               em(i, j, ian) = em(i, j, ian) + en(i) * en(j) / p0
6065  1001      continue
6066  1002   continue
6067 
6068         return
6069         end
6070 
6071 ******************************************************************************
6072 
6073         subroutine zpca2
6074 
6075         implicit double precision (a-h, o-z)
6076         parameter (MAXPTN=400001)
6077 
6078         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
6079 cc      SAVE /para3/
6080         common /para5/ iconfg, iordsc
6081 cc      SAVE /para5/
6082         common /para7/ ioscar,nsmbbbar,nsmmeson
6083 cc      SAVE /para7/
6084         common /ilist6/ t, iopern, icolln
6085 cc      SAVE /ilist6/
6086         common /rndm1/ number
6087 cc      SAVE /rndm1/
6088         common /rndm2/ iff
6089 cc      SAVE /rndm2/
6090         common /rndm3/ iseedp
6091 cc      SAVE /rndm3/
6092         COMMON /AREVT/ IAEVT, IARUN, MISS
6093 cc      SAVE /AREVT/
6094         SAVE   
6095 
6096         if (iconfg .le. 3) then
6097            call zpca2a
6098         else
6099            call zpca2b
6100         end if
6101 
6102         if (ioscar .eq. 1) then
6103            call zpca2c
6104         end if
6105 
6106 cbzdbg2/17/99
6107 c        write (25, *) 'Event', nsevt - 1 + ievt, 
6108 c    &         ', run', isbrun,
6109 c        WRITE (25, *) ' Event ', IAEVT, ', run ', IARUN,
6110 c     &     ',\n\t number of operations = ', iopern,
6111 c     &     ',\n\t number of collisions between particles = ', 
6112 c     &         icolln,
6113 c     &     ',\n\t freezeout time=', t,
6114 c     &     ',\n\t ending at the ', number, 'th random number',
6115 c     &     ',\n\t ending collision iff=', iff
6116         WRITE (25, *) ' Event ', IAEVT, ', run ', IARUN
6117         WRITE (25, *) '    number of operations = ', iopern
6118         WRITE (25, *) '    number of collisions between particles = ', 
6119      &       icolln
6120         WRITE (25, *) '    freezeout time=', t
6121         WRITE (25, *) '    ending at the ', number, 'th random number'
6122         WRITE (25, *) '    ending collision iff=', iff
6123 
6124         return
6125         end
6126 
6127         subroutine zpca2a
6128 
6129         implicit double precision (a-h, o-z)
6130         parameter (MAXPTN=400001)
6131 
6132         common /para1/ mul
6133 cc      SAVE /para1/
6134         common /para2/ xmp, xmu, alpha, rscut2, cutof2
6135 cc      SAVE /para2/
6136         common /para3/ nsevt, nevnt, nsbrun, ievt, isbrun
6137 cc      SAVE /para3/
6138         common /para5/ iconfg, iordsc
6139 cc      SAVE /para5/
6140         common /para6/ centy
6141 cc      SAVE /para6/
6142         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
6143      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
6144      &       xmass(MAXPTN), ityp(MAXPTN)
6145 cc      SAVE /prec2/
6146         common /prec5/ eta(MAXPTN), rap(MAXPTN), tau(MAXPTN)
6147 cc      SAVE /prec5/
6148         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
6149 cc      SAVE /ilist4/
6150         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
6151 cc      SAVE /ilist5/
6152         common /ilist6/ t, iopern, icolln
6153 cc      SAVE /ilist6/
6154         common /rndm1/ number
6155 cc      SAVE /rndm1/
6156         common /rndm2/ iff
6157 cc      SAVE /rndm2/
6158         common /rndm3/ iseedp
6159 cc      SAVE /rndm3/
6160         common /ana1/ ts(12)
6161 cc      SAVE /ana1/
6162         common /ana2/
6163      &     det(12), dn(12), detdy(12), detdn(12), dndy(12),
6164      &     det1(12), dn1(12), detdy1(12), detdn1(12), dndy1(12),
6165      &     det2(12), dn2(12), detdy2(12), detdn2(12), dndy2(12)
6166 cc      SAVE /ana2/
6167         common /ana4/ fdetdy(24), fdndy(24), fdndpt(12)
6168 cc      SAVE /ana4/
6169         SAVE   
6170 
6171         do 1004 i = 1, ichkpt
6172            rapi = rap(i)
6173 c     7/20/01:
6174 c           et = sqrt(px(i) ** 2 + py(i) ** 2 + xmp ** 2)
6175            et = dsqrt(px(i) ** 2 + py(i) ** 2 + xmp ** 2)
6176 
6177            do 1001 j = 1, 24
6178               if (rapi .gt. j + centy - 13d0 
6179      &           .and. rapi .lt. j  + centy - 12d0) then
6180                  fdetdy(j) = fdetdy(j) + et
6181                  fdndy(j) = fdndy(j) + 1d0
6182               end if
6183  1001      continue
6184 
6185            do 1002 j = 1, 12
6186               if (et .gt. 0.5d0 * (j - 1) .and.
6187      &           et .lt. 0.5d0 * j ) then
6188                  fdndpt(j) = fdndpt(j) + 1d0
6189               end if
6190  1002      continue
6191 
6192            if (iconfg .eq. 1) then
6193               t1 = ft(i)
6194               t2 = tlarge
6195               ipic = 11
6196            else
6197               t1 = tau(i)
6198               t2 = tlarge
6199               ipic = 12
6200            end if
6201 
6202            do 1003 ian = 1, ipic
6203               if (t1 .le. ts(ian) .and.
6204      &           t2 .gt. ts(ian)) then
6205                  call zpca1b(rapi, et, ian)
6206               end if
6207  1003      continue
6208 
6209            if (iconfg .eq. 1) then
6210               call zpca1b(rapi, et, 12)
6211            end if
6212  1004   continue
6213 
6214         do 1005 ian = 1, 12
6215            if (dn(ian) .eq. 0d0 .or. dn1(ian) .eq. 0d0 .or.
6216      &        dn2(ian) .eq. 0d0) then
6217 clin-9/2012 suppress output:
6218 c              print *, 'event=', ievt
6219 c              print *, 'dn(', ian, ')=', dn(ian), 'dn1(', ian,
6220 c     &           ')=', dn1(ian), 'dn2(', ian, ')=', dn2(ian)
6221            end if
6222            detdy(ian) = detdy(ian) + det(ian)
6223            if (dn(ian) .ne. 0) then
6224               detdn(ian) = detdn(ian) + det(ian) / dn(ian)
6225            end if
6226            dndy(ian) = dndy(ian) + dn(ian)
6227            detdy1(ian) = detdy1(ian) + det1(ian)
6228            if (dn1(ian) .ne. 0) then
6229               detdn1(ian) = detdn1(ian) + det1(ian) / dn1(ian)
6230            end if
6231            dndy1(ian) = dndy1(ian) + dn1(ian)
6232            detdy2(ian) = detdy2(ian) + det2(ian)
6233            if (dn2(ian) .ne. 0) then
6234               detdn2(ian) = detdn2(ian) + det2(ian) / dn2(ian)
6235            end if
6236            dndy2(ian) = dndy2(ian) + dn2(ian)
6237  1005   continue
6238 
6239         return
6240         end
6241 
6242         subroutine zpca2b
6243 
6244         implicit double precision (a-h, o-z)
6245         parameter (MAXPTN=400001)
6246 
6247         common /prec2/gx(MAXPTN),gy(MAXPTN),gz(MAXPTN),ft(MAXPTN),
6248      &       px(MAXPTN), py(MAXPTN), pz(MAXPTN), e(MAXPTN),
6249      &       xmass(MAXPTN), ityp(MAXPTN)
6250 cc      SAVE /prec2/
6251         common /ilist4/ ifmpt, ichkpt, indx(MAXPTN)
6252 cc      SAVE /ilist4/
6253         common /ilist5/ ct(MAXPTN), ot(MAXPTN), tlarge
6254 cc      SAVE /ilist5/
6255         common /ana1/ ts(12)
6256 cc      SAVE /ana1/
6257         SAVE   
6258 
6259         do 1002 i = 1, ichkpt
6260            t1 = ft(i)
6261            t2 = tlarge
6262            ipic = 12
6263 
6264            do 1001 ian = 1, ipic
6265               if (t1 .le. ts(ian) .and.
6266      &           t2 .gt. ts(ian)) then
6267                  p0 = e(i)
6268                  p1 = px(i)
6269                  p2 = py(i)
6270                  p3 = pz(i)
6271