!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!  Program example of the use of GJR08VFNS GRIDS  !!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! This program calculates and writes to standard output the VFNS pdfs,
!! their error estimates and alpha_s and its error estimate. This is
!! done for the NLO(MSbar) fit (note that the GJR08VFNS package does not
!! contain error estimates for the LO fit) at fixed values of x and Q^2.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      program example
       implicit none
       integer set
       double precision x,Q2
       double precision GJR08VFNSxuv,GJR08VFNSxdv,GJR08VFNSxgl,
     &                  GJR08VFNSxub,GJR08VFNSxdb,GJR08VFNSxsb,
     &                  GJR08VFNSxcb,GJR08VFNSxbb,GJR08VFNSalphas,
     &                  xuv(-13:13),xdv(-13:13),xgl(-13:13),xub(-13:13),
     &                  xdb(-13:13),xsb(-13:13),xcb(-13:13),xbb(-13:13),
     &                  alphas(-13:13),ealphas,
     &                  exuv,exdv,exgl,exub,exdb,exsb,excb,exbb
      !Fixed values of x and Q^2.
       x = 1d-1
       Q2 = 1d4
      !Values of the pdfs (times x) and alpha_s for the different sets.
       do set=-13,13
          call GJR08VFNSinit(set)
          xuv(set) = GJR08VFNSxuv(x,Q2)
          xdv(set) = GJR08VFNSxdv(x,Q2)
          xgl(set) = GJR08VFNSxgl(x,Q2)
          xub(set) = GJR08VFNSxub(x,Q2)
          xdb(set) = GJR08VFNSxdb(x,Q2)
          xsb(set) = GJR08VFNSxsb(x,Q2)
          xcb(set) = GJR08VFNSxcb(x,Q2)
          xbb(set) = GJR08VFNSxbb(x,Q2)
          alphas(set) = GJR08VFNSalphas(Q2)
    1  end do
      !Error estimates.
       exuv = 0d0
       exdv = 0d0
       exgl = 0d0
       exub = 0d0
       exdb = 0d0
       exsb = 0d0
       excb = 0d0
       exbb = 0d0
       ealphas = 0d0
       do set = 1,13
          exuv = exuv + (xuv(set) - xuv(-set))*(xuv(set) - xuv(-set))
          exdv = exdv + (xdv(set) - xdv(-set))*(xdv(set) - xdv(-set))
          exgl = exgl + (xgl(set) - xgl(-set))*(xgl(set) - xgl(-set))
          exub = exub + (xub(set) - xub(-set))*(xub(set) - xub(-set))
          exdb = exdb + (xdb(set) - xdb(-set))*(xdb(set) - xdb(-set))
          exsb = exsb + (xsb(set) - xsb(-set))*(xsb(set) - xsb(-set))
          excb = excb + (xcb(set) - xcb(-set))*(xcb(set) - xcb(-set))
          exbb = exbb + (xbb(set) - xbb(-set))*(xbb(set) - xbb(-set))
          ealphas = ealphas + (alphas(set) - alphas(-set))**2d0
    2  end do
       exuv = 0.5d0*dsqrt(exuv)
       exdv = 0.5d0*dsqrt(exdv)
       exgl = 0.5d0*dsqrt(exgl)
       exub = 0.5d0*dsqrt(exub)
       exdb = 0.5d0*dsqrt(exdb)
       exsb = 0.5d0*dsqrt(exsb)
       excb = 0.5d0*dsqrt(excb)
       exbb = 0.5d0*dsqrt(exbb)
       ealphas = 0.5d0*dsqrt(ealphas)
      !Output.
    3  format(2x,a,1e10.3,a,1e10.3)
       write(*,*)
       write(*,*) "GJR08VFNS GRIDS:"
       write(*,3) "x =",x,", Q^2 =",Q2
       write(*,3) "x*uv =",xuv(0)," +-",exuv
       write(*,3) "x*dv =",xdv(0)," +-",exdv
       write(*,3) "x*gl =",xgl(0)," +-",exgl
       write(*,3) "x*ub =",xub(0)," +-",exub
       write(*,3) "x*db =",xdb(0)," +-",exdb
       write(*,3) "x*sb =",xsb(0)," +-",exsb
       write(*,3) "x*cb =",xcb(0)," +-",excb
       write(*,3) "x*bb =",xbb(0)," +-",exbb
       write(*,3) "alpha_s =",alphas(0)," +-",ealphas
       stop
      end program example