 
C*********************************************************************
 
C...PYRGHM
C...Auxiliary routine to PYVACU for SUSY Higgs calculations.
 
      SUBROUTINE PYRGHM(XMC,XMA,TANB,XMQ,XMUR,XMDL,XMT,AU,AD,XMU,
     &XMHP,HMP,SA,CA,TANBA)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(A-H, O-Z)
      IMPLICIT INTEGER(I-N)
      INTEGER PYK,PYCHGE,PYCOMP
      COMMON/PYHTRI/HHH(7)
 
C...Local variables.
      DIMENSION VH(2,2),XM2(2,2),XM2P(2,2)
 
      XMZ = 91.18D0
      ALP1 = 0.0101D0
      ALP2 = 0.0337D0
      ALP3Z = 0.12D0
      V = 174.1D0
      PI = 3.14159D0
      TANBA = TANB
      TANBT = TANB
 
C...MBOTTOM(XMT) = 3. GEV
      XMB = 3D0
      ALP3 = ALP3Z/(1D0 +(11D0 - 10D0/3D0)/4D0/PI*ALP3Z*
     &LOG(XMT**2/XMZ**2))
 
C...RXMT= RUNNING TOP QUARK MASS
      RXMT = XMT/(1D0+4D0*ALP3/3D0/PI)
      TQ = LOG((XMQ**2+XMT**2)/XMT**2)
      TU = LOG((XMUR**2 + XMT**2)/XMT**2)
      TD = LOG((XMDL**2 + XMT**2)/XMT**2)
      SINB = TANB/((1D0 + TANB**2)**0.5D0)
      COSB = SINB/TANB
      IF(XMA.GT.XMT)
     &TANBA = TANB*(1D0-3D0/32D0/PI**2*
     &(RXMT**2/V**2/SINB**2-XMB**2/V**2/COSB**2)*
     &LOG(XMA**2/XMT**2))
      IF(XMA.LT.XMT.OR.XMA.EQ.XMT) TANBT = TANBA
      SINB = TANBT/((1D0 + TANBT**2)**0.5D0)
      COSB = 1D0/((1D0 + TANBT**2)**0.5D0)
      COS2B = (TANBT**2 - 1D0)/(TANBT**2 + 1D0)
      G1 = (ALP1*4D0*PI)**0.5D0
      G2 = (ALP2*4D0*PI)**0.5D0
      G3 = (ALP3*4D0*PI)**0.5D0
      HU = RXMT/V/SINB
      HD =  XMB/V/COSB
 
      CALL PYGFXX(XMA,TANBA,XMQ,XMUR,XMDL,XMT,AU,AD,
     &XMU,VH,STOP1,STOP2)
 
      IF(XMQ.GT.XMUR) TP = TQ - TU
      IF(XMQ.LT.XMUR.OR.XMQ.EQ.XMUR) TP = TU - TQ
      IF(XMQ.GT.XMUR) TDP = TU
      IF(XMQ.LT.XMUR.OR.XMQ.EQ.XMUR) TDP = TQ
      IF(XMQ.GT.XMDL) TPD = TQ - TD
      IF(XMQ.LT.XMDL.OR.XMQ.EQ.XMDL) TPD = TD - TQ
      IF(XMQ.GT.XMDL) TDPD = TD
      IF(XMQ.LT.XMDL.OR.XMQ.EQ.XMDL) TDPD = TQ
 
      IF(XMQ.GT.XMDL) DLAM1 = 6D0/96D0/PI**2*G1**2*HD**2*TPD
      IF(XMQ.LT.XMDL.OR.XMQ.EQ.XMDL) DLAM1 = 3D0/32D0/PI**2*
     &HD**2*(G1**2/3D0+G2**2)*TPD
 
      IF(XMQ.GT.XMUR) DLAM2 =12D0/96D0/PI**2*G1**2*HU**2*TP
      IF(XMQ.LT.XMUR.OR.XMQ.EQ.XMUR) DLAM2 = 3D0/32D0/PI**2*
     &HU**2*(-G1**2/3D0+G2**2)*TP
 
      DLAM3 = 0D0
      DLAM4 = 0D0
 
      IF(XMQ.GT.XMDL) DLAM3 = -1D0/32D0/PI**2*G1**2*HD**2*TPD
      IF(XMQ.LT.XMDL.OR.XMQ.EQ.XMDL) DLAM3 = 3D0/64D0/PI**2*HD**2*
     &(G2**2-G1**2/3D0)*TPD
 
      IF(XMQ.GT.XMUR) DLAM3 = DLAM3 -
     &1D0/16D0/PI**2*G1**2*HU**2*TP
      IF(XMQ.LT.XMUR.OR.XMQ.EQ.XMUR) DLAM3 = DLAM3 +
     &3D0/64D0/PI**2*HU**2*(G2**2+G1**2/3D0)*TP
 
      IF(XMQ.LT.XMUR) DLAM4 = -3D0/32D0/PI**2*G2**2*HU**2*TP
      IF(XMQ.LT.XMDL) DLAM4 = DLAM4 - 3D0/32D0/PI**2*G2**2*
     &HD**2*TPD
 
      XLAM1 = ((G1**2 + G2**2)/4D0)*
     &(1D0-3D0*HD**2*(TPD + TDPD)/8D0/PI**2)
     &+(3D0*HD**4/16D0/PI**2) *TPD*(1D0
     &+ (3D0*HD**2/2D0 + HU**2/2D0
     &- 8D0*G3**2) * (TPD + 2D0*TDPD)/16D0/PI**2)
     &+(3D0*HD**4/8D0/PI**2) *TDPD*(1D0  + (3D0*HD**2/2D0 + HU**2/2D0
     &- 8D0*G3**2) * TDPD/16D0/PI**2) + DLAM1
      XLAM2 = ((G1**2 + G2**2)/4D0)*(1D0-3D0*HU**2*
     &(TP + TDP)/8D0/PI**2)
     &+(3D0*HU**4/16D0/PI**2) *TP*(1D0
     &+ (3D0*HU**2/2D0 + HD**2/2D0
     &- 8D0*G3**2) * (TP + 2D0*TDP)/16D0/PI**2)
     &+(3D0*HU**4/8D0/PI**2) *TDP*(1D0 + (3D0*HU**2/2D0 + HD**2/2D0
     &- 8D0*G3**2) * TDP/16D0/PI**2) + DLAM2
      XLAM3 = ((G2**2 - G1**2)/4D0)*(1D0-3D0*
     &(HU**2)*(TP + TDP)/16D0/PI**2 -3D0*
     &(HD**2)*(TPD + TDPD)/16D0/PI**2) +DLAM3
      XLAM4 = (- G2**2/2D0)*(1D0
     &-3D0*(HU**2)*(TP + TDP)/16D0/PI**2
     &-3D0*(HD**2)*(TPD + TDPD)/16D0/PI**2) +DLAM4
 
      XLAM5 = 0D0
      XLAM6 = 0D0
      XLAM7 = 0D0
 
C...Defined now in PYSUBH
C      HHH(1)=XLAM1
C      HHH(2)=XLAM2
C      HHH(3)=XLAM3
C      HHH(4)=XLAM4
C      HHH(5)=XLAM5
C      HHH(6)=XLAM6
C      HHH(7)=XLAM7
 
      XM2(1,1) = 2D0*V**2*(XLAM1*COSB**2+2D0*XLAM6*
     &COSB*SINB + XLAM5*SINB**2) + XMA**2*SINB**2
 
      XM2(2,2) = 2D0*V**2*(XLAM5*COSB**2+2D0*XLAM7*
     &COSB*SINB + XLAM2*SINB**2) + XMA**2*COSB**2
      XM2(1,2) = 2D0*V**2*(XLAM6*COSB**2+(XLAM3+XLAM4)*
     &COSB*SINB + XLAM7*SINB**2) - XMA**2*SINB*COSB
 
      XM2(2,1) = XM2(1,2)
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C...THIS IS THE CONTRIBUTION FROM LIGHT CHARGINOS/NEUTRALINOS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
      XMSSU=(0.5D0*(XMQ**2+XMUR**2)+XMT**2)**0.5D0
 
      IF(XMC.GT.XMSSU) GOTO 100
      IF(XMC.LT.XMT) XMC=XMT
 
      TCHAR=LOG(XMSSU**2/XMC**2)
 
      DEL12=(9D0/64D0/PI**2*G2**4+5D0/192D0/PI**2*G1**4)*TCHAR
      DEL3P4=(3D0/64D0/PI**2*G2**4+7D0/192D0/PI**2*G1**4
     &+4D0/32/PI**2*G1**2*G2**2)*TCHAR
 
      DEM112=2D0*DEL12*V**2*COSB**2
      DEM222=2D0*DEL12*V**2*SINB**2
      DEM122=2D0*DEL3P4*V**2*SINB*COSB
 
      XM2(1,1)=XM2(1,1)+DEM112
      XM2(2,2)=XM2(2,2)+DEM222
      XM2(1,2)=XM2(1,2)+DEM122
      XM2(2,1)=XM2(2,1)+DEM122
 
  100 CONTINUE
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C...END OF CHARGINOS/NEUTRALINOS
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
 
      DO 120 I = 1,2
        DO 110 J = 1,2
          XM2P(I,J) = XM2(I,J) + VH(I,J)
  110   CONTINUE
  120 CONTINUE
 
      TRM2P = XM2P(1,1) + XM2P(2,2)
      DETM2P = XM2P(1,1)*XM2P(2,2) - XM2P(1,2)*XM2P(2,1)
 
      XMH2P = (TRM2P - (TRM2P**2 - 4D0* DETM2P)**0.5D0)/2D0
      HM2P = (TRM2P + (TRM2P**2 - 4D0* DETM2P)**0.5D0)/2D0
      HMP = HM2P**0.5D0
      IF(XMH2P.LT.0D0) GOTO 130
      XMHP = XMH2P**0.5D0
      S2ALP = 2D0*XM2P(1,2)/(TRM2P**2-4D0*DETM2P)**0.5D0
      C2ALP = (XM2P(1,1)-XM2P(2,2))/(TRM2P**2-4D0*DETM2P)**0.5D0
      IF(C2ALP.GT.0D0) ALP = ASIN(S2ALP)/2D0
      IF(C2ALP.LT.0D0) ALP = -PI/2D0-ASIN(S2ALP)/2D0
      SA = SIN(ALP)
      CA = COS(ALP)
      SQBMA = (SINB*CA - COSB*SA)**2
  130 XIN = 1D0
  140 CONTINUE
 
      RETURN
      END
