**********************************************************
      program example
**********************************************************
*
*     calculation of BB QCD-NLO(LO) g1p/n/d & errors
*
*     example program
*
**********************************************************
*
      implicit none
*
      integer nq2,nxbj
      parameter (nq2=1 , nxbj=10)
*
      integer iini,iset
      integer iq2,ixbj
*
      double precision q2ref(nq2)
      double precision xbj(nxbj)
      double precision x,q2
      double precision uv,dv,gl,sea,qb,g1p,g1n,g1d
      double precision duv,ddv,dgl,dsea,dqb,dg1p,dg1n,dg1d
*    
      common /intini/ iini
*
*     set data:
*
      data q2ref  / 4.0d0 /
*
      data xbj / 1.0D-1, 2.0D-1, 3.0D-1, 4.0D-1, 5.0D-1, 
     +           6.0D-1, 7.0D-1, 8.0D-1, 9.0D-1, 9.5D-1 /
*
*     start reading data
*
      iini = 0                  ! reading data first time
*
*     choose order of parton sets and fit mode
*
c      iset = 1                  ! LO parton sets , Scenario 1
c      iset = 2                  ! LO parton sets , Scenario 2
      iset = 3                  ! NLO parton sets , Scenario 1
c      iset = 4                  ! NLO parton sets , Scenario 2
*
*     write file header
*
      write(6,2000)
*
*     g1 and pdf vs x at Q2ref
*
      do iq2 = 1,nq2
*
         q2 = q2ref(iq2)
*
         do ixbj = 1,nxbj
*
            x = xbj(ixbj)
*
            call ppdf(iset, x, q2, uv, duv, dv, ddv, gl, dgl,
     1           qb, dqb, g1p, dg1p, g1n, dg1n)

            g1d = 0.5d0 * (g1p+g1n) * (1.0d0-1.5d0*0.05d0)
            dg1d = 0.5d0 * sqrt(dg1p*dg1p+dg1n*dg1n) * 
     1           (1.0d0-1.5d0*0.05d0)
*
*     output to file
*
            write(6,2001) x,q2,g1p,dg1p,g1n,dg1n,g1d,dg1d
*
         enddo
*
      enddo
*
 2000 format('* x,q2,g1p,dg1p,g1n,dg1n,g1d,dg1d')
 2001 format(8(1x,f9.6))
*
      end
*
      include 'ppdf.f'
*
*=======================================================================