Additional protection
[u/mrichter/AliRoot.git] / 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