Project CMSSW displayed by LXR CMSSW_13_3_X_2023-09-29-2300/​GeneratorInterface/​ReggeGribovPartonMCInterface/​src/​epos-int.f	Source navigation Diff markup Identifier search General search
	Version: CMSSW_13_3_X_2023-09-29-2300 CMSSW_13_3_X_2023-09-28-2300 CMSSW_13_3_X_2023-09-27-2300 CMSSW_13_3_X_2023-09-26-2300 CMSSW_13_3_X_2023-09-25-2300 CMSSW_13_3_X_2023-09-24-2300 Revert   Change

File indexing completed on 2023-03-17 11:05:08

 c-----------------------------------------------------------------------
       subroutine bjinta(ier)
 c-----------------------------------------------------------------------
 c  fin. state interactions and decays
 c-----------------------------------------------------------------------
       include 'epos.inc'
       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
       common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat /ctimel/ntc
       common/col3/ncol,kolpt
       double precision ttaun,ttau0,rcproj,rctarg
       common/cttaun/ttaun /cttau0/ttau0 /geom1/rcproj,rctarg
       logical go,lclean
 
       call utpri('bjinta',ish,ishini,4)
 
       ier=0
 
       if(ncol.eq.0.and.iappl.eq.2)goto1000
       if(nevt.ne.1.or.ifrade.eq.0)goto1000
 
       if(iappl.eq.4.or.iappl.eq.9)then
         goto5000
       endif
 
       !if(iappl.eq.1)then
       !  tauxx=0.7+0.94*max(radnuc(maproj),radnuc(matarg))/(0.5*engy)
       !else
       !  tauxx=0.
       !endif
       !tauzz=max(taumin,tauxx)
       !print*,'====',taumin,tauxx,tauzz
       !ttaus=dble(tauzz)
       ttaus=taumin
       ttau0=dsqrt(rcproj*rctarg)
       call jtauin   ! initialize hyperbola
 
       if(iappl.ne.1)goto 5000
 
 
 c     no-secondary-interactions or parton-ladder-fusion
 c     -------------------------------------------------
       if(iorsce.eq.0.and.iorsdf.eq.0.and.iorshh.eq.0
      &      .or.iorsdf.eq.3)then
         if(iorsdf.eq.3)then
           lclean=.false.
           if(nclean.gt.0.and.nptl.gt.mxptl/5)then
             ! if nptl already very big, clean up useless particles in cptl list.
             !(do not use it when gakstr() is called (some information lost)
             nptli=maproj+matarg+1
             do iii=nptli,nptl
               go=.true.
               if(nclean.eq.1.and.istptl(iii).le.istmax)go=.false.
               if(go.and.mod(istptl(iii),10).ne.0)istptl(iii)=99
             enddo
             nptl0=nptl
             call utclea(nptli,nptl0)
             lclean=.true.
           endif
           nptlbpo=nptl
           call jintpo(lclean,iret)  !parton-ladder-fusion
          if(ish.ge.2)call alist('parton-ladder-fusion&',nptlbpo+1,nptl)
          if(iret.eq.1)goto 1001
         endif
         goto 5000
       else
         stop'bjinta: not supported any more (310305).     '
       endif
 
 5000  continue
 
       nptlbd=nptl
 
       call xSpaceTime
 
       if(ifrade.eq.0)goto779  !skip decay
       if(idecay.eq.0)goto779  !skip decay
 
 
       if(ish.ge.2)call alist('final decay&',0,0)
       if(iappl.eq.4.or.iappl.eq.7.or.iappl.eq.9)then
         nptli=1
       else
         nptli=maproj+matarg+1
       endif
       np1=nptli
 41    np2=nptl
       nptli=np1
       ip=np1-1
       do while (ip.lt.np2)
       ip=ip+1
       if(istptl(ip).eq.0)then
       call hdecas(ip,iret)
       if(iret.eq.1)goto 1001
       if(iret.eq.-1)goto 42
 c remove useless particles if not enough space
       if(nclean.gt.0.and.nptl.gt.mxptl/2)then
         nnnpt=0
         do iii=nptli,ip
           go=.true.
           if(nclean.eq.1.and.istptl(iii).le.istmax)go=.false.
           if(go.and.mod(istptl(iii),10).ne.0)then
             istptl(iii)=99
             nnnpt=nnnpt+1
           endif
         enddo
         if(nnnpt.gt.mxptl-nptl)then
           nptl0=nptl
           call utclea(nptli,nptl0)
           np2=np2-nnnpt
           ip=ip-nnnpt
           nptli=ip
         endif
       endif
       endif
 42    continue
       enddo
       nptli=max(nptli,np1)
       np1=np2+1
       if(np1.le.nptl)then
       if(ish.ge.2)then
       if(ish.ge.3)call alist('partial list&',0,0)
       do 6 ip=np1,nptl
         call alist('&',ip,ip)
 6     continue
       endif
       goto 41
       endif
   779 continue
 
 c      if(ish.ge.2)call alist('complete list&',1,nptl)
 
 c     on shell check
 c     --------------
 c      if(iappl.eq.1)call jresc
 
 1000  continue
       call utprix('bjinta',ish,ishini,4)
       return
 
 1001  continue
       ier=1
       goto 1000
 
       end
 
 cc----------------------------------------------------------------------
 c      subroutine jintcs(i,j,ecm,bij,nq,jc,ics)
 cc----------------------------------------------------------------------
 cc compare hadron distance with energy dependent cross section
 cc data taken from particle data group, durham and juelich
 cc input:
 cc   i,j: particle indices
 cc   ecm: center-of-mass energy
 cc   bij: impact parameter
 cc   nq: net quark number of fused object
 cc   jc: jc of fused object
 cc output:
 cc   ics=0 if distance larger than sqrt(sig(E_CMS)/pi)
 cc   ics=1 else
 cc The data are from HEPDATA,
 cc the formulas from Rev. Particle Properties 1995
 cc----------------------------------------------------------------------
 c      include 'epos.inc'
 c      integer jci(nflav,2),jcj(nflav,2),jc(nflav,2),kc(nflav)
 c     *,kci(nflav),kcj(nflav)
 c      common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
 c     *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
 c      parameter(npp=249,napp=205,npn=411,napn=31,npip=441)
 c      parameter(npim=578,nkmp=299,nkmn=41,nkpp=172,nkpn=91)
 cc      parameter(npim=578,nkmp=299,nkmn=41,nkpp=172,nkpn=91,nlp=35)
 c      parameter(npi1=12,npi2=12,npi3=18,npi4=21,npi5=9)
 cc      real ppecm(npp)
 c      real ppbmx(npp)
 c      real appecm(napp),appbmx(napp)
 c      real pnecm(npn),pnbmx(npn)
 c      real apnecm(napn),apnbmx(napn)
 c      real pipecm(npip),pipbmx(npip)
 c      real pimecm(npim),pimbmx(npim)
 c      real kmpecm(nkmp),kmpbmx(nkmp)
 c      real kmnecm(nkmn),kmnbmx(nkmn)
 c      real kppecm(nkpp),kppbmx(nkpp)
 c      real kpnecm(nkpn),kpnbmx(nkpn)
 cc      real lpecm(nlp),lpbmx(nlp)
 c      real pi1ecm(npi1),pi1bmx(npi1)
 c      real pi2ecm(npi2),pi2bmx(npi2)
 c      real pi3ecm(npi3),pi3bmx(npi3)
 c      real pi4ecm(npi4),pi4bmx(npi4)
 c      real pi5ecm(npi5),pi5bmx(npi5)
 c
 cc      data ppecm/
 cc     *    1.8812,   1.8855,   1.8910,   1.8963,   1.9073
 cc     *,   1.9108,   1.9145,   1.9224,   1.9244,   1.9352
 cc     *,   1.9466,   1.9468,   1.9542,   1.9592,   1.9636
 cc     *,   1.9772,   1.9860,   1.9945,   2.0032,   2.0052
 cc     *,   2.0070,   2.0272,   2.0275,   2.0302,   2.0333
 cc     *,   2.0402,   2.0427,   2.0586,   2.0608,   2.0692
 cc     *,   2.0702,   2.0708,   2.0715,   2.0718,   2.0751
 cc     *,   2.0797,   2.0813,   2.0829,   2.0843,   2.0846
 cc     *,   2.0935,   2.1062,   2.1113,   2.1123,   2.1170
 cc     *,   2.1173,   2.1184,   2.1268,   2.1289,   2.1357
 cc     *,   2.1467,   2.1511,   2.1522,   2.1553,   2.1618
 cc     *,   2.1639,   2.1726,   2.1771,   2.1795,   2.1799
 cc     *,   2.1802,   2.1813,   2.1868,   2.2152,   2.2184
 cc     *,   2.2212,   2.2254,   2.2395,   2.2405,   2.2532
 cc     *,   2.2606,   2.2613,   2.2861,   2.2889,   2.2914
 cc     *,   2.2988,   2.3101,   2.3136,   2.3228,   2.3348
 cc     *,   2.3490,   2.3525,   2.3529,   2.3800,   2.3842
 cc     *,   2.3912,   2.4088,   2.4130,   2.4193,   2.4298
 cc     *,   2.4315,   2.4392,   2.4472,   2.4573,   2.5034
 cc     *,   2.5131,   2.5268,   2.5743,   2.5848,   2.5916
 cc     *,   2.6327,   2.6620,   2.6700,   2.6984,   2.7080
 cc     *,   2.7205,   2.7534,   2.7573,   2.7651,   2.7670
 cc     *,   2.7844,   2.8024,   2.8092,   2.8127,   2.8533
 cc     *,   2.8556,   2.8638,   2.9079,   2.9395,   2.9500
 cc     *,   2.9776,   2.9961,   3.0495,   3.0769,   3.0879
 cc     *,   3.1358,   3.1547,   3.2251,   3.2371,   3.3020
 cc     *,   3.3527,   3.3620,   3.4221,   3.4967,   3.5019
 cc     *,   3.5035,   3.5814,   3.5829,   3.6266,   3.8549
 cc     *,   3.8742,   4.0503,   4.0663,   4.0698,   4.0732
 cc     *,   4.0778,   4.1074,   4.1300,   4.4976,   4.5183
 cc     *,   4.5389,   4.5410,   4.5615,   4.6808,   4.9146
 cc     *,   4.9336,   5.0088,   5.2993,   5.3011,   5.4731
 cc     *,   5.6083,   5.6416,   5.6465,   5.9171,   5.9502
 cc     *,   5.9644,   6.1697,   6.1803,   6.2706,   6.3034
 cc     *,   6.3390,   6.5627,   6.7040,   6.8424,   6.9105
 cc     *,   6.9780,   7.1111,   7.1662,   7.6202,   7.8624
 cc     *,   8.2124,   8.7647,   9.0282,   9.2843,   9.5825
 cc     *,   9.7763,   9.9851,  10.2447,  10.6927,  10.8926
 cc     *,  11.4549,  11.5365,  11.7779,  11.8519,  13.6241
 cc     *,  13.6883,  13.7611,  13.8968,  15.0628,  15.1868
 cc     *,  16.6595,  16.8275,  17.9077,  18.0121,  18.1677
 cc     *,  19.2213,  19.4156,  19.6556,  19.7002,  21.2604
 cc     *,  22.9574,  23.3624,  23.4057,  23.4965,  23.4967
 cc     *,  23.5964,  23.7605,  23.8787,  24.1521,  25.2904
 cc     *,  26.3796,  27.5960,  30.5240,  30.5954,  30.5957
 cc     *,  30.6555,  30.6954,  30.7954,  35.1947,  44.6933
 cc     *,  44.6937,  44.6938,  44.7706,  44.8228,  44.8933
 cc     *,  45.1933,  52.6798,  52.7090,  52.7921,  52.7921
 cc     *,  52.7927,  52.8927,  53.1921,  62.2907,  62.3914
 cc     *,  62.4907,  62.6907,  62.6913,  62.7906
 cc     */
 c      data ppbmx/
 c     *    3.1615,   2.2212,   1.7113,   1.4927,   1.1631
 c     *,   1.0911,   1.0388,    .9525,    .9390,    .8885
 c     *,    .8019,    .8956,    .9115,    .8593,    .8840
 c     *,    .9062,    .8444,    .8444,    .8482,    .8713
 c     *,    .8575,    .8795,    .8795,    .8759,    .8630
 c     *,    .8822,    .8593,    .8813,    .9036,    .9150
 c     *,    .8740,    .9219,    .8740,    .9253,    .8722
 c     *,    .8795,    .9288,    .8704,    .9398,    .9271
 c     *,    .9407,    .9373,    .9756,    .9926,   1.0357
 c     *,   1.0037,   1.0408,    .9739,   1.0108,   1.0479
 c     *,   1.0645,   1.0705,   1.1113,   1.0794,   1.0955
 c     *,   1.1085,   1.1256,   1.1535,   1.1731,   1.1424
 c     *,   1.0794,   1.1480,   1.1590,   1.1888,   1.2218
 c     *,   1.2164,   1.2087,   1.2296,   1.2231,   1.2335
 c     *,   1.2322,   1.2309,   1.2114,   1.2296,   1.2293
 c     *,   1.2489,   1.2303,   1.2322,   1.2322,   1.2127
 c     *,   1.2192,   1.2295,   1.2399,   1.2290,   1.2374
 c     *,   1.2205,   1.2278,   1.2284,   1.2061,   1.2296
 c     *,   1.2296,   1.2540,   1.2008,   1.2260,   1.2229
 c     *,   1.2257,   1.2188,   1.2118,   1.2078,   1.1982
 c     *,   1.2039,   1.2012,   1.1991,   1.1808,   1.1969
 c     *,   1.1959,   1.1922,   1.1902,   1.1897,   1.1878
 c     *,   1.1888,   1.1860,   1.1856,   1.1850,   1.2244
 c     *,   1.1782,   1.1790,   1.1718,   1.1696,   1.1726
 c     *,   1.1576,   1.1656,   1.1606,   1.1603,   1.1581
 c     *,   1.1590,   1.1530,   1.1576,   1.1487,   1.1476
 c     *,   1.1447,   1.1794,   1.1448,   1.1507,   1.1507
 c     *,   1.1407,   1.1403,   1.1507,   1.1581,   1.1645
 c     *,   1.1616,   1.1507,   1.1332,   1.1294,   1.1284
 c     *,   1.1227,   1.1284,   1.1298,   1.1261,   1.1199
 c     *,   1.1365,   1.1438,   1.1284,   1.1424,   1.1230
 c     *,   1.1218,   1.1142,   1.1156,   1.1202,   1.1198
 c     *,   1.1099,   1.1099,   1.1175,   1.1241,   1.1168
 c     *,   1.1099,   1.1128,   1.1241,   1.1103,   1.1149
 c     *,   1.1155,   1.1083,   1.1197,   1.1127,   1.1185
 c     *,   1.1113,   1.1128,   1.1113,   1.1083,   1.1056
 c     *,   1.1067,   1.1070,   1.1059,   1.1063,   1.1070
 c     *,   1.1028,   1.0979,   1.1060,   1.1062,   1.1774
 c     *,   1.0805,   1.1039,   1.0940,   1.0955,   1.0940
 c     *,   1.0955,   1.1082,   1.1128,   1.1082,   1.0852
 c     *,   1.1003,   1.1097,   1.1118,   1.0911,   1.1227
 c     *,   1.1070,   1.1206,   1.1142,   1.1213,   1.1174
 c     *,   1.1199,   1.1142,   1.1185,   1.1180,   1.1113
 c     *,   1.1099,   1.1297,   1.1142,   1.1226,   1.1240
 c     *,   1.1251,   1.1368,   1.1403,   1.1319,   1.1296
 c     *,   1.1333,   1.1303,   1.1284,   1.1343,   1.1516
 c     *,   1.1521,   1.1549,   1.1641,   1.1631,   1.1562
 c     *,   1.1631,   1.1697,   1.1726,   1.1756,   1.1659
 c     *,   1.1660,   1.1617,   1.1686,   1.1774,   1.1713
 c     *,   1.1744,   1.1810,   1.1694,   1.1848
 c     */
 c      data appecm/
 c     *    1.9002,   1.9050,   1.9072,   1.9078,   1.9091
 c     *,   1.9129,   1.9157,   1.9162,   1.9174,   1.9176
 c     *,   1.9180,   1.9195,   1.9201,   1.9224,   1.9226
 c     *,   1.9246,   1.9252,   1.9255,   1.9257,   1.9271
 c     *,   1.9282,   1.9293,   1.9301,   1.9310,   1.9319
 c     *,   1.9328,   1.9334,   1.9345,   1.9359,   1.9370
 c     *,   1.9372,   1.9384,   1.9393,   1.9398,   1.9407
 c     *,   1.9426,   1.9430,   1.9433,   1.9452,   1.9454
 c     *,   1.9473,   1.9485,   1.9495,   1.9500,   1.9510
 c     *,   1.9515,   1.9547,   1.9559,   1.9562,   1.9579
 c     *,   1.9610,   1.9615,   1.9644,   1.9680,   1.9718
 c     *,   1.9755,   1.9788,   1.9829,   1.9871,   1.9911
 c     *,   1.9954,   1.9994,   2.0813,   2.0979,   2.1146
 c     *,   2.1180,   2.1316,   2.1487,   2.1660,   2.1834
 c     *,   2.1868,   2.1938,   2.1991,   2.2008,   2.2184
 c     *,   2.2226,   2.2359,   2.2500,   2.2606,   2.2712
 c     *,   2.2889,   2.2995,   2.3066,   2.3243,   2.3419
 c     *,   2.3490,   2.3511,   2.3596,   2.3701,   2.3772
 c     *,   2.3860,   2.3877,   2.3947,   2.4035,   2.4123
 c     *,   2.4298,   2.4472,   2.4629,   2.4820,   2.4993
 c     *,   2.5165,   2.5268,   2.5337,   2.5508,   2.5678
 c     *,   2.5848,   2.5950,   2.6017,   2.6186,   2.6353
 c     *,   2.6377,   2.6520,   2.6654,   2.6687,   2.6852
 c     *,   2.7017,   2.7182,   2.7345,   2.7508,   2.7670
 c     *,   2.7832,   2.7896,   2.7992,   2.8152,   2.8312
 c     *,   2.8439,   2.8470,   2.8565,   2.8628,   2.9377
 c     *,   2.9561,   2.9745,   3.0351,   3.0769,   3.0828
 c     *,   3.1648,   3.2507,   3.2788,   3.3620,   3.4568
 c     *,   3.5492,   3.6266,   3.7893,   3.8501,   3.8597
 c     *,   3.8742,   3.9455,   4.1074,   4.3499,   4.3926
 c     *,   4.5389,   4.5799,   4.6808,   4.7991,   4.9336
 c     *,   4.9901,   5.1742,   5.2993,   5.3520,   5.4731
 c     *,   5.6416,   5.6911,   5.8534,   5.9644,   6.0113
 c     *,   6.2706,   6.3153,   6.7040,   6.9780,   7.3061
 c     *,   7.6202,   7.7422,   7.8624,   7.8743,   8.0393
 c     *,   8.2124,   8.4930,   8.7647,   9.0282,   9.2843
 c     *,   9.5335,   9.7763,  11.5365,  13.7611,  13.8630
 c     *,  15.0628,  16.8275,  17.9077,  19.4156,  21.2604
 c     *,  22.9574,  30.4098,  30.5943,  30.6861,  52.5843
 c     *,  52.7979,  52.7979,  62.2853,  62.4957,  62.6905
 c     *, 539.9198, 546.9191, 899.8658, 900.0000,1803.0007
 c     */
 c      data appbmx/
 c     *    2.7553,   2.6857,   2.6732,   2.6517,   2.6306
 c     *,   2.5376,   2.4398,   2.4779,   2.5212,   2.5143
 c     *,   2.4573,   2.4799,   2.4495,   2.4286,   2.4521
 c     *,   2.3683,   2.4463,   2.4240,   2.4489,   2.3743
 c     *,   2.4162,   2.3561,   2.3903,   2.3796,   2.3221
 c     *,   2.3689,   2.4181,   2.3221,   2.3823,   2.2987
 c     *,   2.3337,   2.3568,   2.3487,   2.2911,   2.3344
 c     *,   2.2701,   2.3097,   2.3207,   2.2631,   2.2603
 c     *,   2.2666,   2.2412,   2.2319,   2.3097,   2.2327
 c     *,   2.2104,   2.2104,   2.2645,   2.2369,   2.1880
 c     *,   2.1756,   2.2162,   2.1491,   2.1484,   2.1365
 c     *,   2.1320,   2.1027,   2.1065,   2.0913,   2.0776
 c     *,   2.0668,   2.0883,   1.9333,   1.9098,   1.8851
 c     *,   1.8932,   1.8721,   1.8623,   1.8520,   1.8409
 c     *,   1.8678,   1.8575,   1.8325,   1.8429,   1.8088
 c     *,   1.7698,   1.7941,   1.7864,   1.7814,   1.7736
 c     *,   1.7645,   1.7623,   1.7576,   1.7523,   1.7445
 c     *,   1.7419,   1.7019,   1.7342,   1.7311,   1.7271
 c     *,   1.7210,   1.7182,   1.7161,   1.7119,   1.7054
 c     *,   1.6947,   1.6816,   1.6780,   1.6678,   1.6599
 c     *,   1.6509,   1.6328,   1.6449,   1.6396,   1.6319
 c     *,   1.6270,   1.5761,   1.6043,   1.6120,   1.6069
 c     *,   1.5978,   1.6018,   1.5891,   1.5948,   1.5905
 c     *,   1.5849,   1.5778,   1.5736,   1.5679,   1.5785
 c     *,   1.5583,   1.5476,   1.5519,   1.5467,   1.5416
 c     *,   1.5233,   1.5368,   1.5492,   1.5313,   1.4895
 c     *,   1.5579,   1.5107,   1.4680,   1.4725,   1.4586
 c     *,   1.3889,   1.4799,   1.4472,   1.4371,   1.3576
 c     *,   1.4228,   1.3922,   1.3762,   1.3854,   1.3669
 c     *,   1.4161,   1.3623,   1.3638,   1.3530,   1.3762
 c     *,   1.3273,   1.2866,   1.2989,   1.3159,   1.2828
 c     *,   1.3086,   1.2797,   1.2704,   1.2940,   1.2741
 c     *,   1.2514,   1.2540,   1.2803,   1.2653,   1.2603
 c     *,   1.2386,   1.2192,   1.2283,   1.2231,   1.2140
 c     *,   1.2048,   1.2114,   1.2099,   1.2052,   1.2048
 c     *,   1.1984,   1.1928,   1.1901,   1.1902,   1.1888
 c     *,   1.1848,   1.1769,   1.1706,   1.1579,   1.1608
 c     *,   1.1521,   1.1534,   1.1520,   1.1495,   1.1549
 c     *,   1.1550,   1.1580,   1.1672,   1.1562,   1.1743
 c     *,   1.1859,   1.1950,   1.1851,   1.1821,   1.1987
 c     *,   1.4740,   1.4037,   1.4483,   1.4417,   1.5149
 c     */
 c      data pnecm/
 c     *    1.87867,   1.87869,   1.87872,   1.87875,   1.87878
 c     *,   1.87881,   1.87884,   1.87887,   1.87890,   1.87893
 c     *,   1.87896,   1.87899,   1.87902,   1.87906,   1.87909
 c     *,   1.87913,   1.87916,   1.87920,   1.87923,   1.87927
 c     *,   1.87931,   1.87934,   1.87938,   1.87942,   1.87946
 c     *,   1.87950,   1.87954,   1.87958,   1.87962,   1.87966
 c     *,   1.87970,   1.87975,   1.87979,   1.87983,   1.87988
 c     *,   1.87992,   1.87997,   1.88001,   1.88006,   1.88011
 c     *,   1.88015,   1.88020,   1.88025,   1.88030,   1.88035
 c     *,   1.88040,   1.88045,   1.88050,   1.88055,   1.88061
 c     *,   1.88066,   1.88071,   1.88077,   1.88082,   1.88087
 c     *,   1.88093,   1.88099,   1.88104,   1.88110,   1.88116
 c     *,   1.88121,   1.88127,   1.88133,   1.88139,   1.88145
 c     *,   1.88151,   1.88157,   1.88163,   1.88170,   1.88176
 c     *,   1.88182,   1.88189,   1.88195,   1.88202,   1.88208
 c     *,   1.88215,   1.88221,   1.88228,   1.88235,   1.88241
 c     *,   1.88248,   1.88255,   1.88262,   1.88269,   1.88276
 c     *,   1.88283,   1.88290,   1.88297,   1.88305,   1.88312
 c     *,   1.88319,   1.88327,   1.88334,   1.88342,   1.88349
 c     *,   1.88357,   1.88364,   1.88372,   1.88380,   1.88388
 c     *,   1.88396,   1.88403,   1.88411,   1.88419,   1.88428
 c     *,   1.88436,   1.88444,   1.88452,   1.88460,   1.88469
 c     *,   1.88477,   1.88485,   1.88494,   1.88502,   1.88511
 c     *,   1.88519,   1.88528,   1.88537,   1.88546,   1.88554
 c     *,   1.88563,   1.88572,   1.88581,   1.88590,   1.88599
 c     *,   1.88608,   1.88618,   1.88627,   1.88636,   1.88645
 c     *,   1.88655,   1.88664,   1.88674,   1.88683,   1.88693
 c     *,   1.88702,   1.88712,   1.88722,   1.88731,   1.88741
 c     *,   1.88751,   1.88761,   1.88771,   1.88781,   1.88791
 c     *,   1.88801,   1.88811,   1.88822,   1.88832,   1.88842
 c     *,   1.88852,   1.88863,   1.88873,   1.88884,   1.88894
 c     *,   1.88905,   1.88916,   1.88926,   1.88937,   1.88948
 c     *,   1.88959,   1.88970,   1.88980,   1.88991,   1.89003
 c     *,   1.89014,   1.89025,   1.89036,   1.89047,   1.89058
 c     *,   1.89070,   1.89081,   1.89093,   1.89104,   1.89116
 c     *,   1.89127,   1.89139,   1.89150,   1.89162,   1.89174
 c     *,   1.89186,   1.89198,   1.89209,   1.89221,   1.89233
 c     *,   1.89245,   1.89258,   1.89270,   1.89282,   1.89294
 c     *,   1.89306,   1.89319,   1.89331,   1.89344,   1.89356
 c     *,   1.89369,   1.89381,   1.89394,   1.89406,   1.89419
 c     *,   1.89432,   1.89445,   1.89458,   1.89496,   1.89549
 c     *,   1.89602,   1.89615,   1.89656,   1.89697,   1.89724
 c     *,   1.89766,   1.89780,   1.89850,   1.89893,   1.89922
 c     *,   1.89995,   1.90068,   1.90098,   1.90143,   1.90174
 c     *,   1.90235,   1.90250,   1.90312,   1.90407,   1.90503
 c     *,   1.90519,   1.90616,   1.90732,   1.90749,   1.90833
 c     *,   1.90936,   1.91216,   1.91379,   1.91545,   1.91714
 c     *,   1.91905,   1.91924,   1.92100,   1.92160,   1.92179
 c     *,   1.92259,   1.92319,   1.92421,   1.92502,   1.92543
 c     *,   1.92605,   1.92646,   1.92792,   1.93046,   1.93241
 c     *,   1.93306,   1.93570,   1.93704,   1.93953,   1.94021
 c     *,   1.94183,   1.94699,   1.94723,   1.95206,   1.95329
 c     *,   1.95452,   1.95651,   1.96029,   1.96080,   1.97927
 c     *,   1.98533,   1.99806,   2.00884,   2.01269,   2.02779
 c     *,   2.02963,   2.04329,   2.04643,   2.05915,   2.05947
 c     *,   2.05979,   2.07207,   2.07337,   2.07533,   2.09178
 c     *,   2.10847,   2.11385,   2.12061,   2.12536,   2.14243
 c     *,   2.15308,   2.15343,   2.15964,   2.17072,   2.17698
 c     *,   2.18185,   2.19442,   2.21193,   2.21369,   2.22353
 c     *,   2.22846,   2.22952,   2.24715,   2.25597,   2.28851
 c     *,   2.29134,   2.29382,   2.31293,   2.31328,   2.32673
 c     *,   2.34935,   2.35501,   2.36207,   2.38251,   2.39729
 c     *,   2.41135,   2.41591,   2.42817,   2.44006,   2.45995
 c     *,   2.48220,   2.49952,   2.50643,   2.52951,   2.56750
 c     *,   2.58756,   2.63546,   2.64719,   2.66487,   2.67285
 c     *,   2.71088,   2.72336,   2.75178,   2.75634,   2.76802
 c     *,   2.76867,   2.76996,   2.78741,   2.80542,   2.81567
 c     *,   2.85639,   2.85860,   2.86681,   2.91102,   2.94265
 c     *,   2.96543,   3.05273,   3.08729,   3.09115,   3.15805
 c     *,   3.16821,   3.22857,   3.24053,   3.29551,   3.35628
 c     *,   3.42580,   3.50727,   3.55297,   3.58520,   3.58675
 c     *,   3.59580,   3.63049,   3.72425,   3.75633,   4.05460
 c     *,   4.06219,   4.07412,   4.11175,   4.37426,   4.52314
 c     *,   4.54378,   4.74537,   4.93885,   5.12513,   5.30495
 c     *,   5.40999,   5.47892,   5.61424,   5.64757,   5.93444
 c     *,   5.97071,   6.27732,   6.42517,   6.51227,   6.56970
 c     *,   6.71114,   6.87703,   6.98545,   7.18446,   7.23904
 c     *,   7.24942,   7.62830,   8.10603,   8.22114,   8.54948
 c     *,   8.77406,   9.29418,   9.78673,  10.15712,  10.25566
 c     *,  10.70407,  11.13446,  11.54882,  12.33587,  13.77577
 c     *,  15.07880,  15.74960,  16.84544,  17.92675,  18.18703
 c     *,  18.44365,  19.43624,  20.14863,  21.28302,  22.69375
 c     *,  22.98187
 c     */
 c      data pnbmx/
 c     *   10.7613,  10.6195,  10.5290,  10.4010,  10.2901
 c     *,  10.1970,  10.1029,  10.0088,   9.9009,   9.7973
 c     *,   9.7125,   9.6198,   9.5613,   9.4281,   9.3628
 c     *,   9.2801,   9.2068,   9.1103,   9.0239,   8.9609
 c     *,   8.8562,   8.8011,   8.7574,   8.6332,   8.5487
 c     *,   8.4788,   8.4157,   8.3606,   8.2746,   8.2604
 c     *,   8.1583,   8.0778,   8.0127,   7.9567,   7.8691
 c     *,   7.8176,   7.8190,   7.7028,   7.6634,   7.5956
 c     *,   7.5808,   7.4956,   7.4144,   7.3600,   7.3071
 c     *,   7.2947,   7.1947,   7.1389,   7.1333,   7.0404
 c     *,   6.9856,   6.9468,   6.9171,   6.8369,   6.7869
 c     *,   6.7829,   6.6804,   6.6483,   6.5982,   6.5464
 c     *,   6.5060,   6.4541,   6.4539,   6.3684,   6.3219
 c     *,   6.2777,   6.2832,   6.1881,   6.1494,   6.1048
 c     *,   6.0642,   6.0207,   6.0255,   5.9756,   5.9161
 c     *,   5.8717,   5.8276,   5.7878,   5.7481,   5.7319
 c     *,   5.6977,   5.6492,   5.6142,   5.5796,   5.5428
 c     *,   5.4939,   5.4659,   5.4549,   5.3933,   5.3601
 c     *,   5.3311,   5.3037,   5.2760,   5.2429,   5.2381
 c     *,   5.2384,   5.1470,   5.1112,   5.0812,   5.0483
 c     *,   5.0350,   5.0330,   4.9640,   4.9410,   4.9050
 c     *,   4.8714,   4.8408,   4.8482,   4.8518,   4.7611
 c     *,   4.7157,   4.7049,   4.6834,   4.6939,   4.6206
 c     *,   4.5931,   4.5723,   4.5498,   4.5210,   4.5337
 c     *,   4.4726,   4.4898,   4.4251,   4.4013,   4.3784
 c     *,   4.3514,   4.3264,   4.3530,   4.2846,   4.2567
 c     *,   4.2339,   4.2171,   4.1919,   4.1757,   4.1503
 c     *,   4.1329,   4.0924,   4.0797,   4.0624,   4.0446
 c     *,   4.0210,   4.0010,   3.9812,   3.9627,   3.9402
 c     *,   3.9221,   3.9029,   3.8925,   3.8680,   3.8439
 c     *,   3.8244,   3.8094,   3.7816,   3.7719,   3.7396
 c     *,   3.7271,   3.7011,   3.6849,   3.6804,   3.6640
 c     *,   3.6496,   3.5954,   3.6042,   3.5697,   3.5552
 c     *,   3.5463,   3.5265,   3.5212,   3.4992,   3.1432
 c     *,   3.4630,   3.4483,   3.4484,   3.4240,   3.4127
 c     *,   3.3595,   3.3767,   3.3582,   3.3059,   3.3258
 c     *,   3.3182,   3.2582,   3.2740,   3.2621,   3.1984
 c     *,   3.2408,   3.2043,   3.2087,   3.1652,   3.1720
 c     *,   3.1476,   3.0754,   3.1236,   3.1494,   3.0653
 c     *,   3.0757,   3.1030,   3.0692,   3.0466,   3.0392
 c     *,   3.0346,   2.9933,   3.0146,   2.8774,   2.7972
 c     *,   2.6463,   2.8068,   2.6869,   2.6643,   2.6673
 c     *,   2.6684,   2.6409,   2.5520,   2.4341,   2.4978
 c     *,   2.4534,   2.4417,   2.4293,   2.3283,   2.4978
 c     *,   2.3042,   2.3111,   2.1996,   2.1996,   2.1960
 c     *,   2.1469,   2.1223,   2.0791,   2.0599,   1.9891
 c     *,   2.0027,   1.8678,   1.7850,   1.7868,   1.7362
 c     *,   1.6478,   1.6254,   1.6166,   1.6478,   1.6516
 c     *,   1.5554,   1.5727,   1.5656,   1.5290,   1.5149
 c     *,   1.5554,   1.5615,   1.5522,   1.4863,   1.4745
 c     *,   1.4745,   1.3991,   1.4139,   1.1686,   1.3458
 c     *,   1.3587,   1.2425,   1.3171,   1.2828,   1.2153
 c     *,   1.2679,   1.2766,   1.1410,   1.2514,   1.0852
 c     *,   1.1445,   1.1672,   1.0624,   1.1089,   1.0899
 c     *,   1.0171,   1.0399,   1.0645,   1.0417,    .9934
 c     *,   1.0403,   1.0029,   1.0357,   1.0388,   1.0337
 c     *,   1.0428,   1.0555,   1.0663,   1.0546,   1.0626
 c     *,   1.0013,   1.0705,   1.0711,   1.0852,   1.0830
 c     *,   1.1119,   1.0940,   1.0976,   1.0675,   1.1205
 c     *,   1.0464,   1.0995,    .9508,   1.0972,   1.1170
 c     *,   1.1015,   1.1251,   1.1296,   1.1027,   1.1028
 c     *,    .9271,   1.1363,   1.1120,   1.1455,   1.1192
 c     *,   1.1498,   1.1491,   1.0108,   1.1282,   1.1550
 c     *,   1.1617,   1.1354,   1.1586,   1.1631,   1.1378
 c     *,   1.1656,   1.1684,   1.1389,   1.1694,   1.1690
 c     *,   1.1702,   1.1705,   1.1390,   1.1714,   1.1716
 c     *,   1.1326,   1.1517,   1.1697,   1.1731,   1.1716
 c     *,   1.0867,   1.1537,   1.1698,   1.1642,   1.1638
 c     *,   1.1199,   1.1635,   1.1713,   1.1631,   1.1601
 c     *,   1.1340,   1.0823,   1.1598,   1.1754,   1.1572
 c     *,   1.1565,   1.1567,   1.1631,   1.1538,   1.0852
 c     *,   1.0342,   1.1645,   1.1452,   1.1099,   1.0940
 c     *,   1.1185,   1.1470,   1.1388,   1.1424,   1.0705
 c     *,   1.1374,   1.1199,   1.1262,   1.1142,   1.1205
 c     *,   1.0867,   1.1095,   1.0734,   1.1234,   1.0925
 c     *,   1.1135,   1.1076,   1.1070,   1.0955,   1.1027
 c     *,   1.1073,   1.0630,   1.1007,   1.1185,   1.1133
 c     *,   1.1128,   1.1028,   1.1033,   1.1086,   1.1120
 c     *,   1.1013,   1.0991,   1.1086,   1.1028,   1.1027
 c     *,   1.1064,   1.1036,   1.1135,   1.1139,   1.1136
 c     *,   1.1145,   1.1166,   1.1152,   1.1168,   1.1172
 c     *,   1.1230,   1.1187,   1.1254,   1.1224,   1.1329
 c     *,   1.1216
 c     */
 c      data apnecm/
 c     *    2.1288,   2.2007,   2.2500,   2.3044,   2.3529
 c     *,   2.9284,   3.5136,   3.6305,   3.7933,   4.1117
 c     *,   4.5438,   4.9389,   5.1797,   5.3049,   5.4789
 c     *,   5.6476,   6.2773,   6.9855,   7.6283,   8.2211
 c     *,   8.7741,   9.2942,   9.7867,  11.5488,  13.7758
 c     *,  15.0788,  16.8454,  17.9267,  19.4362,  21.2830
 c     *,  22.9819
 c     */
 c      data apnbmx/
 c     *    1.9462,   1.7481,   1.8881,   1.8019,   1.8627
 c     *,   1.4745,   1.3231,   1.3762,   1.3351,   1.3505
 c     *,   1.2841,   1.3086,   1.2565,   1.3038,   1.2218
 c     *,   1.2952,   1.2127,   1.2035,   1.1968,   1.1951
 c     *,   1.1604,   1.1794,   1.1747,   1.1724,   1.1595
 c     *,   1.1637,   1.1477,   1.1520,   1.1510,   1.1478
 c     *,   1.1466
 c     */
 c      data pipecm/
 c     *    1.1050,   1.1154,   1.1165,   1.1256,   1.1273
 c     *,   1.1333,   1.1370,   1.1382,   1.1394,   1.1438
 c     *,   1.1495,   1.1579,   1.1592,   1.1677,   1.1691
 c     *,   1.1697,   1.1757,   1.1771,   1.1777,   1.1784
 c     *,   1.1798,   1.1838,   1.1892,   1.1906,   1.1926
 c     *,   1.1933,   1.1953,   1.1960,   1.1967,   1.1981
 c     *,   1.1994,   1.2008,   1.2015,   1.2022,   1.2028
 c     *,   1.2056,   1.2063,   1.2069,   1.2097,   1.2104
 c     *,   1.2111,   1.2124,   1.2131,   1.2138,   1.2152
 c     *,   1.2166,   1.2172,   1.2193,   1.2200,   1.2214
 c     *,   1.2221,   1.2255,   1.2262,   1.2269,   1.2276
 c     *,   1.2283,   1.2303,   1.2310,   1.2317,   1.2352
 c     *,   1.2358,   1.2365,   1.2372,   1.2400,   1.2407
 c     *,   1.2421,   1.2434,   1.2462,   1.2476,   1.2503
 c     *,   1.2517,   1.2524,   1.2538,   1.2545,   1.2565
 c     *,   1.2613,   1.2627,   1.2696,   1.2716,   1.2730
 c     *,   1.2751,   1.2799,   1.2819,   1.2833,   1.2860
 c     *,   1.2867,   1.2915,   1.2922,   1.2929,   1.2990
 c     *,   1.3010,   1.3024,   1.3119,   1.3126,   1.3180
 c     *,   1.3200,   1.3220,   1.3254,   1.3261,   1.3321
 c     *,   1.3382,   1.3388,   1.3415,   1.3422,   1.3449
 c     *,   1.3469,   1.3522,   1.3608,   1.3615,   1.3622
 c     *,   1.3655,   1.3661,   1.3767,   1.3786,   1.3813
 c     *,   1.3898,   1.3917,   1.3950,   1.4047,   1.4080
 c     *,   1.4112,   1.4157,   1.4163,   1.4176,   1.4208
 c     *,   1.4215,   1.4285,   1.4298,   1.4304,   1.4336
 c     *,   1.4349,   1.4425,   1.4470,   1.4495,   1.4526
 c     *,   1.4539,   1.4652,   1.4677,   1.4702,   1.4764
 c     *,   1.4777,   1.4857,   1.4863,   1.4870,   1.4882
 c     *,   1.4919,   1.4993,   1.4999,   1.5030,   1.5097
 c     *,   1.5103,   1.5121,   1.5146,   1.5212,   1.5249
 c     *,   1.5285,   1.5327,   1.5351,   1.5429,   1.5435
 c     *,   1.5465,   1.5513,   1.5531,   1.5548,   1.5566
 c     *,   1.5637,   1.5655,   1.5673,   1.5708,   1.5720
 c     *,   1.5732,   1.5761,   1.5790,   1.5849,   1.5866
 c     *,   1.5965,   1.5977,   1.5994,   1.6023,   1.6063
 c     *,   1.6121,   1.6133,   1.6144,   1.6264,   1.6287
 c     *,   1.6327,   1.6344,   1.6406,   1.6417,   1.6429
 c     *,   1.6502,   1.6519,   1.6536,   1.6625,   1.6664
 c     *,   1.6670,   1.6676,   1.6687,   1.6692,   1.6714
 c     *,   1.6792,   1.6825,   1.6853,   1.6913,   1.6935
 c     *,   1.6941,   1.6963,   1.6985,   1.6990,   1.7012
 c     *,   1.7170,   1.7175,   1.7208,   1.7240,   1.7267
 c     *,   1.7391,   1.7417,   1.7423,   1.7428,   1.7503
 c     *,   1.7535,   1.7609,   1.7635,   1.7651,   1.7714
 c     *,   1.7735,   1.7767,   1.7793,   1.7798,   1.7908
 c     *,   1.7923,   1.7949,   1.7960,   1.8032,   1.8063
 c     *,   1.8099,   1.8197,   1.8207,   1.8238,   1.8253
 c     *,   1.8289,   1.8304,   1.8320,   1.8386,   1.8436
 c     *,   1.8497,   1.8507,   1.8537,   1.8643,   1.8673
 c     *,   1.8678,   1.8703,   1.8713,   1.8738,   1.8777
 c     *,   1.8792,   1.8812,   1.8827,   1.8881,   1.8911
 c     *,   1.8916,   1.8936,   1.8980,   1.9039,   1.9073
 c     *,   1.9108,   1.9151,   1.9161,   1.9181,   1.9186
 c     *,   1.9249,   1.9283,   1.9317,   1.9331,   1.9350
 c     *,   1.9384,   1.9423,   1.9519,   1.9571,   1.9614
 c     *,   1.9652,   1.9681,   1.9756,   1.9780,   1.9804
 c     *,   1.9893,   1.9921,   1.9968,   1.9987,   2.0034
 c     *,   2.0085,   2.0196,   2.0344,   2.0358,   2.0385
 c     *,   2.0545,   2.0636,   2.0654,   2.0704,   2.0816
 c     *,   2.0839,   2.0933,   2.1106,   2.1151,   2.1257
 c     *,   2.1279,   2.1388,   2.1545,   2.1606,   2.1761
 c     *,   2.1770,   2.1804,   2.1817,   2.1920,   2.2022
 c     *,   2.2060,   2.2187,   2.2480,   2.2646,   2.2934
 c     *,   2.3056,   2.3339,   2.3499,   2.3535,   2.3538
 c     *,   2.3737,   2.3894,   2.4050,   2.4128,   2.4283
 c     *,   2.4398,   2.4513,   2.4703,   2.4892,   2.5080
 c     *,   2.5266,   2.5451,   2.5561,   2.5634,   2.5816
 c     *,   2.5997,   2.6069,   2.6177,   2.6355,   2.6568
 c     *,   2.6708,   2.6987,   2.7057,   2.7195,   2.7401
 c     *,   2.7469,   2.7606,   2.7741,   2.7977,   2.8010
 c     *,   2.8077,   2.8343,   2.8409,   2.8672,   2.8769
 c     *,   2.8997,   2.9093,   2.9350,   2.9414,   2.9636
 c     *,   2.9731,   3.0013,   3.0045,   3.0107,   3.0355
 c     *,   3.0570,   3.0662,   3.0876,   3.0967,   3.1268
 c     *,   3.1328,   3.1566,   3.1862,   3.2067,   3.2155
 c     *,   3.2445,   3.2733,   3.3018,   3.3216,   3.3329
 c     *,   3.3609,   3.3887,   3.4163,   3.4190,   3.4327
 c     *,   3.4436,   3.4707,   3.4869,   3.4976,   3.5244
 c     *,   3.5509,   3.6033,   3.6550,   3.7059,   3.7462
 c     *,   3.9247,   3.9887,   3.9981,   4.4001,   4.4341
 c     *,   4.8193,   4.8387,   5.2120,   5.2246,   5.3889
 c     *,   5.5535,   5.5603,   5.7265,   5.8880,   5.8912
 c     *,   5.9671,   6.1984,   6.2421,   6.5084,   6.6369
 c     *,   6.9138,   7.5617,   8.1584,   8.7144,   8.9794
 c     *,   9.2369,   9.7314,   9.9412,  10.2020,  10.6517
 c     *,  11.4987,  13.7295,  15.0339,  16.6334,  16.8018
 c     *,  17.8835,  18.1439,  21.2400,  22.9386,  24.1342
 c     *,  25.2733
 c     */
 c      data pipbmx/
 c     *     .4424,    .5585,    .6180,    .7485,    .7092
 c     *,    .8058,    .7777,    .9113,    .8974,    .9934
 c     *,   1.0896,   1.2653,   1.3078,   1.5076,   1.5472
 c     *,   1.5656,   1.7019,   1.7504,   1.7934,   1.7626
 c     *,   1.8168,   1.9706,   2.0163,   2.1140,   2.1749
 c     *,   2.0576,   2.1851,   2.1148,   2.1924,   2.1851
 c     *,   2.3104,   2.3480,   2.2883,   2.1924,   2.3602
 c     *,   2.4570,   2.3262,   2.2708,   2.4722,   2.3296
 c     *,   2.5124,   2.3194,   2.4476,   2.4033,   2.5497
 c     *,   2.5313,   2.5514,   2.5181,   2.5125,   2.5156
 c     *,   2.5105,   2.4398,   2.4201,   2.5074,   2.4978
 c     *,   2.4463,   2.4863,   2.4856,   2.4469,   2.5231
 c     *,   2.3534,   2.5357,   2.4482,   2.3796,   2.3823
 c     *,   2.3582,   2.3650,   2.3870,   2.1705,   2.1185
 c     *,   2.2341,   2.2057,   2.2284,   2.1178,   2.1705
 c     *,   2.0622,   2.0490,   2.0122,   1.9041,   1.8873
 c     *,   1.9091,   1.7897,   1.7499,   1.7535,   1.8797
 c     *,   1.8455,   1.6087,   1.6381,   1.6361,   1.6737
 c     *,   1.5329,   1.5327,   1.4879,   1.4351,   1.4461
 c     *,   1.3854,   1.3517,   1.3334,   1.3695,   1.2741
 c     *,   1.2989,   1.2183,   1.2361,   1.2074,   1.1456
 c     *,   1.1997,   1.1439,   1.1381,   1.0786,   1.1185
 c     *,   1.0693,   1.1339,    .9958,   1.0077,   1.0491
 c     *,    .9997,    .9611,    .9452,    .9141,    .9190
 c     *,    .8806,    .8885,    .9200,    .8618,    .9680
 c     *,    .8702,    .8643,    .8263,    .8234,    .7150
 c     *,    .8412,    .8101,    .8704,    .7767,    .8766
 c     *,    .7722,    .7694,    .7956,    .7590,    .7128
 c     *,    .7159,    .6894,    .6894,    .7139,    .6898
 c     *,    .7412,    .6784,    .7181,    .6947,    .7695
 c     *,    .6704,    .7110,    .6910,    .6699,    .7199
 c     *,    .6935,    .6759,    .6843,    .6933,    .7190
 c     *,    .7159,    .6865,    .6958,    .7017,    .6794
 c     *,    .7081,    .7004,    .7128,    .7569,    .7269
 c     *,    .7356,    .7230,    .7332,    .7548,    .7464
 c     *,    .7695,    .7703,    .7707,    .8288,    .8117
 c     *,    .8170,    .8054,    .7990,    .7878,    .8556
 c     *,    .8423,    .8917,    .8664,    .8253,    .8649
 c     *,    .8627,    .8813,    .8461,    .8938,    .9288
 c     *,    .8759,    .8345,    .8972,    .8625,    .8818
 c     *,    .8983,    .8953,    .9080,    .9277,    .9117
 c     *,    .9080,    .8579,    .8974,    .8970,    .9184
 c     *,    .9267,    .9184,    .9092,    .8704,    .9146
 c     *,    .9184,    .9591,    .8649,    .9296,    .9071
 c     *,    .9449,    .9407,    .9536,    .9721,    .9608
 c     *,    .9837,    .9233,    .9542,   1.0004,    .9853
 c     *,    .9934,   1.0005,   1.0061,    .9358,   1.0180
 c     *,   1.0280,   1.0573,    .9982,   1.0357,   1.0555
 c     *,    .9877,   1.0682,   1.0694,   1.0817,   1.0720
 c     *,   1.1046,   1.1041,   1.0596,   1.1013,   1.1217
 c     *,   1.1217,   1.1291,   1.1113,   1.1270,   1.1354
 c     *,   1.0799,   1.1199,   1.1439,   1.1459,   1.1213
 c     *,   1.0907,   1.1368,   1.1488,   1.1008,   1.1500
 c     *,   1.1156,   1.1502,   1.0947,   1.1410,   1.1427
 c     *,   1.1480,   1.0729,   1.1634,   1.1319,   1.1480
 c     *,   1.1295,   1.1213,   1.0783,   1.1044,   1.1027
 c     *,   1.1217,   1.0998,   1.0418,   1.1031,   1.0712
 c     *,   1.0720,   1.0600,   1.0540,   1.1180,   1.0490
 c     *,   1.0464,   1.0380,   1.0124,   1.0202,    .9821
 c     *,   1.0013,    .9945,    .9997,   1.0187,    .9805
 c     *,    .9873,    .9796,    .9608,   1.0045,    .9756
 c     *,    .9695,    .9674,    .9646,    .9662,    .9659
 c     *,    .9654,    .9491,    .9643,    .9566,    .9654
 c     *,    .9636,    .9619,    .9676,    .9774,    .9747
 c     *,    .9916,    .9843,    .9950,    .9932,    .9641
 c     *,    .9895,    .9861,    .9871,    .9905,    .9916
 c     *,    .9657,    .9872,    .9770,    .9796,    .9751
 c     *,    .9637,    .9689,    .9633,    .9635,    .9601
 c     *,    .9580,    .9707,    .9581,    .9540,    .9505
 c     *,    .9509,    .9422,    .9494,    .9657,    .9478
 c     *,    .9513,    .9541,    .9472,    .9439,    .9640
 c     *,    .9452,    .9397,    .9440,    .9390,    .9416
 c     *,    .9424,    .9394,    .9334,    .9366,    .9348
 c     *,    .9338,    .9414,    .9312,    .9115,    .9286
 c     *,    .9373,    .9287,    .9190,    .9248,    .9232
 c     *,    .9339,    .9214,    .9201,    .9183,    .9181
 c     *,    .9170,    .9150,    .9138,    .9071,    .9117
 c     *,    .9106,    .9092,    .9084,    .9253,    .9053
 c     *,    .9061,    .9052,    .9132,    .9032,    .9029
 c     *,    .9008,    .8973,    .8954,    .8928,    .9115
 c     *,    .9021,    .8938,    .8997,    .8907,    .8921
 c     *,    .8834,    .8831,    .8795,    .8774,    .8755
 c     *,    .8745,    .8682,    .8849,    .8649,    .8705
 c     *,    .8616,    .8684,    .8691,    .8634,    .8687
 c     *,    .8636,    .8616,    .8581,    .8571,    .8575
 c     *,    .8582,    .8551,    .8575,    .8582,    .8618
 c     *,    .8597,    .8619,    .8634,    .8693,    .8645
 c     *,    .8682,    .8538,    .8759,    .8818,    .8831
 c     *,    .8853/
 c      data (pimecm(i),i=1,400)/
 c     *    1.1046,   1.1133,   1.1394,   1.1425,   1.1495
 c     *,   1.1579,   1.1585,   1.1592,   1.1598,   1.1677
 c     *,   1.1691,   1.1731,   1.1777,   1.1784,   1.1798
 c     *,   1.1831,   1.1858,   1.1879,   1.1892,   1.1906
 c     *,   1.1940,   1.1967,   1.1994,   1.2008,   1.2015
 c     *,   1.2028,   1.2069,   1.2076,   1.2090,   1.2097
 c     *,   1.2111,   1.2124,   1.2131,   1.2159,   1.2166
 c     *,   1.2179,   1.2186,   1.2200,   1.2214,   1.2234
 c     *,   1.2269,   1.2283,   1.2296,   1.2303,   1.2317
 c     *,   1.2324,   1.2352,   1.2358,   1.2365,   1.2372
 c     *,   1.2407,   1.2421,   1.2441,   1.2462,   1.2476
 c     *,   1.2510,   1.2517,   1.2524,   1.2545,   1.2551
 c     *,   1.2579,   1.2586,   1.2593,   1.2607,   1.2613
 c     *,   1.2627,   1.2634,   1.2641,   1.2662,   1.2668
 c     *,   1.2682,   1.2696,   1.2710,   1.2716,   1.2730
 c     *,   1.2758,   1.2778,   1.2806,   1.2819,   1.2826
 c     *,   1.2833,   1.2847,   1.2888,   1.2915,   1.2922
 c     *,   1.2929,   1.2963,   1.3004,   1.3024,   1.3038
 c     *,   1.3058,   1.3065,   1.3072,   1.3112,   1.3119
 c     *,   1.3126,   1.3146,   1.3180,   1.3187,   1.3200
 c     *,   1.3207,   1.3220,   1.3227,   1.3254,   1.3261
 c     *,   1.3301,   1.3308,   1.3321,   1.3335,   1.3362
 c     *,   1.3375,   1.3388,   1.3415,   1.3422,   1.3455
 c     *,   1.3482,   1.3515,   1.3522,   1.3555,   1.3562
 c     *,   1.3615,   1.3628,   1.3641,   1.3655,   1.3681
 c     *,   1.3694,   1.3760,   1.3773,   1.3786,   1.3898
 c     *,   1.3917,   1.3963,   1.3995,   1.4002,   1.4015
 c     *,   1.4047,   1.4099,   1.4112,   1.4163,   1.4176
 c     *,   1.4183,   1.4196,   1.4208,   1.4228,   1.4234
 c     *,   1.4272,   1.4292,   1.4298,   1.4304,   1.4393
 c     *,   1.4425,   1.4470,   1.4476,   1.4489,   1.4507
 c     *,   1.4539,   1.4589,   1.4596,   1.4608,   1.4652
 c     *,   1.4658,   1.4683,   1.4727,   1.4739,   1.4758
 c     *,   1.4764,   1.4777,   1.4795,   1.4808,   1.4833
 c     *,   1.4851,   1.4857,   1.4876,   1.4882,   1.4894
 c     *,   1.4901,   1.4919,   1.4925,   1.4938,   1.4968
 c     *,   1.4974,   1.4981,   1.4993,   1.5030,   1.5054
 c     *,   1.5060,   1.5066,   1.5072,   1.5097,   1.5103
 c     *,   1.5109,   1.5146,   1.5152,   1.5182,   1.5206
 c     *,   1.5212,   1.5224,   1.5231,   1.5249,   1.5309
 c     *,   1.5327,   1.5345,   1.5357,   1.5387,   1.5405
 c     *,   1.5411,   1.5429,   1.5435,   1.5441,   1.5465
 c     *,   1.5489,   1.5507,   1.5513,   1.5519,   1.5548
 c     *,   1.5590,   1.5614,   1.5637,   1.5655,   1.5673
 c     *,   1.5702,   1.5708,   1.5761,   1.5802,   1.5825
 c     *,   1.5831,   1.5843,   1.5849,   1.5866,   1.5890
 c     *,   1.5936,   1.5948,   1.5977,   1.5994,   1.6017
 c     *,   1.6052,   1.6058,   1.6087,   1.6092,   1.6138
 c     *,   1.6173,   1.6178,   1.6190,   1.6236,   1.6258
 c     *,   1.6276,   1.6287,   1.6298,   1.6315,   1.6327
 c     *,   1.6367,   1.6372,   1.6378,   1.6400,   1.6406
 c     *,   1.6423,   1.6480,   1.6513,   1.6519,   1.6525
 c     *,   1.6547,   1.6553,   1.6614,   1.6631,   1.6637
 c     *,   1.6648,   1.6653,   1.6659,   1.6681,   1.6692
 c     *,   1.6703,   1.6720,   1.6726,   1.6731,   1.6742
 c     *,   1.6792,   1.6798,   1.6803,   1.6820,   1.6825
 c     *,   1.6853,   1.6858,   1.6864,   1.6875,   1.6935
 c     *,   1.6957,   1.6963,   1.6979,   1.6990,   1.7067
 c     *,   1.7072,   1.7121,   1.7127,   1.7132,   1.7186
 c     *,   1.7202,   1.7208,   1.7213,   1.7224,   1.7235
 c     *,   1.7240,   1.7267,   1.7278,   1.7294,   1.7342
 c     *,   1.7348,   1.7391,   1.7407,   1.7428,   1.7476
 c     *,   1.7503,   1.7535,   1.7545,   1.7609,   1.7614
 c     *,   1.7667,   1.7714,   1.7735,   1.7741,   1.7767
 c     *,   1.7772,   1.7798,   1.7824,   1.7835,   1.7871
 c     *,   1.7929,   1.7960,   1.8001,   1.8032,   1.8058
 c     *,   1.8130,   1.8182,   1.8212,   1.8218,   1.8238
 c     *,   1.8253,   1.8258,   1.8289,   1.8320,   1.8370
 c     *,   1.8386,   1.8436,   1.8446,   1.8507,   1.8512
 c     *,   1.8557,   1.8588,   1.8638,   1.8678,   1.8713
 c     *,   1.8762,   1.8777,   1.8792,   1.8807,   1.8827
 c     *,   1.8857,   1.8881,   1.8886,   1.8936,   1.8975
 c     *,   1.9010,   1.9063,   1.9132,   1.9186,   1.9205
 c     *,   1.9249,   1.9254,   1.9273,   1.9283,   1.9317
 c     *,   1.9326,   1.9374,   1.9423,   1.9451,   1.9495
 c     *,   1.9519,   1.9538,   1.9614,   1.9709,   1.9733
 c     *,   1.9747,   1.9799,   1.9884,   1.9903,   1.9987
 c     *,   2.0034,   2.0117,   2.0164,   2.0182,   2.0196
 c     *,   2.0330,   2.0335,   2.0358,   2.0540,   2.0581
 c     *,   2.0636,   2.0654,   2.0798,   2.0812,   2.0816
 c     */
 c      data (pimecm(i),i=401,578)/
 c     *    2.0933,   2.1013,   2.1075,   2.1221,   2.1261
 c     *,   2.1305,   2.1375,   2.1445,   2.1588,   2.1658
 c     *,   2.1740,   2.1804,   2.1877,   2.2026,   2.2119
 c     *,   2.2166,   2.2187,   2.2305,   2.2510,   2.2626
 c     *,   2.2688,   2.2717,   2.3040,   2.3121,   2.3315
 c     *,   2.3483,   2.3538,   2.3737,   2.3744,   2.3925
 c     *,   2.4050,   2.4105,   2.4128,   2.4267,   2.4302
 c     *,   2.4513,   2.4627,   2.4892,   2.5069,   2.5266
 c     *,   2.5561,   2.5634,   2.5802,   2.5925,   2.5997
 c     *,   2.6284,   2.6355,   2.6557,   2.6708,   2.6987
 c     *,   2.7057,   2.7401,   2.7673,   2.7741,   2.7966
 c     *,   2.8077,   2.8343,   2.8409,   2.8672,   2.8769
 c     *,   2.8997,   2.9093,   2.9318,   2.9341,   2.9414
 c     *,   2.9636,   2.9731,   3.0045,   3.0262,   3.0355
 c     *,   3.0570,   3.0656,   3.0662,   3.0876,   3.0967
 c     *,   3.1268,   3.1566,   3.1774,   3.1862,   3.2067
 c     *,   3.2155,   3.2445,   3.2733,   3.3018,   3.3216
 c     *,   3.3329,   3.3609,   3.3887,   3.4163,   3.4190
 c     *,   3.4436,   3.4675,   3.4707,   3.4869,   3.4976
 c     *,   3.5509,   3.6033,   3.6550,   3.6575,   3.7059
 c     *,   3.7312,   3.7462,   3.8402,   3.8935,   3.9887
 c     *,   4.0004,   4.3873,   4.4341,   4.6811,   4.7210
 c     *,   4.7408,   4.8387,   4.8406,   5.0288,   5.0844
 c     *,   5.2120,   5.2353,   5.3889,   5.4254,   5.5603
 c     *,   5.6123,   5.7265,   5.7787,   5.8880,   5.9451
 c     *,   5.9671,   5.9953,   6.0792,   6.1984,   6.2241
 c     *,   6.3479,   6.5070,   6.6369,   6.8140,   6.9138
 c     *,   7.0734,   7.2450,   7.3962,   7.5617,   7.7092
 c     *,   7.9841,   7.9865,   8.1584,   8.1814,   8.3628
 c     *,   8.4410,   8.6582,   8.6593,   8.7144,   8.9794
 c     *,   9.2369,   9.2845,   9.4874,   9.7314,   9.7775
 c     *,   9.9694,  10.2020,  10.3580,  10.4293,  10.5988
 c     *,  10.6517,  10.8697,  11.0833,  11.4987,  13.7295
 c     *,  15.0339,  16.6334,  16.8018,  17.8835,  18.1439
 c     *,  19.3933,  19.6336,  21.2400,  22.9386,  24.1342
 c     *,  25.2733,  26.0050,  26.3632
 c     */
 c      data (pimbmx(i),i=1,400)/
 c     *     .5314,    .5314,    .6580,    .7092,    .7506
 c     *,    .8215,    .8579,    .8406,    .8349,    .9441
 c     *,    .9657,   1.0363,   1.0686,   1.1085,   1.1185
 c     *,   1.1899,   1.1658,   1.2218,   1.2322,   1.2653
 c     *,   1.3273,   1.2374,   1.3611,   1.3716,   1.3267
 c     *,   1.3986,   1.4150,   1.3399,   1.4701,   1.4516
 c     *,   1.4723,   1.4955,   1.4494,   1.4161,   1.4558
 c     *,   1.4625,   1.5006,   1.5128,   1.5160,   1.4584
 c     *,   1.4991,   1.4791,   1.4217,   1.4799,   1.4799
 c     *,   1.4588,   1.4844,   1.4172,   1.4273,   1.4439
 c     *,   1.4166,   1.4127,   1.3739,   1.3624,   1.3971
 c     *,   1.3297,   1.3231,   1.3285,   1.2878,   1.2952
 c     *,   1.2679,   1.2641,   1.2976,   1.2386,   1.2463
 c     *,   1.3025,   1.2489,   1.2239,   1.1902,   1.2114
 c     *,   1.1955,   1.2083,   1.1658,   1.1686,   1.1603
 c     *,   1.1424,   1.1185,   1.1160,   1.1013,   1.1070
 c     *,   1.1041,   1.0823,   1.0645,   1.0779,   1.0295
 c     *,   1.0311,    .9950,   1.0155,    .9950,   1.0029
 c     *,    .9767,   1.0249,    .9853,    .9657,    .9906
 c     *,    .9608,    .9591,    .9805,    .9458,    .9575
 c     *,    .9422,    .9407,    .9558,    .9473,    .9271
 c     *,    .9132,    .9236,    .9219,    .9219,    .9575
 c     *,    .9097,    .9224,    .9094,    .9202,    .8903
 c     *,    .9224,    .9160,    .9122,    .9094,    .8956
 c     *,    .9167,    .9271,    .9398,    .9184,    .9277
 c     *,    .9310,    .9184,    .9354,    .9219,    .9300
 c     *,    .9303,    .9508,    .9451,    .9680,    .9399
 c     *,    .9420,    .9387,    .9667,    .9248,    .9915
 c     *,    .9698,    .9626,   1.0187,    .9659,    .9441
 c     *,    .9766,    .9694,   1.0061,    .9789,    .9757
 c     *,    .9803,   1.0495,   1.0218,   1.0026,   1.0706
 c     *,   1.0277,   1.0302,   1.0690,   1.0552,   1.0406
 c     *,   1.0801,   1.0522,   1.1284,   1.1226,   1.1270
 c     *,   1.0934,   1.0968,   1.1381,   1.0915,   1.1543
 c     *,   1.1775,   1.2101,   1.1247,   1.1506,   1.1781
 c     *,   1.1546,   1.1968,   1.1604,   1.1697,   1.2072
 c     *,   1.2058,   1.2489,   1.1912,   1.2141,   1.1928
 c     *,   1.1986,   1.1978,   1.2256,   1.2565,   1.2029
 c     *,   1.2127,   1.1997,   1.2035,   1.1791,   1.1853
 c     *,   1.1914,   1.1927,   1.2616,   1.2303,   1.1848
 c     *,   1.1507,   1.2399,   1.1510,   1.1898,   1.1499
 c     *,   1.1337,   1.1433,   1.1089,   1.1517,   1.2008
 c     *,   1.1598,   1.1030,   1.0798,   1.1165,   1.1213
 c     *,   1.1190,   1.0861,   1.0776,   1.0653,   1.0730
 c     *,   1.1106,   1.0908,   1.0945,   1.0650,   1.0882
 c     *,   1.0720,   1.0540,   1.0925,   1.0718,   1.1377
 c     *,   1.0867,   1.0802,   1.1090,   1.0936,   1.0915
 c     *,   1.1316,   1.1575,   1.1202,   1.1275,   1.0908
 c     *,   1.1142,   1.1507,   1.1890,   1.1863,   1.2114
 c     *,   1.1142,   1.2001,   1.1312,   1.1930,   1.2563
 c     *,   1.2489,   1.2919,   1.2652,   1.2425,   1.2149
 c     *,   1.2386,   1.2616,   1.3207,   1.3299,   1.2816
 c     *,   1.3296,   1.3562,   1.2957,   1.3831,   1.3028
 c     *,   1.3806,   1.3704,   1.3202,   1.4082,   1.3482
 c     *,   1.3730,   1.3351,   1.3886,   1.3348,   1.3957
 c     *,   1.3929,   1.3790,   1.3791,   1.3946,   1.3641
 c     *,   1.3877,   1.3566,   1.3564,   1.3762,   1.3253
 c     *,   1.3739,   1.3252,   1.3234,   1.3423,   1.3159
 c     *,   1.2643,   1.2700,   1.2584,   1.2466,   1.2841
 c     *,   1.2191,   1.2052,   1.2794,   1.2835,   1.2374
 c     *,   1.2061,   1.1936,   1.1947,   1.2476,   1.1481
 c     *,   1.1496,   1.2187,   1.1427,   1.2101,   1.1192
 c     *,   1.1445,   1.1371,   1.1613,   1.0994,   1.1020
 c     *,   1.1233,   1.1041,   1.0896,   1.0890,   1.0880
 c     *,   1.0802,   1.0946,   1.0940,   1.0670,   1.0761
 c     *,   1.0817,   1.0702,   1.0723,   1.0779,   1.0789
 c     *,   1.0733,   1.0771,   1.0767,   1.0710,   1.0705
 c     *,   1.0720,   1.0729,   1.0472,   1.0959,   1.0788
 c     *,   1.0739,   1.0645,   1.0681,   1.0749,   1.0789
 c     *,   1.0799,   1.0795,   1.0763,   1.0588,   1.0805
 c     *,   1.0724,   1.0660,   1.0663,   1.0801,   1.0988
 c     *,   1.0797,   1.0816,   1.0687,   1.0690,   1.0724
 c     *,   1.0721,   1.0777,   1.0650,   1.0715,   1.0617
 c     *,   1.0690,   1.0678,   1.0714,   1.0323,   1.0672
 c     *,   1.0630,   1.0564,   1.0600,   1.0636,   1.0499
 c     *,   1.0479,   1.0481,   1.0463,   1.0510,   1.0450
 c     *,   1.0171,   1.0517,   1.0498,   1.0510,   1.0412
 c     *,   1.0374,   1.0501,   1.0388,   1.0449,   1.0449
 c     *,   1.0311,   1.0510,   1.0505,   1.0449,   1.0567
 c     *,   1.0572,   1.0495,   1.0632,   1.0461,   1.0403
 c     */
 c      data (pimbmx(i),i=401,578)/
 c     *    1.0629,   1.0642,   1.0525,   1.0709,   1.0585
 c     *,   1.0612,   1.0717,   1.0731,   1.0660,   1.0761
 c     *,   1.0696,   1.0767,   1.0755,   1.0720,   1.0615
 c     *,   1.0665,   1.0660,   1.0714,   1.0670,   1.0650
 c     *,   1.0627,   1.0621,   1.0499,   1.0499,   1.0457
 c     *,   1.0405,   1.0372,   1.0412,   1.0331,   1.0299
 c     *,   1.0299,   1.0241,   1.0321,   1.0306,   1.0232
 c     *,   1.0252,   1.0171,   1.0211,   1.0226,   1.0178
 c     *,   1.0083,   1.0157,   1.0138,    .9918,   1.0140
 c     *,   1.0028,   1.0120,   1.0077,   1.0088,    .9979
 c     *,   1.0052,   1.0025,    .9876,    .9987,    .9990
 c     *,    .9944,    .9845,    .9918,    .9777,    .9892
 c     *,    .9754,    .9856,    .9812,    .9816,    .9831
 c     *,    .9674,    .9800,    .9781,    .9679,    .9756
 c     *,    .9805,    .9759,    .9739,    .9672,    .9718
 c     *,    .9701,    .9682,    .9632,    .9664,    .9538
 c     *,    .9647,    .9628,    .9613,    .9596,    .9449
 c     *,    .9585,    .9572,    .9553,    .9540,    .9624
 c     *,    .9527,    .9468,    .9513,    .9525,    .9503
 c     *,    .9479,    .9459,    .9438,    .9424,    .9422
 c     *,    .9407,    .9457,    .9399,    .9385,    .9368
 c     *,    .9409,    .9248,    .9209,    .9034,    .8974
 c     *,    .9150,    .9089,    .9145,    .9146,    .9080
 c     *,    .9044,    .9082,    .9069,    .9062,    .8938
 c     *,    .9034,    .9045,    .9011,    .8921,    .8979
 c     *,    .8925,    .8982,    .8975,    .8976,    .8947
 c     *,    .8962,    .8932,    .8916,    .8913,    .8888
 c     *,    .8889,    .8892,    .8880,    .8881,    .8842
 c     *,    .8815,    .8826,    .8832,    .8813,    .8826
 c     *,    .8791,    .8835,    .8831,    .8833,    .8782
 c     *,    .8789,    .8813,    .8797,    .8780,    .8744
 c     *,    .8782,    .8842,    .8777,    .8777,    .8789
 c     *,    .8815,    .8800,    .8839,    .8740,    .8732
 c     *,    .8751,    .8784,    .8757,    .8779,    .8630
 c     *,    .8802,    .8775,    .8851,    .8874,    .8903
 c     *,    .8935,    .8965,    .8965
 c     */
 c      data kmpecm/
 c     *    1.4691,   1.4720,   1.4750,   1.4780,   1.4811
 c     *,   1.4837,   1.4843,   1.4860,   1.4876,   1.4910
 c     *,   1.4944,   1.4979,   1.5014,   1.5032,   1.5050
 c     *,   1.5087,   1.5091,   1.5124,   1.5132,   1.5162
 c     *,   1.5170,   1.5200,   1.5220,   1.5239,   1.5278
 c     *,   1.5318,   1.5354,   1.5358,   1.5362,   1.5378
 c     *,   1.5399,   1.5440,   1.5523,   1.5607,   1.5654
 c     *,   1.5688,   1.5775,   1.5784,   1.5863,   1.5916
 c     *,   1.5947,   1.6023,   1.6050,   1.6055,   1.6086
 c     *,   1.6145,   1.6159,   1.6172,   1.6191,   1.6236
 c     *,   1.6319,   1.6328,   1.6332,   1.6420,   1.6461
 c     *,   1.6466,   1.6522,   1.6563,   1.6582,   1.6614
 c     *,   1.6642,   1.6694,   1.6712,   1.6717,   1.6768
 c     *,   1.6806,   1.6811,   1.6839,   1.6843,   1.6867
 c     *,   1.6885,   1.6961,   1.6965,   1.7003,   1.7022
 c     *,   1.7083,   1.7087,   1.7172,   1.7177,   1.7181
 c     *,   1.7229,   1.7243,   1.7276,   1.7342,   1.7374
 c     *,   1.7436,   1.7459,   1.7483,   1.7539,   1.7610
 c     *,   1.7629,   1.7633,   1.7718,   1.7770,   1.7789
 c     *,   1.7793,   1.7817,   1.7840,   1.7873,   1.7892
 c     *,   1.8028,   1.8037,   1.8075,   1.8135,   1.8140
 c     *,   1.8219,   1.8247,   1.8261,   1.8308,   1.8369
 c     *,   1.8401,   1.8406,   1.8410,   1.8475,   1.8480
 c     *,   1.8489,   1.8540,   1.8559,   1.8605,   1.8647
 c     *,   1.8721,   1.8744,   1.8767,   1.8785,   1.8836
 c     *,   1.8951,   1.8955,   1.8983,   1.9001,   1.9029
 c     *,   1.9065,   1.9184,   1.9202,   1.9243,   1.9348
 c     *,   1.9393,   1.9411,   1.9434,   1.9483,   1.9547
 c     *,   1.9628,   1.9637,   1.9659,   1.9699,   1.9798
 c     *,   1.9838,   1.9923,   1.9958,   2.0047,   2.0127
 c     *,   2.0149,   2.0162,   2.0255,   2.0277,   2.0308
 c     *,   2.0417,   2.0430,   2.0487,   2.0579,   2.0653
 c     *,   2.0679,   2.0713,   2.0813,   2.0882,   2.0925
 c     *,   2.1028,   2.1105,   2.1126,   2.1211,   2.1232
 c     *,   2.1258,   2.1351,   2.1444,   2.1507,   2.1528
 c     *,   2.1654,   2.1675,   2.1687,   2.1737,   2.1837
 c     *,   2.1862,   2.1937,   2.2044,   2.2131,   2.2143
 c     *,   2.2274,   2.2356,   2.2478,   2.2546,   2.2659
 c     *,   2.2756,   2.2836,   2.2976,   2.2996,   2.3162
 c     *,   2.3166,   2.3296,   2.3335,   2.3363,   2.3535
 c     *,   2.3562,   2.3725,   2.3729,   2.3748,   2.3887
 c     *,   2.3918,   2.3933,   2.4006,   2.4109,   2.4174
 c     *,   2.4223,   2.4299,   2.4352,   2.4488,   2.4518
 c     *,   2.4675,   2.4705,   2.4861,   2.4887,   2.4936
 c     *,   2.5046,   2.5230,   2.5412,   2.5593,   2.5701
 c     *,   2.5773,   2.5952,   2.6023,   2.6059,   2.6130
 c     *,   2.6306,   2.6447,   2.6482,   2.6656,   2.6795
 c     *,   2.6829,   2.7002,   2.7173,   2.7848,   2.8540
 c     *,   2.9216,   2.9407,   2.9878,   3.0096,   3.0250
 c     *,   3.0526,   3.1783,   3.2191,   3.2993,   3.3887
 c     *,   3.5239,   3.6924,   3.7800,   3.9967,   4.0200
 c     *,   4.1348,   4.4617,   4.4826,   4.7663,   4.8636
 c     *,   4.9779,   5.1080,   5.1263,   5.2349,   5.3237
 c     *,   5.4110,   5.5479,   5.5816,   5.8281,   5.9080
 c     *,   6.0801,   6.2173,   6.3664,   6.6545,   6.9306
 c     *,   7.2610,   7.5770,   7.7605,   7.8207,   7.9985
 c     *,   8.1725,   8.2297,   8.4546,   8.7275,   8.9922
 c     *,   9.2493,   9.4994,   9.7431,   9.9809,  10.2131
 c     *,  11.5086,  13.7378,  15.0415,  16.8085,  17.8898
 c     *,  19.3991,  21.2453,  22.9435,  24.1389
 c     */
 c      data kmpbmx/
 c     *    1.9033,   1.7662,   1.7298,   1.7544,   1.5461
 c     *,   1.6991,   1.6205,   1.5898,   1.5817,   1.5023
 c     *,   1.5554,   1.5086,   1.5065,   1.4948,   1.4799
 c     *,   1.4927,   1.4854,   1.6136,   1.6017,   1.7312
 c     *,   1.4884,   1.7075,   1.5574,   1.5260,   1.5002
 c     *,   1.4571,   1.3991,   1.3219,   1.3434,   1.3635
 c     *,   1.3315,   1.2966,   1.2213,   1.1604,   1.1930
 c     *,   1.1382,   1.1354,   1.1583,   1.1185,   1.1156
 c     *,   1.0645,   1.0247,   1.0969,   1.0779,   1.1083
 c     *,   1.0749,   1.0155,   1.0464,   1.0754,   1.0202
 c     *,   1.0470,   1.0615,   1.0540,   1.0357,   1.0605
 c     *,   1.0428,   1.0265,   1.0187,   1.0611,   1.0615
 c     *,   1.0295,   1.0651,   1.0823,   1.0077,   1.0764
 c     *,   1.1298,   1.1013,   1.1080,   1.1070,   1.1368
 c     *,   1.1382,   1.1410,   1.1418,   1.1312,   1.1466
 c     *,   1.1378,   1.1241,   1.0896,   1.1264,   1.0720
 c     *,   1.1368,   1.1013,   1.1442,   1.1466,   1.1612
 c     *,   1.1726,   1.1755,   1.1576,   1.1726,   1.1848
 c     *,   1.1617,   1.2029,   1.1686,   1.2274,   1.2252
 c     *,   1.1726,   1.2256,   1.2244,   1.2008,   1.2332
 c     *,   1.2828,   1.2548,   1.2412,   1.2887,   1.2853
 c     *,   1.2653,   1.2540,   1.2527,   1.2449,   1.2118
 c     *,   1.1781,   1.1902,   1.1767,   1.1859,   1.1410
 c     *,   1.1594,   1.1800,   1.1185,   1.1256,   1.1354
 c     *,   1.1131,   1.1562,   1.1143,   1.1382,   1.0841
 c     *,   1.0587,   1.0632,   1.0311,   1.0617,   1.0150
 c     *,   1.0309,   1.0174,   1.0110,   1.0171,    .9961
 c     *,    .9719,    .9938,    .9869,    .9953,    .9966
 c     *,   1.0189,    .9977,    .9918,    .9948,   1.0034
 c     *,    .9737,   1.0066,   1.0137,   1.0041,   1.0171
 c     *,   1.0303,   1.0223,   1.0322,    .9751,   1.0331
 c     *,   1.0118,   1.0403,   1.0399,   1.0429,   1.0171
 c     *,   1.0279,   1.0467,   1.0414,   1.0217,   1.0432
 c     *,   1.0379,   1.0280,   1.0351,   1.0171,   1.0278
 c     *,   1.0314,   1.0213,   1.0133,   1.0240,   1.0080
 c     *,   1.0018,    .9964,    .9984,    .9889,    .9910
 c     *,    .9903,    .9801,    .9852,    .9853,    .9847
 c     *,    .9800,    .9832,    .9770,    .9749,    .9754
 c     *,    .9723,    .9741,    .9744,    .9738,    .9751
 c     *,    .9780,    .9738,    .9738,    .9684,    .9712
 c     *,    .9669,    .9680,    .9671,    .9576,    .9619
 c     *,    .9624,    .9585,    .9588,    .9586,    .9525
 c     *,    .9253,    .9518,    .9503,    .9491,    .9463
 c     *,    .9476,    .9420,    .9458,    .9434,    .9341
 c     *,    .9444,    .9412,    .9393,    .9395,    .9210
 c     *,    .9370,    .9358,    .8974,    .9400,    .9342
 c     *,    .9305,    .9141,    .9271,    .9267,    .9228
 c     *,    .9233,    .9219,    .9247,    .9296,    .9253
 c     *,    .9089,    .8992,    .8946,    .8992,    .9062
 c     *,    .9069,    .8874,    .9219,    .8746,    .8795
 c     *,    .8740,    .8704,    .8785,    .8795,    .8667
 c     *,    .8849,    .8510,    .8463,    .8612,    .8292
 c     *,    .8292,    .8387,    .8273,    .8273,    .8292
 c     *,    .8292,    .8349,    .8234,    .8349,    .8176
 c     *,    .8292,    .8180,    .8193,    .8154,    .8130
 c     *,    .8121,    .8145,    .8078,    .7959,    .8088
 c     *,    .8075,    .8064,    .8056,    .8086,    .8048
 c     *,    .8080,    .8068,    .8050,    .8042,    .8050
 c     *,    .8054,    .8064,    .8096,    .8095,    .8107
 c     *,    .8135,    .8234,    .8238,    .8263
 c     */
 c      data kmnecm/
 c     *    1.6212,   1.6781,   1.7053,   1.7863,   1.7882
 c     *,   1.8271,   1.9025,   2.0152,   2.1237,   2.2157
 c     *,   2.4252,   2.6054,   2.6160,   2.9441,   3.5279
 c     *,   3.6965,   4.0245,   4.4666,   4.8690,   5.1136
 c     *,   5.2406,   5.4169,   5.5877,   5.9144,   6.2241
 c     *,   6.9381,   7.5852,   8.1813,   8.7369,   9.2592
 c     *,   9.7536,  10.2241,  11.5209,  13.7525,  15.0575
 c     *,  16.8264,  17.9089,  19.4198,  21.2680,  22.9680
 c     *,  24.1646
 c     */
 c      data kmnbmx/
 c     *     .9132,    .9966,    .9772,    .9966,   1.1382
 c     *,   1.0794,    .9674,    .9167,    .8795,    .8500
 c     *,    .8482,    .8444,    .8330,    .8078,    .8349
 c     *,    .8368,    .7919,    .8146,    .8019,    .8156
 c     *,    .7999,    .7961,    .8038,    .8038,    .7949
 c     *,    .7868,    .7864,    .7910,    .7850,    .7870
 c     *,    .7906,    .7885,    .7945,    .7971,    .8011
 c     *,    .8003,    .8029,    .8082,    .8088,    .8156
 c     *,    .8189
 c     */
 c      data kppecm/
 c     *    1.4447,   1.4516,   1.4522,   1.4585,   1.4663
 c     *,   1.4750,   1.5036,   1.5050,   1.5091,   1.5239
 c     *,   1.5378,   1.5523,   1.5654,   1.5784,   1.5916
 c     *,   1.5929,   1.6014,   1.6037,   1.6050,   1.6150
 c     *,   1.6159,   1.6191,   1.6259,   1.6264,   1.6268
 c     *,   1.6328,   1.6378,   1.6461,   1.6517,   1.6587
 c     *,   1.6605,   1.6652,   1.6792,   1.6843,   1.6853
 c     *,   1.6928,   1.7073,   1.7101,   1.7210,   1.7294
 c     *,   1.7374,   1.7422,   1.7483,   1.7539,   1.7643
 c     *,   1.7662,   1.7704,   1.7789,   1.7793,   1.7864
 c     *,   1.7897,   1.8028,   1.8070,   1.8135,   1.8191
 c     *,   1.8327,   1.8355,   1.8373,   1.8517,   1.8587
 c     *,   1.8605,   1.8679,   1.8725,   1.8813,   1.8836
 c     *,   1.9038,   1.9070,   1.9289,   1.9298,   1.9321
 c     *,   1.9524,   1.9533,   1.9749,   1.9807,   1.9950
 c     *,   1.9972,   1.9994,   2.0074,   2.0193,   2.0435
 c     *,   2.0492,   2.0635,   2.0653,   2.0851,   2.1040
 c     *,   2.1066,   2.1083,   2.1279,   2.1296,   2.1490
 c     *,   2.1507,   2.1717,   2.1908,   2.1924,   2.2110
 c     *,   2.2131,   2.2213,   2.2319,   2.2335,   2.2538
 c     *,   2.2724,   2.2740,   2.2940,   2.3123,   2.3138
 c     *,   2.3375,   2.3531,   2.3725,   2.3903,   2.3918
 c     *,   2.4109,   2.4197,   2.4299,   2.4413,   2.4488
 c     *,   2.4675,   2.4861,   2.5046,   2.5230,   2.5266
 c     *,   2.5412,   2.5521,   2.5593,   2.5773,   2.5952
 c     *,   2.6130,   2.6306,   2.6482,   2.6656,   2.6829
 c     *,   2.7002,   2.7173,   3.0096,   3.0556,   3.1754
 c     *,   3.3887,   3.5239,   3.7800,   4.0200,   4.1348
 c     *,   4.4617,   4.6469,   4.7663,   4.8636,   4.9591
 c     *,   5.1263,   5.2349,   5.4110,   5.5816,   5.7308
 c     *,   5.9080,   6.0646,   6.2173,   6.9306,   7.5770
 c     *,   7.8207,   8.1725,   8.7275,   9.2493,   9.7431
 c     *,  10.2131,  11.5086,  13.7378,  15.0415,  16.6402
 c     *,  16.8085,  17.8898,  18.1501,  19.3991,  21.2453
 c     *,  22.9435,  24.1389
 c     */
 c      data kppbmx/
 c     *     .5412,    .6308,    .6024,    .6050,    .5971
 c     *,    .6037,    .6232,    .6103,    .6482,    .6601
 c     *,    .6386,    .6466,    .6438,    .6204,    .6482
 c     *,    .6358,    .6333,    .6445,    .6443,    .6346
 c     *,    .6410,    .6227,    .6283,    .6308,    .6403
 c     *,    .6290,    .6457,    .5984,    .6333,    .5955
 c     *,    .5944,    .6295,    .6346,    .6090,    .6433
 c     *,    .6383,    .6482,    .6425,    .6543,    .6588
 c     *,    .6652,    .6768,    .6730,    .6723,    .6815
 c     *,    .7040,    .6898,    .7014,    .7001,    .7181
 c     *,    .7130,    .7159,    .7067,    .7440,    .7345
 c     *,    .7365,    .7485,    .7382,    .7474,    .7574
 c     *,    .7588,    .7559,    .7590,    .7582,    .7668
 c     *,    .7592,    .7682,    .7661,    .7697,    .7548
 c     *,    .7661,    .7626,    .7626,    .7563,    .7590
 c     *,    .7578,    .7611,    .7557,    .7555,    .7506
 c     *,    .7498,    .7517,    .7508,    .7540,    .7464
 c     *,    .7538,    .7512,    .7527,    .7534,    .7527
 c     *,    .7565,    .7521,    .7529,    .7525,    .7444
 c     *,    .7517,    .7334,    .7485,    .7491,    .7510
 c     *,    .7466,    .7476,    .7478,    .7472,    .7485
 c     *,    .7378,    .7451,    .7468,    .7474,    .7476
 c     *,    .7459,    .7410,    .7461,    .7457,    .7414
 c     *,    .7464,    .7457,    .7444,    .7444,    .7444
 c     *,    .7442,    .7291,    .7421,    .7429,    .7421
 c     *,    .7397,    .7386,    .7373,    .7389,    .7384
 c     *,    .7384,    .7386,    .7378,    .7569,    .7878
 c     *,    .7548,    .7356,    .7526,    .7421,    .7715
 c     *,    .7568,    .7590,    .7777,    .7421,    .7632
 c     *,    .7464,    .7442,    .7548,    .7367,    .7736
 c     *,    .7378,    .7421,    .7455,    .7502,    .7510
 c     *,    .7653,    .7529,    .7580,    .7544,    .7614
 c     *,    .7605,    .7661,    .7703,    .7805,    .7883
 c     *,    .7850,    .7907,    .7611,    .7961,    .8023
 c     *,    .8068,    .8111
 c     */
 c      data kpnecm/
 c     *    1.6875,   1.7430,   1.7670,   1.7816,   1.7905
 c     *,   1.8144,   1.8383,   1.8615,   1.8749,   1.8846
 c     *,   1.9080,   1.9308,   1.9345,   1.9535,   1.9760
 c     *,   1.9974,   1.9983,   2.0205,   2.0460,   2.0647
 c     *,   2.0678,   2.0864,   2.1066,   2.1079,   2.1109
 c     *,   2.1292,   2.1322,   2.1504,   2.1533,   2.1743
 c     *,   2.1922,   2.1951,   2.2157,   2.2239,   2.2334
 c     *,   2.2362,   2.2565,   2.2739,   2.2767,   2.2967
 c     *,   2.3138,   2.3166,   2.3402,   2.3559,   2.3753
 c     *,   2.3919,   2.3946,   2.4138,   2.4328,   2.4517
 c     *,   2.4704,   2.4891,   2.5076,   2.5259,   2.5442
 c     *,   2.5551,   2.5623,   2.5803,   2.5982,   2.6160
 c     *,   2.6337,   2.6513,   2.6687,   2.6861,   2.7033
 c     *,   2.7204,   3.5279,   4.0245,   4.4666,   4.8690
 c     *,   5.2406,   5.4169,   5.5877,   5.9144,   6.2241
 c     *,   6.9381,   7.5852,   8.1813,   8.7369,   9.2592
 c     *,   9.7536,  10.2241,  11.5209,  13.7525,  15.0575
 c     *,  16.8264,  17.9089,  19.4198,  21.2680,  22.9680
 c     *,  24.1646
 c     */
 c      data kpnbmx/
 c     *     .7024,    .7324,    .7485,    .7527,    .7680
 c     *,    .7758,    .8100,    .8224,    .7611,    .8151
 c     *,    .8031,    .7915,    .7674,    .7842,    .7822
 c     *,    .7590,    .7791,    .7767,    .7758,    .7734
 c     *,    .7754,    .7709,    .7674,    .7713,    .7742
 c     *,    .7752,    .7748,    .7721,    .7680,    .7707
 c     *,    .7674,    .7713,    .7715,    .7695,    .7684
 c     *,    .7734,    .7682,    .7709,    .7672,    .7659
 c     *,    .7653,    .7653,    .7506,    .7626,    .7624
 c     *,    .7701,    .7588,    .7622,    .7592,    .7491
 c     *,    .7588,    .7574,    .7592,    .7582,    .7571
 c     *,    .7464,    .7559,    .7538,    .7529,    .7529
 c     *,    .7534,    .7538,    .7487,    .7487,    .7498
 c     *,    .7474,    .7464,    .7485,    .7464,    .7485
 c     *,    .7464,    .7565,    .7442,    .7485,    .7545
 c     *,    .7555,    .7542,    .7634,    .7653,    .7686
 c     *,    .7671,    .7723,    .7695,    .7785,    .7824
 c     *,    .7905,    .7927,    .7923,    .8052,    .8100
 c     *,    .8137
 c     */
 cc      data lpecm/
 cc     *    2.0577,   2.0583,   2.0593,   2.0595,   2.0609
 cc     *,   2.0617,   2.0629,   2.0642,   2.0643,   2.0647
 cc     *,   2.0666,   2.0671,   2.0683,   2.0692,   2.0709
 cc     *,   2.0720,   2.0783,   2.0818,   2.0935,   2.1022
 cc     *,   2.1058,   2.1117,   2.1326,   2.1559,   2.1731
 cc     *,   2.1811,   2.2079,   2.2360,   2.2505,   2.3107
 cc     *,   2.3733,   2.4534,   2.4695,   2.9278,   5.2476
 cc     */
 cc      data lpbmx/
 cc     *    2.5793,   2.3937,   1.9381,   2.3736,   2.0342
 cc     *,   2.2068,   1.9381,   1.8797,   1.7841,   1.7930
 cc     *,   1.2676,   1.6641,   1.4604,   1.0403,   1.3470
 cc     *,   1.0420,    .8722,    .8368,    .5323,    .5352
 cc     *,    .8867,    .5382,    .8106,    .8043,    .8058
 cc     *,    .8630,    .7031,    .7777,    .9441,   1.0403
 cc     *,    .9772,    .6628,   1.1968,   1.0555,   1.0495
 cc     */
 c      data pi1ecm/
 c     *     .45000,    .55000,    .62500,    .68750,    .73750
 c     *,    .78750,    .83750,    .88750,    .95000,   1.02500
 c     *,   1.10000,   1.17500
 c     */
 c      data pi1bmx/
 c     *     .48533,    .47873,    .55852,    .58087,    .55279
 c     *,    .49185,    .47203,    .46181,    .46181,    .47203
 c     *,    .41459,    .40684
 c     */
 c      data pi2ecm/
 c     *     .45000,    .55000,    .62500,    .68750,    .73750
 c     *,    .78750,    .83750,    .88750,    .95000,   1.02500
 c     *,   1.10000,   1.17500
 c     */
 c      data pi2bmx/
 c     *     .45135,    .80385,    .87767,   1.28779,   1.81771
 c     *,   2.07296,   1.55536,    .95912,    .94407,    .77768
 c     *,    .70467,    .54408
 c     */
 c      data pi3ecm/
 c     *     .45000,    .55000,    .61250,    .65000,    .67500
 c     *,    .70000,    .72500,    .75000,    .77500,    .80000
 c     *,    .82500,    .85000,    .87500,    .90000,    .93750
 c     *,   1.02500,   1.10000,   1.17500
 c     */
 c      data pi3bmx/
 c     *     .50463,    .79988,   1.13541,   1.30497,   1.34462
 c     *,   1.52957,   1.63809,   1.94216,   1.98511,   1.66316
 c     *,   1.69163,   1.62933,   1.06152,    .96574,    .83302
 c     *,    .71365,    .66517,    .59173
 c     */
 c      data pi4ecm/
 c     *     .28500,    .29500,    .30500,    .31500,    .37000
 c     *,    .46000,    .54375,    .61250,    .65000,    .67500
 c     *,    .70000,    .72500,    .75000,    .77500,    .80000
 c     *,    .82500,    .85000,    .87500,    .90125,    .92750
 c     *,    .96000
 c     */
 c      data pi4bmx/
 c     *     .29854,    .46181,    .54115,    .54115,    .51709
 c     *,    .72031,    .80976,    .96078,    .90798,   1.07788
 c     *,   1.37389,   1.70101,   2.02716,   1.83687,   1.52644
 c     *,   1.57469,   1.35640,    .93390,    .98531,    .84062
 c     *,    .79188
 c     */
 c      data pi5ecm/0.75,0.80,0.85,0.90,0.95,1.00,1.05,1.10,1.15/
 c      data pi5bmx/0.4,0.5,0.7,0.9,1.1,1.1,0.7,0.6,0.5/
 c
 c      call utpri('jintcs',ish,ishini,7)
 c
 c      ics=0
 c
 c      if(idptl(i).gt.10000)then
 c       call idquad(i,nqi,nai,jci)
 c       if(nqi.eq.0)then
 c        idi=999
 c       else
 c        idi=9999
 c       endif
 c      else
 c       idi=idptl(i)
 c      endif
 c
 c      if(idptl(j).gt.10000)then
 c       call idquad(j,nqj,naj,jcj)
 c       if(nqj.eq.0)then
 c        idj=999
 c       else
 c        idj=9999
 c       endif
 c      else
 c       idj=idptl(j)
 c      endif
 c
 c      do k=1,nflav
 c       kc(k)=jc(k,1)-jc(k,2)
 c       if(nq.lt.0)kc(k)=-kc(k)
 c      enddo
 c      if(nq.lt.0)nq=-nq
 c
 c      if(ish.ge.7)then
 c       write(ifch,*)'i:',i,' id_i:',idptl(i),' j:',j,' id_j:',idptl(j)
 c       write(ifch,*)'E_cm:',ecm,' jc:',(jc(k,1),k=1,6),(jc(k,2),k=1,6)
 c       write(ifch,*)'nq:',nq,' kc:',(kc(k),k=1,6)
 c      endif
 c
 cc check minimal kinetic energy
 c      ekin=ecm-pptl(5,i)-pptl(5,j)
 c      if(ekin.lt.amimel)goto1000
 c
 cc -------------------------------------------------------------------------
 c      if(iabs(idi).gt.1000.or.iabs(idj).gt.1000)then   !baryon involved
 cc -------------------------------------------------------------------------
 c
 c       if(nq.eq.6)then !------------baryon-baryon ----->
 c
 c        if(kc(1).eq.kc(2))then    !isospin_z=0
 cc pn
 c         if(ish.ge.7)write(ifch,*)'sig_pn chosen'
 c         if(ecm.lt.2)then
 c          call utindx(npn,pnecm,ecm,iecm)
 c          bmx=pnbmx(iecm)
 c         else
 c          p=ecm*sqrt((ecm**2/4./.94**2)-1.)
 c          if(p.lt.1.)then
 c           sig=33.+196.*(abs(0.95-p))**2.5
 c          elseif(p.lt.2)then
 c           sig=24.2+8.9*p
 c          else
 c           sig=47.3+.513*(alog(p)**2)-4.27*alog(p)
 c          endif
 c          bmx=sqrt(sig/10./pi)
 c         endif
 c        else                      !isospin_z=+-1
 cc pp
 c         if(ish.ge.7)write(ifch,*)'sig_pp chosen'
 c         if(ecm.le.1.882)then
 c          bmx=ppbmx(1)
 c          if(ish.ge.7)write(ifch,*)'b_mx:',bmx
 c         elseif(ecm.le.1.887)then
 c          bmx=ppbmx(2)
 c          if(ish.ge.7)write(ifch,*)'b_mx:',bmx
 c         elseif(ecm.le.1.893)then
 c          bmx=ppbmx(3)
 c          if(ish.ge.7)write(ifch,*)'b_mx:',bmx
 c         elseif(ecm.le.1.899)then
 c          bmx=ppbmx(4)
 c         elseif(ecm.le.1.909)then
 c          bmx=ppbmx(5)
 c         elseif(ecm.le.1.913)then
 c          bmx=ppbmx(6)
 c         elseif(ecm.le.1.918)then
 c          bmx=ppbmx(7)
 c         else
 c          p=ecm*sqrt((ecm**2/4./.94**2)-1.)
 c          if(p.lt.0.8)then
 c           sig=23.5+1000.*(0.7-p)**4
 c          elseif(p.lt.1.5)then
 c           sig=23.5+24.6/(1+exp(-(p-1.2)/0.1))
 c          elseif(p.lt.5)then
 c           sig=41.+60.*(p-0.9)*exp(-1.2*p)
 c          else
 c           sig=48+.522*(alog(p)**2)-4.51*alog(p)
 c          endif
 c          bmx=sqrt(sig/10./pi)
 c         endif
 c        endif
 c
 c       elseif(nq.eq.0)then !-----------baryon-antibaryon ---->
 c
 c        if(kc(1).eq.0.and.kc(2).eq.0.and.kc(3).eq.0)then  !p-ap, n-an
 cc app
 c         if(ish.ge.7)write(ifch,*)'sig_app chosen'
 c         if(ecm.lt.2)then
 c          call utindx(napp,appecm,ecm,iecm)
 c          bmx=appbmx(iecm)
 c         else
 c          p=ecm*sqrt((ecm**2/4./.94**2)-1.)
 c          sig=38.4+77.6*p**(-.64)+.26*(alog(p)**2)-1.2*alog(p)
 c          bmx=sqrt(sig/10./pi)
 c         endif
 c        else                                              !p-an, n-ap
 cc apn
 c         if(ish.ge.7)write(ifch,*)'sig_apn chosen'
 c         if(ecm.lt.2)then
 c          call utindx(napn,apnecm,ecm,iecm)
 c          bmx=apnbmx(iecm)
 c         else
 c          p=ecm*sqrt((ecm**2/4./.94**2)-1.)
 c          sig=133.6*p**(-.7)-1.22*(alog(p)**2)+13.7*alog(p)
 c          bmx=sqrt(sig/10./pi)
 c         endif
 c        endif
 c
 c       elseif(nq.eq.3)then !----------baryon-meson ---->
 c
 c        if(kc(3).eq.0)then
 cc no kaons involved(except for K-L)
 c         if(kc(1).eq.0.or.kc(2).eq.0)then
 cc pip
 c          if(ish.ge.7)write(ifch,*)'sig_pip chosen'
 c          if(ecm.lt.2.5)then
 c           call utindx(npip,pipecm,ecm,iecm)
 c           bmx=pipbmx(iecm)
 c          else
 c           p=sqrt(((ecm**2-0.94**2-0.14**2)/2./0.94)**2.-0.14**2)
 c           sig=16.4+19.3*p**(-.42)+.19*(alog(p)**2)
 c           bmx=sqrt(sig/10./pi)
 c          endif
 c         else
 cc pim
 c          if(ish.ge.7)write(ifch,*)'sig_pim chosen'
 c          if(ecm.lt.2.5)then
 c           call utindx(npim,pimecm,ecm,iecm)
 c           bmx=pimbmx(iecm)
 c          else
 c           p=sqrt(((ecm**2-0.94**2-0.14**2)/2./0.94)**2.-0.14**2)
 c           sig=33.0+14.0*p**(-1.36)+.456*(alog(p)**2)-4.03*alog(p)
 c           bmx=sqrt(sig/10./pi)
 c          endif
 c         endif
 c        elseif(kc(3).eq.1)then
 cc strange particles involved
 c         if(kc(1).eq.0.or.kc(2).eq.0)then
 cc kmn
 c          if(ish.ge.7)write(ifch,*)'sig_kmn chosen'
 c          if(ecm.lt.2.5)then
 c           call utindx(nkmn,kmnecm,ecm,iecm)
 c           bmx=kmnbmx(iecm)
 c          else
 c           p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2)
 c           sig=25.2+.38*(alog(p)**2)-2.9*alog(p)
 c           bmx=sqrt(sig/10./pi)
 c          endif
 c         else
 cc kmp
 c          if(ish.ge.7)write(ifch,*)'sig_kmp chosen'
 c          if(ecm.lt.2.5)then
 c           call utindx(nkmp,kmpecm,ecm,iecm)
 c           bmx=kmpbmx(iecm)
 c          else
 c           p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2)
 c           sig=32.1+.66*(alog(p)**2)-5.6*alog(p)
 c           bmx=sqrt(sig/10./pi)
 c          endif
 c         endif
 c        elseif(kc(3).eq.-1)then
 cc strange particles involved
 c         if(kc(1).eq.3.or.kc(2).eq.3)then
 cc kpp
 c          if(ish.ge.7)write(ifch,*)'sig_kpp chosen'
 c          if(ecm.lt.2.5)then
 c           call utindx(nkpp,kppecm,ecm,iecm)
 c           bmx=kppbmx(iecm)
 c          else
 c           p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2)
 c           sig=18.1+.26*(alog(p)**2)-alog(p)
 c           bmx=sqrt(sig/10./pi)
 c          endif
 c         else
 cc kpn
 c          if(ish.ge.7)write(ifch,*)'sig_kpn chosen'
 c           if(ecm.lt.2.5)then
 c           call utindx(nkpn,kpnecm,ecm,iecm)
 c           bmx=kpnbmx(iecm)
 c          else
 c           p=sqrt(((ecm**2-0.94**2-0.49**2)/2./0.94)**2.-0.49**2)
 c           sig=18.7+.21*(alog(p)**2)-.89*alog(p)
 c           bmx=sqrt(sig/10./pi)
 c          endif
 c         endif
 c        elseif(kc(3).ge.2)then
 cc two strange particles involved
 c         bmx=1.
 c        elseif(kc(3).le.-2)then
 cc two strange particles involved
 c        bmx=1.
 c        endif
 c
 c       else !-----------baryon_cluster-anything---->
 c
 c        write(ifch,*)'i:',i,' id_i:',idptl(i),' j:',j,' id_j:',idptl(j)
 c        write(ifch,*)'r_i:',radptl(i),' r_j:',radptl(j)
 c        write(ifch,*)'E_cm:',ecm,' jc:',(jc(k,1),k=1,6),(jc(k,2),k=1,6)
 c        write(ifch,*)'nq:',nq,' kc:',(kc(k),k=1,6)
 c        bmx=radptl(i)+radptl(j)
 c
 c       endif ! <--------------
 c
 cc -------------------------------------
 c      else  !   meson-meson
 cc -------------------------------------
 c
 c       call idquad(i,nqi,nai,jci)
 c       call idquad(j,nqj,naj,jcj)
 c       do l=1,nflav
 c        kci(l)=jci(l,1)-jci(l,2)
 c        kcj(l)=jcj(l,1)-jcj(l,2)
 c       enddo
 c
 c       if(kci(3).eq.0.and.pptl(5,i).le.1.0
 c     * .and.kcj(3).eq.0.and.pptl(5,j).le.1.0)then
 c
 c        if(kci(1).eq.0.and.kci(2).eq.0.and.pptl(5,i).le.0.5)then
 c         if(kcj(1).eq.0.and.kcj(2).eq.0.and.pptl(5,j).le.0.5)then
 cc pi0 pi0
 c          goto104
 c         elseif(kcj(1).eq.1.and.kcj(2).eq.-1)then
 cc pi0 pi+
 c          goto102
 c         elseif(kcj(1).eq.-1.and.kcj(2).eq.1)then
 cc pi0 pi-
 c          goto103
 c         else
 cc pi0 eta
 c          goto105
 c         endif
 c        elseif(kci(1).eq.1.and.kci(2).eq.-1)then
 c         if(kcj(1).eq.0.and.kcj(2).eq.0.and.pptl(5,j).le.0.5)then
 cc pi+ pi0
 c          goto102
 c         elseif(kcj(1).eq.1.and.kcj(2).eq.-1)then
 cc pi+ pi+
 c          goto101
 c         elseif(kcj(1).eq.-1.and.kcj(2).eq.1)then
 cc pi+ pi-
 c          goto104
 c         else
 cc pi+ eta
 c          goto105
 c         endif
 c        elseif(kci(1).eq.-1.and.kci(2).eq.1)then
 c         if(kcj(1).eq.0.and.kcj(2).eq.0.and.pptl(5,j).le.0.5)then
 cc pi- pi0
 c          goto103
 c         elseif(kcj(1).eq.1.and.kcj(2).eq.-1)then
 cc pi- pi+
 c          goto104
 c         elseif(kcj(1).eq.-1.and.kcj(2).eq.1)then
 cc pi- pi-
 c          goto101
 c         else
 cc pi- eta
 c          goto105
 c         endif
 c        else
 c         if(pptl(5,j).le.0.5)then
 cc eta pi
 c          goto105
 c         else
 cc eta eta
 c          goto106
 c         endif
 c        endif
 c101     continue
 c        if(ish.ge.7)write(ifch,*)'sig_pi+pi+  chosen'
 c        call utindx(npi1,pi1ecm,ecm,iecm)
 c        bmx=pi1bmx(iecm)
 c        goto110
 c102     continue
 c        if(ish.ge.7)write(ifch,*)'sig_pi+pi0  chosen'
 c        call utindx(npi2,pi2ecm,ecm,iecm)
 c        bmx=pi2bmx(iecm)
 c        goto110
 c103     continue
 c        if(ish.ge.7)write(ifch,*)'sig_pi-pi0  chosen'
 c        call utindx(npi3,pi3ecm,ecm,iecm)
 c        bmx=pi3bmx(iecm)
 c        goto110
 c104     continue
 c        if(ish.ge.7)write(ifch,*)'sig_pi-pi+  chosen'
 c        call utindx(npi4,pi4ecm,ecm,iecm)
 c        bmx=pi4bmx(iecm)
 c        goto110
 c105     continue
 c        if(ish.ge.7)write(ifch,*)'sig_pi_eta  chosen'
 c        call utindx(npi5,pi5ecm,ecm,iecm)
 c        bmx=pi5bmx(iecm)
 c        goto110
 c106     continue
 c        if(ish.ge.7)write(ifch,*)'sig_eta_eta  chosen'
 c        bmx=0.4   ! approx.  sig=5mb
 c110     continue
 c
 c       else !something else involved, strange etc.
 c
 c        bmx=0.6  ! approx.  sig=10mb
 c
 c       endif
 c
 cc --------------------------------
 c      endif
 cc --------------------------------
 c
 c      if(bij.le.bmx)ics=1
 c      if(ish.ge.7)then
 c      write(ifch,*)'b_mx:',bmx,' b_ij:',bij,' ics:',ics
 c      endif
 c
 c1000  continue
 c      call utprix('jintcs',ish,ishini,7)
 c      return
 c      end
 c
 cc----------------------------------------------------------------------
 c      subroutine jintel(i,j,p,amim,xaver)
 cc----------------------------------------------------------------------
 cc  elastic scattering of ptls i,j
 cc----------------------------------------------------------------------
 c      include 'epos.inc'
 c      real xaver(4)
 c      real p(5),u(5),pei(5),pej(5)
 c
 cc     determine momenta of outgoing particles (pei,pej)
 cc     -------------------------------------------------
 c           if(p(5).le.(pptl(5,i)+pptl(5,j))*.99)then
 c      if(ish.ge.1)then
 c      call utmsg('jintel')
 c      write(ifch,132)p(5),pptl(5,i)+pptl(5,j)
 c132   format(1x,'*****  m_fus < m_i+m_j ---> qcm set zero    ( '
 c     *,2f10.3,' )')
 c      write(ifch,133)'p_i:  ',(pptl(k,i),k=1,5)
 c      write(ifch,133)'p_j:  ',(pptl(k,j),k=1,5)
 c      write(ifch,133)'p_fus:',(p(k),k=1,5)
 c133   format(1x,a6,3x,5f10.3)
 c      call utmsgf
 c      endif
 c      qcm=0.
 c           elseif(p(5).le.pptl(5,i)+pptl(5,j))then
 c      qcm=0.
 c           else
 c      qcm=utpcm(p(5),pptl(5,i),pptl(5,j))
 c           endif
 c
 cc isotropic
 c      u(3)=2.*rangen()-1.
 c      phi=2.*pi*rangen()
 c      u(1)=sqrt(1.-u(3)**2)*cos(phi)
 c      u(2)=sqrt(1.-u(3)**2)*sin(phi)
 c      do 47 k=1,3
 c      pei(k)= qcm*u(k)
 c47    pej(k)=-qcm*u(k)
 c
 cc nonisotropic
 cc-c   pt=sqrt(2./pi)*ptq*sqrt(-2*alog(rangen()) !2-dim Gauss
 cc-c   if(pt.ge.qcm)pt=rangen()*qcm
 cc-c   qpl=sqrt(qcm**2-pt**2)
 cc-c   u(3)=qpl
 cc-c   phi=2.*pi*rangen()
 cc-c   u(1)=pt*cos(phi)
 cc-c   u(2)=pt*sin(phi)
 cc-c   call utaxis(i,j,a1,a2,a3)
 cc-c   ivt=1
 cc-c   if(a3.lt.0.)then
 cc-c   a1=-a1
 cc-c   a2=-a2
 cc-c   a3=-a3
 cc-c   ivt=-1
 cc-c   endif
 cc-c   call utrota(-1,a1,a2,a3,u(1),u(2),u(3))
 cc-c   do 47 k=1,3
 cc-c   pei(k)= u(k)*ivt
 cc-c47 pej(k)=-u(k)*ivt
 c
 c      pei(4)=sqrt(qcm**2+pptl(5,i)**2)
 c      pej(4)=sqrt(qcm**2+pptl(5,j)**2)
 c      pei(5)=pptl(5,i)
 c      pej(5)=pptl(5,j)
 c      call utlobo(-1,p(1),p(2),p(3),p(4),p(5)
 c     *,pei(1),pei(2),pei(3),pei(4))
 c      call utlobo(-1,p(1),p(2),p(3),p(4),p(5)
 c     *,pej(1),pej(2),pej(3),pej(4))
 c
 cc     fill /cptl/
 cc     -----------
 c      do 49 lo=1,2
 c      nptl=nptl+1
 c      if(lo.eq.1)ij=i
 c      if(lo.eq.2)ij=j
 c      do 48 k=1,5
 c      if(lo.eq.1)pptl(k,nptl)=pei(k)
 c      if(lo.eq.2)pptl(k,nptl)=pej(k)
 c48    continue
 c      istptl(nptl)=0
 c      idptl(nptl)=idptl(ij)
 c      ibptl(1,nptl)=ibptl(1,ij)
 c      ibptl(2,nptl)=ibptl(2,ij)
 c      ibptl(3,nptl)=ibptl(3,ij)
 c      ibptl(4,nptl)=ibptl(4,ij)
 c      xorptl(1,nptl)=xaver(1)
 c      xorptl(2,nptl)=xaver(2)
 c      xorptl(3,nptl)=xaver(3)
 c      xorptl(4,nptl)=xaver(4)
 c      iorptl(nptl)=i
 c      jorptl(nptl)=j
 c      tivptl(1,nptl)=xaver(4)
 c      call idtau(idptl(nptl),pptl(4,nptl),pptl(5,nptl),taugm)
 c      tivptl(2,nptl)=tivptl(1,nptl)+taugm*(-alog(rangen()))
 c      ifrptl(1,nptl)=0
 c      ifrptl(2,nptl)=0
 c      ityptl(nptl)=ityptl(ij)
 c49    continue
 c
 c1000  return
 c      end
 c
 cc----------------------------------------------------------------------
 c      subroutine jintfs(i,j,p,nq,jc,amim,iret)
 cc----------------------------------------------------------------------
 cc  input:
 cc    i,j: particle indices
 cc  output:
 cc    p: 5-momentum of fused ptl
 cc    nq: net quark number of fused ptl
 cc    jc: jc of fused ptl
 cc    amim : minimum mass of fused ptl
 cc    iret: return code
 cc            0=ok
 cc            1=mass p(5) less than amim
 cc----------------------------------------------------------------------
 c      include 'epos.inc'
 c      integer jci(nflav,2),jcj(nflav,2),jc(nflav,2)
 c      real p(5)
 c      double precision ppfus(4),pp52
 c
 c      iret=0
 c
 c      do 35 k=1,4
 c       p(k)=pptl(k,i)+pptl(k,j)
 c35    ppfus(k)=dble(p(k))
 c      pp52=ppfus(4)**2-ppfus(3)**2-ppfus(2)**2-ppfus(1)**2
 c      if(pp52.le.0.)then
 c       if(ish.ge.1)then
 c        call utmsg('jintfs')
 c        write(ifch,*)'*****  mfus**2 < 0    (',pp52,' )'
 c        write(ifch,*)(ppfus(m),m=1,4)
 c        call utmsgf
 c       endif
 c       goto1001
 c      endif
 c      p(5)=sngl(dsqrt(pp52))
 c
 c      call idquad(i,nqi,nai,jci)
 c      call idquad(j,nqj,naj,jcj)
 c      do 29 n=1,nflav
 c       do 29 k=1,2
 c29    jc(n,k)=jci(n,k)+jcj(n,k)
 c      nq=0
 c      do 54 n=1,nflav
 c54    nq=nq+jc(n,1)-jc(n,2)
 c
 c      call idcomj(jc)
 c      amim=utamnz(jc,5)+amimfs
 c      if(p(5).lt.amim)goto1001
 c      goto1000
 c1001  iret=1
 c1000  return
 c      end
 c
 cc----------------------------------------------------------------------
 c      subroutine jintfu(i,j,p,jc)
 cc----------------------------------------------------------------------
 cc input:
 cc   i,j: particle indices
 cc   p,jc: momentum and jc of fused object
 cc outout:
 cc   id, ib, ist, ior, jor, ifr, ity of fused particle
 cc   written to /cptl/ after increasing nptl by 1
 cc----------------------------------------------------------------------
 c
 c      include 'epos.inc'
 c      parameter (mxlook=10000,mxdky=2000)
 c      common/dkytab/look(mxlook),cbr(mxdky),mode(5,mxdky)
 c      real p(5)
 c      integer jc(nflav,2),ic(2),ib(4)
 c
 c      nptl=nptl+1
 c
 cc momentum
 c      do k=1,5
 c       pptl(k,nptl)=p(k)
 c      enddo
 c      amf=p(5)
 c
 cc determine idr, ib(1-4)
 c      idr=0
 c      do 40 nf=1,nflav
 c      do 40 ij=1,2
 c      if(jc(nf,ij).ge.10)idr=7*10**8
 c40    continue
 c           if(idr/10**8.ne.7)then
 c      call idenco(jc,ic,ireten)
 c      if(ireten.eq.1)call utstop('jintfu: idenco ret code = 1&')
 c      id=idtra(ic,0,0,3)
 c43    amc=amf
 c      call idres(id,amc,idr,iadj)
 c          if(idr.ne.0)then
 c      lid=look(iabs(idr))
 c      if((lid.le.0.or.lid.gt.0.and.mode(2,lid).eq.0)
 c     *.and.pptl(5,nptl).gt.amc+1e-3)then
 c      amf=amf+0.010
 c      goto43
 c      endif
 c      if((lid.le.0.or.lid.gt.0.and.mode(2,lid).eq.0)
 c     *.and.abs(amc-pptl(5,nptl)).gt.1e-3)then
 c      if(ish.ge.1)then
 c      call utmsg('jintfu')
 c      write(ifch,*)'*****  not on mass shell after fusion: '
 c     *,pptl(5,nptl),amc
 c      call utmsgf
 c      endif
 c      endif
 c           endif
 c           if(idr.eq.0)then
 c      if(mod(ic(1),100).ne.0.or.mod(ic(2),100).ne.0)then
 c      idr=9*10**8
 c      else
 c      idr=8*10**8+ic(1)*100+ic(2)/100
 c      endif
 c           endif
 c           else
 c      call idtrbi(jc,ib(1),ib(2),ib(3),ib(4))
 c      idr=idr
 c     *+mod(jc(1,1)+jc(2,1)+jc(3,1)+jc(4,1),10**4)*10**4
 c     *+mod(jc(1,2)+jc(2,2)+jc(3,2)+jc(4,2),10**4)
 c      if(ish.ge.7)write(ifch,*)'ib:',(ib(kk),kk=1,4)
 c           endif
 c
 c      n=nptl
 c
 cc     fill /cptl/
 cc     -----------
 c      idptl(n)=idr
 c      do k=1,4
 c      ibptl(k,n)=ib(k)
 c      enddo
 c      istptl(n)=0
 c      iorptl(n)=i
 c      jorptl(n)=j
 c      ifrptl(1,n)=0
 c      ifrptl(2,n)=0
 c      ityptl(n)=50
 c
 c      return
 c      end
 c
 c----------------------------------------------------------------------
       subroutine jintpo(lclean,iret)
 c----------------------------------------------------------------------
 c  parton-ladder-fusion -- 3D grid -- based on string segments
 c----------------------------------------------------------------------
       include 'epos.inc'
       include 'epos.incems'
       include 'epos.incico'
       common/cxyzt/xptl(mxptl),yptl(mxptl),zptl(mxptl),tptl(mxptl)
      *,optl(mxptl),uptl(mxptl),sptl(mxptl),rptl(mxptl,3)
       parameter(m1grid=10,kgrid=3,kegrid=3)
       parameter(m3xgrid=21*kgrid*kegrid)
       parameter(mxcl=4000,mxcli=1000)
       integer idropgrid(m1grid,m1grid,m3xgrid)
      &       ,jdropgrid(m1grid,m1grid,m3xgrid)
      &       ,jclu(m3xgrid),nsegp4(mxcl)
      &       ,irep(mxcl),mmji(mxcl,mxcli)
      &       ,jccl(mxcl,nflav,2),nseg(mxcl),mseg(mxcl),kclu(mxcl)
      &       ,naseg(0:mxcl),nfseg(mxcl),nsegmx(mxcl)
      &       ,jc(nflav,2),ke(6),jcjj(nflav,2)
      &       ,nsegmt(mxptl),kmean(2,1001:1000+mxcl)
      &       ,ncluj(1001:1000+mxptl)   !,ic(2),nclk(m3xgrid)
       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
      &                ,pptld(5,mxcl),p4tmp
       double precision enp,pzp,ena,pza,ccc,sss,enew,pznew,taa2,pta2
      .,pold,pnew,pnew2,p5a2,eeta,enf,pzf 
       common   /cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat
       common/cdelzet/delzet,delsgr /cvocell/vocell
       common/cranphi/ranphi
       parameter(maxp=40*mxcl)
       dimension yrad(maxp),phirad(maxp),pe(5),yrad2(maxp),nfrag(mxptl)
       logical first,lnew(m1grid,m1grid),lold(m1grid,m1grid),lcont,lpass
      &,lclean
       double precision ptest(5),ttest,p52,p4mean(4,mxcl),zor,tor,xmxms
      &,ppp(5),be(4),energ,bp,rmean(2,mxcl)         !,qqq(5),www(5)
       dimension mpair(mamx,mamx)
       save ncntpo!,ntrr,ntrt
       data ncntpo /0/!,ntrr/0/,ntrt/0/
       ncntpo=ncntpo+1
       call utpri('jintpo',ish,ishini,4)
       tauico=ttaus
 
       if(kigrid.gt.kgrid)then
         write(ifmt,*)'kigrid,kgrid:',kigrid,kgrid
         stop'jintpo: kigrid too big.          '
       endif
 
       fegrid=yhaha/5.36                 !rapidity range/rap range at RHIC
       if(fegrid.gt.kegrid)then
         write(ifmt,*)'fegrid,kegrid:',fegrid,kegrid
         stop'jintpo: kegrid too small.       '
       endif
 
       m3grid=m3xgrid /kgrid *kigrid /kegrid *fegrid 
 c      if(iLHC.eq.1)m3grid=min(m3xgrid,int(m3grid*0.75))
       if(mod(m3grid,2).eq.0)m3grid=m3grid+1
 
       xgrid=8
       sgrid=fsgrid *ttaus  *fegrid *6.0 
       delxgr=2*xgrid/m1grid
       delsgr=2*sgrid/m3grid
       shiftx=(1.-2.*rangen())*rangen()*delxgr/2.
       shifts=(1.-2.*rangen())*rangen()*delsgr/2.
       vocell= delxgr * delxgr * delsgr
       delzet=delsgr/ttaus
       nptla=nptl
       iflhc=1
       if(iLHC.eq.1)iflhc=min(3,max(1,nsegce-1))
       nsegsuj=max(nsegce,nint(float(nsegsu/iflhc)*fegrid))
       xmxms=150d0      !maximum mass for a subcluster
 
       if(ncntpo.eq.1)then
         write(ifmt,*)'EPOS used with FUSION option'
         if(ish.ge.1)then
       print*,'+',('-',i=1,57)
       print*,'| ttaus sgrid nsegce:',ttaus,fsgrid,'  ',nsegce
       print*,'| sgrid  delxgr delsgr vocell:',sgrid,delxgr,delsgr,vocell
       print*,'+',('-',i=1,57)
         endif
       endif
 c      print*,'+',('-',i=1,57)
 
 c...compute x,y,z
 
       if(ish.ge.6)write(ifch,*)'compute x,y,z'
       do n=1,nptla
         iaaptl(n)=1
         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         ! the following check is important for  ioclude= 4 or 5 
         !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         if(ioclude.ne.3.and.istptl(n).eq.10)then
         stop'\n\n remnant cluster detected in jintpo\n\n'
         endif
         if(istptl(n).eq.0.or.istptl(n).eq.10)then
           call jtain(n,xptl(n),yptl(n),zptl(n),tptl(n),nnn,0)
           call idtau(idptl(n),pptl(4,n),pptl(5,n),taurem)
           if(abs(taurem).gt.1e-6.and.
      .    ityptl(n).ge.40.and.ityptl(n).le.59.and.nnn.eq.1)then
            iaaptl(n)=0         !!!!nnn=1: ptl lives later
 c           print*,n,ityptl(n),idptl(iorptl(n)),istptl(iorptl(n)),nnn 
           endif
           call jtaus(zptl(n),tz,sptl(n))
           strap=1e10
           xpl=tptl(n)+zptl(n)
           xmi=tptl(n)-zptl(n)
           if(xmi.gt.0.0.and.xpl.gt.0.0)then
             strap=0.5*log(xpl/xmi)
           else           !particle at eta=inf
             iaaptl(n)=0    
           endif
           dezptl(n)=strap       !space-time-rapidity
 c          write(ifch,*)'dez ???',n,nnn,taurem,ityptl(n)
 c     & ,pptl(4,n),pptl(3,n),idptl(n),xptl(n),yptl(n),zptl(n),sptl(n)
         else
           iaaptl(n)=0    
           xptl(n)=0.
           yptl(n)=0.
           zptl(n)=0.
           sptl(n)=0.
           dezptl(n)=1e10
         endif
       enddo
       ntry=0
 
 c...valid particles
 
 c      print *,'start -----------------------------------------'
       ptmax=ptclu**2
       nptlo=1
       if(iLHC.eq.1)nptlo=maproj+matarg+1
       do n=1,nptla
        !random number calls before if's
        !to keep random number sequence for testing
         rdm=rangen()
         if(istptl(n).eq.0)then
           call idquac(n,idum1,idum2,idum3,jc)
         else
           jc(4,1)=0
           jc(4,2)=0
         endif
        if(iaaptl(n).ne.0)then
         if(istptl(n).ne.10)then
           pt2=pptl(1,n)**2+pptl(2,n)**2
           am2tmp=(pptl(4,n)+pptl(3,n))*(pptl(4,n)-pptl(3,n))-pt2
           if(n.le.nptlo)then         !spectators
             iaaptl(n)=0
           elseif(ityptl(n).eq.47.or.ityptl(n).eq.57)then
             iaaptl(n)=0        !no active spectators
           elseif(istptl(n).ne.0.or.ityptl(n).ge.60)then
             iaaptl(n)=0
           elseif(abs(idptl(n)).le.100)then
             iaaptl(n)=0         !no fundamental particle (electron...)
          !elseif(ityptl(n).eq.41.or.ityptl(n).eq.51)then
          !  iaaptl(n)=0     !avoid particles coming from remnant droplet decay
          !                       !(already droplet decay products !)
           elseif(abs(am2tmp-pptl(5,n)*pptl(5,n)).gt.30.)then
             iaaptl(n)=0         !to discard off shell particles
           elseif(jc(4,1).ne.0.or.jc(4,2).ne.0)then
             iaaptl(n)=0         !to discard charmed particles
           elseif(iLHC.eq.1)then
             if(pt2.lt.1.e-3)then
               iaaptl(n)=0           !to discard slow proton (spectators)
             elseif(ptclu.gt.0.)then
              if(pt2.gt.(1.5*ptclu)**2)then
               ptmax=max(ptmax,pt2)
               if(ioquen.eq.0.or.fploss.le.0.)then  !high pt particle escape
                 iaaptl(n)=0
               elseif(fploss.gt.1.e10)then   !high pt particle absorbed
                 iaaptl(n)=1
               else                       
                 iaaptl(n)=-1                !high pt particle lose energy
               endif
              elseif(pt2.gt.(0.5*ptclu)**2)then
               if(rdm.lt.(sqrt(pt2)-0.5*ptclu)/ptclu)then
                 ptmax=max(ptmax,pt2)
                 if(ioquen.eq.0.or.fploss.le.0.)then !high pt particle escape
                   iaaptl(n)=0
                 elseif(fploss.gt.1.e10)then !high pt particle absorbed
                   iaaptl(n)=1
                 else                       
                   iaaptl(n)=-1  !high pt particle lose energy
                 endif
               endif
              endif
             else
               if(ioquen.eq.0.or.fploss.le.0.)then  !high pt particle escape
                 iaaptl(n)=0
               elseif(fploss.gt.1.e10)then   !high pt particle absorbed
                 iaaptl(n)=1
               else                       
                 iaaptl(n)=-1                !high pt particle lose energy
               endif
             endif
           elseif(pt2.gt.(1.5*ptclu)**2)then
             ptmax=max(ptmax,pt2)
             if(maproj.lt.20.or.matarg.lt.20.or.ioquen.eq.0)then
               iaaptl(n)=0
             endif
           elseif(pt2.gt.(0.5*ptclu)**2)then
             if(rdm.lt.(sqrt(pt2)-0.5*ptclu)/ptclu)then
               if(maproj.lt.20.or.matarg.lt.20.or.ioquen.eq.0)iaaptl(n)=0
             endif
           endif
         endif
        endif
 c        if(iaaptl(n).eq.0.and.abs(idptl(n)).le.10000.and.abs(idptl(n))
 c     &    .ge.100.and.mod(abs(idptl(n)),100).ne.0
 c     &   .and.ityptl(n).ne.0.and.ityptl(n).lt.50)
 c     &  print*,'not selected',idptl(n),ityptl(n),pptl(4,n)
       enddo
 
 c...Start cluster formation
 
  8888 continue
 
       nsegsuj=max(nsegce,nint(nsegsuj/1.08**ntry))
       ntry=ntry+1
       if(ntry.gt.90)        !do not put more than 100 or change limit for p4mean
      &call utstop('jintpo: cluster formation impossible ! &')
       nptl=nptla
 
       do 1 k=1,m3grid
       do 1 j=1,m1grid
       do 1 i=1,m1grid
       idropgrid(i,j,k)=0
   1   continue
 
 c...count string segments in cell
 
       if(ish.ge.6)write(ifch,*)'count string segments in cell'
       do n=1,nptla
         if(iaaptl(n).ne.0)then
           i=1+(xptl(n)+xgrid+shiftx)/delxgr
           j=1+(yptl(n)+xgrid+shiftx)/delxgr
           k=1+(sptl(n)+sgrid+shifts)/delsgr
           if(  i.ge.1.and.i.le.m1grid
      &    .and.j.ge.1.and.j.le.m1grid
      &    .and.k.ge.1.and.k.le.m3grid)then
             nfrag(n)=1
             if(iLHC.eq.1)then  !count number of quarks to absorbe more baryon than mesons
               call idquac(n,idum1,idum2,idum3,jc)
               nfrag(n)=0
               do nj=1,2
                 do ni=1,nflav
                   nfrag(n)=nfrag(n)+ni*jc(ni,nj)
                 enddo
               enddo
               
             endif
 c            print *,idptl(n),nfrag(n)
             idropgrid(i,j,k)=idropgrid(i,j,k)+nfrag(n)
           endif
         endif
       enddo
 
 
       if(iLHC.eq.1)then
 cc...check high pt segments
 c      !to use this part one has to define:
 c      !...  1 = valid
 c      !... -1 = valid but high pt
 c      !...  0 = not valid
 c
 
 c      esu=0
 c      do i=1,nptla
 c      if(istptl(i).eq.0.or.istptl(i).eq.3)esu=esu+pptl(4,i)
 c      enddo
 c      write(ifmt,'(a,41x,f15.1)')' +++++Eajint+++++',esu
 
       if(ish.ge.6)write(ifch,*)'check high pt segments'
       if(fploss.gt.0.0)then
       ein=0
       elo=0
       do n=1,nptla
         delen=0.
         if(iaaptl(n).eq.-1)then
           i=1+(xptl(n)+xgrid+shiftx)/delxgr
           j=1+(yptl(n)+xgrid+shiftx)/delxgr
           k=1+(sptl(n)+sgrid+shifts)/delsgr
           if(   i.ge.1.and.i.le.m1grid
      &     .and.j.ge.1.and.j.le.m1grid
      &     .and.k.ge.1.and.k.le.m3grid)then
            iescape=0  
            xa=xptl(n)
            ya=yptl(n)
            eta=dezptl(n)
            pz=pptl(3,n)
            en=pptl(4,n)
            taa2=pptl(1,n)**2+pptl(2,n)**2+pptl(5,n)**2
            pza=pz   
            ena=sqrt(taa2+pza**2)
            enaxx=en  
            p5a2=pptl(5,n)**2
            pta2=pptl(1,n)**2+pptl(2,n)**2
            !transform to bjo
            eeta=eta
            ccc=cosh(eeta)  
            sss=sinh(eeta)  
            enp= ena*ccc-pza*sss
            pzp=-ena*sss+pza*ccc
            vz=pzp/enp
            vx=pptl(1,n)/enp
            vy=pptl(2,n)/enp
            !print*,'+++++++',eta,vz,pz/en
            if(vx.ne.0.0.or.vy.ne.0.0)then
              if(abs(vx).ge.abs(vy))then
                ica=1
                rat=vy/vx
                is=sign(1.,vx)
                l=i
                va=xa
                wa=ya
              else
                ica=2
                rat=vx/vy
                is=sign(1.,vy)
                l=j
                va=ya
                wa=xa
              endif
              if(is.eq.-1)then
                imax=1
              else                 !if(is.eq. 1)
                imax=m1grid
              endif
              vr=sqrt(vx**2+vy**2)
              ratx=vz/vr
              qu=fploss*delxgr*sqrt(1+rat**2)*sqrt(1+ratx**2)
              do lu=l,imax,is
               vu=-(xgrid+shiftx)+(lu-0.5)*delxgr
               wu=wa+rat*(vu-va)
               mu=1+(wu+xgrid+shiftx)/delxgr
               if(mu.ge.1.and.mu.le.m1grid)then
                 if(ica.eq.1)then
                   ix=lu
                   jx=mu
                 else
                   ix=mu
                   jx=lu
                 endif 
                 if(idropgrid(i,j,k).ge.nsegce)then
                   dens=float(idropgrid(ix,jx,k))/float(iflhc)
                 else
                   dens=0
                 endif
                 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                 ! DeltaE ~ sqrt(T**3*E) * dL (BDMPS2008)
                 ! -> DeltaE ~ dens^3/8 * sqrt(E) *dL
                 !del=dens**(3./8.)*max(1.,sqrt(pptl(4,n))) 
                 !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                 escale=6
                 xen=enp/escale
                 del=0.02*dens**(3./8.)*max(1.,sqrt(xen)) 
                 delen=delen+del*escale*qu
 c                if(k.ge.m3grid/2-2.and.k.le.m3grid/2+4)
 c     &               write(ifch,'(2i3,4x,2i3,4x,2i3,4x,i3,3f7.3)')
 c     &               k,ica,i,j,ix,jx,idropgrid(ix,jx,k),qu,delen
 c     &               ,pptl(4,n)-delen
                endif
              enddo
            endif
            pold=sqrt(pta2+pzp**2)
            enew=enp-delen
            pnew2=enew**2-p5a2
            if(enew.gt.0..and.pnew2.gt.0.)iescape=1  
            if(iescape.eq.1)then  
               iaaptl(n)=0
               idropgrid(i,j,k)=idropgrid(i,j,k)-nfrag(n)
               if(idropgrid(i,j,k).lt.0)stop'jintpo: not possible.    '
               pnew=sqrt(pnew2)
               p1=pnew*pptl(1,n)/pold
               p2=pnew*pptl(2,n)/pold
               pznew=pnew*pzp   /pold      
               enf= enew*ccc+pznew*sss
               pzf= enew*sss+pznew*ccc
               if(ish.ge.1.and.enf.gt.ena.or.enf.ne.enf)then
                 if(ish.ge.1.and.(enf-ena.gt.0.1d0.or.enf.ne.enf))then
                 write(ifmt,*)'WWWW en gain for eta =',eta
                 write(ifmt,*)'e pz bf (lab)',ena,pza,enaxx
                 write(ifmt,*)'e pz bf (bjo)',enp,pzp   ,enp**2-pzp**2
                 write(ifmt,*)'e pz af (bjo)',enew,pznew,enew**2-pznew**2
                 write(ifmt,*)'e pz af (lab)',enf,pzf   ,enf**2-pzf**2
 c                stop
                 endif
               elseif(delen.gt.1e-3)then
 c create fake particle with energy lost in cluster
                 nptl=nptl+1   
                 nptla=nptla+1
                 istptl(nptl)=3  !daughter part of cluster
                 iaaptl(nptl)=1
                 pptl(1,nptl)=pptl(1,n)-p1
                 pptl(2,nptl)=pptl(2,n)-p2
                 pptl(3,nptl)=pptl(3,n)-pzf
                 pptl(4,nptl)=pptl(4,n)-enf
                 do l=1,3
                   xorptl(l,nptl)=xorptl(l,n)
                 enddo
           ! mother
                 pptl(1,n)=p1
                 pptl(2,n)=p2
                 pptl(3,n)=pzf     
                 tivptl(1,n)=xorptl(4,n)
                 pptl(4,n)=sqrt(pptl(1,n)**2+pptl(2,n)**2
      &                        +pptl(3,n)**2+pptl(5,n)**2)
                 if(abs(pptl(4,n)-enf).gt.0.01*pptl(4,n))
      &   print *,'jintpo eloss',pptl(4,n),enf,(pptl(4,n)-enf)/pptl(4,n)
                 call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm)
                 tivptl(2,n)=tivptl(1,n)+taugm*(-alog(rangen()))
                 ifrptl(1,n)=nptl
           ! daughter
                 pptl(5,nptl)=0.
                 pptl(4,nptl)=sqrt(pptl(1,nptl)**2+pptl(2,nptl)**2
      &                           +pptl(3,nptl)**2+pptl(5,nptl)**2)
                 xorptl(4,nptl)=xorptl(4,n)
                 tivptl(1,nptl)=tivptl(1,n)
                 iorptl(nptl)=n
                 jorptl(nptl)=0
                 idptl(nptl)=idptl(n)*100+sign(99,idptl(n))
                 ityptl(nptl)=ityptl(n)+2
                 xptl(nptl)=xptl(n)
                 yptl(nptl)=yptl(n)
                 zptl(nptl)=zptl(n)
                 sptl(nptl)=sptl(n)
                 dezptl(nptl)=dezptl(n)
                 if(ish.ge.6)write(ifch,*)'--> Virtual particle : ',nptl
      & ,idptl(nptl),iorptl(nptl),ityptl(nptl),(pptl(kk,nptl),kk=1,5)
               endif
               elo= elo+ena-pptl(4,n)  
               !write(ifch,*)n,'      escape'
            else
               iaaptl(n)=1
               ein=ein+pptl(4,n) 
               !write(ifch,*)n,'core'
            endif
           endif
         endif
       enddo
       else
       do n=1,nptla
         if(iaaptl(n).eq.-1)then
           iaaptl(n)=0
         endif
       enddo
       endif
 
 c      eloss=esu
 c      esu=0
 c      do i=1,nptla
 c      if(istptl(i).eq.0.or.istptl(i).eq.3)esu=esu+pptl(4,i)
 c      enddo
 c      write(ifmt,'(a,41x,f15.1)')' +++++Eajint+++++',esu
 c      eloss=eloss-esu
 c      write(ifch,*)'+++++ein,elo,eloss ',ein,elo,eloss
 
 
 c      do n=1,nptla
 c        if(iaaptl(n).eq.-1)then
 c          i=1+(xptl(n)+xgrid+shiftx)/delxgr
 c          j=1+(yptl(n)+xgrid+shiftx)/delxgr
 c          k=1+(sptl(n)+sgrid+shifts)/delsgr
 c          if(  i.ge.1.and.i.le.m1grid
 c     &    .and.j.ge.1.and.j.le.m1grid
 c     &    .and.k.ge.1.and.k.le.m3grid)then
 c            if(idropgrid(i,j,k).lt.2*nsegce)then  !surface segments
 c              iaaptl(n)=0
 c              idropgrid(i,j,k)=idropgrid(i,j,k)-1
 c              if(idropgrid(i,j,k).lt.0)stop'jintpo: not possible.    '
 cc           else
 cc             iaaptl(n)=1
 c           endif
 c          endif
 c        endif
 c      enddo
 
       endif
 
 c...identify clusters
 
       ifirst=0
       if(ish.ge.6)write(ifch,*)'identify clusters'
       do k=1,m3grid !~~~~~~k-loop
         jjj=0
         do j=1,m1grid
           first=.true.
           do i=1,m1grid
             if(idropgrid(i,j,k).ge.nsegce)then
               if(first)then
                 ifirst=i
                 jjj=jjj+1
                 if(jjj.gt.mxcl)stop'jintpo: mxcl too small.   '
                 irep(jjj)=0
                 jj=jjj
                 first=.false.
               endif
               jdropgrid(i,j,k)=jj
               if(j.gt.1)then
                if(jdropgrid(i,j-1,k).ne.0)then
                 jjo=jdropgrid(i,j-1,k)
                 if(jjo.lt.jj)then
                   if(jj.eq.jjj)jjj=jjj-1
                   jjx=jj
                   jj=jjo
                   do ii=ifirst,i
                     jdropgrid(ii,j,k)=jj
                     if(jdropgrid(ii,j-1,k).eq.jjx)then
                       if(jjx.gt.jjj)jjj=jjj+1
                       jja=jjx
                       jjb=jj
   90                  continue
                       if(irep(jja).eq.0.or.irep(jja).eq.jjb)then
                         irep(jja)=jjb
                       else
                         mn=min(irep(jja),jjb)
                         mx=max(irep(jja),jjb)
                         irep(jja)=mn
                         jja=mx
                         jjb=mn
                         goto 90
                       endif
                     endif
                   enddo
                 elseif(jdropgrid(i,j-1,k).gt.jj)then
                   irep(jjo)=jj
                 endif
                endif
               endif
             else
               jdropgrid(i,j,k)=0
               first=.true.
             endif
           enddo
         enddo
         !~~~~cluster jj ---> cluster irep(jj)
         do jj=jjj,1,-1
          if(irep(jj).ne.0)then
            do j=1,m1grid
              do i=1,m1grid
                if(jdropgrid(i,j,k).eq.jj)jdropgrid(i,j,k)=irep(jj)
              enddo
            enddo
          endif
         enddo
         !~~~~~remove empty cluster indices
         jjjx=jjj
         jjj=0
         jj=0
         do jjx=1,jjjx
          if(irep(jjx).eq.0)then
            jj=jj+1
            jjj=jjj+1
          else
            do j=1,m1grid
              do i=1,m1grid
                if(jdropgrid(i,j,k).gt.jj)
      &           jdropgrid(i,j,k)=jdropgrid(i,j,k)-1
              enddo
            enddo
          endif
         enddo
         !~~~~~
         jclu(k)=jjj
       enddo !~~~~~~~~~~~~~~~~~ END k-loop
 
 c...add segments to avoid holes
 
       if(ish.ge.6)write(ifch,*)'add segments from holes'
       if(iohole.eq.1.and.iLHC.eq.1)then
       do 83 n=1,nptla
         if(iaaptl(n).eq.0)goto 83
         ihit=0
         i=1+(xptl(n)+xgrid+shiftx)/delxgr
         j=1+(yptl(n)+xgrid+shiftx)/delxgr
         k=1+(sptl(n)+sgrid+shifts)/delsgr
         if(    i.ge.1.and.i.le.m1grid
      &          .and.j.ge.1.and.j.le.m1grid
      &          .and.k.ge.1.and.k.le.m3grid)then
           jgr=jdropgrid(i,j,k)
           nplus=idropgrid(i,j,k)
           if(jgr.eq.0)then          
             isgi=1
             if(rangen().gt.0.5)isgi=-1
             isgj=1
             if(rangen().gt.0.5)isgj=-1
             isgk=1
             if(rangen().gt.0.5)isgk=-1
             do ii=i-isgi,i+isgi,isgi
             do jj=j-isgj,j+isgj,isgj
             do kk=k-isgk,k+isgk,2*isgk
               if(  ii.ge.1.and.ii.le.m1grid
      .        .and.jj.ge.1.and.jj.le.m1grid
      .        .and.kk.ge.1.and.kk.le.m3grid)then
           if(jdropgrid(ii,jj,kk).gt.0)nplus=nplus+1
 c                 nplus=idropgrid(ii,jj,kk)+iflhc
                  if(nplus.ge.nsegce)then
                    ihit=1
 c                   goto 84
                  endif
               endif
             enddo
             enddo
             enddo
           endif
         endif
 c   84   continue
         if(ihit.eq.1)then
           idropgrid(i,j,k)=nplus!idropgrid(i,j,k)+iflhc
 c add cell randomly to adjacent cluster if any
           isgi=1
           if(rangen().gt.0.5)isgi=-1
           isgj=1
           if(rangen().gt.0.5)isgj=-1
           do ii=i-isgi,i+isgi,2*isgi
             do jj=j-isgj,j+isgj,2*isgj
               if(  ii.ge.1.and.ii.le.m1grid
      .        .and.jj.ge.1.and.jj.le.m1grid)then
                 jjj=jdropgrid(ii,jj,k)
                 if(jjj.gt.0)then
                   jdropgrid(i,j,k)=jjj
                   goto 83
                 endif
               endif
             enddo
           enddo
 c if no adjacent cluster, create a new one
           jclu(k)=jclu(k)+1
           jdropgrid(i,j,k)=jclu(k)
         endif
   83  continue
       endif
 
       if(ish.ge.8)then
         write(ifch,*)'segment positions'
         do k=1,m3grid
         write(ifch,*)'k=',k,'  jclu=',jclu(k)
      &              ,'  s=',(k-1)*delsgr-sgrid-shifts
          do j=m1grid,1,-1
         write(ifch,'(10i4,3x,10i4)')(idropgrid(i,j,k),i=1,m1grid)
      &                ,(jdropgrid(i,j,k),i=1,m1grid)
         enddo
         enddo
       endif
 
 c...absolute clusters numbering (for all k)
 
       if(ish.ge.6)write(ifch,*)'absolute clusters numbering'
       if(iLHC.eq.-1)then    !fuse touching clusters along k  (not used)
 c        midbin=m3grid/2
 c        if(rangen().gt.0.5)midbin=midbin+1
         jjj=1000
         k=1
         jclu(k)=0
         do j=1,m1grid
           do i=1,m1grid
             if(jdropgrid(i,j,k).gt.0)then
               jclu(k)=jclu(k)+1
               jjj=jjj+1
               ncluj(jjj)=0
               kmean(1,jjj)=k
               kmean(2,jjj)=k
               jjx=jdropgrid(i,j,k)
               do jj=j,m1grid
                 do ii=i,m1grid
                   if(jdropgrid(ii,jj,k).eq.jjx)then
                     jdropgrid(ii,jj,k)=jjj
                     ncluj(jjj)=ncluj(jjj)+idropgrid(ii,jj,k)
                   endif
                 enddo
               enddo
             endif
           enddo
         enddo
         do k=2,m3grid
           do j=1,m1grid
             do i=1,m1grid
               lnew(i,j)=.true.
               lold(i,j)=.true.
             enddo
           enddo
           jclu(k)=0
           do j=1,m1grid
             do i=1,m1grid
               if(lnew(i,j).and.jdropgrid(i,j,k).gt.0)then
                 lcont=.false.
                 if(jdropgrid(i,j,k-1).gt.0)then
                   jj0=jdropgrid(i,j,k-1)
                   nlim=0
                   if(jj0.gt.1000)nlim=ncluj(jj0)
                   if(nlim.lt.nsegsuj)then
                     jclu(k)=jclu(k)+1
                     jjx=jdropgrid(i,j,k)
                     if(jjx.gt.jj0)then
                       do jj=jjx,jjj
                         ncluj(jj)=ncluj(jj+1)
                         kmean(1,jj)=kmean(1,jj+1)
                         kmean(2,jj)=kmean(2,jj+1)
                       enddo
                       jjj=jjj-1
                       jclu(k)=jclu(k)-1
                     endif
                     kmean(2,jj0)=k
                     do jj=1,m1grid
                       do ii=1,m1grid
                         if(lnew(ii,jj))then
                           if(jdropgrid(ii,jj,k).eq.jjx)then
                             jdropgrid(ii,jj,k)=jj0
                             ncluj(jj0)=ncluj(jj0)+idropgrid(ii,jj,k)
                             lnew(ii,jj)=.false.
                    elseif(jdropgrid(ii,jj,k).gt.jjx.and.jjx.gt.jj0)then
                             jdropgrid(ii,jj,k)=jdropgrid(ii,jj,k)-1
                           endif
                         endif
                       enddo
                     enddo
                   else
                     lcont=.true.
                   endif
                 else
                   lcont=.true.
                 endif
                 if(lcont.and.lold(i,j))then
                   jclu(k)=jclu(k)+1
                   jjj=jjj+1
                   ncluj(jjj)=0
                   kmean(1,jjj)=k
                   kmean(2,jjj)=k
                   jjx=jdropgrid(i,j,k)
                   do jj=j,m1grid
                     do ii=i,m1grid
                       if(jdropgrid(ii,jj,k).eq.jjx)then
                         jdropgrid(ii,jj,k)=jjj
                         ncluj(jjj)=ncluj(jjj)+idropgrid(ii,jj,k)
                         lold(ii,jj)=.false.
                       endif
                     enddo
                   enddo
                 endif
               endif
             enddo
           enddo
         enddo
         jjs=0
         do jj=1001,jjj
           jjs=jjs+1
           kclu(jjs)=(kmean(1,jj)+kmean(2,jj))/2
           nseg(jjs)=0
           nsegp4(jjs)=0
           p4mean(1,jjs)=0d0   
           p4mean(2,jjs)=0d0   
           p4mean(3,jjs)=0d0   
           p4mean(4,jjs)=0d0   
           rmean(1,jjs)=0d0      !cluster distance to center
           rmean(2,jjs)=0d0      !cluster maximal radius
         enddo
         do k=1,m3grid
           do j=1,m1grid
             do i=1,m1grid
               if(jdropgrid(i,j,k).gt.0)then
                 jdropgrid(i,j,k)=jdropgrid(i,j,k)-1000
               endif
             enddo
           enddo
         enddo
         jjj=jjj-1000
       else
         jjj=jclu(1)
         do k=2,m3grid
           do j=1,m1grid
             do i=1,m1grid
               if(jdropgrid(i,j,k).gt.0)
      &             jdropgrid(i,j,k)=jdropgrid(i,j,k)+jjj
             enddo
           enddo
           jjj=jjj+jclu(k)
         enddo
         jjs=0
         do k=1,m3grid
           do jj=1,jclu(k)
             jjs=jjs+1
             kclu(jjs)=k
             nseg(jjs)=0
             nsegp4(jjs)=0
             p4mean(1,jjs)=0d0   
             p4mean(2,jjs)=0d0   
             p4mean(3,jjs)=0d0   
             p4mean(4,jjs)=0d0   
             rmean(1,jjs)=0d0    !cluster distance to center
             rmean(2,jjs)=0d0    !cluster maximal radius
           enddo
         enddo
       endif
       do 20 i=1,maproj
       do 20 j=1,matarg
  20   mpair(j,i)=0   
 
 
 c...calculate mean energy of segments going into in clusters
 
       if(ish.ge.6)write(ifch,*)
      &'calculate mean energy of segments going into in clusters'
       rmax=0.
       do 95 n=1,nptla
         if(iaaptl(n).eq.0)goto 95
         i=1+(xptl(n)+xgrid+shiftx)/delxgr
         j=1+(yptl(n)+xgrid+shiftx)/delxgr
         k=1+(sptl(n)+sgrid+shifts)/delsgr
         if(    i.ge.1.and.i.le.m1grid
      &          .and.j.ge.1.and.j.le.m1grid
      &          .and.k.ge.1.and.k.le.m3grid)then
           jj=jdropgrid(i,j,k)
           if(jj.gt.0)then
             nsegp4(jj)=nsegp4(jj)+1
             io=iorptl(n)
             if(iLHC.eq.1.and.ityptl(n).lt.40.and.io.gt.0)then !look for corresponding pair of nucleons
               if(iorptl(io).gt.0)then
                 do while(iorptl(iorptl(io)).ge.0)
                   io=iorptl(io)
                 enddo
                 ip=iorptl(io)
                 it=jorptl(io)-maproj
                 if(ip.gt.0.and.it.gt.0)mpair(it,ip)=mpair(it,ip)+1
               endif
             endif
             do kk=1,4
               p4mean(kk,jj)=p4mean(kk,jj)+dble(pptl(kk,n))
             enddo
             rr=sqrt(xptl(n)**2+yptl(n)**2)
             rmean(1,jj)=rmean(1,jj)+rr
             rmax=max(rr,rmax)
 c           if(jj.eq.9)write(ifch,*)'after',n
 c     & ,pptl(4,n),pptl(3,n),idptl(n),jj,nseg(jj),i,j,k
           elseif(istptl(n).eq.3)then   !virtual particle is lost (no cluster here)
             iaaptl(n) = 0   !don't use this particle next time
             istptl(n) = 5
             ior=iorptl(n)
             ifrptl(1,ior) = 0
             ifrptl(2,ior) = 0
             do k=1,3         !restore energy to mother particle
               pptl(k,ior)=pptl(k,ior)+pptl(k,n)
             enddo
             pptl(4,ior)=sqrt(pptl(1,ior)**2+pptl(2,ior)**2
      &                          +pptl(3,ior)**2+pptl(5,ior)**2)
             call idtau(idptl(ior),pptl(4,ior),pptl(5,ior),taugm)
             tivptl(2,ior)=tivptl(1,ior)+taugm*(-alog(rangen()))
           endif
         endif
  95   continue
 
       ectot=0.
       amctot=0.
       maptot=0
       mapmax=0
       npair=0
       do jj=1,jjj
         ectot=ectot+p4mean(4,jj)
         amctot=amctot+(p4mean(4,jj)+p4mean(3,jj))
      &     *(p4mean(4,jj)-p4mean(3,jj))-p4mean(1,jj)**2-p4mean(2,jj)**2
 
         if(nsegp4(jj).ge.nsegce/iflhc)then
           p4mean(4,jj)=p4mean(4,jj)/dble(nsegp4(jj))
           rmean(1,jj)=rmean(1,jj)/dble(nsegp4(jj))
         else
           p4mean(4,jj)=1d50
           rmean(1,jj)=-1d0
         endif
       enddo
       yco=0.
       ycor=1.
       if(iLHC.eq.1.and.amctot.gt.0.)then
         do ip=1,maproj
           do it=1,matarg
             if(mpair(it,ip).gt.0)npair=npair+1
             mapmax=max(mapmax,mpair(it,ip))
             maptot=maptot+mpair(it,ip)
           enddo
         enddo
         amctot=sqrt(amctot)
         if(amctot.gt.amuseg**2)then
           visco=1.
           yrmaxi=yradmx
 c          if(fvisco.gt.0)yrmaxi=yrmaxi*(1.-visco)
           yrmaxi=yrmaxi*log(exp(yradpi)+amctot**2)
 c          yrmaxi=yrmaxi*amctot**yradpi
           if(yrmaxi.gt.1e-2)then
             yyrmax=dble(yrmaxi)
             fradflii=sngl(1d0/
      &         ((sinh(yyrmax)*yyrmax-cosh(yyrmax)+1d0)/(yyrmax**2/2d0)))
           else
             yrmaxi=0.
             fradflii=1.
           endif
         else
           amctot=1.
           visco=0.
           fradflii=1.
           yrmaxi=ainfin         !to define ityptl as 19 if mass is too low (aumin=amuseg+yrmaxi in epos-hnb.f)
         endif
         if(ylongmx.lt.0.)then
           delzet=abs(ylongmx)*log(exp(yradmi)+amctot**2)
 c          delzet=abs(ylongmx)*log(exp(yradmi)+amctot)
 c          delzet=abs(ylongmx)*amctot**yradmi
           yco=delzet            !* 1.75
         else
           yco=ylongmx
         endif
         ycor=yco
         if(fvisco.gt.0.)ycor=0.
         if(ycor.gt.1.e-2)then
           ycor=sqrt(sinh(ycor)/ycor)/fradflii
         else
           ycor=1.
         endif
       else
         if(ylongmx.lt.0)delzet=1.75*delzet
       endif
 
 c...mark segments going into in clusters, count them
 
       if(ish.ge.6)write(ifch,*)'mark segments going into in clusters'
       do 96 n=1,nptla
         if(iaaptl(n).eq.0)goto 96
         i=1+(xptl(n)+xgrid+shiftx)/delxgr
         j=1+(yptl(n)+xgrid+shiftx)/delxgr
         k=1+(sptl(n)+sgrid+shifts)/delsgr
         if(    i.ge.1.and.i.le.m1grid
      &          .and.j.ge.1.and.j.le.m1grid
      &          .and.k.ge.1.and.k.le.m3grid)then
           jj=jdropgrid(i,j,k)
           if(jj.gt.0)then
           ! count only particles with reasonnable energy
             pt2=pptl(1,n)**2+pptl(2,n)**2+pptl(5,n)**2
             am2tmp=(pptl(4,n)+pptl(3,n))*(pptl(4,n)-pptl(3,n))-pt2
             if(abs(am2tmp).lt.0.1.or.
      &         dble(pptl(4,n)).le.100.d0/dble(ntry)*p4mean(4,jj))then
               istptl(n) = 3
               nseg(jj)=nseg(jj)+1
               if(rmean(1,jj).gt.0d0)then
                 rmean(2,jj)=max(rmean(2,jj)
      &                   ,abs(rmean(1,jj)-sqrt(xptl(n)**2+yptl(n)**2)))
 c      print *,'r',jj,rmean(1,jj),sqrt(xptl(n)**2+yptl(n)**2),rmean(2,jj)
               endif
 c           if(jj.eq.9)write(ifch,*)'after',n
 c     & ,pptl(4,n),pptl(3,n),idptl(n),jj,nseg(jj),i,j,k
             else
 c              write(ifmt,*)'Rejected particle : ',n,idptl(n)
              if(ish.ge.1)write(ifch,*)'Rejected particle : ',n,idptl(n)
      & ,ityptl(n),(pptl(kk,n),kk=1,5),am2tmp
      & ,dble(pptl(4,n))/dble(nsegp4(jj)),p4mean(4,jj),nsegp4(jj),jj
               nsegp4(jj)=nsegp4(jj)-1
               iaaptl(n)=0
 c        if(iaaptl(n).eq.0)print*,'rejected',idptl(n),ityptl(n)
 c     &                                         ,pptl(4,n)
               if(nsegp4(jj).ge.nsegce/iflhc)then
                 p4mean(4,jj)=(p4mean(4,jj)*dble(nsegp4(jj)+1)-pptl(4,n))
      &               /dble(nsegp4(jj))
               else
                 p4mean(4,jj)=1d50
               endif
             endif
          endif
         endif
   96  continue
 
 
 c...add segments moving towards clusters
 
       if(ish.ge.6)write(ifch,*)'add segments moving towards clusters'
       if(iocluin.eq.1)then
       do 93 n=1,nptla
         if(iaaptl(n).eq.0)goto 93
         ihit=0
         i=1+(xptl(n)+xgrid+shiftx)/delxgr
         j=1+(yptl(n)+xgrid+shiftx)/delxgr
         k=1+(sptl(n)+sgrid+shifts)/delsgr
         if(    i.ge.1.and.i.le.m1grid
      &          .and.j.ge.1.and.j.le.m1grid
      &          .and.k.ge.1.and.k.le.m3grid)then
           jgr=jdropgrid(i,j,k)
           if(jgr.eq.0)then
            if(i.ge.m1grid/2)then
             isi=-1
            else
             isi=1
            endif
            if(j.ge.m1grid/2)then
             jsi=-1
            else
             jsi=1
            endif
             do ii=i,i+2*isi,isi
              do jj=j,j+2*jsi,jsi
                if(.not.(ii.eq.i.and.jj.eq.j))then
                 if(ii.ge.1.and.ii.le.m1grid
      .          .and.jj.ge.1.and.jj.le.m1grid)then
                  jg=jdropgrid(ii,jj,k)
                  if(jg.gt.0)then
 c                  if(nseg(jg).gt.50)then
                    x=xptl(n)
                    y=yptl(n)
                    vrad=( x*pptl(1,n)/pptl(4,n)+y*pptl(2,n)/pptl(4,n))
                     if(vrad.lt.0.)then
                      ihit=1
                      goto 94
                     endif
 c                  endif
                  endif
                 endif
                endif
              enddo
             enddo
           endif
         endif
    94   continue
         if(ihit.eq.1)then
          delx=delxgr*(ii-i)
          dely=delxgr*(jj-j)
          xn=xptl(n)+delx
          yn=yptl(n)+dely
          ix=1+(xn+xgrid+shiftx)/delxgr
          jx=1+(yn+xgrid+shiftx)/delxgr
          jgrx=jdropgrid(ix,jx,k)
          if(jgrx.gt.0)then
           xptl(n)=xn
           yptl(n)=yn
           istptl(n) = 3
           nseg(jgrx)=nseg(jgrx)+1
           if(rmean(1,jg).gt.0d0)
      &         rmean(2,jgrx)=max(rmean(2,jgrx)
      &         ,abs(rmean(1,jgrx)-sqrt(xptl(n)**2+yptl(n)**2)))
          endif
         endif
   93  continue
       endif
 
 
 c Decay Strings not included in clusters taking into account Z
 c do it here to have this particles before the one from cluster 
 c (to define maxfra properly)
       if(ish.ge.6)write(ifch,*)'decay Strings not included in clusters'
       if(ifrade.ne.0.and.ntry.eq.1.and..not.lclean)then
         nptlx=nptl+1
 c copy decayed strings into new strings if no fragments are used for cluster
         do n=1,nptl
           if(istptl(n).eq.29.and.ityptl(n).lt.40)then   !fragmented central strings
             iused=0
             do nn=ifrptl(1,n),ifrptl(2,n)
               if(istptl(nn).eq.3.or.ifrptl(1,nn).ne.0)iused=iused+1
             enddo
             if(iused.eq.0)then   !if no fragment used, reset string
               istptl(n)=28
               do nn=ifrptl(1,n),ifrptl(2,n)     !suppress product
                 istptl(nn)=4
                 iaaptl(nn)=0
               enddo
               do nn=iorptl(n),jorptl(n)         !reinitialize string ends
                 istptl(nn)=20
               enddo              
             endif
           endif
         enddo
         call gakfra(0,iret)      !fragmentation using Z
         if(iret.gt.0)goto 1000
         do nn=nptlx,nptl  !exclude new particle from cluster formation
           iaaptl(nn)=0
           call jtain(nn,xptl(nn),yptl(nn),zptl(nn),tptl(nn),nnn,0)
           call jtaus(zptl(nn),tz,sptl(nn))
           strap=1e10
           xpl=tptl(nn)+zptl(nn)
           xmi=tptl(nn)-zptl(nn)
           if(xmi.gt.0.0.and.xpl.gt.0.0)strap=0.5*log(xpl/xmi)
           dezptl(nn)=strap       !space-time-rapidity
         enddo
         maxfra=nptl
         nptla=nptl
         if(ish.ge.6.and.nptl.ne.nptlx-1)
      &      call alist('list after second fragmentation&',nptlx,nptl)
         if(irescl.eq.1)then
           call utghost(iret)
           if(iret.gt.0)goto 1000
         endif
       endif
 
 
 c      if(iLHC.eq.1)then
 cc...split high pt segment into one part going out of the cluster and one part included in cluster (number of particle in cluster conserved)
 c        if(ish.ge.6)
 c     &  write(ifch,*)'mark high pt segments going into in clusters'
 c        do 100 n=1,nptla
 c          if(iaaptl(n).ge.0)goto 100
 c          i=1+(xptl(n)+xgrid+shiftx)/delxgr
 c          j=1+(yptl(n)+xgrid+shiftx)/delxgr
 c          k=1+(sptl(n)+sgrid+shifts)/delsgr
 c          if(    i.ge.1.and.i.le.m1grid
 c     &          .and.j.ge.1.and.j.le.m1grid
 c     &          .and.k.ge.1.and.k.le.m3grid)then
 c            jj=jdropgrid(i,j,k)
 c            if(jj.gt.0)then
 c              if(rmean(1,jj).gt.0d0)then
 c          ! path through cluster
 c                rp=sngl(rmean(2,jj)
 c     &               -abs(sqrt(xptl(n)**2+yptl(n)**2)-rmean(1,jj)))
 c          ! energy density (cluster mass approximate as sqrt(E))
 c                rh=sngl(sqrt(p4mean(4,jj)*nsegp4(jj))
 c     &               /(4d0/3d0*pi*rmean(2,jj)**3))
 c          ! Et used as Energy (transverse energy)
 c                pt2=pptl(1,n)**2+pptl(2,n)**2
 cc                ept=sqrt(pt2+pptl(3,n)**2)
 c                ept=sqrt(pt2+pptl(5,n)**2)
 c          ! energy loss based on BDMPS (peigne et al.)
 c                eloss=fploss*rh**(3./8.)*max(1.,sqrt(ept))*rp
 c          ! scaling factor
 c                fscal=eloss/ept
 cc                print *,ept,eloss,nsegp4(jj),rmean(2,jj),rh,rp,fscal
 c                
 c                if(fscal.le.0.)then   !exclude fragment from cluster
 c                  nseg(jj)=nseg(jj)-1
 c                  istptl(n)=0
 c                  iaaptl(n)=0
 c          ! update cluster mean energy
 c                  nsegp4(jj)=nsegp4(jj)-1
 c                  if(nsegp4(jj).ge.nsegce/iflhc)then
 c                    p4mean(4,jj)=(p4mean(4,jj)*dble(nsegp4(jj)+1)
 c     &                         -pptl(4,n))/dble(nsegp4(jj))
 c                  else
 c                    p4mean(4,jj)=1d50
 c                    rmean(1,jj)=-1d0
 c                  endif
 c                elseif(fscal.lt.1..and.(1.-fscal)**2*pt2.gt.
 c     &    (0.5-log(rangen()))**2*(ptclu**2))then
 cc     &    (0.5-log(rangen()))**2*(pptl(5,n)**2+ptclu**2))then
 cc                  print *,ept,eloss,nsegp4(jj),rmean(2,jj),rh,rp,fscal
 c                  if(ish.ge.6)write(ifch,*)'Quenched particle : ',n
 c     & ,idptl(n),ityptl(n),(pptl(kk,n),kk=1,5),ept,eloss,nsegp4(jj),jj
 c                  iaaptl(n)=0   !mother particle leave cluster
 c                  istptl(n)=0
 c                  nptl=nptl+1   !create fake particle with energy lost in cluster
 c                  nptla=nptla+1
 c                  istptl(nptl)=3 !daughter part of cluster
 c                  iaaptl(nptl)=1
 c                  do l=1,3
 c                    pptl(l,nptl)=fscal*pptl(l,n)
 c                    pptl(l,n)=(1.-fscal)*pptl(l,n)
 c                    xorptl(l,nptl)=xorptl(l,n)
 c                  enddo
 c! mother
 c                  tivptl(1,n)=xorptl(4,n)
 c                  pptl(4,n)=sqrt(pptl(1,n)**2+pptl(2,n)**2
 c     &                        +pptl(3,n)**2+pptl(5,n)**2)
 c                  call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm)
 c                  tivptl(2,n)=tivptl(1,n)+taugm*(-alog(rangen()))
 c          ! daughter
 c                  pptl(5,nptl)=0.
 c                  pptl(4,nptl)=sqrt(pptl(1,nptl)**2+pptl(2,nptl)**2
 c     &                           +pptl(3,nptl)**2+pptl(5,nptl)**2)
 c                  xorptl(4,nptl)=xorptl(4,n)
 c                  tivptl(1,nptl)=tivptl(1,n)
 c                  iorptl(nptl)=n
 c                  idptl(nptl)=idptl(n)*100+sign(99,idptl(n))
 c                  call idquac(nptl,idum1,idum2,idum3,jc)
 c                  ityptl(nptl)=ityptl(n)+2
 c                  xptl(nptl)=xptl(n)
 c                  yptl(nptl)=yptl(n)
 c                  zptl(nptl)=zptl(n)
 c                  sptl(nptl)=sptl(n)
 c                  dezptl(nptl)=dezptl(n)
 c                if(ish.ge.6)write(ifch,*)'--> Virtual particle : ',nptl
 c     & ,idptl(nptl),ityptl(nptl),(pptl(kk,nptl),kk=1,5)
 c          ! update cluster mean energy
 c                  p4mean(4,jj)=(p4mean(4,jj)*dble(nsegp4(jj))
 c     &                         -pptl(4,n)+pptl(4,nptl))/dble(nsegp4(jj))
 c                else
 c                  istptl(n)=3
 c                  iaaptl(n)=1
 c                  if(ish.ge.6)write(ifch,*)'Absorbed particle : ',n
 c     & ,idptl(n),ityptl(n),(pptl(kk,n),kk=1,5),ept,eloss,nsegp4(jj),jj
 c                endif
 c              endif
 c            endif
 c          endif
 c 100    continue
 c
 c
 c
 c      endif
 c
 
 c      nptlb0=nptl
 
       nptlb=nptl+jjj
 c...prepare /cptl/ for clusters
 
       if(ish.ge.6)write(ifch,*)'prepare /cptl/ for clusters'
       do jj=1,jjj
          nptl=nptl+1
           istptl(nptl)=12
           do l=1,4
             pptl(l,nptl)=0.
             xorptl(l,nptl)=0
           enddo
           sptl(nptl)=0
           uptl(nptl)=0
           optl(nptl)=0
           desptl(nptl)=0
           iorptl(nptl)=0
           jorptl(nptl)=0
 c limit the maximum number of subcluster to half the number of particle
 c (not to have empty subclusters)
           nsegmx(jj)=max(1,nint(float(nseg(jj))/float(nsegsuj)))
           if(ish.ge.6)write(ifch,*)'nsegmx',jj,nseg(jj),nsegmx(jj)
      &                                                 ,nsegsuj
       enddo
 
 c...prepare /cptl/ for subclusters
 
       if(ish.ge.6)write(ifch,*)'prepare /cptl/ for subclusters'
       mm=0
       do jj=1,jjj
         do ii=1,nsegmx(jj)
           mm=mm+1
           if(mm.gt.mxcl)stop'jintpo: mxcl too small.        '
           mmji(jj,ii)=mm
           mseg(mm)=0
           nptl=nptl+1
           istptl(nptl)=10
           do l=1,4
             pptld(l,mm)=0d0
             pptl(l,nptl)=0.
             xorptl(l,nptl)=0
           enddo
           sptl(nptl)=0
           uptl(nptl)=0
           optl(nptl)=0
           desptl(nptl)=0
           do l=1,nflav
             jccl(mm,l,1)=0
             jccl(mm,l,2)=0
           enddo
           iorptl(nptl)=nptla+jj
           jorptl(nptl)=0
           if(ii.eq.1)ifrptl(1,nptla+jj)=nptlb+mm
           ifrptl(2,nptla+jj)=nptlb+mm
         enddo
       enddo
 
 c...separate string segments, add dense area segments to clusters
 
       if(ish.ge.6)write(ifch,*)'separate string segments'
       do 98 n=1,nptla
         if(istptl(n).ne.3)goto 98
         i=1+(xptl(n)+xgrid+shiftx)/delxgr
         j=1+(yptl(n)+xgrid+shiftx)/delxgr
         k=1+(sptl(n)+sgrid+shifts)/delsgr
         if(    i.ge.1.and.i.le.m1grid
      &          .and.j.ge.1.and.j.le.m1grid
      &          .and.k.ge.1.and.k.le.m3grid)then
           jj=jdropgrid(i,j,k)
           if(jj.gt.0)then
             iimx=nsegmx(jj)
             do ii=1,iimx
               mm=mmji(jj,ii)
               if(mseg(mm).eq.0)goto 10          !not to have an empty cluster
             enddo
             ii=1+rangen()*iimx
             ii=min(ii,iimx)
             iini=ii
             am2tmpmx=1e30
  9          ntmp=mmji(jj,ii)
             am2tmp=(pptld(4,ntmp)+pptld(3,ntmp))
      &            *(pptld(4,ntmp)-pptld(3,ntmp))
             if(am2tmp.gt.xmxms**2/2.)then
               if(am2tmp.lt.am2tmpmx)then
                 mm=mmji(jj,ii)
                 am2tmpmx=am2tmp
               endif
               ii=ii+1
               if(ii.gt.iimx)ii=1
               if(ii.ne.iini)then
                 goto 9
               else
                 goto 10
               endif
             endif
             mm=mmji(jj,ii)
  10         continue
             mseg(mm)=mseg(mm)+1
             ifrptl(1,n)=mm      !local use of ifrptl
 c           write(ifch,*)'end',n
 c     & ,pptl(4,n),pptl(3,n),idptl(n),i,j,k,sptl(n),mm
             p4tmp=0d0
             do l=1,3
              pptld(l,mm)=  pptld(l,mm)  + dble(pptl(l,n))
              p4tmp=p4tmp+dble(pptl(l,n))*dble(pptl(l,n))
             enddo
             p4tmp=sqrt(p4tmp+dble(pptl(5,n)*pptl(5,n)))
             pptld(4,mm)=  pptld(4,mm)  + p4tmp
 c           if(mm.eq.86)write(ifch,*)'other',n
 c     & ,pptl(4,n),pptl(3,n),pptl(5,n),idptl(n),k,p4tmp
 c     & ,pptld(4,mm),pptld(3,mm),pptld(2,mm),pptld(1,mm)
 c     & ,(pptld(4,mm)+pptld(3,mm))*(pptld(4,mm)-pptld(3,mm))
             xorptl(1,nptlb+mm)=xorptl(1,nptlb+mm)+xptl(n)
             xorptl(2,nptlb+mm)=xorptl(2,nptlb+mm)+yptl(n)
             xorptl(3,nptlb+mm)=xorptl(3,nptlb+mm)+zptl(n)
             xorptl(4,nptlb+mm)=xorptl(4,nptlb+mm)+tptl(n)
             sptl(nptlb+mm)=sptl(nptlb+mm)+sptl(n)
             aa=cos(phievt)
             bb=sin(phievt)
             cc=-sin(phievt)
             dd=cos(phievt)
             xrot=xptl(n)*aa+yptl(n)*bb
             yrot=xptl(n)*cc+yptl(n)*dd
             uptl(nptlb+mm)=uptl(nptlb+mm)+xrot**2
             optl(nptlb+mm)=optl(nptlb+mm)+yrot**2
             desptl(nptlb+mm)=desptl(nptlb+mm)+xrot*yrot
             call idquac(n,idum1,idum2,idum3,jc)
 c            id=idptl(n)
 c            ida=iabs(id/10)
 c            ids=id/iabs(id)
 c            if(ida.ne.111.and.ida.ne.222.and.ida.ne.333)id=id/10*10
 c            if(ida.eq.111.or. ida.eq.222.or. ida.eq.333)id=id/10*10+ids
 c            if(ida.eq.213)id=1230*ids
 c            ic(1)=idtrai(1,id,1)
 c            ic(2)=idtrai(2,id,1)
 c            call iddeco(ic,jc)
             do l=1,nflav
               jccl(mm,l,1)=jccl(mm,l,1)+jc(l,1)
               jccl(mm,l,2)=jccl(mm,l,2)+jc(l,2)
             enddo
           else
             idropgrid(i,j,k)=0
           endif
         endif
   98  continue
 
 c...associate segments to clusters
 
       if(ish.ge.6)write(ifch,*)'associate segments to clusters'
       naseg(0)=0
       do jj=1,jjj
         do ii=1,nsegmx(jj)
           mm=mmji(jj,ii)
           naseg(mm)=naseg(mm-1)+mseg(mm)
           nfseg(mm)=0
         enddo
       enddo
       do 97 n=1,nptla
         if(istptl(n).ne.3)goto 97
         istptl(n) = 2
 c        write(ifch,*)'final segments ',n
 c     &          ,ityptl(n),istptl(n),idptl(n),dezptl(n),pptl(3,n)
         mm=ifrptl(1,n)
         nfseg(mm)=nfseg(mm)+1
         nsegmt(naseg(mm-1)+nfseg(mm))=n
  97   continue
       do jj=1,jjj
         nst=0
         do ii=1,nsegmx(jj)
           mm=mmji(jj,ii)
           if(mseg(mm).ne.nfseg(mm))stop'jintpo: mseg.ne.nfseg        '
           nst=nst+mseg(mm)
         enddo
         if(nst.ne.nseg(jj))stop'sum(mseg(mm)).ne.nseg(jj)'
       enddo
 
 c...finish cluster storage to /cptl/
 
       if(ish.ge.6)write(ifch,*)'finish cluster storage to /cptl/'
       xx=0.
       yy=0.
       xy=0.
       mjjsegsum=0
       do jj=1,jjj
        njj=nptla+jj
        mjjseg=0
        do l=1,nflav
          jcjj(l,1)=0
          jcjj(l,2)=0
        enddo
        do ii=1,nsegmx(jj)
         mm=mmji(jj,ii)
         n=nptlb+mm
 
         do l=1,nflav
           jc(l,1)=jccl(mm,l,1)
           jc(l,2)=jccl(mm,l,2)
           ke(l)=jc(l,1)-jc(l,2)
           jcjj(l,1)=jcjj(l,1)+jc(l,1)
           jcjj(l,2)=jcjj(l,2)+jc(l,2)
         enddo
         call idenct(jc,idptl(n)
      *  ,ibptl(1,n),ibptl(2,n),ibptl(3,n),ibptl(4,n))
         ttest=0d0
         do ji=1,4
          ptest(ji)=0d0
           do ns=naseg(mm-1)+1,naseg(mm)
             ni=nsegmt(ns)
             ptest(ji)=ptest(ji)+dble(pptl(ji,ni))
           enddo
          ptest(ji)=abs(ptest(ji)-pptld(ji,mm))
          ttest=ttest+ptest(ji)
         enddo
         amcmin=1.01*utamnu(ke(1),ke(2),ke(3),ke(4),ke(5),ke(6),4)
         p52=((pptld(4,mm)+pptld(3,mm))*(pptld(4,mm)-pptld(3,mm))
      &      -pptld(1,mm)**2-pptld(2,mm)**2)
         if(iLHC.eq.1)then
           amcmin=sqrt(max(amcmin**2,sngl(p52)))
           if(amcmin/ycor.gt.amuseg)amcmin=amcmin*ycor
 c     &                    **max(0.,min(1.,pptld(4,mm)-2.*amcmin*ycor))
         endif
           
         pptld(5,mm)=0d0
         if(p52.gt.0d0)then
          pptld(5,mm)=sqrt(p52)
          jerr(2)=jerr(2)+1
          if(iLHC.eq.1.and.pptld(5,mm).lt.amcmin
      &               .and.pptld(4,mm).gt.amcmin)then
            pptld(5,mm)=dble(amcmin)
            bp=sqrt((pptld(4,mm)+pptld(5,mm))*(pptld(4,mm)-pptld(5,mm))
      &            /(pptld(3,mm)*pptld(3,mm)+pptld(2,mm)*pptld(2,mm)
      &                                      +pptld(1,mm)*pptld(1,mm)))
            pptld(1,mm)=bp*pptld(1,mm)
            pptld(2,mm)=bp*pptld(2,mm)
            pptld(3,mm)=bp*pptld(3,mm)
 c           write(ifch,*)"ici ",n,sqrt(p52),pptld(5,mm),bp,pptld(4,mm),
 c     &  sqrt(pptld(3,mm)*pptld(3,mm)
 c     &                   +pptld(2,mm)*pptld(2,mm)
 c     &                   +pptld(1,mm)*pptld(1,mm)
 c     &                   +pptld(5,mm)*pptld(5,mm))
 c      write(ifch,*)'droplet uds=',ke(1),ke(2),ke(3),'   E=',pptld(5,mm)
          endif
         elseif(p52.le.0d0)then
          jerr(3)=jerr(3)+1
          pptld(5,mm)=dble(amcmin)
          pptld(4,mm)=sqrt(pptld(3,mm)*pptld(3,mm)
      &                   +pptld(2,mm)*pptld(2,mm)
      &                   +pptld(1,mm)*pptld(1,mm)
      &                   +pptld(5,mm)*pptld(5,mm))
         endif
         if(ish.ge.1.and.(abs(ttest).gt.1.d0.or.pptld(5,mm).gt.xmxms))
      &    then
           call utmsg('jintpo&')
           write(ifmt,*)'***** Warning in jintpo !',ntry
           write(ifch,*)'***** jintpo: momenta messed up (ttest > 0)'
           write(ifch,*)'*****',mm,n,mseg(mm),p52,ttest
           write(ifch,*)'*****',jj,nsegp4(jj),p4mean(4,jj)
           write(ifch,'(a,10x,4f15.4)')'*****',(pptld(ji,mm),ji=1,4)
           do ns=naseg(mm-1)+1,naseg(mm)
             ni=nsegmt(ns)
             write(ifch,'(a,i5,i9,5f15.4,f12.4)')'*****',ni,idptl(ni)
      *     ,(pptl(ji,ni),ji=1,4),pptl(5,ni)**2
      *     ,(pptl(4,ni)+pptl(3,ni))*(pptl(4,ni)-pptl(3,ni))
      *       -pptl(1,ni)**2-pptl(2,ni)**2
           enddo
         endif
         if(pptld(5,mm).gt.xmxms)then
           p4max=0.
           nh=0
           do ns=naseg(mm-1)+1,naseg(mm)
             ni=nsegmt(ns)
             if(pptl(4,ni).ge.p4max)then
               nh=ni
               p4max=pptl(4,ni)
             endif
           enddo
           if(nh.le.0)then
             stop'Cannot be in jintpo ...'
           else   !put back nh as normal particle
             iaaptl(nh) = 0      !don't use this fragment next time
 c        if(iaaptl(nh).eq.0)print*,'excluded',idptl(nh),ityptl(nh)
 c     &                                         ,pptl(4,nh)
             if(mod(abs(idptl(nh)),100).eq.99)then !restore lost energy
               istptl(nh) = 5
               ior=iorptl(nh)
               iaaptl(ior)=0      !don't use this fragment next time
               ifrptl(1,ior) = 0
               ifrptl(2,ior) = 0
               if(istptl(ior).eq.0)then                   
                 do k=1,3
                   pptl(k,ior)=pptl(k,ior)+pptl(k,n)
                 enddo
                 pptl(4,ior)=sqrt(pptl(1,ior)**2+pptl(2,ior)**2
      &                          +pptl(3,ior)**2+pptl(5,ior)**2)
                 call idtau(idptl(ior),pptl(4,ior),pptl(5,ior),taugm)
                 tivptl(2,ior)=tivptl(1,ior)+taugm*(-alog(rangen()))
               else
                 istptl(ior) = 0
               endif
             elseif(idptl(nh).lt.1e4)then
               istptl(nh) = 0     !particle
               ifrptl(1,nh) = 0
               ifrptl(2,nh) = 0
             else
               istptl(nh) = 10    !droplet
             endif
           endif
           if(ish.ge.1)
      &    write(ifch,*)'***** Redo cluster without heavy particle :'
      &                 ,nh,ntry
           goto 8888
         endif
         do l=1,5
          pptl(l,n)=sngl(pptld(l,mm))
         enddo
         mjjseg=mjjseg+mseg(mm)
         do l=1,4
           pptl(l,njj)=pptl(l,njj)+pptl(l,n)
           xorptl(l,njj)=xorptl(l,njj)+xorptl(l,n)
           xorptl(l,n)=xorptl(l,n)/float(mseg(mm))
         enddo
         sptl(njj)=sptl(njj)+sptl(n)
         uptl(njj)=uptl(njj)+uptl(n)
         optl(njj)=optl(njj)+optl(n)
         desptl(njj)=desptl(njj)+desptl(n)
         sptl(n)=sptl(n)/float(mseg(mm))
         uptl(n)=uptl(n)/float(mseg(mm))
         optl(n)=optl(n)/float(mseg(mm))
         desptl(n)=desptl(n)/float(mseg(mm))
         radptl(n)=0
         istptl(n)=10
         ifrptl(1,n)=0
         ifrptl(2,n)=0
         tivptl(1,n)=xorptl(4,n)
         tivptl(2,n)=ainfin
         ityptl(n)=60
         radptl(n)=0
         dezptl(n)=0.
        enddo !ii
        do l=1,4
         xorptl(l,njj)=xorptl(l,njj)/float(mjjseg)
        enddo
        mjjsegsum=mjjsegsum+mjjseg
        xx=xx+uptl(njj)
        yy=yy+optl(njj)
        xy=xy+desptl(njj)
        sptl(njj)=sptl(njj)/float(mjjseg)
        uptl(njj)=uptl(njj)/float(mjjseg)
        optl(njj)=optl(njj)/float(mjjseg)
        desptl(njj)=desptl(njj)/float(mjjseg)
        pjj52=(pptl(4,njj)+pptl(3,njj))*(pptl(4,njj)-pptl(3,njj))
      &    -pptl(1,njj)**2-pptl(2,njj)**2
         pptl(5,njj)=0
         if(pjj52.gt.0)then
           pptl(5,njj)=sqrt(pjj52)
         endif
         ityptl(njj)=60
         call idenct(jc,idptl(njj)
      *  ,ibptl(1,njj),ibptl(2,njj),ibptl(3,njj),ibptl(4,njj))
       enddo !jj
 
 
 c...ranphi
 
       ranphi=0
       rini=0.
       if(mjjsegsum.gt.0)then
         xx=xx/float(mjjsegsum)
         yy=yy/float(mjjsegsum)
         xy=xy/float(mjjsegsum)
         dta=0.5*(xx-yy)
 c        eba=0.5*(xx+yy)
 c        ww=-xy
 c        !-------------------------------------------------------
 c        !    inertia tensor:
 c        !----------------------+-------------------------------+
 c        !  <y**2>   -<x*y>     !      with <x**2>=uptl         !
 c        !  -<x*y>   <x**2>     !     <y**2>=optl  <xy>=desptl  !
 c        !----------------------+-------------------------------+
 c        ! Eigenvalues ev1, ev2
 c        !-------------------------------------------------------
 c        ev1=eba+sqrt(dta**2+ww**2)
 c        ev2=eba-sqrt(dta**2+ww**2)
         if(xy.lt.0..and.dta.ne.0.)then
           ranphi=0.5*atan(-xy/dta)
         elseif(xy.gt.0..and.dta.ne.0.)then
           ranphi=-0.5*atan(xy/dta)
         else
           ranphi=0
         endif
 c        if(dta.ne.0.)then
 c         ranphi=0.5*atan(abs(ww)/abs(dta))
 c         if(    ww.gt.0..and.dta.gt.0.)then
 c          ranphi=ranphi
 c         elseif(ww.lt.0..and.dta.gt.0.)then
 c          ranphi=-ranphi
 c         elseif(ww.gt.0..and.dta.lt.0.)then
 c          ranphi=pi-ranphi
 c         elseif(ww.lt.0..and.dta.lt.0.)then
 c          ranphi=ranphi-pi
 c         endif
 c        else
 c          ranphi=2.*pi*rangen()
 c        endif
         rini=max(0.01,sqrt(5./3.*(xx+yy))) !<r**2>=3/5*R**2 for sphere of radius R
         volu=(4./3.*pi*(xx+yy)**1.5)
 c        rho=amctot/volu
         flowpp=0.
         flowaa=0.
 
         if(iLHC.eq.1.and.visco.gt.0.)then
           if(npair.gt.0)then
             fcorr=min(1.,float(mapmax)*abs(fvisco)/float(maptot))
             visco=min(1.,2.*float(maptot)/float(npair)/float(mapmax))**2
 c     &                 *abs(fvisco))
           elseif(lclean)then  !large number of particles, npair can't be calculated
             visco=1e-6
             fcorr=1.
           else          !cluster from remnants only
             visco=1.
             fcorr=1.
           endif
           if(visco.ge.1.)yrmaxi=0. !yrmaxi*(1.-visco)
 c          visco=exp(-min(50.,max(0.,
 c     &         float(koievt)/log(amctot)-abs(fvisco))**yradpi)) !mix flow 
 c         &               max(0.,rmax**2-abs(fvisco))))    !mix flow 
 c          yrmaxi=log(amctot**2)
 c          yrmaxi=yradmx*yrmaxi*(1.-visco)
 c          if(visco.lt.1.and.yrmaxi.gt.1e-2)then
 c            yyrmax=dble(yrmaxi)
 c            fradflii=sngl(1d0/
 c     &         ((sinh(yyrmax)*yyrmax-cosh(yyrmax)+1d0)/(yyrmax**2/2d0)))
 c          else
 c            visco=1.
 c            yrmaxi=0.
 c            fradflii=1.
 c          endif
 c          if(visco.gt.1e-5)then
 c            yrmaxi=yradmx*(yrmaxi
 c     &            +visco*(log(amctot)*yradpx/yradmx-yrmaxi))
 c          else
 c            yrmaxi=yradmx*yrmaxi
 c          endif
           fradflii=1.
           if(yrmaxi.gt.0)then
             flowpp=visco*log(fcorr*amctot)*yradpx
             flowaa=yrmaxi
             if(rangen().lt.flowaa/flowpp)then
               visco=0.
               yrmaxi=flowaa+max(0.,flowpp-flowaa)
               yyrmax=dble(yrmaxi)
               fradflii=sngl(1d0/
      &        ((sinh(yyrmax)*yyrmax-cosh(yyrmax)+1d0)/(yyrmax**2/2d0)))
             else
               yrmaxi=0.
               visco=log(fcorr*amctot)/log(amctot)
             endif
           elseif(fcorr.lt.1.)then
             visco=log(fcorr*amctot)/log(amctot)
           endif
         endif
 
         if(ish.ge.3)write(ifch,*)'yrmaxi,delzet=',yrmaxi,delzet
 c        print *,'->',bimevt,yrmaxi,visco*yradpx*log(amctot),flowaa
 c     &        ,flowpp,maptot,log(ectot),visco,log(amctot),rho
 c     &        ,mapmax,npair
 c     &        ,min(1.,2.*float(maptot)/float(npair)/float(mapmax))**2
       endif
 
 c...print
 
       if(ish.ge.5)then
         write(ifch,*)'print'
         do k=1,m3grid
         write(ifch,*)'k=',k,'  jclu=',jclu(k)
      &              ,'  s=',(k-1)*delsgr-sgrid-shifts
          do j=m1grid,1,-1
         write(ifch,'(10i4,3x,10i4)')(idropgrid(i,j,k),i=1,m1grid)
      &                ,(jdropgrid(i,j,k),i=1,m1grid)
         enddo
         enddo
         write(ifch,'(a,a)')
      &    '    k   jj  nseg      mm  mseg     n      mass'
      &   ,'         s         y         z         t '
         do jj=1,jjj
          do ii=1,max(1,nint(1.*nseg(jj)/nsegsuj))
            mm=mmji(jj,ii)
            n=nptlb+mm
            sg=pptl(3,n)/abs(pptl(3,n))
            tm=sqrt(pptl(5,n)**2+pptl(1,n)**2+pptl(2,n)**2)
            y=sg*alog((pptl(4,n)+sg*pptl(3,n))/tm)
 c           if(kclu(jj).eq.44)print *,tm,pptl(4,n),pptl(3,n),iorptl(n)
            write(ifch,'(2i5,i6,i8,2i6,5f10.3)')
      &       kclu(jj),jj,nseg(jj),mm,mseg(mm),n,pptl(5,n)
      &      ,sptl(n),y,xorptl(3,n),xorptl(4,n)
          enddo
         enddo
       endif
 
 c...decay
       iret=0
       if(jjj.gt.0)then     !decay only if some cluster produced
       if(4-abs(typevt).gt.0.0001)typevt=-typevt    !typevt < 0 if fusion but only if not SD (sign used for something else for SD ... and no fusion produced for SD events)
       if(ish.ge.5)write(ifch,*)'decay ...'
       if(ifrade.eq.0.or.ispherio.gt.0)goto1000
       if(jdecay.eq.0)goto1000
       nptlbcf=nptl
       nptl0=nptl
       if(hydt.ne.'---')then
         call HydroFO(ier)
       else
         nclu=0
         ptest(1)=0d0
         ptest(2)=0d0
         ptest(3)=0d0
         ptest(4)=0d0
         ptest(5)=0d0
         do jj=1,jjj
           do ii=1,max(1,nint(1.*nseg(jj)/nsegsuj))
            mm=mmji(jj,ii)
            np=nptlb+mm
 c           print *,'decay',jj
            if(ioclude.eq.3)then
              call hnbaaa(np,iret)
            else
              call DropletDecay(np,iret)
            endif
            if(iret.eq.1)then
              istptl(np)=istptl(np)+2
              do ns=naseg(mm-1)+1,naseg(mm)
                n=nsegmt(ns)
                if(mod(abs(idptl(n)),100).eq.99)then   !restore lost energy
                  istptl(n) = 5
                  ior=iorptl(n)
                  ifrptl(1,ior) = 0
                  ifrptl(2,ior) = 0
                  if(istptl(ior).eq.0)then                   
                    do k=1,3
                      pptl(k,ior)=pptl(k,ior)+pptl(k,n)
                    enddo
                    pptl(4,ior)=sqrt(pptl(1,ior)**2+pptl(2,ior)**2
      &                             +pptl(3,ior)**2+pptl(5,ior)**2)
                    call idtau(idptl(ior),pptl(4,ior),pptl(5,ior),taugm)
                    tivptl(2,ior)=tivptl(1,ior)+taugm*(-alog(rangen()))
                  else
                    istptl(ior) = 0
                  endif
                elseif(idptl(n).lt.1e4)then
                  istptl(n) = 0     !particle
                  ifrptl(1,n) = 0
                  ifrptl(2,n) = 0
                else
                  istptl(n) = 10    !droplet
                endif
              enddo
            elseif(ioclude.eq.3)then
              do i=nptl0+1,nptl
                if(ityptl(i).eq.60)then
                  nclu=nclu+1
                  ptest(1)=ptest(1)+dble(pptl(1,i))
                  ptest(2)=ptest(2)+dble(pptl(2,i))
                  ptest(3)=ptest(3)+dble(pptl(3,i))
                  ptest(4)=ptest(4)+dble(pptl(4,i))
                endif
              enddo
              nptl0=nptl
            endif
           enddo
         enddo
       endif
       do jj=1,jjj
        do ii=1,max(1,nint(1.*nseg(jj)/nsegsuj))
         mm=mmji(jj,ii)
         np=nptlb+mm
         istptl(np)=istptl(np)+1
         ifrptl(1,np)=nptlbcf+1
         ifrptl(2,np)=nptl
         rinptl(np)=kclu(jj)-m3grid/2
        enddo
       enddo
 c add global flow on all particles of all decayed cluster
       if(iLHC.eq.1)then
         yrmax=0.
         if(fvisco.gt.0.)yrmax=yradpx*visco
       else
         yrmax=yradpx
       endif
 c      print *,bimevt,rini,yrmax,yrmaxi,delzet,np,ptmax,visco
       if(ioclude.eq.3.and.yrmax.gt.1e-3.and.nclu.gt.0)then
 c set angular informations
         fecc=0
         aa=1
         bb=0
         cc=0
         dd=1
         dta=0.5*abs(xx-yy)
         ev1=(xx+yy)/2+sqrt(dta**2+xy**2)
         ev2=(xx+yy)/2-sqrt(dta**2+xy**2)
         ecc=(ev1-ev2)/(ev1+ev2)
 c        fecc=facecc*ecc!/(1.+yrmax)
 c        print*,'pp',ecc,ranphi
         fecc=min(facecc,ecc)   !be careful : fecc change <pt> since it is the elliptical deformation of the sub cluster (give strength of v2)
 
         phiclu=mod(phievt+ranphi,2.*pi) !do not change otherwise v2 is gone
         if(phiclu.lt.-pi)phiclu=phiclu+2*pi
         if(phiclu.gt.pi)phiclu=phiclu-2*pi
         aa=cos(phiclu)
         bb=sin(phiclu)
         cc=-sin(phiclu)
         dd=cos(phiclu)
         errlim=0.00005
 c loop on particles for each main cluster
 c        ncl=nptlb0
         npass=max(1,min(nclu/5,jjj)) !to have the same number of group of particles than original clusters but different repartition of particles
         npart=nclu/npass
         if(npart*npass-nclu.gt.max(5,npart/2))npass=npass+1
 c        print *,'ici',nclu,npass,npart
 c        lcont=.true.
 c        if(nclu.lt.50)lcont=.false.
         lcont=.false.
         ncl=0
         nmin=nptlbcf+1
         nmax=nptl
         idrc=-1+2.*int(0.5+rangen())
         ntot=nclu
 
 c       prepare debug output for flow
         nall=0
         if(ish.ge.5)then
           nall=nmax-nmin+1
           do ii=1,nall
             idptl(nptl+ii)=idptl(nptlbcf+ii)
           enddo
           iorptl(nptl+nall+1)=nptl+1
           jorptl(nptl+nall+1)=nptl+nall
           do k=1,5
             pptl(k,nptl+nall+1)=0
             do ii=1,nall
               pptl(k,nptl+ii)=pptl(k,nptlbcf+ii)
               pptl(k,nptl+nall+1)=pptl(k,nptl+nall+1)+pptl(k,nptlbcf+ii)
             enddo
           enddo
         endif
         
 c initialization for rescaling on ipart particles    
         ipart=0
         tecm0=0.
         tecm=0.
         ini0=0
         ifi0=0
 
 c        do 900 while (ncl.le.nptlb-1)
         do while (ntot.gt.0)
           ncl=ncl+1
           if(iLHC.eq.1)then
             yrmax=yradpx*visco
           else
             yrmax=yradpx
           endif
 c         cms frame of all particles from same cluster
           do k=1,5
             ppp(k)=0d0
           enddo
           ini=nptl
           ifi=1
           if(idrc.gt.0)then
             imax=nmax
             imin=nmin-1
             if(ipart.eq.0)ini0=nmin
           else
             imax=nmin
             imin=nmax+1
             if(ipart.eq.0)ifi0=nmax
           endif
 c          npclu=0
 c 880      ncl=ncl+1
           n=0
           i=imin
           lpass=.true.
           do while ((n.lt.npart.or.ncl.eq.npass)
      &              .and.idrc*i.lt.idrc*imax.and.lpass)
             i=i+idrc
 c            if(jorptl(i).eq.ncl)then
 c            if(jorptl(i).eq.ncl.and.ityptl(i).eq.60)then
             
             if(ityptl(i).eq.60)then
               n=n+1
               ntot=ntot-1
               ini=min(ini,i)
               ifi=max(ifi,i)
 c              if(ityptl(i).eq.60)npclu=npclu+1
               do k=1,4
                 ppp(k)=ppp(k)+dble(pptl(k,i))
               enddo
             elseif(.not.lcont.and.n.gt.0)then
               lpass=.false.
             endif
 c            print *,ityptl(i),i,imin,imax,idrc,n,npart,ncl,nclu,nptl
           enddo
           np=n
           if(idrc.gt.0)then
             nmin=nmin+i-imin
             ifi0=i
           else
             nmax=nmax+i-imin
             ini0=i
           endif
           
 c          if(ncl.lt.nptlb.and.npclu.lt.int(0.2*(nptl-nptlbcf+1)))goto880
           if(ifi.le.ini)goto 900
 
 c record info for rescaling
           ipart=ipart+np
 
 c test mass
           ppp(5)=(ppp(4)-ppp(3))*(ppp(4)+ppp(3))-ppp(2)**2-ppp(1)**2
           if(ppp(5).gt.0d0)then
             ppp(5)=sqrt(ppp(5))
           else
             if(ish.ge.1)write(ifch,*)'Precision problem in jintpo, p:',
      &           (ppp(k),k=1,5)
             ppp(5)=0d0
           endif
             if(ish.ge.4)
      &           write(ifch,*)'Group of particle: ',
      &                        idrc,ini,ifi,ncl,'/',npass,npart,nclu
           if(ppp(5).gt.0d0)then    ! here all particle should have flow
             do i=ini,ifi
              if(ityptl(i).eq.60)then
               tecm=tecm+pptl(4,i)  !use energy in collision frame as reference
               call utlob4(1,ppp(1),ppp(2),ppp(3),ppp(4),ppp(5)
      $             ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
              endif
             enddo
             if(tecm.gt.0.)then
               yrmax=yrmax*log(amctot) !/rini**yradpi
             else
               yrmax=0.
             endif
 c           print *,bimevt,rini,ppp(5),yrmax,yrmaxi,delzet,np,ptmax,visco
             if(ish.ge.4)
      &      write(ifch,*)'Radial flow: ',yrmax,tecm,visco,yradpx,yradpp
             if(yrmax.gt.0.)then
               if(np.gt.maxp)stop'maxp too small in jintpo'
               i=0
               if(ish.ge.8)call clist('list before flow&',ini,ifi,60,60)
               do ii=ini,ifi
                if(ityptl(ii).eq.60)then
                 i=i+1
                 yrad(i)=sqrt(rangen())
                 phirad(i)=2.*pi*rangen()
                 pt2=(pptl(1,ii)**2+pptl(2,ii)**2)!+pptl(5,ii)**2)
                 bex=sinh(yrad(i)*yrmax)*cos(phirad(i))
      &             *(1+fecc/(1.+pt2))
                 bey=sinh(yrad(i)*yrmax)*sin(phirad(i))
      &             *(1-fecc/(1.+pt2))
                 be(1)=aa*bex+cc*bey
                 be(2)=bb*bex+dd*bey
                 be(3)=0d0
                 be(4)=sqrt(1+be(1)**2+be(2)**2)
 c                call utlob4(1,be(1),be(2),be(3),be(4),1d0
 c     *         , pptl(1,ii), pptl(2,ii), pptl(3,ii), pptl(4,ii))
 c mimic boost transformation but protect against to high values of p(3) (p(3)~p(4)) 
                 pt2=pptl(1,ii)**2+pptl(2,ii)**2
                 bet=-(pptl(1,ii)*be(1)+pptl(2,ii)*be(2))/(be(4)+1.)
                 
                 pt=yradpp**max(1.,pptl(5,ii))
                 fac=1./(1.+sqrt(pt2/pptl(5,ii)**2))**pt
                 bet=bet+sqrt(pt2+(pptl(3,ii)**2+pptl(5,ii)**2)*fac)
                 
                 pptl(1,ii)=pptl(1,ii)-bet*be(1)
                 pptl(2,ii)=pptl(2,ii)-bet*be(2)
                 pptl(4,ii)=sqrt(pptl(1,ii)**2+pptl(2,ii)**2
      *                   +pptl(3,ii)**2+pptl(5,ii)**2)
               else
                 yrad(i)=0.
                 phirad(i)=0.
                endif
               enddo
               if(ish.ge.8)call clist('list after flow&',ini,ifi,60,60)
 
 
 c           boost back
               pe(1)=0.
               pe(2)=0.
               pe(4)=0.
               do i=ini,ifi
                 if(ityptl(i).eq.60)then
                   call utlob4(-1,ppp(1),ppp(2),ppp(3),ppp(4),ppp(5)
      $             ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
                   pe(1)=pe(1)+pptl(1,i)
                   pe(2)=pe(2)+pptl(2,i)
                   pe(4)=pe(4)+pptl(4,i)
                 endif
               enddo
 
 
 
 cc random rescaling
 c rescaling has to be done for different bins in eta to conserve
 c energy flow (if everything is done at the end on all particles,
 c the energy flow is concentred at mid-rapidity in contradiction
 c with ATLAS data
 c On the other hand, a sufficient number of particles is necessary
 c to have proper eta distribution in particular at high pt.
 c Since the clusters are ordered in rapidity we can do the rescaling on
 c group of clusters having enough particles for the rescaling but
 c not too much to keep the eta dependence of the energy flow.
 
 c set ipart.ge.1 to do rescaling for each subclusters
 
               if((ipart.ge.1.and.ntot.ge.ipart/2).or.ntot.eq.0)then
 
                 ini=ini0
                 ifi=ifi0
 
                 if(ish.ge.8)
      &          call clist('list before flow rescaling&',ini,ifi,60,60)
 
                 if(fplmin.gt.0)then
 
 c                  ntrt=ntrt+1
 
                   niter=0
  611              energ=0.
                   energ0=0.
                   niter=niter+1
                   i=0
                   ptmax=0.
                   plmax=0.
                   do  j=ini,ifi
                     if(ityptl(j).eq.60)then
                       i=i+1
                       pt2=pptl(1,j)**2+pptl(2,j)**2
                       pz2=pptl(3,j)**2
                       pp2=pt2+pz2
 c                      et2=(pp2+pptl(5,j)**2)*pt2/pp2
                       pt=sqrt(pt2)
 c                      pp=sqrt(pp2)
 c base necessary to avoid peak at pt or pl=0
 c epsi change the shape of eta distributions and pt for
 c identified particles (shift the maximum of the distribution)
 c                      epsi=min(0.99,fplmin*rangen()**0.3)
                       base=0. !sqrt(1./(1.-epsi)**2-1.)*pptl(5,j)/pp
                       finc=2.
                       if(energ0.lt.tecm)then
                         finc=finc*sqrt(float(max(1,niter-300)))
 c                        base=base+min(tecm/pe(4),
 c     &                           (max(0,niter-300))*0.3*pptl(5,j)/pp)
                       else
 c                        base=base/log10(max(10.,float(niter-300)))
                         finc=finc/sqrt(max(1.,float(niter-300)))
                       endif
 c                      if(1.-pptl(5,j)/pptl(4,j).ge.epsi
 c     &                   .and.1.+pptl(5,j)**2*(1./pp2-1./pt2).gt.0.)then
                       if(1.+pptl(5,j)**2*(1./pp2-1./pt2).gt.0.)then
                         ptmax=max(pt,ptmax)
                         plmax=max(abs(pptl(3,j)),plmax)
                         yrad2(i)=rangen()
                         yrad2(i)=yrad2(i)*
      &                       max(0.,min(1.,finc*tecm/pe(4))-base)
 
 
 c necessary even with epsi cut to smooth eta distribution
 c and since sumEt has to be conserved, we can constrain scaling for pt
 c and pz :
 c                        ytmp=yrad2(i)
                         yrad2(i)= yrad2(i)**((1.-(pptl(5,j)
      &    /sqrt((min(1.,yrad2(i)*float(max(1,niter-100))))**2
      &         *(pt**2+pptl(3,j)**2)+pptl(5,j)**2))**1.)
      &    *exp(-fplmin*max(0.,pt**2-(pe(4)/tecm)**2)))
 
                         yrad2(i)=min(1.,base+yrad2(i)) !should be here to avoid peak at eta=0.
                       else
                         yrad2(i)=1.
 c                        ytmp=yrad2(i)
                       endif
                       be(1)= pptl(1,j)*yrad2(i)
                       be(2)= pptl(2,j)*yrad2(i)
                       be(3)= pptl(3,j)*yrad2(i)                
 
 c      print *,niter,i,ntry,yrad2(i),energ/tecm,pt,pptl(3,j),finc,base
 c     *       ,ytmp,(1.-(pptl(5,j)
 c     &    /sqrt((min(1.,yrad2(i)*float(max(1,niter-100))))**2
 c     &         *(pt**2+pptl(3,j)**2)+pptl(5,j)**2))**0.1)
 
                       energ=energ+sqrt(be(1)**2+be(2)**2
      *                       +be(3)**2+pptl(5,j)**2)
                     endif
                   enddo
                   energ0=energ
 c          print *,'fin',niter,energ/tecm
                   if(abs(energ-tecm)/tecm.gt.1..and.niter.lt.1000)then
                     goto 611
                   elseif(niter.ge.1000)then
 c      print *,'Rescaling failed:',pe(4),tecm,ptmax,plmax,ipart,ini0,ifi
                     if(ish.ge.2)write(ifch,*)'Random rescaling failed:'
      &                   ,energ,tecm,ptmax,plmax,ipart,ini0,ifi
 c                    ntrr=ntrr+1
                     goto 200
                   endif
 c      print *,'done',niter,energ/tecm,ipart,ini0,ifi,finc
                   if(ish.ge.5)write(ifch,*)'Rescaling done:'
      &               ,tecm,energ/tecm,niter,ptmax,plmax,ipart,ini0,ifi
                   i=0
                   do  j=ini,ifi
                     if(ityptl(j).eq.60)then
                       i=i+1
                       pptl(1,j)= pptl(1,j)*yrad2(i)
                       pptl(2,j)= pptl(2,j)*yrad2(i)
                       pptl(3,j)= pptl(3,j)*yrad2(i)                
                       pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2
      *                     +pptl(3,j)**2+pptl(5,j)**2)
                     endif
                   enddo
                 endif
 
  200            continue
 
 
               if(ish.ge.8)
      &        call clist('list after flow rescaling&',ini,ifi,60,60)
 
               tecm0=tecm0+tecm
               tecm=0.
               ipart=0
 
             endif
 
           endif                 ! yrmax
 
         endif     !p(5)>0
 
 
  900  continue
       enddo
 
 c rescale momentum precisely and globaly to avoid artefacts for matrix
 c but only at 10% level to keep dependence of sumEt with eta 
       if(tecm0.gt.0)then
         ini=nptlbcf+1
         ifi=nptl
         esoll=tecm0
         scal=1.
         do ipass=1,2000
           sum=0.
           n=0
           do  j=ini,ifi
             if(ityptl(j).eq.60)then
               n=n+1
 c             this part is EXTREMELY important for the pseudorapidity shape at various pt
 c             if nothing special a broad peak appear at eta=0
 c             to avoid that, scal has to be reduced when p3 or pt reach 0
               if(scal.lt.1.)then
                 scal0=scal
                 pt=sqrt(pptl(1,j)**2+pptl(2,j)**2)
                 if(fplmin.le.0.)then
                   pt=abs(fplmin)/sqrt(pptl(5,j))*pt
 c                  pt=abs(fplmin)*pt
                 else
                   pt=5.*pt
                 endif
                 pow=sqrt(pptl(5,j)**2+scal**2*(pt**2
      *                           +pptl(3,j)**2))
                 pow=(1.-(1./sqrt(pptl(5,j))+pptl(5,j))
      *                 /(1./sqrt(pptl(5,j))+pow))
                 pow=rangen()*pow
               else
                 pow=1.-2.*pptl(4,j)/engy  !to avoid particle with energy larger than beam energy
                 if(pow.lt.0.)then
                   scal0=1./((1.-pow)
      *                     *exp(-0.25*max(-4.,log(rangen()))))
                   pow=1.
 c                  print *,j,pptl(4,j),pow,scal0,scal
 c     *                              ,pptl(3,j)*scal0**pow
                 else
                   scal0=scal
                   pow=rangen()*pow
                 endif
               endif
 
               do k=1,3
                 pptl(k,j)=pptl(k,j)*scal0**pow !to smooth distributions   
               enddo
               pptl(4,j)=sqrt(pptl(1,j)**2+pptl(2,j)**2
      *             +pptl(3,j)**2+pptl(5,j)**2)
               sum=sum+pptl(4,j)
             endif
           enddo
           scal=esoll/sum
 c         write(ifmt,*)'ipass,scal,e,esoll:'
 c     $             ,ipass,scal,sum,esoll
           if(abs(scal-1.).le.errlim) goto 300
         enddo
  300    continue
 c         write(ifmt,*)'ipass,scal,e,esoll:'
 c     $             ,ipass,scal,sum,esoll
 
 c adjust pt to have pt conservation in cms of particles having flow
         if(nclu.gt.0)then
           ptest(5)=(ptest(4)-ptest(3))*(ptest(4)+ptest(3))-ptest(2)**2
      &            -ptest(1)**2
           if(ptest(5).gt.0d0)then
             ptest(5)=sqrt(ptest(5))
           else
             if(ish.ge.1)write(ifch,*)'Precision problem in jintpo, p:',
      &           (ptest(k),k=1,5)
             ptest(5)=0d0
           endif
           be(1)=0.d0
           be(2)=0.d0
           do i=ini,ifi
             if(ityptl(i).eq.60)then
               call utlob4(1,ptest(1),ptest(2),ptest(3),ptest(4),ptest(5)
      $             ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
               be(1)=be(1)+dble(pptl(1,i))
               be(2)=be(2)+dble(pptl(2,i))
             endif
           enddo
 c shift nclu particles to have sum_pt=0. and boost back in global cms
           pt1shift=-sngl(be(1)/dble(nclu))
           pt2shift=-sngl(be(2)/dble(nclu))
           do i=ini,ifi
             if(ityptl(i).eq.60)then
               pptl(1,i)=pptl(1,i)+pt1shift
               pptl(2,i)=pptl(2,i)+pt2shift
               pptl(4,i)=sqrt(pptl(1,i)**2+pptl(2,i)**2
      &             +pptl(3,i)**2+pptl(5,i)**2)
              call utlob4(-1,ptest(1),ptest(2),ptest(3),ptest(4),ptest(5)
      &             ,pptl(1,i),pptl(2,i),pptl(3,i),pptl(4,i))
             endif
           enddo
         endif
 
       endif
 
 
 c      if(ntrt.gt.0)print *,"jintpo rescaling :",float(ntrr)/float(ntrt)
 c     define life time
       n=0
       do i=ini,ifi
         if(ityptl(i).eq.60)then
           n=n+1
           r=1.15*rini*yrad(n)   !yrad=y/ymax
           tau=2.25/sqrt(yrad(n)**2+0.04)-0.75
           z=xorptl(3,i)
           t=xorptl(4,i)
 !         zeta=0.5*log((t+z)/(t-z))-0.5*delzet+2*0.5*delzet*rangen()
         test=(pptl(4,i)-pptl(3,i))*(pptl(4,i)+pptl(3,i))
         if(test.gt.0.)then
           zeta=0.5*log((pptl(4,i)+pptl(3,i))
      &         /(pptl(4,i)-pptl(3,i)))
         else    !in case of precision problem (but not always good neither for p<0
           pt=sqrt(pptl(2,i)**2+pptl(1,i)**2)
           zeta=0.5*log(1+2*pptl(3,i)*(pptl(4,i)+pptl(3,i))
      &   /(pt*pt+pptl(5,i)**2))
         endif
           z=tau*sinh(zeta)
           t=tau*cosh(zeta)
           xorptl(1,i)=xorptl(1,i)+r*cos(phirad(n))
           xorptl(2,i)=xorptl(2,i)+r*sin(phirad(n))
           xorptl(3,i)=z
           xorptl(4,i)=t
         endif
       enddo
 
 
       if(ish.ge.5)then
         do k=1,5
           pptl(k,nptl+nall+2)=0
           do ii=nptlbcf+1,nptl
             pptl(k,nptl+nall+2)=pptl(k,nptl+nall+2)+pptl(k,ii)
           enddo
         enddo
         iorptl(nptl+nall+2)=nptlbcf+1
         jorptl(nptl+nall+2)=nptl
         call alist2('longitudinal and radial flow&',nptl+1
      &       ,nptl+nall,nptlbcf+1,nptl)
         call alist2('momentum sum&',nptl+nall+1,nptl+nall+1
      &       ,nptl+nall+2,nptl+nall+2)
         write(ifch,'(1x,50a1/)')('-',k=1,50)
       endif
       
 
       endif      !ioclude=3 and flow
 
 
       do n=nptlbcf+1,nptl
         if(ioclude.ne.3)then
           iorptl(n)=nptlb+1
           jorptl(n)=nptlbcf
           rinptl(n)=rinptl((iorptl(n)+jorptl(n))/2)
         else
           rinptl(n)=rinptl(iorptl(n))
         endif
         istptl(n)=0
         ifrptl(1,n)=0
         ifrptl(2,n)=0
         tivptl(1,n)=xorptl(4,n)
         call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm)
         r=rangen()
         tivptl(2,n)=tivptl(1,n)+taugm*(-alog(r))
         radptl(n)=0.
         dezptl(n)=0.
         itsptl(n)=0
       enddo
 
       endif
 
 c Decay droplets not included in clusters
       iret=0
       do mm=1,nptla
         nptlb=nptl
         if(istptl(mm).eq.10)then
           if(ish.ge.5)write(ifch,*)'Decay remaining droplet :',mm
           if(nptlb.gt.mxptl-10)
      &    call utstop('jintpo: mxptl too small (2)&')
           if(ioclude.eq.3)then
             call hnbaaa(mm,iret)
           else
             call DropletDecay(mm,iret) !Decay remn
             iret=0
           endif
           if(iret.eq.0.and.nptl.ne.nptlb)then ! ---successful decay---
             istptl(mm)=istptl(mm)+1
             ifrptl(1,mm)=nptlb+1
             ifrptl(2,mm)=nptl
             t=tivptl(2,mm)
             x=xorptl(1,mm)+(t-xorptl(4,mm))*pptl(1,mm)/pptl(4,mm)
             y=xorptl(2,mm)+(t-xorptl(4,mm))*pptl(2,mm)/pptl(4,mm)
             z=xorptl(3,mm)+(t-xorptl(4,mm))*pptl(3,mm)/pptl(4,mm)
             do 21 n=nptlb+1,nptl
               iorptl(n)=mm
               jorptl(n)=0
               istptl(n)=0
               ifrptl(1,n)=0
               ifrptl(2,n)=0
               radius=0.8*sqrt(rangen())
               phi=2*pi*rangen()
               ti=t
               zi=z
               xorptl(1,n)=x + radius*cos(phi)
               xorptl(2,n)=y + radius*sin(phi)
               xorptl(3,n)=zi
               xorptl(4,n)=ti
               iioo=mm
               zor=dble(xorptl(3,iioo))
               tor=dble(xorptl(4,iioo))
               r=rangen()
               tauran=-taurea*alog(r)
               call jtaix(n,tauran,zor,tor,zis,tis)
               tivptl(1,n)=amax1(ti,tis)
               call idtau(idptl(n),pptl(4,n),pptl(5,n),taugm)
               r=rangen()
               tivptl(2,n)=t+taugm*(-alog(r))
               radptl(n)=0.
               dezptl(n)=0.
               itsptl(n)=0
               rinptl(nptl)=-9999
    21       continue
           else                  ! Unsuccessful decay
             if(ish.ge.1)write(ifch,*)
      *         '***** Unsuccessful remnant cluster decay'
      *             ,' --> redo event.'
           endif
         endif
       enddo
 
 
 
  1000 continue
       call utprix('jintpo',ish,ishini,4)
       end
 
 cc-----------------------------------------------------------------------
 c      subroutine jrad(i,nq,na,jc,rad)
 cc-----------------------------------------------------------------------
 cc     return hadron radius (data taken from huefner and povh)
 cc-----------------------------------------------------------------------
 c      include 'epos.inc'
 c      integer jc(nflav,2),kc(nflav)
 c
 c      id=iabs(idptl(i))
 c      am=pptl(5,i)
 c      if(id.lt.10000)then
 c       k=mod(id,10)
 c      else
 c       k=1
 c      endif
 c      do l=1,nflav
 c       kc(l)=iabs(jc(l,1)-jc(l,2))
 c      enddo
 c
 c      if(nq.eq.0)then   ! mesons
 c       if(kc(1).eq.0.and.kc(2).eq.0.and.kc(3).eq.0.and.kc(4).eq.0)then
 c        if(k.eq.0)then           ! flavor singlet pseudoscalar mesons
 c         if(am.ge.0.000)then
 c          rad=0.64                 ! pi0
 c          if(am.ge.0.500)then
 c           rad=0.60                ! eta
 c           if(am.ge.0.900)then
 c            rad=0.40               ! eta prime
 c            if(am.ge.2.900)then
 c             rad=0.17              ! eta charm
 c            endif
 c           endif
 c          endif
 c         else
 c          write(ifch,*)
 c     *    'i:',i,' id:',idptl(i),' k:',k,' m:',am
 c          write(ifch,*)'jc:',(jc(l,1),l=1,6),(jc(l,2),l=1,6)
 c          call utstop('jrad: meson radius not defined&')
 c         endif
 c        else                    ! flavor singlet vector mesons
 c         if(am.ge.0.000)then
 c          rad=0.72                ! rho,omega
 c          if(am.ge.1.000)then
 c           rad=0.46               ! phi
 c           if(am.ge.3.000)then
 c            rad=0.20              ! J/psi
 c           endif
 c          endif
 c         else
 c          write(ifch,*)
 c     *    'i:',i,' id:',idptl(i),' k:',k,' m:',am
 c          write(ifch,*)'jc:',(jc(l,1),l=1,6),(jc(l,2),l=1,6)
 c          call utstop('jrad: meson radius not defined&')
 c         endif
 c        endif
 c       elseif(kc(3).eq.0.and.kc(4).eq.0)then  ! nonstrange, noncharmed
 c        if(k.eq.0)then
 c         rad=0.64  ! pi
 c        else
 c         rad=0.72  ! resonances
 c        endif
 c       elseif(kc(3).ne.0.and.kc(4).eq.0)then  ! strange
 c        if(k.eq.0)then
 c         rad=0.59  ! kaons
 c        else
 c         rad=0.68  ! kaon resonances
 c        endif
 c       else                                   ! charmed
 c        write(ifch,*)'i:',i,' id:',idptl(i)
 c        call utstop('jrad: radius of meson not defined&')
 c       endif
 c      else   !baryons
 c       if(kc(4).gt.0)then       ! charmed
 c        write(ifch,*)
 c     *  'i:',i,' id:',idptl(i),' k:',k,' m:',am
 c        write(ifch,*)'i:',i,' id:',idptl(i)
 c        call utstop('jrad: radius of charmed baryon not defined&')
 c       elseif(kc(3).eq.0)then   ! nonstrange
 c        if(k.eq.0)then
 c         rad=0.82  !nucleons
 c        else
 c         rad=1.00  !resonances
 c        endif
 c       elseif(kc(3).eq.1)then   ! strange
 c        if(k.eq.0)then
 c         rad=0.76  !lambda, sigma
 c        else
 c         rad=0.93  !resonances
 c        endif
 c       elseif(kc(3).eq.2)then   ! double strange
 c        if(k.eq.0)then
 c         rad=0.71  !cascades
 c        else
 c         rad=0.87  !resonances
 c        endif
 c       elseif(kc(3).ge.3)then   ! triple strange
 c        rad=0.79  !omega
 c       else
 c        write(ifch,*)
 c     *  'i:',i,' id:',idptl(i),' k:',k,' m:',am
 c        write(ifch,*)
 c     *  'q:',(jc(l,1),l=1,6),' qbar:',(jc(l,2),l=1,6),
 c     *  ' |q-qbar|:',(kc(l),l=1,6)
 c        call utstop('jrad: should not happen&')
 c       endif
 cc string fragments with |#q|>3
 c       if(na.gt.3)then
 c        a=(na/3.)**(1./3.)
 c        if(ish.ge.7)then
 c         call utmsg('jrad ')
 c         write(ifch,*)
 c     *   'i:',i,' id:',idptl(i),' k:',k,' m:',am
 c         write(ifch,*)
 c     *   'q:',(jc(l,1),l=1,6),' qbar:',(jc(l,2),l=1,6),
 c     *   ' |q-qbar|:',(kc(l),l=1,6)
 c         write(ifch,*)'nq:',nq,' na:',na,' r:',rad,' ar:',a*rad
 c         call utmsgf
 c        endif
 c        rad=rad*a
 c       endif
 c      endif
 c
 c      if(ish.ge.7)then
 c       write(ifch,*)
 c     * 'i:',i,' id:',idptl(i),' k:',k,' m:',am,' rad:',rad
 c       write(ifch,*)'jc:',(jc(l,1),l=1,6),(jc(l,2),l=1,6)
 c      endif
 c
 c      return
 c      end
 c
 c-----------------------------------------------------------------------
       subroutine jresc
 c-----------------------------------------------------------------------
       include 'epos.inc'
       double precision pa(5),pj(5)
       integer ipptl(mxptl)
 
       call utpri('jresc ',ish,ishini,4)
 
       iret=0
       nptlpt=maproj+matarg
       np=0
       do i=nptlpt+1,nptl
        if(istptl(i).eq.0
      * .and.idptl(i).lt.10000.and.pptl(5,i).gt.0.01)then
         np=np+1
         ipptl(np)=i
        endif
       enddo
       if(np.lt.2)goto1001
       do ii=1,np
        i=ipptl(ii)
        if(mod(iabs(idptl(i)),10).lt.2)then
         call idmass(idptl(i),ami)
         dm=abs(ami-pptl(5,i))
         if(dm.gt.0.001)then
          ntry=0
 1        continue
          ntry=ntry+1
 2        jj=1+int(rangen()*np)
          j=ipptl(jj)
          if(ish.ge.4)write(ifch,*)i,j,istptl(j)
          if(j.eq.i)goto2
          if(mod(iabs(idptl(j)),10).lt.2)then
           call idmass(idptl(j),amj)
          else
           amj=pptl(5,j)
          endif
          do l=1,5
           pa(l)=dble(pptl(l,i))
           pj(l)=dble(pptl(l,j))
          enddo
          if(ish.ge.4)write(ifch,'(70a1)')('-',l=1,70)
          if(ish.ge.4)write(ifch,11)i,idptl(i),'before:',pa,'want:',ami
          if(ish.ge.4)write(ifch,11)j,idptl(j),'before:',pj,'want:',amj
          call jrescl(pa,dble(ami),pj,dble(amj),iret)
          if(iret.eq.1)then
           if(ntry.le.50)then
            goto1
           else
            goto1001
           endif
          endif
          if(ish.ge.4)write(ifch,11)i,idptl(i),' after:',pa
          if(ish.ge.4)write(ifch,11)j,idptl(j),' after:',pj
          if(ish.ge.4)write(ifch,'(70a1)')('-',l=1,70)
          do l=1,5
           pptl(l,i)=sngl(pa(l))
           pptl(l,j)=sngl(pj(l))
          enddo
         endif
        endif
       enddo
 11    format(i5,1x,i5,1x,a,1x,5(d8.2,1x),a,1x,e8.2)
 
 1000  continue
       call utprix('jresc ',ish,ishini,4)
       return
 
 1001  continue
       if(ish.ge.1)then
         write(ifmt,'(a)')'jresc: could not put on shell'
       endif
       goto1000
 
       end
 
 c-----------------------------------------------------------------------
       subroutine jrescl(p1,am1,p2,am2,iret)
 c-----------------------------------------------------------------------
 c rescale momenta of two particles such that the masses assume given
 c values.
 c input:
 c   p1, p2: momenta of the two particles
 c   am1, am2: desired masses of the two particles
 c output:
 c   p1, p2: new momenta of the two particles
 c-----------------------------------------------------------------------
       include 'epos.inc'
       double precision p1(5),p2(5)
      *                ,p1n(5),p2n(5)
      *                ,a1,a2,a12,am1,am2
      *                ,b1,b2,c,d,e,f,g,p,q,r
 
       call utpri('jrescl',ish,ishini,7)
 
       iret=0
       a1=p1(5)**2
       a2=p2(5)**2
       a12=p1(4)*p2(4)-p1(3)*p2(3)-p1(2)*p2(2)-p1(1)*p2(1)
       if(a12.le.(a1+a2))then
        if(ish.ge.7)write(ifch,*)'a_12 < a_1 + a_2'
        if(ish.ge.7)write(ifch,*)a12,' < ',a1+a2
 c      goto1001
       endif
 
 11    format(5(d9.3,1x))
       if(ish.ge.7)write(ifch,11)p1,a1
       if(ish.ge.7)write(ifch,11)p2,a2
       if(ish.ge.7)write(ifch,*)a12
 
       c=(a1+a12)/(a2+a12)
       d=(a1-am1**2-a2+am2**2)/(a2+a12)*0.5d0
 
       e=a1-2d0*a12*c+a2*c**2
       f=2d0*(a1-a12*(c+d)+a2*c*d)
       g=a1-2d0*a12*d+a2*d**2-am1**2
 
       p=f/e
       q=g/e
       r=p**2-4d0*q
 
       if(ish.ge.7)write(ifch,*)'c:',c,' d:',d
       if(ish.ge.7)write(ifch,*)'e:',e,' f:',f,' g:',g
       if(ish.ge.7)write(ifch,*)'p:',p,' q:',q,' r:',r
       if(r.lt.0d0)goto1001
 
       b1=-0.5d0*(p-dsqrt(r))
 
       b2=b1*c+d
 
       if(ish.ge.7)write(ifch,*)'b_1:',b1,' b_2:',b2
 
       do i=1,4
        p1n(i)=(1d0+b1)*p1(i)-b2*p2(i)
        p2n(i)=(1d0+b2)*p2(i)-b1*p1(i)
       enddo
 
       a1=p1n(4)**2-p1n(3)**2-p1n(2)**2-p1n(1)**2
       a2=p2n(4)**2-p2n(3)**2-p2n(2)**2-p2n(1)**2
       if(a1.gt.0d0.and.a2.gt.0d0)then
        do i=1,4
         p1(i)=p1n(i)
         p2(i)=p2n(i)
        enddo
        p1(5)=dsqrt(a1)
        p2(5)=dsqrt(a2)
        if(ish.ge.7)write(ifch,11)p1,a1
        if(ish.ge.7)write(ifch,11)p2,a2
       else
        goto1001
       endif
 
       if(p1(4).lt.0..or.p2(4).lt.0.)goto1001
 
 1000  continue
       call utprix('jrescl',ish,ishini,7)
       return
 
 1001  continue
       iret=1
       goto1000
       end
 
 c-----------------------------------------------------------------------
       subroutine jtain(i,x,y,z,t,n,iopt)
 c-----------------------------------------------------------------------
 c returns intersection (x,y,z,t) of ptl-i-trajectory with taus-line.
 c input:
 c   i: particle number
 c   iopt: formation time considered (0) or not (1)
 c output:
 c   x,y,z,t: 4-vector of intersection point
 c   n: exit code
 c       n=0: ok
 c       n=1: ptl lives later
 c       n=2: ptl lives earlier
 c       n=9: tiv1>tiv2
 c-----------------------------------------------------------------------
       include 'epos.inc'
       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
       common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat
       double precision vv,zza,zz,tt,xo3,xo4,ti1,ti2,derr,dd
       double precision ttp,zzp,ttt,zzt,vvt,vvp,spt2m2E,p4
       common/ctfi/tin,tfi
       double precision ttau0
       common/cttau0/ttau0
 
       n=0
 
       tin=0
       tfi=0
 
       derr=1d-2
       ttp=tpro*ttaus
       zzp=zpro*ttaus
       ttt=ttar*ttaus
       zzt=ztar*ttaus
       vv=sign(min(1.d0,abs(dble(pptl(3,i)))/dble(pptl(4,i)))
      &                ,dble(pptl(3,i)))
 
 
       if(abs(vv).ge.1.d0)then
         spt2m2E=dble(pptl(1,i)*pptl(1,i)+pptl(2,i)*pptl(2,i)
      &              +pptl(5,i)*pptl(5,i))
 c        if(pptl(4,i).le.0.)then
           p4=sqrt(dble(pptl(3,i)*pptl(3,i))+spt2m2E)
 c        else
 c          p4=dble(pptl(4,i))
 c        endif
 ctp to avoid precision problem, replace abs(p3)/p4 by sqrt(1-(pt2+m2)/E2)
         spt2m2E=min(1.d0,sqrt(spt2m2E)/p4)
         vv=sign(sqrt((1d0+spt2m2E)*(1d0-spt2m2E)),dble(pptl(3,i)))
       endif
       xo3=dble(xorptl(3,i))
       xo4=dble(xorptl(4,i))
       zza=xo3-xo4*vv
       if(iopt.eq.0)then
         ti1=dble(tivptl(1,i))
       elseif(iopt.eq.1)then
         ti1=dble(xo4)
       else
         ti1=0
         call utstop("Wrong iopt in jtain !&")
       endif
       ti2=dble(tivptl(2,i))
 
       if(ti1.gt.ti2)then
         n=9
         goto1
       endif
 
       zfi=sngl(xo3+(ti2-xo4)*vv)
       call jtaus(zfi,tzfi,szfi)
       tfi=tzfi
       if(tfi.ge.sngl(ti2))then
         n=2
         goto1
       endif
 
       zin=sngl(xo3+(ti1-xo4)*vv)
       call jtaus(zin,tzin,szin)
       tin=tzin
       if(tin.le.sngl(ti1))then
         n=1
         goto1
       endif
 
 
     1 continue
 
            if(ttaus.le.ttau0)then
       tt=ttaus
       zz=xo3+(tt-xo4)*vv
       if(tt.lt.ti1.and.n.eq.0)n=1
       if(tt.ge.ti2.and.n.eq.0)n=2
       goto1000
            else
       vvt=zzt/ttt
       vvp=zzp/ttp
       tt=(ttt+(zza-zzt)*vvt)/(1-vv*vvt)
       zz=xo3+(tt-xo4)*vv
       if(zz.le.zzt)then
       if(tt.lt.ti1.and.n.eq.0)n=1
       if(tt.ge.ti2.and.n.eq.0)n=2
       goto1000
       endif
       tt=(ttp+(zza-zzp)*vvp)/(1-vv*vvp)
       zz=xo3+(tt-xo4)*vv
       if(zz.ge.zzp)then
       if(tt.lt.ti1.and.n.eq.0)n=1
       if(tt.ge.ti2.and.n.eq.0)n=2
       goto1000
       endif
       dd=1-vv**2
       if(sngl(dd).eq.0..and.vv.gt.0.)then
       tt=-(ttaus**2+zza**2)/2d0/zza
       elseif(sngl(dd).eq.0..and.vv.lt.0.)then
       tt=(ttaus**2+zza**2)/2d0/zza
       else
       tt=(zza*vv+dsqrt(zza**2+ttaus**2*dd))/dd
       endif
       zz=xo3+(tt-xo4)*vv
       if(tt.lt.ti1.and.n.eq.0)n=1
       if(tt.ge.ti2.and.n.eq.0)n=2
         if(dabs(ttaus**2-(tt+zz)*(tt-zz)).gt.derr*ttaus**2.and.
      *dabs(ttaus**2-(tt+zz)*(tt-zz)).gt.derr)then
       if(ish.ge.1)then
       call utmsg('jtain')
       write(ifch,*)'*****  ttaus**2 .ne. (tt+zz)*(tt-zz)'
       write(ifch,*)sngl(ttaus**2),sngl((tt+zz)*(tt-zz))
       call utmsgf
       endif
       goto1000
         endif
            endif
 
 1000  t=sngl(tt)
       z=sngl(zz)
       x=xorptl(1,i)+(t-xorptl(4,i))*pptl(1,i)/pptl(4,i)
       y=xorptl(2,i)+(t-xorptl(4,i))*pptl(2,i)/pptl(4,i)
       return
       end
 
 c-----------------------------------------------------------------------
       subroutine jtaix(i,tau,zor,tor,z,t)
 c-----------------------------------------------------------------------
 c     returns intersection z,t of ptl-i-trajectory with hyperbola h.
 c        h: (t-tor)**2-(z-zor)**2=tau**2 .
 c        zor, tor double precision.
 c-----------------------------------------------------------------------
       include 'epos.inc'
       double precision tor,zor,tors,zors,vv,cc,dd,ttau,derr,tt,zz
       derr=1d-3
       ttau=dble(tau)
       zors=dble(xorptl(3,i))-zor
       tors=dble(xorptl(4,i))-tor
       vv=dble(pptl(3,i)/pptl(4,i))
       vv=dmin1(vv,1d0)
       vv=dmax1(vv,-1d0)
       cc=zors-tors*vv
       dd=1d0-vv**2
       dd=dmax1(dd,0d0)
            if(dd.eq.0d0.and.cc.eq.0d0)then
       if(tau.eq.0.)tt=0d0
       if(tau.ne.0.)tt=dble(ainfin)
       zz=tt
       goto1000
            elseif(dd.eq.0d0)then
       tt=-(ttau**2+cc**2)/2d0/cc/vv
            elseif(dd.lt.1e-8)then
       tt=-(ttau**2+cc**2)/2d0/cc/vv
       call utmsg('jtaix')
       write(ifch,*)'*****  dd = ',dd,'    treated as zero'
       call utmsgf
            else
       tt=(cc*vv+dsqrt(cc**2+ttau**2*dd))
       tt=tt/dd
            endif
       zz=cc+tt*vv
       if(dabs(ttau**2-(tt+zz)*(tt-zz)).gt.derr*ttau**2
      *.and.dabs(ttau**2-(tt+zz)*(tt-zz)).gt.derr
      *.and.tors**2-zors**2.lt.1e6)then
       if(ish.ge.2)then
       call utmsg('jtaix')
       write(ifch,*)'*****  ttau**2 .ne. (tt+zz)*(tt-zz)'
       write(ifch,*)sngl(ttau**2),sngl((tt+zz)*(tt-zz))
       write(ifch,*)'tau,t,z:'
       write(ifch,*)tau,tt,zz
       write(ifch,*)'#,id(ptl):',i,idptl(i)
       write(ifch,*)'zor,tor(str):',zor,tor
       write(ifch,*)'zors,tors,p,e(ptl):'
       write(ifch,*)sngl(zors),sngl(tors),pptl(3,i),pptl(4,i)
       call utmsgf
       endif
       endif
 1000  z=sngl(zz+zor)
       t=sngl(tt+tor)
       return
       end
 
 c-----------------------------------------------------------------------
       subroutine jtaug(su,so,g,y)
 c-----------------------------------------------------------------------
 c  returns g factor and rapidity y for given su, so
 c-----------------------------------------------------------------------
       include 'epos.inc'
       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
       common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat
       double precision ttp,zzp,ttt,zzt,ssp,sst,ssu,sso,ss1,ss2,gg
      *,ssav,yyav,hh
       double precision ttau0
       common/cttau0/ttau0
 
       ssu=dble(su)
       sso=dble(so)
 
       if(ssu.ge.sso)then
       sso=(ssu+sso)*0.5d0 + dble(abs(dezzer))*ttaus*0.5d0
       ssu=(ssu+sso)*0.5d0 - dble(abs(dezzer))*ttaus*0.5d0
       so=real(sso)
       su=real(ssu)
       endif
       if(ssu.ge.sso)then
         print*,ssu,sso,dble(abs(dezzer))*ttaus*0.5d0
         stop'STOP: sr jtaug: ssu.ge.sso'
       endif
 
       g=1
 
       if(ttaus.le.ttau0)return
 
       ttp=tpro*ttaus
       zzp=zpro*ttaus
       ttt=ttar*ttaus
       zzt=ztar*ttaus
       ssp=ttaus*0.5d0*dlog((ttp+zzp)/(ttp-zzp))
       sst=ttaus*0.5d0*dlog((ttt+zzt)/(ttt-zzt))
 
       ssav=(ssu+sso)/2d0
       yyav=ssav/ttaus
       if(ssav.ge.ssp)yyav=detap
       if(ssav.le.sst)yyav=detat
 
       gg=0
       if(ssu.lt.sst)gg=gg + dcosh(detat-yyav) * (dmin1(sst,sso)-ssu)
       if(sso.gt.ssp)gg=gg + dcosh(detap-yyav) * (sso-dmax1(ssp,ssu))
       if(ssu.lt.ssp.and.sso.gt.sst)then
       ss1=dmax1(ssu,sst)
       ss2=dmin1(sso,ssp)
       gg=gg+ttaus*( dsinh(ss2/ttaus-yyav)-dsinh(ss1/ttaus-yyav) )
       endif
       gg=gg/(sso-ssu)
 
       hh=0
       if(ssu.lt.sst)hh=hh + dsinh(detat-yyav) * (dmin1(sst,sso)-ssu)
       if(sso.gt.ssp)hh=hh + dsinh(detap-yyav) * (sso-dmax1(ssp,ssu))
       if(ssu.lt.ssp.and.sso.gt.sst)then
       ss1=dmax1(ssu,sst)
       ss2=dmin1(sso,ssp)
       hh=hh+ttaus*( dcosh(ss2/ttaus-yyav)-dcosh(ss1/ttaus-yyav) )
       endif
       hh=hh/(sso-ssu)
 
       yyav=yyav+0.5d0*dlog((gg+hh)/(gg-hh))
 
       gg=0
       if(ssu.lt.sst)gg=gg + dcosh(detat-yyav) * (dmin1(sst,sso)-ssu)
       if(sso.gt.ssp)gg=gg + dcosh(detap-yyav) * (sso-dmax1(ssp,ssu))
       if(ssu.lt.ssp.and.sso.gt.sst)then
       ss1=dmax1(ssu,sst)
       ss2=dmin1(sso,ssp)
       gg=gg+ttaus*( dsinh(ss2/ttaus-yyav)-dsinh(ss1/ttaus-yyav) )
       endif
       gg=gg/(sso-ssu)
 
       g=sngl(gg)
       y=sngl(yyav)
       return
       end
 
 c-----------------------------------------------------------------------
       subroutine jtaui(s,ts,zs)
 c-----------------------------------------------------------------------
 c  returns time ts and coord zs corresponding to ttaus and inv. length s
 c-----------------------------------------------------------------------
 
       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
       common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat
       double precision ttau0
       common/cttau0/ttau0
 
       double precision ttp,zzp,ttt,zzt,ssp,sst,ss,zeta
 
       zs=s
       ts=sngl(ttaus)
 
       if(ttaus.le.ttau0)return
 
       ttp=tpro*ttaus
       zzp=zpro*ttaus
       ttt=ttar*ttaus
       zzt=ztar*ttaus
       ssp=ttaus*0.5d0*dlog((ttp+zzp)/(ttp-zzp))
       sst=ttaus*0.5d0*dlog((ttt+zzt)/(ttt-zzt))
       ss=dble(s)
 
            if(ss.le.sst)then
       zs=sngl(zzt+ttar*(ss-sst))
       ts=sngl(ttt+(dble(zs)-zzt)*zzt/ttt)
            elseif(ss.ge.ssp)then
       zs=sngl(zzp+tpro*(ss-ssp))
       ts=sngl(ttp+(dble(zs)-zzp)*zzp/ttp)
            else
       zeta=ss/ttaus
       ts=sngl(ttaus*dcosh(zeta))
       zs=sngl(ttaus*dsinh(zeta))
            endif
 
       return
       end
 
 c-----------------------------------------------------------------------
       subroutine jtauin
 c-----------------------------------------------------------------------
 c initializes equal time hyperbola at ttaus
 c-----------------------------------------------------------------------
       include 'epos.inc'
       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
       common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat
       double precision ttau0,rcproj,rctarg
       common/geom1/rcproj,rctarg
       common/cttau0/ttau0
 
       call utpri('jtauin',ish,ishini,6)
 
       if(ttaus.gt.ttau0)then
        if(rcproj.gt.1d-10)then
         detap=dmin1(dble((ypjtl-yhaha)*etafac),dlog(ttaus/rcproj))
        else
         detap=dble((ypjtl-yhaha)*etafac)
        endif
        if(rctarg.gt.1d-10)then
         detat=dmax1(dble(-yhaha*etafac),dlog(rctarg/ttaus))
        else
         detat=dble(-yhaha*etafac)
        endif
        tpro=dcosh(detap)
        zpro=dsinh(detap)
        ttar=dcosh(detat)
        ztar=dsinh(detat)
       else
        detap=0d0
        detat=0d0
        tpro=0d0
        zpro=0d0
        ttar=0d0
        ztar=0d0
       endif
 
       if(ish.ge.6)then
        write(ifch,*)'hyperbola at tau=',ttaus
        write(ifch,*)'r_p:',rcproj,' r_t:',rctarg
        write(ifch,*)'y_p:',detap,' y_t:',detat
        write(ifch,*)'t_p:',tpro,' z_p:',zpro
        write(ifch,*)'t_t:',ttar,' z_t:',ztar
       endif
 
       call utprix('jtauin',ish,ishini,6)
       return
       end
 
 c-----------------------------------------------------------------------
       subroutine jtaus(z,tz,sz)
 c-----------------------------------------------------------------------
 c  returns time tz and inv length sz corresponding to ttaus and z
 c-----------------------------------------------------------------------
       include 'epos.inc'
       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
       common/cttaus/tpro,zpro,ttar,ztar,ttaus,detap,detat
       double precision ttau0
       common/cttau0/ttau0
 
       double precision ttp,zzp,ttt,zzt,zz,tzz
 
       tz=sngl(ttaus)
       sz=z
 
       if(ttaus.le.ttau0)return
 
       ttp=tpro*ttaus
       zzp=zpro*ttaus
       ttt=ttar*ttaus
       zzt=ztar*ttaus
       zz=dble(z)
 
            if(zz.le.zzt)then
       tz=sngl(ttt+(zz-zzt)*zzt/ttt)
       sz=sngl(ttaus*detat+(zz-zzt)/ttar)
            elseif(zz.ge.zzp)then
       tz=sngl(ttp+(zz-zzp)*zzp/ttp)
       sz=sngl(ttaus*detap+(zz-zzp)/tpro)
            else
       if(sngl(ttaus).ge.ainfin)then
       tz=sngl(ttaus)
       sz=0.
       if(ish.ge.1)then
       call utmsg('jtaus')
       write(ifch,*)'*****  large ttaus; set tz=ttaus, sz=0'
       write(ifch,*)'ttaus=',ttaus,'zz=',zz
       call utmsgf
       endif
       else
       tzz=dsqrt(ttaus**2+zz**2)
       tz=sngl(tzz)
       sz=sngl(ttaus*0.5d0*dlog((tzz+zz)/(tzz-zz)))
       endif
            endif
 
       return
       end
 
 c-----------------------------------------------------------------------
       subroutine jtaux(t,z,ttaux)
 c-----------------------------------------------------------------------
 c  returns ttaux (-> tau-line) for given t and z.
 c  ttaux: double precision.
 c-----------------------------------------------------------------------
       double precision ttaux,tt,zz,rcproj,rctarg,zt1,zp1,zt2,zp2,ttau0
       common/geom1/rcproj,rctarg
       common/cttau0/ttau0
       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
       common/cttaus/   tpro,zpro,ttar,ztar,ttaus,detap,detat
 
       tt=dble(t)
       zz=dble(z)
 
       if(tt.gt.ttau0)then
        zt1=rctarg-tt
        zp1=tt-rcproj
        zt2=ztar/ttar*tt
        zp2=zpro/tpro*tt
        if(zz.le.dmax1(zt1,zt2))then
         if(zt1.gt.zt2)then
          ttaux=rctarg*dsqrt((tt-zz)/(2d0*rctarg-tt-zz))
         else
          ttaux=(ttar*tt-ztar*zz)/(ttar**2-ztar**2)
         endif
        elseif(zz.ge.dmin1(zp1,zp2))then
         if(zp1.lt.zp2)then
          ttaux=rcproj*dsqrt((tt+zz)/(2d0*rcproj-tt+zz))
         else
          ttaux=(tpro*tt-zpro*zz)/(tpro**2-zpro**2)
         endif
        else
         ttaux=dsqrt(tt**2-zz**2)
        endif
       else
        ttaux=tt
       endif
 
       return
       end
 
 c-----------------------------------------------------------------------
       subroutine xjden1(ii,itau,x,y,rad,o,u)
 c-----------------------------------------------------------------------
 c ii=0: initialization
 c ii=1: determining dense regions in space of individual events
 c       x,y,rad: tranverse coordinates and radius of particle i
 c       o,u: specifies long range: u < s < o (s: long coordinate)
 c ii=2: plot of individual event
 c       xpar1: itau ; valid: 1,..,10
 c       xpar2: z-range: -xpar2 < z < xpar2
 c       xpar3, x-range: -xpar3 < x < xpar3
 c       xpar4, y-range: -xpar4 < y < xpar4
 c-----------------------------------------------------------------------
       include "epos.inc"
       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
       common/cttaus/   tpro,zpro,ttar,ztar,ttaus,detap,detat
 
       if(idensi.ne.1)stop'STOP in xjden1: idensi must be set 1'
 
       dlcoox=0.5
       dlcooy=0.5
 
            if(ii.eq.0)then
 
       do i=1,nzeta
       do j=1,mxcoox
       do k=1,mxcooy
       kdensi(itau,i,j,k)=0
       enddo
       enddo
       enddo
 
            elseif(ii.eq.1)then
 
       if(itau.lt.1.or.itau.gt.mxtau)return
 
       tau=sngl(ttaus)
       zu=u/tau
       zo=o/tau
 
             do 1 i=1,nzeta
       zi=-nzeta*dlzeta/2-dlzeta/2+i*dlzeta
       if(zu.gt.zi.or.zo.lt.zi)goto1
       do 2 j=1,mxcoox
       xj=-mxcoox*dlcoox/2-dlcoox/2+j*dlcoox
       do 3 k=1,mxcooy
       yk=-mxcooy*dlcooy/2-dlcooy/2+k*dlcooy
       if((x-xj)**2+(y-yk)**2.gt.rad**2)goto3
       kdensi(itau,i,j,k)=1
     3 continue
     2 continue
     1 continue
 
            elseif(ii.eq.2)then
 
       itaux=nint(xpar1)
       if(itaux.gt.mxtau)stop'STOP in xjden1: itaux too large'
 
       iz=nint(xpar2/dlzeta)
       ix=nint(xpar3/dlcoox)
       iy=nint(xpar4/dlcooy)
       if(iz.gt.nzeta/2)stop'STOP in xjden1: zeta-range too large'
       if(ix.gt.mxcoox/2)stop'STOP in xjden1: x-range too large'
       if(iy.gt.mxcooy/2)stop'STOP in xjden1: y-range too large'
 
       do k=mxcooy/2+1-iy,mxcooy/2+iy
       write(ifhi,'(a,f7.2)')      '! tau: ',tauv(itaux)
       write(ifhi,'(a)')      'openhisto'
       write(ifhi,'(a,2f7.2)')'xrange',-iz*dlzeta,iz*dlzeta
       write(ifhi,'(a,2f7.2)')'yrange',-ix*dlcoox,ix*dlcoox
       write(ifhi,'(a)')      'set ityp2d 3'
       write(ifhi,'(a)') 'txt  "xaxis space-time rapidity [z]"'
       write(ifhi,'(a)') 'txt  "yaxis transverse coordinate x (fm)"'
       write(ifhi,'(a,i4)')   'array2d',2*iz
       do j=mxcoox/2+1-ix,mxcoox/2+ix
       write(ifhi,'(40i2)')    (kdensi(itaux,i,j,k),
      *                        i=nzeta/2+1-iz,nzeta/2+iz)
       enddo
       write(ifhi,'(a)')       '  endarray'
       write(ifhi,'(a)')       'closehisto plot2d'
       enddo
 
            else
 
       stop'STOP in xjden1: wrong option'
 
            endif
 
       return
       end
 
 c-----------------------------------------------------------------------
       subroutine xjden2(ii,itau,x,y,rad,s)
 c-----------------------------------------------------------------------
 c ii=0: initialization
 c ii=1: determining dense regions in space of individual events
 c       x,y,rad: tranverse coordinates and radius of particle i
 c       s: long coordinate
 c ii=2: plot of individual event
 c       xpar1: itau ; valid: 1,..,10
 c       xpar2: s-range: -xpar2 < s < xpar2
 c       xpar3, x-range: -xpar3 < x < xpar3
 c       xpar4, y-range: -xpar4 < y < xpar4
 c-----------------------------------------------------------------------
       include "epos.inc"
       double precision tpro,zpro,ttar,ztar,ttaus,detap,detat
       common/cttaus/   tpro,zpro,ttar,ztar,ttaus,detap,detat
       parameter (mxcoos=60)
       common/cdensh/kdensh(matau,mxcoos,mxcoox,mxcooy),ktot(matau)
       character cy*3
 
       dlcoox=0.5
       dlcooy=0.5
       dlcoos=0.5
 
            if(ii.eq.0)then
 
       do i=1,mxcoos
       do j=1,mxcoox
       do k=1,mxcooy
       kdensh(itau,i,j,k)=0
       enddo
       enddo
       enddo
       ktot(itau)=0
 
            elseif(ii.eq.1)then
 
       if(itau.lt.1.or.itau.gt.mxtau)return
 
       tau=sngl(ttaus)
       z=s/tau
 
       do 1 i=1,mxcoos
       si=-mxcoos*dlcoos/2-dlcoos/2+i*dlcoos
       do 2 j=1,mxcoox
       xj=-mxcoox*dlcoox/2-dlcoox/2+j*dlcoox
       do 3 k=1,mxcooy
       yk=-mxcooy*dlcooy/2-dlcooy/2+k*dlcooy
       if(((x-xj)**2+(y-yk)**2+(z-si)**2).gt.rad**2)goto3
       kdensh(itau,i,j,k)=kdensh(itau,i,j,k)+1
       ktot(itau)=ktot(itau)+1
     3 continue
     2 continue
     1 continue
 
            elseif(ii.eq.2)then
 
       itaux=nint(xpar1)
       if(itaux.gt.mxtau)stop'STOP in xjden2: itaux too large'
 
       is=nint(xpar2/dlcoos)
       ix=nint(xpar3/dlcoox)
       iy=nint(xpar4/dlcooy)
       if(is.gt.mxcoos/2)stop'STOP in xjden2: s-range too large'
       if(ix.gt.mxcoox/2)stop'STOP in xjden2: x-range too large'
       if(iy.gt.mxcooy/2)stop'STOP in xjden2: y-range too large'
 
       do k=mxcooy/2+1-iy,mxcooy/2+iy
       write(cy,'(f3.1)')-mxcooy*dlcooy/2-dlcooy/2+k*dlcooy
       write(ifhi,'(a)')      'openhisto'
       write(ifhi,'(a,2f7.2)')'xrange',-is*dlcoos,is*dlcoos
       write(ifhi,'(a,2f7.2)')'yrange',-ix*dlcoox,ix*dlcoox
       write(ifhi,'(a)')      'set ityp2d 5'
       write(ifhi,'(a)') 'txt  "xaxis [z] "'
       write(ifhi,'(a)')
      *'txt  "yaxis  x (fm), y='//cy//' fm"'
       write(ifhi,'(a,i4)')   'array2d',2*is
       do j=mxcoox/2+1-ix,mxcoox/2+ix
       do i=mxcoos/2+1-is,mxcoos/2+is
       write(ifhi,'(e11.3)')
      *kdensh(itaux,i,j,k)/dlcooy/dlcoos/dlcoox/ktot(itaux)
       enddo
       enddo
       write(ifhi,'(a)')       '  endarray'
       write(ifhi,'(a)')       'closehisto plot2d'
       enddo
 
            else
 
       stop'STOP in xjden2: wrong option'
 
            endif
 
       return
       end
 
 cc-----------------------------------------------------------------------
 c      subroutine postscript(iii,ii,ic)
 cc-----------------------------------------------------------------------
 c      include 'epos.inc'
 c      character*10 color(10)
 c      if(iii.eq.0)then
 c      ifps=21
 c      open(unit=ifps,file='zzz.ps',status='unknown')
 c      WRITE(ifps,'(a)') '%!PS-Adobe-2.0'
 c      WRITE(ifps,'(a)') '%%Title: tt2.fig'
 c      WRITE(ifps,'(a)') '%%Orientation: Portrait'
 c      WRITE(ifps,'(a)') '%%BeginSetup'
 c      WRITE(ifps,'(a)') '%%IncludeFeature: *PageSize A4'
 c      WRITE(ifps,'(a)') '%%EndSetup'
 c      WRITE(ifps,'(a)') '%%EndComments'
 c      WRITE(ifps,*) '/l {lineto} bind def'
 c      WRITE(ifps,*) '/rl {rlineto} bind def'
 c      WRITE(ifps,*) '/m {moveto} bind def'
 c      WRITE(ifps,*) '/rm {rmoveto} bind def'
 c      WRITE(ifps,*) '/s {stroke} bind def'
 c      WRITE(ifps,*) '/gr {grestore} bind def'
 c      WRITE(ifps,*) '/gs {gsave} bind def'
 c      WRITE(ifps,*) '/cp {closepath} bind def'
 c      WRITE(ifps,*) '/tr {translate} bind def'
 c      WRITE(ifps,*) '/sc {scale} bind def'
 c      WRITE(ifps,*) '/sd {setdash} bind def'
 c      WRITE(ifps,*) '/sdo {[.01 .05] 0 sd} bind def'
 c      WRITE(ifps,*) '/sdf {[1 .0] 0 sd} bind def'
 c      WRITE(ifps,*) '/n {newpath} bind def'
 c      WRITE(ifps,*) '/slw {setlinewidth } bind def'
 c      write(ifps,*) '/srgb {setrgbcolor} bind def'
 c      write(ifps,*) '/lgrey      { 0 0.95 0.95 srgb} bind def'
 c      write(ifps,*) '/black      { 0 0 0 srgb} bind def'
 c      write(ifps,*) '/red        { 1 0 0 srgb} bind def  '
 c      write(ifps,*) '/green      { 0 1 0  srgb} bind def  '
 c      write(ifps,*) '/blue       { 0 0 1  srgb} bind def  '
 c      write(ifps,*) '/yellow     { 1 0.5 0  srgb} bind def  '
 c      write(ifps,*) '/turquoise  { 0 1 1  srgb} bind def  '
 c      write(ifps,*) '/purple     { 1 0 1  srgb} bind def  '
 cc      write(ifps,*) '/  {   srgb} bind def  '
 cc      write(ifps,*) '/  {   srgb} bind def  '
 c      write(ifps,*) '/ef {eofill} bind def'
 c      WRITE(ifps,'(a)') '%%EndProlog'
 c      WRITE(ifps,*) 'gsave'
 c      WRITE(ifps,*) '/Helvetica findfont 10 scalefont setfont'
 c      color(9)='lgrey     '
 c      color(1)='black     '
 c      color(2)='red       '
 c      color(3)='green     '
 c      color(4)='blue      '
 c      color(7)='yellow    '
 c      color(5)='turquoise '
 c      color(6)='purple    '
 c      np=0
 c         elseif(iii.eq.1)then
 c      np=np+1
 c      write(ifps,'(a,i4)') '%%Page: number ',np
 c      write(ifps,'(a)') 'gsave'
 c      WRITE(ifps,*) '100 700 tr'
 c      scale=0.125
 c      WRITE(ifps,*) 1./scale,1./scale,' sc'
 c      WRITE(ifps,*) scale/2.,' slw'
 c      WRITE(ifps,*) '/Helvetica findfont ',15.*scale
 c     &     ,' scalefont setfont'
 c      write(ifps,*) color(1),' n ',smin,xmin,' m ( tau:',tau,') show '
 c
 c      WRITE(ifps,*) '/Helvetica findfont ',5.*scale
 c     &     ,' scalefont setfont'
 c
 c
 c      yb=-2.
 c      dy=4./12.
 c      yb=yb-dy/2
 c      do iyb=0,11
 c        yb=yb+dy
 c        WRITE(ifps,*) 'gsave'
 c        WRITE(ifps,*) (xmax-xmin)*1.1*float(int(iyb/4))
 c     &       ,-(xmax-xmin)*1.1*mod(iyb,4),' tr'
 c        write(ifps,*) ' n ',smin,xmin,' m ',smax,xmin,' l '
 c     &       ,smax,xmax,' l ',smin,xmax,' l cp s '
 cc.......particles in layer iyb.............
 c        do i=1,nptl
 c          if(ii.gt.0)then
 c              write(ifps,*)  color(mod(i,5)+2)
 c     &             ,' n ',u,x-r,' m ',o,x-r,' l '
 c     &             ,o,x+r,' l ',u,x+r,' l cp s '
 c              write(ifps,*) ' n ',u,x-r,' m (',i,ior,') show '
 c          else
 c              write(ifps,*) ' n ',s,x,r,0,360,' arc ',color(ic),' s '
 c              write(ifps,*) ' n ',s-r,x,' m (',i,io,') show '
 c          endif
 c 10     enddo
 c        write(ifps,*) color(1),' n ',smin,xmin,' m (',yb,') show '
 c        WRITE(ifps,*) 'grestore'
 c      enddo                    !yb bin
 c      write(ifps,'(a)') 'grestore'
 c      write(ifps,*) 'showpage'
 c          elseif(iii.eq.2)then
 c       write(ifps,*) 'gr'
 c
 c       write(ifps,'(a)') '%%Trailer'
 c       write(ifps,'(a,i4)') '%%Pages: ',np
 c       write(ifps,'(a)') '%%EOF'
 c       close(unit=ifps)
 c          endif
 c
 c      return
 c      end
 c
 
 c------------------------------------------------------------------------------
       subroutine xtauev(iii)
 c------------------------------------------------------------------------------
       jdum=iii
       end
 c------------------------------------------------------------------------------
       subroutine wimi
 c------------------------------------------------------------------------------
       end
 c------------------------------------------------------------------------------
       subroutine wimino
 c------------------------------------------------------------------------------
       end
 c------------------------------------------------------------------------------
       subroutine xspace(iii)
 c------------------------------------------------------------------------------
       jdum=iii
       end
 c------------------------------------------------------------------------------
       subroutine wclu
 c------------------------------------------------------------------------------
       end
 c------------------------------------------------------------------------------
       subroutine wclufi
 c------------------------------------------------------------------------------
       end
 c------------------------------------------------------------------------------
       subroutine wtime(iii)
 c------------------------------------------------------------------------------
       jdum=iii
       end