*------------------------------------------------------------------------*
*    This program calculates all the Higgs decay channels for mH<200 GeV *
*    as described in the text in http:/porthos.ist.utl.pt/CTQFT.         *
*    Various tests are included for comparison with other results in     *
*    the literature.                                                     *
*                                                                        *
*    Needs: LoopTools, CUBA integration package, gauss.f for gaussian    *
*           integration.                                                 *
*                                                                        *
*    Author: Jorge Crispim Romao                                         *
*    email:jorge.romao@ist.utl.pt                                        *
*------------------------------------------------------------------------*

      program hgg
      implicit none


**************  Standard Model Physical Constants ******
*
      INCLUDE 'smpc.h'
*
********************************************************

*
*** Variables in the program
*
      integer i,j,iboson
      real*8 lbd2    ! For LoopTools
      real*8 g,G_F,mW,mZ,thetaW,mH,GW,GZ,pi
      real*8 ee,alpha,alphaz,GammaH(10),BR(10),GammaHTotal
      real*8 mlep(3),mup(3),mdown(3)
      real*8 GammaHVVstar,GammaHVV,GammaHFF,GammaHFermions
      real*8 GammaHGG,GammaHGluonGluon,GammaHGZ
      complex*16 XF,XG

      common/data/g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi,iboson
      common/data2/mlep,mup,mdown

*
*** Initialize LoopTools
*
#include "looptools.h"

      call ltini

      lbd2=1.d0
      call setlambda(lbd2)


      open(50,file='BRHiggs.dat',status='unknown')
      open(60,file='XFXG.dat',status='unknown')

      pi=acos(-1.d0)
      mW=SM_MW
      GW=SM_GAW
      mZ=SM_MZ
      GZ=SM_GAZ
      G_F=SM_GFERMI
      mlep(1)=SM_EMASS
      mlep(2)=SM_MUMASS
      mlep(3)=SM_TAUMASS
      mup(3)=SM_TOPMASS
      mup(2)=SM_CHARMMASS
      mup(1)=SM_UPMASS
      mdown(1)=SM_DOWNMASS
      mdown(2)=SM_STRANGEMASS
      mdown(3)=SM_BOTTOMMASS
      alpha=SM_ALPHA
      alphaz=SM_ALPHAZ
      thetaW=asin(sqrt(SM_S2TW))
      ee=sqrt(4*pi*alpha)
      g=ee/sin(thetaW)

      mW=mZ*cos(thetaW)
      call setmudim(mW**2)

      mH=38d0

      do i=1,162
         mH=mH+1d0

         GammaHFermions=0d0
         do j=1,3
            GammaHFermions=GammaHFermions+GammaHFF(mlep(j))
            GammaHFermions=GammaHFermions+3d0*GammaHFF(mdown(j))
         enddo
         do j=1,2
            GammaHFermions=GammaHFermions+3d0*GammaHFF(mup(j))
         enddo

         GammaH(1)=3d0*GammaHFF(mdown(3))
         GammaH(2)=3d0*GammaHFF(mup(2))
         GammaH(3)=GammaHFF(mlep(3))
         GammaH(4)=GammaHFF(mlep(2))
         GammaH(5)=GammaHVVstar(1)
         GammaH(6)=GammaHVVstar(2)
         GammaH(7)=GammaHGG()
         GammaH(8)=GammaHGluonGluon()
         GammaH(9)=GammaHGZ()

         GammaHTotal=GammaHFermions
         do j=5,9
            GammaHTotal=GammaHTotal+GammaH(j)
         enddo
         do j=1,9
            BR(j)=GammaH(j)/GammaHTotal
         enddo

         
         write(50,98)mH,(BR(j),j=1,9)

         write(60,97)mH,XF(mup(3)),XG(),XF(mlep(3))

   
      enddo

      


 98   format(f8.2,4x,9(e14.5,2x))
 97   format(f8.2,4x,6(e14.5,2x))


      end

*-------------------------------------------------------------------*
      real*8 function GammaHGG()
      implicit none
      integer i,iboson
      real*8 mV,cst,beta,mf
      real*8 mlep(3),mup(3),mdown(3)
      real*8 g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi
      complex*16 X,XF,XG

      common/data/g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi,iboson
      common/data2/mlep,mup,mdown


      cst=G_F*mH**3/8d0/pi/sqrt(2d0)*(alpha**2/16d0/pi**2)
      
      X=0d0
      do i=1,3
         X=X+XF(mlep(i))
         X=X+4d0/3d0*XF(mup(i))
         X=X+1d0/3d0*XF(mdown(i))
      enddo

      X=X+XG()

      GammaHGG=cst*abs(X)**2

      return
      end

*-------------------------------------------------------------------*
      real*8 function GammaHGZ()
      implicit none
      integer i,iboson
      real*8 mV,cst,beta,mf,sw2,Nc,gv,qf
      real*8 mlep(3),mup(3),mdown(3)
      real*8 g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi
      complex*16 Y,YF,YG

      common/data/g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi,iboson
      common/data2/mlep,mup,mdown

      
      if(mH.lt.mZ)then
         GammaHGZ=0d0
         return
      endif
      
      sw2=sin(thetaW)**2
      beta=1d0-mZ**2/mH**2
      cst=alpha**2/(64*pi**3)*G_F/sqrt(2d0)*mH**3*beta**3
      
      Y=0d0
      do i=1,3
         Nc=1d0
         qf=-1d0
         gv=-1d0/4d0-qf*sw2
         Y=Y+Nc*qf*gv*YF(mlep(i))
         qf=2d0/3d0
         Nc=3d0
         gv=1d0/4d0-qf*sw2
         Y=Y+Nc*qf*gv*YF(mup(i))
         qf=-1d0/3d0
         Nc=3d0
         gv=-1d0/4d0-qf*sw2
         Y=Y+Nc*qf*gv*YF(mdown(i))
      enddo

      Y=Y+YG()

      GammaHGZ=cst*abs(Y)**2

      return
      end

*-------------------------------------------------------------------*
      real*8 function GammaHGluonGluon()
      implicit none
      integer i,iboson
      real*8 mV,cst,beta,alpha_s
      real*8 mlep(3),mup(3),mdown(3)
      real*8 g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi
      complex*16 X,XF,XG

      common/data/g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi,iboson
      common/data2/mlep,mup,mdown

      alpha_s=0.1d0
      cst=alpha_s**2*g**2*mH**3/512d0/pi**3/mW**2
      

      X=0d0
      do i=1,3
         X=X+XF(mup(i))
         X=X+XF(mdown(i))
      enddo

      GammaHGluonGluon=cst*abs(X)**2

      return
      end

*-------------------------------------------------------------------*
      real*8 function GammaHFF(mf)
      implicit none
      integer iboson
      real*8 mV,cst,beta,mf
      real*8 g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi

      common/data/g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi,iboson

      cst=g**2/32d0/pi/mW**2*mf**2
      
      beta=sqrt(1d0-4d0*mf**2/mH**2)

      GammaHFF=cst*mH*beta**3

      return
      end

*-------------------------------------------------------------------*
      real*8 function GammaHVV(ib)
      implicit none
      integer iboson,ib
      real*8 mV,cst
      real*8 g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi

      common/data/g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi,iboson

      if(ib.eq.1)then
         mV=mW
         cst=g**2
      else
         mV=mZ
         cst=1d0/2d0*g**2/cos(thetaW)**2
      endif
      
      if(mH.lt.2d0*mV)then
         GammaHVV=0d0
      else
      GammaHVV=1d0/64d0/pi*cst*mH**3/mV**2*sqrt(1d0-4d0*mV**2/mH**2)*
     &     (1d0 + 12d0*(mV/mH)**4-4d0*(mV/mH)**2)
      endif

      return
      end


*-------------------------------------------------------------------*
      real*8 function GammaHVVstar(ib)
      implicit none
      integer iflag,iboson,ib
      real*8 res,eps,eta,test1,test2,cst,mV
      real*8 xmin(1),xmax(1)
      real*8 g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi

      common/boundary/xmin,xmax
      common/data/g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi,iboson
      
      iboson=ib

      GammaHVVstar=0d0

      if(iboson.eq.1)then
         mV=mW
         cst=3d0*g**4*mH/512d0/pi**3
      else
         mV=mZ
         cst=g**4*mH/2048d0/pi**3/cos(thetaW)**4*(7d0-40d0/3d0*
     &        sin(thetaW)**2 +160d0/9d0*sin(thetaW)**4)
      endif

      if(mH.le.mV)then
         return
      endif

      eps=mV/mH
      eta=8d0*(gW/mH)*eps
      test1=2d0*eps
      test2=1+eps**2

      if(test1.gt.(1d0+eta))then
         iflag=1
      endif

      if(test1.gt.(1d0-eta).and.test1.lt.(1d0+eta))then
         iflag=2
      endif

      if(test1.lt.(1d0-eta))then
         if(test2.lt.(1d0+eta))then
            iflag=3
         else
            iflag=4
         endif
      endif
      
      if(iflag.eq.1)then
         xmin(1)=2d0*eps
         xmax(1)=1d0+eps**2
         call Integration(res)
         GammaHVVstar=GammaHVVstar + res
      endif

      if(iflag.eq.2)then
         xmin(1)=2d0*eps
         xmax(1)=1d0+eta
         call Integration(res)
         GammaHVVstar=GammaHVVstar + res
         xmin(1)=1d0+eta
         xmax(1)=1d0+eps**2
         call Integration(res)
         GammaHVVstar=GammaHVVstar + res
      endif
      if(iflag.eq.3)then
         xmin(1)=2d0*eps
         xmax(1)=1d0-eta
         call Integration(res)
         GammaHVVstar=GammaHVVstar + res
         xmin(1)=1d0-eta
         xmax(1)=1d0+eps**2
         call Integration(res)
         GammaHVVstar=GammaHVVstar + res
      endif
      if(iflag.eq.4)then
         xmin(1)=2d0*eps
         xmax(1)=1d0-eta
         call Integration(res)
         GammaHVVstar=GammaHVVstar + res
         xmin(1)=1d0-eta
         xmax(1)=1d0+eta
         call Integration(res)
         GammaHVVstar=GammaHVVstar + res
         xmin(1)=1d0+eta
         xmax(1)=1d0+eps**2
         call Integration(res)
         GammaHVVstar=GammaHVVstar + res
      endif


         GammaHVVstar=cst*GammaHVVstar

      return
      end
     

      subroutine Integration(res)
      implicit none
      real*8 res
*
**   Modify the parameter method to choose the method of integration. 
*
      integer i,j,npoints,method
      parameter (npoints = 128)
      parameter (method = 1)
*
***   Variables for all the methods
*
      integer ndim,ncomp
      parameter (ndim = 1)
      parameter (ncomp = 1)
      double precision integral(ncomp)

*
***   Other Variables for vegas,suave,divonne and cuhre
*

      integer verbose,mineval,maxeval
      parameter (verbose = 0)
      parameter (mineval = 0)
      parameter (maxeval = 50000)

      double precision epsrel, epsabs
      parameter (epsrel = 1D-3)
      parameter (epsabs = 1D-12)

      double precision error(ncomp), prob(ncomp)
      integer nregions, neval, fail
      

*
***   Variables only in vegas
*

      integer nstart, nincrease
      parameter (nstart = 5000)
      parameter (nincrease = 500)

      external integrand

      if(method.eq.1)then
         call gauss(ndim, ncomp, npoints, integrand, integral, fail)
         if(fail.eq.1)then
            write(*,*)'Gauss method only works for ndim <= 6'
            stop
         endif

         if(fail.eq.2)then
            write(*,*)'Gauss method only works for
     &           npoints=4,6,8,10,16,20,24,32,64,128'
            stop
         endif
         
      endif

      if(method.eq.2)then
         call vegas(ndim, ncomp, integrand,
     &        epsrel, epsabs, verbose, mineval, maxeval,
     &        nstart, nincrease,
     &        neval, fail, integral, error, prob)
      endif

      res=integral(1)
      
      return
      end

*-------------------------------------------------------------------*
*
*** The integration has to be mapped from the interval [0,1] into
*** [xmin,xmax] 
*
*-------------------------------------------------------------------*
      subroutine integrand(ndim, xx, ncomp, ff)
      implicit none
      integer i, ndim, ncomp
      double precision xx(*), ff(*)
      double precision x(ndim),xmin(1),xmax(1),Jacob
      double precision dfdx

      common/boundary/xmin,xmax

      Jacob=1.d0

      do i=1,ndim
         Jacob=Jacob*(xmax(i)-xmin(i))
         x(i)=xmin(i)+(xmax(i)-xmin(i))*xx(i)
      enddo
      
      ff(1) = dfdx(x(1))*Jacob

      return
      end

*-------------------------------------------------------------------*
*
*** This is the function to integrate. 
*
*-------------------------------------------------------------------*
      double precision function dfdx(x)
      implicit none
      integer iboson
      real*8 x
      real*8 g,alpha,G_F,mW,MZ,thetaW,mH,GW,GZ,pi
      real*8 eps,delta

      common/data/g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi,iboson

      if(iboson.eq.1)then
         eps=mW/mH
         delta=GW/mH
      else
         eps=mZ/mH
         delta=GZ/mH
      endif            

      dfdx=sqrt(x**2-4d0*eps**2)/((x-1d0)**2+(eps*delta)**2)*
     &     (x**2-12d0*eps**2*x+8d0*eps**2+12d0*eps**4)
      return
      end


*-------------------------------------------------------------------*
      complex*16 function XF(mf)
      implicit none
      integer iboson
      real*8 mf,tau,mf2,mH2
      real*8 g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi
      complex*16 F12


      common/data/g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi,iboson

#include "looptools.h"

      mf2=mf*mf
      mH2=mH*mH
      tau=4d0*mf**2/mH**2

      XF=-4d0*(tau-1d0)*C0(0,0,mH2,mf2,mf2,mf2)*mf2-2d0*tau

      if(.FALSE.)then
         write(*,112)MH,XF,f12(tau)
 112     format(f8.4,1x,2(e13.5),1x,2e13.5)
      endif


      return
      end

*-------------------------------------------------------------------*
      complex*16 function YF(mf)
      implicit none


      integer iboson
      real*8 mf,tau,mf2,mH2,MZ2,lbd,cw,sw
      real*8 g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi
      complex*16 aux,I1,I2,I_F,I_F_2

      common/data/g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi,iboson

#include "looptools.h"

      MZ2=mZ*mZ
      mf2=mf*mf
      mH2=mH*mH
      tau=4d0*mf2/MH2
      lbd=4d0*mf2/MZ2
      cw=cos(thetaW)
      sw=sin(thetaW)

      I_F=(-4*mf2*(-2*MZ2*B0(MH2,mf2,mf2) + 2*MZ2*B0(MZ2,mf2,mf2) + 
     &      (MH2 - MZ2)*(-2 + 
     &      (-4*mf2 + MH2 - MZ2)*C0(MZ2,0,MH2,mf2,mf2,mf2))))/
     &      (MH2 - MZ2)**2


************************************************************************
*    Aqui comparamos com as expressões do HHG em termos dos integrais
*    I1 e I2 (aux). Para fazer o teste passar
*
*    .FALSE.   --> .TRUE.
*
*    O resultado aparece no ecran. 
*
*    Nota: I_F e I_F_2 são iguais. I_F é o output para Fortran direto
*          do Mathematica. I_F_2 e a expressao que aparece mais "bonita"
*          nos meus slides.
*    
************************************************************************



      if(.FALSE.)then
      aux=-4d0*(I1(tau,lbd)-I2(tau,lbd))
      I_F=(8*mf2*mZ2)/(mH2 - mZ2)**2*(B0(mH2,mf2,mf2)-B0(mZ2,mf2,mf2))
     &     - 4*mf2/(mH2 - mZ2)*(
     &     -2 + (-4*mf2 + mH2 - mZ2)*C0(mZ2,0,mH2,mf2,mf2,mf2))
      
      write(*,111)MH,aux,I_F
 111  format(f8.4,1x,2(e13.5),1x,2e13.5)
      endif

      YF=I_F/sw/cw

      return
      end



*-------------------------------------------------------------------*
      complex*16 function XG()
      implicit none
      real*8 tau,mW2,mH2
      real*8 g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi,iboson
      complex*16 aux1,J1,J2,aux2,I1,I2

      common/data/g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi,iboson

#include "looptools.h"

      mW2=mW*mW
      mH2=mH*mH
      tau=4d0*mW2/mH2

      XG=6d0*(tau-2d0)*C0(0,0,mH2,mW2,mW2,mW2)*mW2+3d0*tau+2d0

      return
      end

*-------------------------------------------------------------------*
      complex*16 function YG()
      implicit none
      real*8 tau,lbd,mW2,mH2,tw,mZ2
      real*8 g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi,iboson
      complex*16 aux1,J1,J2,aux2,I1,I2,I_W,I_W_2

      common/data/g,alpha,G_F,mW,mZ,thetaW,mH,GW,GZ,pi,iboson

#include "looptools.h"

      mZ2=mZ*mZ
      mW2=mW*mW
      mH2=mH*mH
      tau=4d0*mW2/mH2
      lbd=4d0*mW2/mZ2
      tw=tan(thetaW)

      I_W= 
     & -((MZ2*(-2*MW2*(-5 + tw**2) - MH2*(-1 + tw**2))*B0(MH2,MW2,MW2) + 
     &     MZ2*(2*MW2*(-5 + tw**2) + MH2*(-1 + tw**2))*B0(MZ2,MW2,MW2) + 
     &    (MH2 - MZ2)*(MH2 + 10*MW2 - MH2*tw**2 - 2*MW2*tw**2 + 
     &      2*MW2*(MH2*(-5 + tw**2) - 
     &        2*(MW2*(-5 + tw**2) + MZ2*(-3 + tw**2)))*
     &      C0(MZ2,0,MH2,MW2,MW2,MW2)))/(MH2 - MZ2)**2)


************************************************************************
*    Aqui comparamos com as expressões do HHG em termos dos integrais
*    I1 e I2 (aux2) e com o paper de 86 com o Augusto Barroso em termos das 
*    funcoes J1 e J2 (aux1). Para fazer o teste passar
*
*    .FALSE.   --> .TRUE.
*
*    O resultado aparece no ecran. 
*
*    Nota: I_W e I_W_2 são iguais. I_W é o output para Fortran direto
*          do Mathematica. I_W_2 e a expressao que aparece mais "bonita"
*          nos meus slides.
*    
************************************************************************

      if(.FALSE.)then
      aux1=-4d0*((3d0-tw**2)*J1(mZ2/mW2,mH2/mW2) + 
     &     (-5d0+tw**2-0.5d0*mH2/mW2*(1d0-tw**2))*J2(mZ2/mW2,mH2/mW2))

      aux2=-4d0*(3d0-tw**2)*I2(tau,lbd)-((1d0+2d0/tau)*tw**2-
     &     (5d0+2d0/tau))*I1(tau,lbd)

      I_W_2=-1d0/(mH2 - mZ2)**2*(mH2*(1-tw**2) - 2*mW2*(-5 +tw**2))*mZ2*
     &     (B0(mH2, mW2, mW2)-B0(mZ2, mW2,mW2))
     &     -1d0/(mH2 - mZ2)*(mH2*(1-tw**2) -2*mW2*(-5 +tw**2)
     &     + 2*mW2*((-5 + tw**2)*(mH2 -2*MW2) - 2*mZ2*(-3 +tw**2))*
     &     C0(mZ2, 0, mH2, mW2, mW2, mW2))

      write(*,111)aux1,aux2,I_W_2
 111  format(2(e13.5),1x,2e13.5,1x,2e13.5)
      endif

      YG=I_W/tw

      return
      end

*-------------------------------------------------------------------*
      complex*16 function J1(x,y)
      implicit none
      real*8 x,y
      complex*16 S

      J1=-(S(x)-S(y))/(x-y)
    
      return
      end

*-------------------------------------------------------------------*
      complex*16 function J2(x,y)
      implicit none
      real*8 x,y
      complex*16 S,L

      J2=0.5d0/(x-y)+1d0/(x-y)**2*(S(x)-S(y)-0.5d0*x*(L(x)-L(y)))
    
      return
      end

*-------------------------------------------------------------------*
      complex*16 function S(x)
      implicit none
      real*8 x,pi,ReS,ImS

      pi=acos(-1d0)

      if(x.lt.4d0)then
         S=-2d0*(asin(sqrt(x)/2d0))**2
      else
         reS=-pi**2/2d0+0.5d0*(log((x-2d0+sqrt(x**2-4d0*x))/2d0))**2
         ImS=-2d0*pi*log(sqrt(x)/2d0+sqrt(x/4d0-1d0))
         S=cmplx(reS,ImS)
      endif
      return
      end



*-------------------------------------------------------------------*
      complex*16 function L(x)
      implicit none
      real*8 x,ReL,ImL,pi
      
      pi=acos(-1d0)

      if(x.lt.4d0)then
         L=-2d0+2d0*sqrt((4d0-x)/x)*asin(sqrt(x)/2d0)
      else
         ReL=-2d0+sqrt((x-4d0)/x)*log((x-2d0+sqrt(x*(x-4d0)))/2d0)
         ImL=-pi*sqrt((x-4d0)/x)
         L=cmplx(ReL,ImL)
      endif
      return
      end



*-------------------------------------------------------------------*
      complex*16 function F12(x)
      implicit none
      real*8 x
      complex*16 ftau

      F12=-2d0*x*(1d0+(1d0-x)*ftau(x))

      return
      end

*-------------------------------------------------------------------*
      complex*16 function ftau(x)
      implicit none
      real*8 x,realftau,imgftau,etaplus,etaminus,pi

      pi=acos(-1d0)

      if(x.gt.1d0)then
         ftau=asin(sqrt(1/x))**2
      else
         etaplus=1d0+sqrt(1d0-x)
         etaminus=1d0-sqrt(1d0-x)
         realftau=pi**2/4d0-1d0/4d0*log(etaplus/etaminus)**2
         imgftau=pi/2d0*log(etaplus/etaminus)
         ftau=cmplx(realftau,imgftau)
      endif
      return
      end
*-------------------------------------------------------------------*

*-------------------------------------------------------------------*
      complex*16 function gtau(x)
      implicit none
      real*8 x,realgtau,imggtau,etaplus,etaminus,pi

      pi=acos(-1d0)

      if(x.gt.1d0)then
         gtau=sqrt(x-1d0)*asin(sqrt(1/x))
      else
         etaplus=1d0+sqrt(1d0-x)
         etaminus=1d0-sqrt(1d0-x)
         realgtau=0.5d0*sqrt(1d0-x)*log(etaplus/etaminus)
         imggtau=-pi/2d0*sqrt(1d0-x)
         gtau=cmplx(realgtau,imggtau)
      endif
      return
      end
*-------------------------------------------------------------------*

*-------------------------------------------------------------------*
      complex*16 function I1(x,y)
      implicit none
      real*8 x,y
      complex*16 ftau,gtau

      I1=x*y/2d0/(x-y)+(x*y)**2/2d0/(x-y)**2*(ftau(x)-ftau(y))+
     &     x**2*y/(x-y)**2*(gtau(x)-gtau(y))
    
      return
      end

*-------------------------------------------------------------------*
      complex*16 function I2(x,y)
      implicit none
      real*8 x,y
      complex*16 ftau

      I2=-x*y/2d0/(x-y)*(ftau(x)-ftau(y))
    
      return
      end

*-------------------------------------------------------------------*
