C *********************************************************************
      PROGRAM SAMPLE_1
C ---------------------------------------------------------------------
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION XPDFA(-4:4),gradpdf(19),E(4),AA(9)
      COMMON/ERRM/EM(9,9)
C User should supply XMAX, XMIN, Q2, ISET, and NX in the following.
      DATA XMIN,XMAX/ 1.D-3,1.D0/  ! x-range for NPDF output
      DATA Q2/5.D0/                ! Q^2 value [GeV^2]
C ISET (Type of the nucleus, see npdf04_2.f), 
      DATA ISET/17/  ! 17:Au
c      DATA ISET/2/  !  2: D
c      DATA ISET/1/  !  1: p
      DATA NX/200/                 ! Number of x-steps for writing NPDFs
    
        JSET=3
        DLMAX=DLOG10(XMAX)
        DLMIN=DLOG10(XMIN)
        XLSTEP=(DLMAX-DLMIN)/DFLOAT(NX)

        DO 10 IX=1,NX+1
          DLOGX=DFLOAT(IX-1)*XLSTEP+DLMIN
          X=10.D0**DLOGX
        
C CAll NPDF04 library
         Call NPDF04(ISET,JSET,Q2,X,XPDFA,GRADPDF)
C XPDFA(I) (0: x*gluon, 1: x*u=x*(uv+usea), 2: x*d=x*(dv+dsea),
C           3: x*strange, 4: x*charm)    
      write(21,1010) X,XPDFA(0)          ! xg
     +                ,XPDFA(1)          ! xu (=xuv+xubar)
     +                ,XPDFA(2)          ! xd (=xdv+xdbar)
     +                ,XPDFA(3)          ! xqbar (=xs=xubar=xdbar)
     +                ,XPDFA(4)          ! xc
 1010    FORMAT (1X,6(1PE16.6))

c---------------- NPDF uncertainty by the Hessian method---------------
C GRADPDF(I) (1..9:gluon, 10..12:uv, 13..15:dv, 16..19:antiquark)
        DO N=1,4
          IF(N.EQ.1)THEN      ! xg^A
            NP=9;NS=0;NT=0
          ELSE IF(N.EQ.2)THEN ! xu^A
            NP=7;NS=2;NT=9
          ELSE IF(N.EQ.3)THEN ! xd^A
            NP=7;NS=2;NT=12
          ELSE IF(N.EQ.4)THEN ! xqb^A
            NP=4;NS=5;NT=15
          END IF

          DO I=1,9
            AA(I)=0.D0
          ENDDO  
          EE=0.D0

          DO I=1,NP     ! NPDF uncertainty by the Hessian method 
            DO J=1,NP
              IF((N.EQ.2).and.(J.GT.3)) THEN
                AA(I)=AA(I)+EM(NS+I,NS+J)*GRADPDF(12+J)
              ELSE
                AA(I)=AA(I)+EM(NS+I,NS+J)*GRADPDF(NT+J)
              ENDIF
            ENDDO
            IF((N.EQ.2).and.(I.GT.3)) THEN
              EE=EE+GRADPDF(12+I)*AA(I)
            ELSE
              EE=EE+GRADPDF(NT+I)*AA(I)
            ENDIF
          ENDDO
          E(N)=dsqrt(EE)

        END DO

c Note there is no uncertainty for the charm distribution because it
c is created through the Q2 evolution. 
c                     eg^A,  eu^A, ed^A, eqb^A
      write(22,1011) X,E(1)  ! uncertainty for xg
     +                ,E(2)  !                 xu
     +                ,E(3)  !                 xd
     +                ,E(4)  !                 xqbar (=xs=xubar=xdbar) 
 1011   FORMAT (1X,5(1PE16.6))
 10   CONTINUE 

      STOP
      END
C *********************************************************************
C THE END OF THE PROGRAM.
C *********************************************************************