c       Example of the use of the subroutine  
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        CALL INI

        MODE=1
        X=0.2
        Q2=2.0
        CALL POLFIT(MODE,X,Q2,DUV,DDV,DUBAR,DDBAR,DSTR,DGLU,
     #              G1P,G1N)
c       (remember it gives x times the pdf or x times g1!!!) 

        WRITE(6,*) X,Q2,DUV,DDV,DUBAR,DDBAR
        END


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C									     C
C       POLARIZED PARTON DISTRIBUTIONS (NLO)			             C
C       FROM  PRD71, 094018 (2005)                                           C
C       hep-ph/0504155                                                       C
C	D. de Florian G.A. Navarro and R. Sassot.                            C
C 									     C
C	MODE: 1 SET 1 NLO (MSbar) (based on KRETZER ffs)		     C
C             2 SET 2 NLO (MSbar) (based on KKP     ffs) 		     C
C             3 SET 1  LO         (based on KRETZER ffs) 		     C
C             4 SET 2  LO         (based on KKP     ffs) 		     C
C		                                                             C
C	             Q2=Q^2						     C
C                    DUV :    X * U VALENCE DISTRIBUTION		     C
C                    DDV :    X * D VALENCE DISTRIBUTION 		     C
C                    DUBAR :  X * UBAR DISTRIBUTION   	                     C
C                    DDBAR :  X * DBAR DISTRIBUTION              	     C
C                    DSTR :   X * STRANGE DISTRIBUTION       		     C
C                    DGLU :   X * GLUON DISTRIBUTION 			     C
C                    G1P :    X * POLARIZED STRUCTURE FUNCTION (PROTON)      C
C                    G1N :    X * POLARIZED STRUCTURE FUNCTION(NEUTRON)      C
C									     C
C       REMEMBER: ALWAYS X*DISTRIBUTION !!!				     C  C									     C
C 	BEFORE CALLING THE SUBRUTINE `POLFIT` FOR THE FIRST TIME, THE        C
C       SUBROUTINE `INI` MUST BE CALLED (ONLY ONCE) TO READ THE GRIDS.       C
C              (CALL INI)	
C	RANGE OF VALIDITY OF THE INTERPOLATION:  			     C
C					10**(-4)< X < 0.9		     C
C					1 < Q**2 < 5*10**4                   C
C									     C
C       IN CASE OF PROBLEMS, DOUBTS, ETC, PLEASE REPORT TO		     C
C        dflorian@df.uba.ar     					     C
C        sassot@df.uba.ar         					     C
C									     C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

       SUBROUTINE POLFIT(MODE,X,Q2,DUV,DDV,DUBAR,DDBAR,DSTR,DGLU,
     #  G1P,G1N)
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	DIMENSION GU(6,26,76),GD(6,26,76),GUB(6,26,76),GDB(6,26,76),
     #               GS(6,26,76), GG(6,26,76), GP(6,26,76), GN(6,26,76)	
	DIMENSION DUVG(76,26), DDVG(76,26), DUBARG(76,26),DDBARG(76,26),        
     #            DSTRG(76,26),DGLUG(76,26),GPG(76,26),GNG(76,26),XQ(2)
	COMMON/ GRID / GU,GD,GUB,GDB,GS,GG,GP,GN         

 		
        DO  K = 1, 26
        DO  J = 1, 76
        DUVG(J,K) = GU(MODE,K,J)  
        DDVG(J,K) = GD(MODE,K,J)  
        DUBARG(J,K) = GUB(MODE,K,J)   
        DDBARG(J,K) = GDB(MODE,K,J)   
        DSTRG(J,K)  = GS(MODE,K,J) 
        DGLUG(J,K)  = GG(MODE,K,J)  
        GPG(J,K)  = GP(MODE,K,J)  
        GNG(J,K)  = GN(MODE,K,J) 
  	END DO
	END DO 

  
        XQ(1) = DLOG(X)
        XQ(2) = DLOG(Q2)
	X3=(1.D0-X)**3.D0
	X4=(1.D0-X)**4.D0 
	X5=X**0.5D0
	X6=X**0.6D0
	X7=X**0.7D0
        DUV = PERINOLA(XQ,DUVG) * X3* X6*X
        DDV = PERINOLA(XQ,DDVG) * X4 * X7*X 
        DUBAR = PERINOLA(XQ,DUBARG) * X3 * X5*X
        DDBAR = PERINOLA(XQ,DDBARG) * X3 * X5*X
        DSTR = PERINOLA(XQ,DSTRG)  * X3 * X5*X
        DGLU = PERINOLA(XQ,DGLUG)  * X3 * X5*X
        G1P = PERINOLA(XQ,GPG)  * X3 * X5/X
        G1N = PERINOLA(XQ,GNG)  * X4 * X5/X
c        G1D=(G1P+G1N)*0.5D0*(1.D0-1.5D0*0.058D0)
        RETURN
        END



 

      FUNCTION PERINOLA(ARG,TABLE)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION ARG(2),NENT(2),ENT(102),TABLE(1976)
      DIMENSION D(7),NCOMB(7),IENT(7)
      COMMON/ XARRAY / ENT
       NARG=2
      NENT(1)=76
      NENT(2)=26
       KD=1
      M=1
      JA=1
         DO 5 I=1,NARG
      NCOMB(I)=1
      JB=JA-1+NENT(I)
         DO 2 J=JA,JB
      IF (ARG(I).LE.ENT(J)) GO TO 3
    2 CONTINUE
      J=JB
    3 IF (J.NE.JA) GO TO 4
      J=J+1
    4 JR=J-1
      D(I)=(ENT(J)-ARG(I))/(ENT(J)-ENT(JR))
      IENT(I)=J-JA
      KD=KD+IENT(I)*M
      M=M*NENT(I)
    5 JA=JB+1
      PERINOLA=0.D0
   10 FAC=1.D0
      IADR=KD
      IFADR=1
         DO 15 I=1,NARG
      IF (NCOMB(I).EQ.0) GO TO 12
      FAC=FAC*(1.D0-D(I))
      GO TO 15
   12 FAC=FAC*D(I)
      IADR=IADR-IFADR
   15 IFADR=IFADR*NENT(I)
      PERINOLA=PERINOLA+FAC*TABLE(IADR)
      IL=NARG
   40 IF (NCOMB(IL).EQ.0) GO TO 80
      NCOMB(IL)=0
      IF (IL.EQ.NARG) GO TO 10
      IL=IL+1
         DO 50  K=IL,NARG
   50 NCOMB(K)=1
      GO TO 10
   80 IL=IL-1
      IF(IL.NE.0) GO TO 40
      RETURN
      END



      SUBROUTINE INI
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
	DIMENSION XARRAY(102), GU(6,26,76), GD(6,26,76), GUB(6,26,76),
     #    GDB(6,26,76),GS(6,26,76),GG(6,26,76),GP(6,26,76),GN(6,26,76) 
	COMMON/ XARRAY / XARRAY
	COMMON/ GRID / GU,GD,GUB,GDB,GS,GG,GP,GN

       OPEN(UNIT=10,FILE='DNS1NLO.GRID',STATUS='OLD')
       OPEN(UNIT=11,FILE='DNS2NLO.GRID',STATUS='OLD')
       OPEN(UNIT=12,FILE='DNS1LO.GRID',STATUS='OLD')
       OPEN(UNIT=13,FILE='DNS2LO.GRID',STATUS='OLD')

 
 	DO IFILE=10,13
        DO  K = 1, 76 
        DO  J = 1, 26
	I=IFILE-9
        READ(IFILE,40) GU(I,J,K), GD(I,J,K), GUB(I,J,K),GDB(I,J,K),  
     #                 GS(I,J,K), GG(I,J,K), GP(I,J,K), GN(I,J,K)
 	END DO
	END DO
	CLOSE(IFILE)
	END DO
  40    FORMAT (8(1PE15.7))

	DO LX1=0,25,1
	XX=10.D0**(-4.D0+LX1/25.D0*3.D0)
        XARRAY(LX1+1) = DLOG(XX)
  	END DO
	DO LX1=1,50,1
	XX=0.1D0+LX1/50.D0*0.9D0
        XARRAY(LX1+26) = DLOG(XX)
	END DO
	DO LQ=0,25,1
	QQ=0.6D0*10**(5.D0*LQ/25.D0)
        XARRAY(LQ+77) = DLOG(QQ)
	END DO
 
	RETURN
	END