ISAJET/isasusy/sugra.F

   1 #include "isajet/pilot.h"
   2 C--------------------------------------------------------------------
   3       SUBROUTINE SUGRA(M0,MHF,A0,TANB,SGNMU,MT,IMODEL)
   4 C--------------------------------------------------------------------
   5 C
   6 C     Calculate supergravity spectra for ISAJET using as inputs
   7 C     M0    = M_0       = common scalar mass at GUT scale
   8 C     MHF   = M_(1/2)   = common gaugino mass at GUT scale
   9 C     A0    = A_0       = trilinear soft breaking parameter at GUT scale
  10 C     TANB  = tan(beta) = ratio of vacuum expectation values v_1/v_2
  11 C     SGNMU = sgn(mu)   = +-1 = sign of Higgsino mass term
  12 C     MT    = M_t       = mass of t quark
  13 C     M0    = Lambda    = ratio of vevs <F>/<S>
  14 C     MHF   = M_Mes     = messenger scale
  15 C     A0    = n_5       = number of messenger fields
  16 C     IMODEL            = 1 for SUGRA model
  17 C                       = 2 for GMSB model
  18 C                       = 7 for AMSB model
  19 C
  20 C     Uses Runge-Kutta method to integrate RGE's from M_Z to M_GUT
  21 C     and back, putting in correct thresholds. For the first iteration
  22 C     only the first 6 couplings are included and a common threshold
  23 C     is used.
  24 C
  25 C     See /SUGMG/ for definitions of couplings and masses.
  26 C
  27 #if defined(CERNLIB_IMPNONE)
  28       IMPLICIT NONE
  29 #endif
  30 #include "isajet/sslun.inc"
  31 #include "isajet/sspar.inc"
  32 #include "isajet/sssm.inc"
  33 #include "isajet/sugxin.inc"
  34 #include "isajet/sugmg.inc"
  35 #include "isajet/sugpas.inc"
  36 #include "isajet/sugnu.inc"
  37 #include "isajet/ssinf.inc"
  38       REAL GY(7),W1(21),G(29),W2(87)
  39       REAL G0(29)
  40       COMPLEX*16 SSB0,SSB1
  41       DOUBLE PRECISION DDILOG,XLM
  42       INTEGER IG(29)
  43       EXTERNAL SURG06,SURG26
  44       REAL M0,MHF,A0,TANB,SGNMU,MT,XLAMGM,XMESGM,XN5GM
  45       INTEGER NSTEP
  46       REAL M2,SUALFE,SUALFS,Q,T,A1I,AGUT,A3I,A2I,MTMT,ASMT,DT,
  47      $TGUT,TZ,GGUT,SIG2,SIG1,MH1S,MH2S,AGUTI,
  48      $MUS,MBMZ,MB,MTAU,MZ,MW,SR2,PI,ALEM,MTAMZ,
  49      $MTAMB,MTAMTA,MBMB,ASMB,BETA,COTB,SINB,COS2B,COSB,XC,
  50      $MSN,MG,MT1,MT2,MB1,MB2,MW1,MW2,AMU,BTHAT,BBHAT,BLHAT,AM2
  51       INTEGER II,I,J,IMODEL
  52       REAL G0SAVE(26),DELG0,DELLIM,THRF,THRG,DY,QOLD
  53       INTEGER MXITER,NSTEP0
  54       COMPLEX*16 ZZZ
  55       REAL*8 REAL8
  56 C
  57       DATA MZ/91.187/,MTAU/1.777/,MB/4.9/,ALEM/.0078186/
  58 C          This choice is a compromise between precision and speed:
  59       DATA MXITER/20/,NSTEP0/200/,DELLIM/2.E-2/
  60 C
  61 C          Define REAL(COMPLEX*16) for g77. This might need to be
  62 C          changed for 64-bit machines?
  63 C
  64       REAL8(ZZZ)=DREAL(ZZZ)
  65 C
  66 C          Save input parameters
  67 C
  68       XSUGIN(1)=M0
  69       XSUGIN(2)=MHF
  70       XSUGIN(3)=A0
  71       XSUGIN(4)=TANB
  72       XSUGIN(5)=SGNMU
  73       XSUGIN(6)=MT
  74       XLAMGM=M0
  75       XMESGM=MHF
  76       XN5GM=A0
  77       XGMIN(1)=XLAMGM
  78       XGMIN(2)=XMESGM
  79       XGMIN(3)=XN5GM
  80       XGMIN(4)=TANB
  81       XGMIN(5)=SGNMU
  82       XGMIN(6)=MT
  83       IF (XGMIN(12).EQ.0.) XGMIN(12)=XN5GM
  84       IF (XGMIN(13).EQ.0.) XGMIN(13)=XN5GM
  85       IF (XGMIN(14).EQ.0.) XGMIN(14)=XN5GM
  86 C
  87 C          Compute gauge mediated threshold functions
  88 C
  89       IF (IMODEL.EQ.2) THEN
  90         XLM=XLAMGM/XMESGM
  91         THRF=((1.D0+XLM)*(LOG(1.D0+XLM)-2*DDILOG(XLM/(1.D0+XLM))+
  92      ,        .5*DDILOG(2*XLM/(1.D0+XLM)))+
  93      ,       (1.D0-XLM)*(LOG(1.D0-XLM)-2*DDILOG(-XLM/(1.D0-XLM))+
  94      ,        .5*DDILOG(-2*XLM/(1.D0-XLM))))/XLM**2
  95         THRG=((1.D0+XLM)*LOG(1.D0+XLM)+(1.D0-XLM)*LOG(1.D0-XLM))/XLM**2
  96       END IF
  97 C
  98 C          Initialize standard model parameters in /SSSM/:
  99 C
 100       AMUP=0.0056
 101       AMDN=0.0099
 102       AMST=0.199
 103       AMCH=1.35
 104       AMBT=5.0
 105       AMTP=MT
 106       AMT=MT
 107       AME=0.511E-3
 108       AMMU=0.105
 109       AMTAU=1.777
 110       AMZ=91.17
 111       GAMW=2.12
 112       GAMZ=2.487
 113       ALFAEM=1./128.
 114       SN2THW=0.232
 115       ALFA2=ALFAEM/SN2THW
 116       ALQCD4=0.177
 117       ALFA3=0.118
 118 C
 119       NOGOOD=0
 120       ITACHY=0
 121       PI=4.*ATAN(1.)
 122       SR2=SQRT(2.)
 123       XW=.2324-1.03E-7*(MT**2-138.**2)
 124       MW=MZ*SQRT(1.-XW)
 125       AMW=MW
 126       A1MZ=5*ALEM/3./(1.-XW)
 127       A2MZ=ALEM/XW
 128       G2=SQRT(4*PI*A2MZ)
 129       GP=SQRT(3./5.*A1MZ*4.*PI)
 130       XTANB=TANB
 131       COTB=1./TANB
 132       BETA=ATAN(TANB)
 133       SINB=SIN(BETA)
 134       COSB=COS(BETA)
 135       SIN2B=SIN(2*BETA)
 136       COS2B=COS(2*BETA)
 137       IF (IMODEL.EQ.1) THEN
 138         MSUSY=SQRT(M0**2+4*MHF**2)
 139       ELSE IF (IMODEL.EQ.2) THEN
 140         MSUSY=XLAMGM/100.
 141       ELSE IF (IMODEL.EQ.7) THEN
 142         MSUSY=SQRT(M0**2+(.01*MHF)**2)
 143       END IF
 144 C     USE PIERCE PRESCRIPTION FOR MAGNITUDE OF VEV
 145 C      VEV=SR2*(248.6+0.9*LOG(MSUSY/AMZ)
 146 C      V=SQRT(VEV**2/(1.+COTB))
 147 C     PREVIOUS PRESCRIPTION
 148       V=SQRT(2*MW**2/G2**2/(1.+COTB**2))
 149       VP=V/TANB
 150       VEV=SQRT(V**2+VP**2)
 151 C
 152 C          Compute m(tau), m(b) at z scale using qcd, qed
 153 C
 154       MTAMTA=MTAU*(1.-SUALFE(MTAU**2)/PI)
 155       MTAMB=MTAMTA*(SUALFE(MB**2)/SUALFE(MTAU**2))**(-27./76.)
 156       MTAMZ=MTAMB*(SUALFE(MZ**2)/SUALFE(MB**2))**(-27./80.)
 157       FTAMZ=MTAMZ/COSB/VEV
 158       ASMB=SUALFS(MB**2,.36,MT,3)
 159       MBMB=MB*(1.-4*ASMB/3./PI)
 160       ASMZ=SUALFS(MZ**2,.36,MT,3)
 161       MBMZ=MBMB*(ASMZ/ASMB)**(12./23.)*
 162      $      (SUALFE(MZ**2)/SUALFE(MB**2))**(-3./80.)
 163       FBMZ=MBMZ/COSB/VEV
 164       ASMT=SUALFS(MT**2,.36,MT,3)
 165       MTMT=MT/(1.+4*ASMT/3./PI+(16.11-1.04*(5.-6.63/MT))*(ASMT/PI)**2)
 166       FTMT=MTMT/SINB/VEV
 167       FNMZ=SQRT(XNRIN(2)*XNRIN(1)/(SINB*VEV)**2)
 168       AMNRMJ=XNRIN(2)
 169 C
 170 C          Run the 3 gauge and 3 Yukawa's up to find M_GUT ,A_GUT and
 171 C          Yukawa_GUT
 172 C
 173 C
 174       NSTEP=NSTEP0
 175       GY(1)=SQRT(4*PI*A1MZ)
 176       GY(2)=SQRT(4*PI*A2MZ)
 177       GY(3)=SQRT(4*PI*ALFA3)
 178       GY(4)=FTAMZ
 179       GY(5)=FBMZ
 180       GY(6)=0.
 181       GY(7)=0.
 182       IF (IMODEL.EQ.1.OR.IMODEL.EQ.7) THEN
 183         IF (XSUGIN(7).EQ.0.) THEN
 184           MGUT=1.E19
 185         ELSE
 186           MGUT=XSUGIN(7)
 187         END IF
 188       ELSE IF (IMODEL.EQ.2) THEN
 189         MGUT=XMESGM
 190       END IF
 191       TZ=LOG(MZ/MGUT)
 192       TGUT=0.
 193       DT=(TGUT-TZ)/FLOAT(NSTEP)
 194       DO 200 II=1,NSTEP
 195         T=TZ+(TGUT-TZ)*FLOAT(II-1)/FLOAT(NSTEP)
 196         Q=MGUT*EXP(T)
 197         IF (Q.GT.MT.AND.GY(6).EQ.0.) GY(6)=FTMT
 198         IF (Q.GT.XNRIN(2).AND.GY(7).EQ.0.) GY(7)=FNMZ
 199         CALL RKSTP(7,DT,T,GY,SURG06,W1)
 200         A1I=4*PI/GY(1)**2
 201         A2I=4*PI/GY(2)**2
 202         A3I=4*PI/GY(3)**2
 203         IF (GY(5).GT.10..OR.GY(6).GT.10..OR.GY(7).GT.10.) THEN
 204           NOGOOD=4
 205           GO TO 100
 206         END IF
 207         IF (A1I.LT.A2I.AND.XSUGIN(7).EQ.0.) GO TO 10
 208 200   CONTINUE
 209       IF (MGUT.EQ.1.E19) THEN
 210       WRITE(LOUT,*) 'SUGRA: NO UNIFICATION FOUND'
 211       GO TO 100
 212       END IF
 213 10    IF (XSUGIN(7).EQ.0.) THEN
 214         MGUT=Q
 215       ELSE
 216         MGUT=XSUGIN(7)
 217       END IF
 218       AGUT=(GY(1)**2/4./PI+GY(2)**2/4./PI)/2.
 219       GGUT=SQRT(4*PI*AGUT)
 220       AGUTI=1./AGUT
 221       FTAGUT=GY(4)
 222       FBGUT=GY(5)
 223       FTGUT=GY(6)
 224       IF (XNRIN(1).EQ.0..AND.XNRIN(2).LT.1.E19) THEN
 225 C       UNIFY FN-FT
 226         FNGUT=GY(6)
 227       ELSE
 228         FNGUT=GY(7)
 229       END IF
 230 C
 231 C          Define parameters at GUT scale
 232 C
 233       DO 210 J=1,3
 234         IF (IMODEL.EQ.1) THEN
 235           G(J)=GY(J)
 236           G(J+6)=MHF
 237           G(J+9)=A0
 238         ELSE IF (IMODEL.EQ.2) THEN
 239           G(J)=GY(J)
 240           G(J+6)=XGMIN(11+J)*XGMIN(8)*THRG*(GY(J)/4./PI)**2*XLAMGM
 241           G(J+9)=0.
 242         END IF
 243 210   CONTINUE
 244 C     OVERWRITE ALFA_3 UNIFICATION TO GET ALFA_3(MZ) RIGHT
 245       IF (IMODEL.EQ.1.AND.IAL3UN.NE.0) G(3)=GGUT
 246       G(4)=FTAGUT
 247       G(5)=FBGUT
 248       G(6)=FTGUT
 249 C     IF NR MAJORANA MASS EXISTS, SET EXTRA NR RGE PARAMETERS
 250       IF (XNRIN(2).LT.1.E19) THEN
 251         G(27)=FNGUT
 252         G(28)=XNRIN(4)**2
 253         G(29)=XNRIN(3)
 254       ELSE
 255         G(27)=0.
 256         G(28)=0.
 257         G(29)=0.
 258       END IF
 259       IF (IMODEL.EQ.1) THEN
 260         DO 220 J=13,24
 261           G(J)=M0**2
 262 220     CONTINUE
 263 C       Set possible non-universal boundary conditions
 264       DO 230 J=1,6
 265         IF (XNUSUG(J).LT.1.E19) THEN
 266           G(J+6)=XNUSUG(J)
 267         END IF
 268 230   CONTINUE
 269       DO 231 J=7,18
 270         IF (XNUSUG(J).LT.1.E19) THEN
 271           G(J+6)=XNUSUG(J)**2
 272         END IF
 273 231   CONTINUE
 274       ELSE IF (IMODEL.EQ.2) THEN
 275        XC=2*THRF*XLAMGM**2
 276        DY=SQRT(3./5.)*GY(1)*XGMIN(11)
 277        G(13)=XC*(.75*XGMIN(13)*(GY(2)/4./PI)**4+.6*.25*
 278      ,XGMIN(12)*(GY(1)/4./PI)**4)+XGMIN(9)-DY
 279        G(14)=XC*(.75*XGMIN(13)*(GY(2)/4./PI)**4+.6*.25*
 280      ,XGMIN(12)*(GY(1)/4./PI)**4)+XGMIN(10)+DY
 281        G(15)=XC*(.6*XGMIN(12)*(GY(1)/4./PI)**4)+2*DY
 282        G(16)=XC*(.75*XGMIN(13)*(GY(2)/4./PI)**4+.6*.25*
 283      ,XGMIN(12)*(GY(1)/4./PI)**4)-DY
 284        G(17)=XC*(4*XGMIN(14)*(GY(3)/4./PI)**4/3.+.6*XGMIN(12)*
 285      ,(GY(1)/4./PI)**4/9.)+2*DY/3.
 286        G(18)=XC*(4*XGMIN(14)*(GY(3)/4./PI)**4/3.+.6*4*XGMIN(12)*
 287      ,(GY(1)/4./PI)**4/9.)-4*DY/3.
 288        G(19)=XC*(4*XGMIN(14)*(GY(3)/4./PI)**4/3.+.75*XGMIN(13)*
 289      ,(GY(2)/4./PI)**4+.6*XGMIN(12)*(GY(1)/4./PI)**4/36.)+DY/3.
 290        G(20)=G(15)
 291        G(21)=G(16)
 292        G(22)=G(17)
 293        G(23)=G(18)
 294        G(24)=G(19)
 295       ELSE IF (IMODEL.EQ.7) THEN
 296        G(1)=GY(1)
 297        G(2)=GY(2)
 298        G(3)=GY(3)
 299        BLHAT=G(4)*(-9*G(1)**2/5.-3*G(2)**2+3*G(5)**2+4*G(4)**2)
 300        BBHAT=G(5)*(-7*G(1)**2/15.-3*G(2)**2-16*G(3)**2/3.+
 301      ,             G(6)**2+6*G(5)**2+G(4)**2)
 302        BTHAT=G(6)*(-13*G(1)**2/15.-3*G(2)**2-16*G(3)**2/3.+
 303      ,             6*G(6)**2+G(5)**2)
 304        G(7)=-33*MHF*G(1)**2/5./16./PI**2
 305        G(8)=-MHF*G(2)**2/16./PI**2
 306        G(9)=3*MHF*G(3)**2/16./PI**2
 307        G(10)=BLHAT*MHF/G(4)/16./PI**2
 308        G(11)=BBHAT*MHF/G(5)/16./PI**2
 309        G(12)=BTHAT*MHF/G(6)/16./PI**2
 310        G(13)=(-99*G(1)**4/50.-3*G(2)**4/2.+3*G(5)*BBHAT+G(4)*BLHAT)*
 311      ,        MHF**2/(16*PI**2)**2
 312        G(14)=(-99*G(1)**4/50.-3*G(2)**4/2.+3*G(6)*BTHAT)*
 313      ,        MHF**2/(16*PI**2)**2
 314        G(15)=(-198*G(1)**4/25.)*MHF**2/(16*PI**2)**2
 315        G(16)=(-99*G(1)**4/50.-3*G(2)**4/2.)*MHF**2/(16*PI**2)**2
 316        G(17)=(-22*G(1)**4/25.+8*G(3)**4)*MHF**2/(16*PI**2)**2
 317        G(18)=(-88*G(1)**4/25.+8*G(3)**4)*MHF**2/(16*PI**2)**2
 318        G(19)=(-11*G(1)**4/50.-3*G(2)**4/2.+8*G(3)**4)*
 319      ,        MHF**2/(16*PI**2)**2
 320        G(20)=(-198*G(1)**4/25.+2*G(4)*BLHAT)*MHF**2/(16*PI**2)**2
 321        G(21)=(-99*G(1)**4/50.-3*G(2)**4/2.+G(4)*BLHAT)*
 322      ,        MHF**2/(16*PI**2)**2
 323        G(22)=(-22*G(1)**4/25.+8*G(3)**4+2*G(5)*BBHAT)*
 324      ,MHF**2/(16*PI**2)**2
 325        G(23)=(-88*G(1)**4/25.+8*G(3)**4+2*G(6)*BTHAT)*
 326      ,MHF**2/(16*PI**2)**2
 327        G(24)=(-11*G(1)**4/50.-3*G(2)**4/2.+8*G(3)**4+G(5)*BBHAT+
 328      ,        G(6)*BTHAT)*MHF**2/(16*PI**2)**2
 329        DO 234 I=13,24
 330 234      G(I)=G(I)+M0**2
 331       END IF
 332       G(25)=0.
 333       G(26)=0.
 334       DO 235 I=1,29
 335         IG(I)=0
 336 235   CONTINUE
 337 C          Check for tachyonic sleptons at GUT scale
 338       IF (G(15).LT.0..OR.G(16).LT.0.) THEN
 339         ITACHY=1
 340       END IF
 341 C
 342 C          Initialize thresholds
 343 C
 344       MSS(1)=MSUSY
 345       MSS(2)=MSUSY
 346       MSS(17)=MSUSY
 347       MSS(27)=MSUSY
 348       MSS(31)=MSUSY
 349       MU=MSUSY
 350 C
 351 C          Evolve parameters from mgut to mz
 352 C
 353       TZ=LOG(MZ/MGUT)
 354       TGUT=0.
 355       DT=(TZ-TGUT)/FLOAT(NSTEP)
 356 C          Freeze Higgs parameters at HIGFRZ = Drees' value
 357 C          AMTLSS, AMTRSS initialized to 0 for later use in HIGFRZ
 358       IF (IMODEL.EQ.1) THEN
 359         HIGFRZ=SQRT(M0**2+3*MHF**2)
 360       ELSE IF (IMODEL.EQ.2) THEN
 361         HIGFRZ=MSUSY
 362       ELSE IF (IMODEL.EQ.7) THEN
 363         HIGFRZ=SQRT(M0**2+(.01*MHF)**2)
 364       END IF
 365       AMTLSS=0
 366       AMTRSS=0
 367       DO 240 II=1,NSTEP+2
 368         T=TGUT+(TZ-TGUT)*FLOAT(II-1)/FLOAT(NSTEP)
 369         QOLD=Q
 370         Q=MGUT*EXP(T)
 371         CALL RKSTP(29,DT,T,G,SURG26,W2)
 372         IF (Q.LT.AMNRMJ.AND.QOLD.GE.AMNRMJ.AND.FNMZ.EQ.0.) THEN
 373           FNMZ=G(27)
 374         END IF
 375         IF (Q.LT.AMNRMJ) THEN
 376           G(27)=0.
 377           G(28)=0.
 378           G(29)=0.
 379         END IF
 380         CALL SUGFRZ(Q,G,G0,IG)
 381         IF (NOGOOD.NE.0) GO TO 100
 382         IF (Q.LT.MZ) GO TO 20
 383 240   CONTINUE
 384 20    CONTINUE
 385       ASMZ=G0(3)**2/4./PI
 386 C          Electroweak breaking constraints; tree level
 387       MUS=(G0(13)-G0(14)*TANB**2)/(TANB**2-1.)-MZ**2/2.
 388       IF (MUS.LT.0.) THEN
 389         NOGOOD=2
 390         GO TO 100
 391       END IF
 392       MU=SQRT(MUS)*SIGN(1.,SGNMU)
 393       B=(G0(13)+G0(14)+2*MUS)*SIN2B/MU/2.
 394 C          Compute tree level masses
 395       CALL SUGMAS(G0,0,IMODEL)
 396       IF (NOGOOD.NE.0) GO TO 100
 397 C          Compute effective potential corrections
 398       CALL SUGEFF(G0,SIG1,SIG2)
 399       MH1S=G0(13)+SIG1
 400       MH2S=G0(14)+SIG2
 401       MUS=(MH1S-MH2S*TANB**2)/(TANB**2-1.)-MZ**2/2.
 402       IF (MUS.LT.0.) THEN
 403         NOGOOD=2
 404         GO TO 100
 405       END IF
 406       MU=SQRT(MUS)*SIGN(1.,SGNMU)
 407       B=(MH1S+MH2S+2*MUS)*SIN2B/MU/2.
 408 C
 409 C     Recompute weak scale Yukawa couplings including SUSY loops
 410 C     Follow formulae of Pierce et al. NPB491, 3 (1997)
 411 C
 412       M2=G0(8)
 413       AM2=ABS(M2)
 414       MSN=MSS(16)
 415       MG=MSS(1)
 416       MT1=MSS(12)
 417       MT2=MSS(13)
 418       MB1=MSS(10)
 419       MB2=MSS(11)
 420       MW1=ABS(MSS(27))
 421       MW2=ABS(MSS(28))
 422       AMU=ABS(MU)
 423       XLAM=LOG(MT**2)
 424 C          Be careful in using our convention vs Pierce et al.
 425       IF (ABS(COS(THETAT)).LT..707107) THEN
 426         MTMT=MT/(1.+4*ASMT/3./PI+(16.11-1.04*(5.-6.63/MT))*(ASMT/PI)**2
 427      $  -ASMT/3./PI*(REAL8(SSB1(MT**2,MG,MT1))+
 428      $  REAL8(SSB1(MT**2,MG,MT2))-SIN(2*THETAT)*MG/MT*
 429      $  (REAL8(SSB0(MT**2,MG,MT1))-REAL8(SSB0(MT**2,MG,MT2)))))
 430       ELSE
 431         MTMT=MT/(1.+4*ASMT/3./PI+(16.11-1.04*(5.-6.63/MT))*(ASMT/PI)**2
 432      $  -ASMT/3./PI*(REAL8(SSB1(MT**2,MG,MT2))+
 433      $  REAL8(SSB1(MT**2,MG,MT1))+SIN(2*THETAT)*MG/MT*
 434      $  (REAL8(SSB0(MT**2,MG,MT2))-REAL8(SSB0(MT**2,MG,MT1)))))
 435       END IF
 436       FTMT=MTMT/SINB/VEV
 437       XLAM=LOG(MZ**2)
 438       IF (ABS(COS(THETAB)).LT..707107) THEN
 439         MBMZ=MBMZ*(1.+ASMZ/3./PI*(REAL8(SSB1(MZ**2,MG,MB1))+
 440      $  REAL8(SSB1(MZ**2,MG,MB2))-SIN(2*THETAB)*MG/MB*
 441      $  (REAL8(SSB0(MZ**2,MG,MB1))-REAL8(SSB0(MZ**2,MG,MB2))))
 442      $  -FTMT**2*MU*(-AAT*TANB+MU)/16./PI**2/(MT1**2-MT2**2)*
 443      $  (REAL8(SSB0(MZ**2,AMU,MT1))-REAL8(SSB0(MZ**2,AMU,MT2)))+
 444      $  G2**2*MU*M2*TANB/16./PI**2/(AMU**2-M2**2)*
 445      $  (SIN(THETAT)**2*(REAL8(SSB0(MZ**2,AM2,MT1))-
 446      $  REAL8(SSB0(MZ**2,AMU,MT1)))+
 447      $  COS(THETAT)**2*(REAL8(SSB0(MZ**2,AM2,MT2))-
 448      $  REAL8(SSB0(MZ**2,AMU,MT2)))))
 449       ELSE
 450         MBMZ=MBMZ*(1.+ASMZ/3./PI*(REAL8(SSB1(MZ**2,MG,MB2))+
 451      $  REAL8(SSB1(MZ**2,MG,MB1))+SIN(2*THETAB)*MG/MB*
 452      $  (REAL8(SSB0(MZ**2,MG,MB2))-REAL8(SSB0(MZ**2,MG,MB1))))
 453      $  -FTMT**2*MU*(-AAT*TANB+MU)/16./PI**2/(MT2**2-MT1**2)*
 454      $  (REAL8(SSB0(MZ**2,AMU,MT2))-REAL8(SSB0(MZ**2,AMU,MT1)))+
 455      $  G2**2*MU*M2*TANB/16./PI**2/(AMU**2-M2**2)*
 456      $  (COS(THETAT)**2*(REAL8(SSB0(MZ**2,AM2,MT2))-
 457      $  REAL8(SSB0(MZ**2,AMU,MT2)))+
 458      $  SIN(THETAT)**2*(REAL8(SSB0(MZ**2,AM2,MT1))-
 459      $  REAL8(SSB0(MZ**2,AMU,MT1)))))
 460       END IF
 461       FBMZ=MBMZ/COSB/VEV
 462       MTAMZ=MTAMZ*(1.+G2**2*MU*M2*TANB/16./PI**2/(MUS-M2**2)*
 463      $(REAL8(SSB0(MZ**2,AM2,MSN))-REAL8(SSB0(MZ**2,AMU,MSN))))
 464       FTAMZ=MTAMZ/COSB/VEV
 465 C
 466 C          Iterate entire process, increasing NSTEP each time
 467 C          This time, freeze out parameters at sqrt(t_l t_r)
 468 C
 469       HIGFRZ=MAX(AMZ,(G0(23)*G0(24))**0.25)
 470       DO 300 I=1,MXITER
 471         DO 310 J=1,26
 472 310     G0SAVE(J)=G0(J)
 473         NSTEP=1.2*NSTEP
 474         CALL SUGRGE(M0,MHF,A0,TANB,SGNMU,MT,G,G0,IG,W2,NSTEP,IMODEL)
 475         IF(NOGOOD.NE.0) GO TO 100
 476         DELG0=0.
 477         DO 320 J=1,24
 478 320     DELG0=MAX(DELG0,ABS((G0(J)-G0SAVE(J))/G0(J)))
 479         IF(DELG0.LT.DELLIM) GO TO 400
 480 300   CONTINUE
 481       WRITE(LOUT,1000) MXITER
 482 1000  FORMAT(/' SUGRA WARNING: NO CONVERGENCE IN',I4,' ITERATIONS')
 483 C
 484 C          Save results
 485 C
 486 400   DO 410 I=1,26
 487         GSS(I)=G0(I)
 488 410   CONTINUE
 489       MGUTSS=MGUT
 490       AGUTSS=AGUT
 491       GGUTSS=GGUT
 492 C
 493 C          Fill XISAIN common block
 494 C
 495       XISAIN(1)=MSS(1)
 496       XISAIN(2)=MU
 497       XISAIN(3)=MSS(31)
 498       XISAIN(4)=TANB
 499       XISAIN(5)=SQRT(G0(19))
 500       XISAIN(6)=SQRT(G0(17))
 501       XISAIN(7)=SQRT(G0(18))
 502       XISAIN(8)=SQRT(G0(16))
 503       XISAIN(9)=SQRT(G0(15))
 504       XISAIN(10)=XISAIN(5)
 505       XISAIN(11)=XISAIN(6)
 506       XISAIN(12)=XISAIN(7)
 507       XISAIN(13)=XISAIN(8)
 508       XISAIN(14)=XISAIN(9)
 509       XISAIN(15)=SQRT(G0(24))
 510       XISAIN(16)=SQRT(G0(22))
 511       XISAIN(17)=SQRT(G0(23))
 512       XISAIN(18)=SQRT(G0(21))
 513       XISAIN(19)=SQRT(G0(20))
 514       XISAIN(20)=G0(12)
 515       XISAIN(21)=G0(11)
 516       XISAIN(22)=G0(10)
 517       XISAIN(23)=G0(7)
 518       XISAIN(24)=G0(8)
 519       M2=G0(8)
 520 100   RETURN
 521       END