]> git.uio.no Git - u/mrichter/AliRoot.git/blame - LHAPDF/lhapdf5.3.1/wrapzeus.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
pdfs for LHAPDF version 5.3.1
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapzeus.f
CommitLineData
4e9e3152 1 subroutine ZEUSevolve(x,Q,pdf)
2 implicit double precision (a-h,o-z)
3 include 'parmsetup.inc'
4 character*64 gridname
5 character*16 name(nmxset)
6 integer nmem(nmxset),ndef(nmxset),mem
7 common/NAME/name,nmem,ndef,mem
8 integer nset,iset,isetlast
9 data isetlast/-1/
10 integer Eorder
11 real*8 mc,mc2,mb,mb2,mt,mt2
12 real*8 f(-6:6),pdf(-6:6)
13 integer qnerr
14 parameter(nstartp=7)
15 DIMENSION QSP(NSTARTP)
16 DATA QSP/10.,20.,30.,40.,50.,80.,100./
17 real*8 xval(45)
18 logical HEAVY,VFN
19 real*8 pwgt(20)
20 save
21*
22 call getnset(iset)
23c print *,'iset=',iset,' now calling get_pdfqcd'
24 if(iset.ne.isetlast) then
25 call get_pdfqcd(iset)
26 isetlast = iset
27 endif
28*
29 Q2=Q*Q
30 UCENT=QPDFXQ('UPVAL',X,Q2,IFAIL)
31 DCENT=QPDFXQ('DNVAL',X,Q2,IFAIL)
32 GCENT=QPDFXQ('GLUON',X,Q2,IFAIL)
33 UBCEN=QPDFXQ('UB',X,Q2,IFAIL)
34 DBCEN=QPDFXQ('DB',X,Q2,IFAIL)
35 STCEN=QPDFXQ('SB',X,Q2,IFAIL)
36 IF(Q2.GE.Q2C)THEN
37 CHCEN=QPDFXQ('CB',X,Q2,IFAIL)
38 ELSE
39 CHCEN=0.0
40 ENDIF
41 IF(Q2.GE.Q2B)THEN
42 BTCEN=QPDFXQ('BB',X,Q2,IFAIL)
43 ELSE
44 BTCEN=0.0
45 ENDIF
46 pdf(1)=dcent+dbcen
47 pdf(2)=ucent+ubcen
48 pdf(3)=stcen
49 pdf(4)=chcen
50 pdf(5)=btcen
51 pdf(6)=0d0
52 pdf(0)=gcent
53 pdf(-1)=dbcen
54 pdf(-2)=ubcen
55 pdf(-3)=stcen
56 pdf(-4)=chcen
57 pdf(-5)=btcen
58 pdf(-6)=0d0
59*
60 return
61*
62c===========================================================================
63 entry ZEUSalfa(alfas,Q)
64 Q2=Q*Q
65 nf=6
66 if (Q2.lt.mt2) nf=5
67 if (Q2.lt.mb2) nf=4
68 if (Q2.lt.mc2) nf=3
69 alfas=QALFAS(Q2,Qlam,nf,iflag)
70c print *,q2,alfas,qlam,nf,iflag
71 return
72*
73c===========================================================================
74 entry ZEUSread(nset)
75 read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
76 return
77*
78c===========================================================================
79 entry ZEUSinit(nset,Eorder,Q2fit)
80c print *,name(nset)
81 if(name(nset).eq.'QCDNUM_ZEUS_TR') then
82 HEAVY = .FALSE.
83 VFN = .TRUE.
84 else if(name(nset).eq.'QCDNUM_ZEUS_FF') then
85 HEAVY = .TRUE.
86 VFN = .FALSE.
87 else if(name(nset).eq.'QCDNUM_ZEUS_ZM') then
88 HEAVY = .FALSE.
89 VFN = .FALSE.
90 else
91 print *,'name/scheme not recognized'
92 stop 1
93 endif
94c--try 3 way logic ffn/zm-vfn/rt-vfn
95 IF(HEAVY)THEN
96 IVFN=1
97 ELSE
98 IVFN=0
99 ENDIF
100 IF(VFN)THEN
101 IVFN=IVFN+2
102 ELSE
103 IVFN=IVFN
104 ENDIF
105C IVFN=0 IS ZM-VFN, 1 IS FFN,2 IS RT-VFN, 3 IS NOT ALLOWED
106
107 IF(IVFN.EQ.3)THEN
108 WRITE(*,*)'IVFN=3 SO STOP',IVFN
109 STOP
110 ENDIF
111
112
113
114c--qcdnum initialisation
115 CALL QNINIT
116c--se thresholds
117 Q0=Q2fit
118 ZM=91.187D0
119 ZM2=ZM*ZM
120 ALPHAS=QNALFA(ZM2)
121
122 call getQmassM(nset,4,mc)
123 mc2=mc*mc
124 call getQmassM(nset,5,mb)
125 mb2=mb*mb
126 call getQmassM(nset,6,mt)
127 mt2=mt*mt
128
129c Q2C=1.8225
130 Q2C=mc2
131c Q2B=18.49
132 Q2B=mb2
133c print *,q2c,q2b
134
135 IF (Q0.LT.Q2C) THEN
136 NACT=3
137 ELSE
138 NACT=4
139 ENDIF
140c--this merely defines nact where we startevolution
141c--namely at q0
142 IF (HEAVY) NACT=3
143
144 CALL QNRSET('MCSTF',SQRT(Q2C))
145 CALL QNRSET('MBSTF',SQRT(Q2B))
146 CALL QNRSET('MCALF',SQRT(Q2C))
147 CALL QNRSET('MBALF',SQRT(Q2B))
148
149 IF (HEAVY) THEN
150 CALL QTHRES(1D10,2D10)
151c CALL QTHRES(1D6,2D6)
152 ELSE
153 CALL QTHRES(Q2C,Q2B)
154 ENDIF
155
156 DO I=1,NSTARTP
157 CALL GRQINP(QSP(I),1)
158 ENDDO
159 CALL GRQINP(Q0,1)
160 CALL GRQINP(Q2C,1)
161 CALL GRQINP(Q2B,1)
162c qcdnum grid not my grid
163
164c CALL GRXLIM(120,97D-8)
165 CALL GRXLIM(nx,xmin)
166
167c CALL GRQLIM(61,29D-2,200D3)
168 CALL GRQLIM(nq,qmin,qmax)
169
170C-- Get final grid definitions and grid indices of Q0, Q2C and Q2B
171
172 CALL GRGIVE(NXGRI,XMI,XMA,NQGRI,QMI,QMA)
173c WRITE(*,*)'NX,XL,XH,NQ,QL,QH',NXGRI,XMI,XMA,NQGRI,QMI,QMA
174 IQ0 = IQFROMQ(Q0)
175 IQC = IQFROMQ(Q2C)
176 IQB = IQFROMQ(Q2B)
177C-- Allow for heavy weights
178
179 IF (HEAVY) THEN
180 CALL QNLSET('WTF2C',.TRUE.)
181 CALL QNLSET('WTF2B',.TRUE.)
182 CALL QNLSET('CLOWQ',.FALSE.)
183 CALL QNLSET('WTFLC',.TRUE.)
184 CALL QNLSET('WTFLB',.TRUE.)
185 ENDIF
186
187C-- Compute weights and dump, or read in
188c
189c IF (READIN) THEN
190c OPEN(UNIT=24,FILE='weights.dat',FORM='UNFORMATTED',
191c . STATUS='UNKNOWN')
192c CALL QNREAD(24,ISTOP,IERR)
193c ELSE
194c CALL QNFILW(0,0)
195c IF (HEAVY) THEN
196c OPEN(UNIT=24,FILE='weights.dat',FORM='UNFORMATTED',
197c . STATUS='UNKNOWN')
198c CALL QNDUMP(24)
199c ENDIF
200c ENDIF
201
202
203 if (index(gridname,'none').eq.1) then
204 call qnfilw(0,0)
205 else
206 qnerr=-1
207 open(unit=2,status='old',file=gridname,
208 . form='unformatted',err=1)
209 call QNREAD(2,1,qnerr)
210 1 close(2)
211 if (qnerr.ne.0) then
212 write(*,*) 'Grid file problem: ',gridname
213 if (qnerr.lt.0) then
214 write(*,*) 'Grid file does not exist'
215 write(*,*) 'Calculating and creating grid file'
216 call qnfilw(0,0)
217 open(unit=2,status='unknown',file=gridname,
218 . form='unformatted')
219 call QNDUMP(2)
220 close(2)
221 else
222 write(*,*) 'Existing grid file is inconsistent'
223 if (qnerr.eq.1)
224 . write(*,*) 'Defined grid different'
225 if (qnerr.eq.2)
226 . write(*,*) 'Heavy quark weight table different'
227 if (qnerr.eq.3)
228 . write(*,*) 'Charm mass different'
229 if (qnerr.eq.4)
230 . write(*,*) 'Bottom mass different'
231 stop
232 endif
233 endif
234 endif
235
236C-- Apply cuts to grid
237c--taking away the s cut at 600d0
238 CALL GRCUTS(-1D0,-1D0,-1D0,-1D0)
239
240
241
242
243C-- Choose renormalisation and factorisation scales
244
245 CALL QNRSET('AAAR2',1D0) ! renormalisation
246 CALL QNRSET('BBBR2',0D0)
247 CALL QNRSET('AAM2L',1D0) ! factorisation (light)
248 CALL QNRSET('BBM2L',0D0)
249 CALL QNRSET('AAM2H',1D0) ! factorisation (heavy)
250 CALL QNRSET('BBM2H',0D0)
251
252c ZM=91.187D0
253 imem=0
254c print *,imem
255c -- only need call to listPDF here to get alphas
256 call listPDF(nset,imem,xval)
257c print *,xval
258 AS=XVAL(1)
259c AS=0.118d0
260 CALL QNRSET('ALFQ0',ZM*ZM)
261 CALL QNRSET('ALFAS',AS)
262
263c ZM2=ZM*ZM
264 ALPHAS=QNALFA(ZM2)
265c WRITE(*,*)'ALPHAS AT Mz2',ALPHAS
266
267C-- Book non-singlet distributions
268
269 CALL QNBOOK(2,'UPLUS')
270 CALL QNBOOK(3,'DPLUS')
271 CALL QNBOOK(4,'SPLUS')
272 CALL QNBOOK(5,'CPLUS')
273 CALL QNBOOK(6,'BPLUS')
274 CALL QNBOOK(7,'UPVAL')
275 CALL QNBOOK(8,'DNVAL')
276
277C-- Book linear combinations for proton for f = 3,4,5 flavours
278
279c--define some quark pdfs
280 CALL dVZERO(PWGT,20)
281 PWGT(2) = 0.5
282
283 PWGT(7) = -0.5
284
285 PWGT(1) = 0.5/3.
286 CALL QNLINC(17,'UB',3,PWGT)
287 PWGT(1) = 0.5/4.
288 CALL QNLINC(17,'UB',4,PWGT)
289 PWGT(1) = 0.5/5.
290 CALL QNLINC(17,'UB',5,PWGT)
291 CALL dVZERO(PWGT,20)
292
293 PWGT(4) = 0.5
294 PWGT(1) = 0.5/3.
295 CALL QNLINC(18,'SB',3,PWGT)
296 PWGT(1) = 0.5/4.
297 CALL QNLINC(18,'SB',4,PWGT)
298 PWGT(1) = 0.5/5.
299 CALL QNLINC(18,'SB',5,PWGT)
300 CALL dVZERO(PWGT,20)
301 CALL QNLINC(19,'CB',3,PWGT)
302 PWGT(5) = 0.5
303 PWGT(1) = 0.5/4.
304
305 CALL QNLINC(19,'CB',4,PWGT)
306 PWGT(1) = 0.5/5.
307 CALL QNLINC(19,'CB',5,PWGT)
308 CALL dVZERO(PWGT,20)
309 PWGT(3) = 0.5
310
311 PWGT(8) = -0.5
312
313 PWGT(1) = 0.5/3.
314 CALL QNLINC(20,'DB',3,PWGT)
315 PWGT(1) = 0.5/4.
316 CALL QNLINC(20,'DB',4,PWGT)
317 PWGT(1) = 0.5/5.
318 CALL QNLINC(20,'DB',5,PWGT)
319 CALL dVZERO(PWGT,20)
320 CALL QNLINC(21,'BB',3,PWGT)
321 CALL QNLINC(21,'BB',4,PWGT)
322 PWGT(6) = 0.5
323
324 PWGT(1) = 0.5/5.
325 CALL QNLINC(21,'BB',5,PWGT)
326c---
327
328 return
329*
330c==========================================================================
331 entry ZEUSpdf(nset)
332c ZM = 91.187d0
333c zm2 = zm*zm
334
335c ALPHAS=QNALFA(ZM2)
336
337
338c imem=mem
339 call getnmem(nset,imem)
340 call listPDF(nset,imem,xval)
341c print *,nset,imem
342c print *,xval
343c print *,imem,xval
344
345 UA=XVAL(3)
346 UB=XVAL(4)
347 UE=0.0d0
348 UC=XVAL(5)
349
350 DA=XVAL(7)
351 DB=XVAL(8)
352 DE=0.0d0
353 DC=XVAL(9)
354
355 GA=XVAL(11)
356 GB=XVAL(12)
357 GE=0.0d0
358 GC=XVAL(13)
359
360 SN=XVAL(14)
361 SA=XVAL(15)
362 SB=XVAL(16)
363 SE=0.0d0
364 SC=XVAL(17)
365
366 DLN=XVAL(18)
367 DLA=XVAL(19)
368c DLB=XVAL(20)
369 DLB=XVAL(16)+2.0d0
370 DLE=0.0d0
371 DLC=XVAL(21)
372 AS=XVAL(1)
373 CALL QNRSET('ALFAS',AS)
374C-- Input quark distributions at Q2 = Q02 GeV2
375C-- WRITE IN ACTUAL VALUES TO SAVE USING dgamma
376C UN=2./AREA(UA-1,UB,UE,UC)
377C DN=1./AREA(DA-1,DB,DE,DC)
378 UN = XVAL(2)
379 DN=XVAL(6)
380c UBDBN=DLN/AREA(DLA-1,DLB,DLE,DLC)
381c call grgive(nxgri,xxmin,xxmax,nqgri,qqmin,qqmax)
382 nxgri=nx
383 DO IX = 1,NXGRI
384 xX = XFROMIX(IX)
385 UPVAL=UN*FF_LHA(Xx,UA,UB,UE,UC)
386
387
388 DNVAL=DN*FF_LHA(Xx,DA,DB,DE,DC)
389
390
391 SEA=SN*FF_LHA(Xx,SA,SB,SE,SC)
392
393C GN=(1-UN*AREA(UA,UB,UE,UC)-
394C . DN*AREA(DA,DB,DE,DC)-
395C . SN*AREA(SA,SB,SE,SC))/AREA(GA,GB,GE,GC)
396 GN=XVAL(10)
397 GLUON=GN*FF_LHA(Xx,GA,GB,GE,GC)
398
399
400
401c UMD=UBDBN*FF_LHA(X,DLA,DLB,DLE,DLC)
402 UMD=DLN*FF_LHA(Xx,DLA,DLB,DLE,DLC)
403 SINGL=UPVAL+DNVAL+SEA
404
405c print *,un,dn,sn,gn,dln
406
407 CSEA=0.0
408 SSEA=0.2*SEA
409 USEA=(SEA-SSEA-CSEA)/2.-UMD
410 DSEA=(SEA-SSEA-CSEA)/2.+UMD
411 UPLUS=UPVAL+USEA-1./NACT*SINGL
412 DPLUS=DNVAL+DSEA-1./NACT*SINGL
413 SPLUS=SSEA-1./NACT*SINGL
414
415
416 CALL QNPSET('GLUON',IX,IQ0,GLUON)
417 CALL QNPSET('SINGL',IX,IQ0,SINGL)
418 CALL QNPSET('UPLUS',IX,IQ0,UPLUS)
419 CALL QNPSET('DPLUS',IX,IQ0,DPLUS)
420 CALL QNPSET('SPLUS',IX,IQ0,SPLUS)
421 CALL QNPSET('UPVAL',IX,IQ0,UPVAL)
422 CALL QNPSET('DNVAL',IX,IQ0,DNVAL)
423 ENDDO
424C--THINGS ARE FINE FOR HEAVY SO DO IT
425 IF (HEAVY) THEN
426
427C-- Evolve over whole grid
428
429 CALL EVOLSG(IQ0,1,NQGRI)
430 CALL EVPLUS('UPLUS',IQ0,1,NQGRI)
431 CALL EVPLUS('DPLUS',IQ0,1,NQGRI)
432 CALL EVPLUS('SPLUS',IQ0,1,NQGRI)
433 CALL EVOLNM('UPVAL',IQ0,1,NQGRI)
434 CALL EVOLNM('DNVAL',IQ0,1,NQGRI)
435
436 ELSE
437
438C-- Evolve gluon, singlet and valence over whole grid
439
440 CALL EVOLSG(IQ0,1,NQGRI)
441 CALL EVOLNM('UPVAL',IQ0,1,NQGRI)
442 CALL EVOLNM('DNVAL',IQ0,1,NQGRI)
443
444C-- Be more careful with plus distributions
445
446 IF (NACT.EQ.3) THEN
447C--THINGS ARE ALSO FINE IF 1Q0 IS BELOW 1QC THEN CLEARLY CSEA=0. IS OK
448C--SITUATION CD BE 1<Q0<Q2C<Q2B ETC
449C--GO DOWN QO TO 1 UP Q0 TO Q2C
450 CALL EVPLUS('UPLUS',IQ0,1,IQC)
451 CALL EVPLUS('DPLUS',IQ0,1,IQC)
452 CALL EVPLUS('SPLUS',IQ0,1,IQC)
453C--DEAL WITH CHARM THRESH
454 FACTOR=1./3.-1./4.
455 CALL QADDSI('UPLUS',IQC,FACTOR)
456 CALL QADDSI('DPLUS',IQC,FACTOR)
457 CALL QADDSI('SPLUS',IQC,FACTOR)
458 CALL QNPNUL('CPLUS')
459 FACTOR=-1/4.
460 CALL QADDSI('CPLUS',IQC,FACTOR)
461C--GO UP Q2C TO Q2B
462 CALL EVPLUS('UPLUS',IQC,IQC,IQB)
463 CALL EVPLUS('DPLUS',IQC,IQC,IQB)
464 CALL EVPLUS('SPLUS',IQC,IQC,IQB)
465 CALL EVPLUS('CPLUS',IQC,IQC,IQB)
466
467 ELSEIF(NACT.EQ.4)THEN
468C--1<1QC<1Q0<1QB<ETC
469C--NOW WE NEED A RETHINK OF THE INITIAL CONDITIONS
470C--FIRST DEAL WITH CPLUS TRUNING ON AT IQC
471C--AND GOING UP TO IQB (MIDDLE AGAIN)
472 CALL QNPNUL('CPLUS')
473 FACTOR=-1/4.
474 CALL QADDSI('CPLUS',IQC,FACTOR)
475 CALL EVPLUS('CPLUS',IQC,IQC,IQB)
476C--REDIFINE THE PLUS DUSTNS TAKING INTO ACCOUNT CPLUS
477 CALL QNPNUL('SPLUS')
478 CALL QNPNUL('UPLUS')
479 CALL QNPNUL('DPLUS')
480C--NOW YOU NEED A DO LOOP AGAIN
481 DO IX = 1,NXGRI
482 Xx = XFROMIX(IX)
483 UPVAL=UN*FF_LHA(Xx,UA,UB,UE,UC)
484
485 DNVAL=DN*FF_LHA(Xx,DA,DB,DE,DC)
486 SEA=SN*FF_LHA(Xx,SA,SB,SE,SC)
487 SINGL=UPVAL+DNVAL+SEA
488 UMD=DLN*FF_LHA(Xx,DLA,DLB,DLE,DLC)
489 SSEA=0.2*SEA
490 CSEA=QPDFIJ('CPLUS',IX,IQ0,IFL) + 1./NACT*SINGL
491 USEA=(SEA-SSEA-CSEA)/2.-UMD
492 DSEA=(SEA-SSEA-CSEA)/2.+UMD
493 UPLUS=UPVAL+USEA-1./NACT*SINGL
494 DPLUS=DNVAL+DSEA-1./NACT*SINGL
495 SPLUS=SSEA-1./NACT*SINGL
496 CALL QNPSET('UPLUS',IX,IQ0,UPLUS)
497 CALL QNPSET('DPLUS',IX,IQ0,DPLUS)
498 CALL QNPSET('SPLUS',IX,IQ0,SPLUS)
499 ENDDO
500C-NOW DO MIDDLE FOR UPLUS,DPLUS,SPLUS
501 CALL EVPLUS('UPLUS',IQ0,IQC,IQB)
502 CALL EVPLUS('DPLUS',IQ0,IQC,IQB)
503 CALL EVPLUS('SPLUS',IQ0,IQC,IQB)
504C--THEN GO DOWN IQC TO 1
505 IF(IQC.GT.1)THEN
506 FACTOR=1./4.-1./3.
507 CALL QADDSI('UPLUS',IQC,FACTOR)
508 CALL QADDSI('DPLUS',IQC,FACTOR)
509 CALL QADDSI('SPLUS',IQC,FACTOR)
510 CALL EVPLUS('UPLUS',IQC,1,IQC)
511 CALL EVPLUS('DPLUS',IQC,1,IQC)
512
513 CALL EVPLUS('SPLUS',IQC,1,IQC)
514 ENDIF
515
516 ENDIF
517C--THEN DEAL WITH B THRESHOLD FOR ALL4
518 FACTOR=1./4.-1./5.
519 CALL QADDSI('UPLUS',IQB,FACTOR)
520 CALL QADDSI('DPLUS',IQB,FACTOR)
521 CALL QADDSI('SPLUS',IQB,FACTOR)
522 CALL QADDSI('CPLUS',IQB,FACTOR)
523C--THEN GO UP
524 IF(IQB.LT.NQGRI)THEN
525 CALL EVPLUS('UPLUS',IQB,IQB,NQGRI)
526 CALL EVPLUS('DPLUS',IQB,IQB,NQGRI)
527 CALL EVPLUS('SPLUS',IQB,IQB,NQGRI)
528 CALL EVPLUS('CPLUS',IQB,IQB,NQGRI)
529 ENDIF
530cC--THEN DEAL WITH B TURNING ON AT IQB AND GO UP
531 CALL QNPNUL('BPLUS')
532 FACTOR=-1./5.
533 CALL QADDSI('BPLUS',IQB,FACTOR)
534 CALL EVPLUS('BPLUS',IQB,IQB,NQGRI)
535 ENDIF
536c
537 call getnset(iset)
538 call save_pdfqcd(iset)
539
540 return
541*
542 end
543C------------------------------------------------------------------------
544c DOUBLE PRECISION FUNCTION AREA(A,B,E,C)
545C------------------------------------------------------------------------
546c
547c IMPLICIT DOUBLE PRECISION (A-H,O-Z)
548c
549c IF (A.LE.-0.99.OR.B.LE.-0.99) THEN
550c AREA=1D6
551c RETURN
552c ENDIF
553c
554c AR1=(DGAMMA_LHA(A+1)*DGAMMA_LHA(B+1))/DGAMMA_LHA(A+B+2)
555c AR2=E*(DGAMMA_LHA(A+1.5)*DGAMMA_LHA(B+1))/DGAMMA_LHA(A+B+2.5)
556c AR3=C*(DGAMMA_LHA(A+2)*DGAMMA_LHA(B+1))/DGAMMA_LHA(A+B+3)
557c AREA=AR1+AR2+AR3
558c
559c IF (AREA.LE.1D-6) AREA=1D-6
560c
561c RETURN
562c END
563c
564C------------------------------------------------------------------------
565 DOUBLE PRECISION FUNCTION FF_LHA(X,A,B,E,C)
566C------------------------------------------------------------------------
567
568 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
569
570 FF_LHA=X**A*(1D0-X)**B*(1+E*SQRT(X)+C*X)
571
572 RETURN
573 END
574C------------------------------------------------------------------------
575 SUBROUTINE DVZERO (A,N)
576C
577C CERN PROGLIB# F121 VZERO .VERSION KERNFOR 4.40 940929
578C ORIG. 01/07/71, modif. 24/05/87 to set integer zero
579C modif. 25/05/94 to depend on QINTZERO
580C
581 implicit real*8 (a-h,o-z)
582 DIMENSION A(*)
583 IF (N.LE.0) RETURN
584 DO 9 I= 1,N
585 9 A(I)= 0d0
586 RETURN
587 END
588
589
590
591