C-- This example program demonstrates how to obtain the
C-- MRW 2006 diffractive parton distributions (DPDFs)
C-- for any (xPom,z,Qsq).  These were obtained from
C-- from a "pQCD" fit to H1 LRG data: see hep-ph/0609273.
C-- The diffractive structure functions (DSFs),
C-- F2D3, FLD3, F2D3(charm) and FLD3(charm), are also
C-- provided for any (xPom,beta,Qsq).  The DSFs have three
C-- components: resolved Pomeron, direct Pomeron and
C-- secondary Reggeon.  Add all three together to get the
C-- total DSF.  CERNLIB is required.  Note that all grids
C-- are normalised to M_Y = m_p.  To normalise to
C-- M_Y < 1.6 GeV (H1 LRG data), multiply by a factor 1.23.
C-- Comments to watt(at)hep.ucl.ac.uk

      PROGRAM example

C--   Declare all variables:
      IMPLICIT NONE
      INTEGER iSet,I
      PARAMETER(iSet=1)         ! MRW 2006 (only one set at present)
      DOUBLE PRECISION xPom,z,Qsq,zSingletPom,zGluonPom,
     &     DPDFPom(-6:6),DPDFReg(-6:6),DSF(3),beta
      
C--   Read grid containing DPDFs from file and store in memory:
      CALL InitialiseDPDF(iSet)

C--   Read grid containing DSFs from file and store in memory:
      CALL InitialiseDSF(iSet,'F2D3     ')
      CALL InitialiseDSF(iSet,'FLD3     ')
      CALL InitialiseDSF(iSet,'F2D3charm')
      CALL InitialiseDSF(iSet,'FLD3charm')
C--   Note the five blank spaces after F2D3 and FLD3.

C--   N.B. Initialisation is slow, so don't initialise DPDFs or
C--        DSFs if you are not going to call them.

      WRITE(6,*)
      WRITE(6,*) "Calculate DPDFs ..."

C--   Specify xPom, z and Qsq:
      xPom = 0.01D0
      z = 0.5D0
      Qsq = 25.D0
      WRITE(6,*) "xPom,z,Qsq =",xPom,z,Qsq

C--   Get diffractive quark singlet distribution (Pomeron only):
      WRITE(6,*) "zSingletPom =",zSingletPom(xPom,z,Qsq)

C--   Get diffractive gluon distribution (Pomeron only):
      WRITE(6,*) "zGluonPom =",zGluonPom(xPom,z,Qsq)

C--   Or get DPDFs in PDFLIB-like format:
      CALL GetDPDF(xPom,z,Qsq,'Pom',DPDFPom) ! Pomeron
      CALL GetDPDF(xPom,z,Qsq,'Reg',DPDFReg) ! Reggeon      
      DO I = -6,6
C--   Notation: I = -6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6.
C--               = tbar,bbar,cbar,sbar,ubar,dbar,g,d,u,s,c,b,t.
         WRITE(6,*) "DPDFPom(",I,") = ",DPDFPom(I)
         WRITE(6,*) "DPDFReg(",I,") = ",DPDFReg(I)
      END DO

      WRITE(6,*)
      WRITE(6,*) "Calculate DSFs ..."

C--   Specify xPom, beta and Qsq:
      xPom = 0.01D0
      beta = 0.5D0
      Qsq = 25.D0
      WRITE(6,*) "xPom,beta,Qsq =",xPom,beta,Qsq
      WRITE(6,*)

C--   Get F2D3 (three separate contributions):
      CALL GetDSF(xPom,beta,Qsq,'F2D3     ',DSF)
      WRITE(6,*) "F2D3(Resolved Pomeron) = ",DSF(1)
      WRITE(6,*) "F2D3(Direct Pomeron) = ",DSF(2)
      WRITE(6,*) "F2D3(Secondary Reggeon) = ",DSF(3)
      WRITE(6,*) "F2D3(Total) = ",DSF(1)+DSF(2)+DSF(3)
      WRITE(6,*)

C--   Get FLD3 (three separate contributions):
      CALL GetDSF(xPom,beta,Qsq,'FLD3     ',DSF)
      WRITE(6,*) "FLD3(Resolved Pomeron) = ",DSF(1)
      WRITE(6,*) "FLD3(Direct Pomeron) = ",DSF(2)
      WRITE(6,*) "FLD3(Secondary Reggeon) = ",DSF(3)
      WRITE(6,*) "FLD3(Total) = ",DSF(1)+DSF(2)+DSF(3)
      WRITE(6,*)

C--   Get F2D3charm (three separate contributions):
      CALL GetDSF(xPom,beta,Qsq,'F2D3charm',DSF)
      WRITE(6,*) "F2D3charm(Resolved Pomeron) = ",DSF(1)
      WRITE(6,*) "F2D3charm(Direct Pomeron) = ",DSF(2)
      WRITE(6,*) "F2D3charm(Secondary Reggeon) = ",DSF(3)
      WRITE(6,*) "F2D3charm(Total) = ",DSF(1)+DSF(2)+DSF(3)
      WRITE(6,*)

C--   Get FLD3charm (three separate contributions):
      CALL GetDSF(xPom,beta,Qsq,'FLD3charm',DSF)
      WRITE(6,*) "FLD3charm(Resolved Pomeron) = ",DSF(1)
      WRITE(6,*) "FLD3charm(Direct Pomeron) = ",DSF(2)
      WRITE(6,*) "FLD3charm(Secondary Reggeon) = ",DSF(3)
      WRITE(6,*) "FLD3charm(Total) = ",DSF(1)+DSF(2)+DSF(3)
      WRITE(6,*)

      STOP
      END