program gluon_grid
c          implicit double precision (A-H,O-Z)
          common/ZEUS/val(0:160,0:160)
          dimension f2(0:160,0:160),f2up(0:160,0:160),f2dn(0:160,0:160)
 
c---read in gluon grid at intialisation
          open(unit=05,file='gluon_tr_tot.dat',status='old')
          do i=0,160
          do j=0,160
          read(05,53)f2(j,i),f2up(j,i),f2dn(j,i)
          enddo
          enddo
          close(05)

c-- read the central vlaue into val array to do interpolation
          do i=0,160
          do j=0,160
          val(j,i)=f2(j,i)
          enddo
          enddo
c---------------------------------
c--input what x,q2 you want here
          x=1.e-3
          q2=20.
c--get result
          gq=znlo(x,q2)
          write(*,*)'gluon value, x,q2',gq,x,q2

c---carry on if you want errors
          do i=0,160
          do j=0,160
          val(j,i)=f2up(j,i)
          enddo
          enddo
          gqup=znlo(x,q2)

          do i=0,160
          do j=0,160

          val(j,i)=f2dn(j,i)
          enddo
          enddo
          gqdn=znlo(x,q2)

          write(*,*)'gluon central,up,down  x,q2',gq,gqup,gqdn,x,q2

 53       format(3(e12.4,1x))
          stop
          end
C     ========================
      REAL FUNCTION ZNLO(X,Q2)
C     ========================
C
C     Provide a simple interface for plotting ZEUS NLO fit results
C
C     INPUT:  X      Bjorken x
C             Q2     Q2
C            
C     OUTPUT: ZNLO   Value for PDF/Xsect loaded
C
C     Alex Tapper 18/09/03 based on routines by E. Tassi & A. Cooper-Sarkar
 
      COMMON /ZEUS/F2(0:160,0:160)

C --- Return 1 if x or Q2 out of bounds
      IF (Q2.LT.0.39.OR.X.LT.1E-6.OR.X.GT.1.OR.Q2.GT.101200.) THEN
        ZNLO=1.
        RETURN
      ENDIF

C --- Calculate grid points
      DO I=0,159
        IF(I.LE.71)THEN
          Q2G1=10**(5.*I/120.)-0.7
          Q2G2=10**(5.*(I+1)/120.)-0.7
        ELSEIF(I.EQ.72)THEN
          Q2G1=10**(5.*I/120.)-0.7
          Q2G2=10**(2.*(I-71)/88. +3.)
        ELSE
          Q2G1=10**(2.*(I-72)/88. +3.)
          Q2G2=10**(2.*(I-71)/88. +3.)
        ENDIF
        IF (Q2.GE.Q2G1.AND.Q2.LT.Q2G2) THEN
          IQ2L=I+1
          IQ2H=I+2
          TQ=(Q2-Q2G1)/(Q2G2-Q2G1)
        ENDIF
      ENDDO

      DO I=0,159
      IF(I.LE.79)THEN
        XG1=10**(6.*I/120.-6.)
        XG2=10**(6.*(I+1)/120.-6.)
      ELSEIF(I.EQ.80)THEN
        XG1=10**(6.*I/120.-6.)
        XG2=10**(2./80.*(I-79)-2.)
      ELSE
        XG1=10**(2./80.*(I-80)-2.)
        XG2=10**(2./80.*(I-79)-2.)
      ENDIF
        IF (X.GE.XG1.AND.X.LT.XG2) THEN
          IXL=I+1
          IXH=I+2
          TX=(X-XG1)/(XG2-XG1)
        ENDIF
      ENDDO
 
      ZF2H=(1-TX)*F2(IXL,IQ2H)+TX*F2(IXH,IQ2H)
      ZF2L=(1-TX)*F2(IXL,IQ2L)+TX*F2(IXH,IQ2L)
      ZNLO=(1-TQ)*ZF2L+TQ*ZF2H

      RETURN 
      END