!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Program example of the use of JR09FFNNLO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! This program calculates and writes to standard output the pdfs and
!! alpha_s together with their error estimates. This is done for the
!! dynamical FFNS NNLO(MSbar) distributions at fixed x and Q^2.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      program example
       implicit none
       integer i,set
       double precision x,Q2
       double precision JR09FFNNLOxuv,JR09FFNNLOxdv,JR09FFNNLOxgl,
     &                  JR09FFNNLOxub,JR09FFNNLOxdb,JR09FFNNLOxsb,
     &                  JR09FFNNLOalphas,
     &                  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
      !Ininialization.
      call JR09FFNNLOinit
      !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
          xuv(set) = JR09FFNNLOxuv(x,Q2,set)
          xdv(set) = JR09FFNNLOxdv(x,Q2,set)
          xgl(set) = JR09FFNNLOxgl(x,Q2,set)
          xub(set) = JR09FFNNLOxub(x,Q2,set)
          xdb(set) = JR09FFNNLOxdb(x,Q2,set)
          xsb(set) = JR09FFNNLOxsb(x,Q2,set)
          alphas(set) = JR09FFNNLOalphas(Q2,set)
    1  end do
      !Error estimates.
       exuv = 0d0
       exdv = 0d0
       exgl = 0d0
       exub = 0d0
       exdb = 0d0
       exsb = 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
          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)
       ealphas = 0.5d0*dsqrt(ealphas)
      !Output.
    3  format(2x,a,1E11.4,a,1E11.4)
    4  format(2x,a,1E11.4)
       write(*,*)
       write(*,*) "JR09FFNNLO:"
       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
       stop
      end program example