C-- Public access to MRW 2006 diffractive parton distributions (DPDFs) C-- and diffractive structure functions (DSFs). See example.f for C-- usage. CERNLIB is required for 3-D linear interpolation (FINT) C-- and to access the GRV-P-HO pion PDFs (from PDFLIB). C-- Comments to watt(at)hep.ucl.ac.uk C-- ============================================================== SUBROUTINE InitialiseDPDF(iSet) C-- Purpose : Read DPDF grids from file into memory. C-- Input : iSet,file containing grids. C-- Output : Grids stored in COMMON blocks. C-- ============================================================== C-- Declare all variables: IMPLICIT NONE LOGICAL DPDFinitialised INTEGER io,nXPomGrid,nZGrid,nQsqGrid,iXPom,iZ,iQsq,nCoord, & nXPomGridFile,nZGridFile,nQsqGridFile,iSet PARAMETER(nXPomGrid=80,nZGrid=80,nQsqGrid=60,nCoord=3) INTEGER nXPomZQsqGrid(nCoord) CHARACTER dummy CHARACTER*30 filename CHARACTER*20 PARM(20) DOUBLE PRECISION xPomMin,xPomMax,zMin,zMax,QsqMin,QsqMax, & xPom,z,Qsq,cReg,VALU(20), & xPomGrid(nXPomGrid),zGrid(nZGrid+1),QsqGrid(nQsqGrid) REAL xPomZQsqGrid(nXPomGrid+nZGrid+1+nQsqGrid), & zSingletPomGrid(nXPomGrid,nZGrid+1,nQsqGrid), & zGluonPomGrid(nXPomGrid,nZGrid+1,nQsqGrid) COMMON / commonDPDFGridPts / nXPomZQsqGrid,xPomZQsqGrid COMMON / commonSingletPom / zSingletPomGrid COMMON / commonGluonPom / zGluonPomGrid COMMON / commonCReg / cReg COMMON / commonDPDFinitialised / DPDFinitialised DATA DPDFinitialised / .FALSE. / DATA nXPomZQsqGrid / nXPomGrid, nZGrid+1, nQsqGrid / C-- Check if already initialised: IF (DPDFinitialised) THEN WRITE(6,*) "Warning: already initialised DPDFs." RETURN END IF C-- Open file containing grid: IF (iSet.EQ.1) THEN ! MRW 2006 fit to H1 LRG DDIS data filename = 'mrw2006dpdfgrid.dat' ELSE WRITE(6,*) "Error: iSet = ",iSet STOP END IF WRITE(6,*) "Reading DPDF grid from ",filename OPEN (UNIT=10,FILE=filename,IOSTAT=io,STATUS='OLD') IF (io.NE.0) THEN WRITE(6,*) "Error: FILE = ",filename,", IOSTAT = ",io STOP END IF C-- Read secondary Reggeon normalisation from file: READ(10,*); READ(10,*); READ(10,*); READ(10,*) READ(10,*); READ(10,*) dummy,cReg; READ(10,*) WRITE(6,'(" Read secondary Reggeon normalisation: cReg =", & D13.5)') cReg C-- Use GRV-P-HO pion PDFs (from PDFLIB) for secondary Reggeon: PARM(1) = 'NPTYPE' VALU(1) = 2 PARM(2) = 'NGROUP' VALU(2) = 5 PARM(3) = 'NSET' VALU(3) = 1 WRITE(6,*) " Initialising GRV-P-HO pion PDFS from PDFLIB ..." CALL PDFSET(PARM,VALU) ! Set PDFLIB parameters WRITE(6,*) " ... pion PDFs initialised!" C-- Read (xPom,z,Qsq) grid dimensions from file: READ(10,*); READ(10,*) dummy,nXPomGridFile,xPomMin,xPomMax READ(10,*); READ(10,*) dummy,nZGridFile,zMin,zMax READ(10,*); READ(10,*) dummy,nQsqGridFile,QsqMin,QsqMax C-- Check grid dimensions are correct: IF ((nXPomGridFile.NE.nXPomGrid).OR.(nZGridFile.NE.nZGrid) & .OR.(nQsqGridFile.NE.nQsqGrid)) THEN WRITE(*,*) "Error: Grid from ",filename,"has wrong dimensions" WRITE(*,*) nXPomGridFile, nZGridFile, nQsqGridFile WRITE(*,*) nXPomGrid, nZGrid, nQsqGrid STOP END IF C WRITE(6,*) " nXPomGrid,xPomMin,xPomMax =", C & nXPomGrid,xPomMin,xPomMax C WRITE(6,*) " nZGrid,zMin,zMax =", C & nZGrid,zMin,zMax C WRITE(6,*) " nQsqGrid,QsqMin,QsqMax =",nQsqGrid,QsqMin,QsqMax C-- Add extra grid point at z = 1: zGrid(nZGrid+1) = 1.D0 C-- Loop over xPom, Qsq values: DO iXPom = 1, nXPomGrid DO iQsq = 1, nQsqGrid C-- Store xPom,Qsq values of grid points: READ(10,*); READ(10,*); READ(10,*) dummy,xPom,Qsq xPomGrid(iXPom) = xPom IF (iXPom.EQ.1) QsqGrid(iQsq) = Qsq READ(10,*) C-- Loop over z values: DO iZ = 1, nZGrid C-- Read DPDFs of grid points: IF ((iXPom.EQ.1).AND.(iQsq.EQ.1)) THEN READ(10,*) z,zSingletPomGrid(iXPom,iZ,iQsq), & zGluonPomGrid(iXPom,iZ,iQsq) ELSE READ(10,*) zSingletPomGrid(iXPom,iZ,iQsq), & zGluonPomGrid(iXPom,iZ,iQsq) END IF C-- Store z values of grid points: IF ((iXPom.EQ.1).AND.(iQsq.EQ.1)) zGrid(iZ) = z END DO ! end loop over iZ C-- DPDFs are zero at z = 1: zSingletPomGrid(iXPom,nZGrid+1,iQsq) = 0.D0 zGluonPomGrid(iXPom,nZGrid+1,iQsq) = 0.D0 END DO ! end loop over iXPom END DO ! end loop over iQsq C-- Close file: CLOSE(UNIT=10) C WRITE(6,*) "xPomGrid =",xPomGrid C WRITE(6,*) "zGrid =",zGrid C WRITE(6,*) "QsqGrid =",QsqGrid C-- Store all grid points sequentially in a 1-D REAL array: DO iXPom = 1, nXPomGrid xPomZQsqGrid(iXPom) = REAL(xPomGrid(iXPom)) END DO DO iZ = 1, nZGrid+1 xPomZQsqGrid(nXPomGrid+iZ) = REAL(zGrid(iZ)) END DO DO iQsq = 1, nQsqGrid xPomZQsqGrid(nXPomGrid+nZGrid+1+iQsq) = & REAL(QsqGrid(iQsq)) END DO C WRITE(6,*) "nXPomZQsqGrid =",nXPomZQsqGrid C WRITE(6,*) "xPomZQsqGrid =",xPomZQsqGrid DPDFinitialised = .TRUE. WRITE(6,*) "... DPDFs initialised!" WRITE(6,*) RETURN END C-- ============================================================== SUBROUTINE InitialiseDSF(iSet,whichDSF) C-- Purpose : Read DSF grids from file into memory. C-- Input : iSet,whichDSF,file containing grids. C-- Output : Grids stored in COMMON blocks. C-- ============================================================== C-- Declare all variables: IMPLICIT NONE LOGICAL F2D3initialised,FLD3initialised, & F2D3charminitialised,FLD3charminitialised CHARACTER*9 whichDSF INTEGER io,nXPomGrid,nBetaGrid,nQsqGrid,iXPom,iBeta,iQsq,nCoord, & nXPomGridFile,nBetaGridFile,nQsqGridFile,iSet PARAMETER(nXPomGrid=80,nBetaGrid=80,nQsqGrid=60,nCoord=3) INTEGER nXPomBetaQsqGrid(nCoord) CHARACTER dummy CHARACTER*30 filename DOUBLE PRECISION xPomMin,xPomMax,betaMin,betaMax,QsqMin,QsqMax, & xPom,beta,Qsq, & xPomGrid(nXPomGrid),betaGrid(nBetaGrid+1),QsqGrid(nQsqGrid) REAL xPomBetaQsqGrid(nXPomGrid+nBetaGrid+1+nQsqGrid), & F2D3ResPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & F2D3DirPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & F2D3RegGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & FLD3ResPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & FLD3DirPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & FLD3RegGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & F2D3charmResPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & F2D3charmDirPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & F2D3charmRegGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & FLD3charmResPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & FLD3charmDirPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & FLD3charmRegGrid(nXPomGrid,nBetaGrid+1,nQsqGrid) COMMON / commonDSFGridPts / nXPomBetaQsqGrid,xPomBetaQsqGrid COMMON / commonF2D3 / F2D3ResPomGrid,F2D3DirPomGrid,F2D3RegGrid COMMON / commonFLD3 / FLD3ResPomGrid,FLD3DirPomGrid,FLD3RegGrid COMMON / commonF2D3charm / F2D3charmResPomGrid, & F2D3charmDirPomGrid,F2D3charmRegGrid COMMON / commonFLD3charm / FLD3charmResPomGrid, & FLD3charmDirPomGrid,FLD3charmRegGrid COMMON / commonDSFinitialised / F2D3initialised,FLD3initialised, & F2D3charminitialised,FLD3charminitialised DATA F2D3initialised,FLD3initialised, & F2D3charminitialised,FLD3charminitialised / 4*.FALSE. / DATA nXPomBetaQsqGrid / nXPomGrid, nBetaGrid+1, nQsqGrid / C-- Check if DSFs are already initialised: IF ((whichDSF.EQ.'F2D3').AND.F2D3initialised) THEN WRITE(6,*) "Warning: already initialised F2D3." RETURN ELSE IF ((whichDSF.EQ.'FLD3').AND.FLD3initialised) THEN WRITE(6,*) "Warning: already initialised FLD3." RETURN ELSE IF ((whichDSF.EQ.'F2D3charm').AND.F2D3charminitialised) THEN WRITE(6,*) "Warning: already initialised F2D3charm." RETURN ELSE IF ((whichDSF.EQ.'FLD3charm').AND.FLD3charminitialised) THEN WRITE(6,*) "Warning: already initialised FLD3charm." RETURN END IF C-- Open file containing grid: IF (iSet.EQ.1) THEN ! MRW 2006 fit to H1 LRG DDIS data IF (whichDSF.EQ.'F2D3') THEN F2D3initialised = .TRUE. filename = "mrw2006f2d3grid.dat" ELSE IF (whichDSF.EQ.'FLD3') THEN FLD3initialised = .TRUE. filename = "mrw2006fld3grid.dat" ELSE IF (whichDSF.EQ.'F2D3charm') THEN F2D3charminitialised = .TRUE. filename = "mrw2006f2d3charmgrid.dat" ELSE IF (whichDSF.EQ.'FLD3charm') THEN FLD3charminitialised = .TRUE. filename = "mrw2006fld3charmgrid.dat" ELSE WRITE(6,*) "Error in InitialiseDSF: whichDSF = ",whichDSF STOP END IF ELSE WRITE(6,*) "Error: iSet = ",iSet STOP END IF WRITE(6,*) "Reading DSF grid from ",filename OPEN (UNIT=10,FILE=filename,IOSTAT=io,STATUS='OLD') IF (io.NE.0) THEN WRITE(6,*) "Error: FILE = ",filename,", IOSTAT = ",io STOP END IF READ(10,*); READ(10,*); READ(10,*); READ(10,*) C-- Read (xPom,z,Qsq) grid dimensions from file: READ(10,*); READ(10,*) dummy,nXPomGridFile,xPomMin,xPomMax READ(10,*); READ(10,*) dummy,nBetaGridFile,betaMin,betaMax READ(10,*); READ(10,*) dummy,nQsqGridFile,QsqMin,QsqMax C-- Check grid dimensions are correct: IF ((nXPomGridFile.NE.nXPomGrid).OR.(nBetaGridFile.NE.nBetaGrid) & .OR.(nQsqGridFile.NE.nQsqGrid)) THEN WRITE(*,*) "Error: Grid from ",filename,"has wrong dimensions" WRITE(*,*) nXPomGridFile, nBetaGridFile, nQsqGridFile WRITE(*,*) nXPomGrid, nBetaGrid, nQsqGrid STOP END IF C WRITE(6,*) " nXPomGrid,xPomMin,xPomMax =", C & nXPomGrid,xPomMin,xPomMax C WRITE(6,*) " nBetaGrid,betaMin,betaMax =", C & nBetaGrid,betaMin,betaMax C WRITE(6,*) " nQsqGrid,QsqMin,QsqMax =",nQsqGrid,QsqMin,QsqMax C-- Loop over xPom, Qsq values: DO iXPom = 1, nXPomGrid DO iQsq = 1, nQsqGrid C-- Store xPom,Qsq values of grid points: READ(10,*); READ(10,*); READ(10,*) dummy,xPom,Qsq xPomGrid(iXPom) = xPom IF (iXPom.EQ.1) QsqGrid(iQsq) = Qsq READ(10,*) C-- Loop over beta values: DO iBeta = 1, nBetaGrid+1 C-- Read DSFs of grid points: IF (whichDSF.EQ.'F2D3') THEN IF ((iXPom.EQ.1).AND.(iQsq.EQ.1)) THEN READ(10,*) beta,F2D3ResPomGrid(iXPom,iBeta,iQsq), & F2D3DirPomGrid(iXPom,iBeta,iQsq), & F2D3RegGrid(iXPom,iBeta,iQsq) ELSE READ(10,*) F2D3ResPomGrid(iXPom,iBeta,iQsq), & F2D3DirPomGrid(iXPom,iBeta,iQsq), & F2D3RegGrid(iXPom,iBeta,iQsq) END IF ELSE IF (whichDSF.EQ.'FLD3') THEN IF ((iXPom.EQ.1).AND.(iQsq.EQ.1)) THEN READ(10,*) beta,FLD3ResPomGrid(iXPom,iBeta,iQsq), & FLD3DirPomGrid(iXPom,iBeta,iQsq), & FLD3RegGrid(iXPom,iBeta,iQsq) ELSE READ(10,*) FLD3ResPomGrid(iXPom,iBeta,iQsq), & FLD3DirPomGrid(iXPom,iBeta,iQsq), & FLD3RegGrid(iXPom,iBeta,iQsq) END IF ELSE IF (whichDSF.EQ.'F2D3charm') THEN IF ((iXPom.EQ.1).AND.(iQsq.EQ.1)) THEN READ(10,*) beta, & F2D3charmResPomGrid(iXPom,iBeta,iQsq), & F2D3charmDirPomGrid(iXPom,iBeta,iQsq), & F2D3charmRegGrid(iXPom,iBeta,iQsq) ELSE READ(10,*) F2D3charmResPomGrid(iXPom,iBeta,iQsq), & F2D3charmDirPomGrid(iXPom,iBeta,iQsq), & F2D3charmRegGrid(iXPom,iBeta,iQsq) END IF ELSE IF (whichDSF.EQ.'FLD3charm') THEN IF ((iXPom.EQ.1).AND.(iQsq.EQ.1)) THEN READ(10,*) beta, & FLD3charmResPomGrid(iXPom,iBeta,iQsq), & FLD3charmDirPomGrid(iXPom,iBeta,iQsq), & FLD3charmRegGrid(iXPom,iBeta,iQsq) ELSE READ(10,*) FLD3charmResPomGrid(iXPom,iBeta,iQsq), & FLD3charmDirPomGrid(iXPom,iBeta,iQsq), & FLD3charmRegGrid(iXPom,iBeta,iQsq) END IF ELSE WRITE(6,*) "Error in InitialiseDSF: whichDSF = ", & whichDSF STOP END IF C-- Store beta values of grid points: IF ((iXPom.EQ.1).AND.(iQsq.EQ.1)) betaGrid(iBeta) = beta END DO ! end loop over iBeta END DO ! end loop over iXPom END DO ! end loop over iQsq C-- Close file: CLOSE(UNIT=10) C WRITE(6,*) "xPomGrid =",xPomGrid C WRITE(6,*) "betaGrid =",betaGrid C WRITE(6,*) "QsqGrid =",QsqGrid C-- Store all grid points sequentially in a 1-D REAL array: DO iXPom = 1, nXPomGrid xPomBetaQsqGrid(iXPom) = REAL(xPomGrid(iXPom)) END DO DO iBeta = 1, nBetaGrid+1 xPomBetaQsqGrid(nXPomGrid+iBeta) = REAL(betaGrid(iBeta)) END DO DO iQsq = 1, nQsqGrid xPomBetaQsqGrid(nXPomGrid+nBetaGrid+1+iQsq) = & REAL(QsqGrid(iQsq)) END DO C WRITE(6,*) "nXPomBetaQsqGrid =",nXPomBetaQsqGrid C WRITE(6,*) "xPomBetaQsqGrid =",xPomBetaQsqGrid WRITE(6,*) "... DSFs initialised!" WRITE(6,*) RETURN END C-- ============================================================== DOUBLE PRECISION FUNCTION zSingletPom(xPom,z,Qsq) C-- Purpose : Interpolate diffractive quark singlet distribution C-- (Pomeron contribution). C-- Input : xPom,z,Qsq C-- Output : z*Sigma^D(xPom,z,Qsq) [Pomeron only] C-- ============================================================== IMPLICIT NONE LOGICAL DPDFinitialised INTEGER nXPomGrid,nZGrid,nQsqGrid,nCoord PARAMETER(nXPomGrid=80,nZGrid=80,nQsqGrid=60,nCoord=3) INTEGER nXPomZQsqGrid(nCoord) DOUBLE PRECISION xPom,z,Qsq REAL FINT,xPomZQsq(nCoord), & xPomZQsqGrid(nXPomGrid+nZGrid+1+nQsqGrid), & zSingletPomGrid(nXPomGrid,nZGrid+1,nQsqGrid) COMMON / commonDPDFGridPts / nXPomZQsqGrid,xPomZQsqGrid COMMON / commonSingletPom / zSingletPomGrid COMMON / commonDPDFInitialised / DPDFinitialised IF (.NOT.DPDFinitialised) THEN WRITE(6,*) "Error in zSingletPom: DPDFs not initialised." STOP END IF xPomZQsq(1) = REAL(xPom) xPomZQsq(2) = REAL(z) xPomZQsq(3) = REAL(Qsq) zSingletPom = FINT(nCoord,xPomZQsq,nXPomZQsqGrid, & xPomZQsqGrid,zSingletPomGrid) RETURN END C-- ============================================================== DOUBLE PRECISION FUNCTION zGluonPom(xPom,z,Qsq) C-- Purpose : Interpolate diffractive gluon distribution C-- (Pomeron contribution). C-- Input : xPom,z,Qsq C-- Output : z*g^D(xPom,z,Qsq) [Pomeron only] C-- ============================================================== IMPLICIT NONE LOGICAL DPDFinitialised INTEGER nXPomGrid,nZGrid,nQsqGrid,nCoord PARAMETER(nXPomGrid=80,nZGrid=80,nQsqGrid=60,nCoord=3) INTEGER nXPomZQsqGrid(nCoord) DOUBLE PRECISION xPom,z,Qsq REAL FINT,xPomZQsq(nCoord), & xPomZQsqGrid(nXPomGrid+nZGrid+1+nQsqGrid), & zGluonPomGrid(nXPomGrid,nZGrid+1,nQsqGrid) COMMON / commonDPDFGridPts / nXPomZQsqGrid,xPomZQsqGrid COMMON / commonGluonPom / zGluonPomGrid COMMON / commonDPDFInitialised / DPDFinitialised IF (.NOT.DPDFinitialised) THEN WRITE(6,*) "Error in zGluonPom: DPDFs not initialised." STOP END IF xPomZQsq(1) = REAL(xPom) xPomZQsq(2) = REAL(z) xPomZQsq(3) = REAL(Qsq) zGluonPom = FINT(nCoord,xPomZQsq,nXPomZQsqGrid, & xPomZQsqGrid,zGluonPomGrid) RETURN END C-- ============================================================== SUBROUTINE GetDPDF(xPom,z,Qsq,PomOrReg,DPDF) C-- Purpose : Provide DPDFs at (xPom,z,Qsq). C-- Input : xPom,z,Qsq,PomOrReg. C-- Output : DPDF C-- ============================================================== IMPLICIT NONE INTEGER I CHARACTER*3 PomOrReg DOUBLE PRECISION xPom,z,Qsq,DPDF(-6:6),cReg, & zSingletPom,zGluonPom,zSinglet,zGluon,ReggeonFluxFactor,fReg, & pionPDF(-6:6),XMIN,XMAX,Q2MIN,Q2MAX COMMON / commonCReg / cReg COMMON / W50513 / XMIN,XMAX,Q2MIN,Q2MAX ! PDFLIB variables C-- Initialise DPDFs to zero: DO I = -6, 6 DPDF(I) = 0.D0 END DO IF (PomOrReg.EQ.'Pom') THEN zSinglet = zSingletPom(xPom,z,Qsq) zGluon = zGluonPom(xPom,z,Qsq) DO I = -3, 3 IF (I.EQ.0) THEN DPDF(I) = zGluon ELSE DPDF(I) = zSinglet/6.D0 END IF END DO ELSE IF (PomOrReg.EQ.'Reg') THEN IF ((z.GE.XMIN).AND.(z.LE.XMAX) & .AND.(Qsq.GE.Q2MIN).AND.(Qsq.LE.Q2MAX)) THEN fReg = ReggeonFluxFactor(xPom) CALL PFTOPDG(z,sqrt(Qsq),pionPDF) ! Get pion PDFs DO I = -6, 6 DPDF(I) = cReg*fReg*pionPDF(I) END DO END IF ELSE WRITE(6,*) "Error in GetDPDF: PomOrReg = ",PomOrReg END IF RETURN END C-- ====================================================== DOUBLE PRECISION FUNCTION ReggeonFluxFactor(xPom) C-- Purpose: Calculate non-perturbative Reggeon flux factor C-- Input : xPom C-- Output : ReggeonFluxFactor C-- ===================================================== IMPLICIT NONE DOUBLE PRECISION xPom,alphaRegPrime,BReg,alphaRegZero, + tmin,tcut,mp,xPomNorm,normFactor,tminNorm C-- Parameters same as H1 2006 DPDF fits. PARAMETER(alphaRegPrime=0.30D0,BReg=1.6D0,alphaRegZero=0.50D0) PARAMETER(mp=0.938272029D0,tcut=-1.D0) tmin = -(mp*xPom)**2/(1.D0-xPom) IF ((xPom.GT.0.99999D0).OR.(tmin.LE.tcut)) THEN ReggeonFluxFactor = 0.D0 RETURN END IF ReggeonFluxFactor = & (exp(BReg*tmin)*xPom**(-2.D0*tmin*alphaRegPrime) - & exp(BReg*tcut)*xPom**(-2.D0*tcut*alphaRegPrime)) * & xPom**(1.D0-2.D0*alphaRegZero)/ & (BReg-2.D0*alphaRegPrime*log(xPom)) C-- Normalise xPom*fReg = 1 at xPom = 0.003 (same as H1 2006 fit): xPomNorm = 0.003D0 tminNorm = -(mp*xPomNorm)**2/(1.D0-xPomNorm) normFactor = xPomNorm * & (exp(BReg*tminNorm)*xPomNorm**(-2.D0*tminNorm*alphaRegPrime) - & exp(BReg*tcut)*xPomNorm**(-2.D0*tcut*alphaRegPrime)) * & xPomNorm**(1.D0-2.D0*alphaRegZero)/ & (BReg-2.D0*alphaRegPrime*log(xPomNorm)) ReggeonFluxFactor=ReggeonFluxFactor/normFactor C-- At large xPom, above ReggeonFluxFactor goes negative. C-- When this is the case, simply set to zero. IF (ReggeonFluxFactor.LT.0.D0) ReggeonFluxFactor = 0.D0 RETURN END C-- ============================================================== SUBROUTINE GetDSF(xPom,beta,Qsq,whichDSF,DSF) C-- Purpose : Interpolate diffractive structure function. C-- Input : xPom,beta,Qsq,whichDSF C-- Output : DSF(3) C-- ============================================================== IMPLICIT NONE LOGICAL F2D3initialised,FLD3initialised, & F2D3charminitialised,FLD3charminitialised CHARACTER*9 whichDSF INTEGER nXPomGrid,nBetaGrid,nQsqGrid,nCoord PARAMETER(nXPomGrid=80,nBetaGrid=80,nQsqGrid=60,nCoord=3) INTEGER nXPomBetaQsqGrid(nCoord) DOUBLE PRECISION xPom,beta,Qsq,DSF(3) REAL FINT,xPomBetaQsq(nCoord), & xPomBetaQsqGrid(nXPomGrid+nBetaGrid+1+nQsqGrid), & F2D3ResPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & F2D3DirPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & F2D3RegGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & FLD3ResPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & FLD3DirPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & FLD3RegGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & F2D3charmResPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & F2D3charmDirPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & F2D3charmRegGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & FLD3charmResPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & FLD3charmDirPomGrid(nXPomGrid,nBetaGrid+1,nQsqGrid), & FLD3charmRegGrid(nXPomGrid,nBetaGrid+1,nQsqGrid) COMMON / commonDSFGridPts / nXPomBetaQsqGrid,xPomBetaQsqGrid COMMON / commonF2D3 / F2D3ResPomGrid,F2D3DirPomGrid,F2D3RegGrid COMMON / commonFLD3 / FLD3ResPomGrid,FLD3DirPomGrid,FLD3RegGrid COMMON / commonF2D3charm / F2D3charmResPomGrid, & F2D3charmDirPomGrid,F2D3charmRegGrid COMMON / commonFLD3charm / FLD3charmResPomGrid, & FLD3charmDirPomGrid,FLD3charmRegGrid COMMON / commonDSFinitialised / F2D3initialised,FLD3initialised, & F2D3charminitialised,FLD3charminitialised IF ((whichDSF.EQ.'F2D3').AND.(.NOT.F2D3initialised)) THEN WRITE(6,*) "Error in GetDSF: F2D3 not initialised." STOP ELSE IF ((whichDSF.EQ.'FLD3').AND.(.NOT.FLD3initialised)) THEN WRITE(6,*) "Error in GetDSF: FLD3 not initialised." STOP ELSE IF ((whichDSF.EQ.'F2D3charm').AND. & (.NOT.F2D3charminitialised)) THEN WRITE(6,*) "Error in GetDSF: F2D3charm not initialised." STOP ELSE IF ((whichDSF.EQ.'FLD3charm').AND. & (.NOT.FLD3charminitialised)) THEN WRITE(6,*) "Error in GetDSF: FLD3charm not initialised." STOP ELSE IF ((whichDSF.NE.'F2D3').AND.(whichDSF.NE.'FLD3').AND. & (whichDSF.NE.'F2D3charm').AND.(whichDSF.NE.'FLD3charm')) THEN WRITE(6,*) "Error in GetDSF: whichDSF = ",whichDSF STOP END IF xPomBetaQsq(1) = REAL(xPom) xPomBetaQsq(2) = REAL(beta) xPomBetaQsq(3) = REAL(Qsq) IF (whichDSF.EQ.'F2D3') THEN IF (F2D3initialised) THEN DSF(1) = FINT(nCoord,xPomBetaQsq,nXPomBetaQsqGrid, & xPomBetaQsqGrid,F2D3ResPomGrid) DSF(2) = FINT(nCoord,xPomBetaQsq,nXPomBetaQsqGrid, & xPomBetaQsqGrid,F2D3DirPomGrid) DSF(3) = FINT(nCoord,xPomBetaQsq,nXPomBetaQsqGrid, & xPomBetaQsqGrid,F2D3RegGrid) ELSE WRITE(6,*) "Error in GetDSF: F2D3 not initialised." STOP END IF ELSE IF (whichDSF.EQ.'FLD3') THEN IF (FLD3initialised) THEN DSF(1) = FINT(nCoord,xPomBetaQsq,nXPomBetaQsqGrid, & xPomBetaQsqGrid,FLD3ResPomGrid) DSF(2) = FINT(nCoord,xPomBetaQsq,nXPomBetaQsqGrid, & xPomBetaQsqGrid,FLD3DirPomGrid) DSF(3) = FINT(nCoord,xPomBetaQsq,nXPomBetaQsqGrid, & xPomBetaQsqGrid,FLD3RegGrid) ELSE WRITE(6,*) "Error in GetDSF: FLD3 not initialised." STOP END IF ELSE IF (whichDSF.EQ.'F2D3charm') THEN IF (F2D3charminitialised) THEN DSF(1) = FINT(nCoord,xPomBetaQsq,nXPomBetaQsqGrid, & xPomBetaQsqGrid,F2D3charmResPomGrid) DSF(2) = FINT(nCoord,xPomBetaQsq,nXPomBetaQsqGrid, & xPomBetaQsqGrid,F2D3charmDirPomGrid) DSF(3) = FINT(nCoord,xPomBetaQsq,nXPomBetaQsqGrid, & xPomBetaQsqGrid,F2D3charmRegGrid) ELSE WRITE(6,*) "Error in GetDSF: F2D3charm not initialised." STOP END IF ELSE IF (whichDSF.EQ.'FLD3charm') THEN IF (FLD3charminitialised) THEN DSF(1) = FINT(nCoord,xPomBetaQsq,nXPomBetaQsqGrid, & xPomBetaQsqGrid,FLD3charmResPomGrid) DSF(2) = FINT(nCoord,xPomBetaQsq,nXPomBetaQsqGrid, & xPomBetaQsqGrid,FLD3charmDirPomGrid) DSF(3) = FINT(nCoord,xPomBetaQsq,nXPomBetaQsqGrid, & xPomBetaQsqGrid,FLD3charmRegGrid) ELSE WRITE(6,*) "Error in GetDSF: FLD3charm not initialised." STOP END IF ELSE WRITE(6,*) "Error in GetDSF: whichDSF = ",whichDSF STOP END IF RETURN END