text:
      SUBROUTINE HABER(TANB,SA,CA)

C  PROGRAM HMSUSY
C      WRITTEN BY HOWARD E. HABER
C      LATEST REVISION: AUGUST 28, 1995
C      COMMENTS OR QUESTIONS: SEND E-MAIL TO HABER@SCIPP.UCSC.EDU
C      BASED ON WORK IN COLLABORATION WITH R. HEMPFLING AND A. HOANG
C
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/ANG/ SINA,COSA,SIN2A,COS2A,SINBPA,COSBPA,SINBMA,COSBMA
      COMMON/HINT3/HLAA,HHAA,HLHLHL,HHHLHL,HHHHHL,HHHHHH,
     1      HHHPHM,HLHPHM
      COMMON/SQPARM/ SQM,SQK,SQU,SQD,XMU,STA2,STB2,SBA2,SBB2
      COMMON/QMASS/ TMPOLE
      DIMENSION H(4)

C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      COMMON/HMASS_HDEC/AMSM,AMA,AML,AMH,AMCH,AMAR
      COMMON/BREAK_HDEC/AMEL,AMER,AMSQ,AMUR,AMDR,AL,AU,AD,AMU,AM2
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

C
C   IRC=0 TREE LEVEL ONLY
C   IRC=-1 RADIATIVE CORRECTIONS INCLUDED WITHOUT RGE IMPROVEMENT
C   IRC=1 RADIATIVE CORRECTIONS INCLUDED WITH RGE IMPROVEMENT
C   IRC=2 SLIGHT IMPROVEMENT OF TREATMENT OF STOP-SBOTTOM SECTOR
C              WITH (IRC=-2: WITHOUT) RGE IMPROVEMENT
C   IRC=3 ELLIS, RIDOLFI, AND ZWIRNER LEADING MT**4 AND MB**4
C              CORRECTIONS WITH (IRC=-3: WITHOUT) RGE IMPROVEMENT
C
C     DO 838 IRC= -1,1,2
C
C   INPUT THE MASS OF THE CP-ODD SCALAR
C
C     DO 838 IH=2,2
C     IF (IH .EQ. 1) AM=   50.D0
C     IF (IH .EQ. 2) AM= 1000.D0
C
C   INPUT THE IMPORTANT PARAMETER TANB
C
C     DO 838 ITB= 1,1
C     IF (ITB .EQ. 1) TANB= 1.5D0
C     IF (ITB .EQ. 2) TANB= 20.D0
C     TANB= 10.D0*DFLOAT(ITB)
C
C   INPUT THE SQUARK MIXING A-PARAMETERS (B AND T SECTORS ONLY)
C   IN UNITS OF THE SUPERSYMMETRY BREAKING SCALE PARAMETER (SUSY)
C
C     DO 838 IAT=  0,0
C     AA= 0.4*DFLOAT(IAT)
C     DO 838 ISQ=0,9
C     SUSY= 200.D0+100.D0*DFLOAT(ISQ)
C     IF (ISQ .EQ. 9) SUSY= 2000.D0
C     AT= AA*SUSY
C     AB= AT
C
C   MU PARAMETER SET TO SIGNMU*SUSY
C   IN PRINCIPLE, CAN CHOOSE ANY MU, ALTHOUGH IF DABS(SIGNMU) IS
C         NOT 1, THEN CONTRIBUTIONS OF GAUGINOS AND HIGGSINOS
C         TO THE RADIATIVE CORRECTIONS MUST BE MODIFIED (SLIGHTLY)
C
C     DO  838 IMU= 1,1
C     SIGNMU= -1.D0*DFLOAT(IMU)
C
C   SQUARK MASS PARAMETERS SET IN SETUP SUBROUTINE
C   PARAMETERS DEPOSITED IN COMMON/SQPARM/
C      SQM= COMMON SQUARK/SLEPTON MASS (APART FROM B AND T SECTORS)
C      SQK= SQUARK SU(2)-DOUBLET SOFT-SUSY-BREAKING MASS
C      SQU= TOP-SQUARK SU(2)-SINGLET SOFT-SUSY-BREAKING MASS
C      SQD= BOTTOM-SQUARK SU(2)-SINGLET SOFT-SUSY-BREAKING MASS
C      SUSY= COMMON GAUGINO/HIGGSINO MASS SCALE
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      IRC= 2
      AM = AMA
      SUSY = AMSQ
      AT = SUSY
      AB = SUSY
      SIGNMU = 1
      IF(AMU.LT.0.D0)SIGNMU = -1
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      CALL SETUP(IRC,TANB,TM,SIGNMU,SUSY)
C
C   INITIALLY, ALL FIVE SUSY-BREAKING MASSES SET EQUAL.
C   THIS CAN BE CHANGED EITHER RIGHT HERE OR DIRECTLY IN THE
C   SETUP SUBROUTINE.
C
C   SUBROUTINE HSUSY COMPUTES HIGGS MASSES AND COUPLINGS
C   IT COMPUTES STOP AND SBOTTOM SPECTRUM FIRST, AND DEPOSITS THE
C   RESULTS IN COMMON/SQPARM/.  OUTPUT IS:
C        STA2,STB2= TOP-SQUARK SQUARED-MASSES
C        SBA2,SBB2= BOTTOM-SQUARK SQUARED-MASSES
C   NEXT IT COMPUTES THE HIGGS MASSES AND PUTS THEM IN THE ARRAY H()
C        H(1)= HEAVY CP-EVEN MASS
C        H(2)= LIGHT CP-EVEN MASS
C        H(3)= CP-ODD MASS (INPUT AS "AM" ABOVE)
C        H(4)= CHARGED HIGGS MASS
C   MIXING ANGLE FACTORS THAT APPEAR IN HIGGS FEYNMAN RULES
C   DEPOSITED IN  COMMON/ANG/.  TANB IS INPUT, WHILE OUTPUT INCLUDES:
C        SINA= SIN(ALPHA) [ALPHA: CP-EVEN HIGGS MIXING ANGLE]
C        COSA= COS(ALPHA)
C        SIN2A= SIN(2*ALPHA)
C        COS2A= COS(2*ALPHA)
C        SINBPA= SIN(BETA+ALPHA)
C        COSBPA= COS(BETA+ALPHA)
C        SINBMA= SIN(BETA-ALPHA)
C        COSBMA= COS(BETA-ALPHA)
C   THREE-HIGGS COUPLING DEPOSITED IN COMMON/HINT3/.  NOTATION:
C        HL= LIGHT CP-EVEN HIGGS
C        HH= HEAVY CP-EVEN HIGGS
C         A= CP-ODD HIGGS
C        HPHM= CHARGED HIGGS PAIR
C   FOR EXAMPLE, HLHPHM IS THE FEYNMAN RULE FOR THE INTERACTION OF
C   THE LIGHT CP-EVEN HIGGS WITH A CHARGED HIGGS PAIR (WITHOUT A
C   FACTOR OF I), HLHLHL IS THE THREE LIGHT CP-HIGGS COUPLING, ETC.
C
C   IERR=0 IF THE PROGRAM SUCCEEDS.  IF ONE FINDS A NEGATIVE
C   SQUARED MASS FOR ANY SCALAR (USUALLY, SCALAR-TOP, SCALAR-BOTTOM,
C   THE LIGHT CP-EVEN HIGGS OR THE CHARGED HIGGS), THE PROGRAM
C   RETURNS IERR>0.
C
      CALL HMSUSY(IRC,AM,TM,TANB,SUSY,AT,AB,H,ALF,IERR)
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      AML = H(2)
      AMH = H(1)
C     AMA = H(3)
      AMCH = H(4)
      AMAR = H(3)
      SA = SINA
      CA = COSA
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C
C   SLIGHT INACCURACIES OF THE PROGRAM:
C      1. IRC=2 SLIGHTLY IMPROVES NEUTRAL SCALAR MASSES, BUT
C            NO IMPROVEMENT YET FOR CHARGED SCALAR MASSES.
C      2. LEADING LOGS FOR CHARGED SCALAR MASS CORRECT ONLY
C            IF M(A)=O(ZMASS).
C      3. LEADING LOGS FOR HEAVY CP-EVEN SCALAR MASS CORRECT ONLY
C            IF M(A)=O(ZMASS).
C   IN CASES 2 AND 3, THE RELATIVE ERROR MADE IS RATHER MINOR,
C   SINCE HEAVY SCALARS ARE ROUGHLY DEGENERATE WITH THE CP-ODD SCALAR.
C      4. NO DETAILED STUDY YET OF ACCURACY OF RGE-IMPROVEMENT VIA
C            THE USE OF MT(Q) FOR WIDELY SPLIT SQK, SQU, AND SQD
C
      IF (STB2 .LT. 0.D0) GO TO 838
      IF (SBB2 .LT. 0.D0) GO TO 838
      STOP1= DSQRT(STA2)
      STOP2= DSQRT(STB2)
      SBOT1= DSQRT(SBA2)
      SBOT2= DSQRT(SBB2)
      IF (IERR .GT. 0) GO TO 838
C 500 WRITE(6,565) H(2),H(1),H(3),H(4),TANB,AA,SUSY,XMU
C 500 WRITE(6,565) H(2),H(3),TANB,AA,SUSY,XMU,H(4),STOPB
C 500 WRITE(6,564) IRC,H(2),H(3),TANB,TM,AT,SUSY,XMU
C 500 WRITE(6,565) XMU,H(3),H(2) ,H(4),TANB,AT,AB,SUSY
C     WRITE(6,565) SQK,SQU,SQD,XMU,STOP1,STOP2,SBOT1,SBOT2
  564 FORMAT(I2,7(F9.3,1X))
  565 FORMAT(8(F9.3,1X))
  838 CONTINUE
C     STOP
      RETURN
      END

      SUBROUTINE SETUP(IRC,TANB,TMASS,SIGNMU,SUSY)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON ZMASS,WMASS,SINW2,ALPHA,BMASS,G2,GP2,SINB,COSB,COTB
      COMMON/SQPARM/ SQM,SQK,SQU,SQD,XMU,STA2,STB2,SBA2,SBB2
      COMMON/QMASS/ TMPOLE
      DATA PI/3.1415926535D0/
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      COMMON/PARAM_HDEC/GF0,ALPH0,AMTAU0,AMMUON0,AMZ0,AMW0
      COMMON/MASSES_HDEC/AMS0,AMC0,AMB0,AMT0
      COMMON/BREAK_HDEC/AMEL,AMER,AMSQ,AMUR,AMDR,AL,AU,AD,AMU,AM2
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      TMPOLE= AMT0
      BMASS= AMB0
      ALS=ALPHAS_HDEC(TMPOLE,3)
      SINW2 = 1-AMW0**2/AMZ0**2
      ALPHA = DSQRT(2.D0)*GF0/PI*AMW0**2*SINW2
      ZMASS= AMZ0
      WMASS= AMW0
      SQM= AMEL
      SQK= AMSQ
      SQU= AMUR
      SQD= AMDR
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C     TMPOLE= 175.D0
C     ALS=ALFS(TMPOLE)
      ASCORR= 1.D0+4.D0*ALS/(3.D0*PI)
      TMASS= TMPOLE/ASCORR
C     TMASST= TMPOLE/ASCORR
C     TMASS= TMPOLE
C     IF (IRC .GT. 0.D0) TMASS= TMASST
C     BMASS= 4.7D0
C     ALPHA= 1.D0/128.D0
C     SINW2= 0.23D0
      COSW2= 1.D0-SINW2
      COSW= DSQRT(COSW2)
C     ZMASS= 91.1888D0
C     WMASS= ZMASS*COSW
      W2= WMASS**2
      E2= 4.D0*PI*ALPHA
      G2= E2/SINW2
      GP2= E2/COSW2
      SINB= TANB/DSQRT(1.D0+TANB**2)
      COSB= 1.D0/DSQRT(1.D0+TANB**2)
      COTB= 1.D0/TANB
C     SQM= SUSY
C     SQK= SUSY
C     SQU= SUSY
C     SQD= SUSY
      XMU= SIGNMU*SUSY
      RETURN
      END
C
C   IF CHARGED HIGGS MASS SQUARED IS NEGATIVE, PROGRAM RETURNS ZERO
C
      SUBROUTINE HMSUSY(IRC,AMASS,TMASS,TANB,SUSY,AT,AB,H,ALF,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON ZMASS,WMASS,SINW2,ALPHA,BMASS,G2,GP2,SINB,COSB,COTB
      COMMON/ANG/ SINA,COSA,SIN2A,COS2A,SINBPA,COSBPA,SINBMA,COSBMA
      COMMON/SQPARM/ SQM,SQK,SQU,SQD,XMU,STA2,STB2,SBA2,SBB2
      COMMON/HINT3/HLAA,HHAA,HLHLHL,HHHLHL,HHHHHL,HHHHHH,
     1      HHHPHM,HLHPHM
      COMMON/QMASS/ TMPOLE
      COMMON/TEST/ D11,D12,D22,DP
      DIMENSION H(4),DL(7)
      DATA PI/3.1415926535D0/, ZERO/0.D0/
      IERR= 0
      STEP= 0.D0
      IF (AMASS .GE. ZMASS) STEP= 1.D0
      Z2= ZMASS**2
      W2= WMASS**2
      A2= AMASS**2
      B2= BMASS**2
      COSW2= 1.D0-SINW2
      COS2W= COSW2-SINW2
      SINW= DSQRT(SINW2)
      COSW= DSQRT(COSW2)
      EU= 2.D0/3.D0
      ED= -1.D0/3.D0
      SINB2= SINB**2
      COSB2= COSB**2
      COS2B= COSB**2-SINB**2
      COS2B2= COS2B**2
      SIN2B= 2.D0*SINB*COSB
      TANB2= TANB**2
      COTB2= COTB**2
      S11= A2*SINB2+Z2*COSB2
      S22= A2*COSB2+Z2*SINB2
      S12= -(A2+Z2)*SINB*COSB
      SPM= A2+W2
      IF (IRC .NE. 0) GO TO 25
      CALL CPEVEN(S11,S22,S12,H1R,H2R,SIN2A,COS2A)
      DPM= SPM
      GO TO 50
C
C   HERE WE ALLOW FOR
C   MIXING OF Q(L) AND Q(R) ONLY FOR THE S-TOP AND S-BOTTOM
C   OTHERWISE, THERE IS NO MIXING.
C
   25 SUSY2= SUSY**2
      SQM2= SQM*SQM
      SQBL2= SQK**2+BMASS**2-Z2*COS2B*(0.5D0+ED*SINW2)
      SQBR2= SQD**2+BMASS**2+Z2*COS2B*ED*SINW2
      SQTL2= SQK**2+TMPOLE**2+Z2*COS2B*(0.5D0-EU*SINW2)
      SQTR2= SQU**2+TMPOLE**2+Z2*COS2B*EU*SINW2
      XB= AB-XMU*TANB
      XT= AT-XMU*COTB
      AB3= -XB*BMASS/(2.D0*WMASS)
      AT3= -XT*TMPOLE/(2.D0*WMASS)
      DIFFT= SQTL2-SQTR2
      DIFFB= SQBL2-SQBR2
      TMIX= -4.D0*AT3*WMASS
      BMIX= -4.D0*AB3*WMASS
      STA2= 0.5D0*(SQTL2+SQTR2+DSQRT(DIFFT**2+TMIX**2))
      STB2= 0.5D0*(SQTL2+SQTR2-DSQRT(DIFFT**2+TMIX**2))
      SBA2= 0.5D0*(SQBL2+SQBR2+DSQRT(DIFFB**2+BMIX**2))
      SBB2= 0.5D0*(SQBL2+SQBR2-DSQRT(DIFFB**2+BMIX**2))
      IF (SBB2 .LT. ZERO) GO TO 291
      IF (STB2 .LT. ZERO) GO TO 292
      SQT1= DSQRT(STA2)
      SQT2= DSQRT(STB2)
      SQB1= DSQRT(SBA2)
      SQB2= DSQRT(SBB2)
      SQTM2= SQT1*SQT2
      SQBM2= SQB1*SQB2
      SQTBM2= DSQRT(SQTM2*SQBM2)
      SMT2= 0.5D0*(SQK**2+SQU**2)
      SMT= DSQRT(SMT2)
      QP= DSQRT(SMT*TMPOLE)
      TMASS1= TMASSR(QP,TMASS)
      IF (IRC .LE. 0) TMASS1= TMASS
      T2= TMASS1**2
      TMASS2= TMASSR(SMT,TMASS)
      IF (IRC .LE. 0) TMASS2= TMASS
      COLOR= 3.D0
      GENS= 3.D0
      PT= COLOR*(1.D0-4.D0*EU*SINW2+8.D0*EU**2*SINW2**2)
      PB= COLOR*(1.D0+4.D0*ED*SINW2+8.D0*ED**2*SINW2**2)
      PF= (GENS-1.D0)*(PT+PB)+2.D0*GENS*(1.D0-2.D0*SINW2+4.D0*SINW2**2)
      P22= -2.D0*(32.D0*SINW2**2-54.D0*SINW2+27.D0)
      P12= -2.D0*(8.D0*SINW2**2-6.D0*SINW2-9.D0)
      P2H= -10.D0+2.D0*SINW2-2*SINW2**2
      P2HP= 8.D0-22.D0*SINW2+10.D0*SINW2**2
      P1H= -9.D0*COS2B2**2+(1.D0-2.D0*SINW2+2.D0*SINW2**2)*COS2B2
      PH12= (P1H-P2H)*STEP
      PH12P= (P1H+P2HP)*STEP
      FACT= COLOR*G2/(16.D0*PI**2*W2)
      FACT1= COLOR*G2*T2/(32.D0*PI**2*COSW2*SINB2)
      FACT2= COLOR*G2*B2/(32.D0*PI**2*COSW2*COSB2)
      FACTT= G2*Z2*SINB2/(96.D0*PI**2*COSW2)
      FACBB= G2*Z2*COSB2/(96.D0*PI**2*COSW2)
      FACTBT= G2*Z2/(96.D0*PI**2*COSW2)
      DLOGW= DLOG(SUSY2/W2)
      DLOGZ= DLOG(SUSY2/Z2)
      DLOGF= DLOG(SQM2/Z2)
      DLOGFW= DLOG(SQM2/W2)
C     DLOGT= DLOG(SUSY2/T2)
C     DLOGTB= DLOG(SUSY2/T2)
C     DLOGB= DLOG(SUSY2/Z2)
      DLOGT= DLOG(SQTM2/T2)
      DLOGTB= DLOG(SQTBM2/T2)
      DLOGB= DLOG(SQBM2/Z2)
      DLOGA= DLOG(A2/Z2)
      IF (IABS(IRC) .NE. 3) GO TO 30
      XB2= XB**2
      XT2= XT**2
      AB2= AB**2
      AT2= AT**2
      XMU2= XMU**2
      DLBF= FH(SBA2,SBB2)
      DLTF= FH(STA2,STB2)
      DLBG= FG(SBA2,SBB2)
      DLTG= FG(STA2,STB2)
      T4= TMASS2**4/SINB2
      B4= BMASS**4/COSB2
      DH11= FACT*(B4*(AB2*XB2*DLBG+2.D0*AB*XB*DLBF)+T4*XMU2*XT2*DLTG)
      DH22= FACT*(T4*(AT2*XT2*DLTG+2.D0*AT*XT*DLTF)+B4*XMU2*XB2*DLBG)
      DH12= -XMU*FACT*(B4*XB*(DLBF+AB*XB*DLBG)+T4*XT*(DLTF+AT*XT*DLTG))
      DM11= S11+DH11+2.D0*FACT*B4*DLOGB
      DM22= S22+DH22+2.D0*FACT*T2**2*DLOGT/SINB2
      DM12= S12+DH12
      GO TO 40
   30 DNL22= 0.D0
C
C   DNL22: MT**2 NON-LEADING LOG PIECE
C
      DNL22= FACT*T2*Z2/(3.D0*SINB2)
      CALL DELUD(IRC,TANB,TMASS1,AT,AB,DUD11,DUD22,DUD12,DUDPM)
      CALL DELHM(IRC,TANB,TMASS2,AT,AB,DHZ11,DHZ22,DHZ12,DHZPM)
      DM11= S11+DHZ11+DUD11
     1        +(FACT*(2.D0*B2**2/COSB2-Z2*B2)+FACBB*PB)*DLOGB
     2        +FACBB*(PT*DLOGT+PF*DLOGF+P22*DLOGZ+PH12*DLOGA)
      DM22= S22+DHZ22+DUD22+DNL22
     1        +(FACT*(2.D0*T2**2/SINB2-Z2*T2)+FACTT*PT)*DLOGT
     2        +FACTT*(PB*DLOGB+PF*DLOGF+P22*DLOGZ+PH12*DLOGA)
      DM12= S12+DHZ12+DUD12
     1        +SINB*COSB*((FACT1-FACTBT*PT)*DLOGT
     2        +(FACT2-FACTBT*PB)*DLOGB-FACTBT*(PF*DLOGF+P12*DLOGZ)
     3        +FACTBT*PH12P*DLOGA)
   40 CALL CPEVEN(DM11,DM22,DM12,H1R,H2R,SIN2A,COS2A)
      FACTG= COLOR*G2/(32.D0*PI**2)
      FACTG2= G2/(48.D0*PI**2)
      FACTGP= 5.D0*GP2*W2/(16.D0*PI**2)
      XCG= COLOR*(GENS-1.D0)+GENS
      DPML= SPM+FACTG*(2.D0*T2*B2/(W2*SINB2*COSB2)-T2/SINB2-B2/COSB2
     1         +2.D0*W2/3.D0)*DLOGTB + FACTGP*DLOGW
     2         +FACTG2*(XCG*DLOGFW-9.D0*DLOGW)
C     CORRQ=-0.5D0*T2*(W2+A2*COSB2-4.D0*B2/COSB2)/(W2*SINB2)
C    1     -(W2*(3.D0-5.D0*SINB2)/9.D0+2.D0*A2*COSB2/3.D0
C    2     -0.5D0*B2*(4.D0-5.D0*SINB2)/COSB2)/SINB2
C    3     -0.5D0*A2*(A2*COSB2/3.D0
C    4         -B2*(4.D0*COSB2+3.D0*SINB2**2)/(2.D0*COSB2))/(W2*SINB2)
C    5     +B2*(A2*SINB2*DLOG(T2/A2)-2.D0*B2*DLOG(T2/B2)/SINB2)
C    6         /(W2*COSB2)
C     CORRSQ=T2**3/(2.D0*W2*SINB2*SMT2)
C    1    +T2**2*(W2*(8.D0-5.D0*SINB2)-Z2*COS2B+A2*COSB2
C    2    -3.D0*B2*(2.D0+SINB2)/COSB2)/(6.D0*W2*SINB2*SMT2)
C     DPMNL= FACTG*(CORRQ+CORRSQ)
C     DPM= DPML+DPMNL+DHZPM
      DPM= DPML+DHZPM
   50 SIGN= 1.D0
      IF (SIN2A .LT. 0.D0) SIGN= -1.D0
C     COSA= DSQRT((1.D0+COS2A)/2.D0)
C     SINA= SIGN*DSQRT((1.D0-COS2A)/2.D0)
C     ALF= SIGN*DACOS(COSA)
      ALF= 0.5D0*ATAN2(SIN2A,COS2A)
      COSA= COS(ALF)
      SINA= SIN(ALF)
      BET= DATAN(TANB)
      SINBPA= DSIN(BET+ALF)
      SINBMA= DSIN(BET-ALF)
      COSBPA= DCOS(BET+ALF)
      COSBMA= DCOS(BET-ALF)
      H(1)= DSQRT(H1R)
      IF (H2R .LT. ZERO) GO TO 293
      H(2)= DSQRT(H2R)
      H(3)= AMASS
      H(4)= 0.D0
      IF (DPM .LT. ZERO) GO TO 294
      H(4)= DSQRT(DPM)
      CALL HLAMBG(IRC,TANB,AMASS,TMASS1,TMASS2,SUSY,AT,AB,DL)
      CALL HIGGS3(IRC,H,DL)
      RETURN
  291 IERR=1
      WRITE(6,991)
      RETURN
  292 IERR=2
      WRITE(6,992)
      RETURN
  293 IERR= 3
      WRITE(6,993)
      RETURN
  294 IERR= 4
      WRITE(6,994)
  991 FORMAT(1X,'ERROR: S-BOTTOM MASS SQUARED EIGENVALUE IS NEGATIVE')
  992 FORMAT(1X,'ERROR: S-TOP MASS SQUARED EIGENVALUE IS NEGATIVE')
  993 FORMAT(1X,'ERROR: LIGHT CP-EVEN SCALAR SQUARED-MASS IS NEGATIVE')
  994 FORMAT(1X,'ERROR: CHARGED HIGGS SQUARED MASS IS NEGATIVE')
      RETURN
      END

      SUBROUTINE CPEVEN(S11,S22,S12,H1,H2,SIN2A,COS2A)
      IMPLICIT REAL*8(A-H,O-Z)
      DISCR= DSQRT((S11-S22)**2+4.D0*S12**2)
      H1= 0.5D0*(S11+S22+DISCR)
      H2= 0.5D0*(S11+S22-DISCR)
      SIN2A= 2.D0*S12/DISCR
      COS2A= (S11-S22)/DISCR
      RETURN
      END

      SUBROUTINE DELHM(IRC,TANB,TMASS,AT,AB,DHZ11,DHZ22,DHZ12,DHZPM)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON ZMASS,WMASS,SINW2,ALPHA,BMASS,G2,GP2,SINB,COSB,COTB
      COMMON/SQPARM/ SQM,SQK,SQU,SQD,XMU,STA2,STB2,SBA2,SBB2
      DATA PI/3.1415926535D0/
      W2= WMASS*WMASS
      Z2= ZMASS*ZMASS
      AT2= AT*AT
      AB2= AB*AB
      XMU2= XMU*XMU
      GSQ= G2+GP2
      EU= 2.D0/3.D0
      ED= -1.D0/3.D0
      COLOR= 3.D0
      CB2= COSB**2
      SB2= SINB**2
      T2= TMASS**2
      B2= BMASS**2
      T4= TMASS**4/SB2
      B4= BMASS**4/CB2
      T2D= TMASS**2/SB2
      B2D= BMASS**2/CB2
      XB= AB-XMU*TANB
      XT= AT-XMU*COTB
      YB= AB+XMU*COTB
      YT= AT+XMU*TANB
      XB2= XB**2
      XT2= XT**2
      Q2= SQK**2
      U2= SQU**2
      D2= SQD**2
      UF= FF_HABER(U2,Q2)
      DF= FF_HABER(D2,Q2)
      UG= FG(Q2,U2)
      DG= FG(Q2,D2)
      UKP= UF+2.D0*EU*SINW2*(Q2-U2)*UG
      DKP= DF-2.D0*ED*SINW2*(Q2-D2)*DG
      UH= FH(Q2,U2)
      DH= FH(Q2,D2)
      UB= BP(Q2,U2)
      DB= BP(Q2,D2)
      FPBPU= UF+UB
      FPBPD= DF+DB
      FMBPU= UF-UB
      FMBPD= DF-DB
      HUD= HP(Q2,U2,D2)
      FUD= FP(Q2,U2,D2)
      FACT= COLOR*G2/(16.D0*PI**2*W2)
      DHZPM= 0.5D0*FACT*(XMU2*(T4*UF+B4*DF+2.D0*T2D*B2D*HUD)
     1       -T2D*B2D*(AB2*DF+AT2*UF+2.D0*AT*AB*HUD+FUD*(XMU2-AT*AB)**2)
     2       -W2*B2D*(XMU2*FPBPD-AB2*FMBPD)
     3       -W2*T2D*(XMU2*FPBPU-AT2*FMBPU))
      IF (IABS(IRC) .NE. 2) GO TO 10
      TQ2= STA2
      BQ2= SBA2
      U2= STB2
      D2= SBB2
      UF= FF_HABER(U2,TQ2)
      DF= FF_HABER(D2,BQ2)
      UG= FG(TQ2,U2)
      DG= FG(BQ2,D2)
      UKP= UF+2.D0*EU*SINW2*(TQ2-U2)*UG
      DKP= DF-2.D0*ED*SINW2*(BQ2-D2)*DG
      UH= FH(TQ2,U2)
      DH= FH(BQ2,D2)
      UB= BP(TQ2,U2)
      DB= BP(BQ2,D2)
   10 DHZ11= FACT*(B4*AB*XB*(2.D0*DH+AB*XB*DG)+T4*XMU2*XT2*UG
     1      +Z2*(B2*AB*(XB*DKP-AB*DB)+T2*XMU*COTB*(XT*UKP-XMU*COTB*UB)))
      DHZ22= FACT*(T4*AT*XT*(2.D0*UH+AT*XT*UG)+B4*XMU2*XB2*DG
     1      +Z2*(T2*AT*(XT*UKP-AT*UB)+B2*XMU*TANB*(XB*DKP-XMU*TANB*DB)))
      DHZ12= -FACT*(B4*XMU*XB*(DH+AB*XB*DG)+T4*XMU*XT*(UH+AT*XT*UG)
     1       -0.5D0*Z2*(T2*COTB*((XMU2+AT2)*UB-XT*YT*UKP)
     2                 +B2*TANB*((XMU2+AB2)*DB-XB*YB*DKP)))
      RETURN
      END

      SUBROUTINE DELUD(IRC,TANB,TMASS,AT,AB,DHZ11,DHZ22,DHZ12,DHZPM)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON ZMASS,WMASS,SINW2,ALPHA,BMASS,G2,GP2,SINB,COSB,COTB
      COMMON/SQPARM/ SQM,SQK,SQU,SQD,XMU,STA2,STB2,SBA2,SBB2
      DATA PI/3.1415926535D0/
      SINB2= SINB**2
      COSB2= COSB**2
      W2= WMASS*WMASS
      Z2= ZMASS*ZMASS
      T2= TMASS**2
      B2= BMASS**2
      Q2= SQK**2
      U2= SQU**2
      D2= SQD**2
      EU= 2.D0/3.D0
      ED= -1.D0/3.D0
      COLOR= 3.D0
      FACT= COLOR*G2*Z2/(16.D0*PI**2*W2)
      FACTG= COLOR*G2/(64.D0*PI**2)
      DHZPM= FACTG*(T2/SINB2+B2/COSB2-2.D0*W2/3.D0)*DLOG(Q2**2/(U2*D2))
      IF (IABS(IRC) .EQ. 2) GO TO 50
      DHZ11= -FACT*(B2*(0.5D0+2.D0*ED*SINW2)*DLOG(Q2/D2))
      DHZ22= -FACT*(T2*(0.5D0-2.D0*EU*SINW2)*DLOG(Q2/U2))
      DHZ12= 0.5D0*FACT*(B2*TANB*(0.5D0+2.D0*ED*SINW2)*DLOG(Q2/D2)
     1              +T2*COTB*(0.5D0-2.D0*EU*SINW2)*DLOG(Q2/U2))
      RETURN
   50 DHZ11= -FACT*(B2*(0.5D0+2.D0*ED*SINW2)*DLOG(SBA2/SBB2))
      DHZ22= -FACT*(T2*(0.5D0-2.D0*EU*SINW2)*DLOG(STA2/STB2))
      DHZ12= 0.5D0*FACT*(B2*TANB*(0.5D0+2.D0*ED*SINW2)*DLOG(SBA2/SBB2)
     1              +T2*COTB*(0.5D0-2.D0*EU*SINW2)*DLOG(STA2/STB2))
      END

      SUBROUTINE HLAMBG(IRC,TANB,AMASS,TM1,TM2,SUSY,AT,AB,DLT)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON ZMASS,WMASS,SINW2,ALPHA,BMASS,G2,GP2,SINB,COSB,COTB
      COMMON/SQPARM/ SQM,SQK,SQU,SQD,XMU,STA2,STB2,SBA2,SBB2
      COMMON/TEST/ D11,D12,D22,DP
      DIMENSION DL1(7),DL2(7),DL3(7),DL4(7),DLT(7)
      DATA PI/3.1415926535D0/
C
C   DLN(3)= LAMBDA(3)+LAMBDA(4),   N=1,...,4
C   DLN(4)= LAMBDA(5)-LAMBDA(4),   N=1,...,4
C
      GSQ= G2+GP2
      DLT(1)= 0.25D0*GSQ
      DLT(2)= 0.25D0*GSQ
      DLT(3)= -0.25D0*GSQ
      DLT(4)= 0.5D0*G2
      DLT(5)=0
      DLT(6)=0
      DLT(7)=0
      IF (IRC .EQ. 0) RETURN
      STEP= 0.D0
      IF (AMASS .GE. ZMASS) STEP= 1.D0
      A2= AMASS**2
      COSW2= 1.D0-SINW2
      T2= TM1**2
      B2= BMASS**2
      SUSY2= SUSY*SUSY
      SQM2= SQM*SQM
      SQTM2= DSQRT(STA2*STB2)
      SQBM2= DSQRT(SBA2*SBB2)
      W2= WMASS*WMASS
      Z2= ZMASS**2
      AT2= AT*AT
      AB2= AB*AB
      XMU2= XMU*XMU
      V2= 4.D0*W2/G2
      CB2= COSB**2
      SB2= SINB**2
      COS2B= COSB**2-SINB**2
      COS2B2= COS2B**2
      EU= 2.D0/3.D0
      ED= -1.D0/3.D0
      COLOR= 3.D0
      GENS= 3.D0
      PT= COLOR*(1.D0-4.D0*EU*SINW2+8.D0*EU**2*SINW2**2)
      PB= COLOR*(1.D0+4.D0*ED*SINW2+8.D0*ED**2*SINW2**2)
      PF= (GENS-1.D0)*(PT+PB)+2.D0*GENS*(1.D0-2.D0*SINW2+4.D0*SINW2**2)
      P22= -2.D0*(32.D0*SINW2**2-54.D0*SINW2+27.D0)
      P12= -2.D0*(8.D0*SINW2**2-6.D0*SINW2-9.D0)
      P2H= -10.D0+2.D0*SINW2-2*SINW2**2
      P2HP= 8.D0-22.D0*SINW2+10.D0*SINW2**2
      P1H= -9.D0*COS2B2**2+(1.D0-2.D0*SINW2+2.D0*SINW2**2)*COS2B2
      PH12= (P1H-P2H)*STEP
      PH12P= (P1H+P2HP)*STEP
      FACT= COLOR*G2/(16.D0*PI**2*W2*Z2)
      FACT1= COLOR*G2*T2/(32.D0*PI**2*COSW2*SB2*Z2)
      FACT2= COLOR*G2*B2/(32.D0*PI**2*COSW2*CB2*Z2)
      FACT3= G2/(96.D0*PI**2*COSW2)
      FACTC= COLOR*G2/(32.D0*PI**2)
      FACTC2= G2/(48.D0*PI**2)
      FACTPC= 5.D0*GP2/(16.D0*PI**2)
      DLOGW= DLOG(SUSY2/W2)
      DLOGZ= DLOG(SUSY2/Z2)
      DLOGF= DLOG(SQM2/Z2)
      DLOGFW= DLOG(SQM2/W2)
      DLOGT= DLOG(SQTM2/T2)
      DLOGB= DLOG(SQBM2/Z2)
      DLOGA= DLOG(A2/Z2)
      XCG= COLOR*(GENS-1.D0)+GENS
      DL1(1)= 0.25D0*GSQ*(1.D0
     1        +(FACT*(2.D0*B2**2/CB2-Z2*B2)/CB2+FACT3*PB)*DLOGB
     2        +FACT3*(PF*DLOGF+P22*DLOGZ+PT*DLOGT+PH12*DLOGA))
      DL1(2)= 0.25D0*GSQ*(1.D0
     1        +(FACT*(2.D0*T2**2/SB2-Z2*T2)/SB2+FACT3*PT)*DLOGT
     2        +FACT3*(PF*DLOGF+P22*DLOGZ+PB*DLOGB+PH12*DLOGA))
      DL1(3)= -0.25D0*GSQ*(1.D0
     1        -(FACT1-FACT3*PT)*DLOGT-(FACT2-FACT3*PB)*DLOGB
     2        +FACT3*(PF*DLOGF+P12*DLOGZ-PH12P*DLOGA))
      DL1(4)= 0.5D0*G2*(1.D0
     1        +FACTC*(2.D0*T2*B2/(W2**2*SB2*CB2)+2.D0/3.D0
     2        -(T2/SB2+B2/CB2)/W2)*DLOGT+FACTPC*DLOGW
     3        +FACTC2*(XCG*DLOGFW-9.D0*DLOGW))
      DL1(5)=0
      DL1(6)=0
      DL1(7)=0
      AFACT= COLOR/(16.D0*PI**2)
      HB2= 0.5D0*G2*B2/(W2*CB2)
      HT2= 0.5D0*G2*TM2**2/(W2*SB2)
      HB4= HB2*HB2
      HT4= HT2*HT2
      Q2= SQK**2
      U2= SQU**2
      D2= SQD**2
      UF= FF_HABER(U2,Q2)
      DF= FF_HABER(D2,Q2)
      UG= FG(Q2,U2)
      DG= FG(Q2,D2)
      UK= 0.25D0*(GSQ*UF+2.D0*EU*GP2*(Q2-U2)*UG)
      DK= 0.25D0*(GSQ*DF-2.D0*ED*GP2*(Q2-D2)*DG)
      UH= FH(Q2,U2)
      DH= FH(Q2,D2)
      UB= BP(Q2,U2)
      DB= BP(Q2,D2)
      A11= -AFACT*(HT2*XMU2*UB+HB2*AB2*DB)
      A22= -AFACT*(HT2*AT2*UB+HB2*XMU2*DB)
      DL2(1)= 0.5D0*GSQ*A11
      DL2(2)= 0.5D0*GSQ*A22
      DL2(3)= -0.25D0*GSQ*(A11+A22)
      DL2(4)= 0.5D0*G2*(A11+A22)
      DL2(5)= 0.D0
      DL2(6)= 0.D0
      DL2(7)= 0.D0
      DL3(1)= 2.D0*AFACT*(AB2*HB4*DH-HT2*XMU2*UK+HB2*AB2*DK)
      DL3(2)= 2.D0*AFACT*(AT2*HT4*UH-HB2*XMU2*DK+HT2*AT2*UK)
      DL3(3)= AFACT*(XMU2*(HB4*DH+HT4*UH)+HB2*(XMU2-AB2)*DK
     1           +HT2*(XMU2-AT2)*UK)
      DL3(4)= -0.5D0*AFACT*(4.D0*HB2*HT2*(AT*AB-XMU2)*HP(Q2,U2,D2)
     1           +HB2*DF*(2.D0*(HT2*AB2-HB2*XMU2)+G2*(XMU2-AB2))
     2           +HT2*UF*(2.D0*(HB2*AT2-HT2*XMU2)+G2*(XMU2-AT2)))
      DL3(5)= 0.D0
      DL3(6)= -AFACT*XMU*(HB4*AB*DH+HB2*AB*DK-HT2*AT*UK)
      DL3(7)= -AFACT*XMU*(HT4*AT*UH+HT2*AT*UK-HB2*AB*DK)
      DL4(1)= AFACT*(DG*(AB2*HB2)**2+UG*(XMU2*HT2)**2)
      DL4(2)= AFACT*(UG*(AT2*HT2)**2+DG*(XMU2*HB2)**2)
      DL4(3)= 2.D0*AFACT*XMU2*(AT2*HT4*UG+AB2*HB4*DG)
      DL4(4)= -AFACT*(HT2*HB2*(XMU2-AT*AB)**2*FP(Q2,U2,D2))
      DL4(5)= AFACT*XMU2*(AT2*HT4*UG+AB2*HB4*DG)
      DL4(6)= -AFACT*XMU*(AT*XMU2*HT4*UG+AB*AB2*HB4*DG)
      DL4(7)= -AFACT*XMU*(AB*XMU2*HB4*DG+AT*AT2*HT4*UG)
      DO 10 I=1,7
   10 DLT(I)= DL1(I)+DL2(I)+DL3(I)+DL4(I)
      D11= A2*SB2+V2*(DLT(1)*CB2+DLT(5)*SB2+2.D0*SINB*COSB*DLT(6))
      D22= A2*CB2+V2*(DLT(2)*SB2+DLT(5)*CB2+2.D0*SINB*COSB*DLT(7))
      D12= -A2*SINB*COSB+V2*(DLT(3)*SINB*COSB+DLT(6)*CB2+DLT(7)*SB2)
      DP= A2+0.5D0*V2*DLT(4)
      RETURN
      END
C
      FUNCTION FF_HABER(X2,Y2)
      IMPLICIT REAL*8(A-H,O-Z)
      DATA SMALL/1.D-4/
      FF_HABER= 0.D0
      IF (X2 .LE. 0.D0) GO TO 100
      IF (Y2 .LE. 0.D0) GO TO 100
      XYD= (X2-Y2)/Y2
      DXYD= DABS(XYD)
      IF (DXYD .LT. SMALL) GO TO 50
      FF_HABER= 1.D0/(X2-Y2)-X2*DLOG(X2/Y2)/(X2-Y2)**2
      RETURN
   50 FF_HABER= (-1.D0/2.D0+XYD/6.D0-XYD**2/12.D0+XYD**3/20.D0)/Y2
      RETURN
  100 WRITE(6,200)
  200 FORMAT(1X,'ERROR: ARGUMENTS OF FF-FUNCTION MUST BE POSITIVE')
      END
C
      FUNCTION FFD(X2,Y2)
      IMPLICIT REAL*8(A-H,O-Z)
      DATA SMALL/1.D-4/
      FFD= 0.D0
      IF (X2 .LE. 0.D0) GO TO 100
      IF (Y2 .LE. 0.D0) GO TO 100
      XYD= (X2-Y2)/Y2
      DXYD= DABS(XYD)
      IF (DXYD .LT. SMALL) GO TO 50
      FFD= 1.D0-X2*DLOG(X2/Y2)/(X2-Y2)
      RETURN
   50 FFD= -XYD/2.D0+XYD**2/6.D0-XYD**3/12.D0+XYD**4/20.D0
      RETURN
  100 WRITE(6,200)
  200 FORMAT(1X,'ERROR: ARGUMENTS OF FFD-FUNCTION MUST BE POSITIVE')
      END
C
      FUNCTION FG(X2,Y2)
      IMPLICIT REAL*8(A-H,O-Z)
      DATA SMALL/1.D-4/
      FG= 0.D0
      IF (X2 .LE. 0.D0) GO TO 100
      IF (Y2 .LE. 0.D0) GO TO 100
      XYD= (X2-Y2)/Y2
      DXYD= DABS(XYD)
      IF (DXYD .LT. SMALL) GO TO 50
      FG= 2.D0/(X2-Y2)**2-(X2+Y2)*DLOG(X2/Y2)/(X2-Y2)**3
      RETURN
   50 FG= (-1.D0/6.D0+XYD/6.D0-3.D0*XYD**2/20.D0
     1           +2.D0*XYD**3/15.D0)/Y2**2
      RETURN
  100 WRITE(6,200)
  200 FORMAT(1X,'ERROR: ARGUMENTS OF FG-FUNCTION MUST BE POSITIVE')
      END
C
      FUNCTION FH(X2,Y2)
      IMPLICIT REAL*8(A-H,O-Z)
      DATA SMALL/1.D-4/
      FH= 0.D0
      IF (X2 .LE. 0.D0) GO TO 100
      IF (Y2 .LE. 0.D0) GO TO 100
      XYD= (X2-Y2)/Y2
      DXYD= DABS(XYD)
      IF (DXYD .LT. SMALL) GO TO 50
      FH= DLOG(X2/Y2)/(X2-Y2)
      RETURN
   50 FH= (1.D0-XYD/2.D0+XYD**2/3.D0-XYD**3/4.D0+XYD**4/5.D0)/Y2
      RETURN
  100 WRITE(6,200)
  200 FORMAT(1X,'ERROR: ARGUMENTS OF FH-FUNCTION MUST BE POSITIVE')
      END
C
      FUNCTION BP(X2,Y2)
      IMPLICIT REAL*8(A-H,O-Z)
      DATA SMALL/1.D-4/
      BP= 0.D0
      IF (X2 .LE. 0.D0) GO TO 100
      IF (Y2 .LE. 0.D0) GO TO 100
      XYD= (X2-Y2)/Y2
      DXYD= DABS(XYD)
      IF (DXYD .LT. SMALL) GO TO 50
      BP= 0.5D0*(X2+Y2-2.D0*X2*Y2*DLOG(X2/Y2)/(X2-Y2))/(X2-Y2)**2
      RETURN
   50 BP= (1.D0/6.D0-XYD/12.D0+XYD**2/20.D0-XYD**3/30.D0)/Y2
      RETURN
  100 WRITE(6,200)
  200 FORMAT(1X,'ERROR: ARGUMENTS OF BP-FUNCTION MUST BE POSITIVE')
      END
C
      FUNCTION FP(Z2,X2,Y2)
      IMPLICIT REAL*8(A-H,O-Z)
      DATA SMALL/1.D-4/
      FP= 0.D0
      IF (X2 .LE. 0.D0) GO TO 100
      IF (Y2 .LE. 0.D0) GO TO 100
      XYD= (X2-Y2)/Y2
      DXYD= DABS(XYD)
      IF (DXYD .LT. SMALL) GO TO 50
      FP= (FF_HABER(X2,Z2)-FF_HABER(Y2,Z2))/(X2-Y2)
      RETURN
   50 FP= -FG(Y2,Z2)+XYD*FP1(Y2,Z2)+XYD**2*FP2(Y2,Z2)+XYD**3*FP3(Y2,Z2)
      RETURN
  100 WRITE(6,200)
  200 FORMAT(1X,'ERROR: ARGUMENTS OF FP-FUNCTION MUST BE POSITIVE')
      END
C
      FUNCTION HP(Z2,X2,Y2)
      IMPLICIT REAL*8(A-H,O-Z)
      DATA SMALL/1.D-4/
      HP= 0.D0
      IF (X2 .LE. 0.D0) GO TO 100
      IF (Y2 .LE. 0.D0) GO TO 100
      XYD= (X2-Y2)/Y2
      DXYD= DABS(XYD)
      IF (DXYD .LT. SMALL) GO TO 50
      HP= (FFD(X2,Z2)-FFD(Y2,Z2))/(X2-Y2)
      RETURN
   50 HP= FF_HABER(Z2,Y2)+XYD*BP(Y2,Z2)+XYD**2*HP2(Y2,Z2)
     .    +XYD**3*HP3(Y2,Z2)
      RETURN
  100 WRITE(6,200)
  200 FORMAT(1X,'ERROR: ARGUMENTS OF HP-FUNCTION MUST BE POSITIVE')
      END
C
      FUNCTION FP1(Y2,Z2)
      IMPLICIT REAL*8(A-H,O-Z)
      DATA SMALL/1.D-4/
      YZD= (Y2-Z2)/Z2
      DYZD= DABS(YZD)
      IF (DYZD .LT. SMALL) GO TO 50
      FP1= 0.5D0*(5.D0*Y2+Z2)/(Y2-Z2)**3
     1         -Y2*(2.D0*Z2+Y2)*DLOG(Y2/Z2)/(Y2-Z2)**4
      RETURN
   50 FP1=(-1.D0/12.D0+YZD/15.D0-YZD**2/20.D0+4.D0*YZD**3/105.D0)/Z2**2
      RETURN
      END
C
      FUNCTION FP2(Y2,Z2)
      IMPLICIT REAL*8(A-H,O-Z)
      DATA SMALL/1.D-4/
      YZD= (Y2-Z2)/Z2
      DYZD= DABS(YZD)
      IF (DYZD .LT. SMALL) GO TO 50
      FP2= (Z2**2-8.D0*Z2*Y2-17.D0*Y2**2)/(6.D0*(Y2-Z2)**4)
     1         +Y2**2*(3.D0*Z2+Y2)*DLOG(Y2/Z2)/(Y2-Z2)**5
      RETURN
   50 FP2=(1.D0/20.D0-YZD/30.D0+3.D0*YZD**2/140.D0+YZD**3/70.D0)/Z2**2
      RETURN
      END

      FUNCTION FP3(Y2,Z2)
      IMPLICIT REAL*8(A-H,O-Z)
      DATA SMALL/1.D-4/
      YZD= (Y2-Z2)/Z2
      DYZD= DABS(YZD)
      IF (DYZD .LT. SMALL) GO TO 50
      POLY= Z2**3-7.D0*Z2**2*Y2+29.D0*Z2*Y2**2+37.D0*Y2**3
      FP3= POLY/(12.D0*(Y2-Z2)**5)
     1         -Y2**3*(4.D0*Z2+Y2)*DLOG(Y2/Z2)/(Y2-Z2)**6
      RETURN
   50 FP3=(-1.D0/30.D0+2.D0*YZD/105.D0-3.D0*YZD**2/280.D0
     1         +2.D0*YZD**3/315.D0)/Z2**2
      RETURN
      END
C
      FUNCTION HP2(Y2,Z2)
      IMPLICIT REAL*8(A-H,O-Z)
      DATA SMALL/1.D-4/
      YZD= (Y2-Z2)/Z2
      DYZD= DABS(YZD)
      IF (DYZD .LT. SMALL) GO TO 50
      HP2= (Z2**2-5.D0*Z2*Y2-2.D0*Y2**2)/(6.D0*(Y2-Z2)**3)
     1         +Y2**2*Z2*DLOG(Y2/Z2)/(Y2-Z2)**4
      RETURN
   50 HP2=(-1.D0/12.D0+YZD/30.D0-YZD**2/60.D0+YZD**3/105.D0)/Z2
      RETURN
      END
C
      FUNCTION HP3(Y2,Z2)
      IMPLICIT REAL*8(A-H,O-Z)
      DATA SMALL/1.D-4/
      YZD= (Y2-Z2)/Z2
      DYZD= DABS(YZD)
      IF (DYZD .LT. SMALL) GO TO 50
      POLY= Z2**3-5.D0*Z2**2*Y2+13.D0*Z2*Y2**2+3.D0*Y2**3
      HP3= POLY/(12.D0*(Y2-Z2)**4)
     1         -Y2**3*Z2*DLOG(Y2/Z2)/(Y2-Z2)**5
      RETURN
   50 HP3=(1.D0/20.D0-YZD/60.D0+YZD**2/140.D0-YZD**3/280.D0)/Z2
      RETURN
      END

      FUNCTION ALFS(QK)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON/QMASS/ TMPOLE
      DATA PI/3.1415926535D0/
C
C   ONE-LOOP STRONG COUPLING FOR QK>BMASS
C
      AST= 0.11D0
      ALFS= AST
      IF (QK .GE. TMPOLE) F= 6.D0
      IF (QK .LT. TMPOLE) F= 5.D0
      IF (QK .LT. 5.D0) GO TO 10
      DEN= 1.D0+AST*(11.D0-2.D0*F/3.D0)*DLOG(QK/TMPOLE)/(2.D0*PI)
      ALFS= AST/DEN
      RETURN
   10 WRITE (6,100)
  100 FORMAT(1X,'ERROR: STRONG COUPLING CONSTANT EVALUATED BELOW M(B)')
      RETURN
      END

      FUNCTION TMASSR(Q,TMASS)
      IMPLICIT REAL*8 (A-H,O-Z)
      COMMON ZMASS,WMASS,SINW2,ALPHA,BMASS,G2,GP2,SINB,COSB,COTB
      COMMON/QMASS/ TMPOLE
      DATA PI/3.1415926535/
      QP= TMPOLE
      HB2= 0.5D0*G2*BMASS**2/(WMASS*COSB)**2
      HT2= 0.5D0*G2*TMASS**2/(WMASS*SINB)**2
      ALS= ALFS(QP)
      ALST= ALS-3.D0*HT2/(64.D0*PI)
     1          +GP2/(48.D0*PI)
     2          -HB2/(64.D0*PI)
      TMASSR= TMASS*(1.D0-ALST*DLOG(Q**2/TMASS**2)/PI)
      RETURN
      END

      SUBROUTINE HIGGS3(IRC,H,DL)
      IMPLICIT REAL*8(A-H,O-Z)
      COMMON ZMASS,WMASS,SINW2,ALPHA,BMASS,G2,GP2,SINB,COSB,COTB
      COMMON/ANG/ SINA,COSA,SIN2A,COS2A,SINBPA,COSBPA,SINBMA,COSBMA
      COMMON/HINT3/HLAA,HHAA,HLHLHL,HHHLHL,HHHHHL,HHHHHH,
     1      HHHPHM,HLHPHM
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      DOUBLE PRECISION LA1,LA2,LA3,LA4,LA5,LA6,LA7
      COMMON/HSELF_HDEC/LA1,LA2,LA3,LA4,LA5,LA6,LA7
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      DIMENSION H(4),DL(7)
      DATA PI/3.1415926535D0/
C
C   DL(3)= LAMBDA(3)+LAMBDA(4)
C   DL(4)= LAMBDA(5)-LAMBDA(4)
C
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      LA1 = DL(1)
      LA2 = DL(2)
      LA4 = DL(5)-DL(4)
      LA3 = DL(3)-LA4
      LA5 = DL(5)
      LA6 = DL(6)
      LA7 = DL(7)
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      V2= 4.D0*WMASS**2/G2
      DL4NL= 2.D0*(H(4)**2-H(3)**2)/V2
      COSW= DSQRT(1.D0-SINW2)
      G= DSQRT(G2)
      SINB2= SINB**2
      COSB2= COSB**2
      SINA2= SINA**2
      COSA2= COSA**2
      SINA3= SINA**3
      COSA3= COSA**3
      COS2B= COSB2-SINB2
      IF (IRC .EQ. 0) GO TO 50
      FACN= 2.D0*WMASS/G
      HLAA= FACN*(DL(1)*SINA*COSB*SINB2-DL(2)*COSA*SINB*COSB2
     1        -(DL(3)+DL(5))*(COSA*SINB**3-SINA*COSB**3)
     2        +2.D0*DL(5)*SINBMA-DL(6)*SINB*(COSB*SINBPA+SINA*COS2B)
     3        -DL(7)*COSB*(COSA*COS2B-SINB*SINBPA))
      HHAA= -FACN*(DL(1)*COSA*COSB*SINB2+DL(2)*SINA*SINB*COSB2
     1        +(DL(3)+DL(5))*(SINA*SINB**3+COSA*COSB**3)
     2        -2.D0*DL(5)*COSBMA-DL(6)*SINB*(COSB*COSBPA+COSA*COS2B)
     3        +DL(7)*COSB*(SINB*COSBPA+SINA*COS2B))
      HHHLHL= -3.D0*FACN*(DL(1)*COSA*SINA2*COSB+DL(2)*SINA*COSA2*SINB
     1        +(DL(3)+DL(5))*(SINA3*SINB+COSA3*COSB-2.D0*COSBMA/3.D0)
     2        -DL(6)*SINA*(COSB*COS2A+COSA*COSBPA)
     3        +DL(7)*COSA*(SINB*COS2A+SINA*COSBPA))
      HHHHHL= 3.D0*FACN*(DL(1)*SINA*COSB*COSA2-DL(2)*COSA*SINB*SINA2
     1        +(DL(3)+DL(5))*(SINA3*COSB-COSA3*SINB+2.D0*SINBMA/3.D0)
     2        -DL(6)*COSA*(COSB*COS2A-SINA*SINBPA)
     3        -DL(7)*SINA*(SINB*COS2A+COSA*SINBPA))
      HHHHHH= -3.D0*FACN*(DL(1)*COSB*COSA3+DL(2)*SINB*SINA3
     1        +(DL(3)+DL(5))*SINA*COSA*SINBPA
     2        +DL(6)*COSA2*(3.D0*SINA*COSB+COSA*SINB)
     3        +DL(7)*SINA2*(3.D0*COSA*SINB+SINA*COSB))
      HLHLHL= 3.D0*FACN*(DL(1)*COSB*SINA3-DL(2)*SINB*COSA3
     1        +(DL(3)+DL(5))*SINA*COSA*COSBPA
     2        -DL(6)*SINA2*(3.D0*COSA*COSB-SINA*SINB)
     3        +DL(7)*COSA2*(3.D0*SINA*SINB-COSA*COSB))
      HHHPHM= HHAA-FACN*COSBMA*DL(4)
      HLHPHM= HLAA-FACN*SINBMA*DL(4)
C
C   THESE COUPLINGS INCLUDE NON-LEADING LOG CONTRIBUTIONS
C   WHICH ARE TAKEN TO BE THE SAME AS IN THE CORRECTIONS
C   TO THE CHARGED HIGGS MASS SQUARED.
C
C     HHHPM= HHAA-FACN*COSBMA*DL4NL
C     HLHPM= HLAA-FACN*SINBMA*DL4NL
      RETURN
   50 FAC= -G*WMASS
      FAN= -0.5D0*G*ZMASS/COSW
      HHHPHM= FAC*COSBMA-FAN*COS2B*COSBPA
      HLHPHM= FAC*SINBMA+FAN*COS2B*SINBPA
      HHHHHH= 3.D0*FAN*COS2A*COSBPA
      HLHLHL= 3.D0*FAN*COS2A*SINBPA
      HHHHHL= -FAN*(2.D0*SIN2A*COSBPA+COS2A*SINBPA)
      HHHLHL= FAN*(2.D0*SIN2A*SINBPA-COS2A*COSBPA)
      HHAA= -FAN*COS2B*COSBPA
      HLAA= FAN*COS2B*SINBPA
      RETURN
      END


text:
program main
use iso_fortran_env, only: iostat_end
implicit none
integer :: iu,i,ierr,j
character (len=*), parameter :: fmt_i = "(i0)"
character (len=100)          :: text
open (newunit=iu,file="temp.txt")
write (iu,fmt_i) (i,i=1,3)
rewind (iu) ! move position to beginning of file
write (iu,fmt_i) (2*i,i=1,3)
! deletes previous data and writes 2 4 6 on successive lines
! position is now after last line written
backspace (iu) ! move file position to last line written
read (iu,*) j
print*,j ! 6
backspace (iu) ; backspace (iu) ; backspace (iu)
read (iu,*) j
print*,j ! 2
read (iu,*) ! skip a line -- move file position by 1 line
read (iu,*) j
print*,j ! 6
rewind (iu)
write (iu,"(i0/i0/a)") 10,20,"abc"
rewind (iu)
do
   read (iu,*,iostat=ierr) j
   if (ierr == 0) then
      print*,"j =",j
   else if (ierr == iostat_end) then
      stop "reached end of file"
   else
      backspace (iu)
      read (iu,"(a)") text
      stop "could not read integer from '" // trim(text) // "'"
   end if
end do
end program main
! output:
!            6
!            2
!            6
!  j =          10
!  j =          20
! STOP could not read integer from 'abc'

max_stars_repo_name: awvwgk/FortranTip

max_issues_repo_name: awvwgk/FortranTip

max_forks_repo_name: awvwgk/FortranTip


SUBMODULE (vtk_io) vtk_io_implementation
    !! author: Ian Porter
    !! date: 12/1/2017
    !!
    !! This module implements the ability to read/write VTK formatted files
    !!
    CONTAINS

        MODULE PROCEDURE vtk_legacy_write
        USE Misc,     ONLY : to_uppercase
        USE vtk_vars, ONLY : default_fn, default_title, filetype, vtkfilename, vtktitle, ascii, binary, &
          &                  version, fcnt, file_extension
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/1/2017
        !!
        !! This subroutines writes the legacy vtk output file
        !!
        INTEGER(i4k)     :: i, inputstat
        LOGICAL          :: file_is_open
        CHARACTER(LEN=8) :: fcnt_char = ''
        CHARACTER(LEN=:), ALLOCATABLE :: form, filetype_text

        IF (PRESENT(data_type)) filetype    = data_type    !! Calling program provided what file type to use for vtk file
        IF (PRESENT(filename)) THEN
            vtkfilename = filename                         !! Calling program provided a filename
        ELSE
            vtkfilename = default_fn                       !! Calling program did not provide a filename. Use default
        END IF
        IF (PRESENT(multiple_io)) THEN
            IF (multiple_io) THEN
                WRITE (fcnt_char,FMT='(i8)') fcnt
                vtkfilename = vtkfilename(1:INDEX(to_uppercase(vtkfilename),to_uppercase(file_extension))-1) // "_" // &
                  &           TRIM(ADJUSTL(fcnt_char)) // file_extension
                fcnt = fcnt + 1
            END IF
        END IF
        IF (PRESENT(title)) THEN
            vtktitle    = title                            !! Calling program provided a title
        ELSE
            vtktitle    = default_title                    !! Calling program did not provide a title. Use default
        END IF

        INQUIRE(unit = unit, opened = file_is_open)        !! Check to see if file is already open
        IF (.NOT. file_is_open) THEN                       !! File is not yet open. Determine format to open file
            SELECT CASE (filetype)
            CASE (ascii)
                form          = 'formatted'
                filetype_text = 'ASCII'
            CASE (binary)
                form          = 'unformatted'
                filetype_text = 'BINARY'
            CASE DEFAULT
                WRITE(*,*) 'Warning: filetype is incorrectly defined. Will default to ASCII'
                form          = 'formatted'
                filetype_text = 'ASCII'
            END SELECT
            OPEN(unit=unit, file=vtkfilename, iostat=inputstat, status='unknown', form=form)
                                                           !! Open the VTK file
        END IF

        WRITE(unit,100) version                            !! VTK version (currently, 3.0)
        WRITE(unit,100) vtktitle                           !! VTK title card
        WRITE(unit,100) filetype_text                      !! VTK file type

        CALL geometry%write(unit)                          !! Write the geometry information
        IF (PRESENT(celldatasets)) THEN
            WRITE(unit,101) celldatasets(1)%n
            DO i = 1, SIZE(celldatasets)
                CALL celldatasets(i)%attribute%write(unit) !! Write the cell data values
            END DO
        ELSE IF (PRESENT(celldata)) THEN
            WRITE(unit,101) celldatasets(1)%n
            CALL celldata%write(unit)                      !! Write the cell data values
        END IF
        IF (PRESENT(pointdatasets)) THEN
            WRITE(unit,102) pointdatasets(1)%n
            DO I = 1, SIZE(pointdatasets)
                CALL pointdatasets(i)%attribute%write(unit)
            END DO
        ELSE IF (PRESENT(pointdata)) THEN
            WRITE(unit,102) pointdatasets(1)%n
            CALL pointdata%write(unit)                     !! Write the point data values
        END IF

        CLOSE(unit)                                        !! Close the VTK file
100     FORMAT(a)
101     FORMAT('CELL_DATA ',i0)
102     FORMAT('POINT_DATA ',i0)
        END PROCEDURE vtk_legacy_write

        MODULE PROCEDURE vtk_legacy_read
        USE vtk_vars, ONLY : default_fn, default_title, filetype, vtkfilename, vtktitle, ascii, binary, version
        IMPLICIT NONE
        !! author: Ian Porter
        !! date: 12/20/2017
        !!
        !! This subroutines reads the legacy vtk output file
        !!
        INTEGER(i4k) :: i, inputstat
        LOGICAL      :: file_is_open
        CHARACTER(LEN=:), ALLOCATABLE :: form, filetype_text, vtk_version, line

        INQUIRE(unit = unit, opened = file_is_open)        !! Check to see if file is already open
        IF (.NOT. file_is_open) THEN                       !! File is not yet open. Determine format to open file
            OPEN(unit=unit, file=filename, iostat=inputstat, status='old')
                                                           !! Open the VTK file
        END IF

        READ(unit,100) vtk_version                         !! VTK version (currently, 3.0)
        IF (vtk_version /= version) THEN
                                                           !! If the file is not version 3.0, abort read
            ERROR STOP 'ERROR: vtkmofo only works with version 3.0'
        END IF

        READ(unit,100) title                               !! VTK title card
        READ(unit,100) filetype_text                       !! VTK file type

        CLOSE(unit)                                        !! Close the file to re-open it in the proper format

        SELECT CASE (filetype_text)
        CASE ('ASCII')
            data_type = ascii
            form      = 'formatted'
        CASE ('BINARY')
            data_type = binary
            form      = 'unformatted'
        END SELECT

        OPEN(unit=unit, file=filename, iostat=inputstat, status='old', form=form)
                                                           !! Open the VTK file in the proper format
        READ(unit,100) line                                !! Skip over this line
        READ(unit,100) line                                !! Skip over this line
        READ(unit,100) line                                !! Skip over this line

        CALL geometry%read(unit)                           !! Read the information from the file

        IF (PRESENT(celldatasets)) THEN
            READ(unit,*)
            DO i = 1, SIZE(celldatasets)
                CALL celldatasets(i)%attribute%read(unit)  !! Write the cell data values
            END DO
        ELSE IF (PRESENT(celldata)) THEN
            READ(unit,*)
            CALL celldata%write(unit)                      !! Write the cell data values
        END IF
        IF (PRESENT(pointdatasets)) THEN
            READ(unit,*)
            DO I = 1, SIZE(pointdatasets)
                CALL pointdatasets(i)%attribute%read(unit)
            END DO
        ELSE IF (PRESENT(pointdata)) THEN
            READ(unit,*)
            CALL pointdata%write(unit)                     !! Write the point data values
        END IF

        CLOSE(unit)                                        !! Close the VTK file

        vtkfilename = filename                             !! Save the filename for future internal use
        vtktitle    = title                                !! Save the title for future internal use
        filetype    = data_type                            !! Save the file type for future internal use

100     FORMAT(a)
        END PROCEDURE vtk_legacy_read

END SUBMODULE vtk_io_implementation


text:
Chef Debra Personal Chef Services
After a free inhome consultation, Chef Debra will create a master Recipes menu especially for you from which you can order each week.

As a personal chef, I look forward to cultivating relationships with my clients and making someone’s life a little easier and more delicious.
There are several packages available, depending on your needs.
Gift Certificates are available.
Cost for all services include:
Initial client interview  Menu planning  Grocery Shopping  Meal preparation   Packaging  Labeling and Cleanup.  Grocery costs are separate.

I offer special services such as:
 Vegetarian, Vagan and Raw Vegan diets, special diet needs, Pantry Makeovers, I can help you clean up your pantry and your eating habits for a healthier way of life, by teaching you which staples to keep on hand and how to find healthy substitutions for everyday ingredients.  You will feel confident making conscientious food choices and have an easier time creating your own meals that are satisfying in every way.
Dinner Parties, Special Events, Cooking classes and Cooking Parties.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
C
C
      SUBROUTINE DLAEIG(N, NBAND, NL, NR, A, EIGVAL, LDE,
     1   EIGVEC, BOUND, ATEMP, D, VTEMP, EPS, TMIN, TMAX)
C
C  THIS IS A SPECIALIZED VERSION OF THE SUBROUTINE BNDEIG TAILORED
C  SPECIFICALLY FOR USE BY THE LASO PACKAGE.
C
      INTEGER N, NBAND, NL, NR, LDE
      DOUBLE PRECISION A(NBAND,1), EIGVAL(1),
     1   EIGVEC(LDE,1), BOUND(2,1), ATEMP(1), D(1), VTEMP(1),
     2   EPS, TMIN, TMAX
C
C  LOCAL VARIABLES
C
      INTEGER I, M, NVAL
      DOUBLE PRECISION ARTOL, ATOL
C
C  FUNCTIONS CALLED
C
      DOUBLE PRECISION DMAX1
C
C  SUBROUTINES CALLED
C
C     DLABCM, DLABFC, DLAGER, DCOPY
C
C  SET PARAMETERS
C
      ATOL = DBLE(FLOAT(N))*EPS*DMAX1(TMAX,-TMIN)
      ARTOL = ATOL/DSQRT(EPS)
      NVAL = NR - NL + 1
C
C   CHECK FOR SPECIAL CASE OF N = 1
C
      IF(N .NE. 1) GO TO 30
      EIGVAL(1) = A(1,1)
      EIGVEC(1,1) = 1.0D0
      RETURN
C
C   SET UP INITIAL EIGENVALUE BOUNDS
C
   30 M = NVAL + 1
      DO 50 I = 2, M
         BOUND(1,I) = TMIN
         BOUND(2,I) = TMAX
   50 CONTINUE
      BOUND(2,1) = TMAX
      BOUND(1,NVAL + 2) = TMIN
      IF(NL .EQ. 1) BOUND(2,1) = TMIN
      IF(NR .EQ. N) BOUND(1,NVAL + 2) = TMAX
C
   60 CALL DLABCM(N, NBAND, NL, NR, A, EIGVAL, LDE,
     1  EIGVEC, ATOL, ARTOL, BOUND, ATEMP, D, VTEMP)
      RETURN
      END


text:
! { dg-do compile }
      SUBROUTINE A11_3(AA, BB, CC, DD, N)
      INTEGER N
      REAL AA(N,N), BB(N,N), CC(N,N), DD(N,N)
      REAL R
        R=0
!$OMP PARALLEL
!$OMP WORKSHARE
             AA = BB
!$OMP ATOMIC
               R = R + SUM(AA)
             CC = DD
!$OMP END WORKSHARE
!$OMP END PARALLEL
      END SUBROUTINE A11_3

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
Sunset is the occurrence of the sun receding beneath the western horizon. When viewed from Davis the sun disappears beneath the Coastal Range between the hours of 4 PM and 9 PM, depending on the time of year.

Sunset is also a great time for wildlife viewing, particularly birds and bats. For Bats, try looking beneath the Bike Overpasses pedestrian walkway over Covell, near Davis Senior High School, under the Poleline Pedestrian overpass, or under the Causeway. For birds check out anywhere there is water and or food as birds tend to congregate there at dusk and dawn.

 Places to watch the sunset
Freeway Bridges Noisy but expansive views.
Bermuda Avenue This is the end of residential Davis and has views to the Coastal Range.
Putah Creek Albeit lower than the surrounding areas, still offers nice reflections during sunset.
West Pond This spot is particularly nice during Winter when the pond is full of water and Town Fauna wildlife.
The top of Sproul Hall on the UCD campus offers a great silhouette of the Blue Ridge.

 Pictures of sunsets


text:
program main
    implicit none
    include 'omp_lib.h'
    integer :: i
    integer, parameter :: S = 2000
    double precision, pointer, dimension(:, :) :: a, b, c, d, e

    double precision t

    allocate(a(S, S), b(S, S), c(S, S), d(S, S), e(S, S))
    a = 0.0
    b = 0.0
    c = 0.0
    d = 0.0
    e = 0.0

    t = omp_get_wtime()
    do i = 1, 100
        call expr1(a, b)
    end do
    print *, "expr1", omp_get_wtime() - t

    t = omp_get_wtime()
    do i = 1, 100
        call expr2(a, b, c, d, e)
    end do
    print *, "expr2", omp_get_wtime() - t

    t = omp_get_wtime()
    do i = 1, 100
        call expr3(a, b, c, d, e)
    end do
    print *, "expr3", omp_get_wtime() - t

    print *, sum(a)

contains
    subroutine expr1(a, b)
        implicit none
        double precision, dimension(:, :) :: a, b
        a(:, :) = a(:, :) + transpose(b(:, :))
    end subroutine

    subroutine expr2(a, b, c, d, e)
        implicit none
        double precision, dimension(:, :) :: a, b, c, d, e
        a(:, :) = a + transpose(b) + c + transpose(d) + e
    end subroutine

    subroutine expr3(a, b, c, d, e)
        implicit none
        double precision, dimension(:, :) :: a, b, c, d, e
        a(:, :) = transpose(a) + b + transpose(c) + &
                  d + transpose(e)
    end subroutine
 end program

max_stars_repo_name: markflorisson/minivect

max_issues_repo_name: markflorisson/minivect

max_forks_repo_name: markflorisson/minivect


text:
	SUBROUTINE WP_OUT ( lunf, wtype, isutim, strtim, stptim, wnum,
     +			    nmpt, alat, alon, iret )
C************************************************************************
C* WP_OUT								*
C*									*
C* This subroutine writes the decoded watch corner point information	*
C* (watch type, start time, ending time, watch number, number of lat/lon*
C* pairs, and the lat/lon pairs) to an ASCII file.			*
C*									*
C* WO_OUT ( LUNF, ISUTIM, STRTIM, STPTIM, WNUM, NUMPT, ALAT, ALON, IRET)*
C*									*
C* Input parameters:							*
C*	LUNF		INTEGER		ASCII output file number	*
C*	WTYPE		CHAR*		Watch type (TN, TS, NA)		*
C*	ISUTIM		CHAR*		Issue time--GEMPAK format	*
C*	STRTIM		CHAR*		Watch start time--GEMPAK format	*
C*	STPTIM		CHAR*		Watch ending time--GEMPAK format*
C*      WNUM            CHAR*		Watch number			*
C*      NMPT            INTEGER		Number of lat/lon pairs decoded *
C*      ALAT (*)	REAL            Latitudes in degree-decimal     *
C*      ALON (*)	REAL            Longitude in degree-decimal	*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code 			*
C*					  0 = normal return		*
C*									*
C**									*
C* Log:									*
C* F. J. Yen/NCEP	 4/05						*
C* M. James/UCAR        05/16   gfortran 4.6+ fix for writing files     *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER*(*)	wtype, isutim, strtim, stptim, wnum
        REAL            alat (*), alon (*)
C*
	CHARACTER       string*55, anmpt*3
C-----------------------------------------------------------------------
	iret = 0
C
C*	Setting up the information string.
C
	CALL ST_LSTR ( wtype, lnwtyp, ier )
	CALL ST_LSTR ( isutim, lnissu, ier )
	CALL ST_LSTR ( strtim, lnstrt, ier )
	CALL ST_LSTR ( stptim, lnstpt, ier )
	CALL ST_LSTR ( wnum, lnwnum, ier )
	CALL ST_INCH ( nmpt, anmpt, ier )
	CALL ST_LSTR ( anmpt, lnnmpt, ier )

 	string = '|' // wtype(:lnwtyp)  // '|' // isutim(:lnissu) //
     +		 '|' // strtim(:lnstrt) // '|' // stptim(:lnstpt) //
     +		 '|' // wnum(:lnwnum)   // '|' // anmpt(:lnnmpt)
C
C*	Write the WCP information to the ASCII file.
C
	BACKSPACE(lunf)
	WRITE ( lunf, 20 ) string
 20     FORMAT ( A )
C
C*	Write the lat/lon pairs to the ASCII file.
C
	DO ii = 1, nmpt
	    WRITE ( lunf, 100 ) alat(ii), alon(ii)
 100	    FORMAT ( 2F9.2 )
	END DO
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
NotTires has left the building.  Not full graduated, I am also fully not in Davis.  I will rarely be checking this site, although might see your comment eventually.  I have moved on to another, distressingly smaller wikispot.

NotTires is a graduate students graduate student who will remain anonymous except to those who already suspect who he is!   This is partially because he might be your Teaching Assistants TA one day (and already was the TA for one major wikiuser that everyone has heard from).  He will eventually notice comments left here and respond, especially about Town Flora and Town Wildlife, the Climate Change environment, biodiesel biofuels, and Live Theater.  He is Judaism Jewish, in case you are curious, and happy to answer questions about religion, since talking is good.  Frequently seen in the outskirts driving around looking for plants for his research, he loves sharing his opinions, whether asked or not.   of course, everyone here does or we wouldnt edit a wiki, but I admit it up front.

p.s. NotTires is a pun on his name

Here is the comments section, in case anyone actually bothers to read this:

20080621 10:27:50 nbsp Hey... theres no template for departed businesses.  You just put {{{Include(Departed Business)}}} at the top of the entry.  Its just a notice to let people know when they are reading the entry.  Users/JabberWokky

20080826 01:53:28 nbsp Someone had actually posted that on my page at one point.  And it is true.  I am indeed not a pipe.  There is, however, a puzzle embedded into the items that Im listed as not being one of.  Figure out why I chose those particular items, and you solve the puzzle.  :) Users/IDoNotExist

20081113 05:41:33 nbsp Re: Your edit to Recycling: You dont have to refer to the other editor as commentor... you are both full editors and can write your positions as content, not as dialog. In other words, they stated your position, and you can state your position with equal footing (and give them the same footing).  Thus the entry starts to reflect the myriad of different viewpoints within the actual community of Davis. Users/JabberWokky

20081127 17:10:06 nbsp I will guess that your name is Jim. Users/ThanhVu


text:
This editor can click the Edit button to tell the community about him/herself.

20100609 21:47:16 nbsp Hi skurniaw, and welcome to the wiki!  You asked why no one had created a page for the Fox Squirrel Contraceptive Research Project.  All of the pages on the wiki are created by volunteer editors such as yourself.  Do you have information about the project?  If so, you can click on the dotted link in this comment and that will create a page.  Type in what you know. If you have trouble with the formatting, others can help out.  It would be great to have a page for that! Users/CovertProfessor


text:
      subroutine bmpfixed

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine applies fixed removal eff. from the .ops to upland loads

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name          |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ihru          |none          |HRU number
!!    bmp_sed(:)    |%             | Sediment removal by BMP
!!    bmp_pp(:)     |%             | Particulate P removal by BMP
!!    bmp_sp(:)     |%             | Soluble P removal by BMP
!!    bmp_pn(:)     |%             | Particulate N removal by BMP
!!    bmp_sn(:)     |%             | Soluble N removal by BMP
!!    bmp_bac(:)    |%             | Bacteria removal by BMP
!!    sedyld(:)     |metric tons   |daily soil loss caused by water erosion
!!    minpgw(:)     |kg P/ha       |soluble P loading to reach in groundwater
!!    sedminpa(:)   |kg P/ha       |amount of active mineral phosphorus sorbed
!!                                 |to sediment in surface runoff in HRU for day
!!    sedminps(:)   |kg P/ha       |amount of stable mineral phosphorus sorbed
!!                                 |to sediment in surface runoff in HRU for day
!!    sedorgp(:)    |kg P/ha       |amount of organic phosphorus in surface
!!                                 |runoff in HRU for the day
!!    surqsolp(:)   |kg P/ha       |amount of soluble phosphorus in surface
!!                                 |runoff in HRU for the day
!!    latno3(:)     |kg N/ha       |amount of NO3-N in lateral flow in HRU for
!!    surqno3(:)    |kg N/ha       |amount of NO3-N in surface runoff in HRU for
!!                                 |the day
!!    sedorgn(:)   |kg N/ha        |amount of organic nitrogen in surface runoff
!!    no3gw(:)     |kg N/ha        |nitrate loading to reach in groundwater
!!                                 |in HRU for the day
!!    bactrolp      |# cfu/m^2     |less persistent bacteria transported to main
!!                                 |channel with surface runoff
!!    bactrop       |# cfu/m^2     |persistent bacteria transported to main
!!                                 |channel with surface runoff
!!    bactsedlp     |# cfu/m^2     |less persistent bacteria transported with
!!                                 |sediment in surface runoff
!!    bactsedp      |# cfu/m^2     |persistent bacteria transported with
!!                                 |sediment in surface runoff

!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    sedyld(:)     |metric tons   |daily soil loss caused by water erosion
!!    minpgw(:)     |kg P/ha       |soluble P loading to reach in groundwater
!!    sedminpa(:)   |kg P/ha       |amount of active mineral phosphorus sorbed
!!                                 |to sediment in surface runoff in HRU for day
!!    sedminps(:)   |kg P/ha       |amount of stable mineral phosphorus sorbed
!!                                 |to sediment in surface runoff in HRU for day
!!    sedorgp(:)    |kg P/ha       |amount of organic phosphorus in surface
!!                                 |runoff in HRU for the day
!!    surqsolp(:)   |kg P/ha       |amount of soluble phosphorus in surface
!!                                 |runoff in HRU for the day
!!    latno3(:)     |kg N/ha       |amount of NO3-N in lateral flow in HRU for
!!    surqno3(:)    |kg N/ha       |amount of NO3-N in surface runoff in HRU for
!!                                 |the day
!!    sedorgn(:)   |kg N/ha        |amount of organic nitrogen in surface runoff
!!    no3gw(:)     |kg N/ha        |nitrate loading to reach in groundwater
!!                                 |in HRU for the day
!!    bactrolp      |# cfu/m^2     |less persistent bacteria transported to main
!!                                 |channel with surface runoff
!!    bactrop       |# cfu/m^2     |persistent bacteria transported to main
!!                                 |channel with surface runoff
!!    bactsedlp     |# cfu/m^2     |less persistent bacteria transported with
!!                                 |sediment in surface runoff
!!    bactsedp      |# cfu/m^2     |persistent bacteria transported with
!!                                 |sediment in surface runoff
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use parm

!!	set variables
      j = ihru

!! Subtract reductions from sediment, nutrients, bacteria, NOT SURFACE RUNOFF to protect water balance
!! Sediment
	sedyld(j) = sedyld(j) * (1. - bmp_sed(j)/100)

!! Phosphorus
      !! Particulate
	sedminpa(j) = sedminpa(j) * (1. - bmp_pp(j)/100.)
	sedminps(j) = sedminps(j) * (1. - bmp_pp(j)/100.)
	sedorgp(j) = sedorgp(j) * (1. - bmp_pp(j)/100.)
      !! Soluble
	surqsolp(j) = surqsolp(j) * (1. - bmp_sp(j)/100.)
	minpgw(j)= minpgw(j) * (1. - bmp_sp(j)/100.)

!! Nitrogen
	!! Particulate
	sedorgn(j) = sedorgn(j) * (1. - bmp_pn(j)/100.)
      !! Soluble
      surqno3(j) = surqno3(j) * (1. - bmp_sn(j)/100.)
	latno3(j) = latno3(j) * (1. - bmp_sn(j)/100.)
	no3gw(j) = no3gw(j) * (1. - bmp_sn(j)/100.)

!! Bacteria
      bactrop = bactrop * (1. - bmp_bac(j)/100.)
      bactrolp = bactrolp * (1. - bmp_bac(j)/100.)
      bactsedp = bactsedp * (1. - bmp_bac(j)/100.)
      bactsedlp = bactsedlp * (1. - bmp_bac(j)/100.)

      return
      end

max_stars_repo_name: GISWAT/erosion-sediment

max_issues_repo_name: GISWAT/erosion-sediment

max_forks_repo_name: GISWAT/erosion-sediment


text:
Okay, so were going to have a little ice breaker..Say your name, your major, and something funny about yourself!

Discussion sections are designed to counter the large lecture arrangement of most lowerdivision courses.  They are almost always taught by Teaching Assistants.


text:
MODULE kinetic_variables

  implicit none

  PRIVATE
  PUBLIC :: Md, coef_lap0, coef_lap, coef_nab &
           ,coef_nabk, const_k2, zcoef_kin, ggg &
           ,flag_nab, flag_n12, flag_n23, flag_n31
  PUBLIC :: wk
  PUBLIC :: SYStype
  PUBLIC :: read_kinetic, read_oldformat_kinetic
  PUBLIC :: coef_kin

  integer :: SYStype=0
  integer :: Md
  logical :: flag_nab, flag_n12, flag_n23, flag_n31
  real(8) :: coef_lap0, ggg(6)
  real(8),allocatable :: coef_lap(:,:), coef_nab(:,:)
  real(8),allocatable :: coef_nabk(:,:,:), const_k2(:)
  real(8),allocatable :: coef_kin(:)
  complex(8),allocatable :: zcoef_kin(:,:,:)
#ifdef _DRSDFT_
  real(8),allocatable :: wk(:,:,:,:)
#else
  complex(8),allocatable :: wk(:,:,:,:)
#endif

CONTAINS

  SUBROUTINE read_kinetic(rank,unit)
    implicit none
    integer,intent(IN) :: rank,unit
    integer :: i
    character(7) :: cbuf,ckey
    Md = 6
    SYStype = 0
    if ( rank == 0 ) then
       rewind unit
       do i=1,10000
          read(unit,*,END=999) cbuf
          call convert_capital(cbuf,ckey)
          if ( ckey(1:2) == "MD" ) then
             backspace(unit)
             read(unit,*) cbuf,Md
          else if ( ckey(1:7) == "SYSTYPE" ) then
             backspace(unit)
             read(unit,*) cbuf,SYStype
          end if
       end do
999    continue
       write(*,*) "Md =",Md
       write(*,*) "SYStype =",SYStype
    end if
    call send_kinetic(0)
  END SUBROUTINE read_kinetic

  SUBROUTINE read_oldformat_kinetic(rank,unit)
    implicit none
    integer,intent(IN) :: rank,unit
    if ( rank == 0 ) then
       read(unit,*) Md, SYStype
       write(*,*) "Md =",Md
       write(*,*) "SYStype =",SYStype
    end if
    call send_kinetic(0)
  END SUBROUTINE read_oldformat_kinetic

  SUBROUTINE send_kinetic(rank)
    implicit none
    integer,intent(IN) :: rank
    integer :: ierr
    include 'mpif.h'
    call mpi_bcast(Md,1,MPI_INTEGER,rank,MPI_COMM_WORLD,ierr)
    call mpi_bcast(SYStype,1,MPI_INTEGER,rank,MPI_COMM_WORLD,ierr)
  END SUBROUTINE send_kinetic

END MODULE kinetic_variables

max_stars_repo_name: j-iwata/RSDFT

max_issues_repo_name: j-iwata/RSDFT

max_forks_repo_name: j-iwata/RSDFT


text:
! -*- F90 -*-

      subroutine GRVGevolvep0(xin,qin,p2in,ip2in,pdf)
      include 'parmsetup.inc'
      real*8 xin,qin,q2in,p2in,pdf(-6:6),xval(45),qcdl4,qcdl5
      real*8 upv,dnv,usea,dsea,str,chm,bot,top,glu,zbot,zchm
      character*16 name(nmxset)
      integer nmem(nmxset),ndef(nmxset),mmem
      common/NAME/name,nmem,ndef,mmem
      integer nset

      save

      call getnset(iset)
      call getnmem(iset,imem)
      if(imem.eq.1) then
        call GRVGALO (xin,qin,upv,dnv,usea,dsea,str,chm,bot,glu)
      elseif(imem.eq.2.or.imem.eq.0) then
        q2in = qin*qin
! calls GRVGALO for charm and bottom, rest from GRSGALO
        call GRVGALO(xin,qin,upv,dnv,usea,dsea,str,chm,bot,glu)
        call GRSGALO(xin,q2in,p2in,                                     &
     &               upv,dnv,usea,dsea,str,zchm,zbot,glu)
      else
        CONTINUE
      endif
      pdf(-6)= 0.0d0
      pdf(6)= 0.0d0
      pdf(-5)= bot
      pdf(5 )= bot
      pdf(-4)= chm
      pdf(4 )= chm
      pdf(-3)= str
      pdf(3 )= str
      pdf(-2)= usea
      pdf(2 )= upv
      pdf(-1)= dsea
      pdf(1 )= dnv
      pdf(0 )= glu

      return
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry GRVGevolvep1(xin,qin,p2in,ip2in,pdf)

      if(imem.eq.1) then
        call GRVGAH0 (xin,qin,upv,dnv,usea,dsea,str,chm,bot,glu)
      elseif(imem.eq.2 .or. imem.eq.0) then
        call GRVGAHO (xin,qin,upv,dnv,usea,dsea,str,chm,bot,glu)
      else
        CONTINUE
      endif

      pdf(-6)= 0.0d0
      pdf(6)= 0.0d0
      pdf(-5)= bot
      pdf(5 )= bot
      pdf(-4)= chm
      pdf(4 )= chm
      pdf(-3)= str
      pdf(3 )= str
      pdf(-2)= usea
      pdf(2 )= upv
      pdf(-1)= dsea
      pdf(1 )= dnv
      pdf(0 )= glu

      return
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry GRVGread(nset)
      read(1,*)nmem(nset),ndef(nset)
      return
!
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry GRVGalfa(alfas,qalfa)
        call getnset(iset)
        call GetOrderAsM(iset,iord)
        call Getlam4M(iset,imem,qcdl4)
        call Getlam5M(iset,imem,qcdl5)
        call aspdflib(alfas,Qalfa,iord,qcdl5)
      return
!
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry GRVGinit(Eorder,Q2fit)
      return
!
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      entry GRVGpdf(mem)
      call getnset(iset)
      call setnmem(iset,mem)

!      imem = mem
      return
!
 1000 format(5e13.5)
      END
!
       SUBROUTINE GRVGAH0 (ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZBB,ZGL)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!                                                                 *
!      G R V - P H O T O N - P A R A M E T R I Z A T I O N S      *
!                                                                 *
!                 FOR A DETAILED EXPLANATION SEE :                *
!              M. GLUECK, E.REYA, A.VOGT: DO-TH 91/31             *
!                                                                 *
!    THE OUTPUT IS ALWAYS   1./ ALPHA(EM) * X * PARTON DENSITY    *
!    output modified by HPB to be always    X * PARTON DENSITY    *
!                                                                 *
!   THE PARAMETRIZATIONS ARE FITTED TO THE PARTON DISTRIBUTIONS   *
!   FOR Q ** 2 BETWEEN MU ** 2 (=  0.25 / 0.30  GEV ** 2  IN LO   *
!   / HO) AND  1.E6 GEV ** 2  AND FOR X BETWEEN  1.E-5  AND  1.   *
!                                                                 *
!              HEAVY QUARK THRESHOLDS  Q(H) = M(H) :              *
!         M(C)  =  1.5,  M(B)  =  4.5,  M(T)  =  100  GEV         *
!                                                                 *
!      CORRESPONDING LAMBDA(F) VALUES FOR F ACTIVE FLAVOURS :     *
!      LO :   LAMBDA(3)  =  0.232,   LAMBDA(4)  =  0.200,         *
!             LAMBDA(5)  =  0.153,   LAMBDA(6)  =  0.082  GEV     *
!      HO :   LAMBDA(3)  =  0.248,   LAMBDA(4)  =  0.200,         *
!             LAMBDA(5)  =  0.131,   LAMBDA(6)  =  0.053  GEV     *
!                                                                 *
!      HO DISTRIBUTIONS REFER TO THE DIS(GAMMA) SCHEME, SEE :     *
!              M. GLUECK, E.REYA, A.VOGT: DO-TH 91/26             *
!                                                                 *
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
       IMPLICIT REAL (A - Y)
       double precision                                                 &
     &        ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZBB,ZGL
       REAL  X, Q
       DATA ALPHEM/7.29927D-3/
       X = ZX
       Q = ZQ
       MU2  = 0.3
       LAM2 = 0.248 * 0.248
       Q2 = Q*Q
       S  = ALOG (ALOG(Q2/LAM2) / ALOG(MU2/LAM2))
       SS = SQRT (S)
       S2 = S * S
!...X * U = X * UBAR :
       AL =  1.447
       BE =  0.848
       AK =  0.527 + 0.200 * S  - 0.107 * S2
       BK =  7.106 - 0.310 * SS - 0.786 * S2
       AG =  0.197 + 0.533 * S
       BG =  0.062 - 0.398 * S  + 0.109 * S2
       C  =          0.755 * S  - 0.112 * S2
       D  =  0.318 - 0.059 * S
       E  =  4.225 + 1.708 * S
       ES =  1.752 + 0.866 * S
       U0 =  GRVGF (X, S, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZUV = U0 * ALPHEM
       ZUB = ZUV
!...X * D = X * DBAR :
       AL =  1.424
       BE =  0.770
       AK =  0.500 + 0.067 * SS - 0.055 * S2
       BK =  0.376 - 0.453 * SS + 0.405 * S2
       AG =  0.156 + 0.184 * S
       BG =   0.0  - 0.528 * S  + 0.146 * S2
       C  =  0.121 + 0.092 * S
       D  =  0.379 - 0.301 * S  + 0.081 * S2
       E  =  4.346 + 1.638 * S
       ES =  1.645 + 1.016 * S
       D0  =  GRVGF (X, S, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZDV = D0 * ALPHEM
       ZDB = ZDV
!...X * G :
       AL =  0.661
       BE =  0.793
       AK =  0.537 - 0.600 * SS
       BK =  6.389              - 0.953 * S2
       AG =  0.558 - 0.383 * SS + 0.261 * S2
       BG =   0.0  - 0.305 * S
       C  = -0.222              + 0.078 * S2
       D  =  0.153 + 0.978 * S  - 0.209 * S2
       E  =  1.429 + 1.772 * S
       ES =  3.331 + 0.806 * S
       G0 =  GRVGF (X, S, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZGL = G0 * ALPHEM
!...X * S = X * SBAR :
       SF =   0.0
       AL =  1.578
       BE =  0.863
       AK =  0.622 + 0.332 * S  - 0.300 * S2
       BK =  2.469
       AG =  0.211 - 0.064 * SS - 0.018 * S2
       BG = -0.215 + 0.122 * S
       C  =  0.153
       D  =   0.0  + 0.253 * S  - 0.081 * S2
       E  =  3.990 + 2.014 * S
       ES =  1.720 + 0.986 * S
       S0 =  GRVGFS (X, S, SF, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZSB = S0 * ALPHEM
!...X * C = X * CBAR :
       SF =  0.820
       AL =  0.929
       BE =  0.381
       AK =  1.228 - 0.231 * S
       BK =  3.806             - 0.337 * S2
       AG =  0.932 + 0.150 * S
       BG = -0.906
       C  =  1.133
       D  =   0.0  + 0.138 * S  - 0.028 * S2
       E  =  5.588 + 0.628 * S
       ES =  2.665 + 1.054 * S
       C0 =  GRVGFS (X, S, SF, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZCB = C0 * ALPHEM
!...X * B = X * BBAR :
       SF =  1.297
       AL =  0.970
       BE =  0.207
       AK =  1.719 - 0.292 * S
       BK =  0.928 + 0.096 * S
       AG =  0.845 + 0.178 * S
       BG = -2.310
       C  =  1.558
       D  = -0.191 + 0.151 * S
       E  =  6.089 + 0.282 * S
       ES =  3.379 + 1.062 * S
       B0 =  GRVGFS (X, S, SF, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZBB = B0 * ALPHEM
!
       RETURN
      END
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       SUBROUTINE GRVGAHO (ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZBB,ZGL)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!                                                                 *
!      G R V - P H O T O N - P A R A M E T R I Z A T I O N S      *
!                                                                 *
!                 FOR A DETAILED EXPLANATION SEE :                *
!              M. GLUECK, E.REYA, A.VOGT: DO-TH 91/31             *
!                                                                 *
!    THE OUTPUT IS ALWAYS   1./ ALPHA(EM) * X * PARTON DENSITY    *
!    output modified by HPB to be always    X * PARTON DENSITY    *
!                                                                 *
!   THE PARAMETRIZATIONS ARE FITTED TO THE PARTON DISTRIBUTIONS   *
!   FOR Q ** 2 BETWEEN MU ** 2 (=  0.25 / 0.30  GEV ** 2  IN LO   *
!   / HO) AND  1.E6 GEV ** 2  AND FOR X BETWEEN  1.E-5  AND  1.   *
!                                                                 *
!              HEAVY QUARK THRESHOLDS  Q(H) = M(H) :              *
!         M(C)  =  1.5,  M(B)  =  4.5,  M(T)  =  100  GEV         *
!                                                                 *
!      CORRESPONDING LAMBDA(F) VALUES FOR F ACTIVE FLAVOURS :     *
!      LO :   LAMBDA(3)  =  0.232,   LAMBDA(4)  =  0.200,         *
!             LAMBDA(5)  =  0.153,   LAMBDA(6)  =  0.082  GEV     *
!      HO :   LAMBDA(3)  =  0.248,   LAMBDA(4)  =  0.200,         *
!             LAMBDA(5)  =  0.131,   LAMBDA(6)  =  0.053  GEV     *
!                                                                 *
!      HO DISTRIBUTIONS REFER TO THE DIS(GAMMA) SCHEME, SEE :     *
!              M. GLUECK, E.REYA, A.VOGT: DO-TH 91/26             *
!                                                                 *
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
       IMPLICIT REAL (A - Y)
       double precision                                                 &
     &        ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZBB,ZGL
       DATA ALPHEM/7.29927D-3/
       REAL  X, Q
       X = ZX
       Q = ZQ
       MU2  = 0.3
       LAM2 = 0.248 * 0.248
       Q2 = Q*Q
       S  = ALOG (ALOG(Q2/LAM2) / ALOG(MU2/LAM2))
       SS = SQRT (S)
       S2 = S * S
!...X * U = X * UBAR :
       AL =  0.583
       BE =  0.688
       AK =  0.449 - 0.025 * S  - 0.071 * S2
       BK =  5.060 - 1.116 * SS
       AG =  0.103
       BG =  0.319 + 0.422 * S
       C  =  1.508 + 4.792 * S  - 1.963 * S2
       D  =  1.075 + 0.222 * SS - 0.193 * S2
       E  =  4.147 + 1.131 * S
       ES =  1.661 + 0.874 * S
       UH =  GRVGF (X, S, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZUV = UH * ALPHEM
       ZUB = ZUV
!...X * D = X * DBAR :
       AL =  0.591
       BE =  0.698
       AK =  0.442 - 0.132 * S  - 0.058 * S2
       BK =  5.437 - 1.916 * SS
       AG =  0.099
       BG =  0.311 - 0.059 * S
       C  =  0.800 + 0.078 * S  - 0.100 * S2
       D  =  0.862 + 0.294 * SS - 0.184 * S2
       E  =  4.202 + 1.352 * S
       ES =  1.841 + 0.990 * S
       DH  =  GRVGF (X, S, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZDV = DH * ALPHEM
       ZDB = ZDV
!...X * G :
       AL =  1.161
       BE =  1.591
       AK =  0.530 - 0.742 * SS + 0.025 * S2
       BK =  5.662
       AG =  0.533 - 0.281 * SS + 0.218 * S2
       BG =  0.025 - 0.518 * S  + 0.156 * S2
       C  = -0.282              + 0.209 * S2
       D  =  0.107 + 1.058 * S  - 0.218 * S2
       E  =   0.0  + 2.704 * S
       ES =  3.071 - 0.378 * S
       GH =  GRVGF (X, S, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZGL = GH * ALPHEM
!...X * S = X * SBAR :
       SF =   0.0
       AL =  0.635
       BE =  0.456
       AK =  1.770 - 0.735 * SS - 0.079 * S2
       BK =  3.832
       AG =  0.084 - 0.023 * S
       BG =  0.136
       C  =  2.119 - 0.942 * S  + 0.063 * S2
       D  =  1.271 + 0.076 * S  - 0.190 * S2
       E  =  4.604 + 0.737 * S
       ES =  1.641 + 0.976 * S
       SH =  GRVGFS (X, S, SF, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZSB = SH * ALPHEM
!...X * C = X * CBAR :
       SF =  0.820
       AL =  0.926
       BE =  0.152
       AK =  1.142 - 0.175 * S
       BK =  3.276
       AG =  0.504 + 0.317 * S
       BG = -0.433
       C  =  3.334
       D  =  0.398 + 0.326 * S  - 0.107 * S2
       E  =  5.493 + 0.408 * S
       ES =  2.426 + 1.277 * S
       CH =  GRVGFS (X, S, SF, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZCB = CH * ALPHEM
!...X * B = X * BBAR :
       SF =  1.297
       AL =  0.969
       BE =  0.266
       AK =  1.953 - 0.391 * S
       BK =  1.657 - 0.161 * S
       AG =  1.076 + 0.034 * S
       BG = -2.015
       C  =  1.662
       D  =  0.353 + 0.016 * S
       E  =  5.713 + 0.249 * S
       ES =  3.456 + 0.673 * S
       BH =  GRVGFS (X, S, SF, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZBB = BH * ALPHEM
!
       RETURN
      END
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       SUBROUTINE GRVGALO (ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZBB,ZGL)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!                                                                 *
!      G R V - P H O T O N - P A R A M E T R I Z A T I O N S      *
!                                                                 *
!                 FOR A DETAILED EXPLANATION SEE :                *
!              M. GLUECK, E.REYA, A.VOGT: DO-TH 91/31             *
!                                                                 *
!    THE OUTPUT IS ALWAYS   1./ ALPHA(EM) * X * PARTON DENSITY    *
!    output modified by HPB to be always    X * PARTON DENSITY    *
!                                                                 *
!   THE PARAMETRIZATIONS ARE FITTED TO THE PARTON DISTRIBUTIONS   *
!   FOR Q ** 2 BETWEEN MU ** 2 (=  0.25 / 0.30  GEV ** 2  IN LO   *
!   / HO) AND  1.E6 GEV ** 2  AND FOR X BETWEEN  1.E-5  AND  1.   *
!                                                                 *
!              HEAVY QUARK THRESHOLDS  Q(H) = M(H) :              *
!         M(C)  =  1.5,  M(B)  =  4.5,  M(T)  =  100  GEV         *
!                                                                 *
!      CORRESPONDING LAMBDA(F) VALUES FOR F ACTIVE FLAVOURS :     *
!      LO :   LAMBDA(3)  =  0.232,   LAMBDA(4)  =  0.200,         *
!             LAMBDA(5)  =  0.153,   LAMBDA(6)  =  0.082  GEV     *
!      HO :   LAMBDA(3)  =  0.248,   LAMBDA(4)  =  0.200,         *
!             LAMBDA(5)  =  0.131,   LAMBDA(6)  =  0.053  GEV     *
!                                                                 *
!      HO DISTRIBUTIONS REFER TO THE DIS(GAMMA) SCHEME, SEE :     *
!              M. GLUECK, E.REYA, A.VOGT: DO-TH 91/26             *
!                                                                 *
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
       IMPLICIT REAL (A - Y)
       double precision                                                 &
     &        ZX,ZQ,ZUV,ZDV,ZUB,ZDB,ZSB,ZCB,ZBB,ZGL
       REAL  X, Q
       DATA ALPHEM/7.29927D-3/
       X = ZX
       Q = ZQ
       MU2  = 0.25
       LAM2 = 0.232 * 0.232
       Q2 = Q*Q
       S  = ALOG (ALOG(Q2/LAM2) / ALOG(MU2/LAM2))
       SS = SQRT (S)
       S2 = S * S
!...X * U = X * UBAR :
       AL =  1.717
       BE =  0.641
       AK =  0.500 - 0.176 * S
       BK = 15.00  - 5.687 * SS - 0.552 * S2
       AG =  0.235 + 0.046 * SS
       BG =  0.082 - 0.051 * S  + 0.168 * S2
       C  =   0.0  + 0.459 * S
       D  =  0.354 - 0.061 * S
       E  =  4.899 + 1.678 * S
       ES =  2.046 + 1.389 * S
       UL =  GRVGF (X, S, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZUV = UL * ALPHEM
       ZUB = ZUV
!...X * D = X * DBAR :
       AL =  1.549
       BE =  0.782
       AK =  0.496 + 0.026 * S
       BK =  0.685 - 0.580 * SS + 0.608 * S2
       AG =  0.233 + 0.302 * S
       BG =   0.0  - 0.818 * S  + 0.198 * S2
       C  =  0.114 + 0.154 * S
       D  =  0.405 - 0.195 * S  + 0.046 * S2
       E  =  4.807 + 1.226 * S
       ES =  2.166 + 0.664 * S
       DL  =  GRVGF (X, S, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZDV = DL * ALPHEM
       ZDB = ZDV
!...X * G :
       AL =  0.676
       BE =  1.089
       AK =  0.462 - 0.524 * SS
       BK =  5.451              - 0.804 * S2
       AG =  0.535 - 0.504 * SS + 0.288 * S2
       BG =  0.364 - 0.520 * S
       C  = -0.323              + 0.115 * S2
       D  =  0.233 + 0.790 * S  - 0.139 * S2
       E  =  0.893 + 1.968 * S
       ES =  3.432 + 0.392 * S
       GL =  GRVGF (X, S, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZGL = GL * ALPHEM
!...X * S = X * SBAR :
       SF =   0.0
       AL =  1.609
       BE =  0.962
       AK =  0.470              - 0.099 * S2
       BK =  3.246
       AG =  0.121 - 0.068 * SS
       BG = -0.090 + 0.074 * S
       C  =  0.062 + 0.034 * S
       D  =   0.0  + 0.226 * S  - 0.060 * S2
       E  =  4.288 + 1.707 * S
       ES =  2.122 + 0.656 * S
       SL =  GRVGFS (X, S, SF, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZSB = SL * ALPHEM
!...X * C = X * CBAR :
       SF =  0.888
       AL =  0.970
       BE =  0.545
       AK =  1.254 - 0.251 * S
       BK =  3.932              - 0.327 * S2
       AG =  0.658 + 0.202 * S
       BG = -0.699
       C  =  0.965
       D  =   0.0  + 0.141 * S  - 0.027 * S2
       E  =  4.911 + 0.969 * S
       ES =  2.796 + 0.952 * S
       CL =  GRVGFS (X, S, SF, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZCB = CL * ALPHEM
!...X * B = X * BBAR :
       SF =  1.351
       AL =  1.016
       BE =  0.338
       AK =  1.961 - 0.370 * S
       BK =  0.923 + 0.119 * S
       AG =  0.815 + 0.207 * S
       BG = -2.275
       C  =  1.480
       D  = -0.223 + 0.173 * S
       E  =  5.426 + 0.623 * S
       ES =  3.819 + 0.901 * S
       BL =  GRVGFS (X, S, SF, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       ZBB = BL * ALPHEM
!
       RETURN
      END
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!-----------------------------------------------------------------------
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!                                                                 *
!           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.19
!
!*******************************************************
      SUBROUTINE GRSGALO(DX,DQ2,DP2,                                    &
     &                     DUPV,DDNV,DUSEA,DDSEA,DSTR,DCHM,DBOT,DGL)
!      subroutine grsgalo(x,q2,p2,ugam,dgam,sgam,ggam)
      implicit real*8 (a-h,o-z)
      double precision                                                  &
     &          x, q2, p2, mu2, lam2,                                   &
     &          ugam, dgam, sgam, ggam,                                 &
     &          DUPV,DDNV,DUSEA,DDSEA,DSTR,DCHM,DBOT,DGL
!
      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/
!
      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/
!
      save u1,ds1,g1,ud2,s2,g2,up0,dsp0,gp0
!
      x  = DX
      q  = SQRT(DQ2)
      q2 = DQ2
      p2 = DP2
      mu2=0.25d0
      lam2=0.232d0*0.232d0
!
      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
!
      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)
!
      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)
!
      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)
!
      s=log(log(q2/lam2)/log(mu2/lam2))
      suppr=1.d0/(1.d0+p2/0.59d0)**2
!
      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)
!
      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)
!
      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)
!
      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
!
      DCHM = 0.D0
      DBOT = 0.D0
!
      return
      END
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       FUNCTION GRVGF (X, S, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       IMPLICIT REAL (A - Z)
       SX = SQRT (X)
       LX = ALOG (1./X)
       GRVGF  = (X**AK * (AG + BG * SX + C * X**BK)  +  S**AL           &
     &       * EXP (-E + SQRT (ES * S**BE * LX))) * (1.- X)**D
       RETURN
      END
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
       FUNCTION GRVGFS (X, S, SF, AL, BE, AK, BK, AG, BG, C, D, E, ES)
       IMPLICIT REAL (A - Z)
       IF (S .LE. SF) THEN
          GRVGFS = 0.0
       ELSE
          SX = SQRT (X)
          LX = ALOG (1./X)
          DS = S - SF
          GRVGFS = (DS * X**AK * (AG + BG * SX + C * X**BK) + DS**AL    &
     &         * EXP (-E + SQRT (ES * S**BE * LX))) * (1.- X)**D
       END IF
       RETURN
      END
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      double precision function grsf1(x,s,alp,bet,a,b,ga,gb,gc,gd,      &
     &                                ge,gep)
      implicit real*8 (a-h,o-z)
!
      grsf1=(x**a*(ga+gb*sqrt(x)+gc*x**b)+                              &
     &      s**alp*exp(-ge+sqrt(gep*s**bet*log(1.d0/x))))*              &
     &      (1.d0-x)**gd
      return
      END
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      double precision function grsf2(x,s,alp,bet,a,b,ga,gb,gc,gd,      &
     &                                ge,gep)
      implicit real*8 (a-h,o-z)
!
      grsf2=(s*x**a*(ga+gb*sqrt(x)+gc*x**b)+                            &
     &      s**alp*exp(-ge+sqrt(gep*s**bet*log(1.d0/x))))*              &
     &      (1.d0-x)**gd
      return
      END

max_stars_repo_name: AllaMaevskaya/AliRoot

max_issues_repo_name: AllaMaevskaya/AliRoot

max_forks_repo_name: AllaMaevskaya/AliRoot