HIJING/hipyset1_35/luindf_hijing.F

   1 * $Id$
   2
   3 C*********************************************************************
   4
   5       SUBROUTINE LUINDF_HIJING(IP)
   6
   7 C...Purpose: to handle the fragmentation of a jet system (or a single
   8 C...jet) according to independent fragmentation models.
   9       IMPLICIT DOUBLE PRECISION(D)
  10 #include "lujets_hijing.inc"
  11 #include "ludat1_hijing.inc"
  12 #include "ludat2_hijing.inc"
  13       DIMENSION DPS(5),PSI(4),NFI(3),NFL(3),IFET(3),KFLF(3),
  14      &KFLO(2),PXO(2),PYO(2),WO(2)
  15
  16 C...Reset counters. Identify parton system and take copy. Check flavour.
  17       NSAV=N
  18       NJET=0
  19       KQSUM=0
  20       DO 100 J=1,5
  21   100 DPS(J)=0.
  22       I=IP-1
  23   110 I=I+1
  24       IF(I.GT.MIN(N,MSTU(4)-MSTU(32))) THEN
  25          CALL LUERRM_HIJING(12
  26      $        ,'(LUINDF_HIJING:) failed to reconstruct jet system')
  27         IF(MSTU(21).GE.1) RETURN
  28       ENDIF
  29       IF(K(I,1).NE.1.AND.K(I,1).NE.2) GOTO 110
  30       KC=LUCOMP_HIJING(K(I,2))
  31       IF(KC.EQ.0) GOTO 110
  32       KQ=KCHG(KC,2)*ISIGN(1,K(I,2))
  33       IF(KQ.EQ.0) GOTO 110
  34       NJET=NJET+1
  35       IF(KQ.NE.2) KQSUM=KQSUM+KQ
  36       DO 120 J=1,5
  37       K(NSAV+NJET,J)=K(I,J)
  38       P(NSAV+NJET,J)=P(I,J)
  39   120 DPS(J)=DPS(J)+P(I,J)
  40       K(NSAV+NJET,3)=I
  41       IF(K(I,1).EQ.2.OR.(MSTJ(3).LE.5.AND.N.GT.I.AND.
  42      &K(I+1,1).EQ.2)) GOTO 110
  43       IF(NJET.NE.1.AND.KQSUM.NE.0) THEN
  44          CALL LUERRM_HIJING(12
  45      $        ,'(LUINDF_HIJING:) unphysical flavour combination')
  46         IF(MSTU(21).GE.1) RETURN
  47       ENDIF
  48
  49 C...Boost copied system to CM frame. Find CM energy and sum flavours.
  50       IF(NJET.NE.1) CALL LUDBRB_HIJING(NSAV+1,NSAV+NJET,0.,0.,-DPS(1)
  51      $     /DPS(4),-DPS(2)/DPS(4),-DPS(3)/DPS(4))
  52       PECM=0.
  53       DO 130 J=1,3
  54   130 NFI(J)=0
  55       DO 140 I=NSAV+1,NSAV+NJET
  56       PECM=PECM+P(I,4)
  57       KFA=IABS(K(I,2))
  58       IF(KFA.LE.3) THEN
  59         NFI(KFA)=NFI(KFA)+ISIGN(1,K(I,2))
  60       ELSEIF(KFA.GT.1000) THEN
  61         KFLA=MOD(KFA/1000,10)
  62         KFLB=MOD(KFA/100,10)
  63         IF(KFLA.LE.3) NFI(KFLA)=NFI(KFLA)+ISIGN(1,K(I,2))
  64         IF(KFLB.LE.3) NFI(KFLB)=NFI(KFLB)+ISIGN(1,K(I,2))
  65       ENDIF
  66   140 CONTINUE
  67
  68 C...Loop over attempts made. Reset counters.
  69       NTRY=0
  70   150 NTRY=NTRY+1
  71       N=NSAV+NJET
  72       IF(NTRY.GT.200) THEN
  73          CALL LUERRM_HIJING(14
  74      $        ,'(LUINDF_HIJING:) caught in infinite loop')
  75         IF(MSTU(21).GE.1) RETURN
  76       ENDIF
  77       DO 160 J=1,3
  78       NFL(J)=NFI(J)
  79       IFET(J)=0
  80   160 KFLF(J)=0
  81
  82 C...Loop over jets to be fragmented.
  83       DO 230 IP1=NSAV+1,NSAV+NJET
  84       MSTJ(91)=0
  85       NSAV1=N
  86
  87 C...Initial flavour and momentum values. Jet along +z axis.
  88       KFLH=IABS(K(IP1,2))
  89       IF(KFLH.GT.10) KFLH=MOD(KFLH/1000,10)
  90       KFLO(2)=0
  91       WF=P(IP1,4)+SQRT(P(IP1,1)**2+P(IP1,2)**2+P(IP1,3)**2)
  92
  93 C...Initial values for quark or diquark jet.
  94   170 IF(IABS(K(IP1,2)).NE.21) THEN
  95         NSTR=1
  96         KFLO(1)=K(IP1,2)
  97         CALL LUPTDI_HIJING(0,PXO(1),PYO(1))
  98         WO(1)=WF
  99
 100 C...Initial values for gluon treated like random quark jet.
 101       ELSEIF(MSTJ(2).LE.2) THEN
 102         NSTR=1
 103         IF(MSTJ(2).EQ.2) MSTJ(91)=1
 104         KFLO(1)=INT(1.+(2.+PARJ(2))*RLU_HIJING(0))*(-1)
 105      $       **INT(RLU_HIJING(0)+0.5)
 106         CALL LUPTDI_HIJING(0,PXO(1),PYO(1))
 107         WO(1)=WF
 108
 109 C...Initial values for gluon treated like quark-antiquark jet pair,
 110 C...sharing energy according to Altarelli-Parisi splitting function.
 111       ELSE
 112         NSTR=2
 113         IF(MSTJ(2).EQ.4) MSTJ(91)=1
 114         KFLO(1)=INT(1.+(2.+PARJ(2))*RLU_HIJING(0))*(-1)
 115      $       **INT(RLU_HIJING(0)+0.5)
 116         KFLO(2)=-KFLO(1)
 117         CALL LUPTDI_HIJING(0,PXO(1),PYO(1))
 118         PXO(2)=-PXO(1)
 119         PYO(2)=-PYO(1)
 120  111    WO(1)=WF*RLU_HIJING(0)**(1./3.)
 121         WO(2)=WF-WO(1)
 122         IF (ABS(WO(2)).LE.0..OR.ABS(WO(1)).LE.0.) GOTO 111
 123       ENDIF
 124
 125 C...Initial values for rank, flavour, pT and W+.
 126       DO 220 ISTR=1,NSTR
 127   180 I=N
 128       IRANK=0
 129       KFL1=KFLO(ISTR)
 130       PX1=PXO(ISTR)
 131       PY1=PYO(ISTR)
 132       W=WO(ISTR)
 133
 134 C...New hadron. Generate flavour and hadron species.
 135   190 I=I+1
 136       IF(I.GE.MSTU(4)-MSTU(32)-NJET-5) THEN
 137          CALL LUERRM_HIJING(11
 138      $        ,'(LUINDF_HIJING:) no more memory left in LUJETS_HIJING')
 139         IF(MSTU(21).GE.1) RETURN
 140       ENDIF
 141       IRANK=IRANK+1
 142       K(I,1)=1
 143       K(I,3)=IP1
 144       K(I,4)=0
 145       K(I,5)=0
 146   200 CALL LUKFDI_HIJING(KFL1,0,KFL2,K(I,2))
 147       IF(K(I,2).EQ.0) GOTO 180
 148       IF(MSTJ(12).GE.3.AND.IRANK.EQ.1.AND.IABS(KFL1).LE.10.AND.
 149      &IABS(KFL2).GT.10) THEN
 150         IF(RLU_HIJING(0).GT.PARJ(19)) GOTO 200
 151       ENDIF
 152
 153 C...Find hadron mass. Generate four-momentum.
 154       P(I,5)=ULMASS_HIJING(K(I,2))
 155       CALL LUPTDI_HIJING(KFL1,PX2,PY2)
 156       P(I,1)=PX1+PX2
 157       P(I,2)=PY1+PY2
 158       PR=P(I,5)**2+P(I,1)**2+P(I,2)**2
 159       CALL LUZDIS_HIJING(KFL1,KFL2,PR,Z)
 160       P(I,3)=0.5*(Z*W-PR/(Z*W))
 161       P(I,4)=0.5*(Z*W+PR/(Z*W))
 162       IF(MSTJ(3).GE.1.AND.IRANK.EQ.1.AND.KFLH.GE.4.AND.
 163      &P(I,3).LE.0.001) THEN
 164         IF(W.GE.P(I,5)+0.5*PARJ(32)) GOTO 180
 165         P(I,3)=0.0001
 166         P(I,4)=SQRT(PR)
 167         Z=P(I,4)/W
 168       ENDIF
 169
 170 C...Remaining flavour and momentum.
 171       KFL1=-KFL2
 172       PX1=-PX2
 173       PY1=-PY2
 174       W=(1.-Z)*W
 175       DO 210 J=1,5
 176   210 V(I,J)=0.
 177
 178 C...Check if pL acceptable. Go back for new hadron if enough energy.
 179       IF(MSTJ(3).GE.0.AND.P(I,3).LT.0.) I=I-1
 180       IF(W.GT.PARJ(31)) GOTO 190
 181   220 N=I
 182       IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) WF=WF+0.1*PARJ(32)
 183       IF(MOD(MSTJ(3),5).EQ.4.AND.N.EQ.NSAV1) GOTO 170
 184
 185 C...Rotate jet to new direction.
 186       THE=ULANGL_HIJING(P(IP1,3),SQRT(P(IP1,1)**2+P(IP1,2)**2))
 187       PHI=ULANGL_HIJING(P(IP1,1),P(IP1,2))
 188       CALL LUDBRB_HIJING(NSAV1+1,N,THE,PHI,0D0,0D0,0D0)
 189       K(K(IP1,3),4)=NSAV1+1
 190       K(K(IP1,3),5)=N
 191
 192 C...End of jet generation loop. Skip conservation in some cases.
 193   230 CONTINUE
 194       IF(NJET.EQ.1.OR.MSTJ(3).LE.0) GOTO 470
 195       IF(MOD(MSTJ(3),5).NE.0.AND.N-NSAV-NJET.LT.2) GOTO 150
 196
 197 C...Subtract off produced hadron flavours, finished if zero.
 198       DO 240 I=NSAV+NJET+1,N
 199       KFA=IABS(K(I,2))
 200       KFLA=MOD(KFA/1000,10)
 201       KFLB=MOD(KFA/100,10)
 202       KFLC=MOD(KFA/10,10)
 203       IF(KFLA.EQ.0) THEN
 204         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))*(-1)**KFLB
 205         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(I,2))*(-1)**KFLB
 206       ELSE
 207         IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)-ISIGN(1,K(I,2))
 208         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)-ISIGN(1,K(I,2))
 209         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISIGN(1,K(I,2))
 210       ENDIF
 211   240 CONTINUE
 212       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
 213      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
 214       IF(NREQ.EQ.0) GOTO 320
 215
 216 C...Take away flavour of low-momentum particles until enough freedom.
 217       NREM=0
 218   250 IREM=0
 219       P2MIN=PECM**2
 220       DO 260 I=NSAV+NJET+1,N
 221       P2=P(I,1)**2+P(I,2)**2+P(I,3)**2
 222       IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) IREM=I
 223   260 IF(K(I,1).EQ.1.AND.P2.LT.P2MIN) P2MIN=P2
 224       IF(IREM.EQ.0) GOTO 150
 225       K(IREM,1)=7
 226       KFA=IABS(K(IREM,2))
 227       KFLA=MOD(KFA/1000,10)
 228       KFLB=MOD(KFA/100,10)
 229       KFLC=MOD(KFA/10,10)
 230       IF(KFLA.GE.4.OR.KFLB.GE.4) K(IREM,1)=8
 231       IF(K(IREM,1).EQ.8) GOTO 250
 232       IF(KFLA.EQ.0) THEN
 233         ISGN=ISIGN(1,K(IREM,2))*(-1)**KFLB
 234         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISGN
 235         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)-ISGN
 236       ELSE
 237         IF(KFLA.LE.3) NFL(KFLA)=NFL(KFLA)+ISIGN(1,K(IREM,2))
 238         IF(KFLB.LE.3) NFL(KFLB)=NFL(KFLB)+ISIGN(1,K(IREM,2))
 239         IF(KFLC.LE.3) NFL(KFLC)=NFL(KFLC)+ISIGN(1,K(IREM,2))
 240       ENDIF
 241       NREM=NREM+1
 242       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
 243      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
 244       IF(NREQ.GT.NREM) GOTO 250
 245       DO 270 I=NSAV+NJET+1,N
 246   270 IF(K(I,1).EQ.8) K(I,1)=1
 247
 248 C...Find combination of existing and new flavours for hadron.
 249   280 NFET=2
 250       IF(NFL(1)+NFL(2)+NFL(3).NE.0) NFET=3
 251       IF(NREQ.LT.NREM) NFET=1
 252       IF(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)).EQ.0) NFET=0
 253       DO 290 J=1,NFET
 254       IFET(J)=1+(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3)))*RLU_HIJING(0)
 255       KFLF(J)=ISIGN(1,NFL(1))
 256       IF(IFET(J).GT.IABS(NFL(1))) KFLF(J)=ISIGN(2,NFL(2))
 257   290 IF(IFET(J).GT.IABS(NFL(1))+IABS(NFL(2))) KFLF(J)=ISIGN(3,NFL(3))
 258       IF(NFET.EQ.2.AND.(IFET(1).EQ.IFET(2).OR.KFLF(1)*KFLF(2).GT.0))
 259      &GOTO 280
 260       IF(NFET.EQ.3.AND.(IFET(1).EQ.IFET(2).OR.IFET(1).EQ.IFET(3).OR.
 261      &IFET(2).EQ.IFET(3).OR.KFLF(1)*KFLF(2).LT.0.OR.KFLF(1)*KFLF(3).
 262      &LT.0.OR.KFLF(1)*(NFL(1)+NFL(2)+NFL(3)).LT.0)) GOTO 280
 263       IF(NFET.EQ.0) KFLF(1)=1+INT((2.+PARJ(2))*RLU_HIJING(0))
 264       IF(NFET.EQ.0) KFLF(2)=-KFLF(1)
 265       IF(NFET.EQ.1) KFLF(2)=ISIGN(1+INT((2.+PARJ(2))*RLU_HIJING(0)),
 266      $     -KFLF(1))
 267       IF(NFET.LE.2) KFLF(3)=0
 268       IF(KFLF(3).NE.0) THEN
 269         KFLFC=ISIGN(1000*MAX(IABS(KFLF(1)),IABS(KFLF(3)))+
 270      &  100*MIN(IABS(KFLF(1)),IABS(KFLF(3)))+1,KFLF(1))
 271         IF(KFLF(1).EQ.KFLF(3).OR.(1.+3.*PARJ(4))*RLU_HIJING(0).GT.1.)
 272      &  KFLFC=KFLFC+ISIGN(2,KFLFC)
 273       ELSE
 274         KFLFC=KFLF(1)
 275       ENDIF
 276       CALL LUKFDI_HIJING(KFLFC,KFLF(2),KFLDMP,KF)
 277       IF(KF.EQ.0) GOTO 280
 278       DO 300 J=1,MAX(2,NFET)
 279   300 NFL(IABS(KFLF(J)))=NFL(IABS(KFLF(J)))-ISIGN(1,KFLF(J))
 280
 281 C...Store hadron at random among free positions.
 282       NPOS=MIN(1+INT(RLU_HIJING(0)*NREM),NREM)
 283       DO 310 I=NSAV+NJET+1,N
 284       IF(K(I,1).EQ.7) NPOS=NPOS-1
 285       IF(K(I,1).EQ.1.OR.NPOS.NE.0) GOTO 310
 286       K(I,1)=1
 287       K(I,2)=KF
 288       P(I,5)=ULMASS_HIJING(K(I,2))
 289       P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
 290   310 CONTINUE
 291       NREM=NREM-1
 292       NREQ=(IABS(NFL(1))+IABS(NFL(2))+IABS(NFL(3))-IABS(NFL(1)+
 293      &NFL(2)+NFL(3)))/2+IABS(NFL(1)+NFL(2)+NFL(3))/3
 294       IF(NREM.GT.0) GOTO 280
 295
 296 C...Compensate for missing momentum in global scheme (3 options).
 297   320 IF(MOD(MSTJ(3),5).NE.0.AND.MOD(MSTJ(3),5).NE.4) THEN
 298         DO 330 J=1,3
 299         PSI(J)=0.
 300         DO 330 I=NSAV+NJET+1,N
 301   330   PSI(J)=PSI(J)+P(I,J)
 302         PSI(4)=PSI(1)**2+PSI(2)**2+PSI(3)**2
 303         PWS=0.
 304         DO 340 I=NSAV+NJET+1,N
 305         IF(MOD(MSTJ(3),5).EQ.1) PWS=PWS+P(I,4)
 306         IF(MOD(MSTJ(3),5).EQ.2) PWS=PWS+SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
 307      &  PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
 308   340   IF(MOD(MSTJ(3),5).EQ.3) PWS=PWS+1.
 309         DO 360 I=NSAV+NJET+1,N
 310         IF(MOD(MSTJ(3),5).EQ.1) PW=P(I,4)
 311         IF(MOD(MSTJ(3),5).EQ.2) PW=SQRT(P(I,5)**2+(PSI(1)*P(I,1)+
 312      &  PSI(2)*P(I,2)+PSI(3)*P(I,3))**2/PSI(4))
 313         IF(MOD(MSTJ(3),5).EQ.3) PW=1.
 314         DO 350 J=1,3
 315   350   P(I,J)=P(I,J)-PSI(J)*PW/PWS
 316   360   P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
 317
 318 C...Compensate for missing momentum withing each jet separately.
 319       ELSEIF(MOD(MSTJ(3),5).EQ.4) THEN
 320         DO 370 I=N+1,N+NJET
 321         K(I,1)=0
 322         DO 370 J=1,5
 323   370   P(I,J)=0.
 324         DO 390 I=NSAV+NJET+1,N
 325         IR1=K(I,3)
 326         IR2=N+IR1-NSAV
 327         K(IR2,1)=K(IR2,1)+1
 328         PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
 329      &  (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
 330         DO 380 J=1,3
 331   380   P(IR2,J)=P(IR2,J)+P(I,J)-PLS*P(IR1,J)
 332         P(IR2,4)=P(IR2,4)+P(I,4)
 333   390   P(IR2,5)=P(IR2,5)+PLS
 334         PSS=0.
 335         DO 400 I=N+1,N+NJET
 336   400   IF(K(I,1).NE.0) PSS=PSS+P(I,4)/(PECM*(0.8*P(I,5)+0.2))
 337         DO 420 I=NSAV+NJET+1,N
 338         IR1=K(I,3)
 339         IR2=N+IR1-NSAV
 340         PLS=(P(I,1)*P(IR1,1)+P(I,2)*P(IR1,2)+P(I,3)*P(IR1,3))/
 341      &  (P(IR1,1)**2+P(IR1,2)**2+P(IR1,3)**2)
 342         DO 410 J=1,3
 343   410   P(I,J)=P(I,J)-P(IR2,J)/K(IR2,1)+(1./(P(IR2,5)*PSS)-1.)*PLS*
 344      &  P(IR1,J)
 345   420   P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
 346       ENDIF
 347
 348 C...Scale momenta for energy conservation.
 349       IF(MOD(MSTJ(3),5).NE.0) THEN
 350         PMS=0.
 351         PES=0.
 352         PQS=0.
 353         DO 430 I=NSAV+NJET+1,N
 354         PMS=PMS+P(I,5)
 355         PES=PES+P(I,4)
 356   430   PQS=PQS+P(I,5)**2/P(I,4)
 357         IF(PMS.GE.PECM) GOTO 150
 358         NECO=0
 359   440   NECO=NECO+1
 360         PFAC=(PECM-PQS)/(PES-PQS)
 361         PES=0.
 362         PQS=0.
 363         DO 460 I=NSAV+NJET+1,N
 364         DO 450 J=1,3
 365   450   P(I,J)=PFAC*P(I,J)
 366         P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
 367         PES=PES+P(I,4)
 368   460   PQS=PQS+P(I,5)**2/P(I,4)
 369         IF(NECO.LT.10.AND.ABS(PECM-PES).GT.2E-6*PECM) GOTO 440
 370       ENDIF
 371
 372 C...Origin of produced particles and parton daughter pointers.
 373   470 DO 480 I=NSAV+NJET+1,N
 374       IF(MSTU(16).NE.2) K(I,3)=NSAV+1
 375   480 IF(MSTU(16).EQ.2) K(I,3)=K(K(I,3),3)
 376       DO 490 I=NSAV+1,NSAV+NJET
 377       I1=K(I,3)
 378       K(I1,1)=K(I1,1)+10
 379       IF(MSTU(16).NE.2) THEN
 380         K(I1,4)=NSAV+1
 381         K(I1,5)=NSAV+1
 382       ELSE
 383         K(I1,4)=K(I1,4)-NJET+1
 384         K(I1,5)=K(I1,5)-NJET+1
 385         IF(K(I1,5).LT.K(I1,4)) THEN
 386           K(I1,4)=0
 387           K(I1,5)=0
 388         ENDIF
 389       ENDIF
 390   490 CONTINUE
 391
 392 C...Document independent fragmentation system. Remove copy of jets.
 393       NSAV=NSAV+1
 394       K(NSAV,1)=11
 395       K(NSAV,2)=93
 396       K(NSAV,3)=IP
 397       K(NSAV,4)=NSAV+1
 398       K(NSAV,5)=N-NJET+1
 399       DO 500 J=1,4
 400       P(NSAV,J)=DPS(J)
 401   500 V(NSAV,J)=V(IP,J)
 402       P(NSAV,5)=SQRT(MAX(0D0,DPS(4)**2-DPS(1)**2-DPS(2)**2-DPS(3)**2))
 403       V(NSAV,5)=0.
 404       DO 510 I=NSAV+NJET,N
 405       DO 510 J=1,5
 406       K(I-NJET+1,J)=K(I,J)
 407       P(I-NJET+1,J)=P(I,J)
 408   510 V(I-NJET+1,J)=V(I,J)
 409       N=N-NJET+1
 410
 411 C...Boost back particle system. Set production vertices.
 412       IF(NJET.NE.1) CALL LUDBRB_HIJING(NSAV+1,N,0.,0.,DPS(1)/DPS(4),
 413      &DPS(2)/DPS(4),DPS(3)/DPS(4))
 414       DO 520 I=NSAV+1,N
 415       DO 520 J=1,4
 416   520 V(I,J)=V(IP,J)
 417
 418       RETURN
 419       END