!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!  Program example of the use of GJR08FFNS GRIDS  !!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! This program calculates and writes to standard output the pdfs and
!! their error estimates, alpha_s and its error estimate and the
!! predictions for F2_light, F2_charm and F2_bottom. This is done for 
!! the NLO(MSbar) fit (note that the GJR08FFNS package does not contain
!! error estimates or structure functions for the other fits presented
!! in the paper) at fixed values of x and Q^2.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      program example
       implicit none
       integer set
       double precision x,Q2
       double precision GJR08FFNSxuv,GJR08FFNSxdv,GJR08FFNSxgl,
     &                  GJR08FFNSxub,GJR08FFNSxdb,GJR08FFNSxsb,
     &                  GJR08FFNSalphas,
     &                  GJR08FFNSF2l,GJR08FFNSF2c,GJR08FFNSF2b,
     &                  xuv(-13:13),xdv(-13:13),xgl(-13:13),xub(-13:13),
     &                  xdb(-13:13),xsb(-13:13),alphas(-13:13),
     &                  exuv,exdv,exgl,exub,exdb,exsb,ealphas
      !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 GJR08FFNSinit(set)
          xuv(set) = GJR08FFNSxuv(x,Q2)
          xdv(set) = GJR08FFNSxdv(x,Q2)
          xgl(set) = GJR08FFNSxgl(x,Q2)
          xub(set) = GJR08FFNSxub(x,Q2)
          xdb(set) = GJR08FFNSxdb(x,Q2)
          xsb(set) = GJR08FFNSxsb(x,Q2)
          alphas(set) = GJR08FFNSalphas(Q2)
    1  end do
      !Error estimates.
       exuv = 0d0
       exdv = 0d0
       exgl = 0d0
       exub = 0d0
       exdb = 0d0
       exsb = 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))
          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)
       ealphas = 0.5d0*dsqrt(ealphas)
      !Reloading the NLO(MSbar) set for the calculation of the structure functions.
       call GJR08FFNSinit(0)
      !Output.
    3  format(2x,a,1e10.3,a,1e10.3)
    4  format(2x,a,1e10.3)
       write(*,*)
       write(*,*) "GJR08FFNS 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) "alpha_s =",alphas(0)," +-",ealphas
       write(*,4) "F2_light =",GJR08FFNSF2l(x,Q2)
       write(*,4) "F2_charm =",GJR08FFNSF2c(x,Q2)
       write(*,4) "F2_bottom =",GJR08FFNSF2b(x,Q2)
       stop
      end program example