[u/mrichter/AliRoot.git] / PDF / spdf / grsgalo.F

#include "pdf/pilot.h"
C-----------------------------------------------------------------------
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
*                                                                 *
*           G R S - LO - VIRTUAL PHOTON PARAMETRIZATIONS          *
*                                                                 *
*                 FOR A DETAILED EXPLANATION SEE                  *
*                M. GLUECK, E.REYA, M. STRATMANN :                *
*                    PHYS. REV. D51 (1995) 3220                   *
*                                                                 *
*   THE PARAMETRIZATIONS ARE FITTED TO THE EVOLVED PARTONS FOR    *
*        Q**2 / GEV**2  BETWEEN   0.6   AND  5.E4                 *
*                       AND (!)  Q**2 > 5 P**2                    *
*        P**2 / GEV**2  BETWEEN   0.0   AND  10.                  *
*                       P**2 = 0  <=> REAL PHOTON                 *
*             X         BETWEEN  1.E-4  AND   1.                  *
*                                                                 *
*   HEAVY QUARK THRESHOLDS  Q(H) = M(H)  IN THE BETA FUNCTION :   *
*                   M(C)  =  1.5,  M(B)  =  4.5                   *
*   CORRESPONDING LAMBDA(F) VALUES IN GEV FOR  Q**2 > M(H)**2 :   *
*      LO :   LAMBDA(3)  =  0.232,   LAMBDA(4)  =  0.200,         *
*             LAMBDA(5)  =  0.153,                                *
*   THE NUMBER OF ACTIVE QUARK FLAVOURS IS  NF = 3  EVERYWHERE    *
*   EXCEPT IN THE BETA FUNCTION, I.E. THE HEAVY QUARKS C,B,...    *
*   ARE NOT PRESENT AS PARTONS IN THE Q2-EVOLUTION.               *
*                                                                 *
*   PLEASE REPORT ANY STRANGE BEHAVIOUR TO :                      *
*          STRAT@HAL1.PHYSIK.UNI-DORTMUND.DE                      *
* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
*
*...INPUT PARAMETERS : 
*
*    X   = MOMENTUM FRACTION 
*    Q2  = SCALE Q**2 IN GEV**2
*    P2  = VIRTUALITY OF THE PHOTON IN GEV**2 
*
*...OUTPUT (ALWAYS X TIMES THE DISTRIBUTION DIVIDED BY ALPHA_EM) :
*...OUTPUT (ALWAYS X TIMES THE DISTRIBUTION) :  modified H.P.-B. 10.9.1996
*
********************************************************
      SUBROUTINE GRSGALO(DX,DQ2,DP2,
     +                     DUPV,DDNV,DUSEA,DDSEA,DSTR,DCHM,DBOT,DGL)
C      subroutine grsgalo(x,q2,p2,ugam,dgam,sgam,ggam)
#include "pdf/impdp.inc"
#if defined(CERNLIB_DOUBLE)
      DOUBLE PRECISION
#endif
#if defined(CERNLIB_SINGLE)
      REAL
#endif
     +          x, q2, p2, mu2, lam2,
     +          ugam, dgam, sgam, ggam,
     +          DUPV,DDNV,DUSEA,DDSEA,DSTR,DCHM,DBOT,DGL
C
      dimension u1(40),ds1(40),g1(40)
      dimension ud2(20),s2(20),g2(20)
      dimension up0(20),dsp0(20),gp0(20)
      DATA ALPHEM/7.29927D-3/ 
c
      data u1/-0.139d0,0.783d0,0.132d0,0.087d0,0.003d0,-0.0134d0,
     +   0.009d0,-0.017d0,0.092d0,-0.516d0,-0.085d0,0.439d0,
     +   0.013d0,0.108d0,-0.019d0,-0.272d0,-0.167d0,0.138d0,
     +   0.076d0,0.026d0,-0.013d0,0.27d0,0.107d0,-0.097d0,0.04d0,
     +   0.064d0,0.011d0,0.002d0,0.057d0,-0.057d0,0.162d0,
     +   -0.172d0,0.124d0,-0.016d0,-0.065d0,0.044d0,-1.009d0,
     +   0.622d0,0.227d0,-0.184d0/
      data ds1/0.033d0,0.007d0,-0.0516d0,0.12d0,0.001d0,-0.013d0,
     +   0.018d0,-0.028d0,0.102d0,-0.595d0,-0.114d0,0.669d0,
     +   0.022d0,0.001d0,-0.003d0,-0.0583d0,-0.041d0,0.035d0,
     +   0.009d0,0.009d0,0.004d0,0.054d0,0.025d0,-0.02d0,
     +   0.007d0,0.021d0,0.01d0,0.004d0,-0.067d0,0.06d0,-0.148d0,
     +   0.13d0,0.032d0,-0.009d0,-0.06d0,0.036d0,-0.39d0,0.033d0,
     +   0.245d0,-0.171d0/
      data g1/0.025d0,0.d0,-0.018d0,0.112d0,-0.025d0,0.177d0, 
     +   -0.022d0,0.024d0,0.001d0,-0.0104d0,0.d0,0.d0,-1.082d0,
     +   -1.666d0,0.d0,0.086d0,0.d0,0.053d0,0.005d0,-0.058d0,
     +   0.034d0,0.073d0,1.08d0,1.63d0,-0.0256d0,-0.088d0,0.d0,
     +   0.d0,-0.004d0,0.016d0,0.007d0,-0.012d0,0.01d0,-0.673d0,
     +   0.126d0,-0.167d0,0.032d0,-0.227d0,0.086d0,-0.159d0/
      data ud2/0.756d0,0.187d0,0.109d0,-0.163d0,0.002d0,0.004d0,
     +   0.054d0,-0.039d0,22.53d0,-21.02d0,5.608d0,0.332d0,
     +   -0.008d0,-0.021d0,0.381d0,0.572d0,4.774d0,1.436d0,
     +   -0.614d0,3.548d0/ 
      data s2/0.902d0,0.182d0,0.271d0,-0.346d0,0.017d0,-0.01d0,
     +   -0.011d0,0.0065d0,17.1d0,-13.29d0,6.519d0,0.031d0,
     +   -0.0176d0,0.003d0,1.243d0,0.804d0,4.709d0,1.499d0,
     +   -0.48d0,3.401d0/
      data g2/0.364d0,1.31d0,0.86d0,-0.254d0,0.611d0,0.008d0,
     +   -0.097d0,-2.412d0,-0.843d0,2.248d0,-0.201d0,1.33d0,
     +   0.572d0,0.44d0,1.233d0,0.009d0,0.954d0,1.862d0,3.791d0,
     +   -0.079d0/
      data up0/1.551d0,0.105d0,1.089d0,-0.172d0,3.822d0,-2.162d0,
     +   0.533d0,-0.467d0,-0.412d0,0.2d0,0.377d0,0.299d0,0.487d0,
     +   0.0766d0,0.119d0,0.063d0,7.605d0,0.234d0,-0.567d0,
     +   2.294d0/ 
      data dsp0/2.484d0,1.214d0,1.088d0,-0.1735d0,4.293d0,
     +   -2.802d0,0.5975d0,-0.1193d0,-0.0872d0,0.0418d0,0.128d0,
     +   0.0337d0,0.127d0,0.0135d0,0.14d0,0.0423d0,6.946d0,
     +   0.814d0,1.531d0,0.124d0/
      data gp0/1.682d0,1.1d0,0.5888d0,-0.4714d0,0.5362d0,0.0127d0,
     +   -2.438d0,0.03399d0,0.07825d0,0.05842d0,0.08393d0,2.348d0,
     +   -0.07182d0,1.084d0,0.3098d0,-0.07514d0,3.327d0,1.1d0,
     +   2.264d0,0.2675d0/
c
      save u1,ds1,g1,ud2,s2,g2,up0,dsp0,gp0
c
      x  = DX
      q  = SQRT(DQ2)
      q2 = DQ2
      p2 = DP2
      mu2=0.25d0
      lam2=0.232d0*0.232d0
c
      if(p2.le.0.25d0) then
         s=log(log(q2/lam2)/log(mu2/lam2))
         lp1=0.d0
         lp2=0.d0
      else
         if(q2.lt.p2) then
            write(*,1000)
 1000       format
     +      ('  WARNING: GRSGALO has been called with Q2 < P2 !',/,
     +       '           GRSGALO is about to blow up, therefore',/,
     +       '           Q2 is set equal to P2')
            q2=p2
         endif
         s=log(log(q2/lam2)/log(p2/lam2))
         lp1=log(p2/mu2)*log(p2/mu2)
         lp2=log(p2/mu2+log(p2/mu2))
      endif
c
      alp=up0(1)+lp1*u1(1)+lp2*u1(2)
      bet=up0(2)+lp1*u1(3)+lp2*u1(4)
      a=up0(3)+lp1*u1(5)+lp2*u1(6)+
     +  (up0(4)+lp1*u1(7)+lp2*u1(8))*s
      b=up0(5)+lp1*u1(9)+lp2*u1(10)+
     +  (up0(6)+lp1*u1(11)+lp2*u1(12))*s**0.5+
     +  (up0(7)+lp1*u1(13)+lp2*u1(14))*s**2
      gb=up0(8)+lp1*u1(15)+lp2*u1(16)+
     +  (up0(9)+lp1*u1(17)+lp2*u1(18))*s+
     +  (up0(10)+lp1*u1(19)+lp2*u1(20))*s**2
      ga=up0(11)+lp1*u1(21)+lp2*u1(22)+
     +  (up0(12)+lp1*u1(23)+lp2*u1(24))*s**0.5
      gc=up0(13)+lp1*u1(25)+lp2*u1(33)+
     +  (up0(14)+lp1*u1(26)+lp2*u1(34))*s
      gd=up0(15)+lp1*u1(27)+lp2*u1(35)+
     +  (up0(16)+lp1*u1(28)+lp2*u1(36))*s
      ge=up0(17)+lp1*u1(29)+lp2*u1(37)+
     +  (up0(18)+lp1*u1(30)+lp2*u1(38))*s
      gep=up0(19)+lp1*u1(31)+lp2*u1(39)+
     +  (up0(20)+lp1*u1(32)+lp2*u1(40))*s
      upart1=grsf2(x,s,alp,bet,a,b,ga,gb,gc,gd,ge,gep)
c
      alp=dsp0(1)+lp1*ds1(1)+lp2*ds1(2)
      bet=dsp0(2)+lp1*ds1(3)+lp2*ds1(4)
      a=dsp0(3)+lp1*ds1(5)+lp2*ds1(6)+
     +  (dsp0(4)+lp1*ds1(7)+lp2*ds1(8))*s
      b=dsp0(5)+lp1*ds1(9)+lp2*ds1(10)+
     +  (dsp0(6)+lp1*ds1(11)+lp2*ds1(12))*s**0.5+
     +  (dsp0(7)+lp1*ds1(13)+lp2*ds1(14))*s**2
      gb=dsp0(8)+lp1*ds1(15)+lp2*ds1(16)+
     +  (dsp0(9)+lp1*ds1(17)+lp2*ds1(18))*s+
     +  (dsp0(10)+lp1*ds1(19)+lp2*ds1(20))*s**2
      ga=dsp0(11)+lp1*ds1(21)+lp2*ds1(22)+
     +  (dsp0(12)+lp1*ds1(23)+lp2*ds1(24))*s
      gc=dsp0(13)+lp1*ds1(25)+lp2*ds1(33)+
     +  (dsp0(14)+lp1*ds1(26)+lp2*ds1(34))*s
      gd=dsp0(15)+lp1*ds1(27)+lp2*ds1(35)+
     +  (dsp0(16)+lp1*ds1(28)+lp2*ds1(36))*s
      ge=dsp0(17)+lp1*ds1(29)+lp2*ds1(37)+
     +  (dsp0(18)+lp1*ds1(30)+lp2*ds1(38))*s
      gep=dsp0(19)+lp1*ds1(31)+lp2*ds1(39)+
     +  (dsp0(20)+lp1*ds1(32)+lp2*ds1(40))*s
      dspart1=grsf2(x,s,alp,bet,a,b,ga,gb,gc,gd,ge,gep)
c
      alp=gp0(1)+lp1*g1(1)+lp2*g1(2)
      bet=gp0(2)+lp1*g1(3)+lp2*g1(4)
      a=gp0(3)+lp1*g1(5)+lp2*g1(6)+
     +  (gp0(4)+lp1*g1(7)+lp2*g1(8))*s**0.5
      b=gp0(5)+lp1*g1(9)+lp2*g1(10)+
     +  (gp0(6)+lp1*g1(11)+lp2*g1(12))*s**2
      gb=gp0(7)+lp1*g1(13)+lp2*g1(14)+
     +  (gp0(8)+lp1*g1(15)+lp2*g1(16))*s
      ga=gp0(9)+lp1*g1(17)+lp2*g1(18)+
     +  (gp0(10)+lp1*g1(19)+lp2*g1(20))*s**0.5+
     +  (gp0(11)+lp1*g1(21)+lp2*g1(22))*s**2
      gc=gp0(12)+lp1*g1(23)+lp2*g1(24)+
     +  (gp0(13)+lp1*g1(25)+lp2*g1(26))*s**2
      gd=gp0(14)+lp1*g1(27)+lp2*g1(28)+
     +  (gp0(15)+lp1*g1(29)+lp2*g1(30))*s+
     +  (gp0(16)+lp1*g1(31)+lp2*g1(32))*s**2
      ge=gp0(17)+lp1*g1(33)+lp2*g1(34)+
     +  (gp0(18)+lp1*g1(35)+lp2*g1(36))*s
      gep=gp0(19)+lp1*g1(37)+lp2*g1(38)+
     +  (gp0(20)+lp1*g1(39)+lp2*g1(40))*s
      gpart1=grsf2(x,s,alp,bet,a,b,ga,gb,gc,gd,ge,gep)
c
      s=log(log(q2/lam2)/log(mu2/lam2))
      suppr=1.d0/(1.d0+p2/0.59d0)**2
c
      alp=ud2(1)
      bet=ud2(2)
      a=ud2(3)+ud2(4)*s
      ga=ud2(5)+ud2(6)*s**0.5
      gc=ud2(7)+ud2(8)*s
      b=ud2(9)+ud2(10)*s+ud2(11)*s**2
      gb=ud2(12)+ud2(13)*s+ud2(14)*s**2
      gd=ud2(15)+ud2(16)*s
      ge=ud2(17)+ud2(18)*s
      gep=ud2(19)+ud2(20)*s
      udpart2=suppr*grsf1(x,s,alp,bet,a,b,ga,gb,gc,gd,ge,gep)
c
      alp=s2(1)
      bet=s2(2)
      a=s2(3)+s2(4)*s
      ga=s2(5)+s2(6)*s**0.5
      gc=s2(7)+s2(8)*s
      b=s2(9)+s2(10)*s+s2(11)*s**2
      gb=s2(12)+s2(13)*s+s2(14)*s**2
      gd=s2(15)+s2(16)*s
      ge=s2(17)+s2(18)*s
      gep=s2(19)+s2(20)*s
      spart2=suppr*grsf2(x,s,alp,bet,a,b,ga,gb,gc,gd,ge,gep)
c
      alp=g2(1)
      bet=g2(2)
      a=g2(3)+g2(4)*s**0.5
      b=g2(5)+g2(6)*s**2
      gb=g2(7)+g2(8)*s
      ga=g2(9)+g2(10)*s**0.5+g2(11)*s**2
      gc=g2(12)+g2(13)*s**2
      gd=g2(14)+g2(15)*s+g2(16)*s**2
      ge=g2(17)+g2(18)*s
      gep=g2(19)+g2(20)*s
      gpart2=suppr*grsf1(x,s,alp,bet,a,b,ga,gb,gc,gd,ge,gep)
c
      ugam=upart1+udpart2
      DUPV = UGAM * ALPHEM
      DUSEA = DUPV
      dgam=dspart1+udpart2
      DDNV = DGAM * ALPHEM
      DDSEA = DDNV
      sgam=dspart1+spart2
      DSTR = SGAM * ALPHEM
      ggam=gpart1+gpart2
      DGL = GGAM * ALPHEM
C
      DCHM = 0.D0
      DBOT = 0.D0
c
      return
      end
Commit	Line	Data
21886bb6	1	#include "pdf/pilot.h"
	2	C-----------------------------------------------------------------------
	3	* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
	4	* *
	5	* G R S - LO - VIRTUAL PHOTON PARAMETRIZATIONS *
	6	* *
	7	* FOR A DETAILED EXPLANATION SEE *
	8	* M. GLUECK, E.REYA, M. STRATMANN : *
	9	* PHYS. REV. D51 (1995) 3220 *
	10	* *
	11	* THE PARAMETRIZATIONS ARE FITTED TO THE EVOLVED PARTONS FOR *
	12	* Q2 / GEV2 BETWEEN 0.6 AND 5.E4 *
	13	* AND (!) Q2 > 5 P2 *
	14	* P2 / GEV2 BETWEEN 0.0 AND 10. *
	15	* P*2 = 0 <=> REAL PHOTON
	16	* X BETWEEN 1.E-4 AND 1. *
	17	* *
	18	* HEAVY QUARK THRESHOLDS Q(H) = M(H) IN THE BETA FUNCTION : *
	19	* M(C) = 1.5, M(B) = 4.5 *
	20	* CORRESPONDING LAMBDA(F) VALUES IN GEV FOR Q2 > M(H)2 : *
	21	* LO : LAMBDA(3) = 0.232, LAMBDA(4) = 0.200, *
	22	* LAMBDA(5) = 0.153, *
	23	* THE NUMBER OF ACTIVE QUARK FLAVOURS IS NF = 3 EVERYWHERE *
	24	* EXCEPT IN THE BETA FUNCTION, I.E. THE HEAVY QUARKS C,B,... *
	25	* ARE NOT PRESENT AS PARTONS IN THE Q2-EVOLUTION. *
	26	* *
	27	* PLEASE REPORT ANY STRANGE BEHAVIOUR TO : *
	28	* STRAT@HAL1.PHYSIK.UNI-DORTMUND.DE *
	29	* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
	30	*
	31	*...INPUT PARAMETERS :
	32	*
	33	* X = MOMENTUM FRACTION
	34	* Q2 = SCALE Q2 IN GEV2
	35	* P2 = VIRTUALITY OF THE PHOTON IN GEV**2
	36	*
	37	*...OUTPUT (ALWAYS X TIMES THE DISTRIBUTION DIVIDED BY ALPHA_EM) :
	38	*...OUTPUT (ALWAYS X TIMES THE DISTRIBUTION) : modified H.P.-B. 10.9.1996
	39	*
	40	********************************************************
	41	SUBROUTINE GRSGALO(DX,DQ2,DP2,
	42	+ DUPV,DDNV,DUSEA,DDSEA,DSTR,DCHM,DBOT,DGL)
	43	C subroutine grsgalo(x,q2,p2,ugam,dgam,sgam,ggam)
	44	#include "pdf/impdp.inc"
7ef50f50	45	#if defined(CERNLIB_DOUBLE)
	46	DOUBLE PRECISION
	47	#endif
	48	#if defined(CERNLIB_SINGLE)
	49	REAL
	50	#endif
21886bb6	51	+ x, q2, p2, mu2, lam2,
	52	+ ugam, dgam, sgam, ggam,
	53	+ DUPV,DDNV,DUSEA,DDSEA,DSTR,DCHM,DBOT,DGL
	54	C
	55	dimension u1(40),ds1(40),g1(40)
	56	dimension ud2(20),s2(20),g2(20)
	57	dimension up0(20),dsp0(20),gp0(20)
	58	DATA ALPHEM/7.29927D-3/
	59	c
	60	data u1/-0.139d0,0.783d0,0.132d0,0.087d0,0.003d0,-0.0134d0,
	61	+ 0.009d0,-0.017d0,0.092d0,-0.516d0,-0.085d0,0.439d0,
	62	+ 0.013d0,0.108d0,-0.019d0,-0.272d0,-0.167d0,0.138d0,
	63	+ 0.076d0,0.026d0,-0.013d0,0.27d0,0.107d0,-0.097d0,0.04d0,
	64	+ 0.064d0,0.011d0,0.002d0,0.057d0,-0.057d0,0.162d0,
	65	+ -0.172d0,0.124d0,-0.016d0,-0.065d0,0.044d0,-1.009d0,
	66	+ 0.622d0,0.227d0,-0.184d0/
	67	data ds1/0.033d0,0.007d0,-0.0516d0,0.12d0,0.001d0,-0.013d0,
	68	+ 0.018d0,-0.028d0,0.102d0,-0.595d0,-0.114d0,0.669d0,
	69	+ 0.022d0,0.001d0,-0.003d0,-0.0583d0,-0.041d0,0.035d0,
	70	+ 0.009d0,0.009d0,0.004d0,0.054d0,0.025d0,-0.02d0,
	71	+ 0.007d0,0.021d0,0.01d0,0.004d0,-0.067d0,0.06d0,-0.148d0,
	72	+ 0.13d0,0.032d0,-0.009d0,-0.06d0,0.036d0,-0.39d0,0.033d0,
	73	+ 0.245d0,-0.171d0/
	74	data g1/0.025d0,0.d0,-0.018d0,0.112d0,-0.025d0,0.177d0,
	75	+ -0.022d0,0.024d0,0.001d0,-0.0104d0,0.d0,0.d0,-1.082d0,
	76	+ -1.666d0,0.d0,0.086d0,0.d0,0.053d0,0.005d0,-0.058d0,
	77	+ 0.034d0,0.073d0,1.08d0,1.63d0,-0.0256d0,-0.088d0,0.d0,
	78	+ 0.d0,-0.004d0,0.016d0,0.007d0,-0.012d0,0.01d0,-0.673d0,
	79	+ 0.126d0,-0.167d0,0.032d0,-0.227d0,0.086d0,-0.159d0/
	80	data ud2/0.756d0,0.187d0,0.109d0,-0.163d0,0.002d0,0.004d0,
	81	+ 0.054d0,-0.039d0,22.53d0,-21.02d0,5.608d0,0.332d0,
	82	+ -0.008d0,-0.021d0,0.381d0,0.572d0,4.774d0,1.436d0,
	83	+ -0.614d0,3.548d0/
	84	data s2/0.902d0,0.182d0,0.271d0,-0.346d0,0.017d0,-0.01d0,
	85	+ -0.011d0,0.0065d0,17.1d0,-13.29d0,6.519d0,0.031d0,
	86	+ -0.0176d0,0.003d0,1.243d0,0.804d0,4.709d0,1.499d0,
	87	+ -0.48d0,3.401d0/
	88	data g2/0.364d0,1.31d0,0.86d0,-0.254d0,0.611d0,0.008d0,
	89	+ -0.097d0,-2.412d0,-0.843d0,2.248d0,-0.201d0,1.33d0,
	90	+ 0.572d0,0.44d0,1.233d0,0.009d0,0.954d0,1.862d0,3.791d0,
	91	+ -0.079d0/
	92	data up0/1.551d0,0.105d0,1.089d0,-0.172d0,3.822d0,-2.162d0,
	93	+ 0.533d0,-0.467d0,-0.412d0,0.2d0,0.377d0,0.299d0,0.487d0,
	94	+ 0.0766d0,0.119d0,0.063d0,7.605d0,0.234d0,-0.567d0,
	95	+ 2.294d0/
	96	data dsp0/2.484d0,1.214d0,1.088d0,-0.1735d0,4.293d0,
	97	+ -2.802d0,0.5975d0,-0.1193d0,-0.0872d0,0.0418d0,0.128d0,
	98	+ 0.0337d0,0.127d0,0.0135d0,0.14d0,0.0423d0,6.946d0,
	99	+ 0.814d0,1.531d0,0.124d0/
	100	data gp0/1.682d0,1.1d0,0.5888d0,-0.4714d0,0.5362d0,0.0127d0,
	101	+ -2.438d0,0.03399d0,0.07825d0,0.05842d0,0.08393d0,2.348d0,
	102	+ -0.07182d0,1.084d0,0.3098d0,-0.07514d0,3.327d0,1.1d0,
	103	+ 2.264d0,0.2675d0/
	104	c
	105	save u1,ds1,g1,ud2,s2,g2,up0,dsp0,gp0
	106	c
	107	x = DX
	108	q = SQRT(DQ2)
	109	q2 = DQ2
	110	p2 = DP2
	111	mu2=0.25d0
	112	lam2=0.232d0*0.232d0
	113	c
	114	if(p2.le.0.25d0) then
115	s=log(log(q2/lam2)/log(mu2/lam2))
116	lp1=0.d0
117	lp2=0.d0
118	else
119	if(q2.lt.p2) then
120	write(*,1000)
121	1000 format
122	+ (' WARNING: GRSGALO has been called with Q2 < P2 !',/,
123	+ ' GRSGALO is about to blow up, therefore',/,
124	+ ' Q2 is set equal to P2')
125	q2=p2
126	endif
127	s=log(log(q2/lam2)/log(p2/lam2))
128	lp1=log(p2/mu2)*log(p2/mu2)
129	lp2=log(p2/mu2+log(p2/mu2))
130	endif
131	c
132	alp=up0(1)+lp1u1(1)+lp2u1(2)
133	bet=up0(2)+lp1u1(3)+lp2u1(4)
134	a=up0(3)+lp1u1(5)+lp2u1(6)+
135	+ (up0(4)+lp1u1(7)+lp2u1(8))*s
136	b=up0(5)+lp1u1(9)+lp2u1(10)+
137	+ (up0(6)+lp1u1(11)+lp2u1(12))s*0.5+
138	+ (up0(7)+lp1u1(13)+lp2u1(14))s*2
139	gb=up0(8)+lp1u1(15)+lp2u1(16)+
140	+ (up0(9)+lp1u1(17)+lp2u1(18))*s+
141	+ (up0(10)+lp1u1(19)+lp2u1(20))s*2
142	ga=up0(11)+lp1u1(21)+lp2u1(22)+
143	+ (up0(12)+lp1u1(23)+lp2u1(24))s*0.5
144	gc=up0(13)+lp1u1(25)+lp2u1(33)+
145	+ (up0(14)+lp1u1(26)+lp2u1(34))*s
146	gd=up0(15)+lp1u1(27)+lp2u1(35)+
147	+ (up0(16)+lp1u1(28)+lp2u1(36))*s
148	ge=up0(17)+lp1u1(29)+lp2u1(37)+
149	+ (up0(18)+lp1u1(30)+lp2u1(38))*s
150	gep=up0(19)+lp1u1(31)+lp2u1(39)+
151	+ (up0(20)+lp1u1(32)+lp2u1(40))*s
152	upart1=grsf2(x,s,alp,bet,a,b,ga,gb,gc,gd,ge,gep)
153	c
154	alp=dsp0(1)+lp1ds1(1)+lp2ds1(2)
155	bet=dsp0(2)+lp1ds1(3)+lp2ds1(4)
156	a=dsp0(3)+lp1ds1(5)+lp2ds1(6)+
157	+ (dsp0(4)+lp1ds1(7)+lp2ds1(8))*s
158	b=dsp0(5)+lp1ds1(9)+lp2ds1(10)+
159	+ (dsp0(6)+lp1ds1(11)+lp2ds1(12))s*0.5+
160	+ (dsp0(7)+lp1ds1(13)+lp2ds1(14))s*2
161	gb=dsp0(8)+lp1ds1(15)+lp2ds1(16)+
162	+ (dsp0(9)+lp1ds1(17)+lp2ds1(18))*s+
163	+ (dsp0(10)+lp1ds1(19)+lp2ds1(20))s*2
164	ga=dsp0(11)+lp1ds1(21)+lp2ds1(22)+
165	+ (dsp0(12)+lp1ds1(23)+lp2ds1(24))*s
166	gc=dsp0(13)+lp1ds1(25)+lp2ds1(33)+
167	+ (dsp0(14)+lp1ds1(26)+lp2ds1(34))*s
168	gd=dsp0(15)+lp1ds1(27)+lp2ds1(35)+
169	+ (dsp0(16)+lp1ds1(28)+lp2ds1(36))*s
170	ge=dsp0(17)+lp1ds1(29)+lp2ds1(37)+
171	+ (dsp0(18)+lp1ds1(30)+lp2ds1(38))*s
172	gep=dsp0(19)+lp1ds1(31)+lp2ds1(39)+
173	+ (dsp0(20)+lp1ds1(32)+lp2ds1(40))*s
174	dspart1=grsf2(x,s,alp,bet,a,b,ga,gb,gc,gd,ge,gep)
175	c
176	alp=gp0(1)+lp1g1(1)+lp2g1(2)
177	bet=gp0(2)+lp1g1(3)+lp2g1(4)
178	a=gp0(3)+lp1g1(5)+lp2g1(6)+
179	+ (gp0(4)+lp1g1(7)+lp2g1(8))s*0.5
180	b=gp0(5)+lp1g1(9)+lp2g1(10)+
181	+ (gp0(6)+lp1g1(11)+lp2g1(12))s*2
182	gb=gp0(7)+lp1g1(13)+lp2g1(14)+
183	+ (gp0(8)+lp1g1(15)+lp2g1(16))*s
184	ga=gp0(9)+lp1g1(17)+lp2g1(18)+
185	+ (gp0(10)+lp1g1(19)+lp2g1(20))s*0.5+
186	+ (gp0(11)+lp1g1(21)+lp2g1(22))s*2
187	gc=gp0(12)+lp1g1(23)+lp2g1(24)+
188	+ (gp0(13)+lp1g1(25)+lp2g1(26))s*2
189	gd=gp0(14)+lp1g1(27)+lp2g1(28)+
190	+ (gp0(15)+lp1g1(29)+lp2g1(30))*s+
191	+ (gp0(16)+lp1g1(31)+lp2g1(32))s*2
192	ge=gp0(17)+lp1g1(33)+lp2g1(34)+
193	+ (gp0(18)+lp1g1(35)+lp2g1(36))*s
194	gep=gp0(19)+lp1g1(37)+lp2g1(38)+
195	+ (gp0(20)+lp1g1(39)+lp2g1(40))*s
196	gpart1=grsf2(x,s,alp,bet,a,b,ga,gb,gc,gd,ge,gep)
197	c
198	s=log(log(q2/lam2)/log(mu2/lam2))
199	suppr=1.d0/(1.d0+p2/0.59d0)**2
200	c
201	alp=ud2(1)
202	bet=ud2(2)
203	a=ud2(3)+ud2(4)*s
204	ga=ud2(5)+ud2(6)s*0.5
205	gc=ud2(7)+ud2(8)*s
206	b=ud2(9)+ud2(10)s+ud2(11)s**2
207	gb=ud2(12)+ud2(13)s+ud2(14)s**2
208	gd=ud2(15)+ud2(16)*s
209	ge=ud2(17)+ud2(18)*s
210	gep=ud2(19)+ud2(20)*s
211	udpart2=suppr*grsf1(x,s,alp,bet,a,b,ga,gb,gc,gd,ge,gep)
212	c
213	alp=s2(1)
214	bet=s2(2)
215	a=s2(3)+s2(4)*s
216	ga=s2(5)+s2(6)s*0.5
217	gc=s2(7)+s2(8)*s
218	b=s2(9)+s2(10)s+s2(11)s**2
219	gb=s2(12)+s2(13)s+s2(14)s**2
220	gd=s2(15)+s2(16)*s
221	ge=s2(17)+s2(18)*s
222	gep=s2(19)+s2(20)*s
223	spart2=suppr*grsf2(x,s,alp,bet,a,b,ga,gb,gc,gd,ge,gep)
224	c
225	alp=g2(1)
226	bet=g2(2)
227	a=g2(3)+g2(4)s*0.5
228	b=g2(5)+g2(6)s*2
229	gb=g2(7)+g2(8)*s
230	ga=g2(9)+g2(10)s0.5+g2(11)s**2
231	gc=g2(12)+g2(13)s*2
232	gd=g2(14)+g2(15)s+g2(16)s**2
233	ge=g2(17)+g2(18)*s
234	gep=g2(19)+g2(20)*s
235	gpart2=suppr*grsf1(x,s,alp,bet,a,b,ga,gb,gc,gd,ge,gep)
236	c
237	ugam=upart1+udpart2
238	DUPV = UGAM * ALPHEM
239	DUSEA = DUPV
240	dgam=dspart1+udpart2
241	DDNV = DGAM * ALPHEM
242	DDSEA = DDNV
243	sgam=dspart1+spart2
244	DSTR = SGAM * ALPHEM
245	ggam=gpart1+gpart2
246	DGL = GGAM * ALPHEM
247	C
248	DCHM = 0.D0
249	DBOT = 0.D0
250	c
251	return
252	end