]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - ISAJET/isasusy/ssmhn.F
First commit.
[u/mrichter/AliRoot.git] / ISAJET / isasusy / ssmhn.F
diff --git a/ISAJET/isasusy/ssmhn.F b/ISAJET/isasusy/ssmhn.F
new file mode 100644 (file)
index 0000000..5801f82
--- /dev/null
@@ -0,0 +1,454 @@
+#include "isajet/pilot.h"
+      SUBROUTINE SSMHN(MHLNEG)
+C-----------------------------------------------------------------------
+C
+C          Calculate HL, HH masses and ALFAH 
+C          (scalar Higgs mixing angle) using radiative
+C          corrections calculated by M. Bisset
+C          and save results in /SSPAR/.
+C
+C          Both top and bottom couplings are now 
+C          included.  Non-degenerate mixed squark
+C          masses and A-terms are also included.
+C          The D-terms from the squark mass matrix
+C          (terms prop. to g**2 * Yukawa coupling)
+C          are included as an option: 
+C                 INRAD = 1 ==> D-TERMS ON
+C                 INRAD = 2 ==> D-TERMS OFF    .
+C
+C         10/18/93 D-terms are now turned on.
+C                     INRAD = 1 
+C
+C         There is an arbitrary mass scale that must
+C         chosen to avoid dimensionful logarithms.
+C         The choice does not matter if D-terms are
+C         not included, but it does matter if D-terms
+C         are included. 
+C     
+C         Arbitrary mass scale updated to 
+C              QQQ = HIGFRZ = SQRT(AMTLSS*AMTRSS)
+C         with running masses to include dominant 2-loop 
+C         effects. 12/10/96 H. Baer
+C
+C         It is assumed that the A-terms are real.
+C         Complex A-terms are allowed 
+C         (unless RTT=0 or RBB=0 --see below) in 
+C         this subroutine, but the imaginary parts
+C         are now set to zero. 
+C
+C-----------------------------------------------------------------------
+#if defined(CERNLIB_IMPNONE)
+      IMPLICIT NONE
+#endif
+#include "isajet/sslun.inc"
+#include "isajet/sssm.inc"
+#include "isajet/sspar.inc"
+C
+      REAL PI,PI2,SR2,G2,GP2,GGP,GG1,GG2
+      REAL TANB,COTB,COSB,SINB,BE
+      REAL SINB2,COSB2,COS2B,SIN2B
+      REAL V2,VP2,V,VP,VVP,VPVM,VVPP
+      REAL MT2,MB2,FT2,FB2,MW2,ZAP,QQQ2
+      REAL EP,EP2,RR,MHP2
+      REAL ATI,ABI,ATR,ABR,AT2,AB2
+      REAL TLRM,BLRM,TLRP,BLRP
+      REAL MST1SQ,MST2SQ,MSB1SQ,MSB2SQ
+
+      REAL RTT,TTT1,TEMPT,TM1BT
+      REAL TEMPS,T1RD,T2RD,T1RPD,T2RPD
+      REAL CT1,A1,A2,T1RR,T2RR
+      REAL CT5,A5,A6,T1RPRP,T2RPRP
+      REAL A9,T1RRP,T2RRP
+      REAL TEMPSQ,DT1,DT2,VRRT,VRPRPT,VRRPT
+      REAL ALPHAT,LAT
+C
+      REAL RBB,BBB1,TEMPB,TM1BB
+      REAL B1RD,B2RD,B1RPD,B2RPD
+      REAL CB3,A3,A4,B1RR,B2RR
+      REAL CB7,A7,A8,B1RPRP,B2RPRP
+      REAL A10,B1RRP,B2RRP
+      REAL DB1,DB2,VRRB,VRPRPB,VRRPB
+      REAL ALPHAB,LAB
+C
+      REAL DVRR,DVRPRP,DVRRP,TEMPH
+      REAL MHL2,MHH2,TRACEM,TPAL,TANAH
+      REAL ASMB,MBMB,MBQ,ASMT,MTMT,MTQ,SUALFS,HIGFRZ
+      DOUBLE PRECISION SSMQCD
+      INTEGER INRAD,MHLNEG
+C
+      MHLNEG=0
+      PI=4.*ATAN(1.)
+      PI2 = PI**2
+      SR2=SQRT(2.)
+      G2=4.*PI*ALFAEM/SN2THW
+      GP2=G2*SN2THW/(1.-SN2THW)
+      HIGFRZ=SQRT(AMTLSS*AMTRSS)
+      ASMB=SUALFS(AMBT**2,.36,AMTP,3)
+      MBMB=AMBT*(1.-4*ASMB/3./PI)
+      MBQ=SSMQCD(DBLE(MBMB),DBLE(HIGFRZ))
+      ASMT=SUALFS(AMTP**2,.36,AMTP,3)
+      MTMT=AMTP/(1.+4*ASMT/3./PI+(16.11-1.04*(5.-6.63/AMTP))*
+     $(ASMT/PI)**2)
+      MTQ=SSMQCD(DBLE(MTMT),DBLE(HIGFRZ))
+      MT2=MTQ**2
+      MB2=MBQ**2
+      MW2=AMW**2
+      EP=TWOM1
+      EP2=EP**2
+      RR=RV2V1
+      MHP2=AMHA**2
+      TANB=1.0/RR
+      COTB=RR
+      BE=ATAN(1./RV2V1)
+      SINB=SIN(BE)
+      COSB=COS(BE)
+      SINB2=SINB**2
+      COSB2=COSB**2
+      SIN2B=SIN(2.*BE)
+      COS2B=COS(2.*BE)
+      V2=2.0*MW2*SINB2/G2
+      VP2=2.0*MW2*COSB2/G2
+      V=SQRT(V2)
+      VP=SQRT(VP2)
+      VVP=SQRT(V2*VP2)
+      VPVM=VP2-V2
+      GGP=G2+GP2
+      GG1=G2-5.0*GP2/3.0
+      GG2=G2-GP2/3.0
+      VVPP=2.0*AMZ**2/GGP
+      FT2=MT2/V2
+      FB2=MB2/VP2
+C
+      TLRM=AMTLSS**2-AMTRSS**2
+      BLRM=AMBLSS**2-AMBRSS**2 
+      TLRP=AMTLSS**2+AMTRSS**2
+      BLRP=AMBLSS**2+AMBRSS**2 
+C
+C          Higgs mass matrix
+C
+C     (AAT and AAB are also assumed to be real)
+C          
+      ATR=AAT
+      ABR=AAB
+      ATI=0.0
+      ABI=0.0
+      AT2=ATR**2+ATI**2
+      AB2=ABR**2+ABI**2
+C
+      MST1SQ=AMT1SS**2
+      MST2SQ=AMT2SS**2
+      MSB1SQ=AMB1SS**2
+      MSB2SQ=AMB2SS**2
+      INRAD=1
+      QQQ2=HIGFRZ**2
+C
+      ZAP = 1.0
+C
+C                  STOP TERMS
+C          
+      RTT=(TLRM+VPVM*ZAP*GG1/4.0)**2
+     $      +4.0*MT2*(EP*COTB+ATR)**2+4.0*MT2*ATI**2
+      RTT=SQRT(RTT)
+C
+C     calculate 2M1*B term
+C
+      TTT1=0.5*TLRP+MT2+VPVM*ZAP*GGP/8.0
+      IF(RTT.NE.0.0) THEN
+        TEMPT=4.0*EP*FT2*VVP*ATI**2/(RTT**2)
+        TM1BT=-2.0*FT2*(TEMPT+ATR)*TTT1
+     $               *LOG(MST2SQ/MST1SQ)/RTT
+        TM1BT=TM1BT-FT2*ATR
+     $               *LOG(MST1SQ*MST2SQ/QQQ2/QQQ2)
+        TM1BT=TM1BT+FT2*(2.0*TEMPT-ATR)
+        TM1BT=3.0*EP*TM1BT/32.0/PI2
+C
+C        calculate first derivatives w.r.t H_R
+C           divided by sqrt(2) * v
+C        
+         TEMPS=-ZAP*GG1*(TLRM+ZAP*GG1*VPVM/4.0)/2.0 
+         TEMPS=TEMPS+4.0*FT2*(AT2+EP*COTB*ATR)
+         TEMPS=TEMPS/RTT/4.0 
+         T1RD=FT2-ZAP*GGP/8.0-TEMPS
+         T2RD=FT2-ZAP*GGP/8.0+TEMPS
+C
+C        calculate first derivatives w.r.t H_R'
+C           divided by sqrt(2) * v'
+C        
+         TEMPS=ZAP*GG1*(TLRM+ZAP*GG1*VPVM/4.0)/2.0
+         TEMPS=TEMPS+4.0*FT2*EP*(EP+TANB*ATR)
+         TEMPS=TEMPS/RTT/4.0
+         T1RPD=ZAP*GGP/8.0-TEMPS
+         T2RPD=ZAP*GGP/8.0+TEMPS
+C
+C        calculate second derivatives w.r.t. H_R
+C
+         CT1=-V*ZAP*GG1*(TLRM+ZAP*GG1*VPVM/4.0)/SR2
+         CT1=CT1+4.0*SR2*FT2*V*(EP*COTB*ATR+AT2)
+         A1=-CT1**2/(RTT**3)/8.0
+         A2=-ZAP*GG1*(TLRM+ZAP*GG1*VPVM/4.0)/2.0
+         A2=A2+V2*ZAP*GG1**2/4.0+4.0*FT2*AT2
+         A2=A2/RTT/4.0
+         T1RR=FT2-ZAP*GGP/8.0-A1-A2
+         T2RR=FT2-ZAP*GGP/8.0+A1+A2
+C
+C        calculate second derivatives w.r.t. H_R'
+C
+         CT5=VP*ZAP*GG1*(TLRM+ZAP*GG1*VPVM/4.0)/SR2
+         CT5=CT5+4.0*SR2*FT2*VP*EP*(EP+TANB*ATR)
+         A5=-CT5**2/(RTT**3)/8.0
+         A6=ZAP*GG1*(TLRM+ZAP*GG1*VPVM/4.0)/2.0 
+         A6=A6+VP2*ZAP*GG1**2/4.0+4.0*FT2*EP2
+         A6=A6/RTT/4.0
+         T1RPRP=ZAP*GGP/8.0-A5-A6
+         T2RPRP=ZAP*GGP/8.0+A5+A6
+C
+C        calculate second derivatives w.r.t. H_R and H_R'
+C
+         A9=-VVP*ZAP*(GG1**2)/4.0+4.0*FT2*EP*ATR
+         A9=A9/RTT/4.0
+         T1RRP=CT1*CT5/(RTT**3)/8.0-A9
+         T2RRP=-CT1*CT5/(RTT**3)/8.0+A9
+C
+C        calculate D^2 V / D^2 H_R
+C
+         TEMPSQ=MST1SQ*(T1RR-T1RD)
+         DT1=2.0*(2.0*V2*T1RD**2+TEMPSQ)*LOG(MST1SQ/QQQ2)
+         DT1=DT1+6.0*V2*T1RD**2+TEMPSQ
+         TEMPSQ=MST2SQ*(T2RR-T2RD)
+         DT2=2.0*(2.0*V2*T2RD**2+TEMPSQ)*LOG(MST2SQ/QQQ2)
+         DT2=DT2+6.0*V2*T2RD**2+TEMPSQ
+         VRRT=DT1+DT2-8.0*FT2*MT2*LOG(MT2/QQQ2)-12.0*FT2*MT2
+         VRRT=-TM1BT*COTB+3.0*VRRT/32.0/PI2
+C
+C        calculate D^2 V / D^2 H'_R
+C
+         TEMPSQ=MST1SQ*(T1RPRP-T1RPD)
+         DT1=2.0*(2.0*VP2*T1RPD**2+TEMPSQ)*LOG(MST1SQ/QQQ2)
+         DT1=DT1+6.0*VP2*T1RPD**2+TEMPSQ
+         TEMPSQ=MST2SQ*(T2RPRP-T2RPD)
+         DT2=2.0*(2.0*VP2*T2RPD**2+TEMPSQ)*LOG(MST2SQ/QQQ2)
+         DT2=DT2+6.0*VP2*T2RPD**2+TEMPSQ
+         VRPRPT=-TM1BT*TANB+3.0*(DT1+DT2)/32.0/PI2
+C
+C        calculate D^2 V / D^H_R  D^H_R'
+C
+         DT1=2.0*VVP*T1RD*T1RPD+MST1SQ*T1RRP
+         DT1=2.0*DT1*LOG(MST1SQ/QQQ2)
+         DT1=DT1+6.0*VVP*T1RD*T1RPD+MST1SQ*T1RRP
+         DT2=2.0*VVP*T2RD*T2RPD+MST2SQ*T2RRP
+         DT2=2.0*DT2*LOG(MST2SQ/QQQ2)
+         DT2=DT2+6.0*VVP*T2RD*T2RPD+MST2SQ*T2RRP
+         VRRPT=TM1BT+3.0*(DT1+DT2)/32.0/PI2
+C
+      ELSE IF(RTT.EQ.0.0) THEN
+C
+         ALPHAT=TLRP/2.0+MT2+ZAP*GGP*VPVM/8.0
+         LAT=2.0*LOG(ALPHAT/QQQ2)+3.0
+C
+C        calculate D^2 V / D^2 H_R
+C
+         VRRT=V2*(GGP**2+GG1**2)/16.0-MT2*GGP
+         VRRT=ZAP*VRRT*LAT+8.0*FT2*MT2*LOG(ALPHAT/MT2)
+         VRRT=3.0*VRRT/32.0/PI2
+C
+C        calculate D^2 V / D^2 H_R'
+C
+         VRPRPT=ZAP*VP2*(GGP**2+GG1**2)/16.0
+         VRPRPT=3.0*(VRPRPT*LAT)/32.0/PI2
+C
+C        calculate D^2 V / D^H_R D^H_R'
+C
+         VRRPT=FT2*GGP-(GGP**2+GG1**2)/8.0
+         VRRPT=ZAP*VVP*VRRPT*LAT/2.0
+         VRRPT=3.0*VRRPT/32.0/PI2
+C
+C
+      ENDIF
+C
+C     SBOTTOM TERMS
+C
+      RBB=(BLRM-VPVM*ZAP*GG2/4.0)**2
+     $      +4.0*MB2*(EP*TANB+ABR)**2+4.0*MB2*ABI**2
+      RBB=SQRT(RBB)
+C      IF(RBB.EQ.0.0.AND.ABI.NE.0.0) THEN
+C        WRITE(6,*) 'RBB=0, ABI NOT 0'
+C        WRITE(6,*) 'ERROR: THIS CASE NOT COVERED YET'
+C        GO TO 1000
+C      ENDIF
+C
+      IF(RBB.NE.0.0) THEN
+C
+C     calculate 2M1*B term
+C
+        BBB1=0.5*BLRP+MB2-VPVM*ZAP*GGP/8.0
+        TEMPB=4.0*EP*FB2*VVP*ABI**2/(RBB**2)       
+        TM1BB=-2.0*FB2*(TEMPB+ABR)*BBB1
+     $               *LOG(MSB2SQ/MSB1SQ)/RBB
+        TM1BB=TM1BB-FB2*ABR
+     $               *LOG(MSB1SQ*MSB2SQ/QQQ2/QQQ2)
+        TM1BB=TM1BB+FB2*(2.0*TEMPB-ABR)
+        TM1BB=3.0*EP*TM1BB/32.0/PI2
+C
+C        calculate first derivatives w.r.t H_R
+C           divided by sqrt(2) * v
+C        
+        TEMPS=ZAP*GG2*(BLRM-ZAP*GG2*VPVM/4.0)/2.0
+        TEMPS=TEMPS+4.0*FB2*EP*(EP+COTB*ABR)
+        TEMPS=TEMPS/RBB/4.0
+        B1RD=ZAP*GGP/8.0-TEMPS
+        B2RD=ZAP*GGP/8.0+TEMPS
+
+C        calculate first derivatives w.r.t H_R'
+C           divided by sqrt(2) * v'
+C
+        TEMPS=-ZAP*GG2*(BLRM-ZAP*GG2*VPVM/4.0)/2.0
+        TEMPS=TEMPS+4.0*FB2*(AB2+EP*TANB*ABR)
+        TEMPS=TEMPS/RBB/4.0
+        B1RPD=FB2-ZAP*GGP/8.0-TEMPS
+        B2RPD=FB2-ZAP*GGP/8.0+TEMPS
+C
+C        calculate second derivatives w.r.t. H_R
+C
+        CB3=V*ZAP*GG2*(BLRM-ZAP*GG2*VPVM/4.0)/SR2
+        CB3=CB3+4.0*SR2*FB2*V*EP*(EP+COTB*ABR)
+        A3=-CB3**2/(RBB**3)/8.0
+        A4=ZAP*GG2*(BLRM-ZAP*GG2*VPVM/4.0)/2.0
+        A4=A4+V2*ZAP*GG2**2/4.0+4.0*FB2*EP2
+        A4=A4/RBB/4.0
+        B1RR=ZAP*GGP/8.0-A3-A4
+        B2RR=ZAP*GGP/8.0+A3+A4
+C
+C       calculate second derivatives w.r.t. H_R'
+C
+        CB7=-VP*ZAP*GG2*(BLRM-ZAP*GG2*VPVM/4.0)/SR2
+        CB7=CB7+4.0*SR2*FB2*VP*(AB2+EP*TANB*ABR)
+        A7=-CB7**2/(RBB**3)/8.0
+        A8=-ZAP*GG2*(BLRM-ZAP*GG2*VPVM/4.0)/2.0
+        A8=A8+VP2*ZAP*GG2**2/4.0+4.0*FB2*AB2
+        A8=A8/RBB/4.0
+        B1RPRP=FB2-ZAP*GGP/8.0-A7-A8
+        B2RPRP=FB2-ZAP*GGP/8.0+A7+A8
+C
+C       calculate second derivatives w.r.t. H_R and H_R'
+C
+        A10=-VVP*ZAP*(GG2**2)/4.0+4.0*FB2*EP*ABR
+        A10=A10/RBB/4.0
+        B1RRP=CB3*CB7/(RBB**3)/8.0-A10
+        B2RRP=-CB3*CB7/(RBB**3)/8.0+A10
+C
+C       calculate  D^2 V / D^2 H_R
+C
+        TEMPSQ=MSB1SQ*(B1RR-B1RD)
+        DB1=2.0*(2.0*V2*B1RD**2+TEMPSQ)*LOG(MSB1SQ/QQQ2)
+        DB1=DB1+6.0*V2*B1RD**2+TEMPSQ
+        TEMPSQ=MSB2SQ*(B2RR-B2RD)
+        DB2=2.0*(2.0*V2*B2RD**2+TEMPSQ)*LOG(MSB2SQ/QQQ2)
+        DB2=DB2+6.0*V2*B2RD**2+TEMPSQ
+        VRRB=-TM1BB*COTB+3.0*(DB1+DB2)/32.0/PI2
+C
+C       calculate  D^2 V / D^2 H'_R
+C
+        TEMPSQ=MSB1SQ*(B1RPRP-B1RPD)
+        DB1=2.0*(2.0*VP2*B1RPD**2+TEMPSQ)*LOG(MSB1SQ/QQQ2)
+        DB1=DB1+6.0*VP2*B1RPD**2+TEMPSQ
+        TEMPSQ=MSB2SQ*(B2RPRP-B2RPD)
+        DB2=2.0*(2.0*VP2*B2RPD**2+TEMPSQ)*LOG(MSB2SQ/QQQ2)
+        DB2=DB2+6.0*VP2*B2RPD**2+TEMPSQ
+        VRPRPB=DB1+DB2
+        VRPRPB=DB1+DB2-8.0*FB2*MB2*LOG(MB2/QQQ2)-12.0*FB2*MB2
+        VRPRPB=-TM1BB*TANB+3.0*VRPRPB/32.0/PI2
+C
+C       calculate  D^2 V / D H_R  D H'_R
+C
+        DB1=2.0*VVP*B1RD*B1RPD+MSB1SQ*B1RRP
+        DB1=2.0*DB1*LOG(MSB1SQ/QQQ2)
+        DB1=DB1+6.0*VVP*B1RD*B1RPD+MSB1SQ*B1RRP
+        DB2=2.0*VVP*B2RD*B2RPD+MSB2SQ*B2RRP
+        DB2=2.0*DB2*LOG(MSB2SQ/QQQ2)
+        DB2=DB2+6.0*VVP*B2RD*B2RPD+MSB2SQ*B2RRP
+        VRRPB=TM1BB+3.0*(DB1+DB2)/32.0/PI2
+      ELSE IF(RBB.EQ.0.0) THEN
+C
+        ALPHAB=BLRP/2.0+MB2-ZAP*GGP*VPVM/8.0
+        LAB=2.0*LOG(ALPHAB/QQQ2)+3.0
+C
+C       calculate  D^2 V / D^2 H_R
+C
+        VRRB=ZAP*V2*(GGP**2 + GG2**2)/16.0
+        VRRB=3.0*(VRRB*LAB)/32.0/PI2
+C
+C       calculate  D^2 V / D^2 H_R'
+C
+        VRPRPB=VP2*(GGP**2+GG2**2)/16.0-MB2*GGP         
+        VRPRPB=ZAP*VRPRPB*LAB+8.0*FB2*MB2*LOG(ALPHAB/MB2)
+        VRPRPB=3.0*VRPRPB/32.0/PI2
+C
+C       calculate  D^2 V / D^H_R D^H_R'
+C
+        VRRPB=FB2*GGP-(GGP**2+GG2**2)/8.0
+        VRRPB=ZAP*VVP*VRRPB*LAB/2.0
+        VRRPB=3.0*VRRPB/32.0/PI2
+C
+      ENDIF
+C
+      DVRR=VRRT+VRRB+VP2*MHP2/VVPP + V2*GGP/2.0
+      DVRPRP=VRPRPT+VRPRPB+V2*MHP2/VVPP + VP2*GGP/2.0
+      DVRRP=VRRPT+VRRPB-VVP*MHP2/VVPP - VVP*GGP/2.0
+C          TEMPH is always non-negative:
+      TEMPH=(DVRR-DVRPRP)**2+4*DVRRP**2
+      TEMPH=0.5*SQRT(TEMPH)
+      MHL2=0.5*(DVRR+DVRPRP)-TEMPH
+      MHH2=0.5*(DVRR+DVRPRP)+TEMPH
+      IF(MHL2.LT.0.0) THEN
+        MHLNEG=1
+C        WRITE(LOUT,*) 'SSMHN: ERROR:  MHL**2 < 0.0 FOR PARAMETERS:'
+C        WRITE(LOUT,*) 'MHP =', AMHA, 'TANB =', 1.0/RR
+C        WRITE(LOUT,*) 'MSTL=', AMTLSS, 'MSBL=', AMBLSS 
+C        WRITE(LOUT,*) 'MSTR=', AMTRSS, 'MSBR=', AMBRSS
+C        WRITE(LOUT,*) 'AT=', AAT, 'AB=', AAB
+C        WRITE(LOUT,*) 'MU=-2M1=', -EP
+C        WRITE(LOUT,*) 'MT=', AMTP, 'MB=', AMBT
+C        WRITE(LOUT,*) 'D-TERMS? 1=YES 2=NO :', INRAD
+C        WRITE(LOUT,*) 'MASS SCALE (QQQ)=', SQRT(QQQ2)
+        AMHH=SQRT(MHH2)
+        AMHL=SQRT(ABS(MHL2))
+        GO TO 1000
+      ENDIF
+      AMHL=SQRT(MHL2)
+      AMHH=SQRT(MHH2)
+
+C
+C     Now calculate mixing angle ALFAH
+C
+      TRACEM=DVRR-DVRPRP
+      TPAL=TRACEM**2 + 4.0*DVRRP**2
+      TANAH=TRACEM+SQRT(TPAL)
+      IF(DVRRP.EQ.0.0) THEN
+        WRITE(LOUT,*) 'SSMHN: OFF-DIAGONAL TERM OF SCALAR HIGGS',
+     $  ' MASS MATRIX IS ZERO '
+        IF(TANAH.NE.0.0) THEN
+          WRITE(LOUT,*) 'SSMHN: WARNING: TAN(ALFAH) FORMULA',
+     $    ' YIELDS INFINITY'
+        ELSE IF(TANAH.EQ.0.0) THEN
+          WRITE(LOUT,*) 'SSMHN: WARNING: TAN(ALFAH) FORMULA',
+     $    ' YIELDS 0/0 '
+        ENDIF
+        IF(DVRR.GT.DVRPRP) THEN
+          WRITE(LOUT,*) 'SSMHN: DVRR > DVRPRP ==> SET ALFAH=PI/2'
+          ALFAH = PI/2.0
+        ELSE IF (DVRR .LT. DVRPRP) THEN
+          WRITE(LOUT,*) 'SSMHN: DVRR < DVRPRP ==> SET ALFAH=0'
+          ALFAH = 0.0
+        ELSE IF (DVRR .EQ. DVRPRP) THEN
+          WRITE(LOUT,*) 'SSMHN: DVRR = DVRPRP ==> ALFAH INDETERMINANT'
+          WRITE(LOUT,*) 'SETTING SCALAR MIXING ANGLE ALPHA=PI/4'
+          ALFAH=PI/4.0
+        ENDIF
+        GO TO 1000
+      ENDIF
+      TANAH = -0.5*TANAH/DVRRP
+      ALFAH = ATAN(TANAH)
+C
+1000  RETURN
+      END