!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Program example of the use of JR09VFNNLO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! This program calculates and writes to standard output the pdfs and
!! alpha_s together with their error estimates. This is done for the
!! dynamical VFNS NLO(MSbar) distributions at fixed x and Q^2.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      program example
       implicit none
       integer i,set
       double precision x,Q2
       double precision JR09VFNNLOxuv,JR09VFNNLOxdv,JR09VFNNLOxgl,
     &                  JR09VFNNLOxub,JR09VFNNLOxdb,JR09VFNNLOxsb,
     &                  JR09VFNNLOxcb,JR09VFNNLOxbb,
     &                  JR09VFNNLOalphas,
     &                  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),
     &                  exuv,exdv,exgl,exub,exdb,exsb,excb,exbb,ealphas
      !Ininialization.
      call JR09VFNNLOinit
      !Fixed values of x and Q^2.
       x = 1.308d-4
       Q2 = 8314d0
c       Q2 = 1d2
      !Values of the pdfs (times x) and alpha_s for the different sets.
       do set=-13,13
          xuv(set) = JR09VFNNLOxuv(x,Q2,set)
          xdv(set) = JR09VFNNLOxdv(x,Q2,set)
          xgl(set) = JR09VFNNLOxgl(x,Q2,set)
          xub(set) = JR09VFNNLOxub(x,Q2,set)
          xdb(set) = JR09VFNNLOxdb(x,Q2,set)
          xsb(set) = JR09VFNNLOxsb(x,Q2,set)
          xcb(set) = JR09VFNNLOxcb(x,Q2,set)
          xbb(set) = JR09VFNNLOxbb(x,Q2,set)
          alphas(set) = JR09VFNNLOalphas(Q2,set)
    1  end do
      !Error estimates.
       exuv = 0d0
       exdv = 0d0
       exgl = 0d0
       exub = 0d0
       exdb = 0d0
       exsb = 0d0
       excb = 0d0
       exbb = 0d0
       ealphas = 0d0
       do i = 1,13
          exuv = exuv + (xuv(i) - xuv(-i))**2d0
          exdv = exdv + (xdv(i) - xdv(-i))**2d0
          exgl = exgl + (xgl(i) - xgl(-i))**2d0
          exub = exub + (xub(i) - xub(-i))**2d0
          exdb = exdb + (xdb(i) - xdb(-i))**2d0
          exsb = exsb + (xsb(i) - xsb(-i))**2d0
          excb = excb + (xcb(i) - xcb(-i))**2d0
          exbb = exbb + (xbb(i) - xbb(-i))**2d0
          ealphas = ealphas + (alphas(i) - alphas(-i))**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,1E11.4,a,1E11.4)
    4  format(2x,a,1E11.4)
       write(*,*)
       write(*,*) "JR09VFNNLO:"
       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