C     Calculates LO alphaS using same input parameters
C     that were used to generate GS 2009 LO grids, i.e.
C     mcharm = 1.40GeV, mbot = 4.75GeV and alphaS(1Gev)
C     = 0.68183. Q must be input in GeV.

      FUNCTION GSALPS(Q)

      IMPLICIT NONE
      
      DOUBLE PRECISION BQCD
      DOUBLE PRECISION GSALPS, T
      DOUBLE PRECISION LLS
      DOUBLE PRECISION Q, MCharm, Mbot
      DOUBLE PRECISION B1
      EXTERNAL B1 
      DOUBLE PRECISION LLSINI
      DOUBLE PRECISION, SAVE :: LGLQSQ(3)
      LOGICAL, SAVE :: FSTTIM = .TRUE.
      INTEGER NF
      DOUBLE PRECISION LGMCSQ, LGMBSQ
      DOUBLE PRECISION PIFAC

      Mcharm = 1.40
      Mbot = 4.75

      T = LOG(Q*Q)

      IF(Q .LE. Mcharm) THEN
         NF = 3
      ELSEIF(Q .GT. Mcharm .AND. Q .LT. Mbot) THEN
         NF = 4
      ELSE 
         NF = 5
      END IF      

      IF(FSTTIM) THEN

         LGMCSQ = LOG(Mcharm**2)
         LGMBSQ = LOG(Mbot**2)
         PIFAC = ACOS(-1d0)

         LLSINI = LOG(1D0)-1.0/(B1(3.0D0)*0.68183d0)

         LGLQSQ(1) = LLSINI
         LGLQSQ(2) = (-LGMCSQ/(6.0*PIFAC)+B1(3.0D0)*LGLQSQ(1))/B1(4.0D0)
         LGLQSQ(3) = (-LGMBSQ/(6.0*PIFAC)+B1(4.0D0)*LGLQSQ(2))/B1(5.0D0)
 
         FSTTIM = .FALSE.

      END IF

      LLS = LGLQSQ(NF-2)

      BQCD = B1(NF*1d0)

      GSALPS = 1.0/(BQCD*(T-LLS))

      END FUNCTION GSALPS

      FUNCTION B1(NFREE)

      IMPLICIT NONE

      DOUBLE PRECISION B1, PIFAC, NFREE

      PIFAC = 3.14159265358979323846264338

      B1 = (33.0 - NFREE*2.0)/(12.0*PIFAC)

      END FUNCTION B1