      SUBROUTINE PYINDF(IP)
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
      INTEGER PYK,PYCHGE,PYCOMP
C...Commonblocks.
      COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
      COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
      COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
      SAVE /PYJETS/,/PYDAT1/,/PYDAT2/
C...Local arrays.
      DIMENSION DPS(5),PSI(4),NFI(3),NFL(3),IFET(3),KFLF(3),
     &KFLO(2),PXO(2),PYO(2),WO(2)
C.. MOPS error message
      IF(MSTJ(12).GT.3) CALL PYERRM(9,'(PYINDF:) MSTJ(12)>3 options'//
     &' are not treated as expected in independent fragmentation')
C...Reset counters. Identify parton system and take copy. Check flavour.
      NSAV=N
      MSTU90=MSTU(90)
      NJET=0
      KQSUM=0
      DO 100 J=1,5
        DPS(J)=0D0
  100 CONTINUE
      I=IP-1
  110 I=I+1
      IF(I.GT.MIN(N,MSTU(4)-MSTU(32))) THEN
        CALL PYERRM(12,'(PYINDF:) failed to reconstruct jet system')
        IF(MSTU(21).GE.1) RETURN
      ENDIF
      IF(K(I,1).NE.1.AND.K(I,1).NE.2) GOTO 110
      KC=PYCOMP(K(I,2))
      IF(KC.EQ.0) GOTO 110
      KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
      IF(KQ.EQ.0) GOTO 110
      NJET=NJET+1
      IF(KQ.NE.2) KQSUM=KQSUM+KQ
      DO 120 J=1,5
        K(NSAV+NJET,J)=K(I,J)
        P(NSAV+NJET,J)=P(I,J)
        DPS(J)=DPS(J)+P(I,J)
  120 CONTINUE
      K(NSAV+NJET,3)=I
      IF(K(I,1).EQ.2.OR.(MSTJ(3).LE.5.AND.N.GT.I.AND.
     &K(I+1,1).EQ.2)) GOTO 110
      IF(NJET.NE.1.AND.KQSUM.NE.0) THEN
        CALL PYERRM(12,'(PYINDF:) unphysical flavour combination')
        IF(MSTU(21).GE.1) RETURN
      ENDIF
C...Boost copied system to CM frame. Find CM energy and sum flavours.
      IF(NJET.NE.1) THEN
        MSTU(33)=1
        CALL PYROBO(NSAV+1,NSAV+NJET,0D0,0D0,-DPS(1)/DPS(4),
     &  -DPS(2)/DPS(4),-DPS(3)/DPS(4))
      ENDIF
      PECM=0D0
      DO 130 J=1,3
        NFI(J)=0
  130 CONTINUE
      DO 140 I=NSAV+1,NSAV+NJET
        PECM=PECM+P(I,4)
        KFA=IABS(K(I,2))
        IF(KFA.LE.3) THEN
          NFI(KFA)=NFI(KFA)+ISIGN(1,K(I,2))
        ELSEIF(KFA.GT.1000) THEN
          KFLA=MOD(KFA/1000,10)
          KFLB=MOD(KFA/100,10)
          IF(KFLA.LE.3) NFI(KFLA)=NFI(KFLA)+ISIGN(1,K(I,2))
          IF(KFLB.LE.3) NFI(KFLB)=NFI(KFLB)+ISIGN(1,K(I,2))
        ENDIF
  140 CONTINUE
C...Loop over attempts made. Reset counters.
      NTRY=0
  150 NTRY=NTRY+1
      IF(NTRY.GT.200) THEN
        CALL PYERRM(14,'(PYINDF:) caught in infinite loop')
        IF(MSTU(21).GE.1) RETURN
      ENDIF
      N=NSAV+NJET
      MSTU(90)=MSTU90
      DO 160 J=1,3
        NFL(J)=NFI(J)
        IFET(J)=0
        KFLF(J)=0
  160 CONTINUE
C...Loop over jets to be fragmented.
      DO 230 IP1=NSAV+1,NSAV+NJET
        MSTJ(91)=0
        NSAV1=N
        MSTU91=MSTU(90)
C...Initial flavour and momentum values. Jet along +z axis.
        KFLH=IABS(K(IP1,2))
        IF(KFLH.GT.10) KFLH=MOD(KFLH/1000,10)
        KFLO(2)=0
        WF=P(IP1,4)+SQRT(P(IP1,1)**2+P(IP1,2)**2+P(IP1,3)**2)
C...Initial values for quark or diquark jet.
  170   IF(IABS(K(IP1,2)).NE.21) THEN
          NSTR=1
          KFLO(1)=K(IP1,2)
          CALL PYPTDI(0,PXO(1),PYO(1))
          WO(1)=WF
C...Initial values for gluon treated like random quark jet.
        ELSEIF(MSTJ(2).LE.2) THEN
          NSTR=1
          IF(MSTJ(2).EQ.2) MSTJ(91)=1
          KFLO(1)=INT(1D0+(2D0+PARJ(2))*PYR(0))*(-1)**INT(PYR(0)+0.5D0)
          CALL PYPTDI(0,PXO(1),PYO(1))
          WO(1)=WF
C...Initial values for gluon treated like quark-antiquark jet pair,
C...sharing energy according to Altarelli-Parisi splitting function.
        ELSE
          NSTR=2
          IF(MSTJ(2).EQ.4) MSTJ(91)=1
          KFLO(1)=INT(1D0+(2D0+PARJ(2))*PYR(0))*(-1)**INT(PYR(0)+0.5D0)
          KFLO(2)=-KFLO(1)
          CALL PYPTDI(0,PXO(1),PYO(1))
          PXO(2)=-PXO(1)
          PYO(2)=-PYO(1)
          WO(1)=WF*PYR(0)**(1D0/3D0)
          WO(2)=WF-WO(1)
        ENDIF
C...Initial values for rank, flavour, pT and W+.
        DO 220 ISTR=1,NSTR
  180     I=N
          MSTU(90)=MSTU91
          IRANK=0
          KFL1=KFLO(ISTR)
          PX1=PXO(ISTR)
          PY1=PYO(ISTR)
          W=WO(ISTR)
C...New hadron. Generate flavour and hadron species.
  190     I=I+1
          IF(I.GE.MSTU(4)-MSTU(32)-NJET-5) THEN
            CALL PYERRM(11,'(PYINDF:) no more memory left in PYJETS')
            IF(MSTU(21).GE.1) RETURN
          ENDIF
          IRANK=IRANK+1
          K(I,1)=1
          K(I,3)=IP1
          K(I,4)=0
          K(I,5)=0
  200     CALL PYKFDI(KFL1,0,KFL2,K(I,2))
          IF(K(I,2).EQ.0) GOTO 180
          IF(IRANK.EQ.1.AND.IABS(KFL1).LE.10.AND.IABS(KFL2).GT.10) THEN
            IF(PYR(0).GT.PARJ(19)) GOTO 200
          ENDIF
C...Find hadron mass. Generate four-momentum.
          P(I,5)=PYMASS(K(I,2))
          CALL PYPTDI(KFL1,PX2,PY2)
          P(I,1)=PX1+PX2
          P(I,2)=PY1+PY2
          PR=P(I,5)**2+P(I,1)**2+P(I,2)**2
          CALL PYZDIS(KFL1,KFL2,PR,Z)
          MZSAV=0
          IF(IABS(KFL1).GE.4.AND.IABS(KFL1).LE.8.AND.MSTU(90).LT.8) THEN
            MZSAV=1
            MSTU(90)=MSTU(90)+1
            MSTU(90+MSTU(90))=I
            PARU(90+MSTU(90))=Z
          ENDIF
          P(I,3)=0.5D0*(Z*W-PR/MAX(1D-4,Z*W))
          P(I,4)=0.5D0*(Z*W+PR/MAX(1D-4,Z*W))
          IF(MSTJ(3).GE.1.AND.IRANK.EQ.1.AND.KFLH.GE.4.AND.
     &    P(I,3).LE.0.001D0) THEN
            IF(W.GE.P(I,5)+0.5D0*PARJ(32)) GOTO 180
            P(I,3)=0.0001D0
            P(I,4)=SQRT(PR)
            Z=P(I,4)/W
          ENDIF
C...Remaining flavour and momentum.
          KFL1=-KFL2
          PX1=-PX2
          PY1=-PY2
          W=(1D0-Z)*W
          DO 210 J=1,5
            V(I,J)=0D0
  210     CONTINUE
C...Check if pL acceptable. Go back for new hadron if enough energy.
          IF(MSTJ(3).GE.0.AND.P(I,3).LT.0D0) THEN
            I=I-1
            IF(MZSAV.EQ.1) MSTU(90)=MSTU(90)-1
          ENDIF
          IF(W.GT.PARJ(31)) GOTO 190
          N=I
  220   CONTINUE
        IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) WF=WF+0.1D0*PARJ(32)
        IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) GOTO 170
C...Rotate jet to new direction.
        THE=PYANGL(P(IP1,3),SQRT(P(IP1,1)**2+P(IP1,2)**2))
        PHI=PYANGL(P(IP1,1),P(IP1,2))
        MSTU(33)=1
        CALL PYROBO(NSAV1+1,N,THE,PHI,0D0,0D0,0D0)
        K(K(IP1,3),4)=NSAV1+1
        K(K(IP1,3),5)=N
C...End of jet generation loop. Skip conservation in some cases.
  230 CONTINUE
      IF(NJET.EQ.1.OR.MSTJ(3).LE.0) GOTO 490
      IF(MOD(MSTJ(3),5).NE.0.AND.N-NSAV-NJET.LT.2) GOTO 150
C...Subtract off produced hadron flavours, finished if zero.
      DO 240 I=NSAV+NJET+1,N
        KFA=IABS(K(I,2))
        KFLA=MOD(KFA/1000,10)
        KFLB=MOD(KFA/100,10)
        KFLC=MOD(KFA/10,10)
        IF(KFLA.EQ.0) THEN
          IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))*(-1)**KFLB
          IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(I,2))*(-1)**KFLB
        ELSE
          IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)-ISIGN(1,K(I,2))
          IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))
          IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISIGN(1,K(I,2))
        ENDIF
  240 CONTINUE
      NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
     &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
      IF(NREQ.EQ.0) GOTO 320
C...Take away flavour of low-momentum particles until enough freedom.
      NREM=0
  250 IREM=0
      P2MIN=PECM**2
      DO 260 I=NSAV+NJET+1,N
        P2=P(I,1)**2+P(I,2)**2+P(I,3)**2
        IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) IREM=I
        IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) P2MIN=P2
  260 CONTINUE
      IF(IREM.EQ.0) GOTO 150
      K(IREM,1)=7
      KFA=IABS(K(IREM,2))
      KFLA=MOD(KFA/1000,10)
      KFLB=MOD(KFA/100,10)
      KFLC=MOD(KFA/10,10)
      IF(KFLA.GE.4.OR.KFLB.GE.4) K(IREM,1)=8
      IF(K(IREM,1).EQ.8) GOTO 250
      IF(KFLA.EQ.0) THEN
        ISGN=ISIGN(1,K(IREM,2))*(-1)**KFLB
        IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISGN
        IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISGN
      ELSE
        IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)+ISIGN(1,K(IREM,2))
        IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISIGN(1,K(IREM,2))
        IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(IREM,2))
      ENDIF
      NREM=NREM+1
      NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
     &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
      IF(NREQ.GT.NREM) GOTO 250
      DO 270 I=NSAV+NJET+1,N
        IF(K(I,1).EQ.8) K(I,1)=1
  270 CONTINUE
C...Find combination of existing and new flavours for hadron.
  280 NFET=2
      IF(NFL(1)+NFL(2)+NFL(3).NE.0) NFET=3
      IF(NREQ.LT.NREM) NFET=1
      IF(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)).EQ.0) NFET=0
      DO 290 J=1,NFET
        IFET(J)=1+(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)))*PYR(0)
        KFLF(J)=ISIGN(1,NFL(1))
        IF(IFET(J).GT.IABS(NFL(1))) KFLF(J)=ISIGN(2,NFL(2))
        IF(IFET(J).GT.IABS(NFL(1))+IABS(NFL(2))) KFLF(J)=ISIGN(3,NFL(3))
  290 CONTINUE
      IF(NFET.EQ.2.AND.(IFET(1).EQ.IFET(2).OR.KFLF(1)*KFLF(2).GT.0))
     &GOTO 280
      IF(NFET.EQ.3.AND.(IFET(1).EQ.IFET(2).OR.IFET(1).EQ.IFET(3).OR.
     &IFET(2).EQ.IFET(3).OR.KFLF(1)*KFLF(2).LT.0.OR.KFLF(1)*KFLF(3)
     &.LT.0.OR.KFLF(1)*(NFL(1)+NFL(2)+NFL(3)).LT.0)) GOTO 280
      IF(NFET.EQ.0) KFLF(1)=1+INT((2D0+PARJ(2))*PYR(0))
      IF(NFET.EQ.0) KFLF(2)=-KFLF(1)
      IF(NFET.EQ.1) KFLF(2)=ISIGN(1+INT((2D0+PARJ(2))*PYR(0)),-KFLF(1))
      IF(NFET.LE.2) KFLF(3)=0
      IF(KFLF(3).NE.0) THEN
        KFLFC=ISIGN(1000*MAX(IABS(KFLF(1)),IABS(KFLF(3)))+
     &  100*MIN(IABS(KFLF(1)),IABS(KFLF(3)))+1,KFLF(1))
        IF(KFLF(1).EQ.KFLF(3).OR.(1D0+3D0*PARJ(4))*PYR(0).GT.1D0)
     &  KFLFC=KFLFC+ISIGN(2,KFLFC)
      ELSE
        KFLFC=KFLF(1)
      ENDIF
      CALL PYKFDI(KFLFC,KFLF(2),KFLDMP,KF)
      IF(KF.EQ.0) GOTO 280
      DO 300 J=1,MAX(2,NFET)
        NFL(IABS(KFLF(J)))=NFL(IABS(KFLF(J)))-ISIGN(1,KFLF(J))
  300 CONTINUE
C...Store hadron at random among free positions.
      NPOS=MIN(1+INT(PYR(0)*NREM),NREM)
      DO 310 I=NSAV+NJET+1,N
        IF(K(I,1).EQ.7) NPOS=NPOS-1
        IF(K(I,1).EQ.1.OR.NPOS.NE.0) GOTO 310
        K(I,1)=1
        K(I,2)=KF
        P(I,5)=PYMASS(K(I,2))
        P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
  310 CONTINUE
      NREM=NREM-1
      NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
     &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
      IF(NREM.GT.0) GOTO 280
C...Compensate for missing momentum in global scheme (3 options).
  320 IF(MOD(MSTJ(3),5).NE.0.AND.MOD(MSTJ(3),5).NE.4) THEN
        DO 340 J=1,3
          PSI(J)=0D0
          DO 330 I=NSAV+NJET+1,N
            PSI(J)=PSI(J)+P(I,J)
  330     CONTINUE
  340   CONTINUE
        PSI(4)=PSI(1)**2+PSI(2)**2+PSI(3)**2
        PWS=0D0
        DO 350 I=NSAV+NJET+1,N
          IF(MOD(MSTJ(3),5).EQ.1) PWS=PWS+P(I,4)
          IF(MOD(MSTJ(3),5).EQ.2) PWS=PWS+SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
     &    PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
          IF(MOD(MSTJ(3),5).EQ.3) PWS=PWS+1D0
  350   CONTINUE
        DO 370 I=NSAV+NJET+1,N
          IF(MOD(MSTJ(3),5).EQ.1) PW=P(I,4)
          IF(MOD(MSTJ(3),5).EQ.2) PW=SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
     &    PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
          IF(MOD(MSTJ(3),5).EQ.3) PW=1D0
          DO 360 J=1,3
            P(I,J)=P(I,J)-PSI(J)*PW/PWS
  360     CONTINUE
          P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
  370   CONTINUE
C...Compensate for missing momentum withing each jet separately.
      ELSEIF(MOD(MSTJ(3),5).EQ.4) THEN
        DO 390 I=N+1,N+NJET
          K(I,1)=0
          DO 380 J=1,5
            P(I,J)=0D0
  380     CONTINUE
  390   CONTINUE
        DO 410 I=NSAV+NJET+1,N
          IR1=K(I,3)
          IR2=N+IR1-NSAV
          K(IR2,1)=K(IR2,1)+1
          PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
     &    (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
          DO 400 J=1,3
            P(IR2,J)=P(IR2,J)+P(I,J)-PLS*P(IR1,J)
  400     CONTINUE
          P(IR2,4)=P(IR2,4)+P(I,4)
          P(IR2,5)=P(IR2,5)+PLS
  410   CONTINUE
        PSS=0D0
        DO 420 I=N+1,N+NJET
          IF(K(I,1).NE.0) PSS=PSS+P(I,4)/(PECM*(0.8D0*P(I,5)+0.2D0))
  420   CONTINUE
        DO 440 I=NSAV+NJET+1,N
          IR1=K(I,3)
          IR2=N+IR1-NSAV
          PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
     &    (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
          DO 430 J=1,3
            P(I,J)=P(I,J)-P(IR2,J)/K(IR2,1)+(1D0/(P(IR2,5)*PSS)-1D0)*
     &      PLS*P(IR1,J)
  430     CONTINUE
          P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
  440   CONTINUE
      ENDIF
C...Scale momenta for energy conservation.
      IF(MOD(MSTJ(3),5).NE.0) THEN
        PMS=0D0
        PES=0D0
        PQS=0D0
        DO 450 I=NSAV+NJET+1,N
          PMS=PMS+P(I,5)
          PES=PES+P(I,4)
          PQS=PQS+P(I,5)**2/P(I,4)
  450   CONTINUE
        IF(PMS.GE.PECM) GOTO 150
        NECO=0
  460   NECO=NECO+1
        PFAC=(PECM-PQS)/(PES-PQS)
        PES=0D0
        PQS=0D0
        DO 480 I=NSAV+NJET+1,N
          DO 470 J=1,3
            P(I,J)=PFAC*P(I,J)
  470     CONTINUE
          P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
          PES=PES+P(I,4)
          PQS=PQS+P(I,5)**2/P(I,4)
  480   CONTINUE
        IF(NECO.LT.10.AND.ABS(PECM-PES).GT.2D-6*PECM) GOTO 460
      ENDIF
C...Origin of produced particles and parton daughter pointers.
  490 DO 500 I=NSAV+NJET+1,N
        IF(MSTU(16).NE.2) K(I,3)=NSAV+1
        IF(MSTU(16).EQ.2) K(I,3)=K(K(I,3),3)
  500 CONTINUE
      DO 510 I=NSAV+1,NSAV+NJET
        I1=K(I,3)
        K(I1,1)=K(I1,1)+10
        IF(MSTU(16).NE.2) THEN
          K(I1,4)=NSAV+1
          K(I1,5)=NSAV+1
        ELSE
          K(I1,4)=K(I1,4)-NJET+1
          K(I1,5)=K(I1,5)-NJET+1
          IF(K(I1,5).LT.K(I1,4)) THEN
            K(I1,4)=0
            K(I1,5)=0
          ENDIF
        ENDIF
  510 CONTINUE
C...Document independent fragmentation system. Remove copy of jets.
      NSAV=NSAV+1
      K(NSAV,1)=11
      K(NSAV,2)=93
      K(NSAV,3)=IP
      K(NSAV,4)=NSAV+1
      K(NSAV,5)=N-NJET+1
      DO 520 J=1,4
        P(NSAV,J)=DPS(J)
        V(NSAV,J)=V(IP,J)
  520 CONTINUE
      P(NSAV,5)=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))
      V(NSAV,5)=0D0
      DO 540 I=NSAV+NJET,N
        DO 530 J=1,5
          K(I-NJET+1,J)=K(I,J)
          P(I-NJET+1,J)=P(I,J)
          V(I-NJET+1,J)=V(I,J)
  530   CONTINUE
  540 CONTINUE
      N=N-NJET+1
      DO 550 IZ=MSTU90+1,MSTU(90)
        MSTU(90+IZ)=MSTU(90+IZ)-NJET+1
  550 CONTINUE
C...Boost back particle system. Set production vertices.
      IF(NJET.NE.1) CALL PYROBO(NSAV+1,N,0D0,0D0,DPS(1)/DPS(4),
     &DPS(2)/DPS(4),DPS(3)/DPS(4))
      DO 570 I=NSAV+1,N
        DO 560 J=1,4
          V(I,J)=V(IP,J)
  560   CONTINUE
  570 CONTINUE
      RETURN
      END
