CDECK  ID>, HWUALF.

*CMZ :-        -15/07/92  14.08.45  by  Mike Seymour

*-- Author :    Bryan Webber

C-----------------------------------------------------------------------

      FUNCTION HWUALF(IOPT,SCALE)

C-----------------------------------------------------------------------

C     STRONG COUPLING CONSTANT

C     IOPT.EQ.0  INITIALIZES

C         .EQ.1  TWO-LOOP, FLAVOUR THRESHOLDS

C         .EQ.2  RATIO OF ABOVE TO ONE-LOOP

C                WITH 5-FLAVOUR BETA, LAMBDA=QCDL3

C         .EQ.3  ONE-LOOP WITH 5-FLAVOUR BETA, LAMBDA=QCDL3

C-----------------------------------------------------------------------

      INCLUDE 'HERWIG61.INC'

      DOUBLE PRECISION HWUALF,SCALE,KAFAC,B3,B4,B5,B6,C3,C4,C5,C6,C35,

     & C45,C65,D35,RHO,RAT,RLF,DRH,EPS

      INTEGER IOPT,ITN

      SAVE B3,B4,B5,B6,C3,C4,C5,C6,C35,C45,C65,D35

      DATA EPS/1.D-6/

      IF (IOPT.EQ.0) THEN

C---INITIALIZE CONSTANTS

        CAFAC=FLOAT(NCOLO)

        CFFAC=FLOAT(NCOLO**2-1)/(2.*CAFAC)

        B3=((11.*CAFAC)- 6.)/(12.*PIFAC)

        B4=((11.*CAFAC)- 8.)/(12.*PIFAC)

        B5=((11.*CAFAC)-10.)/(12.*PIFAC)

        B6=((11.*CAFAC)-12.)/(12.*PIFAC)

        BETAF=6.*PIFAC*B5

        C3=((17.*CAFAC**2)-(5.*CAFAC+3.*CFFAC)*3.)/(24.*PIFAC**2)/B3**2

        C4=((17.*CAFAC**2)-(5.*CAFAC+3.*CFFAC)*4.)/(24.*PIFAC**2)/B4**2

        C5=((17.*CAFAC**2)-(5.*CAFAC+3.*CFFAC)*5.)/(24.*PIFAC**2)/B5**2

        C6=((17.*CAFAC**2)-(5.*CAFAC+3.*CFFAC)*6.)/(24.*PIFAC**2)/B6**2

        KAFAC=CAFAC*(67./18.-PIFAC**2/6.)-25./9.

C---QCDLAM IS 5-FLAVOUR LAMBDA-MS-BAR AT LARGE X OR Z

C---QCDL5  IS 5-FLAVOUR LAMBDA-MC

        QCDL5=QCDLAM*EXP(KAFAC/(4.*PIFAC*B5))/SQRT(2.D0)

C---COMPUTE THRESHOLD MATCHING

        RHO=2.*LOG(RMASS(6)/QCDL5)

        RAT=LOG(RHO)/RHO

        C65=(B5/(1.-C5*RAT)-B6/(1.-C6*RAT))*RHO

        RHO=2.*LOG(RMASS(5)/QCDL5)

        RAT=LOG(RHO)/RHO

        C45=(B5/(1.-C5*RAT)-B4/(1.-C4*RAT))*RHO

        RHO=2.*LOG(RMASS(4)/QCDL5)

        RAT=LOG(RHO)/RHO

        C35=(B4/(1.-C4*RAT)-B3/(1.-C3*RAT))*RHO+C45

C---FIND QCDL3

        D35=-1./(B3*C35)

        DO 10 ITN=1,100

          RAT=LOG(D35)/D35

          RLF=B3*D35/(1.-C3*RAT)

          DRH=B3*(RLF+C35)*D35**2/((1.-2.*C3*RAT+C3/D35)*RLF**2)

          D35=D35-DRH

          IF (ABS(DRH).LT.EPS*D35) GOTO 20

   10   CONTINUE

   20   QCDL3=QCDL5*EXP(0.5*D35)

      ENDIF

      IF (SCALE.LE.QCDL5) CALL HWWARN('HWUALF',51,*999)

      RHO=2.*LOG(SCALE/QCDL5)

      IF (IOPT.EQ.3) THEN

        IF (RHO.LE.D35) CALL HWWARN('HWUALF',52,*999)

        HWUALF=1./(B5*(RHO-D35))

        RETURN

      ENDIF

      RAT=LOG(RHO)/RHO

      IF (SCALE.GT.RMASS(6)) THEN

        RLF=B6*RHO/(1.-C6*RAT)+C65

      ELSEIF (SCALE.GT.RMASS(5)) THEN

        RLF=B5*RHO/(1.-C5*RAT)

      ELSEIF (SCALE.GT.RMASS(4)) THEN

        RLF=B4*RHO/(1.-C4*RAT)+C45

      ELSE

        RLF=B3*RHO/(1.-C3*RAT)+C35

      ENDIF

      IF (RLF.LE.ZERO) CALL HWWARN('HWUALF',53,*999)

      IF (IOPT.EQ.1) THEN

        HWUALF=1./RLF

      ELSE

        HWUALF=B5*(RHO-D35)/RLF

        IF (HWUALF.GT.ONE) CALL HWWARN('HWUALF',54,*999)

      ENDIF

      RETURN

 999  HWUALF=ZERO

      END