The following changes are done:
[u/mrichter/AliRoot.git] / EPOS / epos167 / epos-tim-155.f
CommitLineData
9ef1c2d9 1c---------------------------------------------------------------------
60ed90b3 2 subroutine timann
9ef1c2d9 3c---------------------------------------------------------------------
4c electron-positron
5c---------------------------------------------------------------------
6 include 'epos.inc'
7 include 'epos.incsem'
60ed90b3 8 common/nfla/nfla
9ef1c2d9 9 parameter (ntim=1000)
10 common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim)
11 &,idaprt(2,ntim)
12 integer jcp(6,2),jcm(6,2)
13 dimension jorprt(ntim),p1(5)
14
15 call utpri('timann',ish,ishini,5)
16 egyevt=engy
17 qsqevt=engy**2
18
19 123 nprtj=1
20 en=engy
21 q20=0.
22 nptl=0
23 nptl=nptl+1 !electron
24 pptl(1,nptl)=0.
25 pptl(2,nptl)=0.
26 pptl(3,nptl)=sqrt(en**2/4-2.61121e-7)
27 pptl(4,nptl)=en/2.
28 pptl(5,nptl)=0.000511
29 idptl(nptl)=12
30 istptl(nptl)=1
31
32 nptl=nptl+1 !positron
33 pptl(1,nptl)=0.
34 pptl(2,nptl)=0.
35 pptl(3,nptl)=-sqrt(en**2/4-2.61121e-7)
36 pptl(4,nptl)=en/2.
37 pptl(5,nptl)=0.000511
38 idptl(nptl)=-12
39 istptl(nptl)=1
40
41 nptl=nptl+1 !virtual gamma
42 pptl(1,nptl)=0.
43 pptl(2,nptl)=0.
44 pptl(3,nptl)=0.
45 pptl(4,nptl)=en
46 pptl(5,nptl)=en
47 idptl(nptl)=10
48 iorptl(nptl)=nptl-2
49 jorptl(nptl)=nptl-1
50 istptl(nptl)=1
51
52 pprt(1,nprtj)=0.
53 pprt(2,nprtj)=0.
54 pprt(3,nprtj)=0.
60ed90b3 55 pprt(4,nprtj)=en
56 pprt(5,nprtj)=en
9ef1c2d9 57 idaprt(1,nprtj)=2
58 idaprt(2,nprtj)=3
59 if(q20.gt.0.)then
60 pprt(5,nprtj)=sqrt(q20)
61 else
62 pprt(5,nprtj)=en
63 endif
64 nfla=0
65 do i=1,naflav
66 call idmass(i,am)
67 if (2.*am.lt.en) nfla=i
68 enddo
69 if(itflav.eq.0)then
60ed90b3 70 s=engy**2
9ef1c2d9 71 dlz=2.4
72 amz=91.1885
73 al=1./real(137.035989d0)
74 gf=1.16639e-5
75 ak=sqrt(2.)*gf*amz**2/(16*pi*al)
76 chi1=ak*s*(s-amz**2)/((s-amz**2)**2+dlz**2*amz**2)
77 chi2=ak**2*s**2/((s-amz**2)**2+dlz**2*amz**2)
78 qe=-1.
79 ve=0.25-2.*qe*0.232
80 ae=-0.5
81 qf=2./3.
60ed90b3 82 vf=sign(.5,qf)-2.*qf*0.232
9ef1c2d9 83 af=sign(.5,qf)
84 dsmax1=
60ed90b3 85 $ 2.*(qf**2-2.*qf*ve*vf*chi1
9ef1c2d9 86 $ +(ae**2+ve**2)*(af**2+vf**2)*chi2)
87 $ + abs(-4*qf*ae*af*chi1+8*ae*ve*af*vf*chi2)
60ed90b3 88
9ef1c2d9 89 qf=-1./3.
60ed90b3 90 vf=sign(.5,qf)-2.*qf*0.232
9ef1c2d9 91 af=sign(.5,qf)
92 dsmax2=
60ed90b3 93 $ 2.*(qf**2-2.*qf*ve*vf*chi1
9ef1c2d9 94 $ +(ae**2+ve**2)*(af**2+vf**2)*chi2)
95 $ + abs(-4*qf*ae*af*chi1+8*ae*ve*af*vf*chi2)
96
60ed90b3 97 100 iq1=1+INT(nfla*rangen())
9ef1c2d9 98 call idchrg(iq1,qf)
99 ct=-1.+2.*rangen()
60ed90b3 100 vf=sign(.5,qf)-2.*qf*0.232
9ef1c2d9 101 af=sign(.5,qf)
102
103 dsigma=
60ed90b3 104 $ (1.+ct**2)*(qf**2-2.*qf*ve*vf*chi1
9ef1c2d9 105 $ +(ae**2+ve**2)*(af**2+vf**2)*chi2)
106 $ + ct*(-4*qf*ae*af*chi1+8*ae*ve*af*vf*chi2)
107 if(rangen().gt.dsigma/max(dsmax1,dsmax2)) goto 100
108 else
109 iq1=itflav
110 endif
111 if(rangen().lt.0.5)iq1=-iq1
60ed90b3 112
9ef1c2d9 113 nprtj=nprtj+1
114 idprt(nprtj)=iq1
115 pprt(4,nprtj)=en/2.
116 pprt(5,nprtj)=en
117 iorprt(nprtj)=1
60ed90b3 118
9ef1c2d9 119 nprtj=nprtj+1
120 idprt(nprtj)=-iq1
121 pprt(4,nprtj)=en/2.
122 pprt(5,nprtj)=en
123 iorprt(nprtj)=1
124
125 call timsho(2,3)
126
127 jorprt(1)=0 !!color-connection, no origin!!
128 jt=1
129 if(idprt(idaprt(1,1)).lt.0)jt=2
60ed90b3 130 do i=1,nprtj
9ef1c2d9 131 if(idaprt(1,i).ne.0) then
132 js=jt
133 if(idprt(i).lt.0.and.
134 & ((idprt(idaprt(2,i)).eq.9.and.jt.eq.1).or.
60ed90b3 135 & (idprt(idaprt(1,i)).eq.9.and.jt.eq.2)))then
9ef1c2d9 136 js=3-jt
137 elseif(idprt(i).gt.0.and.idprt(i).ne.9.and.
138 & ((idprt(idaprt(2,i)).eq.9.and.jt.eq.2).or.
60ed90b3 139 & (idprt(idaprt(1,i)).eq.9.and.jt.eq.1)))then
9ef1c2d9 140 js=3-jt
141 elseif(idprt(i).eq.9.and.idprt(idaprt(1,i)).ne.9.and.
142 & ((idprt(idaprt(1,i)).lt.0.and.jt.eq.2).or.
143 & (idprt(idaprt(1,i)).gt.0.and.jt.eq.1)))then
144 js=3-jt
145 endif
146 jorprt(idaprt(3-js,i))=jorprt(i)
147 jorprt(i)=idaprt(js,i)
148 jorprt(idaprt(js,i))=idaprt(3-js,i)
149 endif
150 enddo
60ed90b3 151
9ef1c2d9 152 if(ish.ge.5)then
60ed90b3 153 i=1
9ef1c2d9 154 do while(i.ne.0)
60ed90b3 155 if(idaprt(1,i) .eq. 0 ) then
9ef1c2d9 156 write(ifch,*) idprt(i)
157 write(ifch,*) '|'
158 endif
159 i=jorprt(i)
160 enddo
161 endif
162
163 iptl=nptl
60ed90b3 164 i=1
9ef1c2d9 165 do while(i.gt.0)
60ed90b3 166 if(idaprt(1,i) .eq. 0) then
9ef1c2d9 167 nptl=nptl+1
168 do j=1,5
169 pptl(j,nptl)=pprt(j,i)
170 enddo
171 idptl(nptl)=idprt(i)
172 istptl(nptl)=20
173 iorptl(nptl)=iptl
174 qsqptl(nptl)=en**2
175 endif
176 i=jorprt(i)
177 enddo
178 ifrptl(1,iptl)=iptl+1
179 ifrptl(2,iptl)=nptl
180
181 nk1=iptl+1
182 441 nk2=nk1+1
183 do while (idptl(nk2).eq.9)
184 nk2=nk2+1
185 enddo
186 do i=1,4
187 p1(i)=0.
188 do j=nk1,nk2
189 p1(i)=p1(i)+pptl(i,j)
190 enddo
191 enddo
192 p1(5)=sqrt(max(0.,p1(4)**2-p1(3)**2-p1(2)**2-p1(1)**2))
193 do i=1,2
194 do j=1,6
195 jcm(j,i)=0
196 jcp(j,i)=0
197 enddo
198 enddo
199 ii=1
200 if(idptl(nk1).lt.0)ii=2
201 jcp(abs(idptl(nk1)),ii)=1
202 jcm(abs(idptl(nk2)),3-ii)=1
203 amm= utamnx(jcp,jcm)
204 if(amm.gt.p1(5))goto 123
205 nk1=nk2+1
206 if(nk1.lt.nptl)goto 441
207
208
209 if(ish.ge.5)then
210 write(ifch,*)
211 do i=1,nprtj
212 write(ifch,98) i,(pprt(j,i),j=1,5),idprt(i)
213 & ,iorprt(i),idaprt(1,i),idaprt(2,i),jorprt(i)
214 enddo
215 write(ifch,*)
216 do i=1,nprtj
60ed90b3 217 if(pprt(5,i).eq.0.)
9ef1c2d9 218 & write(ifch,99) i,(pprt(j,i),j=1,5),idprt(i)
219 enddo
220 write(ifch,*)
221 write(ifch,*)
222 endif
60ed90b3 223
9ef1c2d9 224 99 format(i4,5g10.3,1i4)
225 98 format(i4,5g10.3,5i4)
60ed90b3 226
9ef1c2d9 227 call utprix('timann',ish,ishini,5)
228 return
229 end
60ed90b3 230
9ef1c2d9 231c---------------------------------------------------------------------
232 subroutine timsh1(q20,en,idfla)
233c---------------------------------------------------------------------
234 include 'epos.inc'
235 include 'epos.incsem'
236 common/nfla/nfla
237 parameter (ntim=1000)
238 common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim)
239 &,idaprt(2,ntim)
240 common/cprtx/nprtjx,pprtx(5,2)
241 nfla=naflav
242 nprtj=1
243 pprt(1,nprtj)=0.
244 pprt(2,nprtj)=0.
245 pprt(3,nprtj)=0.
60ed90b3 246 pprt(4,nprtj)=en
9ef1c2d9 247 pprt(5,nprtj)=sqrt(q20)
248 idprt(nprtj)=idfla
249 call timsho(1,0)
250 nprtjx=0
251 return
252 end
253
254
255c---------------------------------------------------------------------
60ed90b3 256 subroutine timsh2(q20,q21,en,idfla1,idfla2)
9ef1c2d9 257c---------------------------------------------------------------------
258 include 'epos.inc'
259 include 'epos.incsem'
260 common/nfla/nfla
261 parameter (ntim=1000)
262 common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim)
263 &,idaprt(2,ntim)
264 common/cprtx/nprtjx,pprtx(5,2)
265 nfla=naflav
266
267 nprtj=1
268 pprt(1,nprtj)=0.
269 pprt(2,nprtj)=0.
270 pprt(3,nprtj)=0.
60ed90b3 271 pprt(4,nprtj)=en/2.
9ef1c2d9 272 pprt(5,nprtj)=sqrt(q20)
273 idprt(nprtj)=idfla1
274
275 nprtj=2
276 pprt(1,nprtj)=0.
277 pprt(2,nprtj)=0.
278 pprt(3,nprtj)=0.
279 pprt(4,nprtj)=en/2.
280 pprt(5,nprtj)=sqrt(q21)
281 idprt(nprtj)=idfla2
282
283 call timsho(1,2)
284
285 nprtjx=1
286 pprtx(1,nprtjx)=0.
287 pprtx(2,nprtjx)=0.
60ed90b3 288 pprtx(3,nprtjx)=en/2.
289 pprtx(4,nprtjx)=en/2.
9ef1c2d9 290 pprtx(5,nprtjx)=0
291 nprtjx=2
292 pprtx(1,nprtjx)=0.
293 pprtx(2,nprtjx)=0.
60ed90b3 294 pprtx(3,nprtjx)=-en/2.
9ef1c2d9 295 pprtx(4,nprtjx)=en/2.
296 pprtx(5,nprtjx)=0
297
298 return
299 end
300
301c---------------------------------------------------------------------
302 subroutine timsho(j1,j2)
303c---------------------------------------------------------------------
304
305
306 include 'epos.inc'
307 include 'epos.incsem'
308 parameter (ntim=1000)
309 common/cprt/nprtj,pprt(5,ntim),idprt(ntim),iorprt(ntim)
310 &,idaprt(2,ntim)
311 dimension pz(ntim),id(2,ntim)
312 dimension ij(2),ee(2),amm2(-6:10)
313 logical first
314 common/nfla/nfla
315 call utpri('timsho',ish,ishini,5)
316 if(iappl.eq.6)then
317 do i=1,6
318 call idmass(i,am)
319 amm2(i)=am**2
320 amm2(-i)=am**2
321 enddo
322 else
323 do i=1,6
324 amm2(i)=0.
325 amm2(-i)=0.
326 enddo
327 endif
328 amm2(9)=0.
329 amm2(0)=0.
330 amm2(10)=0.
331 nn=nprtj
332 ij(1)=j1
333 ij(2)=j2
334 ii2=2
335 if(j2.eq.0)ii2=1
336
337 ii=1
338 first=.true.
339 n10=0
340 10 n10=n10+1
341 if(float(n10).gt.1e7)then
342 goto9999
343 endif
344 io=iorprt(ij(ii))
345 idfl=idprt(ij(ii))
346 q2start=pprt(5,ij(ii))**2
347 E=pprt(4,ij(ii))
348 if(ij(2).eq.j2.and.ii2.eq.2)E=pprt(4,ij(1))+pprt(4,ij(2))
349 zetamx=0.
350 if(ij(1).ne.j1)then
60ed90b3 351 zetamx=pprt(5,io)/pprt(4,io)/sqrt(pz(io)*(1.-pz(io)))
9ef1c2d9 352 endif
353
354c call timdev(idfl,q2start,E,zetamx,idfla,idflb,qa2,z)
355c......................................................................
356 380 q2=q2start
357 z=0.5
358 PT2MIN=max(qcdlam*1.1,q2fin)
60ed90b3 359 ALFM=LOG(PT2MIN/qcdlam)
9ef1c2d9 360 if(ish.ge.9)then
361 write (ifch,*) '---------------------',ii
362 $ ,pprt(5,ij(1)) !,pprt(5,ij(2))
363 write(ifch,*) ' idfl,q2start,zetamx:',idfl,q2start,zetamx
364 endif
365 if (q2.lt.4.*q2fin+2.*amm2(idfl) )then
366 if(ish.ge.9)then
367 write(ifch,'(a,i4,i4,2f15.8)') 'null:',0.
368 endif
369 q2=amm2(idfl)
370 idfla=0
371 idflb=0
372 goto 999
373 endif
60ed90b3 374
9ef1c2d9 375 390 zc=.5*(1.-sqrt(max(0.000001,1.-4.*q2fin/q2)))
376 if(ish.ge.9)then
377 write(ifch,*) 'zc=',zc
378 endif
60ed90b3 379 IF(idfl.EQ.9) THEN
380 FBR=6.*LOG((1.-ZC)/ZC)+nfla*(0.5-ZC)
9ef1c2d9 381 ELSE
60ed90b3 382 FBR=(8./3.)*LOG((1.-ZC)/ZC)
9ef1c2d9 383 endif
60ed90b3 384
9ef1c2d9 385 B0=(33.-2.*nfla)/6.
60ed90b3 386
387
9ef1c2d9 388c.....select new q2
389 r=rangen()
390 q2=q2*exp(log(r)*B0*ALFM/FBR)
391 if(ish.ge.9)then
392 write(ifch,*) 'q^2=',q2
393 endif
394 if (q2.lt.4.*q2fin+2.*amm2(idfl))then
395 q2=amm2(idfl)
396 if(ish.ge.9)then
397 write(ifch,'(a,i4,i4,2f15.8)') 'null:',0.
398 endif
399 idfla=0
400 idflb=0
401 goto 999
402 endif
60ed90b3 403
9ef1c2d9 404c.....select flavor and z-value .....................................
60ed90b3 405 IF(idfl.EQ.9) THEN
9ef1c2d9 406 if(rangen()*FBR.lt.nfla*(0.5-ZC))then
407 ! .................g -> qqbar
408 Z=ZC+(1.-2.*ZC)*rangen()
60ed90b3 409 IF(Z**2+(1.-Z)**2.LT.rangen()) GOTO 390
9ef1c2d9 410 idfla=int(1.+rangen()*real(nfla))
411 idflb=-idfla
412 else !..................g -> gg
413 Z=(1.-ZC)*(ZC/(1.-ZC))**rangen()
60ed90b3 414 IF(rangen().GT.0.5) Z=1.-Z
415 IF((1.-Z*(1.-Z))**2.lt.rangen()) GOTO 390
9ef1c2d9 416 idfla=9
417 idflb=9
418 endif
419 ELSE
420 Z=1.-(1.-ZC)*(ZC/(1.-ZC))**rangen() !!........q -> qg
60ed90b3 421 IF(1.+Z**2.LT.2.*rangen()) GOTO 390
9ef1c2d9 422 idfla=idfl
423 idflb=9
424 endif
425
60ed90b3 426 if(alfm/log(q2*z*(1.-z)/qcdlam).lt.rangen()) goto 390
427
9ef1c2d9 428 if(ij(1).ne.j1.or.ij(2).eq.0)then
429 if(E.lt.sqrt(q2))goto 390
430 pzz=sqrt((E-sqrt(q2))*(E+sqrt(q2)))
431 pt2=(E**2*(z*(1.-z)*q2-z*amm2(idflb)-(1.-z)*amm2(idfla))
432 $ -.25*(amm2(idflb)-amm2(idfla)-q2)**2+q2*amm2(idfla))/pzz**2
433 if (pt2.le.0.) then
434 if(ish.ge.9)then
435 write(ifch,*) 'z not good for pt2:',z,pt2
436 endif
437 goto 390
438 endif
439 endif
440
441 zeta = sqrt(q2)/E/sqrt(z*(1.-z))
442 if(zetamx.gt.0)then
443 if (zeta.gt.zetamx)then
444 if(ish.ge.9)then
445 write(ifch,*) zeta,' > ',zetamx,'zeta-Ablehnung'
446 endif
447 goto 390
448 endif
449 endif
450 999 continue
451c......................................................................
60ed90b3 452
9ef1c2d9 453 pprt(5,ij(ii))=sqrt(q2)
454 id(1,ii)=idfla
455 id(2,ii)=idflb
456 pz(ij(ii))=z
457 if(first)then
458 ii=2
459 first=.false.
460 if(ii2.eq.2)goto 10
461 endif
462
463 if(ij(2).eq.0)then
464 z1=1.
465 pt=0.
466 alpha=0.
467 pprt(3,ij(1))=sqrt(E**2-pprt(5,ij(1))**2)
468 elseif(ij(1).eq.j1.and.ij(2).eq.j2)then
469 E=pprt(4,ij(1))+pprt(4,ij(2))
470 ee(1)=E*.5+(pprt(5,ij(1))**2-pprt(5,ij(2))**2)/2./E
471 ee(2)=E-ee(1)
472 ii=1
473 do while (ii.le.2)
474 if(ee(ii)-pprt(5,ij(ii)).lt.0.) then
475 if(ish.ge.5)then
476 write(ifch,*) 'goto 11'
477 endif
478 ii=1
479 if ( pprt(5,ij(1))**2-amm2(idprt(ij(1))).lt.
480 $ pprt(5,ij(2))**2-amm2(idprt(ij(2))) ) ii=2
481 goto 10
482 endif
483c zc=.5*(1.-sqrt(1.-pprt(5,ij(ii))**2/ee(ii)**2))
484c if(pz(ij(ii)).lt.zc.or.pz(ij(ii)).gt.1.-zc)then
485c if(ish.ge.7)then
486c write(ifch,*) 'first branching rejected'
487c endif
488c goto 10
489c endif
490 z=pz(ij(ii))
491 q2=pprt(5,ij(ii))**2
492 pzz=sqrt((ee(ii)-pprt(5,ij(ii)))*(ee(ii)+pprt(5,ij(ii))))
493 if(pzz.gt.0.)then
494 pt2=(ee(ii)**2*(z*(1.-z)*q2-z*amm2(id(2,ii))
495 $ -(1.-z)*amm2(id(1,ii)))
496 $ -.25*(amm2(id(2,ii))-amm2(id(1,ii))-q2)**2
497 $ +q2*amm2(id(1,ii)))/pzz**2
498 else
499 pt2=0.
500 endif
501 if(id(1,ii).ne.0.and.pt2.le.0.)then
502 if(ish.ge.7)then
503 write(ifch,*) 'first branching rejected for pt2',ii
504 $ ,z1,q2,ee(ii),pprt(5,ij(ii)),id(1,ii),id(2,ii)
505 endif
506 goto 10
507 endif
508 ii=ii+1
509 enddo
510 z1=ee(1)/E
511 z2=ee(2)/E
512 if(ish.ge.7)then
513 write(ifch,*) 'z of first branching',z1
514 endif
515 pprt(3,ij(1))= sqrt(max(0.,ee(1)**2-pprt(5,ij(1))**2))
516 pprt(3,ij(2))=-sqrt(max(0.,ee(2)**2-pprt(5,ij(2))**2))
517 pt=0.
518 alpha=0.
519 else
520 E=pprt(4,io)
521 z1=pz(io)
522 z2=1.-z1
523 am0=pprt(5,io)
524 am1=pprt(5,ij(1))
525 am2=pprt(5,ij(2))
526 aM=am2**2-am0**2-am1**2
527 pzz=sqrt((E-am0)*(E+am0))
528 pprt(3,ij(1))=.5*(aM+2.*z1*E**2)/pzz
529 pprt(3,ij(2))=pzz-pprt(3,ij(1))
530 pt2=(E**2*(z1*z2*am0**2-z1*am2**2-z2*am1**2)
531 $ -.25*aM**2+am0**2*am1**2)/pzz**2
532 if(ish.ge.8)then
533 write(ifch,*) 'pt2,pzz=',pt2,pzz,z1**2*E*pprt(5,io),z1,E
534 $ ,pprt(5,io)
535 endif
60ed90b3 536 111 if(pt2.lt.0.) then
9ef1c2d9 537 ii=1
538 if ( pprt(5,ij(1))**2-amm2(idprt(ij(1))).lt.
539 $ pprt(5,ij(2))**2-amm2(idprt(ij(2))) ) ii=2
60ed90b3 540 goto 10
9ef1c2d9 541 endif
542 pt=sqrt(pt2)
543 alpha=2.*pi*rangen()
544 endif
545 pprt(4,ij(1))=z1*E
546 pprt(1,ij(1))=cos(alpha)*pt
547 pprt(2,ij(1))=sin(alpha)*pt
548 if(ii2.eq.2)then
549 pprt(4,ij(2))=z2*E
550 pprt(1,ij(2))=-cos(alpha)*pt
551 pprt(2,ij(2))=-sin(alpha)*pt
552 endif
6801ca52 553 if(ij(1).ne.j1)then
554 if(pprt(1,io).ne.0..or.pprt(2,io).ne.0.)then
9ef1c2d9 555 do ii=1,2
556 call utrota(-1,pprt(1,io),pprt(2,io),pprt(3,io)
557 & ,pprt(1,ij(ii)),pprt(2,ij(ii)),pprt(3,ij(ii)))
558 enddo
559 endif
6801ca52 560 endif
561 if(ij(1).ne.j1)then
562 if(pprt(3,io).lt.0.)then
9ef1c2d9 563 do k=1,3
564 do ii=1,2
565 pprt(k,ij(ii)) = -pprt(k,ij(ii))
566 enddo
567 enddo
568 endif
6801ca52 569 endif
9ef1c2d9 570 do ii=1,ii2
571 if(id(1,ii).ne.0)then
572 idprt(nprtj+1)=id(1,ii)
573 idprt(nprtj+2)=id(2,ii)
574 pprt(4,nprtj+1)=pz(ij(ii))*pprt(4,ij(ii))
575 pprt(5,nprtj+1)=pz(ij(ii))*pprt(5,ij(ii))
576 pprt(4,nprtj+2)=(1.-pz(ij(ii)))*pprt(4,ij(ii))
577 pprt(5,nprtj+2)=(1.-pz(ij(ii)))*pprt(5,ij(ii))
578 iorprt(nprtj+1)=ij(ii)
579 iorprt(nprtj+2)=ij(ii)
580 idaprt(1,ij(ii))=nprtj+1
581 idaprt(2,ij(ii))=nprtj+2
582 idaprt(1,nprtj+1)=0
583 idaprt(2,nprtj+1)=0
584 idaprt(1,nprtj+2)=0
585 idaprt(2,nprtj+2)=0
586 nprtj=nprtj+2
587 else
588 idaprt(1,ij(ii))=0
589 idaprt(2,ij(ii))=0
590 endif
591 enddo
592 if(ish.ge.5)then
593 if(ij(1).ne.j1)then
594 write(ifch,98) io,(pprt(j,io),j=1,5),pz(io)
595 & ,idprt(io)
596 write(ifch,*) '->'
597 endif
598 write(ifch,99) ij(1),(pprt(j,ij(1)),j=1,5),pz(ij(1))
599 & ,idprt(ij(1)),'->',id(1,1),id(2,1)
600 if(ij(2).ne.0)write(ifch,99) ij(2),
601 & (pprt(j,ij(2)),j=1,5),pz(ij(2))
602 & ,idprt(ij(2)),'->',id(1,2),id(2,2)
603 write(ifch,*)
604 endif
605
606 98 format(i4,6g10.3,1i4)
607 99 format(i4,6g10.3,1i4,a,2i4)
608
609
610 if(ij(1).le.nn)ij(1)=nn-1
611 ij(1)=ij(1)+2
612 ij(2)=ij(1)+1
613 ii=1
614 ii2=2
615 first=.true.
616 if(ij(1).le.nprtj)goto10
617
618 if(ish.ge.5)then
619 write(ifch,*)
620 do i=1,nprtj
60ed90b3 621 write(ifch,'(i4,1x,5g10.4,2i4,a,2i4)')
9ef1c2d9 622 & i,(pprt(j,i),j=1,5),idprt(i)
623 & ,iorprt(i),' -->',idaprt(1,i),idaprt(2,i)
624 enddo
625 write(ifch,*)
626 endif
627 9999 call utprix('timsho',ish,ishini,5)
628 return
629 end
630