**********************************************************
      program example
**********************************************************
*
*    calculation of LSS'06 NLO PDFs and g1(p,n,d)
*
*     example program
*
**********************************************************
*
c      implicit none
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
*
      integer nq2,nxbj
      parameter (nq2=29 , nxbj=48)
*
      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,g1pLT,g1p,g1nLT,g1n,g1d
      double precision duv,ddv,dgl,dsea,dqb,dg1pLT,dg1p,dg1nLT,dg1n,dg1d
*    
      common /intini/ iini

*
*     set data:
*
c        data q2ref  / 3.0d0 /


c       DATA xbj /
c     1           
c     3           1.D-2, 1.5D-2, 2.2D-2, 3.2D-2, 5.0D-2, 7.5D-2,
c     4           0.1, 0.125, 0.15, 0.175, 0.2, 0.225, 0.25, 0.275,
c     5           0.3, 0.325, 0.35, 0.375, 0.4, 0.45,  0.5, 0.55,
c     6           0.6, 0.65,0.7,0.75,0.8 /

       DATA q2ref / 1.0D0, 1.25d0, 1.5D0, 2.d0, 2.5D0,
     1           4.0D0, 6.4D0, 1.0D1, 1.5D1, 2.5D1, 4.0D1, 6.4D1,
     2           1.0D2, 1.8D2, 3.2D2, 5.8D2, 1.0D3, 1.8D3,
     3           3.2D3, 5.8D3, 1.0D4, 1.8D4, 3.2D4, 5.8D4,
     4           1.0D5, 1.8D5, 3.2D5, 5.8D5, 1.0D6  /
       DATA xbj /
     1           1.D-5, 1.5D-5, 2.2D-5, 3.2D-5, 4.8D-5, 7.D-5,
     1           1.D-4, 1.5D-4, 2.2D-4, 3.2D-4, 4.8D-4, 7.D-4,
     2           1.D-3, 1.5D-3, 2.2D-3, 3.2D-3, 4.8D-3, 7.D-3,
     3           1.D-2, 1.5D-2, 2.2D-2, 3.2D-2, 5.0D-2, 7.5D-2,
     4           0.1, 0.125, 0.15, 0.175, 0.2, 0.225, 0.25, 0.275,
     5           0.3, 0.325, 0.35, 0.375, 0.4, 0.45,  0.5, 0.55,
     6           0.6, 0.65,  0.7,  0.75,  0.8, 0.85,  0.9, 1.0 /

      OPEN(6,FILE='Print_grids_iset123a.dat'
     1 ,STATUS='NEW',FORM='FORMATTED')

*

                                       do iiset=1,3
*     start reading data
*
      iini = 0                  ! reading data first time
*
*     choose order of parton sets and fit mode
*   INPUT:   ISET = number of the parton set                         *
*              (TO BE DEFINED BY THE USER ):                         *
*            ISET = 1   NEXT-TO-LEADING ORDER (xDelta G > 0)         *
*                      (DATA FILE 'NLO_MS_delGpos.grid' UNIT=11)     *
*            ISET = 2   NEXT-TO-LEADING ORDER (xDelta G < 0)         *
*                      (DATA FILE 'NLO_MS_delGneg.grid' UNIT=22)     *
*            ISET = 3   NEXT-TO-LEADING ORDER (changing sign         *
*                       xDelta G)                                    *
*                      (DATA FILE 'NLO_MS_chsign_delG.grid' UNIT=33) *
*     write file header
*


*
*     g1 and pdf vs x at Q2ref
*


*         MS-bar
                              	iset=iiset
  
      PRINT '(''iset = '' ,i3)'
     1 ,iset	  
         write(6,2000)	 
      do iq2 = 1,nq2-1
*
         q2 = q2ref(iq2)

*
         do ixbj = 1,nxbj

*
            x = xbj(ixbj)

c     To obtain the values of PPDs and g1 outside the 
c     resonanse region use the cut W**2>4 GeV**2  
    
      call LSS2006 (ISET,X,Q2,UUB,DDB,SSB,GL,UV,DV,UB,DB,ST,
     1g1pLT,g1p,g1nLT,g1n)

          g1d = 0.5d0 * (g1p+g1n) * (1.0d0-1.5d0*0.05d0)

        write(6,2001) q2,x,UUB,DDB,SSB,GL,UV,DV,UB,DB,ST
	1,g1pLT,g1p,g1nLT,g1n
*
      enddo

      enddo
 
         enddo


 2000 format('* Q2 x UUB DDB SSB GL UV DV UB DB ST g1pLT g1p g1nLT g1n')
 2001   FORMAT (2e9.3,13(1pe12.4))
*
      end

*


      include 'LSS2006pdf_g1.f'        

*=======================================================================