]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HIJING/hijing1_36/quench.F
Fix long line
[u/mrichter/AliRoot.git] / HIJING / hijing1_36 / quench.F
CommitLineData
e74335a4 1* $Id$
2C
3C
4C
5C
6 SUBROUTINE QUENCH(JPJT,NTP)
7 DIMENSION RDP(300),LQP(300),RDT(300),LQT(300)
8#include "hijcrdn.inc"
9#include "hiparnt.inc"
10C
11#include "hijjet1.inc"
12#include "hijjet2.inc"
13#include "histrng.inc"
14C
15 SAVE
16 IF(NTP.EQ.2) GO TO 400
17 IF(NTP.EQ.3) GO TO 2000
18C*******************************************************
19C Jet interaction for proj jet in the direction PHIP
20C******************************************************
21C
22 IF(NFP(JPJT,7).NE.1) RETURN
23
24 JP=JPJT
25 DO 290 I=1,NPJ(JP)
26 PTJET0=SQRT(PJPX(JP,I)**2+PJPY(JP,I)**2)
27 IF(PTJET0.LE.HIPR1(11)) GO TO 290
28 PTOT=SQRT(PTJET0*PTJET0+PJPZ(JP,I)**2)
29 IF(PTOT.LT.HIPR1(8)) GO TO 290
30 PHIP=ULANGL_HIJING(PJPX(JP,I),PJPY(JP,I))
31C******* find the wounded proj which can interact with jet***
32 KP=0
33 DO 100 I2=1,IHNT2(1)
34 IF(NFP(I2,5).NE.3 .OR. I2.EQ.JP) GO TO 100
35 DX=YP(1,I2)-YP(1,JP)
36 DY=YP(2,I2)-YP(2,JP)
37 PHI=ULANGL_HIJING(DX,DY)
38 DPHI=ABS(PHI-PHIP)
39 IF(DPHI.GE.HIPR1(40)/2.0) GO TO 100
40 RD0=SQRT(DX*DX+DY*DY)
41 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 100
42 KP=KP+1
43 LQP(KP)=I2
44 RDP(KP)=COS(DPHI)*RD0
45 100 CONTINUE
46C******* rearrange according decending rd************
47 DO 110 I2=1,KP-1
48 DO 110 J2=I2+1,KP
49 IF(RDP(I2).LT.RDP(J2)) GO TO 110
50 RD=RDP(I2)
51 LQ=LQP(I2)
52 RDP(I2)=RDP(J2)
53 LQP(I2)=LQP(J2)
54 RDP(J2)=RD
55 LQP(J2)=LQ
56 110 CONTINUE
57C****** find wounded targ which can interact with jet********
58 KT=0
59 DO 120 I2=1,IHNT2(3)
60 IF(NFT(I2,5).NE.3) GO TO 120
61 DX=YT(1,I2)-YP(1,JP)
62 DY=YT(2,I2)-YP(2,JP)
63 PHI=ULANGL_HIJING(DX,DY)
64 DPHI=ABS(PHI-PHIP)
65 IF(DPHI.GT.HIPR1(40)/2.0) GO TO 120
66 RD0=SQRT(DX*DX+DY*DY)
67 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 120
68 KT=KT+1
69 LQT(KT)=I2
70 RDT(KT)=COS(DPHI)*RD0
71 120 CONTINUE
72C******* rearrange according decending rd************
73 DO 130 I2=1,KT-1
74 DO 130 J2=I2+1,KT
75 IF(RDT(I2).LT.RDT(J2)) GO TO 130
76 RD=RDT(I2)
77 LQ=LQT(I2)
78 RDT(I2)=RDT(J2)
79 LQT(I2)=LQT(J2)
80 RDT(J2)=RD
81 LQT(J2)=LQ
82 130 CONTINUE
83
84 MP=0
85 MT=0
86 R0=0.0
87 NQ=0
88 DP=0.0
89 PTOT=SQRT(PJPX(JP,I)**2+PJPY(JP,I)**2+PJPZ(JP,I)**2)
90 V1=PJPX(JP,I)/PTOT
91 V2=PJPY(JP,I)/PTOT
92 V3=PJPZ(JP,I)/PTOT
93
94 200 RN=RLU_HIJING(0)
95 210 IF(MT.GE.KT .AND. MP.GE.KP) GO TO 290
96 IF(MT.GE.KT) GO TO 220
97 IF(MP.GE.KP) GO TO 240
98 IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 240
99 220 MP=MP+1
100 DRR=RDP(MP)-R0
101 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 210
102 DP=DRR*HIPR1(14)
103 IF(KFPJ(JP,I).NE.21) DP=0.5*DP
104C ********string tension of quark jet is 0.5 of gluon's
105 IF(DP.LE.0.2) GO TO 210
106 IF(PTOT.LE.0.4) GO TO 290
107 IF(PTOT.LE.DP) DP=PTOT-0.2
108 DE=DP
109
110 IF(KFPJ(JP,I).NE.21) THEN
111 PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
112 & +PP(LQP(MP),3)**2
113 DE=SQRT(PJPM(JP,I)**2+PTOT**2)
114 & -SQRT(PJPM(JP,I)**2+(PTOT-DP)**2)
115 ERSHU=(PP(LQP(MP),4)+DE-DP)**2
116 AMSHU=ERSHU-PRSHU
117 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 210
118 PP(LQP(MP),4)=SQRT(ERSHU)
119 PP(LQP(MP),5)=SQRT(AMSHU)
120 ENDIF
121C ********reshuffle the energy when jet has mass
122 R0=RDP(MP)
123 DP1=DP*V1
124 DP2=DP*V2
125 DP3=DP*V3
126C ********momentum and energy transfer from jet
127
128 NPJ(LQP(MP))=NPJ(LQP(MP))+1
129 KFPJ(LQP(MP),NPJ(LQP(MP)))=21
130 PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
131 PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
132 PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
133 PJPE(LQP(MP),NPJ(LQP(MP)))=DP
134 PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
135 GO TO 260
136
137 240 MT=MT+1
138 DRR=RDT(MT)-R0
139 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 210
140 DP=DRR*HIPR1(14)
141 IF(DP.LE.0.2) GO TO 210
142 IF(PTOT.LE.0.4) GO TO 290
143 IF(PTOT.LE.DP) DP=PTOT-0.2
144 DE=DP
145
146 IF(KFPJ(JP,I).NE.21) THEN
147 PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
148 & +PT(LQT(MT),3)**2
149 DE=SQRT(PJPM(JP,I)**2+PTOT**2)
150 & -SQRT(PJPM(JP,I)**2+(PTOT-DP)**2)
151 ERSHU=(PT(LQT(MT),4)+DE-DP)**2
152 AMSHU=ERSHU-PRSHU
153 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 210
154 PT(LQT(MT),4)=SQRT(ERSHU)
155 PT(LQT(MT),5)=SQRT(AMSHU)
156 ENDIF
157C ********reshuffle the energy when jet has mass
158
159 R0=RDT(MT)
160 DP1=DP*V1
161 DP2=DP*V2
162 DP3=DP*V3
163C ********momentum and energy transfer from jet
164 NTJ(LQT(MT))=NTJ(LQT(MT))+1
165 KFTJ(LQT(MT),NTJ(LQT(MT)))=21
166 PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
167 PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
168 PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
169 PJTE(LQT(MT),NTJ(LQT(MT)))=DP
170 PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
171
172 260 PJPX(JP,I)=(PTOT-DP)*V1
173 PJPY(JP,I)=(PTOT-DP)*V2
174 PJPZ(JP,I)=(PTOT-DP)*V3
175 PJPE(JP,I)=PJPE(JP,I)-DE
176
177 PTOT=PTOT-DP
178 NQ=NQ+1
179 GO TO 200
180 290 CONTINUE
181
182 RETURN
183
184C*******************************************************
185C Jet interaction for target jet in the direction PHIT
186C******************************************************
187C
188C******* find the wounded proj which can interact with jet***
189
190 400 IF(NFT(JPJT,7).NE.1) RETURN
191 JT=JPJT
192 DO 690 I=1,NTJ(JT)
193 PTJET0=SQRT(PJTX(JT,I)**2+PJTY(JT,I)**2)
194 IF(PTJET0.LE.HIPR1(11)) GO TO 690
195 PTOT=SQRT(PTJET0*PTJET0+PJTZ(JT,I)**2)
196 IF(PTOT.LT.HIPR1(8)) GO TO 690
197 PHIT=ULANGL_HIJING(PJTX(JT,I),PJTY(JT,I))
198 KP=0
199 DO 500 I2=1,IHNT2(1)
200 IF(NFP(I2,5).NE.3) GO TO 500
201 DX=YP(1,I2)-YT(1,JT)
202 DY=YP(2,I2)-YT(2,JT)
203 PHI=ULANGL_HIJING(DX,DY)
204 DPHI=ABS(PHI-PHIT)
205 IF(DPHI.GT.HIPR1(40)/2.0) GO TO 500
206 RD0=SQRT(DX*DX+DY*DY)
207 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 500
208 KP=KP+1
209 LQP(KP)=I2
210 RDP(KP)=COS(DPHI)*RD0
211 500 CONTINUE
212C******* rearrange according to decending rd************
213 DO 510 I2=1,KP-1
214 DO 510 J2=I2+1,KP
215 IF(RDP(I2).LT.RDP(J2)) GO TO 510
216 RD=RDP(I2)
217 LQ=LQP(I2)
218 RDP(I2)=RDP(J2)
219 LQP(I2)=LQP(J2)
220 RDP(J2)=RD
221 LQP(J2)=LQ
222 510 CONTINUE
223C****** find wounded targ which can interact with jet********
224 KT=0
225 DO 520 I2=1,IHNT2(3)
226 IF(NFT(I2,5).NE.3 .OR. I2.EQ.JT) GO TO 520
227 DX=YT(1,I2)-YT(1,JT)
228 DY=YT(2,I2)-YT(2,JT)
229 PHI=ULANGL_HIJING(DX,DY)
230 DPHI=ABS(PHI-PHIT)
231 IF(DPHI.GT.HIPR1(40)/2.0) GO TO 520
232 RD0=SQRT(DX*DX+DY*DY)
233 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 520
234 KT=KT+1
235 LQT(KT)=I2
236 RDT(KT)=COS(DPHI)*RD0
237 520 CONTINUE
238C******* rearrange according to decending rd************
239 DO 530 I2=1,KT-1
240 DO 530 J2=I2+1,KT
241 IF(RDT(I2).LT.RDT(J2)) GO TO 530
242 RD=RDT(I2)
243 LQ=LQT(I2)
244 RDT(I2)=RDT(J2)
245 LQT(I2)=LQT(J2)
246 RDT(J2)=RD
247 LQT(J2)=LQ
248 530 CONTINUE
249
250 MP=0
251 MT=0
252 NQ=0
253 DP=0.0
254 R0=0.0
255 PTOT=SQRT(PJTX(JT,I)**2+PJTY(JT,I)**2+PJTZ(JT,I)**2)
256 V1=PJTX(JT,I)/PTOT
257 V2=PJTY(JT,I)/PTOT
258 V3=PJTZ(JT,I)/PTOT
259
260 600 RN=RLU_HIJING(0)
261 610 IF(MT.GE.KT .AND. MP.GE.KP) GO TO 690
262 IF(MT.GE.KT) GO TO 620
263 IF(MP.GE.KP) GO TO 640
264 IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 640
265620 MP=MP+1
266 DRR=RDP(MP)-R0
267 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 610
268 DP=DRR*HIPR1(14)
269 IF(KFTJ(JT,I).NE.21) DP=0.5*DP
270C ********string tension of quark jet is 0.5 of gluon's
271 IF(DP.LE.0.2) GO TO 610
272 IF(PTOT.LE.0.4) GO TO 690
273 IF(PTOT.LE.DP) DP=PTOT-0.2
274 DE=DP
275C
276 IF(KFTJ(JT,I).NE.21) THEN
277 PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
278 & +PP(LQP(MP),3)**2
279 DE=SQRT(PJTM(JT,I)**2+PTOT**2)
280 & -SQRT(PJTM(JT,I)**2+(PTOT-DP)**2)
281 ERSHU=(PP(LQP(MP),4)+DE-DP)**2
282 AMSHU=ERSHU-PRSHU
283 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 610
284 PP(LQP(MP),4)=SQRT(ERSHU)
285 PP(LQP(MP),5)=SQRT(AMSHU)
286 ENDIF
287C ********reshuffle the energy when jet has mass
288C
289 R0=RDP(MP)
290 DP1=DP*V1
291 DP2=DP*V2
292 DP3=DP*V3
293C ********momentum and energy transfer from jet
294 NPJ(LQP(MP))=NPJ(LQP(MP))+1
295 KFPJ(LQP(MP),NPJ(LQP(MP)))=21
296 PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
297 PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
298 PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
299 PJPE(LQP(MP),NPJ(LQP(MP)))=DP
300 PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
301
302 GO TO 660
303
304640 MT=MT+1
305 DRR=RDT(MT)-R0
306 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 610
307 DP=DRR*HIPR1(14)
308 IF(DP.LE.0.2) GO TO 610
309 IF(PTOT.LE.0.4) GO TO 690
310 IF(PTOT.LE.DP) DP=PTOT-0.2
311 DE=DP
312
313 IF(KFTJ(JT,I).NE.21) THEN
314 PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
315 & +PT(LQT(MT),3)**2
316 DE=SQRT(PJTM(JT,I)**2+PTOT**2)
317 & -SQRT(PJTM(JT,I)**2+(PTOT-DP)**2)
318 ERSHU=(PT(LQT(MT),4)+DE-DP)**2
319 AMSHU=ERSHU-PRSHU
320 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 610
321 PT(LQT(MT),4)=SQRT(ERSHU)
322 PT(LQT(MT),5)=SQRT(AMSHU)
323 ENDIF
324C ********reshuffle the energy when jet has mass
325
326 R0=RDT(MT)
327 DP1=DP*V1
328 DP2=DP*V2
329 DP3=DP*V3
330C ********momentum and energy transfer from jet
331 NTJ(LQT(MT))=NTJ(LQT(MT))+1
332 KFTJ(LQT(MT),NTJ(LQT(MT)))=21
333 PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
334 PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
335 PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
336 PJTE(LQT(MT),NTJ(LQT(MT)))=DP
337 PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
338
339660 PJTX(JT,I)=(PTOT-DP)*V1
340 PJTY(JT,I)=(PTOT-DP)*V2
341 PJTZ(JT,I)=(PTOT-DP)*V3
342 PJTE(JT,I)=PJTE(JT,I)-DE
343
344 PTOT=PTOT-DP
345 NQ=NQ+1
346 GO TO 600
347690 CONTINUE
348 RETURN
349C********************************************************
350C Q-QBAR jet interaction
351C********************************************************
3522000 ISG=JPJT
353 IF(IASG(ISG,3).NE.1) RETURN
354C
355 JP=IASG(ISG,1)
356 JT=IASG(ISG,2)
357 XJ=(YP(1,JP)+YT(1,JT))/2.0
358 YJ=(YP(2,JP)+YT(2,JT))/2.0
359 DO 2690 I=1,NJSG(ISG)
360 PTJET0=SQRT(PXSG(ISG,I)**2+PYSG(ISG,I)**2)
361 IF(PTJET0.LE.HIPR1(11).OR.PESG(ISG,I).LT.HIPR1(1))
362 & GO TO 2690
363 PTOT=SQRT(PTJET0*PTJET0+PZSG(ISG,I)**2)
364 IF(PTOT.LT.MAX(HIPR1(1),HIPR1(8))) GO TO 2690
365 PHIQ=ULANGL_HIJING(PXSG(ISG,I),PYSG(ISG,I))
366 KP=0
367 DO 2500 I2=1,IHNT2(1)
368 IF(NFP(I2,5).NE.3.OR.I2.EQ.JP) GO TO 2500
369 DX=YP(1,I2)-XJ
370 DY=YP(2,I2)-YJ
371 PHI=ULANGL_HIJING(DX,DY)
372 DPHI=ABS(PHI-PHIQ)
373 IF(DPHI.GT.HIPR1(40)/2.0) GO TO 2500
374 RD0=SQRT(DX*DX+DY*DY)
375 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 2500
376 KP=KP+1
377 LQP(KP)=I2
378 RDP(KP)=COS(DPHI)*RD0
379 2500 CONTINUE
380C******* rearrange according to decending rd************
381 DO 2510 I2=1,KP-1
382 DO 2510 J2=I2+1,KP
383 IF(RDP(I2).LT.RDP(J2)) GO TO 2510
384 RD=RDP(I2)
385 LQ=LQP(I2)
386 RDP(I2)=RDP(J2)
387 LQP(I2)=LQP(J2)
388 RDP(J2)=RD
389 LQP(J2)=LQ
390 2510 CONTINUE
391C****** find wounded targ which can interact with jet********
392 KT=0
393 DO 2520 I2=1,IHNT2(3)
394 IF(NFT(I2,5).NE.3 .OR. I2.EQ.JT) GO TO 2520
395 DX=YT(1,I2)-XJ
396 DY=YT(2,I2)-YJ
397 PHI=ULANGL_HIJING(DX,DY)
398 DPHI=ABS(PHI-PHIQ)
399 IF(DPHI.GT.HIPR1(40)/2.0) GO TO 2520
400 RD0=SQRT(DX*DX+DY*DY)
401 IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 2520
402 KT=KT+1
403 LQT(KT)=I2
404 RDT(KT)=COS(DPHI)*RD0
405 2520 CONTINUE
406C******* rearrange according to decending rd************
407 DO 2530 I2=1,KT-1
408 DO 2530 J2=I2+1,KT
409 IF(RDT(I2).LT.RDT(J2)) GO TO 2530
410 RD=RDT(I2)
411 LQ=LQT(I2)
412 RDT(I2)=RDT(J2)
413 LQT(I2)=LQT(J2)
414 RDT(J2)=RD
415 LQT(J2)=LQ
416 2530 CONTINUE
417
418 MP=0
419 MT=0
420 NQ=0
421 DP=0.0
422 R0=0.0
423 PTOT=SQRT(PXSG(ISG,I)**2+PYSG(ISG,I)**2
424 & +PZSG(ISG,I)**2)
425 V1=PXSG(ISG,I)/PTOT
426 V2=PYSG(ISG,I)/PTOT
427 V3=PZSG(ISG,I)/PTOT
428
429 2600 RN=RLU_HIJING(0)
430 2610 IF(MT.GE.KT .AND. MP.GE.KP) GO TO 2690
431 IF(MT.GE.KT) GO TO 2620
432 IF(MP.GE.KP) GO TO 2640
433 IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 2640
434 2620 MP=MP+1
435 DRR=RDP(MP)-R0
436 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 2610
437 DP=DRR*HIPR1(14)/2.0
438 IF(DP.LE.0.2) GO TO 2610
439 IF(PTOT.LE.0.4) GO TO 2690
440 IF(PTOT.LE.DP) DP=PTOT-0.2
441 DE=DP
442C
443 IF(K2SG(ISG,I).NE.21) THEN
444 IF(PTOT.LT.DP+HIPR1(1)) GO TO 2690
445 PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
446 & +PP(LQP(MP),3)**2
447 DE=SQRT(PMSG(ISG,I)**2+PTOT**2)
448 & -SQRT(PMSG(ISG,I)**2+(PTOT-DP)**2)
449 ERSHU=(PP(LQP(MP),4)+DE-DP)**2
450 AMSHU=ERSHU-PRSHU
451 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 2610
452 PP(LQP(MP),4)=SQRT(ERSHU)
453 PP(LQP(MP),5)=SQRT(AMSHU)
454 ENDIF
455C ********reshuffle the energy when jet has mass
456C
457 R0=RDP(MP)
458 DP1=DP*V1
459 DP2=DP*V2
460 DP3=DP*V3
461C ********momentum and energy transfer from jet
462 NPJ(LQP(MP))=NPJ(LQP(MP))+1
463 KFPJ(LQP(MP),NPJ(LQP(MP)))=21
464 PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
465 PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
466 PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
467 PJPE(LQP(MP),NPJ(LQP(MP)))=DP
468 PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
469
470 GO TO 2660
471
472 2640 MT=MT+1
473 DRR=RDT(MT)-R0
474 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 2610
475 DP=DRR*HIPR1(14)
476 IF(DP.LE.0.2) GO TO 2610
477 IF(PTOT.LE.0.4) GO TO 2690
478 IF(PTOT.LE.DP) DP=PTOT-0.2
479 DE=DP
480
481 IF(K2SG(ISG,I).NE.21) THEN
482 IF(PTOT.LT.DP+HIPR1(1)) GO TO 2690
483 PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
484 & +PT(LQT(MT),3)**2
485 DE=SQRT(PMSG(ISG,I)**2+PTOT**2)
486 & -SQRT(PMSG(ISG,I)**2+(PTOT-DP)**2)
487 ERSHU=(PT(LQT(MT),4)+DE-DP)**2
488 AMSHU=ERSHU-PRSHU
489 IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 2610
490 PT(LQT(MT),4)=SQRT(ERSHU)
491 PT(LQT(MT),5)=SQRT(AMSHU)
492 ENDIF
493C ********reshuffle the energy when jet has mass
494
495 R0=RDT(MT)
496 DP1=DP*V1
497 DP2=DP*V2
498 DP3=DP*V3
499C ********momentum and energy transfer from jet
500 NTJ(LQT(MT))=NTJ(LQT(MT))+1
501 KFTJ(LQT(MT),NTJ(LQT(MT)))=21
502 PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
503 PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
504 PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
505 PJTE(LQT(MT),NTJ(LQT(MT)))=DP
506 PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
507
508 2660 PXSG(ISG,I)=(PTOT-DP)*V1
509 PYSG(ISG,I)=(PTOT-DP)*V2
510 PZSG(ISG,I)=(PTOT-DP)*V3
511 PESG(ISG,I)=PESG(ISG,I)-DE
512
513 PTOT=PTOT-DP
514 NQ=NQ+1
515 GO TO 2600
516 2690 CONTINUE
517 RETURN
518 END