Reading QA thresholds from external file II
[u/mrichter/AliRoot.git] / LHAPDF / lhapdf5.3.1 / wrapzeus.f
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)
23 c      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 *
62 c===========================================================================
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)
70 c      print *,q2,alfas,qlam,nf,iflag
71       return
72 *
73 c===========================================================================
74       entry ZEUSread(nset)
75          read(1,*) gridname,nx,xmin,xmax,nq,qmin,qmax
76       return
77 *
78 c===========================================================================
79       entry ZEUSinit(nset,Eorder,Q2fit)
80 c      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
94 c--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
105 C 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       
114 c--qcdnum initialisation
115       CALL QNINIT
116 c--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     
129 c      Q2C=1.8225
130       Q2C=mc2
131 c      Q2B=18.49
132       Q2B=mb2 
133 c      print  *,q2c,q2b    
134
135       IF (Q0.LT.Q2C) THEN
136         NACT=3
137       ELSE
138         NACT=4
139       ENDIF
140 c--this merely defines nact where we startevolution
141 c--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)
151 c        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)
162 c      qcdnum grid not my grid
163
164 c      CALL GRXLIM(120,97D-8)
165       CALL GRXLIM(nx,xmin)
166
167 c      CALL GRQLIM(61,29D-2,200D3) 
168       CALL GRQLIM(nq,qmin,qmax) 
169
170 C--   Get final grid definitions and grid indices of Q0, Q2C and Q2B
171
172       CALL GRGIVE(NXGRI,XMI,XMA,NQGRI,QMI,QMA)
173 c      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)
177 C--   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
187 C--   Compute weights and dump, or read in
188 c
189 c      IF (READIN) THEN 
190 c        OPEN(UNIT=24,FILE='weights.dat',FORM='UNFORMATTED',
191 c     .                                  STATUS='UNKNOWN')
192 c        CALL QNREAD(24,ISTOP,IERR)
193 c      ELSE
194 c        CALL QNFILW(0,0)
195 c        IF (HEAVY) THEN
196 c          OPEN(UNIT=24,FILE='weights.dat',FORM='UNFORMATTED',
197 c     .                                    STATUS='UNKNOWN')
198 c          CALL QNDUMP(24)
199 c        ENDIF
200 c      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
236 C--   Apply cuts to grid
237 c--taking away the s cut at 600d0
238       CALL GRCUTS(-1D0,-1D0,-1D0,-1D0)
239
240
241
242
243 C--   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  
252 c       ZM=91.187D0
253       imem=0
254 c      print *,imem
255 c -- only need call to listPDF here to get alphas
256       call listPDF(nset,imem,xval)
257 c      print *,xval
258         AS=XVAL(1)
259 c        AS=0.118d0
260       CALL QNRSET('ALFQ0',ZM*ZM)
261       CALL QNRSET('ALFAS',AS)
262
263 c      ZM2=ZM*ZM
264       ALPHAS=QNALFA(ZM2)
265 c      WRITE(*,*)'ALPHAS AT Mz2',ALPHAS
266
267 C--   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  
277 C--   Book linear combinations for proton for f = 3,4,5 flavours
278
279 c--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) 
326 c---
327
328       return
329 *
330 c==========================================================================
331       entry ZEUSpdf(nset)
332 c      ZM = 91.187d0
333 c      zm2 = zm*zm
334
335 c      ALPHAS=QNALFA(ZM2)
336       
337       
338 c      imem=mem
339       call getnmem(nset,imem)    
340       call listPDF(nset,imem,xval)
341 c      print *,nset,imem
342 c      print *,xval
343 c      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)
368 c      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)
374 C--   Input quark distributions at Q2 = Q02 GeV2
375 C-- WRITE IN ACTUAL VALUES TO SAVE USING dgamma
376 C       UN=2./AREA(UA-1,UB,UE,UC)
377 C       DN=1./AREA(DA-1,DB,DE,DC)
378       UN = XVAL(2)
379       DN=XVAL(6)
380 c       UBDBN=DLN/AREA(DLA-1,DLB,DLE,DLC)
381 c      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
393 C        GN=(1-UN*AREA(UA,UB,UE,UC)-
394 C     .        DN*AREA(DA,DB,DE,DC)-
395 C     .        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
401 c        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
405 c        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
424 C--THINGS ARE FINE FOR HEAVY SO DO IT
425       IF (HEAVY) THEN
426
427 C--     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
438 C--     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         
444 C--     Be more careful with plus distributions
445
446         IF (NACT.EQ.3) THEN
447 C--THINGS ARE ALSO FINE IF 1Q0 IS BELOW 1QC THEN CLEARLY CSEA=0. IS OK
448 C--SITUATION CD BE 1<Q0<Q2C<Q2B ETC
449 C--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)
453 C--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)
461 C--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
468 C--1<1QC<1Q0<1QB<ETC
469 C--NOW WE NEED A RETHINK OF THE INITIAL CONDITIONS
470 C--FIRST DEAL WITH CPLUS TRUNING ON AT IQC
471 C--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)
476 C--REDIFINE THE PLUS DUSTNS TAKING INTO ACCOUNT CPLUS
477           CALL QNPNUL('SPLUS')
478           CALL QNPNUL('UPLUS')
479           CALL QNPNUL('DPLUS')
480 C--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
500 C-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)
504 C--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
517 C--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)
523 C--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
530 cC--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  
536 c
537       call getnset(iset)
538       call save_pdfqcd(iset)
539
540       return
541 *
542       end
543 C------------------------------------------------------------------------ 
544 c      DOUBLE PRECISION FUNCTION AREA(A,B,E,C)
545 C------------------------------------------------------------------------ 
546 c
547 c      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
548
549 c      IF (A.LE.-0.99.OR.B.LE.-0.99) THEN
550 c        AREA=1D6
551 c        RETURN
552 c      ENDIF
553 c
554 c      AR1=(DGAMMA_LHA(A+1)*DGAMMA_LHA(B+1))/DGAMMA_LHA(A+B+2)
555 c      AR2=E*(DGAMMA_LHA(A+1.5)*DGAMMA_LHA(B+1))/DGAMMA_LHA(A+B+2.5)
556 c      AR3=C*(DGAMMA_LHA(A+2)*DGAMMA_LHA(B+1))/DGAMMA_LHA(A+B+3)
557 c      AREA=AR1+AR2+AR3
558 c     
559 c      IF (AREA.LE.1D-6) AREA=1D-6
560 c
561 c      RETURN
562 c      END
563 c
564 C------------------------------------------------------------------------ 
565       DOUBLE PRECISION FUNCTION FF_LHA(X,A,B,E,C)
566 C------------------------------------------------------------------------ 
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
574 C------------------------------------------------------------------------ 
575       SUBROUTINE DVZERO (A,N)
576 C
577 C CERN PROGLIB# F121    VZERO           .VERSION KERNFOR  4.40  940929
578 C ORIG. 01/07/71, modif. 24/05/87 to set integer zero
579 C                 modif. 25/05/94 to depend on QINTZERO
580 C
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