10 IMPLICIT DOUBLE PRECISION(a-
h, o-
z)
15 common/
pydat1/mstu(200),paru(200),mstj(200),parj(200)
16 common/
pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
17 common/
pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
20 dimension vdcy(4),kflo(4),kfl1(4),pv(10,5),rord(10),ue(3),be(3),
21 &wtcor(10),ptau(4),pcmtau(4),dbetau(3)
23 DATA wtcor/2d0,5d0,15d0,60d0,250d0,1500d0,1.2d4,1.2d5,150d0,16d0/
26 pawt(a,b,
c)=sqrt((a**2-(b+
c)**2)*(a**2-(b-
c)**2))/(2d0*a)
40 ELSEIF(
k(ip,1).NE.4)
THEN
44 vdcy(
j)=
v(ip,
j)+
v(ip,5)*
p(ip,
j)/
p(ip,5)
49 IF(mstj(22).EQ.2)
THEN
50 IF(pmas(kc,4).GT.parj(71)) mout=1
51 ELSEIF(mstj(22).EQ.3)
THEN
52 IF(vdcy(1)**2+vdcy(2)**2+vdcy(3)**2.GT.parj(72)**2) mout=1
53 ELSEIF(mstj(22).EQ.4)
THEN
54 IF(vdcy(1)**2+vdcy(2)**2.GT.parj(73)**2) mout=1
55 IF(abs(vdcy(3)).GT.parj(74)) mout=1
57 IF(mout.EQ.1.AND.
k(ip,1).NE.5)
THEN
63 IF(kfa.EQ.15.AND.mstj(28).GE.1)
THEN
81 ELSEIF(
k(imtau,2).EQ.
k(itau,2))
THEN
83 IF(
k(
k(imtau,4),2).EQ.22)
THEN
85 pcmtau(
j)=pcmtau(
j)+
p(
k(imtau,4),
j)
87 ELSEIF(
k(
k(imtau,5),2).EQ.22)
THEN
89 pcmtau(
j)=pcmtau(
j)+
p(
k(imtau,5),
j)
94 ELSEIF(iabs(
k(imtau,2)).GT.100)
THEN
97 kforig=-isign(24,
k(itau,2))
99 DO 160 ii=
k(imtau,4),
k(imtau,5)
100 IF(
k(ii,2)*isign(1,
k(itau,2)).EQ.-16)
THEN
102 pcmtau(
j)=pcmtau(
j)+
p(ii,
j)
112 DO 170 ii=imtau+1,ip-1
113 IF(
k(ii,2).EQ.kforig.AND.
k(ii,3).EQ.iorig.AND.
114 & abs(
p(ii,5)-
p(iorig,5)).LT.1d-5*
p(iorig,5)) iorig=ii
124 dbetau(
j)=pcmtau(
j)/pcmtau(4)
126 IF(kforig.NE.0) CALL
pyrobo(itau,itau,0d0,0d0,-dbetau(1),
127 & -dbetau(2),-dbetau(3))
129 CALL
pyrobo(itau,itau,0d0,-phitau,0d0,0d0,0d0)
131 CALL
pyrobo(itau,itau,-thetau,0d0,0d0,0d0,0d0)
134 IF(kforig.NE.0.OR.mstj(28).EQ.2)
THEN
135 CALL
pytaud(itau,iorig,kforig,ndecay)
136 DO 200 ii=nsav+1,nsav+ndecay
149 IF(kforig.NE.0.OR.mstj(28).EQ.2)
THEN
150 CALL
pyrobo(nsav+1,
n,thetau,phitau,0d0,0d0,0d0)
151 IF(kforig.NE.0) CALL
pyrobo(nsav+1,
n,0d0,0d0,dbetau(1),
152 & dbetau(2),dbetau(3))
164 IF((kfa.EQ.511.OR.kfa.EQ.531).AND.mstj(26).GE.1)
THEN
166 IF(kfa.EQ.531) xbbmix=parj(77)
167 IF(sin(0.5d0*xbbmix*
v(ip,5)/pmas(kc,4))**2.GT.
pyr(0)) mmix=1
168 IF(mmix.EQ.1) kfs=-kfs
173 IF(mdcy(kc,2).GT.0)
THEN
174 mdmdcy=mdme(mdcy(kc,2),2)
175 IF(mdmdcy.GT.80.AND.mdmdcy.LE.90) kca=mdmdcy
177 IF(mdcy(kca,2).LE.0.OR.mdcy(kca,3).LE.0)
THEN
178 CALL
pyerrm(9,
'(PYDECY:) no decay channel defined')
181 IF(mod(kfa/1000,10).EQ.0.AND.kca.EQ.85) kfs=-kfs
182 IF(kchg(kc,3).EQ.0)
THEN
185 IF(
pyr(0).GT.0.5d0) kfs=-kfs
186 ELSEIF(kfs.GT.0)
THEN
197 DO 230 idl=mdcy(kca,2),mdcy(kca,2)+mdcy(kca,3)-1
198 IF(mdme(idl,1).NE.1.AND.kfsp*mdme(idl,1).NE.2.AND.
199 & kfsn*mdme(idl,1).NE.3) goto 230
200 IF(mdme(idl,2).GT.100) goto 230
205 CALL
pyerrm(2,
'(PYDECY:) all decay channels closed by user')
213 IF(mdme(idl,1).NE.1.AND.kfsp*mdme(idl,1).NE.2.AND.
214 &kfsn*mdme(idl,1).NE.3)
THEN
215 IF(idl.LT.mdcy(kca,2)+mdcy(kca,3)-1) goto 250
216 ELSEIF(mdme(idl,2).GT.100)
THEN
217 IF(idl.LT.mdcy(kca,2)+mdcy(kca,3)-1) goto 250
221 IF(idl.LT.mdcy(kca,2)+mdcy(kca,3)-1.AND.rbr.GT.0d0) goto 250
227 IF(mod(
ntry,200).EQ.0)
THEN
228 WRITE(cidc,
'(I4)') idc
230 IF(kfa.NE.113.AND.kfa.NE.115.AND.kfa.NE.215)
231 & CALL
pyerrm(4,
'(PYDECY:) caught in loop for decay channel'//
235 IF(
ntry.GT.1000)
THEN
236 CALL
pyerrm(14,
'(PYDECY:) caught in infinite loop')
237 IF(mstu(21).GE.1)
RETURN
243 IF(mmat.GE.11.AND.
p(ip,4).GT.20d0*
p(ip,5)) mbst=1
246 IF(mbst.EQ.0) pv(1,
j)=
p(ip,
j)
248 IF(mbst.EQ.1) pv(1,4)=
p(ip,5)
254 IF(kfa.GT.80) mhaddy=1
263 IF(mdme(idc+1,2).EQ.101) jtmax=10
265 IF(jt.LE.5) kp=kfdp(idc,jt)
266 IF(jt.GE.6) kp=kfdp(idc+1,jt-5)
270 IF(kpa.GT.80) mhaddy=1
271 IF(kchg(kcp,3).EQ.0.AND.kpa.NE.81.AND.kpa.NE.82)
THEN
273 ELSEIF(kpa.NE.81.AND.kpa.NE.82)
THEN
275 ELSEIF(kpa.EQ.81.AND.mod(kfa/1000,10).EQ.0)
THEN
276 kfp=-kfs*mod(kfa/10,10)
277 ELSEIF(kpa.EQ.81.AND.mod(kfa/100,10).GE.mod(kfa/10,10))
THEN
278 kfp=kfs*(100*mod(kfa/10,100)+3)
279 ELSEIF(kpa.EQ.81)
THEN
280 kfp=kfs*(1000*mod(kfa/10,10)+100*mod(kfa/100,10)+1)
281 ELSEIF(kp.EQ.82)
THEN
282 CALL
pydcyk(-kfs*int(1d0+(2d0+parj(2))*
pyr(0)),0,kfp,kdump)
283 IF(kfp.EQ.0) goto 260
287 IF(pv(1,5).LT.parj(32)+2d0*
pymass(kfp)) goto 260
288 ELSEIF(kp.EQ.-82)
THEN
291 IF(kpa.EQ.81.OR.kpa.EQ.82) kcp=
pycomp(kfp)
296 IF(mmat.GE.11.AND.mmat.LE.30.AND.kqp.NE.0)
THEN
300 IF(kp.EQ.82.AND.mstu(121).GT.0) jtmo=nq
303 ELSEIF((mmat.EQ.42.OR.mmat.EQ.43.OR.mmat.EQ.48).AND.
np.EQ.3.AND.
304 & mod(nq,2).EQ.1)
THEN
310 IF(
k(
i,2).EQ.0) goto 260
317 IF(mmat.NE.33.AND.kqp.NE.0) nq=nq+1
318 IF(mmat.EQ.33.AND.kqp.NE.0.AND.kqp.NE.2) nq=nq+1
320 IF(mmat.EQ.4.AND.jt.LE.2.AND.kfp.EQ.21)
k(
i,1)=2
321 IF(mmat.EQ.4.AND.jt.EQ.3)
k(
i,1)=1
333 IF(
ps+parj(64).GT.pv(1,5)) goto 240
337 290
IF(mmat.GE.11.AND.mmat.LE.30)
THEN
339 cnde=parj(61)*
log(
max((pv(1,5)-
ps-psq)/parj(62),1.1d0))
340 IF(mmat.EQ.12) cnde=cnde+parj(63)
346 ELSEIF(irndmo.EQ.1)
THEN
351 IF(
ntry.GT.1000)
THEN
352 CALL
pyerrm(14,
'(PYDECY:) caught in infinite loop')
353 IF(mstu(21).GE.1)
RETURN
356 gauss=sqrt(-2d0*cnde*
log(
max(1d-10,
pyr(0))))*
357 & sin(paru(2)*
pyr(0))
358 nd=0.5d0+0.5d0*
np+0.25d0*nq+cnde+gauss
359 IF(nd.LT.
np+nq/2.OR.nd.LT.2.OR.nd.GT.10) goto 300
360 IF(mmat.EQ.13.AND.nd.EQ.2) goto 300
361 IF(mmat.EQ.14.AND.nd.LE.3) goto 300
362 IF(mmat.EQ.15.AND.nd.LE.4) goto 300
368 IF(mstu(121).GT.mstu(125)) goto 300
374 IF(nd.EQ.
np+nq/2) goto 330
375 DO 320
i=
n+
np+1,
n+nd-nq/2
378 IF(jt.EQ.0) jt=1+int((nq-1)*
pyr(0))
380 IF(
k(
i,2).EQ.0) goto 300
381 mstu(125)=mstu(125)-1
383 IF(mstu(121).GT.0) jtmo=jt
389 IF(nq.EQ.4.AND.
pyr(0).LT.parj(66)) jt=4
390 IF(jt.EQ.4.AND.isign(1,kfl1(1)*(10-iabs(kfl1(1))))*
391 & isign(1,kfl1(jt)*(10-iabs(kfl1(jt)))).GT.0) jt=3
394 CALL
pydcyk(kfl1(1),kfl1(jt),kfldmp,
k(
n+nd-nq/2+1,2))
395 IF(
k(
n+nd-nq/2+1,2).EQ.0) goto 300
396 IF(nq.EQ.4) CALL
pydcyk(kfl1(jt2),kfl1(jt3),kfldmp,
k(
n+nd,2))
397 IF(nq.EQ.4.AND.
k(
n+nd,2).EQ.0) goto 300
409 IF(
ps+parj(64).GT.pv(1,5)) goto 300
412 ELSEIF((mmat.EQ.31.OR.mmat.EQ.33.OR.mmat.EQ.44)
415 pqt=(
p(
n+
np,5)+parj(65))/pv(1,5)
418 pv(1,
j)=(1d0-pqt)*pv(1,
j)
420 IF(
ps+parj(64).GT.pv(1,5)) goto 260
426 IF(
np.GE.2.AND.
ps+parj(64).GT.pv(1,5)) goto 260
436 IF(im.LT.0.OR.im.GE.ip) im=0
437 IF(im.NE.0) kfam=iabs(
k(im,2))
439 DO 360 il=
max(ip-2,im+1),min(ip+2,
n)
440 IF(
k(il,3).EQ.im) nm=nm+1
441 IF(
k(il,3).EQ.im.AND.il.NE.ip) isis=il
443 IF(nm.NE.2.OR.kfam.LE.100.OR.mod(kfam,10).NE.1.OR.
444 & mod(kfam/1000,10).NE.0) nm=0
447 IF((kfas.LE.100.OR.mod(kfas,10).NE.1.OR.
448 & mod(kfas/1000,10).NE.0).AND.kfas.NE.22) nm=0
464 wtmax=1d0/wtcor(nd-2)
465 pmax=pv(1,5)-
ps+
p(
n+nd,5)
469 pmin=pmin+
p(
n+il+1,5)
470 wtmax=wtmax*pawt(pmax,pmin,
p(
n+il,5))
476 ELSEIF(mmat.EQ.2)
THEN
477 pmes=4d0*pmas(11,1)**2
478 pmrho2=pmas(131,1)**2
479 pgrho2=pmas(131,2)**2
480 400 pmst=pmes*(
p(ip,5)**2/pmes)**
pyr(0)
481 wt=(1+0.5d0*pmes/pmst)*sqrt(
max(0d0,1d0-pmes/pmst))*
482 & (1d0-pmst/
p(ip,5)**2)**3*(1d0+pgrho2/pmrho2)/
483 & ((1d0-pmst/pmrho2)**2+pgrho2/pmrho2)
484 IF(wt.LT.
pyr(0)) goto 400
485 pv(2,5)=
max(2.00001d0*pmas(11,1),sqrt(pmst))
492 DO 420 il2=il1-1,1,-1
493 IF(rsav.LE.rord(il2)) goto 430
494 rord(il2+1)=rord(il2)
501 pv(il,5)=pv(il+1,5)+
p(
n+il,5)+(rord(il)-rord(il+1))*
503 wt=wt*pawt(pv(il,5),pv(il+1,5),
p(
n+il,5))
505 IF(wt.LT.
pyr(0)*wtmax) goto 410
510 pa=pawt(pv(il,5),pv(il+1,5),
p(
n+il,5))
513 ue(1)=sqrt(1d0-ue(3)**2)*cos(
phi)
514 ue(2)=sqrt(1d0-ue(3)**2)*sin(
phi)
519 p(
n+il,4)=sqrt(pa**2+
p(
n+il,5)**2)
520 pv(il+1,4)=sqrt(pa**2+pv(il+1,5)**2)
529 be(
j)=pv(il,
j)/pv(il,4)
533 bep=be(1)*
p(
i,1)+be(2)*
p(
i,2)+be(3)*
p(
i,3)
535 p(
i,
j)=
p(
i,
j)+ga*(ga*bep/(1d0+ga)+
p(
i,4))*be(
j)
537 p(
i,4)=ga*(
p(
i,4)+bep)
543 IF(
ntry.GT.800) goto 560
547 wt=(
p(
n+1,5)*
p(
n+2,5)*
p(
n+3,5))**2-(
p(
n+1,5)*four(
n+2,
n+3))**2
548 & -(
p(
n+2,5)*four(
n+1,
n+3))**2-(
p(
n+3,5)*four(
n+1,
n+2))**2
549 & +2d0*four(
n+1,
n+2)*four(
n+1,
n+3)*four(
n+2,
n+3)
550 IF(
max(wt*wtcor(9)/
p(ip,5)**6,0.001d0).LT.
pyr(0)) goto 390
553 ELSEIF(mmat.EQ.2)
THEN
556 wt=(pmst-0.5d0*pmes)*(four12**2+four13**2)+
557 & pmes*(four12*four13+four12**2+four13**2)
558 IF(wt.LT.
pyr(0)*0.25d0*pmst*(
p(ip,5)**2-pmst)**2) goto 460
563 ELSEIF(mmat.EQ.3.AND.nm.EQ.2)
THEN
570 IF(kfas.NE.22) hnum=(four10*four12-pms1*four02)**2
571 IF(kfas.EQ.22) hnum=pms1*(2d0*four10*four12*four02-
572 & pms1*four02**2-pms0*four12**2-pms2*four10**2+pms1*pms0*pms2)
573 hnum=
max(1d-6*pms1**2*pms0*pms2,hnum)
574 hden=(four10**2-pms1*pms0)*(four12**2-pms1*pms2)
575 IF(hnum.LT.
pyr(0)*hden) goto 460
578 ELSEIF(mmat.EQ.4)
THEN
579 hx1=2d0*four(ip,
n+1)/
p(ip,5)**2
580 hx2=2d0*four(ip,
n+2)/
p(ip,5)**2
581 hx3=2d0*four(ip,
n+3)/
p(ip,5)**2
582 wt=((1d0-hx1)/(hx2*hx3))**2+((1d0-hx2)/(hx1*hx3))**2+
583 & ((1d0-hx3)/(hx1*hx2))**2
584 IF(wt.LT.2d0*
pyr(0)) goto 390
585 IF(
k(ip+1,2).EQ.22.AND.(1d0-hx1)*
p(ip,5)**2.LT.4d0*parj(32)**2)
589 ELSEIF(mmat.EQ.41)
THEN
590 IF(mbst.EQ.0) hx1=2d0*four(ip,
n+1)/
p(ip,5)**2
591 IF(mbst.EQ.1) hx1=2d0*
p(
n+1,4)/
p(ip,5)
592 hxm=min(0.75d0,2d0*(1d0-
ps/
p(ip,5)))
593 IF(hx1*(3d0-2d0*hx1).LT.
pyr(0)*hxm*(3d0-2d0*hxm)) goto 390
596 ELSEIF((mmat.EQ.42.OR.mmat.EQ.43.OR.mmat.EQ.44.OR.mmat.EQ.48)
598 IF(mbst.EQ.0) wt=four(ip,
n+1)*four(
n+2,
n+3)
599 IF(mbst.EQ.1) wt=
p(ip,5)*
p(
n+1,4)*four(
n+2,
n+3)
600 IF(wt.LT.
pyr(0)*
p(ip,5)*pv(1,5)**3/wtcor(10)) goto 390
601 ELSEIF(mmat.EQ.42.OR.mmat.EQ.43.OR.mmat.EQ.44.OR.mmat.EQ.48)
THEN
608 IF(mbst.EQ.0) wt=four(ip,
n+1)*four(
n+2,
n+
np+1)
609 IF(mbst.EQ.1) wt=
p(ip,5)*
p(
n+1,4)*four(
n+2,
n+
np+1)
610 IF(wt.LT.
pyr(0)*
p(ip,5)*pv(1,5)**3/wtcor(10)) goto 390
614 560
IF(mrem.EQ.1)
THEN
616 pv(1,
j)=pv(1,
j)/(1d0-pqt)
624 IF(mmat.EQ.31.AND.nd.EQ.3)
THEN
629 IF(
p(
n+2,5)**2+
p(
n+3,5)**2+2d0*four(
n+2,
n+3).GE.
630 & (parj(32)+pm2+pm3)**2) goto 630
634 IF(
k(
n+2,2).EQ.0) goto 260
641 ELSEIF(mmat.EQ.44)
THEN
646 IF(
p(
n+3,5)**2+
p(
n+4,5)**2+2d0*four(
n+3,
n+4).GE.
647 & (parj(32)+pm3+pm4)**2) goto 600
651 IF(
k(
n+3,2).EQ.0) goto 260
656 p(
n+3,4)=sqrt(
p(
n+3,1)**2+
p(
n+3,2)**2+
p(
n+3,3)**2+
p(
n+3,5)**2)
657 ha=
p(
n+1,4)**2-
p(
n+2,4)**2
658 hb=ha-(
p(
n+1,5)**2-
p(
n+2,5)**2)
659 hc=(
p(
n+1,1)-
p(
n+2,1))**2+(
p(
n+1,2)-
p(
n+2,2))**2+
660 & (
p(
n+1,3)-
p(
n+2,3))**2
661 hd=(pv(1,4)-
p(
n+3,4))**2
662 he=ha**2-2d0*hd*(
p(
n+1,4)**2+
p(
n+2,4)**2)+hd**2
665 hh=(sqrt(hg**2+he*hf)-hg)/(2d0*hf)
667 pcor=hh*(
p(
n+1,
j)-
p(
n+2,
j))
671 p(
n+1,4)=sqrt(
p(
n+1,1)**2+
p(
n+1,2)**2+
p(
n+1,3)**2+
p(
n+1,5)**2)
672 p(
n+2,4)=sqrt(
p(
n+2,1)**2+
p(
n+2,2)**2+
p(
n+2,3)**2+
p(
n+2,5)**2)
677 600
IF((mmat.EQ.42.OR.mmat.EQ.43.OR.mmat.EQ.44.OR.mmat.EQ.48)
678 &.AND.iabs(
k(
n+1,2)).LT.10)
THEN
679 pmr=sqrt(
max(0d0,
p(
n+1,5)**2+
p(
n+2,5)**2+2d0*four(
n+1,
n+2)))
684 IF(pmr.GT.parj(32)+pm1+pm2) goto 610
685 kfldum=int(1.5d0+
pyr(0))
686 CALL
pykfdi(
k(
n+1,2),-isign(kfldum,
k(
n+1,2)),kfldmp,kf1)
687 CALL
pykfdi(
k(
n+2,2),-isign(kfldum,
k(
n+2,2)),kfldmp,kf2)
688 IF(kf1.EQ.0.OR.kf2.EQ.0) goto 260
690 IF((mmat.EQ.42.OR.mmat.EQ.48).AND.pmr.GT.parj(64)+psm) goto 610
691 IF(mmat.GE.43.AND.pmr.GT.0.2d0*parj(32)+psm) goto 610
692 IF(mmat.EQ.48) goto 390
693 IF(nd.EQ.4.OR.kfa.EQ.15) goto 260
697 IF(
k(
n+1,2).EQ.0) goto 260
702 IF(
ps+parj(64).GT.pv(1,5)) goto 260
710 610
IF((mmat.EQ.42.OR.mmat.EQ.48).AND.iabs(
k(
n+1,2)).LT.10)
THEN
736 be(
j)=
p(ip,
j)/
p(ip,4)
740 bep=be(1)*
p(
i,1)+be(2)*
p(
i,2)+be(3)*
p(
i,3)
742 p(
i,
j)=
p(
i,
j)+ga*(ga*bep/(1d0+ga)+
p(
i,4))*be(
j)
744 p(
i,4)=ga*(
p(
i,4)+bep)
757 IF(mstj(23).GE.1.AND.mmat.EQ.4.AND.
k(nsav+1,2).EQ.21)
THEN
761 k(nsav+1,4)=mstu(5)*(nsav+2)
762 k(nsav+1,5)=mstu(5)*(nsav+3)
763 k(nsav+2,4)=mstu(5)*(nsav+3)
764 k(nsav+2,5)=mstu(5)*(nsav+1)
765 k(nsav+3,4)=mstu(5)*(nsav+1)
766 k(nsav+3,5)=mstu(5)*(nsav+2)
768 ELSEIF(mstj(23).GE.1.AND.mmat.EQ.4)
THEN
771 k(nsav+2,4)=mstu(5)*(nsav+3)
772 k(nsav+2,5)=mstu(5)*(nsav+3)
773 k(nsav+3,4)=mstu(5)*(nsav+2)
774 k(nsav+3,5)=mstu(5)*(nsav+2)
776 ELSEIF(mstj(23).GE.1.AND.(mmat.EQ.32.OR.mmat.EQ.44).AND.
777 & iabs(
k(nsav+1,2)).LE.10.AND.iabs(
k(nsav+2,2)).LE.10)
THEN
780 k(nsav+1,4)=mstu(5)*(nsav+2)
781 k(nsav+1,5)=mstu(5)*(nsav+2)
782 k(nsav+2,4)=mstu(5)*(nsav+1)
783 k(nsav+2,5)=mstu(5)*(nsav+1)
785 ELSEIF(mstj(23).GE.1.AND.(mmat.EQ.32.OR.mmat.EQ.44).AND.
786 & iabs(
k(nsav+1,2)).LE.20.AND.iabs(
k(nsav+2,2)).LE.20)
THEN
788 ELSEIF(mstj(23).GE.1.AND.mmat.EQ.33.AND.iabs(
k(nsav+2,2)).EQ.21)
794 kqp=kchg(kcp,2)*isign(1,
k(nsav+1,2))
797 k(nsav+1,jcon)=mstu(5)*(nsav+2)
798 k(nsav+2,9-jcon)=mstu(5)*(nsav+1)
799 k(nsav+2,jcon)=mstu(5)*(nsav+3)
800 k(nsav+3,9-jcon)=mstu(5)*(nsav+2)
802 ELSEIF(mstj(23).GE.1.AND.mmat.EQ.33)
THEN
805 k(nsav+1,4)=mstu(5)*(nsav+3)
806 k(nsav+1,5)=mstu(5)*(nsav+3)
807 k(nsav+3,4)=mstu(5)*(nsav+1)
808 k(nsav+3,5)=mstu(5)*(nsav+1)
813 IF(
k(ip,1).EQ.5)
k(ip,1)=15
814 IF(
k(ip,1).LE.10)
k(ip,1)=11
815 IF(mmix.EQ.1.AND.mstj(26).EQ.2.AND.
k(ip,1).EQ.11)
k(ip,1)=12