19 SUBROUTINE pyptis(MODE,PT2NOW,PT2CUT,PT2,IFAIL)
22 IMPLICIT DOUBLE PRECISION(a-
h, o-
z)
26 parameter(maxnur=1000)
28 common/pypart/
npart,npartd,ipart(maxnur),ptpart(maxnur)
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)
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)
43 & /
pyint2/,/pyintm/,/pyismx/,/pyctag/,/pyisjn/
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
66 aem2pi=paru(101)/paru(2)
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)
81 IF (
mstp(61).GE.1)
THEN
90 IF(kfbeam(js).NE.22.AND.(kflcb.EQ.4.OR.kflcb.EQ.5))
THEN
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.')
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 '//
124 ELSEIF(mode.EQ.0)
THEN
130 IF (
mint(44+js).EQ.1)
RETURN
135 IF (nisgen(1,mi).EQ.0.AND.nisgen(2,mi).EQ.0) shtnow(mi)=shat
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
143 IF(mecor.GE.1) CALL
pymemx(mecor,wtff,wtgf,wtfg,wtgg)
146 kflb=
k(imi(js,mi,1),2)
150 ksvcb=
max(-1,imi(js,mi,2))
153 IF(imi(js,mi,2).GT.imisep(mi)) ksvcb=-1
155 xb=xmi(js,mi)/
vint(142+js)
159 IF (kflba.EQ.4.OR.kflba.EQ.5)
THEN
161 IF (kflba.EQ.5) rmq2=rmb2
163 IF (kfbeam(js).NE.22) mqmass=kflba
174 ELSEIF(mqmass.EQ.0)
THEN
181 xq0=
max(1d-10,xpsvc(kflb,ksvcb))
189 xfb(kflb)=xpsvc(kflb,ksvcb)
194 IF (imi(js,misea,2).NE.imi(js,mi,1)) goto 120
206 130
IF (
pt2.GT.tmin*rmb2)
THEN
208 pt2mne=
max(tmin*rmb2,pt2cut)
211 ELSEIF(
pt2.GT.tmin*rmc2)
THEN
213 pt2mne=
max(tmin*rmc2,pt2cut)
225 IF (mqmass.EQ.0)
THEN
227 zmax=1d0-0.5d0*(pt2mne/shtnow(mi))*(sqrt(1d0+4d0*shtnow(mi)
231 fmq=sqrt(rmq2/shtnow(mi))
237 IF(
pt2.LT.pt2cut.OR.zmax.LE.zmin)
RETURN
247 IF (
mstp(96).NE.0)
THEN
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)
263 wtap(kflb)=(8d0/3d0)*
log((1d0-zmin)/(1d0-zmax))
265 rmin=(1+sqrt(zmin))/(1-sqrt(zmin))
266 rmax=(1+sqrt(zmax))/(1-sqrt(zmax))
267 wtap(kflb)=(8d0/3d0)*
log(rmax/rmin)
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)
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!")
288 IF (
mstp(96).NE.0.AND.mjoind(js,mi).EQ.0)
THEN
292 IF (imi(js,mj,2).NE.imi(js,mi,1)) goto 160
293 IF (mjoind(js,mj).EQ.0)
THEN
296 wtapj(mj)=
z*(1d0-
z)*0.5d0*(
z**2+(1d0-
z)**2)
297 IF (wtapj(mj).GT.1d-6)
THEN
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
316 ELSEIF (imi(js,mi,2).GE.0)
THEN
319 ELSEIF (mqmass.EQ.0)
THEN
321 wtap(21)=wtap(21)*1.25d0
325 ELSEIF(kflb.EQ.21)
THEN
329 wtapq=(16d0/3d0)*(sqrt(1d0/zmin)-sqrt(1d0/zmax))
331 DO 180 kfl=1,min(3,
mstp(58))
335 wtap(21)=6d0*
log(zmax*(1d0-zmin)/(zmin*(1d0-zmax)))
336 IF(mecor.GE.1.AND.nisgen(js,mi).EQ.0)
THEN
338 wtap(21)=wtgg*wtap(21)
341 IF (
mstp(96).NE.0.AND.
mint(31).GE.2.AND.mjoind(js,mi).EQ.0)
344 IF (mj.EQ.mi.OR.mjoind(js,mj).NE.0) goto 190
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)
350 IF (mj.GT.mi.AND.kflc.EQ.21) goto 190
351 z=xmi(js,mj)/(xmi(js,mi)+xmi(js,mj))
353 wtapj(mj)=6d0*(
z**2+(1d0-
z)**2)
355 wtapj(mj)=
z*4d0/3d0*(1d0+
z**2)
357 IF (wtapj(mj).GT.1d-6)
THEN
367 IF (mqmass.NE.0)
THEN
368 rml=(rmq2+
vint(18))/alam2
381 CALL
pyerrm(9,
'(PYPTIS:) failed to evolve shower.')
388 xfbo=
max(1d-10,xfb(kflb))
390 wtpdf(kfl)=xfb(kfl)/xfbo
391 wtsum=wtsum+wtap(kfl)*wtpdf(kfl)
395 wtpdf(21)=xfb(21)/xfbo
396 wtsum=wtsum+wtap(21)*wtpdf(21)
398 wtsum=
max(0.0001d0,wtsum)
402 IF (
mstp(96).NE.0.AND.njn.NE.0)
THEN
404 IF (wtapj(mj).LT.1d-3) goto 220
407 kflc=
k(imi(js,mj,1),2)
409 ksvcc=
max(-1,imi(js,mj,2))
410 IF (imi(js,mj,2).GT.imisep(mj)) ksvcc=-1
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!'
420 wtpdfj(mj)=wtpdfj(mj)/xfj(kflc)
424 IF (kflca.EQ.21.AND.kflba.LE.5)
THEN
427 ELSEIF (kflba.EQ.21.AND.kflca.LE.5)
THEN
435 IF (kflba.EQ.21)
mint(36)=mj
438 IF (iabs(kfla).LE.6) xfj(kfla)=xpsvc(kfla,ksvca)
443 wtpdfj(mj)=xfj(kfla)*wtpdfj(mj)
444 wtjoin=wtjoin+wtapj(mj)*wtpdfj(mj)
456 pt2qed=(pt2old+
vint(18))*
pyr(0)**(1d0/(aem2pi*wtape))-
vint(18)
457 IF(pt2qed.GT.
pt2)
THEN
466 IF (mqmass.NE.0)
THEN
467 pt2cr=(rmq2+
vint(18))*(rml**(tpm/(tpl*
pyr(0)**(-tml/wn)-tpm)))
470 IF (pt2cr.LT.tmin*rmq2)
THEN
472 IF (nthres.GT.50)
THEN
473 CALL
pyerrm(9,
'(PYPTIS:) no phase space left for '//
474 &
'massive quark creation. Gave up trying.')
483 IF (pt2cr.GT.
pt2)
THEN
500 IF (
mstp(96).NE.0.AND.njn.NE.0)
THEN
501 pt2jn=alam2*((pt2old+
vint(18))/alam2)**(
pyr(0)**(b0/wtjoin))
503 IF (pt2jn.GE.
pt2)
THEN
510 IF(izrg.EQ.3.AND.
pt2.LT.rmb2)
THEN
513 ELSEIF(izrg.EQ.2.AND.
pt2.LT.rmc2)
THEN
522 IF (
pt2.LT.pt2mx.OR.
pt2.LT.pt2cut)
RETURN
525 IF (mqmass.EQ.0.AND.kflc.NE.22.AND.mjoin.EQ.0)
THEN
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
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.'//
548 IF (kflba.LE.6) xfb(kflb)=xpsvc(kflb,ksvcb)
553 wtveto=1d0/wtpdfj(mj)/xfb(kflb)
554 kflc=
k(imi(js,mj,1),2)
556 ksvcc=
max(-1,imi(js,mj,2))
557 IF (ksvcb.GE.1) ksvcc=-1
563 IF (kflca.LE.6.AND.ksvcc.LE.0) xfj(kflc)=xpsvc(kflc,ksvcc)
564 wtveto=wtveto/xfj(kflc)
568 IF (kflca.EQ.21.AND.kflba.LE.5)
THEN
571 ELSEIF (kflba.EQ.21.AND.kflca.LE.5)
THEN
577 IF (kflb.EQ.21)
mint(36)=mj
580 IF (iabs(kfla).LE.6) xfj(kfla)=xpsvc(kfla,ksvca)
584 wtveto=wtveto*xfj(kfla)
586 IF (wtveto.LT.
pyr(0)) goto 200
588 IF (
pt2.GT.pt2mx)
THEN
607 IF (kflaa.LE.5.AND.kflba.LE.5)
THEN
609 z=1d0-(1d0-zmin)*((1d0-zmax)/(1d0-zmin))**
pyr(0)
611 zfac=rmin*(rmax/rmin)**
pyr(0)
612 z=((1-zfac)/(1+zfac))**2
616 IF (kflba.GE.4) wtz=wtz-
z*(1d0-
z)**2*rmq2/
pt2
618 IF (ksvcb.GE.0) wtz=wtz*sqrt(
z)
622 ELSEIF (kflaa.LE.5.AND.kflb.EQ.21)
THEN
624 z=zmax/(1d0+
pyr(0)*(sqrt(zmax/zmin)-1d0))**2
625 wtz=0.5d0*(1d0+(1d0-
z)**2)*sqrt(
z)
628 ELSEIF (kfla.EQ.21.AND.kflba.LE.5)
THEN
630 z=zmin+
pyr(0)*(zmax-zmin)
633 IF (mqmass.NE.0)
THEN
634 wtz=wtz+2d0*
z*(1d0-
z)*rmq2/
pt2
636 ELSEIF (ksvcb.LT.0)
THEN
641 ELSEIF (kfla.EQ.21.AND.kflb.EQ.21)
THEN
643 z=1d0/(1d0+((1d0-zmin)/zmin)*((1d0-zmax)*zmin/
644 & (zmax*(1d0-zmin)))**
pyr(0))
645 wtz=(1d0-
z*(1d0-
z))**2
650 IF (kflba.GE.4) q2b=q2b-rmq2
654 pt2adj=q2b-
z*(shtnow(mi)+q2b)*(q2b+rm2c)/shtnow(mi)
655 IF (pt2adj.LT.1d-6) goto 230
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))
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)
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)
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)
677 ELSEIF(kfla.EQ.21.AND.kflb.EQ.21)
THEN
678 CALL
pymewt(mecor,4,q2b*shat/shtnow(mi),
z,
phi,wtme)
687 IF (kflba.LE.6.AND.ksvcb.LE.0) xfn(kflb)=xpsvc(kflb,ksvcb)
691 IF(xfbn.LT.1d-20)
THEN
692 IF(kfla.EQ.kflb)
THEN
709 IF (kflba.LE.5.AND.kflaa.LE.5)
THEN
712 xfa(kfla)=xpsvc(kfla,ksvcb)
719 IF(xfan.LT.1d-20)
THEN
726 wtpdfa=1d0/wtpdf(kfla)
727 wttot=wtz*xfan/xfbn*wtpdfa
728 ELSEIF(mcrqq.EQ.1)
THEN
730 wttot=wtcrqq*wtz*xfan/xfbn*wtpdfa
732 ELSEIF(mcrqq.EQ.2)
THEN
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
747 IF(mstu(29).EQ.0) mstu(23)=mstu(23)-1
749 &
'(PYPTIS:) Weight '//chwt//
' above unity')
750 IF (
pt2.GT.ptemax) ptemax=
pt2
751 IF (wttot.GT.wtemax) wtemax=wttot
754 &
'(PYPTIS:) Weight '//chwt//
' above unity')
765 IF(
pt2.GT.pt2mx)
THEN
780 ELSEIF (mode.EQ.1)
THEN
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
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
820 IF (imi(js,ji,2).EQ.is) imi(js,ji,2)=im
823 IF (ipart(ifs).GE.
it) ipart(ifs)=ipart(ifs)+2
844 k(im,3)=
mint(83)+js+2
845 IF(mi.GE.2)
k(im,3)=
mint(83)+js
851 IF(
k(
it,2).EQ.22)
THEN
856 ELSEIF(
k(im,2).GT.0.AND.
k(im,2).LE.5.AND.
k(
it,2).EQ.21)
THEN
860 ELSEIF(
k(im,2).GT.0.AND.
k(im,2).LE.5)
THEN
864 ELSEIF(
k(im,2).LT.0.AND.
k(im,2).GE.-5.AND.
k(
it,2).EQ.21)
THEN
868 ELSEIF(
k(im,2).LT.0.AND.
k(im,2).GE.-5)
THEN
872 ELSEIF((
k(
it,2).EQ.21.AND.
pyr(0).GT.0.5d0).OR.
k(
it,2).LT.0)
THEN
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
884 k(id1,5)=
k(id1,5)+mstu(5)*id2
885 k(id2,4)=
k(id2,4)+mstu(5)*id1
887 IF(
k(
it,1).EQ.1)
THEN
900 IF (
k(is,kcs)/mstu(5)**2.LE.1) CALL
pycttr(is,-kcs,im)
901 IF(
mint(51).NE.0)
RETURN
906 IF (
k(
it,kcs)/mstu(5)**2.LE.1) CALL
pycttr(
it,kcs,im)
907 IF(
mint(51).NE.0)
RETURN
911 betaz=
side*(1d0-(1d0+q2bmx/shat)**2)/
912 & (1d0+(1d0+q2bmx/shat)**2)
914 CALL
pyrobo(ir,ir,0d0,0d0,0d0,0d0,betaz)
920 IF (mi.EQ.1) imin=
mint(83)+5
924 CALL
pyrobo(imin,imax,0d0,-phimx,0d0,0d0,0d0)
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)
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)
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
944 IF (
p(is,5).LT.0d0)
THEN
945 p(is,5)=-sqrt(abs(
p(is,5)))
947 p(is,5)=sqrt(
p(is,5))
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')
960 CALL
pyrobo(imin,imax,0d0,0d0,-betax,0d0,-betaz)
963 CALL
pyrobo(imin,imax,-theta,phimx,0d0,0d0,0d0)
971 IF (
k(
it,2).NE.22)
THEN
974 ptpart(
npart)=sqrt(pt2amx)
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
986 IF (ksa.EQ.21.AND.kma.GE.1.AND.kma.LE.5)
THEN
990 CALL
pyptmi(2,pt2now,ptdum1,ptdum2,ifail)
991 IF(
mint(51).NE.0)
RETURN
993 IF(ksa.GE.1.AND.ksa.LE.5.AND.kma.EQ.21)
THEN
998 CALL
pyerrm(9,
'(PYPTIS:) Sorry, companion quark radiated'
999 & //
' away. Cannot handle that yet. Giving up.')
1002 ELSEIF(icmp.LT.0)
THEN
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)
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
1017 nvc(js,ifl)=nvc(js,ifl)-1
1021 ELSEIF(ksa.GE.1.AND.ksa.LE.5.AND.kma.NE.21)
THEN
1025 IF (imi(js,mi,2).LT.0)
THEN
1028 xassoc(js,ifl,icmp)=xmi(jsmx,mimx)