Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyrvch.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyrvch.f
1 
2 C*********************************************************************
3 
4 C...PYRVCH
5 C...Calculates R-violating chargino decay widths.
6 C...P. Z. Skands
7 
8  SUBROUTINE pyrvch(KFIN,XLAM,IDLAM,LKNT)
9 
10 C...Double precision and integer declarations.
11  IMPLICIT DOUBLE PRECISION(a-h, o-z)
12  IMPLICIT INTEGER(i-n)
13 C...Parameter statement to help give large particle numbers.
14  parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
15  &kexcit=4000000,kdimen=5000000)
16 C...Commonblocks.
17  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
18  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
19  common/pymssm/imss(0:99),rmss(0:99)
20  common/pyssmt/zmix(4,4),umix(2,2),vmix(2,2),smz(4),smw(2),
21  &sfmix(16,4),zmixi(4,4),umixi(2,2),vmixi(2,2)
22  common/pymsrv/rvlam(3,3,3), rvlamp(3,3,3), rvlamb(3,3,3)
23 C...Local variables.
24  DOUBLE PRECISION xlam(0:400)
25  INTEGER idlam(400,3), pycomp
26 C...Information from main routine to PYRVGW
27  common/pyrvnv/ab(2,16,2),rms(0:3),res(6,2),intres(6,3),idr,idr2
28  & ,dcmass,kfr(3)
29 C...Auxiliary variables needed for BV (RV Gauge STOre)
30  common/rvgsto/xresi,xresj,xresk,xresij,xresik,xresjk,rvlijk,rvlkij
31  & ,rvljki,rvljik
32 C...Running quark masses
33  DOUBLE PRECISION rmq(6)
34 C...Decay product masses on/off
35  LOGICAL dcmass
36  SAVE /pydat1/,/pydat2/,/pymssm/,/pyssmt/,/pymsrv/,/pyrvnv/,
37  & /rvgsto/
38 
39 
40 C...IF R-VIOLATION ON.
41  IF ((imss(51).GE.1).OR.(imss(52).GE.1).OR.(imss(53).GE.1)) THEN
42  kfsm=kfin-ksusy1
43  IF(kfsm.EQ.24.OR.kfsm.EQ.37) THEN
44 C...WHICH CHARGINO ?
45  nchi = 1
46  IF (kfsm.EQ.37) nchi = 2
47 
48 C...Useful parameters for calculating the A and B constants.
49 C...SIGN OF MASS (Opposite convention as HERWIG)
50  ism = 1
51  IF (smw(nchi).LT.0d0) ism = -1
52  wmass = pmas(pycomp(24),1)
53  cosb = 1/(sqrt(1+rmss(5)**2))
54  sinb = rmss(5)/sqrt(1+rmss(5)**2)
55  gw2 = 4*paru(103)*paru(1)/paru(102)
56  c1u = umix(nchi,2)/(sqrt(2d0)*cosb*wmass)
57  c1v = vmix(nchi,2)/(sqrt(2d0)*sinb*wmass)
58  c2 = umix(nchi,1)
59  c3 = vmix(nchi,1)
60 C...Running masses at Q^2=MCHI^2.
61  sqmchi = pmas(pycomp(kfsm),1)**2
62  DO 100 i=1,6
63  rmq(i)=pymrun(i,sqmchi)
64  100 CONTINUE
65 
66 C... AB(x,y,z) coefficients:
67 C x=1-2 : A or B coefficient (1:A ; 2:B)
68 C y=1-16 : Sparticle's SM code (1-6:d,u,s,c,b,t ;
69 C 11-16:e,nu_e,mu,...)
70 C z=1-2 : Mass eigenstate number
71  DO 110 i = 11,15,2
72 C...Intermediate sleptons
73  ab(1,i,1) = 0d0
74  ab(1,i,2) = 0d0
75  ab(2,i,1) = -pmas(pycomp(i),1)*c1u*sfmix(i,2) +
76  & sfmix(i,1)*c2
77  ab(2,i,2) = -pmas(pycomp(i),1)*c1u*sfmix(i,4) +
78  & sfmix(i,3)*c2
79 C...Intermediate sneutrinos
80  ab(1,i+1,1) = -pmas(pycomp(i),1)*c1u
81  ab(1,i+1,2) = 0d0
82  ab(2,i+1,1) = ism*c3
83  ab(2,i+1,2) = 0d0
84 C...Intermediate sdown
85  j=i-10
86  ab(1,j,1) = -rmq(j+1)*c1v*sfmix(j,1)
87  ab(1,j,2) = -rmq(j+1)*c1v*sfmix(j,3)
88  ab(2,j,1) = -ism*(rmq(j)*c1u*sfmix(j,2) - sfmix(j,1)*c2)
89  ab(2,j,2) = -ism*(rmq(j)*c1u*sfmix(j,4) - sfmix(j,3)*c2)
90 C...Intermediate sup
91  j=j+1
92  ab(1,j,1) = -rmq(j-1)*c1u*sfmix(j,1)
93  ab(1,j,2) = -rmq(j-1)*c1u*sfmix(j,3)
94  ab(2,j,1) = -ism*(rmq(j)*c1v*sfmix(j,2) - sfmix(j,1)*c3)
95  ab(2,j,2) = -ism*(rmq(j)*c1v*sfmix(j,4) - sfmix(j,3)*c3)
96  110 CONTINUE
97 
98 C...LLE TYPE R-VIOLATION
99  IF (imss(51).GE.1) THEN
100 C...LOOP OVER DECAY MODES
101  DO 140 isc=0,26
102 
103 C...CHI+ -> NUBAR_I + LEPTON+_J + NU_K.
104  IF(mod(isc/9,3).NE.mod(isc/3,3)) THEN
105  lknt = lknt+1
106  idlam(lknt,1) = -12 -2*mod(isc/9,3)
107  idlam(lknt,2) = -11 -2*mod(isc/3,3)
108  idlam(lknt,3) = 12 +2*mod(isc,3)
109  xlam(lknt) = 0d0
110 C...Set coupling, and decay product masses on/off
111  rvlamc = gw2 * 5d-1 *
112  & rvlam(mod(isc/9,3)+1,mod(isc/3,3)+1,mod(isc,3)+1)
113  & **2
114  dcmass=.false.
115  IF (idlam(lknt,2).EQ.-15) dcmass = .true.
116 C...Resonance KF codes (1=I,2=J,3=K).
117  kfr(1) = 0
118  kfr(2) = 0
119  kfr(3) = -idlam(lknt,3)+1
120 C...Calculate width.
121  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),
122  & idlam(lknt,3),xlam(lknt))
123  xlam(lknt)=xlam(lknt)*rvlamc/((2*paru(1)*rms(0))**3*32)
124 C...KINEMATICS CHECK
125  IF (xlam(lknt).EQ.0d0) THEN
126  lknt=lknt-1
127  ENDIF
128 
129 C * CHI+ -> NU_I + NU_J + LEPTON+_K. (NOTE: SYMM. IN I AND J)
130  120 IF (mod(isc/9,3).LT.mod(isc/3,3)) THEN
131  lknt = lknt+1
132  idlam(lknt,1) = 12 +2*mod(isc/9,3)
133  idlam(lknt,2) = 12 +2*mod(isc/3,3)
134  idlam(lknt,3) =-11 -2*mod(isc,3)
135  xlam(lknt) = 0d0
136 C...Set coupling, and decay product masses on/off
137  rvlamc = gw2 * 5d-1 *
138  & rvlam(mod(isc/9,3)+1,mod(isc/3,3)+1,mod(isc,3)+1)**2
139 C...I,J SYMMETRY => FACTOR 2
140  rvlamc=2*rvlamc
141  dcmass=.false.
142  IF (idlam(lknt,3).EQ.-15) dcmass = .true.
143 C...Resonance KF codes (1=I,2=J,3=K)
144  kfr(1)=idlam(lknt,1)-1
145  kfr(2)=idlam(lknt,2)-1
146  kfr(3)=0
147 C...Calculate width.
148  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),
149  & idlam(lknt,3),xlam(lknt))
150  xlam(lknt)=xlam(lknt)*rvlamc/((2*paru(1)*rms(0))**3*32)
151 C...KINEMATICS CHECK
152  IF (xlam(lknt).EQ.0d0) THEN
153  lknt=lknt-1
154  ENDIF
155  130 ENDIF
156 
157 C * CHI+ -> LEPTON+_I + LEPTON+_J + LEPTON-_K
158  lknt = lknt+1
159  idlam(lknt,1) =-11 -2*mod(isc/9,3)
160  idlam(lknt,2) =-11 -2*mod(isc/3,3)
161  idlam(lknt,3) = 11 +2*mod(isc,3)
162  xlam(lknt) = 0d0
163 C...Set coupling, and decay product masses on/off
164  rvlamc = gw2 * 5d-1 *
165  & rvlam(mod(isc/9,3)+1,mod(isc/3,3)+1,mod(isc,3)+1)**2
166 C...I,J SYMMETRY => FACTOR 2
167  rvlamc=2*rvlamc
168  dcmass=.false.
169  IF (idlam(lknt,1).EQ.-15.OR.idlam(lknt,2).EQ.-15
170  & .OR.idlam(lknt,3).EQ.15) dcmass = .true.
171 C...Resonance KF codes (1=I,2=J,3=K)
172  kfr(1) =-idlam(lknt,1)+1
173  kfr(2) =-idlam(lknt,2)+1
174  kfr(3) = 0
175 C...Calculate width.
176  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),
177  & idlam(lknt,3),xlam(lknt))
178  xlam(lknt)=xlam(lknt)*rvlamc/((2*paru(1)*rms(0))**3*32)
179 C...KINEMATICS CHECK
180  IF (xlam(lknt).EQ.0d0) THEN
181  lknt=lknt-1
182  ENDIF
183  ENDIF
184  140 CONTINUE
185  ENDIF
186 
187 C...LQD TYPE R-VIOLATION
188  IF (imss(52).GE.1) THEN
189 C...LOOP OVER DECAY MODES
190  DO 180 isc=0,26
191 
192 C...CHI+ -> NUBAR_I + DBAR_J + U_K
193  lknt = lknt+1
194  idlam(lknt,1) =-12 -2*mod(isc/9,3)
195  idlam(lknt,2) = -1 -2*mod(isc/3,3)
196  idlam(lknt,3) = 2 +2*mod(isc,3)
197  xlam(lknt) = 0d0
198 C...Set coupling, and decay product masses on/off
199  rvlamc = 3. * gw2 * 5d-1 *
200  & rvlamp(mod(isc/9,3)+1,mod(isc/3,3)+1,mod(isc,3)+1)**2
201  dcmass=.false.
202  IF (idlam(lknt,2).EQ.-5.OR.idlam(lknt,3).EQ.6)
203  & dcmass = .true.
204 C...Resonance KF codes (1=I,2=J,3=K)
205  kfr(1)=0
206  kfr(2)=0
207  kfr(3)=-idlam(lknt,3)+1
208 C...Calculate width.
209  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),idlam(lknt,3)
210  & ,xlam(lknt))
211  xlam(lknt)=xlam(lknt)*rvlamc/((2*paru(1)*rms(0))**3*32)
212 C...KINEMATICS CHECK
213  IF (xlam(lknt).EQ.0d0) THEN
214  lknt=lknt-1
215  ENDIF
216 
217 C * CHI+ -> LEPTON+_I + UBAR_J + U_K.
218  150 lknt = lknt+1
219  idlam(lknt,1) =-11 -2*mod(isc/9,3)
220  idlam(lknt,2) = -2 -2*mod(isc/3,3)
221  idlam(lknt,3) = 2 +2*mod(isc,3)
222  xlam(lknt) = 0d0
223 C...Set coupling, and decay product masses on/off
224  rvlamc = 3. * gw2 * 5d-1 *
225  & rvlamp(mod(isc/9,3)+1,mod(isc/3,3)+1,mod(isc,3)+1)**2
226  dcmass=.false.
227  IF (idlam(lknt,1).EQ.-11.OR.idlam(lknt,2).EQ.-6
228  & .OR.idlam(lknt,3).EQ.6) dcmass = .true.
229 C...Resonance KF codes (1=I,2=J,3=K)
230  kfr(1)=0
231  kfr(2)=0
232  kfr(3)=-idlam(lknt,3)+1
233 C...Calculate width.
234  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),idlam(lknt,3)
235  & ,xlam(lknt))
236  xlam(lknt)=xlam(lknt)*rvlamc/((2*paru(1)*rms(0))**3*32)
237 C...KINEMATICS CHECK
238  IF (xlam(lknt).EQ.0d0) THEN
239  lknt=lknt-1
240  ENDIF
241 
242 C * CHI+ -> LEPTON+_I + DBAR_J + D_K.
243  160 lknt = lknt+1
244  idlam(lknt,1) =-11 -2*mod(isc/9,3)
245  idlam(lknt,2) = -1 -2*mod(isc/3,3)
246  idlam(lknt,3) = 1 +2*mod(isc,3)
247  xlam(lknt) = 0d0
248 C...Set coupling, and decay product masses on/off
249  rvlamc = 3. * gw2 * 5d-1 *
250  & rvlamp(mod(isc/9,3)+1,mod(isc/3,3)+1,mod(isc,3)+1)**2
251  dcmass = .false.
252  IF (idlam(lknt,1).EQ.-15.OR.idlam(lknt,2).EQ.-5
253  & .OR.idlam(lknt,3).EQ.5) dcmass = .true.
254 C...Resonance KF codes (1=I,2=J,3=K)
255  kfr(1)=-idlam(lknt,1)+1
256  kfr(2)=-idlam(lknt,2)+1
257  kfr(3)=0
258 C...Calculate width.
259  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),idlam(lknt,3)
260  & ,xlam(lknt))
261  xlam(lknt)=xlam(lknt)*rvlamc/((2*paru(1)*rms(0))**3*32)
262 C...KINEMATICS CHECK
263  IF (xlam(lknt).EQ.0d0) THEN
264  lknt=lknt-1
265  ENDIF
266 
267 C * CHI+ -> NU_I + U_J + DBAR_K.
268  170 lknt = lknt+1
269  idlam(lknt,1) = 12 +2*mod(isc/9,3)
270  idlam(lknt,2) = 2 +2*mod(isc/3,3)
271  idlam(lknt,3) = -1 -2*mod(isc,3)
272  xlam(lknt) = 0d0
273 C...Set coupling, and decay product masses on/off
274  dcmass = .false.
275  rvlamc = 3. * gw2 * 5d-1 *
276  & rvlamp(mod(isc/9,3)+1,mod(isc/3,3)+1,mod(isc,3)+1)**2
277  IF (idlam(lknt,2).EQ.6.OR.idlam(lknt,3).EQ.-5)
278  & dcmass = .true.
279 C...Resonance KF codes (1=I,2=J,3=K)
280  kfr(1)=idlam(lknt,1)-1
281  kfr(2)=idlam(lknt,2)-1
282  kfr(3)=0
283 C...Calculate width.
284  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),idlam(lknt,3)
285  & ,xlam(lknt))
286  xlam(lknt)=xlam(lknt)*rvlamc/((2*paru(1)*rms(0))**3*32)
287 C...KINEMATICS CHECK
288  IF (xlam(lknt).EQ.0d0) THEN
289  lknt=lknt-1
290  ENDIF
291 
292  180 CONTINUE
293  ENDIF
294 
295 C...UDD TYPE R-VIOLATION
296 C...These decays need special treatment since more than one BV coupling
297 C...contributes (with interference). Consider e.g. (symbolically)
298 C |M|^2 = |l''_{ijk}|^2*(PYRVI1(RES_I) + PYRVI2(RES_I))
299 C +|l''_{jik}|^2*(PYRVI1(RES_J) + PYRVI2(RES_J))
300 C +l''_{ijk}*l''_{jik}*PYRVI3(PYRVI4(RES_I,RES_J))
301 C...The problem is that a single call to PYRVGW would evaluate all
302 C...these terms and sum them, but without the different couplings. The
303 C...way out is to call PYRVGW three times, once for the first line, once
304 C...for the second line, and then once for all the lines (it is
305 C...impossible to get just the last line out) without multiplying by
306 C...couplings. The last line is then obtained as the result of the third
307 C...call minus the results of the two first calls. Each term is then
308 C...multiplied by its respective coupling before the whole thing is
309 C...summed up in XLAM.
310 C...Note that with three interfering resonances, this procedure becomes
311 C...more complicated, as can be seen in the CHI+ -> 3*DBAR mode.
312 
313  IF (imss(53).GE.1) THEN
314 C...LOOP OVER DECAY MODES
315  DO 190 isc=1,25
316 
317 C...CHI+ -> U_I + U_J + D_K
318 C...Decay mode I<->J symmetric.
319  IF (mod(isc/9,3).LE.mod(isc/3,3).AND.isc.NE.13) THEN
320  lknt = lknt+1
321  idlam(lknt,1) = 2 +2*mod(isc/9,3)
322  idlam(lknt,2) = 2 +2*mod(isc/3,3)
323  idlam(lknt,3) = 1 +2*mod(isc,3)
324  xlam(lknt) = 0d0
325 C...Set coupling, and decay product masses on/off
326  rvlamc= 6. * gw2 * 5d-1
327  rvljik= rvlamb(mod(isc/3,3)+1,mod(isc/9,3)+1,mod(isc,3)
328  & +1)
329  rvlijk= rvlamb(mod(isc/9,3)+1,mod(isc/3,3)+1,mod(isc,3)
330  & +1)
331  IF (mod(isc/9,3).EQ.mod(isc/3,3)) rvlamc = 5d-1
332  & * rvlamc
333  dcmass=.false.
334  IF (idlam(lknt,1).EQ.6.OR.idlam(lknt,2).EQ.6
335  & .OR.idlam(lknt,3).EQ.5) dcmass =.true.
336 C...Resonance KF codes (1=I,2=J,3=K)
337  kfr(1) = -idlam(lknt,1)+1
338  kfr(2) = 0
339  kfr(3) = 0
340 C...Calculate width.
341  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),
342  & idlam(lknt,3),xresi)
343 C...Resonance KF codes (1=I,2=J,3=K)
344  kfr(1) = 0
345  kfr(2) = -idlam(lknt,2)+1
346  kfr(3) = 0
347 C...Calculate width.
348  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),
349  & idlam(lknt,3),xresj)
350 C...Resonance KF codes (1=I,2=J,3=K)
351  kfr(1) = -idlam(lknt,1)+1
352  kfr(2) = -idlam(lknt,2)+1
353  kfr(3) = 0
354 C...Calculate width.
355  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),
356  & idlam(lknt,3),xresij)
357  IF (abs((xresi+xresj)/xresij-1.).GT.1d-4) THEN
358  xresij = xresij-xresi-xresj
359  ELSE
360  xresij = 0d0
361  ENDIF
362 C...CALCULATE TOTAL WIDTH
363  xlam(lknt) = rvljik**2 * xresi + rvlijk**2 * xresj
364  & + rvljik*rvlijk * xresij
365  xlam(lknt)=xlam(lknt)*rvlamc/((2*paru(1)*rms(0))**3*32)
366 C...KINEMATICS CHECK
367  IF (xlam(lknt).EQ.0d0) THEN
368  lknt=lknt-1
369  ENDIF
370  ENDIF
371 C...CHI+ -> DBAR_I + DBAR_J + DBAR_K
372 C...Symmetry I<->J<->K.
373  IF ((mod(isc/9,3).LE.mod(isc/3,3)).AND.(mod(isc/3,3).le
374  & .mod(isc,3)).AND.isc.NE.13) THEN
375  lknt = lknt+1
376  idlam(lknt,1) = -1 -2*mod(isc/9,3)
377  idlam(lknt,2) = -1 -2*mod(isc/3,3)
378  idlam(lknt,3) = -1 -2*mod(isc,3)
379  xlam(lknt) = 0d0
380 C...Set coupling, and decay product masses on/off
381  rvlamc = 6. * gw2 * 5d-1
382  rvlijk = rvlamb(mod(isc/9,3)+1,mod(isc/3,3)+1,mod(isc,3)
383  & +1)
384  rvlkij = rvlamb(mod(isc,3)+1,mod(isc/9,3)+1,mod(isc/3,3)
385  & +1)
386  rvljki = rvlamb(mod(isc/3,3)+1,mod(isc,3)+1,mod(isc/9,3)
387  & +1)
388  dcmass = .false.
389  IF (idlam(lknt,1).EQ.-5.OR.idlam(lknt,2).EQ.-5
390  & .OR.idlam(lknt,3).EQ.-5) dcmass = .true.
391 C...Collect symmetry factors
392  IF (mod(isc/9,3).EQ.mod(isc/3,3).OR.mod(isc/3,3).eq
393  & .mod(isc,3).OR.mod(isc/9,3).EQ.mod(isc,3))
394  & rvlamc = 5d-1 * rvlamc
395 C...Resonance KF codes (1=I,2=J,3=K)
396  kfr(1) = idlam(lknt,1)-1
397  kfr(2) = 0
398  kfr(3) = 0
399 C...Calculate width.
400  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),
401  & idlam(lknt,3),xresi)
402 C...Resonance KF codes (1=I,2=J,3=K)
403  kfr(1) = 0
404  kfr(2) = idlam(lknt,2)-1
405  kfr(3) = 0
406 C...Calculate width.
407  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),
408  & idlam(lknt,3),xresj)
409 C...Resonance KF codes (1=I,2=J,3=K)
410  kfr(1) = 0
411  kfr(2) = 0
412  kfr(3) = idlam(lknt,3)-1
413 C...Calculate width.
414  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),
415  & idlam(lknt,3),xresk)
416 C...Resonance KF codes (1=I,2=J,3=K)
417  kfr(1) = idlam(lknt,1)-1
418  kfr(2) = idlam(lknt,2)-1
419  kfr(3) = 0
420 C...Calculate width.
421  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),
422  & idlam(lknt,3),xresij)
423  IF (abs(xresij/(xresi+xresj)-1.).GT.1d-4) THEN
424  xresij = xresi+xresj-xresij
425  ELSE
426  xresij = 0d0
427  ENDIF
428 C...Resonance KF codes (1=I,2=J,3=K)
429  kfr(1) = 0
430  kfr(2) = idlam(lknt,2)-1
431  kfr(3) = idlam(lknt,3)-1
432 C...Calculate width.
433  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),
434  & idlam(lknt,3),xresjk)
435  IF (abs(xresjk/(xresj+xresk)-1.).GT.1d-4) THEN
436  xresjk = xresj+xresk-xresjk
437  ELSE
438  xresjk = 0d0
439  ENDIF
440 C...Resonance KF codes (1=I,2=J,3=K)
441  kfr(1) = idlam(lknt,1)-1
442  kfr(2) = 0
443  kfr(3) = idlam(lknt,3)-1
444 C...Calculate width.
445  CALL pyrvgw(kfin,idlam(lknt,1),idlam(lknt,2),
446  & idlam(lknt,3),xresik)
447  IF (abs(xresik/(xresi+xresk)-1.).GT.1d-4) THEN
448  xresik = xresi+xresk-xresik
449  ELSE
450  xresik = 0d0
451  ENDIF
452 C...CALCULATE TOTAL WIDTH
453  xlam(lknt) =
454  & rvlijk**2 * xresi
455  & + rvljki**2 * xresj
456  & + rvlkij**2 * xresk
457  & + rvlijk*rvljki * xresij
458  & + rvlijk*rvlkij * xresik
459  & + rvljki*rvlkij * xresjk
460  xlam(lknt)=xlam(lknt)*rvlamc/((2.*paru(1)*rms(0))**3*32)
461 C...KINEMATICS CHECK
462  IF (xlam(lknt).EQ.0d0) THEN
463  lknt=lknt-1
464  ENDIF
465  ENDIF
466  190 CONTINUE
467  ENDIF
468  ENDIF
469  ENDIF
470 
471  RETURN
472  END