**********************************************************************
*                                                                    *
*   Program for generating Electromagnetic Structure Functions using *
*   Consistent Charm/Bottom  evolution according to MRRS:            *
*       A.D. Martin, R.G. Roberts, M. Ryskin and W.J. Stirling       *
*       University of Durham preprint DTP/96/102 (1996)              *
*   Setting MODE=1 corresponds to charm mass = 1.35 GeV              *
*   Setting MODE=2 corresponds to charm mass = 1.50 GeV              *
*   Setting MODE=3 corresponds to charm mass = 1.20 GeV              *
*   All modes have Lambda(MSbar,4) = 300MeV ==> alpha_S(M_Z) = 0.118 *
*                                                                    *
**********************************************************************

      implicit real*8(a-h,o-z)
      COMMON/GAUS96/XI(96),WI(96),XX(97),NTERMS
      COMMON/INPUT/alambda,flavor,qsct,qsdt,iord
      common/alfscl/sclfac

      call wate96
      iord=1
      flavor=3
      qsdt=7.29
      sclfac=0.0
      qsct=74.
      alambda=0.3
   99 continue
      write(6,*) 'input mode,x,q2:'
      read(5,*) mode,x,q2
      if(mode.eq.1) qsdt=7.29
      if(mode.eq.2) qsdt=9.00
      if(mode.eq.3) qsdt=5.76

      call sfun(x,q2,mode,f2p,flp,f1p,rp,f2n,fln,f1n,rn,
     xf2c,flc,f1c,f2b,flb,f1b)
c      scale=dsqrt(q2)
c      call mrschm(x,scale,mode,upv,dnv,usea,dsea,str,chm,bot,glu)

      print 100,x,q2,f2p,flp,f1p,f2n,fln,f1n
      print 101,f2c,flc,f1c,f2b,flb,f1b

      go to 99
    1 continue
  100 format(1x,'x= ',f10.6,' Q2= ',f10.2/
     11x,' Proton Structure Functions '/
     25x,' F2= ',f8.4,' FL= ',f8.4,' 2xF1=',f8.4/
     31x,' Neutron Structure Functions '/
     45x,' F2= ',f8.4,' FL= ',f8.4,' 2xF1=',f8.4/)
  101 format(1x,' Contributions from Charm'/
     25x,' F2(c)= ',f8.4,' FL(c)= ',f8.4,' 2xF1(c)=',f8.4/
     31x,' Contributions from Bottom'/
     45x,' F2(b)= ',f8.4,' FL(b)= ',f8.4,' 2xF1(b)=',f8.4/)

      end

c -- this is the corrected version of sfun (12/2/97)

      subroutine sfun(x,q2,mode,f2p,flp,f1p,rp,f2n,fln,f1n,rn,
     xf2c,flc,f1c,f2b,flb,f1b)
      implicit real*8(a-h,o-z)
      COMMON/GAUS96/XI(96),WI(96),XX(97),NTERMS
      COMMON/INPUT/alambda,flavor,qsct,qsdt,iord
      common/alfscl/sclfac
      data pi,pi2/3.14159,9.8696/
      XLAM=alambda
      IFL=4
      xlam2=xlam*xlam
      t=dlog(q2/xlam2)
      al=alpha(t)/(4.*pi)
      tt=t+sclfac
      argmin=qsdt/4./xlam2
      ttmin=dlog(argmin)
      if(tt.lt.ttmin) tt=ttmin
      altt=alpha(tt)/(4.*pi)
      scale=dsqrt(q2)
      q24=dexp(sclfac)*q2
      if(q24.le.qsdt/4.) q24=qsdt/4.
      scale4=dsqrt(q24)

      cf=4./3.
      ca=3.
      enf=flavor
      iorder=iord
      dpsi2=2./9.

      ww2=(1.-x)*q2/x
      epsc4=qsdt/q2
      fpsc4=qsdt/ww2
      epsc=epsc4/4.
      xcc=(1.+epsc)*x
      if(xcc.ge.0.95) xcc=0.95
      thcq=1.
      if(epsc.gt.1.) thcq=0.
      thcq4=1.
      if(epsc4.gt.1.) thcq4=0.
      threshc=thcq4*(1.-epsc4)
      thcw=1.
      if(fpsc4.gt.1.) thcw=0.
      epsb4=qsct/q2
      fpsb4=qsct/ww2
      epsb=epsb4/4.
      xbb=(1.+epsb)*x
      if(xbb.ge.0.95) xbb=0.95
      thbq=1.
      if(epsb.gt.1.) thbq=0.
      thbq4=1.
      if(epsb4.gt.1.) thbq4=0.
      threshb=thbq4*(1.-epsb4)
      thbw=1.
      if(fpsb4.gt.1.) thbw=0.

      call mrschm(x,scale,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
      fp=(4.*upv+dnv+8.*usea+2.*dsea+2.*str)/9.
      fn=(4.*dnv+upv+2.*usea+8.*dsea+2.*str)/9.
      fc=8.*chm/9.*threshc
      fb=2.*bot/9.*threshb

      call mrschm(xcc,scale,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
      fcc=(8.*chm/9.)*(1.+epsc4)*threshc
      if(q2.lt.qsdt/4.) fcc=0.
      call mrschm(xbb,scale,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
      fbb=(2.*bot/9.)*(1.+epsb4)*threshb
      if(q2.lt.qsct/4.) fbb=0.
      ffp=0.
      ffn=0.
      ffc=0.
      ffb=0.
      chmq=0.
      chmg=0.
      chmgsw=0.
      chmdcg=0.
      botq=0.
      botg=0.
      botgsw=0.
      botdcg=0.
      fflp=0.
      ffln=0.
      fflc=epsc4*fcc
      fflb=epsb4*fbb

      IF(IORD.LE.0.) THEN 
        GO TO 27
      ELSE
        GO TO 22
      ENDIF
   22 CONTINUE
      IF3=0

      FAC=FLAVOR
      facc=1.
      facb=1.

      IF(ifl.EQ.3.OR.ifl.EQ.4) then
      FAC=6./9.
      facc=4./9.
      facb=1./9.
      endif

      IF(IFL.EQ.1) IF3=1
      CF=4./3.
      AL1=dLOG(1.-X)

      ffp=Fp+Fp*AL*CF*(-9.-2.*PI2/3.+AL1*(-3.+2.*AL1))
      ffn=Fn+Fn*AL*CF*(-9.-2.*PI2/3.+AL1*(-3.+2.*AL1))
      ffc=fcc+fc*AL*CF*(-9.-2.*PI2/3.+AL1*(-3.+2.*AL1))
      ffb=fbb+fb*AL*CF*(-9.-2.*PI2/3.+AL1*(-3.+2.*AL1))

      DO 23 I=1,NTERMS
      Y=0.5*(1.-X)*XI(I)+0.5*(1.+X)
      XY=X/Y
      yc=y*(1.+epsc)
      yb=y*(1.+epsb)
      AL1=dLOG(1.-Y)
      call mrschm(xy,scale,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
      fpxy=(4.*upv+dnv+8.*usea+2.*dsea+2.*str)/9.
      fnxy=(4.*dnv+upv+2.*usea+8.*dsea+2.*str)/9.
      fcxy=8.*chm/9.*threshc
      fbxy=2.*bot/9.*threshb
      gluxy=glu
      call mrschm(xy,scale4,mode,upv,dnv,usea,dsea,str,chm,bot,glu)
      gluxy4=glu
      C22=CF*(6.+4.*Y-2.*(1.+Y*Y)/(1.-Y)*dLOG(Y)-2.*(1.+Y)*AL1
     2-IF3*2.*(1.+Y))
      C23=CF*(-3.+4.*AL1)/(1.-Y)
      CG2=2.*FAC*(-1.+8.*Y*(1.-Y)+(1.-2.*Y+2.*Y*Y)*dLOG(1./Y-1.))
      f1lq=4.*cf*y
      f1lg=8.*enf*y*(1.-y)

      wwz2=(1.-y)*q2/y
      fzpsc4=qsdt/wwz2
      thcwz=1.
      if(fzpsc4.gt.1.) thcwz=0.
      fzpsb4=qsct/wwz2
      thbwz=1.
      if(fzpsb4.gt.1.) thbwz=0.

      if(thcw.eq.0..or.thcwz.eq.0.) then
      cg21c=0.
      clg21c=0.
      go to 531
      endif
      betac2=1.-epsc4*y/(1.-y)
      betac=dsqrt(betac2)
      argc=(1.+betac)/(1.-betac)
      alc=dlog(argc)
      cg21c=2.*facc*(1.-2.*y+2.*y*y+epsc4*y*(1.-3.*y)
     x-8.*epsc*epsc*y*y)*alc+2.*facc*(8.*y*(1.-y)-1.
     x-epsc4*y*(1.-y))*betac
      clg21c=2.*facc*(-8.*epsc*y*y*alc+4.*y*(1.-y)*betac)
  531 argc=1./epsc
      if(argc.lt.1..or.yc.gt.1.) then
      cg22c=0.
      clg22c=0.
      go to 53
      endif
      alc=dlog(argc)
      cg22c=2.*facc*((1.-2.*yc+2.*yc*yc)*alc
     x+2.*yc*(1.-yc)*(1.-epsc))
      clg22c=cg22c
   53 cg22c=cg22c*thcq*(1.+epsc4)*threshc
      clg22c=clg22c*thcq*epsc4*threshc

      if(thbw.eq.0..or.thbwz.eq.0.) then
      cg21b=0.
      clg21b=0.
      go to 541
      endif
      betab2=1.-epsb4*y/(1.-y)
      betab=dsqrt(betab2)
      argb=(1.+betab)/(1.-betab)
      alb=dlog(argb)
      cg21b=2.*facb*(1.-2.*y+2.*y*y+epsb4*y*(1.-3.*y)
     x-8.*epsb*epsb*y*y)*alb+2.*facb*(8.*y*(1.-y)-1.
     x-epsb4*y*(1.-y))*betab
      clg21b=2.*facb*(-8.*epsb*y*y*alb+4.*y*(1.-y)*betab)
  541 argb=1./epsb     
      if(argb.lt.1..or.yb.gt.1.) then
      cg22b=0.
      clg22b=0.
      go to 54
      endif
      alb=dlog(argb)
      cg22b=2.*facb*((1.-2.*yb+2.*yb*yb)*alb
     x+2.*yb*(1.-yb)*(1.-epsb))
      clg22b=cg22b
   54 cg22b=cg22b*thbq*(1.+epsb4)*threshb
      clg22b=clg22b*thbq*epsb4*threshb

      ffp=ffp+.5*(1.-X)*WI(I)*AL*(C22*fpxy+C23*(fpxy-fp))
      ffn=ffn+.5*(1.-X)*WI(I)*AL*(C22*fnxy+C23*(fnxy-fn))
      ffc=ffc+.5*(1.-X)*WI(I)*AL*(C22*fcxy+C23*(fcxy-fc))
      ffb=ffb+.5*(1.-X)*WI(I)*AL*(C22*fbxy+C23*(fbxy-fb))

      fflp=fflp+.5*(1.-x)*wi(i)*al*f1lq*fpxy
      ffln=ffln+.5*(1.-x)*wi(i)*al*f1lq*fnxy
      fflc=fflc+.5*(1.-x)*wi(i)*al*f1lq*fcxy
      fflb=fflb+.5*(1.-x)*wi(i)*al*f1lq*fbxy

      IF(IFL-1) 23,23,24
   24 CONTINUE

      ffp=ffp+.5*(1.-X)*WI(I)*AL*CG2*gluxy
      ffn=ffn+.5*(1.-X)*WI(I)*AL*CG2*gluxy
      ffc=ffc+.5*(1.-X)*WI(I)*(ALtt*cg21c-altt*cg22c)*gluxy4
      ffb=ffb+.5*(1.-X)*WI(I)*(ALtt*cg21b-altt*cg22b)*gluxy4

      fflp=fflp+.5*(1.-x)*wi(i)*al*dpsi2*f1lg*gluxy
      ffln=ffln+.5*(1.-x)*wi(i)*al*dpsi2*f1lg*gluxy
      fflc=fflc+.5*(1.-x)*wi(i)*(altt*clg21c-altt*clg22c)*gluxy4
      fflb=fflb+.5*(1.-x)*wi(i)*(altt*clg21b-altt*clg22b)*gluxy4

   23 CONTINUE
   21 CONTINUE
      f2p=ffp+ffc+ffb
      f2n=ffn+ffc+ffb
      f2c=ffc
      f2b=ffb
      chmq=ffc-chmg
      botq=ffb-botg

      flp=fflp+fflc+fflb
      fln=ffln+fflc+fflb
      flc=fflc
      flb=fflb
      f1p=f2p-flp
      f1n=f2n-fln
      f1c=f2c-flc
      f1b=f2b-flb
      rp=flp/f1p
      rn=fln/f1n

   27 RETURN
      END

      DOUBLE PRECISION FUNCTION	ALPHA(T)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON/INPUT/alambda,flavor,qsct,qsdt,iord
      DATA PI/3.14159/
      DATA TOL/.0005/
      ITH=0
      TT=T
      qsdtt=qsdt/4.
      qsctt=qsct/4.
      AL=ALAMBDA
      AL2=AL*AL
      FLAV=4.
      QS=AL2*dEXP(T)

      if(qs.lt.0.5d0) then   !!  running stops below 0.5
          qs=0.5d0
          t=dlog(qs/al2)
          tt=t
      endif

      IF(QS.gt.QSCTT) GO	TO 12  
      IF(QS.lt.QSDTT) GO	TO 312  
   11 CONTINUE
      B0=11-2.*FLAV/3. 
      IF(IORD)1,1,2
c     IF(IORD)2,2,2	!TAKE CARE !!
    1 CONTINUE
      ALPHA=4.*PI/B0/T
      RETURN
    2 CONTINUE
      X1=4.*PI/B0
      B1=102.-38.*FLAV/3.
      X2=B1/B0**2
      AS=X1/T*(1.-X2*dLOG(T)/T)
    5 CONTINUE
      F=-T+X1/AS-X2*dLOG(X1/AS+X2)
      FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
      AS2=AS-F/FP
      DEL=ABS(F/FP/AS)
      IF(DEL-TOL)3,3,4
    3 CONTINUE
      ALPHA=AS2
      IF(ITH.EQ.0) RETURN
      GO TO (13,14,15) ITH
    4 CONTINUE
      AS=AS2
      GO TO 5
   12 ITH=1
      T=dLOG(QSCTT/AL2)
      GO TO 11
   13 ALFQC4=ALPHA
      FLAV=5.   
      ITH=2
      GO TO 11
   14 ALFQC5=ALPHA
      ITH=3
      T=TT
      GO TO 11
   15 ALFQS5=ALPHA
      ALFINV=1./ALFQS5+1./ALFQC4-1./ALFQC5
      ALPHA=1./ALFINV
      RETURN

  311 CONTINUE
      B0=11-2.*FLAV/3. 
      IF(IORD)31,31,32
c     IF(IORD)32,32,32	!TAKE CARE !!
   31 CONTINUE
      ALPHA=4.*PI/B0/T
      RETURN
   32 CONTINUE
      X1=4.*PI/B0
      B1=102.-38.*FLAV/3.
      X2=B1/B0**2
      AS=X1/T*(1.-X2*dLOG(T)/T)
   35 CONTINUE
      F=-T+X1/AS-X2*dLOG(X1/AS+X2)
      FP=-X1/AS**2*(1.-X2/(X1/AS+X2))
      AS2=AS-F/FP
      DEL=ABS(F/FP/AS)
      IF(DEL-TOL)33,33,34
   33 CONTINUE
      ALPHA=AS2
      IF(ITH.EQ.0) RETURN
      GO TO (313,314,315) ITH
   34 CONTINUE
      AS=AS2
      GO TO 35
  312 ITH=1
      T=dLOG(QSDTT/AL2)
      GO TO 311
  313 ALFQC4=ALPHA
      FLAV=3.   
      ITH=2
      GO TO 311
  314 ALFQC3=ALPHA
      ITH=3
      T=TT
      GO TO 311
  315 ALFQS3=ALPHA
      ALFINV=1./ALFQS3+1./ALFQC4-1./ALFQC3
      ALPHA=1./ALFINV
      RETURN
      END


      SUBROUTINE WATE96
C*******************************************************************
C*****							       *****
C***** THE X(I)	AND W(I) ARE THE DIRECT	OUTPUT FROM A PROGRAM  *****
C***** USING NAG ROUTINE D01BCF	TO CALCULATE THE	       *****
C***** GAUSS-LEGENDRE WEIGHTS FOR 96 POINT INTEGRATION.	       *****
C***** THEY AGREE TO TYPICALLY 14 DECIMAL PLACES WITH THE      *****
C***** TABLE IN	ABRAMOWITZ & STEGUN, PAGE 919.		       *****
C*****							       *****
C***** ---->   PETER HARRIMAN, APRIL 3RD 1990.		       *****
C*****							       *****
C*******************************************************************
      IMPLICIT REAL*8 (A-H,O-Z)
      DIMENSION	X(48),W(48)
      COMMON/GAUS96/XI(96),WI(96),XX(97),NTERMS
      NTERMS=96

      X( 1)=   0.01627674484960183561
      X( 2)=   0.04881298513604856015
      X( 3)=   0.08129749546442434360
      X( 4)=   0.11369585011066471632
      X( 5)=   0.14597371465489567682
      X( 6)=   0.17809688236761733390
      X( 7)=   0.21003131046056591064
      X( 8)=   0.24174315616383866556
      X( 9)=   0.27319881259104774468
      X(10)=   0.30436494435449495954
      X(11)=   0.33520852289262397655
      X(12)=   0.36569686147231213885
      X(13)=   0.39579764982890709712
      X(14)=   0.42547898840729897474
      X(15)=   0.45470942216774136446
      X(16)=   0.48345797392059470382
      X(17)=   0.51169417715466604391
      X(18)=   0.53938810832435567233
      X(19)=   0.56651041856139533470
      X(20)=   0.59303236477757022282
      X(21)=   0.61892584012546672523
      X(22)=   0.64416340378496526886
      X(23)=   0.66871831004391424358
      X(24)=   0.69256453664216964528
      X(25)=   0.71567681234896561582
      X(26)=   0.73803064374439816819
      X(27)=   0.75960234117664555964
      X(28)=   0.78036904386743123629
      X(29)=   0.80030874413913884180
      X(30)=   0.81940031073792957139
      X(31)=   0.83762351122818502758
      X(32)=   0.85495903343459936363
      X(33)=   0.87138850590929436968
      X(34)=   0.88689451740241818933
      X(35)=   0.90146063531585023110
      X(36)=   0.91507142312089592706
      X(37)=   0.92771245672230655266
      X(38)=   0.93937033975275308073
      X(39)=   0.95003271778443564022
      X(40)=   0.95968829144874048809
      X(41)=   0.96832682846326217918
      X(42)=   0.97593917458513455843
      X(43)=   0.98251726356301274934
      X(44)=   0.98805412632962202890
      X(45)=   0.99254390032376081654
      X(46)=   0.99598184298720747465
      X(47)=   0.99836437586317963722
      X(48)=   0.99968950388322870559
      W( 1)=   0.03255061449236316962
      W( 2)=   0.03251611871386883307
      W( 3)=   0.03244716371406427668
      W( 4)=   0.03234382256857594104
      W( 5)=   0.03220620479403026124
      W( 6)=   0.03203445623199267876
      W( 7)=   0.03182875889441101874
      W( 8)=   0.03158933077072719007
      W( 9)=   0.03131642559686137819
      W(10)=   0.03101033258631386231
      W(11)=   0.03067137612366917839
      W(12)=   0.03029991542082762553
      W(13)=   0.02989634413632842385
      W(14)=   0.02946108995816795100
      W(15)=   0.02899461415055528410
      W(16)=   0.02849741106508543861
      W(17)=   0.02797000761684838950
      W(18)=   0.02741296272602931385
      W(19)=   0.02682686672559184485
      W(20)=   0.02621234073567250055
      W(21)=   0.02557003600534944960
      W(22)=   0.02490063322248370695
      W(23)=   0.02420484179236479915
      W(24)=   0.02348339908592633665
      W(25)=   0.02273706965832950717
      W(26)=   0.02196664443874448477
      W(27)=   0.02117293989219144572
      W(28)=   0.02035679715433347898
      W(29)=   0.01951908114014518992
      W(30)=   0.01866067962741165898
      W(31)=   0.01778250231604547316
      W(32)=   0.01688547986424539715
      W(33)=   0.01597056290256253144
      W(34)=   0.01503872102699521608
      W(35)=   0.01409094177231515264
      W(36)=   0.01312822956696188190
      W(37)=   0.01215160467108866759
      W(38)=   0.01116210209983888144
      W(39)=   0.01016077053500880978
      W(40)=   0.00914867123078384552
      W(41)=   0.00812687692569928101
      W(42)=   0.00709647079115442616
      W(43)=   0.00605854550423662775
      W(44)=   0.00501420274292825661
      W(45)=   0.00396455433844564804
      W(46)=   0.00291073181793626202
      W(47)=   0.00185396078894924657
      W(48)=   0.00079679206555731759
      DO 1 I=1,48
      XI(I)=-X(49-I)
      WI(I)=W(49-I)
      XI(I+48)=X(I)
      WI(I+48)=W(I)
    1 CONTINUE
      DO 2 I=1,96
    2 XX(I)=0.5*(XI(I)+1.)
      XX(97)=1.0
      EXPON=1.0
      DO 3 I=1,96
      YI=2.*(0.5*(1.+XI(I)))**EXPON-1.
      WI(I)=WI(I)/(1.+XI(I))*(1.+YI)*EXPON
      XI(I)=YI
      XX(I)=0.5*(1.+YI)
    3 CONTINUE
      RETURN
      END
