--- /dev/null
+#include "isajet/pilot.h"
+ SUBROUTINE SSMSSM(XMG,XMU,XMHA,XTANB,XMQ1,XMDR,XMUR,
+ $XML1,XMER,XMQ2,XMSR,XMCR,XML2,XMMR,XMQ3,XMBR,XMTR,
+ $XML3,XMLR,XAT,XAB,XAL,XM1,XM2,XMT,IALLOW,IMODEL)
+C-----------------------------------------------------------------------
+C
+C Calculate MSSM masses and decays using parameters:
+C XM1 = U(1) mass
+C > 1e19: use scaling from XMG
+C XM2 = SU(2) mass
+C > 1e19: use scaling from XMG
+C XMG = gluino mass
+C XMQ1,... = 1st gen. su(2) soft squark mass,...
+C XMTL = m(stop-left)
+C XMTR = m(stop-right)
+C XMBR = m(sbot-right)
+C XML1 = left selectron mass
+C XMER = right selectron mass
+C XMN1 = 1st ge. sneutrino mass
+C XTANB = v/v' = ratio of vev's
+C XMU = -2*m_1 = SUSY Higgs mass
+C XMHA = m(pseudo-scalar-Higgs)
+C XMT = m(top)
+C XAT = stop trilinear coupling
+C XAB = sbottom trilinear coupling
+C XAL = stau trilinear coupling
+C IALLOW = 0 for valid point, 1 otherwise
+C IMODEL = 1 for SUGRA or MSSM, 2 for GMSB
+C
+C Program outline:
+C SSMSSM: Initialize standard model parameters in /SSSM/ and
+C SUSY parameters in /SSPAR/.
+C SSMASS: Calculate dependent SUSY masses and mixings.
+C SSTPBF: Calculate top decays; save in /SSMODE/.
+C SSSTBF: Calculate stop decays; save in /SSMODE/.
+C SSGLBF: Calcualte gluino decays; save in /SSMODE/.
+C SSQKBF: Calculate squark decays; save in /SSMODE/.
+C SSWZBF: Calculate gaugino decays; save in /SSMODE/.
+C SSHIBF: Calculate Higgs decays; save in /SSMODE/.
+C
+C Notes:
+C 1) All particle ID codes are defined with symbolic names in
+C /SSTYPE/, making it easy to change them.
+C
+C 2) /SSMODE/ contains the parent, the daughters, the width, and
+C the branching ratio for each mode. Decay modes for a given parent
+C need not be adjacent, so they must be sorted at the end.
+C
+C 3) Some of Baer's original routines used single precision and others
+C double precision. To accomodate this, the variable names used in
+C /SSSM/ and /SSPAR/ have all been changed to longer, more
+C descriptive ones.
+C
+C 4) All routines have been strongly typed.
+C
+C Source: H. Baer, et al.
+C Modified: F. Paige, Aug. 1992
+C-----------------------------------------------------------------------
+#if defined(CERNLIB_IMPNONE)
+ IMPLICIT NONE
+#endif
+#include "isajet/sslun.inc"
+#include "isajet/ssmode.inc"
+#include "isajet/sssm.inc"
+#include "isajet/sspar.inc"
+#include "isajet/dkyss3.inc"
+C
+ REAL XR21,PI,SR2
+ REAL XMG,XMU,XMHA,XTANB,XMQ1,XMDR,XMUR,XML1,XMER,XMQ2,XMSR,
+ $XMCR,XML2,XMMR,XMQ3,XMBR,XMTR,XML3,XMLR,XAT,XAB,XAL,XM1,XM2,
+ $XMT,MU1,MU2,BETA,COS2B
+ INTEGER IALLOW,MHLNEG,MHCNEG,IMODEL
+C
+ NSSMOD=0
+C
+C Standard model and SUSY parameters
+C
+ IALLOW=0
+ XR21=1./XTANB
+ PI=4.*ATAN(1.)
+ SR2=SQRT(2.)
+ AMDN=0.0099
+ AMUP=0.0056
+ AMST=0.199
+ AMCH=1.35
+ AMBT=5.0
+ AMTP=XMT
+ AME=0.511E-3
+ AMMU=0.105
+ AMTAU=1.777
+ AMW=80.0
+ AMZ=91.17
+ GAMW=2.12
+ GAMZ=2.487
+ ALFAEM=1./128.
+ SN2THW=0.232
+ ALFA2=ALFAEM/SN2THW
+ BETA=ATAN(XTANB)
+ COS2B=COS(2*BETA)
+C
+C SU(2) and U(1) gaugino masses are reset in SSMASS if
+C they are > 1e19.
+C
+ MU2=XM2
+ MU1=XM1
+C Set 2nd gen soft terms equal to 1st gen. soft terms
+c unless previously set by user.
+ IF (XMQ2.GE.1.E19) THEN
+ XMQ2=XMQ1
+ XMSR=XMDR
+ XMCR=XMUR
+ XML2=XML1
+ XMMR=XMER
+ END IF
+C
+C The results can be quite sensitive to the choice of the
+C 4-flavor QCD scale ALQCD4 and the expression for the QCD
+C coupling ALFA3. Select among the following lines:
+C
+ ALQCD4=0.177
+ ALFA3=0.12
+C
+C Calculate simple masses; other masses via SSMASS
+ AMGLSS=XMG
+ AMULSS=SQRT(XMQ1**2+AMUP**2+(.5-2.*SN2THW/3.)*AMZ**2*COS2B)
+ AMURSS=SQRT(XMUR**2+AMUP**2+2./3.*SN2THW*AMZ**2*COS2B)
+ AMDLSS=SQRT(XMQ1**2+AMDN**2+(-.5+SN2THW/3.)*AMZ**2*COS2B)
+ AMDRSS=SQRT(XMDR**2+AMDN**2-1./3.*SN2THW*AMZ**2*COS2B)
+ AMCLSS=SQRT(XMQ2**2+AMCH**2+(.5-2.*SN2THW/3.)*AMZ**2*COS2B)
+ AMCRSS=SQRT(XMCR**2+AMCH**2+2./3.*SN2THW*AMZ**2*COS2B)
+ AMSLSS=SQRT(XMQ2**2+AMST**2+(-.5+SN2THW/3.)*AMZ**2*COS2B)
+ AMSRSS=SQRT(XMSR**2+AMST**2-1./3.*SN2THW*AMZ**2*COS2B)
+ AMELSS=SQRT(XML1**2+AME**2-(.5-SN2THW)*AMZ**2*COS2B)
+ AMERSS=SQRT(XMER**2+AME**2-SN2THW*AMZ**2*COS2B)
+ AMMLSS=SQRT(XML2**2+AMMU**2-(.5-SN2THW)*AMZ**2*COS2B)
+ AMMRSS=SQRT(XMMR**2+AMMU**2-SN2THW*AMZ**2*COS2B)
+ AMN1SS=SQRT(XML1**2+.5*AMZ**2*COS2B)
+ AMN2SS=SQRT(XML2**2+.5*AMZ**2*COS2B)
+ AMN3SS=SQRT(XML3**2+.5*AMZ**2*COS2B)
+ AMTLSS=XMQ3
+ AMTRSS=XMTR
+ AMBLSS=XMQ3
+ AMBRSS=XMBR
+ AMLLSS=XML3
+ AMLRSS=XMLR
+ AMHA=XMHA
+ AAT=XAT
+ AAB=XAB
+ AAL=XAL
+ TWOM1=-XMU
+ RV2V1=XR21
+C
+C Calculate mass eigenstates and check Z1SS = LSP
+C
+ CALL SSMASS(MU1,MU2,IALLOW,1,MHLNEG,MHCNEG,IMODEL)
+ IF (MHLNEG.EQ.1.OR.MHCNEG.EQ.1) IALLOW=10
+C IF(IALLOW.NE.0) RETURN
+C
+C Initialize counters for matrix elements
+C Calculate decay widths and branching rations
+C
+ NMSS3=0
+ NPSS3=0
+ CALL SSTPBF
+ CALL SSGLBF
+ CALL SSQKBF
+ CALL SSSTBF
+ CALL SSLPBF
+ CALL SSWZBF
+ CALL SSHIBF
+C
+ RETURN
+ END