Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyptis.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyptis.f
1 C*********************************************************************
2 
3 C...PYPTIS
4 C...Generates pT-ordered spacelike initial-state parton showers and
5 C...trial joinings.
6 C...MODE=-1: Initialize ISR from scratch, starting from the hardest
7 C... interaction initiators at PT2NOW.
8 C...MODE= 0: Generate a trial branching on interaction MINT(36), side
9 C... MINT(30). Start evolution at PT2NOW, solve Sudakov for PT2.
10 C... Store in /PYISMX/ if PT2 is largest so far. Abort if PT2
11 C... is below PT2CUT.
12 C... (Also generate test joinings if MSTP(96)=1.)
13 C...MODE= 1: Accept stored shower branching. Update event record etc.
14 C...PT2NOW : Starting (max) PT2 scale for evolution.
15 C...PT2CUT : Lower limit for evolution.
16 C...PT2 : Result of evolution. Generated PT2 for trial emission.
17 C...IFAIL : Status return code. IFAIL=0 when all is well.
18 
19  SUBROUTINE pyptis(MODE,PT2NOW,PT2CUT,PT2,IFAIL)
20 
21 C...Double precision and integer declarations.
22  IMPLICIT DOUBLE PRECISION(a-h, o-z)
23  IMPLICIT INTEGER(i-n)
24  INTEGER pyk,pychge,pycomp
25 C...Parameter statement for maximum size of showers.
26  parameter(maxnur=1000)
27 C...Commonblocks.
28  common/pypart/npart,npartd,ipart(maxnur),ptpart(maxnur)
29  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
30  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
31  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
32  common/pypars/mstp(200),parp(200),msti(200),pari(200)
33  common/pyint1/mint(400),vint(400)
34  common/pyint2/iset(500),kfpr(500,2),coef(500,20),icol(40,4,2)
35  common/pyintm/kfival(2,3),nmi(2),imi(2,800,2),nvc(2,-6:6),
36  & xassoc(2,-6:6,240),xpsvc(-6:6,-1:240),pvctot(2,-1:1),
37  & xmi(2,240),pt2mi(240),imisep(0:240)
38  common/pyismx/mimx,jsmx,kflamx,kflcmx,kfbeam(2),nisgen(2,240),
39  & pt2mx,pt2amx,zmx,rm2cmx,q2bmx,phimx
40  common/pyctag/nct,mct(4000,2)
41  common/pyisjn/mjn1mx,mjn2mx,mjoind(2,240)
42  SAVE /pypart/,/pyjets/,/pydat1/,/pydat2/,/pypars/,/pyint1/,
43  & /pyint2/,/pyintm/,/pyismx/,/pyctag/,/pyisjn/
44 C...Local variables
45  dimension zsav(2,240),pt2sav(2,240),
46  & xfb(-25:25),xfa(-25:25),xfn(-25:25),xfj(-25:25),
47  & wtap(-25:25),wtpdf(-25:25),shtnow(240),
48  & wtapj(240),wtpdfj(240),x1(240),y(240)
49  SAVE zsav,pt2sav,xfb,xfa,xfn,wtap,wtpdf,xmxc,shtnow,
50  & rmb2,rmc2,alam3,alam4,alam5,tmin,ptemax,wtemax,aem2pi
51 C...For check on excessive weights.
52  CHARACTER chwt*12
53 
54  DATA ptemax /0d0/
55  DATA wtemax /0d0/
56 
57  ifail=-1
58 
59 C----------------------------------------------------------------------
60 C...MODE=-1: Initialize initial state showers from scratch, i.e.
61 C...starting from the hardest interaction initiators.
62  IF (mode.EQ.-1) THEN
63 C...Set hard scattering SHAT.
64  shtnow(1)=vint(44)
65 C...Mass thresholds and Lambda for QCD evolution.
66  aem2pi=paru(101)/paru(2)
67  rmb=pmas(5,1)
68  rmc=pmas(4,1)
69  alam4=parp(61)
70  IF(mstu(112).LT.4) alam4=parp(61)*(parp(61)/rmc)**(2d0/25d0)
71  IF(mstu(112).GT.4) alam4=parp(61)*(rmb/parp(61))**(2d0/25d0)
72  alam5=alam4*(alam4/rmb)**(2d0/23d0)
73  alam3=alam4*(rmc/alam4)**(2d0/27d0)
74  rmb2=rmb**2
75  rmc2=rmc**2
76 C...Massive quark forced creation threshold (in M**2).
77  tmin=1.01d0
78 C...Set upper limit for X (ensures some X left for beam remnant).
79  xmxc=1d0-2d0*parp(111)/vint(1)
80 
81  IF (mstp(61).GE.1) THEN
82 C...Initial values: flavours, momenta, virtualities.
83  DO 100 js=1,2
84  nisgen(js,1)=0
85 
86 C...Special kinematics check for c/b quarks (that g -> c cbar or
87 C...b bbar kinematically possible).
88  kflb=k(imi(js,1,1),2)
89  kflcb=iabs(kflb)
90  IF(kfbeam(js).NE.22.AND.(kflcb.EQ.4.OR.kflcb.EQ.5)) THEN
91 C...Check PT2MAX > mQ^2
92  IF (vint(56).LT.1.05d0*pmas(pycomp(kflcb),1)**2) THEN
93  CALL pyerrm(9,'(PYPTIS:) PT2MAX < 1.05 * MQ**2. '//
94  & 'No Q creation possible.')
95  mint(51)=1
96  RETURN
97  ELSE
98 C...Check for physical z values (m == MQ / sqrt(s))
99 C...For creation diagram, x < z < (1-m)/(1+m(1-m))
100  fmq=pmas(kflcb,1)/sqrt(shtnow(1))
101  zmxcr=(1d0-fmq)/(1d0+fmq*(1d0-fmq))
102  IF (xmi(js,1).GT.0.9d0*zmxcr) THEN
103  CALL pyerrm(9,'(PYPTIS:) No physical z value for '//
104  & 'Q creation.')
105  mint(51)=1
106  RETURN
107  ENDIF
108  ENDIF
109  ENDIF
110  100 CONTINUE
111  ENDIF
112 
113  mint(354)=0
114 C...Zero joining array
115  DO 110 mj=1,240
116  mjoind(1,mj)=0
117  mjoind(2,mj)=0
118  110 CONTINUE
119 
120 C----------------------------------------------------------------------
121 C...MODE= 0: Generate a trial branching on interaction MINT(36) side
122 C...MINT(30). Store if emission PT2 scale is largest so far.
123 C...Also generate test joinings if MSTP(96)=1.
124  ELSEIF(mode.EQ.0) THEN
125  ifail=-1
126  mecor=0
127  isub=mint(1)
128  js=mint(30)
129 C...No shower for structureless beam
130  IF (mint(44+js).EQ.1) RETURN
131  mi=mint(36)
132  shat=vint(44)
133 C...Absolute shower max scale = VINT(56)
134  pt2=min(pt2now,vint(56))
135  IF (nisgen(1,mi).EQ.0.AND.nisgen(2,mi).EQ.0) shtnow(mi)=shat
136 C...Define for which processes ME corrections have been implemented.
137  IF(mstp(68).EQ.1.OR.mstp(68).EQ.3) THEN
138  IF(isub.EQ.1.OR.isub.EQ.2.OR.isub.EQ.141.OR.isub.eq
139  & .142.OR.isub.EQ.144) mecor=1
140  IF(isub.EQ.102.OR.isub.EQ.152.OR.isub.EQ.157) mecor=2
141  IF(isub.EQ.3.OR.isub.EQ.151.OR.isub.EQ.156) mecor=3
142 C...Calculate preweighting factor for ME-corrected processes.
143  IF(mecor.GE.1) CALL pymemx(mecor,wtff,wtgf,wtfg,wtgg)
144  ENDIF
145 C...Basic info on daughter for which to find mother.
146  kflb=k(imi(js,mi,1),2)
147  kflba=iabs(kflb)
148 C...KSVCB: -1 for sea or first companion, 0 for valence or gluon, >1 for
149 C...second companion.
150  ksvcb=max(-1,imi(js,mi,2))
151 C...Treat "first" companion of a pair like an ordinary sea quark
152 C...(except that creation diagram is not allowed)
153  IF(imi(js,mi,2).GT.imisep(mi)) ksvcb=-1
154 C...X (rescaled to [0,1])
155  xb=xmi(js,mi)/vint(142+js)
156 C...Massive quarks (use physical masses.)
157  rmq2=0d0
158  mqmass=0
159  IF (kflba.EQ.4.OR.kflba.EQ.5) THEN
160  rmq2=rmc2
161  IF (kflba.EQ.5) rmq2=rmb2
162 C...Special threshold treatment for non-photon beams
163  IF (kfbeam(js).NE.22) mqmass=kflba
164  ENDIF
165 
166 C...Flags for parton distribution calls.
167  mint(105)=mint(102+js)
168  mint(109)=mint(106+js)
169  vint(120)=vint(2+js)
170 
171 C...Calculate initial parton distribution weights.
172  IF(xb.GE.xmxc) THEN
173  RETURN
174  ELSEIF(mqmass.EQ.0) THEN
175  CALL pypdfu(kfbeam(js),xb,pt2,xfb)
176  ELSE
177 C...Initialize massive quark PT2 dependent pdf underestimate.
178  pt20=pt2
179  CALL pypdfu(kfbeam(js),xb,pt20,xfb)
180 C.!.Tentative treatment of massive valence quarks.
181  xq0=max(1d-10,xpsvc(kflb,ksvcb))
182  xg0=xfb(21)
183  tpm0=log(pt20/rmq2)
184  wpdf0=tpm0*xg0/xq0
185  ENDIF
186  IF (kflba.LE.6) THEN
187 C...For quarks, only include respective sea, val, or cmp part.
188  IF (ksvcb.LE.0) THEN
189  xfb(kflb)=xpsvc(kflb,ksvcb)
190  ELSE
191 C...Find companion's companion
192  misea=0
193  120 misea=misea+1
194  IF (imi(js,misea,2).NE.imi(js,mi,1)) goto 120
195  xs=xmi(js,misea)
196  xrem=vint(142+js)
197  ys=xs/(xrem+xs)
198 C...Momentum fraction of the companion quark.
199 C...Rescale from XB = x/XREM to YB = x/(1-Sum_rest) -> factor (1-YS).
200  yb=xb*(1d0-ys)
201  xfb(kflb)=pyfcmp(yb/vint(140),ys/vint(140),mstp(87))
202  ENDIF
203  ENDIF
204 
205 C...Determine overestimated z range: switch at c and b masses.
206  130 IF (pt2.GT.tmin*rmb2) THEN
207  izrg=3
208  pt2mne=max(tmin*rmb2,pt2cut)
209  b0=23d0/6d0
210  alam2=alam5**2
211  ELSEIF(pt2.GT.tmin*rmc2) THEN
212  izrg=2
213  pt2mne=max(tmin*rmc2,pt2cut)
214  b0=25d0/6d0
215  alam2=alam4**2
216  ELSE
217  izrg=1
218  pt2mne=pt2cut
219  b0=27d0/6d0
220  alam2=alam3**2
221  ENDIF
222 C...Divide Lambda by PARP(64) (equivalent to mult pT2 by PARP(64))
223  alam2=alam2/parp(64)
224 C...Overestimated ZMAX:
225  IF (mqmass.EQ.0) THEN
226 C...Massless
227  zmax=1d0-0.5d0*(pt2mne/shtnow(mi))*(sqrt(1d0+4d0*shtnow(mi)
228  & /pt2mne)-1d0)
229  ELSE
230 C...Massive (limit for bremsstrahlung diagram > creation)
231  fmq=sqrt(rmq2/shtnow(mi))
232  zmax=1d0/(1d0+fmq)
233  ENDIF
234  zmin=xb/xmxc
235 
236 C...If kinematically impossible then do not evolve.
237  IF(pt2.LT.pt2cut.OR.zmax.LE.zmin) RETURN
238 
239 C...Reset Altarelli-Parisi and PDF weights.
240  DO 140 kfl=-5,5
241  wtap(kfl)=0d0
242  wtpdf(kfl)=0d0
243  140 CONTINUE
244  wtap(21)=0d0
245  wtpdf(21)=0d0
246 C...Zero joining weights and compute X(partner) and X(mother) values.
247  IF (mstp(96).NE.0) THEN
248  njn=0
249  DO 150 mj=1,mint(31)
250  wtapj(mj)=0d0
251  wtpdfj(mj)=0d0
252  x1(mj)=xmi(js,mj)/(vint(142+js)+xmi(js,mj))
253  y(mj)=(xmi(js,mi)+xmi(js,mj))/(vint(142+js)+xmi(js,mj)
254  & +xmi(js,mi))
255  150 CONTINUE
256  ENDIF
257 
258 C...Approximate Altarelli-Parisi weights (integrated AP dz).
259 C...q -> q, g -> q or q -> q + gamma (already set which).
260  IF(kflba.LE.5) THEN
261 C...Val and cmp quarks get an extra sqrt(z) to smooth their bumps.
262  IF (ksvcb.LT.0) THEN
263  wtap(kflb)=(8d0/3d0)*log((1d0-zmin)/(1d0-zmax))
264  ELSE
265  rmin=(1+sqrt(zmin))/(1-sqrt(zmin))
266  rmax=(1+sqrt(zmax))/(1-sqrt(zmax))
267  wtap(kflb)=(8d0/3d0)*log(rmax/rmin)
268  ENDIF
269  wtap(21)=0.5d0*(zmax-zmin)
270  wtape=(2d0/9d0)*log((1d0-zmin)/(1d0-zmax))
271  IF(mod(kflba,2).EQ.0) wtape=4d0*wtape
272  IF(mecor.GE.1.AND.nisgen(js,mi).EQ.0) THEN
273  wtap(kflb)=wtff*wtap(kflb)
274  wtap(21)=wtgf*wtap(21)
275  wtape=wtff*wtape
276  ENDIF
277  IF (ksvcb.GE.1) THEN
278 C...Kill normal creation but add joining diagrams for cmp quark.
279  wtap(21)=0d0
280  IF (kflba.EQ.4.OR.kflba.EQ.5) THEN
281  CALL pyerrm(9,'(PYPTIS:) Sorry, I got a heavy companion'//
282  & " quark here. Not handled yet, giving up!")
283  pt2=0d0
284  mint(51)=1
285  RETURN
286  ENDIF
287 C...Check for possible joinings
288  IF (mstp(96).NE.0.AND.mjoind(js,mi).EQ.0) THEN
289 C...Find companion's companion.
290  mj=0
291  160 mj=mj+1
292  IF (imi(js,mj,2).NE.imi(js,mi,1)) goto 160
293  IF (mjoind(js,mj).EQ.0) THEN
294  y(mi)=yb+ys
295  z=yb/y(mi)
296  wtapj(mj)=z*(1d0-z)*0.5d0*(z**2+(1d0-z)**2)
297  IF (wtapj(mj).GT.1d-6) THEN
298  njn=1
299  ELSE
300  wtapj(mj)=0d0
301  ENDIF
302  ENDIF
303 C...Add trial gluon joinings.
304  DO 170 mj=1,mint(31)
305  kflc=k(imi(js,mj,1),2)
306  IF (kflc.NE.21.OR.mjoind(js,mj).NE.0) goto 170
307  z=xmi(js,mj)/(xmi(js,mi)+xmi(js,mj))
308  wtapj(mj)=6d0*(z**2+(1d0-z)**2)
309  IF (wtapj(mj).GT.1d-6) THEN
310  njn=njn+1
311  ELSE
312  wtapj(mj)=0d0
313  ENDIF
314  170 CONTINUE
315  ENDIF
316  ELSEIF (imi(js,mi,2).GE.0) THEN
317 C...Kill creation diagram for val quarks and sea quarks with companions.
318  wtap(21)=0d0
319  ELSEIF (mqmass.EQ.0) THEN
320 C...Extra safety factor for massless sea quark creation.
321  wtap(21)=wtap(21)*1.25d0
322  ENDIF
323 
324 C... q -> g, g -> g.
325  ELSEIF(kflb.EQ.21) THEN
326 C...Here we decide later whether a quark picked up is valence or
327 C...sea, so we maintain the extra factor sqrt(z) since we deal
328 C...with the *sum* of sea and valence in this context.
329  wtapq=(16d0/3d0)*(sqrt(1d0/zmin)-sqrt(1d0/zmax))
330 C...new: do not allow backwards evol to pick up heavy flavour.
331  DO 180 kfl=1,min(3,mstp(58))
332  wtap(kfl)=wtapq
333  wtap(-kfl)=wtapq
334  180 CONTINUE
335  wtap(21)=6d0*log(zmax*(1d0-zmin)/(zmin*(1d0-zmax)))
336  IF(mecor.GE.1.AND.nisgen(js,mi).EQ.0) THEN
337  wtapq=wtfg*wtapq
338  wtap(21)=wtgg*wtap(21)
339  ENDIF
340 C...Check for possible joinings (companions handled separately above)
341  IF (mstp(96).NE.0.AND.mint(31).GE.2.AND.mjoind(js,mi).EQ.0)
342  & THEN
343  DO 190 mj=1,mint(31)
344  IF (mj.EQ.mi.OR.mjoind(js,mj).NE.0) goto 190
345  ksvcc=imi(js,mj,2)
346  IF (imi(js,mj,2).GT.imisep(mj)) ksvcc=-1
347  IF (ksvcc.GE.1) goto 190
348  kflc=k(imi(js,mj,1),2)
349 C...Only try g -> g + g once.
350  IF (mj.GT.mi.AND.kflc.EQ.21) goto 190
351  z=xmi(js,mj)/(xmi(js,mi)+xmi(js,mj))
352  IF (kflc.EQ.21) THEN
353  wtapj(mj)=6d0*(z**2+(1d0-z)**2)
354  ELSE
355  wtapj(mj)=z*4d0/3d0*(1d0+z**2)
356  ENDIF
357  IF (wtapj(mj).GT.1d-6) THEN
358  njn=njn+1
359  ELSE
360  wtapj(mj)=0d0
361  ENDIF
362  190 CONTINUE
363  ENDIF
364  ENDIF
365 
366 C...Initialize massive quark evolution
367  IF (mqmass.NE.0) THEN
368  rml=(rmq2+vint(18))/alam2
369  tml=log(rml)
370  tpl=log((pt2+vint(18))/alam2)
371  tpm=log((pt2+vint(18))/rmq2)
372  wn=wtap(21)*wpdf0/b0
373  ENDIF
374 
375 
376 C...Loopback point for iteration
377  ntry=0
378  nthres=0
379  200 ntry=ntry+1
380  IF(ntry.GT.500) THEN
381  CALL pyerrm(9,'(PYPTIS:) failed to evolve shower.')
382  mint(51)=1
383  RETURN
384  ENDIF
385 
386 C... Calculate PDF weights and sum for evolution rate.
387  wtsum=0d0
388  xfbo=max(1d-10,xfb(kflb))
389  DO 210 kfl=-5,5
390  wtpdf(kfl)=xfb(kfl)/xfbo
391  wtsum=wtsum+wtap(kfl)*wtpdf(kfl)
392  210 CONTINUE
393 C...Only add gluon mother diagram for massless KFLB.
394  IF(mqmass.EQ.0) THEN
395  wtpdf(21)=xfb(21)/xfbo
396  wtsum=wtsum+wtap(21)*wtpdf(21)
397  ENDIF
398  wtsum=max(0.0001d0,wtsum)
399  wtsums=wtsum
400 C...Add joining diagrams where applicable.
401  wtjoin=0d0
402  IF (mstp(96).NE.0.AND.njn.NE.0) THEN
403  DO 220 mj=1,mint(31)
404  IF (wtapj(mj).LT.1d-3) goto 220
405  wtpdfj(mj)=1d0/xfbo
406 C...x and x*pdf (+ sea/val) for parton C.
407  kflc=k(imi(js,mj,1),2)
408  kflca=iabs(kflc)
409  ksvcc=max(-1,imi(js,mj,2))
410  IF (imi(js,mj,2).GT.imisep(mj)) ksvcc=-1
411  mint(30)=js
412  mint(36)=mj
413  CALL pypdfu(kfbeam(js),x1(mj),pt2,xfj)
414  mint(36)=mi
415  IF (kflca.LE.6.AND.ksvcc.LE.0) THEN
416  xfj(kflc)=xpsvc(kflc,ksvcc)
417  ELSEIF (ksvcc.GE.1) THEN
418  print*, 'error! parton C is companion!'
419  ENDIF
420  wtpdfj(mj)=wtpdfj(mj)/xfj(kflc)
421 C...x and x*pdf (+ sea/val) for parton A.
422  kfla=21
423  ksvca=0
424  IF (kflca.EQ.21.AND.kflba.LE.5) THEN
425  kfla=kflb
426  ksvca=ksvcb
427  ELSEIF (kflba.EQ.21.AND.kflca.LE.5) THEN
428  kfla=kflc
429  ksvca=ksvcc
430  ENDIF
431  mint(30)=js
432  IF (ksvca.LE.0) THEN
433 C...Consider C the "evolved" parton if B is gluon. Val/sea
434 C...counting will then be done correctly in PYPDFU.
435  IF (kflba.EQ.21) mint(36)=mj
436  CALL pypdfu(kfbeam(js),y(mj),pt2,xfj)
437  mint(36)=mi
438  IF (iabs(kfla).LE.6) xfj(kfla)=xpsvc(kfla,ksvca)
439  ELSE
440 C...If parton A is companion, use Y(MI) and YS in call to PYFCMP.
441  xfj(kfla)=pyfcmp(y(mi)/vint(140),ys/vint(140),mstp(87))
442  ENDIF
443  wtpdfj(mj)=xfj(kfla)*wtpdfj(mj)
444  wtjoin=wtjoin+wtapj(mj)*wtpdfj(mj)
445  220 CONTINUE
446  ENDIF
447 
448 C...Pick normal pT2 (in overestimated z range).
449  230 pt2old=pt2
450  wtsum=wtsums
451  pt2=alam2*((pt2+vint(18))/alam2)**(pyr(0)**(b0/wtsum))-vint(18)
452  kflc=21
453 
454 C...Evolve q -> q gamma separately, pick it if larger pT.
455  IF(kflba.LE.5) THEN
456  pt2qed=(pt2old+vint(18))*pyr(0)**(1d0/(aem2pi*wtape))-vint(18)
457  IF(pt2qed.GT.pt2) THEN
458  pt2=pt2qed
459  kflc=22
460  kfla=kflb
461  ENDIF
462  ENDIF
463 
464 C... Evolve massive quark creation separately.
465  mcrqq=0
466  IF (mqmass.NE.0) THEN
467  pt2cr=(rmq2+vint(18))*(rml**(tpm/(tpl*pyr(0)**(-tml/wn)-tpm)))
468  & -vint(18)
469 C... Ensure mininimum PT2CR and force creation near threshold.
470  IF (pt2cr.LT.tmin*rmq2) THEN
471  nthres=nthres+1
472  IF (nthres.GT.50) THEN
473  CALL pyerrm(9,'(PYPTIS:) no phase space left for '//
474  & 'massive quark creation. Gave up trying.')
475  mint(51)=1
476  RETURN
477  ENDIF
478  pt2=0d0
479  pt2cr=tmin*rmq2
480  mcrqq=2
481  ENDIF
482 C... Select largest PT2 (brems or creation):
483  IF (pt2cr.GT.pt2) THEN
484  mcrqq=max(mcrqq,1)
485  wtsum=0d0
486  pt2=pt2cr
487  kfla=21
488  ELSE
489  mcrqq=0
490  kfla=kflb
491  ENDIF
492 C... Compute logarithms for this PT2
493  tpl=log((pt2+vint(18))/alam2)
494  tpm=log((pt2+vint(18))/(rmq2+vint(18)))
495  wtcrqq=tpm/log(pt2/rmq2)
496  ENDIF
497 
498 C...Evolve joining separately
499  mjoin=0
500  IF (mstp(96).NE.0.AND.njn.NE.0) THEN
501  pt2jn=alam2*((pt2old+vint(18))/alam2)**(pyr(0)**(b0/wtjoin))
502  & -vint(18)
503  IF (pt2jn.GE.pt2) THEN
504  mjoin=1
505  pt2=pt2jn
506  ENDIF
507  ENDIF
508 
509 C...Loopback if crossed c/b mass thresholds.
510  IF(izrg.EQ.3.AND.pt2.LT.rmb2) THEN
511  pt2=rmb2
512  goto 130
513  ELSEIF(izrg.EQ.2.AND.pt2.LT.rmc2) THEN
514  pt2=rmc2
515  goto 130
516  ENDIF
517 
518 C...Speed up shower. Skip if higher-PT acceptable branching
519 C...already found somewhere else.
520 C...Also finish if below lower cutoff.
521 
522  IF (pt2.LT.pt2mx.OR.pt2.LT.pt2cut) RETURN
523 
524 C...Select parton A flavour (massive Q handled above.)
525  IF (mqmass.EQ.0.AND.kflc.NE.22.AND.mjoin.EQ.0) THEN
526  wtran=pyr(0)*wtsum
527  kfla=-6
528  240 kfla=kfla+1
529  wtran=wtran-wtap(kfla)*wtpdf(kfla)
530  IF(kfla.LE.5.AND.wtran.GT.0d0) goto 240
531  IF(kfla.EQ.6) kfla=21
532  ELSEIF (mjoin.EQ.1) THEN
533 C...Tentative joining accept/reject.
534  wtran=pyr(0)*wtjoin
535  mj=0
536  250 mj=mj+1
537  wtran=wtran-wtapj(mj)*wtpdfj(mj)
538  IF(mj.LE.mint(31)-1.AND.wtran.GT.0d0) goto 250
539  IF(mjoind(js,mj).NE.0.OR.mjoind(js,mi).NE.0) THEN
540  CALL pyerrm(9,'(PYPTIS:) Attempted double joining.'//
541  & ' Rejected.')
542  goto 230
543  ENDIF
544 C...x*pdf (+ sea/val) at new pT2 for parton B.
545  IF (ksvcb.LE.0) THEN
546  mint(30)=js
547  CALL pypdfu(kfbeam(js),xb,pt2,xfb)
548  IF (kflba.LE.6) xfb(kflb)=xpsvc(kflb,ksvcb)
549  ELSE
550 C...Companion distributions do not evolve.
551  xfb(kflb)=xfbo
552  ENDIF
553  wtveto=1d0/wtpdfj(mj)/xfb(kflb)
554  kflc=k(imi(js,mj,1),2)
555  kflca=iabs(kflc)
556  ksvcc=max(-1,imi(js,mj,2))
557  IF (ksvcb.GE.1) ksvcc=-1
558 C...x*pdf (+ sea/val) at new pT2 for parton C.
559  mint(30)=js
560  mint(36)=mj
561  CALL pypdfu(kfbeam(js),x1(mj),pt2,xfj)
562  mint(36)=mi
563  IF (kflca.LE.6.AND.ksvcc.LE.0) xfj(kflc)=xpsvc(kflc,ksvcc)
564  wtveto=wtveto/xfj(kflc)
565 C...x and x*pdf (+ sea/val) at new pT2 for parton A.
566  kfla=21
567  ksvca=0
568  IF (kflca.EQ.21.AND.kflba.LE.5) THEN
569  kfla=kflb
570  ksvca=ksvcb
571  ELSEIF (kflba.EQ.21.AND.kflca.LE.5) THEN
572  kfla=kflc
573  ksvca=ksvcc
574  ENDIF
575  IF (ksvca.LE.0) THEN
576  mint(30)=js
577  IF (kflb.EQ.21) mint(36)=mj
578  CALL pypdfu(kfbeam(js),y(mj),pt2,xfj)
579  mint(36)=mi
580  IF (iabs(kfla).LE.6) xfj(kfla)=xpsvc(kfla,ksvca)
581  ELSE
582  xfj(kfla)=pyfcmp(y(mj)/vint(140),ys/vint(140),mstp(87))
583  ENDIF
584  wtveto=wtveto*xfj(kfla)
585 C...Monte Carlo veto.
586  IF (wtveto.LT.pyr(0)) goto 200
587 C...If accept, save PT2 of this joining.
588  IF (pt2.GT.pt2mx) THEN
589  pt2mx=pt2
590  jsmx=2+js
591  mjn1mx=mj
592  mjn2mx=mi
593  wtapj(mj)=0d0
594  njn=0
595  ENDIF
596 C...Exit and continue evolution.
597  goto 380
598  ENDIF
599  kflaa=iabs(kfla)
600 
601 C...Choose z value (still in overestimated range) and corrective weight.
602 C...Unphysical z will be rejected below when Q2 has is computed.
603  wtz=0d0
604 
605 C...Note: ME and MQ>0 give corrections to overall weights, not shapes.
606 C...q -> q + g or q -> q + gamma (already set which).
607  IF (kflaa.LE.5.AND.kflba.LE.5) THEN
608  IF (ksvcb.LT.0) THEN
609  z=1d0-(1d0-zmin)*((1d0-zmax)/(1d0-zmin))**pyr(0)
610  ELSE
611  zfac=rmin*(rmax/rmin)**pyr(0)
612  z=((1-zfac)/(1+zfac))**2
613  ENDIF
614  wtz=0.5d0*(1d0+z**2)
615 C...Massive weight correction.
616  IF (kflba.GE.4) wtz=wtz-z*(1d0-z)**2*rmq2/pt2
617 C...Valence quark weight correction (extra sqrt)
618  IF (ksvcb.GE.0) wtz=wtz*sqrt(z)
619 
620 C...q -> g + q.
621 C...NB: MQ>0 not yet implemented. Forced absent above.
622  ELSEIF (kflaa.LE.5.AND.kflb.EQ.21) THEN
623  kflc=kfla
624  z=zmax/(1d0+pyr(0)*(sqrt(zmax/zmin)-1d0))**2
625  wtz=0.5d0*(1d0+(1d0-z)**2)*sqrt(z)
626 
627 C...g -> q + qbar.
628  ELSEIF (kfla.EQ.21.AND.kflba.LE.5) THEN
629  kflc=-kflb
630  z=zmin+pyr(0)*(zmax-zmin)
631  wtz=z**2+(1d0-z)**2
632 C...Massive correction
633  IF (mqmass.NE.0) THEN
634  wtz=wtz+2d0*z*(1d0-z)*rmq2/pt2
635 C...Extra safety margin for light sea quark creation
636  ELSEIF (ksvcb.LT.0) THEN
637  wtz=wtz/1.25d0
638  ENDIF
639 
640 C...g -> g + g.
641  ELSEIF (kfla.EQ.21.AND.kflb.EQ.21) THEN
642  kflc=21
643  z=1d0/(1d0+((1d0-zmin)/zmin)*((1d0-zmax)*zmin/
644  & (zmax*(1d0-zmin)))**pyr(0))
645  wtz=(1d0-z*(1d0-z))**2
646  ENDIF
647 
648 C...Derive Q2 from pT2.
649  q2b=pt2/(1d0-z)
650  IF (kflba.GE.4) q2b=q2b-rmq2
651 
652 C...Loopback if outside allowed z range for given pT2.
653  rm2c=pymass(kflc)**2
654  pt2adj=q2b-z*(shtnow(mi)+q2b)*(q2b+rm2c)/shtnow(mi)
655  IF (pt2adj.LT.1d-6) goto 230
656 
657 C...Loopback if nonordered in angle/rapidity.
658  IF (mstp(62).GE.3.AND.nisgen(js,mi).GE.1) THEN
659  IF(pt2.GT.((1d0-z)/(z*(1d0-zsav(js,mi))))**2*pt2sav(js,mi))
660  & goto 230
661  ENDIF
662 
663 C...Select phi angle of branching at random.
664  phi=paru(2)*pyr(0)
665 
666 C...Matrix-element corrections for some processes.
667  IF (mecor.GE.1.AND.nisgen(js,mi).EQ.0) THEN
668  IF (kflaa.LE.20.AND.kflba.LE.20) THEN
669  CALL pymewt(mecor,1,q2b*shat/shtnow(mi),z,phi,wtme)
670  wtz=wtz*wtme/wtff
671  ELSEIF((kfla.EQ.21.OR.kfla.EQ.22).AND.kflba.LE.20) THEN
672  CALL pymewt(mecor,2,q2b*shat/shtnow(mi),z,phi,wtme)
673  wtz=wtz*wtme/wtgf
674  ELSEIF(kflaa.LE.20.AND.(kflb.EQ.21.OR.kflb.EQ.22)) THEN
675  CALL pymewt(mecor,3,q2b*shat/shtnow(mi),z,phi,wtme)
676  wtz=wtz*wtme/wtfg
677  ELSEIF(kfla.EQ.21.AND.kflb.EQ.21) THEN
678  CALL pymewt(mecor,4,q2b*shat/shtnow(mi),z,phi,wtme)
679  wtz=wtz*wtme/wtgg
680  ENDIF
681  ENDIF
682 
683 C...Parton distributions at new pT2 but old x.
684  mint(30)=js
685  CALL pypdfu(kfbeam(js),xb,pt2,xfn)
686 C...Treat val and cmp separately
687  IF (kflba.LE.6.AND.ksvcb.LE.0) xfn(kflb)=xpsvc(kflb,ksvcb)
688  IF (ksvcb.GE.1)
689  & xfn(kflb)=pyfcmp(yb/vint(140),ys/vint(140),mstp(87))
690  xfbn=xfn(kflb)
691  IF(xfbn.LT.1d-20) THEN
692  IF(kfla.EQ.kflb) THEN
693  wtap(kflb)=0d0
694  goto 200
695  ELSE
696  xfbn=1d-10
697  xfn(kflb)=xfbn
698  ENDIF
699  ENDIF
700  DO 260 kfl=-5,5
701  xfb(kfl)=xfn(kfl)
702  260 CONTINUE
703  xfb(21)=xfn(21)
704 
705 C...Parton distributions at new pT2 and new x.
706  xa=xb/z
707  mint(30)=js
708  CALL pypdfu(kfbeam(js),xa,pt2,xfa)
709  IF (kflba.LE.5.AND.kflaa.LE.5) THEN
710 C...q -> q + g: only consider respective sea, val, or cmp content.
711  IF (ksvcb.LE.0) THEN
712  xfa(kfla)=xpsvc(kfla,ksvcb)
713  ELSE
714  ya=xa*(1d0-ys)
715  xfa(kflb)=pyfcmp(ya/vint(140),ys/vint(140),mstp(87))
716  ENDIF
717  ENDIF
718  xfan=xfa(kfla)
719  IF(xfan.LT.1d-20) THEN
720  goto 200
721  ENDIF
722 
723 C...If weighting fails continue evolution.
724  wttot=0d0
725  IF (mcrqq.EQ.0) THEN
726  wtpdfa=1d0/wtpdf(kfla)
727  wttot=wtz*xfan/xfbn*wtpdfa
728  ELSEIF(mcrqq.EQ.1) THEN
729  wtpdfa=tpm/wpdf0
730  wttot=wtcrqq*wtz*xfan/xfbn*wtpdfa
731  xbest=tpm/tpm0*xq0
732  ELSEIF(mcrqq.EQ.2) THEN
733 C...Force massive quark creation.
734  wttot=1d0
735  ENDIF
736 
737 C...Loop back if trial emission fails.
738  IF(wttot.GE.0d0.AND.wttot.LT.pyr(0)) goto 200
739  wtacc=((1d0+pt2)/(0.25d0+pt2))**2
740  IF(wttot.LT.0d0) THEN
741  WRITE(chwt,'(1P,E12.4)') wttot
742  CALL pyerrm(19,'(PYPTIS:) Weight '//chwt//' negative')
743  ELSEIF(wttot.GT.wtacc) THEN
744  WRITE(chwt,'(1P,E12.4)') wttot
745  IF (pt2.GT.ptemax.OR.wttot.GE.wtemax) THEN
746 C...Too high weight: write out as error, but do not update error counter.
747  IF(mstu(29).EQ.0) mstu(23)=mstu(23)-1
748  CALL pyerrm(19,
749  & '(PYPTIS:) Weight '//chwt//' above unity')
750  IF (pt2.GT.ptemax) ptemax=pt2
751  IF (wttot.GT.wtemax) wtemax=wttot
752  ELSE
753  CALL pyerrm(9,
754  & '(PYPTIS:) Weight '//chwt//' above unity')
755  ENDIF
756 C...Useful for debugging but commented out for distribution:
757 C print*, 'JS, MI',JS, MI
758 C print*, 'PT:',SQRT(PT2), ' MCRQQ',MCRQQ
759 C print*, 'A -> B C',KFLA, KFLB, KFLC
760 C XFAO=XFBO/WTPDFA
761 C print*, 'WT(Z,XFA,XFB)',WTZ, XFAN/XFAO, XFBO/XFBN
762  ENDIF
763 
764 C...Save acceptable branching.
765  IF(pt2.GT.pt2mx) THEN
766  mimx=mint(36)
767  jsmx=js
768  pt2mx=pt2
769  kflamx=kfla
770  kflcmx=kflc
771  rm2cmx=rm2c
772  q2bmx=q2b
773  zmx=z
774  pt2amx=pt2adj
775  phimx=phi
776  ENDIF
777 
778 C----------------------------------------------------------------------
779 C...MODE= 1: Accept stored shower branching. Update event record etc.
780  ELSEIF (mode.EQ.1) THEN
781  mi=mimx
782  js=jsmx
783  shat=shtnow(mi)
784  side=3d0-2d0*js
785 C...Shift down rest of event record to make room for insertion.
786  it=imisep(mi)+1
787  im=it+1
788  is=imi(js,mi,1)
789  DO 280 i=n,it,-1
790  IF (k(i,3).GE.it) k(i,3)=k(i,3)+2
791  kt1=k(i,4)/mstu(5)**2
792  kt2=k(i,5)/mstu(5)**2
793  id1=mod(k(i,4),mstu(5))
794  id2=mod(k(i,5),mstu(5))
795  im1=mod(k(i,4)/mstu(5),mstu(5))
796  im2=mod(k(i,5)/mstu(5),mstu(5))
797  IF (id1.GE.it) id1=id1+2
798  IF (id2.GE.it) id2=id2+2
799  IF (im1.GE.it) im1=im1+2
800  IF (im2.GE.it) im2=im2+2
801  k(i,4)=kt1*mstu(5)**2+im1*mstu(5)+id1
802  k(i,5)=kt2*mstu(5)**2+im2*mstu(5)+id2
803  DO 270 ix=1,5
804  k(i+2,ix)=k(i,ix)
805  p(i+2,ix)=p(i,ix)
806  v(i+2,ix)=v(i,ix)
807  270 CONTINUE
808  mct(i+2,1)=mct(i,1)
809  mct(i+2,2)=mct(i,2)
810  280 CONTINUE
811  n=n+2
812 C...Also update shifted-down pointers in IMI, IMISEP, and IPART.
813  DO 290 ji=1,mint(31)
814  IF (imi(1,ji,1).GE.it) imi(1,ji,1)=imi(1,ji,1)+2
815  IF (imi(1,ji,2).GE.it) imi(1,ji,2)=imi(1,ji,2)+2
816  IF (imi(2,ji,1).GE.it) imi(2,ji,1)=imi(2,ji,1)+2
817  IF (imi(2,ji,2).GE.it) imi(2,ji,2)=imi(2,ji,2)+2
818  IF (ji.GE.mi) imisep(ji)=imisep(ji)+2
819 C...Also update companion pointers to the present mother.
820  IF (imi(js,ji,2).EQ.is) imi(js,ji,2)=im
821  290 CONTINUE
822  DO 300 ifs=1,npart
823  IF (ipart(ifs).GE.it) ipart(ifs)=ipart(ifs)+2
824  300 CONTINUE
825 C...Zero entries dedicated for new timelike and mother partons.
826  DO 320 i=it,it+1
827  DO 310 j=1,5
828  k(i,j)=0
829  p(i,j)=0d0
830  v(i,j)=0d0
831  310 CONTINUE
832  mct(i,1)=0
833  mct(i,2)=0
834  320 CONTINUE
835 
836 C...Define timelike and new mother partons. History.
837  k(it,1)=3
838  k(it,2)=kflcmx
839  k(im,1)=14
840  k(im,2)=kflamx
841  k(is,3)=im
842  k(it,3)=im
843 C...Set mother origin = side.
844  k(im,3)=mint(83)+js+2
845  IF(mi.GE.2) k(im,3)=mint(83)+js
846 
847 C...Define colour flow of branching.
848  im1=im
849  im2=im
850 C...q -> q + gamma.
851  IF(k(it,2).EQ.22) THEN
852  k(it,1)=1
853  id1=is
854  id2=is
855 C...q -> q + g.
856  ELSEIF(k(im,2).GT.0.AND.k(im,2).LE.5.AND.k(it,2).EQ.21) THEN
857  id1=it
858  id2=is
859 C...q -> g + q.
860  ELSEIF(k(im,2).GT.0.AND.k(im,2).LE.5) THEN
861  id1=is
862  id2=it
863 C...qbar -> qbar + g.
864  ELSEIF(k(im,2).LT.0.AND.k(im,2).GE.-5.AND.k(it,2).EQ.21) THEN
865  id1=is
866  id2=it
867 C...qbar -> g + qbar.
868  ELSEIF(k(im,2).LT.0.AND.k(im,2).GE.-5) THEN
869  id1=it
870  id2=is
871 C...g -> g + g; g -> q + qbar..
872  ELSEIF((k(it,2).EQ.21.AND.pyr(0).GT.0.5d0).OR.k(it,2).LT.0) THEN
873  id1=is
874  id2=it
875  ELSE
876  id1=it
877  id2=is
878  ENDIF
879  IF(im1.EQ.im) k(im1,4)=k(im1,4)+id1
880  IF(im2.EQ.im) k(im2,5)=k(im2,5)+id2
881  k(id1,4)=k(id1,4)+mstu(5)*im1
882  k(id2,5)=k(id2,5)+mstu(5)*im2
883  IF(id1.NE.id2) THEN
884  k(id1,5)=k(id1,5)+mstu(5)*id2
885  k(id2,4)=k(id2,4)+mstu(5)*id1
886  ENDIF
887  IF(k(it,1).EQ.1) THEN
888  k(it,4)=0
889  k(it,5)=0
890  ENDIF
891 C...Update IMI and colour tag arrays.
892  imi(js,mi,1)=im
893  DO 330 mc=1,2
894  mct(it,mc)=0
895  mct(im,mc)=0
896  330 CONTINUE
897  DO 340 jcs=4,5
898  kcs=jcs
899 C...If mother flag not yet set for spacelike parton, trace it.
900  IF (k(is,kcs)/mstu(5)**2.LE.1) CALL pycttr(is,-kcs,im)
901  IF(mint(51).NE.0) RETURN
902  340 CONTINUE
903  DO 350 jcs=4,5
904  kcs=jcs
905 C...If mother flag not yet set for timelike parton, trace it.
906  IF (k(it,kcs)/mstu(5)**2.LE.1) CALL pycttr(it,kcs,im)
907  IF(mint(51).NE.0) RETURN
908  350 CONTINUE
909 
910 C...Boost recoiling parton to compensate for Q2 scale.
911  betaz=side*(1d0-(1d0+q2bmx/shat)**2)/
912  & (1d0+(1d0+q2bmx/shat)**2)
913  ir=imi(3-js,mi,1)
914  CALL pyrobo(ir,ir,0d0,0d0,0d0,0d0,betaz)
915 
916 C...Define system to be rotated and boosted
917 C...(not including the 2 just added partons)
918 C...(but including the docu lines for first interaction)
919  imin=imisep(mi-1)+1
920  IF (mi.EQ.1) imin=mint(83)+5
921  imax=imisep(mi)-2
922 
923 C...Rotate back system in phi to compensate for subsequent rotation.
924  CALL pyrobo(imin,imax,0d0,-phimx,0d0,0d0,0d0)
925 
926 C...Define kinematics of new partons in old frame.
927  imax=imisep(mi)
928  p(im,1)=sqrt(pt2amx)*shat/(zmx*(shat+q2bmx))
929  p(im,3)=0.5d0*sqrt(shat)*((shat-q2bmx)/((shat
930  & +q2bmx)*zmx)+(q2bmx+rm2cmx)/shat)*side
931  p(im,4)=sqrt(p(im,1)**2+p(im,3)**2)
932  p(it,1)=p(im,1)
933  p(it,3)=p(im,3)-0.5d0*(shat+q2bmx)/sqrt(shat)*side
934  p(it,4)=sqrt(p(it,1)**2+p(it,3)**2+rm2cmx)
935  p(it,5)=sqrt(rm2cmx)
936 
937 C...Update internal line, now spacelike
938  p(is,1)=p(im,1)-p(it,1)
939  p(is,2)=p(im,2)-p(it,2)
940  p(is,3)=p(im,3)-p(it,3)
941  p(is,4)=p(im,4)-p(it,4)
942  p(is,5)=p(is,4)**2-p(is,1)**2-p(is,2)**2-p(is,3)**2
943 C...Represent spacelike virtualities as -sqrt(abs(Q2)) .
944  IF (p(is,5).LT.0d0) THEN
945  p(is,5)=-sqrt(abs(p(is,5)))
946  ELSE
947  p(is,5)=sqrt(p(is,5))
948  ENDIF
949 
950 C...Boost entire system and rotate to new frame.
951 C...(including docu lines)
952  betax=(p(im,1)+p(ir,1))/(p(im,4)+p(ir,4))
953  betaz=(p(im,3)+p(ir,3))/(p(im,4)+p(ir,4))
954  IF(betax**2+betaz**2.GE.1d0) THEN
955  CALL pyerrm(1,'(PYPTIS:) boost bigger than unity')
956  mint(51)=1
957  ifail=-1
958  RETURN
959  ENDIF
960  CALL pyrobo(imin,imax,0d0,0d0,-betax,0d0,-betaz)
961  i1=imi(1,mi,1)
962  theta=pyangl(p(i1,3),p(i1,1))
963  CALL pyrobo(imin,imax,-theta,phimx,0d0,0d0,0d0)
964 
965 C...Global statistics.
966  mint(352)=mint(352)+1
967  vint(352)=vint(352)+sqrt(p(it,1)**2+p(it,2)**2)
968  IF (mint(352).EQ.1) vint(357)=sqrt(p(it,1)**2+p(it,2)**2)
969 
970 C...Add parton with relevant pT scale for timelike shower.
971  IF (k(it,2).NE.22) THEN
972  npart=npart+1
973  ipart(npart)=it
974  ptpart(npart)=sqrt(pt2amx)
975  ENDIF
976 
977 C...Update saved variables.
978  shtnow(mimx)=shtnow(mimx)/zmx
979  nisgen(jsmx,mimx)=nisgen(jsmx,mimx)+1
980  xmi(jsmx,mimx)=xmi(jsmx,mimx)/zmx
981  pt2sav(jsmx,mimx)=pt2mx
982  zsav(js,mimx)=zmx
983 
984  ksa=iabs(k(is,2))
985  kma=iabs(k(im,2))
986  IF (ksa.EQ.21.AND.kma.GE.1.AND.kma.LE.5) THEN
987 C...Gluon reconstructs to quark.
988 C...Decide whether newly created quark is valence or sea:
989  mint(30)=js
990  CALL pyptmi(2,pt2now,ptdum1,ptdum2,ifail)
991  IF(mint(51).NE.0) RETURN
992  ENDIF
993  IF(ksa.GE.1.AND.ksa.LE.5.AND.kma.EQ.21) THEN
994 C...Quark reconstructs to gluon.
995 C...Now some guy may have lost his companion. Check.
996  icmp=imi(js,mi,2)
997  IF (icmp.GT.0) THEN
998  CALL pyerrm(9,'(PYPTIS:) Sorry, companion quark radiated'
999  & //' away. Cannot handle that yet. Giving up.')
1000  mint(51)=1
1001  RETURN
1002  ELSEIF(icmp.LT.0) THEN
1003 C...A sea quark with companion still in BR was reconstructed to a gluon.
1004 C...Companion should now be removed from the beam remnant.
1005 C...(Momentum integral is automatically updated in next call to PYPDFU.)
1006  icmp=-icmp
1007  ifl=-k(is,2)
1008  DO 370 jcmp=icmp,nvc(js,ifl)-1
1009  xassoc(js,ifl,jcmp)=xassoc(js,ifl,jcmp+1)
1010  DO 360 ji=1,mint(31)
1011  kmi=-imi(js,ji,2)
1012  jfl=-k(imi(js,ji,1),2)
1013  IF (kmi.EQ.jcmp+1.AND.jfl.EQ.ifl) imi(js,ji,2)=imi(js,ji
1014  & ,2)+1
1015  360 CONTINUE
1016  370 CONTINUE
1017  nvc(js,ifl)=nvc(js,ifl)-1
1018  ENDIF
1019 C...Set gluon IMI(JS,MI,2) = 0.
1020  imi(js,mi,2)=0
1021  ELSEIF(ksa.GE.1.AND.ksa.LE.5.AND.kma.NE.21) THEN
1022 C...Quark reconstructing to quark. If sea with companion still in BR
1023 C...then update associated x value.
1024 C...(Momentum integral is automatically updated in next call to PYPDFU.)
1025  IF (imi(js,mi,2).LT.0) THEN
1026  icmp=-imi(js,mi,2)
1027  ifl=-k(is,2)
1028  xassoc(js,ifl,icmp)=xmi(jsmx,mimx)
1029  ENDIF
1030  ENDIF
1031 
1032  ENDIF
1033 
1034 C...If reached this point, normal exit.
1035  380 ifail=0
1036 
1037  RETURN
1038  END