pyslha.f

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyslha.f

C*********************************************************************
 
C...PYSLHA
C...Read/write spectrum or decay data from SLHA standard file(s).
C...P. Skands
 
C...MUPDA=1 : READ SPECTRUM ON LUN=IMSS(21)
C...MUPDA=2 : LOOK FOR DECAY TABLE FOR KF=KFORIG ON LUN=IMSS(22)
C...MUPDA=3 : WRITE SPECTRUM ON LUN=IMSS(23)
C...(MUPDA=4 : WRITE DECAY TABLE FOR KF=KFORIG ON LUN=IMSS(24))
C...MUPDA=5 : READ MASS FOR KF=KFORIG ONLY (WITH DECAY TABLE)
      SUBROUTINE pyslha(MUPDA,KFORIG,IRETRN)
 
C...Double precision and integer declarations.
      IMPLICIT DOUBLE PRECISION(a-h, o-z)
      IMPLICIT INTEGER(i-n)
      INTEGER pyk,pychge,pycomp
      parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
     &kexcit=4000000,kdimen=5000000)
C...Commonblocks.
      common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
      common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
      common/pydat3/mdcy(500,3),mdme(8000,2),brat(8000),kfdp(8000,5)
      common/pydat4/chaf(500,2)
      CHARACTER chaf*16
      CHARACTER*40 isaver,visaje
      common/pyint4/mwid(500),wids(500,5)
      SAVE /pydat1/,/pydat2/,/pydat3/,/pydat4/,/pyint4/
C...SUSY blocks
      common/pymssm/imss(0:99),rmss(0:99)
      common/pyssmt/zmix(4,4),umix(2,2),vmix(2,2),smz(4),smw(2),
     &sfmix(16,4),zmixi(4,4),umixi(2,2),vmixi(2,2)
      common/pymsrv/rvlam(3,3,3), rvlamp(3,3,3), rvlamb(3,3,3)
      SAVE /pymssm/,/pyssmt/,/pymsrv/
 
C...Local arrays, character variables and data.
      common/pylh3p/modsel(200),parmin(100),parext(200),rmsoft(0:100),
     &     au(3,3),ad(3,3),ae(3,3)
      common/pylh3c/cpro(2),cver(2)
      SAVE /pylh3p/,/pylh3c/
      dimension mmod(100),mspc(100),mdec(100)
C...MMOD: flags to set for each block read in.
C... 1: MODSEL     2: MINPAR     3: EXTPAR     4: SMINPUTS
C...MSPC: Flags to set for each block read in.
C... 1: MASS       2: NMIX       3: UMIX       4: VMIX       5: SBOTMIX
C... 6: STOPMIX    7: STAUMIX    8: HMIX       9: GAUGE     10: AU
C...11: AD        12: AE        13: YU        14: YD        15: YE
C...16: SPINFO    17: ALPHA     18: MSOFT     19: QNUMBERS
      CHARACTER cpro*12,cver*12,chnlin*6
      CHARACTER doc*11, chdum*120, chblck*60
      CHARACTER chinl*120,chkf*9,chtmp*16
      INTEGER verbos
      SAVE verbos
C...Date of last Change
      parameter(doc='05 Mar 2007')
C...MQREAD(0): Number of entries I in MQREAD
C...      (I): KF code for which a QNUMBERS block has been read.
      dimension idc(5),kfsusy(50),mqread(0:100)
      SAVE kfsusy,mqread
      DATA verbos /1/
      DATA nhello /0/
      DATA kfsusy/
     &1000001,1000002,1000003,1000004,1000005,1000006,
     &2000001,2000002,2000003,2000004,2000005,2000006,
     &1000011,1000012,1000013,1000014,1000015,1000016,
     &2000011,2000012,2000013,2000014,2000015,2000016,
     &1000021,1000022,1000023,1000025,1000035,1000024,
     &1000037,1000039,     25,     35,     36,     37,
     &      6,     24,     45,     46,1000045, 9*0/
      rmfun(ip)=pmas(pycomp(ip),1)
 
C...Hello World
      IF (nhello.EQ.0) THEN
        WRITE(mstu(11),5000) doc
        nhello=1
      ENDIF

C...SLHA file assumed opened by user on unit LFN, stored in IMSS(20
C...+MUPDA).
      lfn=imss(20+mupda)
      IF (mupda.EQ.5) lfn=imss(21)
      IF (mupda.EQ.0) lfn=imss(21)
C...Flag that we have not yet found whatever we were asked to find.
      iretrn=1
 
C...STOP IF LFN IS ZERO (i.e. if no LFN was given).
      IF (lfn.EQ.0) THEN
        WRITE(mstu(11),*) '* (PYSLHA:) No valid unit given in IMSS'
        goto 9999
      ENDIF
 
C...If told to read spectrum, first zero all previous information.
      IF (mupda.EQ.1) THEN
C...Zero all block read flags
        DO 100 m=1,100
          mmod(m)=0
          mspc(m)=0
          mdec(m)=0
  100   CONTINUE
C...Zero all (MSSM) masses, widths, and lifetimes in PYTHIA
        DO 110 isusy=1,36
          kc=pycomp(kfsusy(isusy))
          pmas(kc,1)=0d0
          pmas(kc,2)=0d0
          pmas(kc,3)=0d0
          pmas(kc,4)=0d0
  110   CONTINUE
C...Zero all (3rd gen sfermion + gaugino/higgsino) mixing matrices.
        DO 130 j=1,4
          sfmix(5,j) =0d0
          sfmix(6,j) =0d0
          sfmix(15,j)=0d0
          DO 120 l=1,4
            zmix(l,j) =0d0
            zmixi(l,j)=0d0
            IF (j.LE.2.AND.l.LE.2) THEN
              umix(l,j) =0d0
              umixi(l,j)=0d0
              vmix(l,j) =0d0
              vmixi(l,j)=0d0
            ENDIF
  120     CONTINUE
C...Zero signed masses.
          smz(j)=0d0
          IF (j.LE.2) smw(j)=0d0
  130   CONTINUE
C...NB: RMSS array not zeroed.
        WRITE(mstu(11),*)
     &       '* (PYSLHA:) Reading in SLHA spectrum from unit ', lfn
 
C...If reading decays, reset PYTHIA decay counters.
      ELSEIF (mupda.EQ.2) THEN
        kcc=100
        ndc=0
        brsum=0d0
        DO 140 kc=1,mstu(6)
          IF(kc.GT.100.AND.kchg(kc,4).GT.100) kcc=kc
          ndc=max(ndc,mdcy(kc,2)+mdcy(kc,3)-1)
  140   CONTINUE
      ELSEIF (mupda.EQ.5) THEN
C...Zero block read flags
        DO 150 m=1,100
          mspc(m)=0
 150    CONTINUE
      ENDIF
 
C............READ
C...(spectrum or look for decays of KF=KFORIG or MASS of KF=KFORIG
      IF(mupda.EQ.0.OR.mupda.EQ.1.OR.mupda.EQ.2.OR.mupda.EQ.5) THEN
C...Initialize program and version strings
        IF(mupda.EQ.1.OR.mupda.EQ.2) THEN
        cpro(mupda)=' '
        cver(mupda)=' '
        ENDIF
 
C...Initialize read loop
        merr=0
        nline=0
        chblck=' '
C...READ NEW LINE INTO CHINL. GOTO 300 AT END-OF-FILE.
  160   chinl=' '
        READ(lfn,'(A120)',end=300) chinl
C...Count which line number we're at.
        nline=nline+1
        WRITE(chnlin,'(I6)') nline
 
C...Skip comment and empty lines without processing.
        IF (chinl(1:1).EQ.'#'.OR.chinl.EQ.' ') goto 160
 
C...We assume all upper case below. Rewrite CHINL to all upper case.
        inl=0
        igood=0
  170   inl=inl+1
        IF (chinl(inl:inl).NE.'#') THEN
          DO 180 ich=97,122
            IF (char(ich).EQ.chinl(inl:inl)) chinl(inl:inl)=char(ich-32)
  180     CONTINUE
C...Extra safety. Chek for sensible input on line
          IF (igood.EQ.0) THEN
            DO 190 ich=48,90
              IF (char(ich).EQ.chinl(inl:inl)) igood=1
  190       CONTINUE
          ENDIF
          IF (inl.LT.120) goto 170
        ENDIF
        IF (igood.EQ.0) goto 160
 
C...Check for BLOCK begin statement (spectrum).
        IF (chinl(1:1).EQ.'B') THEN
          merr=0
          READ(chinl,'(A6,A)',err=460) chdum,chblck
C...Check if another of this type of block was already read.
C...(logarithmic interpolation not yet implemented, so duplicates always
C...give errors)
          IF (chblck(1:6).EQ.'MODSEL'.AND.mmod(1).NE.0) merr=7
          IF (chblck(1:6).EQ.'MINPAR'.AND.mmod(2).NE.0) merr=7
          IF (chblck(1:6).EQ.'EXTPAR'.AND.mmod(3).NE.0) merr=7
          IF (chblck(1:8).EQ.'SMINPUTS'.AND.mmod(4).NE.0) merr=7
          IF (chblck(1:4).EQ.'MASS'.AND.mspc(1).NE.0) merr=7
          IF (chblck(1:4).EQ.'NMIX'.AND.mspc(2).NE.0) merr=7
          IF (chblck(1:4).EQ.'UMIX'.AND.mspc(3).NE.0) merr=7
          IF (chblck(1:4).EQ.'VMIX'.AND.mspc(4).NE.0) merr=7
          IF (chblck(1:7).EQ.'SBOTMIX'.AND.mspc(5).NE.0) merr=7
          IF (chblck(1:7).EQ.'STOPMIX'.AND.mspc(6).NE.0) merr=7
          IF (chblck(1:7).EQ.'STAUMIX'.AND.mspc(7).NE.0) merr=7
          IF (chblck(1:4).EQ.'HMIX'.AND.mspc(8).NE.0) merr=7
          IF (chblck(1:5).EQ.'ALPHA'.AND.mspc(17).NE.0) merr=7
          IF (chblck(1:5).EQ.'AU'.AND.mspc(10).NE.0) merr=7
          IF (chblck(1:5).EQ.'AD'.AND.mspc(11).NE.0) merr=7
          IF (chblck(1:5).EQ.'AE'.AND.mspc(12).NE.0) merr=7
          IF (chblck(1:5).EQ.'MSOFT'.AND.mspc(18).NE.0) merr=7
C...Check for new particles
          IF (chblck(1:8).EQ.'QNUMBERS'.OR.chblck(1:8).EQ.'PARTICLE')
     &        THEN
            mspc(19)=mspc(19)+1
C...Read PDG code
            READ(chblck(9:60),*) kfq

            DO 121 mq=1,mqread(0)
              IF (mqread(mq).EQ.kfq) THEN
                merr=17
                goto 290
              ENDIF
 121        CONTINUE
            WRITE(mstu(11),'(A,I9,A,F12.3)')
     &           ' * (PYSLHA:) Reading in '//chblck(1:8)//
     &           ' for KF =',kfq
            mqread(0)=mqread(0)+1
            mqread(mqread(0))=kfq
            mspc(19)=mspc(19)+1
            kcq=pycomp(kfq)
            IF (kcq.EQ.0) THEN
              DO 123 kct=100,mstu(6)
                IF(kchg(kct,4).GT.100) kcq=kct
 123          CONTINUE
              kcq=kcq+1
              kcc=kcq
              kchg(kcq,4)=kfq              
C...First write PDG code as name
              WRITE(chtmp,*) kfq
C...Then look for real name
              icmt=9
 90           icmt=icmt+1
              IF (chblck(icmt:icmt).NE.'#'.AND.icmt.LT.59) goto 90
              IF (icmt.LT.59) THEN
                READ(chblck(icmt+1:60),'(A)',err=95) chdum
                IF (chdum.NE.' ') chtmp=chdum
              ENDIF
 95           IF (chtmp(1:1).EQ.' ') THEN
                READ(chtmp,'(1x,A)') chaf(kcq,1)
              ELSE
                READ(chtmp,'(A)') chaf(kcq,1)
              ENDIF
              mstu(20)=0
C...Set stable for now
              pmas(kcq,2)=1d-6
              mwid(kcq)=0
              mdcy(kcq,1)=0
              mdcy(kcq,2)=0
              mdcy(kcq,3)=0
            ELSE
              WRITE(mstu(11),*)
     &           '* (PYSLHA:) KF =',kfq,' already exists: ',
     &             chaf(kcq,1), '. Entry ignored.'
              merr=7
            ENDIF
          ENDIF
C...Finalize this line and read next.
          goto 290
C...Check for DECAY begin statement (decays).
        ELSEIF (chinl(1:1).EQ.'D') THEN
          merr=0
          brsum=0d0
          chblck='DECAY'
C...Read KF code and WIDTH
          mpsign=1
          READ(chinl(7:inl),*,err=470) kf, width
          IF (kf.LE.0) THEN
            kf=-kf
            mpsign=-1
          ENDIF
C...If this is not the KF we're looking for...
          IF (kf.NE.kforig.OR.mupda.NE.2) THEN
C...Set block skip flag and read next line.
            merr=16
            goto 290
          ENDIF
 
C...Determine PYTHIA KC code of particle
          kcrep=0
          IF(kf.LE.100) THEN
            kcrep=kf
          ELSE
            DO 200 kcr=101,kcc
              IF(kchg(kcr,4).EQ.kf) kcrep=kcr
  200       CONTINUE
          ENDIF
          kc=kcrep
          IF (kcrep.NE.0) THEN
C...Particle is already known. Don't do anything yet.
          ELSE
C...  Add new particle. Actually, this should not happen.
C...  New particles should be added already when reading the spectrum
C...  information, so go under previously stable category.
            kcc=kcc+1
            kc=kcc
          ENDIF
 
          IF (width.LE.0d0) THEN
C...Stable (i.e. LSP)
            WRITE(mstu(11),*)
     &           '* (PYSLHA:) Reading in SLHA stable particle: ',
     &           chaf(kcrep,1)
            IF (width.LT.0d0) THEN
              CALL pyerrm(19,'(PYSLHA:) Negative width forced to'//
     &             ' zero !')
              width=0d0
            ENDIF
            pmas(kc,2)=1d-6
            mwid(kc)=0
            mdcy(kc,1)=0
C...Ignore any decay lines that may be present for this KF
            merr=16
            mdcy(kc,2)=0
            mdcy(kc,3)=0
C...Return ok
            iretrn=0
          ENDIF
C...Finalize and start reading in decay modes.
          goto 290
        ELSEIF (mod(merr,10).GE.6) THEN
C...If ignore block flag set, skip directly to next line.
          goto 160
        ENDIF
 
C...READ SPECTRUM
        IF (mupda.EQ.0.AND.merr.EQ.0) THEN
          IF (chblck(1:8).EQ.'QNUMBERS'.OR.chblck(1:8).EQ.'PARTICLE') 
     &        THEN
            READ(chinl,*) indx, ival
            IF (indx.EQ.1) kchg(kcq,1)=ival
            IF (indx.EQ.3) kchg(kcq,2)=0
            IF (indx.EQ.3.AND.ival.EQ.3) kchg(kcq,2)=1
            IF (indx.EQ.3.AND.ival.EQ.-3) kchg(kcq,2)=-1
            IF (indx.EQ.3.AND.ival.EQ.8) kchg(kcq,2)=2
            IF (indx.EQ.4) THEN
              kchg(kcq,3)=ival
              IF (ival.EQ.1) THEN 
                chtmp=chaf(kcq,1)
                IF (chtmp.EQ.' ') THEN
                  WRITE(chaf(kcq,1),*) kchg(kcq,4)
                  WRITE(chaf(kcq,2),*) -kchg(kcq,4)
                ELSE
                  ilast=17
 116              ilast=ilast-1
                  IF (chtmp(ilast:ilast).EQ.' ') goto 116
                  IF (chtmp(ilast:ilast).EQ.'+') THEN
                    chtmp(ilast:ilast)='-'
                  ELSE
                    chtmp(ilast+1:min(16,ilast+4))='bar'
                  ENDIF
                  chaf(kcq,2)=chtmp
                ENDIF
              ENDIF
            ENDIF
          ELSE
            merr=8
          ENDIF
        ELSEIF ((mupda.EQ.1.OR.mupda.EQ.5).AND.merr.EQ.0) THEN
C...MASS: Mass spectrum
          IF (chblck(1:4).EQ.'MASS') THEN
            READ(chinl,*) kf, val
            merr=1
            kc=0
            IF (mupda.EQ.1.OR.kf.EQ.kforig) THEN
C...Read in masses for anything
              merr=0
              kc=pycomp(kf)
              IF (kc.NE.0) THEN
                mspc(1)=mspc(1)+1
                pmas(kc,1) = abs(val)
                IF (mupda.EQ.5) THEN
                  WRITE(mstu(11),'(A,I9,A,F12.3)')
     &                 ' * (PYSLHA:) Reading in MASS entry for KF =',
     &                 kf, ', pole mass =', val
                  iretrn=0
                ENDIF
C...  Signed masses
                IF (kf.EQ.1000021.AND.mspc(18).EQ.0) rmss(3)=val
                IF (kf.EQ.1000022) smz(1)=val
                IF (kf.EQ.1000023) smz(2)=val
                IF (kf.EQ.1000025) smz(3)=val
                IF (kf.EQ.1000035) smz(4)=val
                IF (kf.EQ.1000024) smw(1)=val
                IF (kf.EQ.1000037) smw(2)=val
              ENDIF
            ELSEIF (mupda.EQ.5) THEN
              merr=0
            ENDIF
          ELSEIF (mupda.EQ.5) THEN
C...Only read MASS if MUPDA = 5. Skip any other blocks.
            merr=8
          ELSEIF (chblck(1:8).EQ.'QNUMBERS'.OR.
     &          chblck(1:8).EQ.'PARTICLE') THEN
C...Don't print a warning for QNUMBERS when reading spectrum
            merr=8
C...  MODSEL: Model selection and global switches
          ELSEIF (chblck(1:6).EQ.'MODSEL') THEN
            READ(chinl,*) indx, ival
            IF (indx.LE.200.AND.indx.GT.0) THEN
              modsel(indx)=ival
              mmod(1)=mmod(1)+1
              IF (indx.EQ.3.AND.ival.EQ.1) THEN
C...  Switch on NMSSM
                WRITE(mstu(11),*) '* (PYSLHA:) switching on NMSSM'
                imss(13)=max(1,imss(13))
C...  Add NMSSM states if not already done
 
                kfn=25
                kcn=kfn
                chaf(kcn,1)='H_10'
                chaf(kcn,2)=' '
 
                kfn=35
                kcn=kfn
                chaf(kcn,1)='H_20'
                chaf(kcn,2)=' '
 
                kfn=45
                kcn=kfn
                chaf(kcn,1)='H_30'
                chaf(kcn,2)=' '
 
                kfn=36
                kcn=kfn
                chaf(kcn,1)='A_10'
                chaf(kcn,2)=' '
 
                kfn=46
                kcn=kfn
                chaf(kcn,1)='A_20'
                chaf(kcn,2)=' '
 
                kfn=1000045
                kcn=pycomp(kfn)
                IF (kcn.EQ.0) THEN
                  DO 234 kct=100,mstu(6)
                    IF(kchg(kct,4).GT.100) kcn=kct
 234              CONTINUE
                  kcn=kcn+1
                  kchg(kcn,4)=kfn
                  mstu(20)=0
                ENDIF
C...  Set stable for now
                pmas(kcn,2)=1d-6
                mwid(kcn)=0
                mdcy(kcn,1)=0
                mdcy(kcn,2)=0
                mdcy(kcn,3)=0
                chaf(kcn,1)='~chi_50'
                chaf(kcn,2)=' '
              ENDIF
            ELSE
              merr=1
            ENDIF
C...MINPAR: Minimal model parameters
          ELSEIF (chblck(1:6).EQ.'MINPAR') THEN
            IF (modsel(1).NE.0) THEN
              READ(chinl,*) indx, val
              IF (indx.LE.100.AND.indx.GT.0) THEN
                parmin(indx)=val
                mmod(2)=mmod(2)+1
              ELSE
                merr=1
              ENDIF
            ELSEIF (mmod(3).NE.0) THEN
              WRITE(mstu(11),*)
     &             '* (PYSLHA:) MINPAR after EXTPAR !'
              merr=1
            ELSE
              WRITE(mstu(11),*)
     &             '* (PYSLHA:) Reading MINPAR, but no MODSEL !' 
              merr=1
            ENDIF
C...tan(beta)
            IF (indx.EQ.3) rmss(5)=val
C...EXTPAR: non-minimal model parameters.
          ELSEIF (chblck(1:6).EQ.'EXTPAR') THEN
            IF (mmod(1).NE.0) THEN
              READ(chinl,*) indx, val
              IF (indx.LE.200.AND.indx.GT.0) THEN
                parext(indx)=val
                mmod(3)=mmod(3)+1
              ELSE
                merr=1
              ENDIF
            ELSE
              WRITE(mstu(11),*)
     &             '* (PYSLHA:) Reading EXTPAR, but no MODSEL !'
              merr=1
            ENDIF
C...tan(beta)
            IF (indx.EQ.25) rmss(5)=val
          ELSEIF (chblck(1:8).EQ.'SMINPUTS') THEN
            READ(chinl,*) indx, val
            IF (indx.LE.3.OR.indx.EQ.5.OR.indx.GE.7) THEN
              merr=1
            ELSEIF (indx.EQ.4) THEN
              pmas(pycomp(23),1)=val
            ELSEIF (indx.EQ.6) THEN
              pmas(pycomp(6),1)=val
            ENDIF
          ELSEIF (chblck(1:4).EQ.'NMIX'.OR.chblck(1:4).EQ.'VMIX'.or
     $           .chblck(1:4).EQ.'UMIX'.OR.chblck(1:7).EQ.'STOPMIX'.or
     $           .chblck(1:7).EQ.'SBOTMIX'.OR.chblck(1:7).EQ.'STAUMIX')
     $           THEN
C...NMIX,UMIX,VMIX,STOPMIX,SBOTMIX, and STAUMIX. Mixing.
            im=0
            IF (chblck(5:6).EQ.'IM') im=1
  250       READ(chinl,*) indx1, indx2, val
            IF (chblck(1:1).EQ.'N'.AND.indx1.LE.4.AND.indx2.LE.4) THEN
              IF (im.EQ.0) zmix(indx1,indx2) = val
              IF (im.EQ.1) zmixi(indx1,indx2)= val
              mspc(2)=mspc(2)+1
            ELSEIF (chblck(1:1).EQ.'U') THEN
              IF (im.EQ.0) umix(indx1,indx2) = val
              IF (im.EQ.1) umixi(indx1,indx2)= val
              mspc(3)=mspc(3)+1
            ELSEIF (chblck(1:1).EQ.'V') THEN
              IF (im.EQ.0) vmix(indx1,indx2) = val
              IF (im.EQ.1) vmixi(indx1,indx2)= val
              mspc(4)=mspc(4)+1
            ELSEIF (chblck(1:4).EQ.'STOP'.OR.chblck(1:4).EQ.'SBOT'.or
     $             .chblck(1:4).EQ.'STAU') THEN
              IF (chblck(1:4).EQ.'STOP') THEN
                kfsm=6
                ispc=6
              ELSEIF (chblck(1:4).EQ.'SBOT') THEN
                kfsm=5
                ispc=5
              ELSEIF (chblck(1:4).EQ.'STAU') THEN
                kfsm=15
                ispc=7
              ENDIF
C...Set SFMIX element
              sfmix(kfsm,2*(indx1-1)+indx2)=val
              mspc(ispc)=mspc(ispc)+1
            ENDIF
C...Running parameters
          ELSEIF (chblck(1:4).EQ.'HMIX') THEN
            READ(chblck(8:25),*,err=510) q
            READ(chinl,*) indx, val
            mspc(8)=mspc(8)+1
            IF (indx.EQ.1) THEN
              rmss(4) = val
            ELSE
              merr=1
              mspc(8)=mspc(8)-1
            ENDIF
          ELSEIF (chblck(1:5).EQ.'ALPHA') THEN
            READ(chinl,*,err=520) val
            rmss(18)= val
            mspc(17)=mspc(17)+1
C...Higgs parameters set manually or with FeynHiggs.
            imss(4)=max(2,imss(4))
          ELSEIF (chblck(1:2).EQ.'AU'.OR.chblck(1:2).EQ.'AD'.or
     &           .chblck(1:2).EQ.'AE') THEN
            READ(chblck(9:26),*,err=510) q
            READ(chinl,*) indx1, indx2, val
            IF (chblck(2:2).EQ.'U') THEN
              au(indx1,indx2)=val
              IF (indx1.EQ.3.AND.indx2.EQ.3) rmss(16)=val
              mspc(11)=mspc(11)+1
            ELSEIF (chblck(2:2).EQ.'D') THEN
              ad(indx1,indx2)=val
              IF (indx1.EQ.3.AND.indx2.EQ.3) rmss(15)=val
              mspc(10)=mspc(10)+1
            ELSEIF (chblck(2:2).EQ.'E') THEN
              ae(indx1,indx2)=val
              IF (indx1.EQ.3.AND.indx2.EQ.3) rmss(17)=val
              mspc(12)=mspc(12)+1
            ELSE
              merr=1
            ENDIF
          ELSEIF (chblck(1:5).EQ.'MSOFT') THEN
            IF (mspc(18).EQ.0) THEN
              READ(chblck(9:25),*,err=510) q
              rmsoft(0)=q
            ENDIF
            READ(chinl,*) indx, val
            rmsoft(indx)=val
            mspc(18)=mspc(18)+1
          ELSEIF (chblck(1:5).EQ.'GAUGE') THEN
            merr=8
          ELSEIF (chblck(1:2).EQ.'YU'.OR.chblck(1:2).EQ.'YD'.or
     &           .chblck(1:2).EQ.'YE') THEN
            merr=8
          ELSEIF (chblck(1:6).EQ.'SPINFO') THEN
            READ(chinl(1:6),*) indx
            it=0
            mird=0
  260       it=it+1
            IF (chinl(it:it).EQ.' ') goto 260
C...Don't read index
            IF (chinl(it:it).EQ.char(indx+48).AND.mird.EQ.0) THEN
              mird=1
              goto 260
            ENDIF
            IF (indx.EQ.1) cpro(1)=chinl(it:it+12)
            IF (indx.EQ.2) cver(1)=chinl(it:it+12)
          ELSE
C...  Set unrecognized block flag.
            merr=6
          ENDIF
 
C...DECAY TABLES
C...Read in decay information
        ELSEIF (mupda.EQ.2.AND.merr.EQ.0) THEN
C...Read new decay chanel
          IF(chinl(1:1).EQ.' '.AND.chblck(1:5).EQ.'DECAY') THEN
            ndc=ndc+1
C...Read in branching ratio and number of daughters for this mode.
            READ(chinl(4:50),*,err=480) brat(ndc)
            READ(chinl(4:50),*,err=490) dum, nda
            IF (nda.LE.5) THEN
              IF(ndc.GT.mstu(7)) CALL pyerrm(27,
     &             '(PYSLHA:) Decay data arrays full by KF ='
     $             //chaf(kc,1))
C...If first decay chanel, set decays start point in decay table
              IF(brsum.LE.0d0.AND.brat(ndc).NE.0d0) THEN 
                WRITE(mstu(11),*)
     &              '* (PYSLHA:) Reading in SLHA decay table for ',
     &              chaf(kcrep,1)
C...Set particle parameters (mass set when reading BLOCK MASS above)
                pmas(kc,2)=width
                IF (kf.EQ.25.OR.kf.EQ.35.OR.kf.EQ.36) THEN
                  WRITE(mstu(11),*)
     &                '*  Note: the Pythia gg->h/H/A cross section'//
     &                ' is proportional to the h/H/A->gg width'
                ENDIF
                pmas(kc,3)=0d0
                pmas(kc,4)=paru(3)*1d-12/width
                mwid(kc)=2
                mdcy(kc,1)=1
                mdcy(kc,2)=ndc
                mdcy(kc,3)=0
C...Return ok
                iretrn=0
              ENDIF
C...  Count up number of decay modes for this particle
              mdcy(kc,3)=mdcy(kc,3)+1
C...  Read in decay daughters.
              READ(chinl(4:120),*,err=500) dum,idm, (idc(ida),ida=1,nda)
C...  Flip sign if reading antiparticle decays (if antipartner exists)
              DO 270 ida=1,nda
                IF (kchg(pycomp(idc(ida)),3).NE.0)
     &               idc(ida)=mpsign*idc(ida)
  270         CONTINUE
C...Switch on decay channel, with products ordered in decreasing ABS(KF)
              mdme(ndc,1)=1
              IF (brat(ndc).LE.0d0) mdme(ndc,1)=0
              brsum=brsum+abs(brat(ndc))
              brat(ndc)=abs(brat(ndc))
 274          iflip=0
              DO 277 ida=1,nda-1
                IF (iabs(idc(ida+1)).GT.iabs(idc(ida))) THEN
                  itmp=idc(ida)
                  idc(ida)=idc(ida+1)
                  idc(ida+1)=itmp
                  iflip=iflip+1
                ENDIF
 277          CONTINUE
              IF (iflip.GT.0) goto 274
C              WRITE(MSTU(11),7510) BRAT(NDC), NDA, (IDC(IDA),IDA=1,NDA)
C...Treat as ordinary decay, no fancy stuff.
              mdme(ndc,2)=0
              DO 280 ida=1,5
                IF (ida.LE.nda) THEN
                  kfdp(ndc,ida)=idc(ida)
                ELSE
                  kfdp(ndc,ida)=0
                ENDIF
  280         CONTINUE
            ELSE
              CALL pyerrm(7,'(PYSLHA:) Too many daughters on line '//
     &             chnlin)
              merr=11
              ndc=ndc-1
            ENDIF
          ELSEIF(chinl(1:1).EQ.'+') THEN
            merr=11
          ELSEIF(chblck(1:6).EQ.'DCINFO') THEN
            merr=16
          ELSE
            merr=16
          ENDIF
        ENDIF
C...  Error check.
  290   IF (mod(merr,10).EQ.1.AND.(mupda.EQ.1.OR.mupda.EQ.2)) THEN
          WRITE(mstu(11),*) '* (PYSLHA:) Ignoring line '//chnlin//': '
     &         //chinl(1:40)
          merr=0
        ELSEIF (merr.EQ.6.AND.mupda.EQ.1) THEN
          WRITE(mstu(11),*) '* (PYSLHA:) Ignoring BLOCK '//
     &         chblck(1:inl)//'... on line'//chnlin
        ELSEIF (merr.EQ.8.AND.mupda.EQ.1) THEN
          WRITE(mstu(11),*) '* (PYSLHA:) PYTHIA will not use BLOCK '
     &         //chblck(1:inl)//'... on line'//chnlin
        ELSEIF (merr.EQ.16.AND.mupda.EQ.2.AND.imss(21).EQ.0.AND.
     &         chblck(1:1).NE.'D'.AND.verbos.EQ.1) THEN
          WRITE(mstu(11),*) '* (PYSLHA:) Ignoring BLOCK '//chblck(1:inl)
     &         //'... on line'//chnlin
        ELSEIF (merr.EQ.7.AND.mupda.EQ.1) THEN
          WRITE(mstu(11),*) '* (PYSLHA:) Ignoring extra BLOCK '/
     &         /chblck(1:inl)//'... on line'//chnlin
        ELSEIF (merr.EQ.2.AND.mupda.EQ.1) THEN
          WRITE (chtmp,*) kf
          WRITE(mstu(11),*)
     &         '* (PYSLHA:) Ignoring extra MASS entry for KF='//
     &         chtmp(1:9)//' on line'//chnlin
        ENDIF
C...  End of loop
        goto 160
  300   CONTINUE
C...Set flag that KC codes have been rearranged.
        mstu(20)=0
        verbos=0
 
C...Perform possible tests that new information is consistent.
        IF (mupda.EQ.1) THEN
          mstu23=mstu(23)
          mstu27=mstu(27)
C...Check Z and top masses
          IF (abs(pmas(pycomp(23),1)-91.2d0).GT.1d0) THEN
            WRITE(chtmp,*) pmas(pycomp(23),1)
            CALL pyerrm(19,'(PYSLHA:) note Z boson mass, M ='//chtmp)
          ENDIF
          IF (abs(pmas(pycomp(6),1)-175d0).GT.25d0) THEN
            WRITE(chtmp,*) pmas(pycomp(6),1)
            CALL pyerrm(19,'(PYSLHA:) note top quark mass, M ='
     &           //chtmp//'GeV')
          ENDIF
C...Check masses
          DO 310 isusy=1,37
            kf=kfsusy(isusy)
C...Don't complain about right-handed neutrinos
            IF (kf.EQ.ksusy2+12.OR.kf.EQ.ksusy2+14.OR.kf.EQ.ksusy2
     &           +16) goto 310
C...Only check gravitino in GMSB scenarios
            IF (modsel(1).NE.2.AND.kf.EQ.ksusy1+39) goto 310
            kc=pycomp(kf)
            IF (pmas(kc,1).EQ.0d0) THEN
              WRITE(chtmp,*) kf
              CALL pyerrm(9
     &             ,'(PYSLHA:) No mass information found for KF = '
     &             //chtmp)
            ENDIF
  310     CONTINUE
C...Check mixing matrices (MSSM only)
          IF (imss(13).EQ.0) THEN
            IF (mspc(2).NE.16.AND.mspc(2).NE.32) CALL pyerrm(9
     &           ,'(PYSLHA:) Inconsistent # of elements in NMIX')
            IF (mspc(3).NE.4.AND.mspc(3).NE.8) CALL pyerrm(9
     &           ,'(PYSLHA:) Inconsistent # of elements in UMIX')
            IF (mspc(4).NE.4.AND.mspc(4).NE.8) CALL pyerrm(9
     &           ,'(PYSLHA:) Inconsistent # of elements in VMIX')
            IF (mspc(5).NE.4) CALL pyerrm(9
     &           ,'(PYSLHA:) Inconsistent # of elements in SBOTMIX')
            IF (mspc(6).NE.4) CALL pyerrm(9
     &           ,'(PYSLHA:) Inconsistent # of elements in STOPMIX')
            IF (mspc(7).NE.4) CALL pyerrm(9
     &           ,'(PYSLHA:) Inconsistent # of elements in STAUMIX')
            IF (mspc(8).LT.1) CALL pyerrm(9
     &           ,'(PYSLHA:) Too few elements in HMIX')
            IF (mspc(10).EQ.0) CALL pyerrm(9
     &           ,'(PYSLHA:) Missing A_b trilinear coupling')
            IF (mspc(11).EQ.0) CALL pyerrm(9
     &           ,'(PYSLHA:) Missing A_t trilinear coupling')
            IF (mspc(12).EQ.0) CALL pyerrm(9
     &           ,'(PYSLHA:) Missing A_tau trilinear coupling')
            IF (mspc(17).LT.1) CALL pyerrm(9
     &           ,'(PYSLHA:) Missing Higgs mixing angle alpha')
          ENDIF
C...Check wavefunction normalizations.
C...Sfermions
          DO 320 ispc=5,7
            IF (mspc(ispc).EQ.4) THEN
              kfsm=ispc
              IF (ispc.EQ.7) kfsm=15
              check=abs(sfmix(kfsm,1)*sfmix(kfsm,4)-sfmix(kfsm,2)
     &             *sfmix(kfsm,3))
              IF (abs(1d0-check).GT.1d-3) THEN
                kcsm=pycomp(kfsm)
                CALL pyerrm(17
     &               ,'(PYSLHA:) Non-orthonormal mixing matrix for ~'
     &               //chaf(kcsm,1))
              ENDIF
            ENDIF
  320     CONTINUE
C...Neutralinos + charginos
          DO 340 j=1,4
            cn1=0d0
            cn2=0d0
            cu1=0d0
            cu2=0d0
            cv1=0d0
            cv2=0d0
            DO 330 l=1,4
              cn1=cn1+zmix(j,l)**2
              cn2=cn2+zmix(l,j)**2
              IF (j.LE.2.AND.l.LE.2) THEN
                cu1=cu1+umix(j,l)**2
                cu2=cu2+umix(l,j)**2
                cv1=cv1+vmix(j,l)**2
                cv2=cv2+vmix(l,j)**2
              ENDIF
  330       CONTINUE
C...NMIX normalization
            IF (mspc(2).EQ.16.AND.(abs(1d0-cn1).GT.1d-3.OR.abs(1d0-cn2)
     &           .GT.1d-3).AND.imss(13).EQ.0) THEN
              CALL pyerrm(19,
     &             '(PYSLHA:) NMIX: Inconsistent normalization.')
              WRITE(mstu(11),'(7x,I2,1x,":",2(1x,F7.4))') j, cn1, cn2
            ENDIF
C...UMIX, VMIX normalizations
            IF (mspc(3).EQ.4.OR.mspc(4).EQ.4.AND.imss(13).EQ.0) THEN
              IF (j.LE.2) THEN
                IF (abs(1d0-cu1).GT.1d-3.OR.abs(1d0-cu2).GT.1d-3) THEN
                  CALL pyerrm(19
     &                ,'(PYSLHA:) UMIX: Inconsistent normalization.')
                  WRITE(mstu(11),'(7x,I2,1x,":",2(1x,F6.2))') j, cu1,
     &                 cu2
                ENDIF
                IF (abs(1d0-cv1).GT.1d-3.OR.abs(1d0-cv2).GT.1d-3) THEN
                  CALL pyerrm(19,
     &                '(PYSLHA:) VMIX: Inconsistent normalization.')
                  WRITE(mstu(11),'(7x,I2,1x,":",2(1x,F6.2))') j, cv1,
     &                 cv2
                ENDIF
              ENDIF
            ENDIF
  340     CONTINUE
          IF (mstu(27).EQ.mstu27.AND.mstu(23).EQ.mstu23) THEN
            WRITE(mstu(11),'(1x,"*"/1x,A/1x,"*")')
     &           '*  PYSLHA:  No spectrum inconsistencies were found.'
          ELSE
            WRITE(mstu(11),'(1x,"*"/1x,A/1x,"*",A/1x,"*",A/)')
     &           '* (PYSLHA:) INCONSISTENT SPECTRUM WARNING.'
     &           ,'Warning: one or more (serious)'//
     &           ' inconsistencies were found in the spectrum!!!'
     &           ,'Read the error messages above and check your'//
     &           ' input file.'
          ENDIF
C...Increase precision in Higgs sector using FeynHiggs
          IF (imss(4).EQ.3) THEN
C...FeynHiggs needs MSOFT.
            ierr=0
            IF (mspc(18).EQ.0) THEN
              WRITE(mstu(11),'(1x,"*"/1x,A/)')
     &             '* (PYSLHA:) BLOCK MSOFT not found in SLHA file.'//
     &              ' Cannot call FeynHiggs.'
              ierr=-1
            ELSE
              WRITE(mstu(11),'(1x,/1x,A/)')
     &             '* (PYSLHA:) Now calling FeynHiggs.'
              CALL pyfeyn(ierr)
              IF (ierr.NE.0) imss(4)=2
            ENDIF
          ENDIF
        ELSEIF (mupda.EQ.2.AND.iretrn.EQ.0) THEN
          kf=kforig
          kc=pycomp(kf)
          WRITE(chkf,8300) kf
          IF(min(pmas(kc,1),pmas(kc,2),pmas(kc,3),pmas(kc,1)-pmas(kc,3
     $         ),pmas(kc,4)).LT.0d0.OR.mdcy(kc,3).LT.0.OR.(mdcy(kc,3)
     $         .EQ.0.AND.mdcy(kc,1).GE.1)) CALL pyerrm(17
     $         ,'(PYSLHA:) Mass/width/life/(# channels) wrong for KF='
     $         //chkf)
          brsum=0d0
          bropn=0d0
          DO 360 ida=mdcy(kc,2),mdcy(kc,2)+mdcy(kc,3)-1
            IF(mdme(ida,2).GT.80) goto 360
            kq=kchg(kc,1)
            pms=pmas(kc,1)-pmas(kc,3)-parj(64)
            merr=0
            DO 350 j=1,5
              kp=kfdp(ida,j)
              IF(kp.EQ.0.OR.kp.EQ.81.OR.iabs(kp).EQ.82) THEN
                IF(kp.EQ.81) kq=0
              ELSEIF(pycomp(kp).EQ.0) THEN
                merr=3
              ELSE
                kq=kq-pychge(kp)
                kpc=pycomp(kp)
                pms=pms-pmas(kpc,1)
                IF(mstj(24).GT.0) pms=pms+0.5d0*min(pmas(kpc,2),
     &               pmas(kpc,3))
              ENDIF
  350       CONTINUE
            IF(kq.NE.0) merr=max(2,merr)
            IF(mwid(kc).EQ.0.AND.kf.NE.311.AND.pms.LT.0d0)
     &           merr=max(1,merr)
            IF(merr.EQ.3) CALL pyerrm(17,
     &           '(PYSLHA:) Unknown particle code in decay of KF ='
     $           //chkf)
            IF(merr.EQ.2) CALL pyerrm(17,
     &           '(PYSLHA:) Charge not conserved in decay of KF ='
     $           //chkf)
            IF(merr.EQ.1) CALL pyerrm(7,
     &           '(PYSLHA:) Kinematically unallowed decay of KF ='
     $           //chkf)
            brsum=brsum+brat(ida)
            IF (mdme(ida,1).GT.0) bropn=bropn+brat(ida)
  360     CONTINUE
C...Check branching ratio sum.
          IF (bropn.LE.0d0) THEN
C...If zero, set stable. 
             WRITE(chtmp,8500) bropn
             CALL pyerrm(7
     &            ,"(PYSLHA:) Effective BR sum for KF="//chkf//' is '//
     &            chtmp(9:16)//'. Changed to stable.')
             pmas(kc,2)=1d-6
             mwid(kc)=0
C...If BR's > 1, rescale.
          ELSEIF (brsum.GT.(1d0+1d-6)) THEN
             WRITE(chtmp,8500) brsum
             CALL pyerrm(7
     &            ,"(PYSLHA:) Forced rescaling of BR's for KF="//chkf//
     &            ' ; sum was'//chtmp(9:16)//'.')
             fac=1d0/brsum
             DO 370 ida=mdcy(kc,2),mdcy(kc,2)+mdcy(kc,3)-1
                IF(mdme(ida,2).GT.80) goto 370
                brat(ida)=fac*brat(ida)
 370         CONTINUE
          ELSEIF (brsum.LT.(1d0-1d-6)) THEN
C...If BR's < 1, insert dummy mode for proper cross section rescaling.
             WRITE(chtmp,8500) brsum
             CALL pyerrm(7
     &            ,"(PYSLHA:) Sum of BR's for KF="//chkf//' is '//
     &            chtmp(9:16)//'. Dummy mode will be inserted.')
C...  Insert dummy mode
             mdcy(kc,3)=mdcy(kc,3)+1
             ida=mdcy(kc,2)+mdcy(kc,3)-1
             brat(ida)=1d0-brsum
             kfdp(ida,1)=0
             kfdp(ida,2)=0
             kfdp(ida,3)=0
             kfdp(ida,4)=0
             kfdp(ida,5)=0
             mdme(ida,1)=0
             brsum=1d0
          ENDIF
       ENDIF
 
C...WRITE SPECTRUM ON SLHA FILE
      ELSEIF(mupda.EQ.3) THEN
C...If SPYTHIA or ISASUSY runtime was called for SUGRA, update PARMIN.
        IF (imss(1).EQ.2.OR.imss(1).EQ.12) THEN
          modsel(1)=1
          parmin(1)=rmss(8)
          parmin(2)=rmss(1)
          parmin(3)=rmss(5)
          parmin(4)=sign(1d0,rmss(4))
          parmin(5)=rmss(36)
        ENDIF
C...Write spectrum
        WRITE(lfn,7000) 'SLHA MSSM spectrum'
        WRITE(lfn,7000) 'Pythia 6.4: T. Sjostrand, S. Mrenna,'
     &    // ' P. Skands.'
        WRITE(lfn,7010) 'MODSEL',  'Model selection'
        WRITE(lfn,7110) 1, modsel(1)
        WRITE(lfn,7010) 'MINPAR', 'Parameters for minimal model.'
        IF (modsel(1).EQ.1) THEN
          WRITE(lfn,7210) 1, parmin(1), 'm0'
          WRITE(lfn,7210) 2, parmin(2), 'm12'
          WRITE(lfn,7210) 3, parmin(3), 'tan(beta)'
          WRITE(lfn,7210) 4, parmin(4), 'sign(mu)'
          WRITE(lfn,7210) 5, parmin(5), 'a0'
        ELSEIF(modsel(2).EQ.2) THEN
          WRITE(lfn,7210) 1, parmin(1), 'Lambda'
          WRITE(lfn,7210) 2, parmin(2), 'M'
          WRITE(lfn,7210) 3, parmin(3), 'tan(beta)'
          WRITE(lfn,7210) 4, parmin(4), 'sign(mu)'
          WRITE(lfn,7210) 5, parmin(5), 'N5'
          WRITE(lfn,7210) 6, parmin(6), 'c_grav'
        ENDIF
        WRITE(lfn,7000) ' '
        WRITE(lfn,7010) 'MASS', 'Mass spectrum'
        DO 380 i=1,36
          kf=kfsusy(i)
          kc=pycomp(kf)
          IF (kf.EQ.1000039.AND.modsel(1).NE.2) goto 380
          kfsm=kf-ksusy1
          IF (kfsm.GE.22.AND.kfsm.LE.37) THEN
            IF (kfsm.EQ.22)  WRITE(lfn,7220) kf, smz(1), chaf(kc,1)
            IF (kfsm.EQ.23)  WRITE(lfn,7220) kf, smz(2), chaf(kc,1)
            IF (kfsm.EQ.25)  WRITE(lfn,7220) kf, smz(3), chaf(kc,1)
            IF (kfsm.EQ.35)  WRITE(lfn,7220) kf, smz(4), chaf(kc,1)
            IF (kfsm.EQ.24)  WRITE(lfn,7220) kf, smw(1), chaf(kc,1)
            IF (kfsm.EQ.37)  WRITE(lfn,7220) kf, smw(2), chaf(kc,1)
          ELSE
            WRITE(lfn,7220) kf, pmas(kc,1), chaf(kc,1)
          ENDIF
  380   CONTINUE
C...SUSY scale
        rmsusy=sqrt(pmas(pycomp(ksusy1+6),1)*pmas(pycomp(ksusy2+6),1))
        WRITE(lfn,7020) 'HMIX',rmsusy,'Higgs parameters'
        WRITE(lfn,7210) 1, rmss(4),'mu'
        WRITE(lfn,7010) 'ALPHA',' '
        WRITE(lfn,7210) 1, rmss(18), 'alpha'
        WRITE(lfn,7020) 'AU',rmsusy
        WRITE(lfn,7410) 3, 3, rmss(16), 'A_t'
        WRITE(lfn,7020) 'AD',rmsusy
        WRITE(lfn,7410) 3, 3, rmss(15), 'A_b'
        WRITE(lfn,7020) 'AE',rmsusy
        WRITE(lfn,7410) 3, 3, rmss(17), 'A_tau'
        WRITE(lfn,7010) 'STOPMIX','~t mixing matrix'
        WRITE(lfn,7410) 1, 1, sfmix(6,1)
        WRITE(lfn,7410) 1, 2, sfmix(6,2)
        WRITE(lfn,7410) 2, 1, sfmix(6,3)
        WRITE(lfn,7410) 2, 2, sfmix(6,4)
        WRITE(lfn,7010) 'SBOTMIX','~b mixing matrix'
        WRITE(lfn,7410) 1, 1, sfmix(5,1)
        WRITE(lfn,7410) 1, 2, sfmix(5,2)
        WRITE(lfn,7410) 2, 1, sfmix(5,3)
        WRITE(lfn,7410) 2, 2, sfmix(5,4)
        WRITE(lfn,7010) 'STAUMIX','~tau mixing matrix'
        WRITE(lfn,7410) 1, 1, sfmix(15,1)
        WRITE(lfn,7410) 1, 2, sfmix(15,2)
        WRITE(lfn,7410) 2, 1, sfmix(15,3)
        WRITE(lfn,7410) 2, 2, sfmix(15,4)
        WRITE(lfn,7010) 'NMIX','~chi0 mixing matrix'
        DO 400 i1=1,4
          DO 390 i2=1,4
            WRITE(lfn,7410) i1, i2, zmix(i1,i2)
  390     CONTINUE
  400   CONTINUE
        WRITE(lfn,7010) 'UMIX','~chi^+ U mixing matrix'
        DO 420 i1=1,2
          DO 410 i2=1,2
            WRITE(lfn,7410) i1, i2, umix(i1,i2)
  410     CONTINUE
  420   CONTINUE
        WRITE(lfn,7010) 'VMIX','~chi^+ V mixing matrix'
        DO 440 i1=1,2
          DO 430 i2=1,2
            WRITE(lfn,7410) i1, i2, vmix(i1,i2)
  430     CONTINUE
  440   CONTINUE
        WRITE(lfn,7010) 'SPINFO'
        IF (imss(1).EQ.2) THEN
          cpro(1)='PYTHIA'
          cver(1)='6.4'
        ELSEIF (imss(1).EQ.12) THEN
          isaver=visaje()
          cpro(1)='ISASUSY'
          cver(1)=isaver(1:12)
        ENDIF
        WRITE(lfn,7310) 1, cpro(1), 'Spectrum Calculator'
        WRITE(lfn,7310) 2, cver(1), 'Version number'
      ENDIF
 
C...Print user information about spectrum
      IF (mupda.EQ.1.OR.mupda.EQ.3) THEN
        IF (cpro(mod(mupda,2)).NE.' '.AND.cver(mod(mupda,2)).NE.' ')
     &       WRITE(mstu(11),5030) cpro(1), cver(1)
        IF (imss(4).EQ.3) WRITE(mstu(11),5040)
        IF (mupda.EQ.1) THEN
          WRITE(mstu(11),5020) lfn
        ELSE
          WRITE(mstu(11),5010) lfn
        ENDIF
 
        WRITE(mstu(11),5400)
        WRITE(mstu(11),5500) 'Pole masses'
        WRITE(mstu(11),5700) (rmfun(ksusy1+ip),ip=1,6)
     $       ,(rmfun(ksusy2+ip),ip=1,6)
        WRITE(mstu(11),5800) (rmfun(ksusy1+ip),ip=11,16)
     $       ,(rmfun(ksusy2+ip),ip=11,16)
        IF (imss(13).EQ.0) THEN
          WRITE(mstu(11),5900) rmfun(ksusy1+21),rmfun(ksusy1+22)
     $         ,rmfun(ksusy1+23),rmfun(ksusy1+25),rmfun(ksusy1+35),
     $         rmfun(ksusy1+24),rmfun(ksusy1+37)
          WRITE(mstu(11),6000) chaf(25,1),chaf(35,1),chaf(36,1),
     &         chaf(37,1), ' ', ' ',' ',' ',
     &         rmfun(25), rmfun(35), rmfun(36), rmfun(37)
        ELSEIF (imss(13).EQ.1) THEN
          kf1=ksusy1+21
          kf2=ksusy1+22
          kf3=ksusy1+23
          kf4=ksusy1+25
          kf5=ksusy1+35
          kf6=ksusy1+45
          kf7=ksusy1+24
          kf8=ksusy1+37
          WRITE(mstu(11),6000) chaf(pycomp(kf1),1),chaf(pycomp(kf2),1),
     &         chaf(pycomp(kf3),1),chaf(pycomp(kf4),1),
     &         chaf(pycomp(kf5),1),chaf(pycomp(kf6),1),
     &         chaf(pycomp(kf7),1),chaf(pycomp(kf8),1),
     &         rmfun(kf1),rmfun(kf2),rmfun(kf3),rmfun(kf4),
     &         rmfun(kf5),rmfun(kf6),rmfun(kf7),rmfun(kf8)
          WRITE(mstu(11),6000) chaf(25,1), chaf(35,1), chaf(45,1),
     &         chaf(36,1), chaf(46,1), chaf(37,1),' ',' ',
     &         rmfun(25), rmfun(35), rmfun(45), rmfun(36), rmfun(46),
     &         rmfun(37)
        ENDIF
        WRITE(mstu(11),5400)
        WRITE(mstu(11),5500) 'Mixing structure'
        WRITE(mstu(11),6100) ((zmix(i,j), j=1,4),i=1,4)
        WRITE(mstu(11),6200) (umix(1,j), j=1,2),(vmix(1,j),j=1,2)
     &       ,(umix(2,j), j=1,2),(vmix(2,j),j=1,2)
        WRITE(mstu(11),6300) (sfmix(5,j), j=1,2),(sfmix(6,j),j=1,2)
     &       ,(sfmix(15,j), j=1,2),(sfmix(5,j),j=3,4),(sfmix(6,j), j=3,4
     &       ),(sfmix(15,j),j=3,4)
        WRITE(mstu(11),5400)
        WRITE(mstu(11),5500) 'Couplings'
        WRITE(mstu(11),6400) rmss(15),rmss(16),rmss(17)
        WRITE(mstu(11),6450) rmss(18), rmss(5), rmss(4)
        WRITE(mstu(11),5400)
        WRITE(mstu(11),6500)
 
      ENDIF
 
C...Only rewind when reading
      IF (mupda.LE.2.OR.mupda.EQ.5) rewind(lfn)
 
 9999 RETURN
 
C...Serious error catching
  460 write(*,*) '* (PYSLHA:) read BLOCK error on line',nline
      write(*,*) chinl(1:80)
      CALL pystop(106)
  470 WRITE(*,*) '* (PYSLHA:) read DECAY error on line',nline
      WRITE(*,*) chinl(1:72)
      CALL pystop(106)
  480 WRITE(*,*) '* (PYSLHA:) read BR error on line',nline
      WRITE(*,*) chinl(1:80)
      CALL pystop(106)
  490 WRITE(*,*) '* (PYSLHA:) read NDA error on line',nline
      WRITE(*,*) chinl(1:80)
      CALL pystop(106)
  500 WRITE(*,*) '* (PYSLHA:) decay daughter read error on line',nline
      WRITE(*,*) chinl(1:80)
  510 WRITE(*,*) '* (PYSLHA:) read Q error in BLOCK ',chblck
      CALL pystop(106)
  520 WRITE(*,*) '* (PYSLHA:) read error in line ',nline,':'
      WRITE(*,*) chinl(1:80)
      CALL pystop(106)
 
 8300 FORMAT(i9)
 8500 FORMAT(f16.5)
 
C...Formats for user information printout.
 5000 FORMAT(1x,15('*'),1x,'PYSLHA v1.09: SUSY/BSM SPECTRUM '
     &     ,'INTERFACE',1x,15('*')/1x,'*',2x
     &     ,'PYSLHA:  Last Change',1x,a,1x,'-',1x,'P.Z. Skands')
 5010 FORMAT(1x,'*',3x,'Wrote spectrum file on unit: ',i3)
 5020 FORMAT(1x,'*',3x,'Read spectrum file on unit: ',i3)
 5030 FORMAT(1x,'*',3x,'Spectrum Calculator was: ',a,' version ',a)
 5040 FORMAT(1x,'*',3x,'Higgs sector corrected with FeynHiggs')
 5100 FORMAT(1x,'*',1x,'Model parameters:'/1x,'*',1x,'----------------')
 5200 FORMAT(1x,'*',1x,3x,'M_0',6x,'M_1/2',5x,'A_0',3x,'Tan(beta)',
     &     3x,'Sgn(mu)',3x,'M_t'/1x,'*',1x,4(f8.2,1x),i8,2x,f8.2)
 5300 FORMAT(1x,'*'/1x,'*',1x,'Model spectrum :'/1x,'*',1x
     &     ,'----------------')
 5400 FORMAT(1x,'*',1x,a)
 5500 FORMAT(1x,'*',1x,a,':')
 5600 FORMAT(1x,'*',2x,2x,'M_GUT',2x,2x,'g_GUT',2x,1x,'alpha_GUT'/
     &       1x,'*',2x,1p,2(1x,e8.2),2x,e8.2)
 5700 FORMAT(1x,'*',4x,4x,'~d',2x,1x,4x,'~u',2x,1x,4x,'~s',2x,1x,
     &     4x,'~c',2x,1x,1x,'~b(12)',1x,1x,1x,'~t(12)'/1x,'*',2x,'L',1x
     &     ,6(f8.2,1x)/1x,'*',2x,'R',1x,6(f8.2,1x))
 5800 FORMAT(1x,'*'/1x,'*',4x,4x,'~e',2x,1x,2x,'~nu_e',2x,1x,3x,'~mu',2x
     &     ,1x,1x,'~nu_mu',1x,1x,'~tau(12)',1x,1x,'~nu_tau'/1x,'*',2x
     &     ,'L',1x,6(f8.2,1x)/1x,'*',2x,'R',1x,6(f8.2,1x))
 5900 FORMAT(1x,'*'/1x,'*',4x,4x,'~g',2x,1x,1x,'~chi_10',1x,1x,'~chi_20'
     &     ,1x,1x,'~chi_30',1x,1x,'~chi_40',1x,1x,'~chi_1+',1x
     &     ,1x,'~chi_2+'/1x,'*',3x,1x,7(f8.2,1x))
 6000 FORMAT(1x,'*'/1x,'*',3x,1x,8(1x,a7,1x)/1x,'*',3x,1x,8(f8.2,1x))
 6100 FORMAT(1x,'*',11x,'|',3x,'~B',3x,'|',2x,'~W_3',2x,'|',2x
     &     ,'~H_1',2x,'|',2x,'~H_2',2x,'|'/1x,'*',3x,'~chi_10',1x,4('|'
     &     ,1x,f6.3,1x),'|'/1x,'*',3x,'~chi_20',1x,4('|'
     &     ,1x,f6.3,1x),'|'/1x,'*',3x,'~chi_30',1x,4('|'
     &     ,1x,f6.3,1x),'|'/1x,'*',3x,'~chi_40',1x,4('|'
     &     ,1x,f6.3,1x),'|')
 6200 FORMAT(1x,'*'/1x,'*',6x,'L',4x,'|',3x,'~W',3x,'|',3x,'~H',3x,'|'
     &     ,12x,'R',4x,'|',3x,'~W',3x,'|',3x,'~H',3x,'|'/1x,'*',3x
     &     ,'~chi_1+',1x,2('|',1x,f6.3,1x),'|',9x,'~chi_1+',1x,2('|',1x
     &     ,f6.3,1x),'|'/1x,'*',3x,'~chi_2+',1x,2('|',1x,f6.3,1x),'|',9x
     &     ,'~chi_2+',1x,2('|',1x,f6.3,1x),'|')
 6300 FORMAT(1x,'*'/1x,'*',8x,'|',2x,'~b_L',2x,'|',2x,'~b_R',2x,'|',8x
     &     ,'|',2x,'~t_L',2x,'|',2x,'~t_R',2x,'|',10x
     &     ,'|',1x,'~tau_L',1x,'|',1x,'~tau_R',1x,'|'/
     &     1x,'*',3x,'~b_1',1x,2('|',1x,f6.3,1x),'|',3x,'~t_1',1x,2('|'
     &     ,1x,f6.3,1x),'|',3x,'~tau_1',1x,2('|',1x,f6.3,1x),'|'/
     &     1x,'*',3x,'~b_2',1x,2('|',1x,f6.3,1x),'|',3x,'~t_2',1x,2('|'
     &     ,1x,f6.3,1x),'|',3x,'~tau_2',1x,2('|',1x,f6.3,1x),'|')
 6400 FORMAT(1x,'*',3x,'  A_b = ',f8.2,4x,'      A_t = ',f8.2,4x
     &     ,'A_tau = ',f8.2)
 6450 FORMAT(1x,'*',3x,'alpha = ',f8.2,4x,'tan(beta) = ',f8.2,4x
     &     ,'   mu = ',f8.2)
 6500 FORMAT(1x,32('*'),1x,'END OF PYSLHA',1x,31('*'))
 
C...Format to use for comments
 7000 FORMAT('# ',a)
C...Format to use for block statements
 7010 FORMAT('Block',1x,a,3x,'#',1x,a)
 7020 FORMAT('Block',1x,a,1x,'Q=',1p,e16.8,0p,3x,'#',1x,a)
C...Indexed Int
 7110 FORMAT(1x,i4,1x,i4,3x,'#')
C...Non-Indexed Double
 7200 FORMAT(9x,1p,e16.8,0p,3x,'#',1x,a)
C...Indexed Double
 7210 FORMAT(1x,i4,3x,1p,e16.8,0p,3x,'#',1x,a)
C...Long Indexed Double (PDG + double)
 7220 FORMAT(1x,i9,3x,1p,e16.8,0p,3x,'#',1x,a)
C...Indexed Char(12)
 7310 FORMAT(1x,i4,3x,a12,3x,'#',1x,a)
C...Single matrix
 7410 FORMAT(1x,i2,1x,i2,3x,1p,e16.8,0p,3x,'#',1x,a)
C...Double Matrix
 7420 FORMAT(1x,i2,1x,i2,3x,1p,e16.8,3x,e16.8,0p,3x,'#',1x,a)
C...Write Decay Table
 7500 FORMAT('Decay',1x,i9,1x,'WIDTH=',1p,e16.8,0p,3x,'#',1x,a)
 7510 FORMAT(4x,1p,e16.8,0p,3x,i2,3x,'IDA=',1x,5(1x,i9),3x,'#',1x,a)
 
      END