From 6f6ae960d0b14af7375297dc8e0bbd273f8fa6df Mon Sep 17 00:00:00 2001
From: Alexander Voigt <Alexander.Voigt@desy.de>
Date: Fri, 2 Sep 2016 20:06:20 +0200
Subject: [PATCH] Bugfix: correcting Ben's catches for NaNs

Fortran statements of the form

   if (condition) expression

does not have an affect with my gfortran 4.9.2.
---
 src/mssm_twoloophiggs_impl.f | 35 ++++++++++++++++++++++++++---------
 1 file changed, 26 insertions(+), 9 deletions(-)

diff --git a/src/mssm_twoloophiggs_impl.f b/src/mssm_twoloophiggs_impl.f
index fbad2167c..84694e0c2 100644
--- a/src/mssm_twoloophiggs_impl.f
+++ b/src/mssm_twoloophiggs_impl.f
@@ -5569,7 +5569,9 @@ function tauF1ab(t,A0,BL,T1,T2,s2t,c2t,cb,sb,q,mu)
       Nc = 1d0
 
       mu2 = mu**2
-      if(mu2.eq.0d0) mu2 = 1d-10
+      if (mu2.eq.0d0) then
+         mu2 = 1d-10
+      endif
       
       tauF1ab = 
      $     (2.*BL*mu2/t*(mu2-BL+t)/delt(BL,mu2,t)
@@ -5608,7 +5610,9 @@ function tauF1c(t,A0,BL,T1,T2,s2t,c2t,cb,sb,q,mu)
       Nc = 1d0
       
       mu2 = mu**2
-      if(mu2.eq.0d0) mu2 = 1d-10
+      if (mu2.eq.0d0) then
+         mu2 = 1d-10
+      endif
       
       Xt  = s2t*(T1-T2)/2d0/Sqrt(t)
       Yt  = s2t*(T1-T2)/2d0/Sqrt(t) - mu/sb/cb
@@ -5721,7 +5725,9 @@ function tauF2c(t,A0,BL,T1,T2,s2t,c2t,cb,sb,q,mu)
       Nc = 1d0
       
       mu2 = mu**2
-      if(mu2.eq.0d0) mu2 = 1d-10
+      if (mu2.eq.0d0) then
+         mu2 = 1d-10
+      endif
       
       Xt = s2t*(T1-T2)/2d0/Sqrt(t)
       Yt = Xt - mu/cb/sb
@@ -5837,7 +5843,9 @@ function tauF3c(t,A0,BL,T1,T2,s2t,c2t,cb,sb,q,mu)
       Nc = 1d0
       
       mu2 = mu**2
-      if(mu2.eq.0d0) mu2 = 1d-10
+      if (mu2.eq.0d0) then
+         mu2 = 1d-10
+      endif
       
       Xt = s2t*(T1-T2)/2d0/Sqrt(t)
       Yt = Xt - mu/cb/sb
@@ -6316,12 +6324,16 @@ function tauFAc(t,A0,BL,T1,T2,s2t,c2t,cb,sb,q,mu)
       Nc = 1d0
       
       mu2 = mu**2
-      if(mu2.eq.0d0) mu2 = 1d-10
+      if (mu2.eq.0d0) then
+         mu2 = 1d-10
+      endif
       
       Xt = s2t*(T1-T2)/2d0/Sqrt(t)
       Yt = Xt - mu/cb/sb
       At = sb**2*Xt+cb**2*Yt
-      if (At.eq.0d0) At=1d-10
+      if (At.eq.0d0) then
+         At=1d-10
+      endif
       
       ct2 = (1d0+c2t)/2d0
       st2 = (1d0-c2t)/2d0
@@ -6462,11 +6474,13 @@ subroutine tausqtad(t,A0,BL,T1,T2,st,ct,q,mu,tb,vv,
 
 c     ADDED by BEN: guards against NANs when sin theta is zero!
       if (dabs(st).lt.1.0d-10) then
-         if (st.ge.0.0d0) st=1.0d-10
+         if (st.ge.0.0d0) then
+            st=1.0d-10
          else 
             st=-1.0d-10
          endif
-         ct = sqrt(1.d0-st*st)
+      endif
+      ct = sqrt(1.d0-st*st)
 c     end of addition by BEN 13/6/12
       
       s2t = 2d0*ct*st
@@ -6523,12 +6537,15 @@ subroutine tautadfuncs(t,A0,BL,T1,T2,s2t,c2t,q,mu,sb,cb,F2l,G2l)
       cb2 = cb**2
       mu2 = mu**2
 
-      if(mu2.eq.0d0) mu2 = 1d-10
+      if (mu2.eq.0d0) then
+         mu2 = 1d-10
+      endif
 
       Xt = s2t*(T1-T2)/2d0/Sqrt(t)
       Yt = Xt - mu/cb/sb
       At = sb2*Xt+cb2*Yt
       
+
       F2l = (delt(mu2,t,T1)+2*mu2*t)/T1*phi(mu2,t,T1)
      $     -(delt(mu2,t,T2)+2*mu2*t)/T2*phi(mu2,t,T2)
      $     +A0*cb2*(2*Sqrt(t)+s2t*Yt)/4/s2t/T1/(T1-T2)*