C-----------------------------------------------------------------------
#include "pdf/pilot.h"
#if defined(CERNLIB_SINGLE)
      FUNCTION Ctq4PnX (IPRTN, X, Q)
#endif
#if defined(CERNLIB_DOUBLE)
      DOUBLE PRECISION FUNCTION Ctq4PnX (IPRTN, X, Q)
#endif
C
C   Given the parton distribution function in the array Upd in
C   COMMON / CtqPar1 / , this routine fetches u(fl, x, q) at any value of
C   x and q using Mth-order polynomial interpolation for x and Ln(Q/Lambda).
C
#include "pdf/impdp.inc"
C
      PARAMETER (MXX = 105, MXQ = 25, MXF = 6)
      PARAMETER (MXPQX = (MXF *2 +2) * MXQ * MXX)
      PARAMETER (M= 2, M1 = M + 1)
C
      COMMON / W5051IC / Al, XV(0:MXX), QL(0:MXQ), UPD(MXPQX)
      COMMON / W5051ID / Nx, Nt, NfMx
      COMMON / W5051IE / Qini, Qmax, Xmin
C
      Dimension Fq(M1), Df(M1)
C                                                 Work with Log (Q)
      QG  = LOG (Q/AL)

C                           Find lower end of interval containing X
      JL = -1
      JU = Nx+1
 11   If (JU-JL .GT. 1) Then
         JM = (JU+JL) / 2
         If (X .GT. XV(JM)) Then
            JL = JM
         Else
            JU = JM
         Endif
         Goto 11
      Endif

      Jx = JL - (M-1)/2
      If (X .lt. Xmin) Then
C         Print '(A, 2(1pE12.4))', 
C     >     ' WARNING: X < Xmin, extrapolation used; X, Xmin =', X, Xmin
         If (Jx .LT. 0) Jx = 0
      Elseif (Jx .GT. Nx-M) Then
         Jx = Nx - M
      Endif
C                                    Find the interval where Q lies
      JL = -1
      JU = NT+1
 12   If (JU-JL .GT. 1) Then
         JM = (JU+JL) / 2
         If (QG .GT. QL(JM)) Then
            JL = JM
         Else
            JU = JM
         Endif
         Goto 12
      Endif

      Jq = JL - (M-1)/2
      If (Jq .LT. 0) Then
         Jq = 0
C         If (Q .lt. Qini)  Print '(A, 2(1pE12.4))', 
C     >     ' WARNING: Q < Qini, extrapolation used; Q, Qini =', Q, Qini
      Elseif (Jq .GT. Nt-M) Then
         Jq = Nt - M
C         If (Q .gt. Qmax)  Print '(A, 2(1pE12.4))', 
C     >     ' WARNING: Q > Qmax, extrapolation used; Q, Qmax =', Q, Qmax
      Endif

      If (Iprtn .GE. 3) Then
         Ip = - Iprtn
      Else
         Ip = Iprtn
      EndIf
C                             Find the off-set in the linear array Upd
      JFL = Ip + NfMx
      J0  = (JFL * (NT+1) + Jq) * (NX+1) + Jx
C
C                                           Now interpolate in x for M1 Q's
      Do 21 Iq = 1, M1
         J1 = J0 + (Nx+1)*(Iq-1) + 1
         Call Dpolin (XV(Jx), Upd(J1), M1, X, Fq(Iq), Df(Iq))
 21   Continue
C                                          Finish off by interpolating in Q
      Call Dpolin (QL(Jq), Fq(1), M1, QG, Ftmp, Ddf)

      Ctq4PnX = Ftmp
C
      RETURN
C                        ****************************
      END