! -----------------------------------------------------------------------------
! This file was automatically created by SARAH version 4.12.2
! SARAH References: arXiv:0806.0538, 0909.2863, 1002.0840, 1207.0906, 1309.7223
! (c) Florian Staub, 2013
! ------------------------------------------------------------------------------
! File created at 18:49 on 2.10.2017
! ----------------------------------------------------------------------

Module Wrapper_OneLoopDecay_Glu_seesaw1nmssm
Use Model_Data_seesaw1nmssm
Use Kinematics
Use OneLoopDecay_Glu_seesaw1nmssm
Use Control
Use Settings

Contains

Subroutine OneLoopDecay_Glu(MSdOS,MSd2OS,MSvImOS,MSvIm2OS,MSvReOS,MSvRe2OS,           &
& MSuOS,MSu2OS,MSeOS,MSe2OS,MhhOS,Mhh2OS,MAhOS,MAh2OS,MHpmOS,MHpm2OS,MChiOS,             &
& MChi2OS,MFvOS,MFv2OS,MChaOS,MCha2OS,MFeOS,MFe2OS,MFdOS,MFd2OS,MFuOS,MFu2OS,            &
& MGluOS,MGlu2OS,MVZOS,MVZ2OS,MVWmOS,MVWm2OS,ZDOS,ZVIOS,ZVROS,ZUOS,ZEOS,ZHOS,            &
& ZAOS,ZPOS,ZNOS,UVOS,UMOS,UPOS,ZELOS,ZEROS,ZDLOS,ZDROS,ZULOS,ZUROS,MAh,MAh2,            &
& MCha,MCha2,MChi,MChi2,MFd,MFd2,MFe,MFe2,MFu,MFu2,MFv,MFv2,MGlu,MGlu2,Mhh,              &
& Mhh2,MHpm,MHpm2,MSd,MSd2,MSe,MSe2,MSu,MSu2,MSvIm,MSvIm2,MSvRe,MSvRe2,MVWm,             &
& MVWm2,MVZ,MVZ2,pG,TW,UM,UP,UV,v,ZA,ZD,ZDL,ZDR,ZE,ZEL,ZER,ZH,ZN,ZP,ZU,ZUL,              &
& ZUR,ZVI,ZVR,ZW,ZZ,betaH,g1,g2,g3,Yd,Ye,lam,kap,lamN,Yu,Yv,Td,Te,Tlam,Tk,               &
& TLN,Tu,Tv,mq2,ml2,mHd2,mHu2,md2,mu2,me2,mv2,ms2,M1,M2,M3,vd,vu,vS,dg1,dg2,             &
& dg3,dYd,dTd,dYe,dTe,dlam,dTlam,dkap,dTk,dlamN,dTLN,dYu,dTu,dYv,dTv,dmq2,               &
& dml2,dmHd2,dmHu2,dmd2,dmu2,dme2,dmv2,dms2,dM1,dM2,dM3,dvd,dvu,dvS,dpG,dZD,             &
& dZVI,dZVR,dZU,dZE,dZH,dZA,dZP,dZN,dUV,dUM,dUP,dZEL,dZER,dZDL,dZDR,dZUL,dZUR,           &
& dSinTW,dCosTW,dTanTW,ZfVG,ZffG,ZfVP,ZfVZ,ZfVWm,ZfSd,ZfSvIm,ZfSvRe,ZfSu,ZfSe,           &
& Zfhh,ZfAh,ZfHpm,ZfL0,ZfFV,ZfLm,ZfLp,ZfFEL,ZfFER,ZfFDL,ZfFDR,ZfFUL,ZfFUR,               &
& ZfVPVZ,ZfVZVP,cplAhSdcSd,cplAhSucSu,cplcChaFdcSuL,cplcChaFdcSuR,cplcFdChiSdL,          &
& cplcFdChiSdR,cplcFdFdAhL,cplcFdFdAhR,cplcFdFdhhL,cplcFdFdhhR,cplcFdFdVGL,              &
& cplcFdFdVGR,cplcFdFdVPL,cplcFdFdVPR,cplcFdFdVZL,cplcFdFdVZR,cplcFdFuHpmL,              &
& cplcFdFuHpmR,cplcFdFuVWmL,cplcFdFuVWmR,cplcFdGluSdL,cplcFdGluSdR,cplcFuChiSuL,         &
& cplcFuChiSuR,cplcFuFdcHpmL,cplcFuFdcHpmR,cplcFuFdcVWmL,cplcFuFdcVWmR,cplcFuFuAhL,      &
& cplcFuFuAhR,cplcFuFuhhL,cplcFuFuhhR,cplcFuFuVGL,cplcFuFuVGR,cplcFuFuVPL,               &
& cplcFuFuVPR,cplcFuFuVZL,cplcFuFuVZR,cplcFuGluSuL,cplcFuGluSuR,cplChaFucSdL,            &
& cplChaFucSdR,cplChiFdcSdL,cplChiFdcSdR,cplChiFucSuL,cplChiFucSuR,cplGluFdcSdL,         &
& cplGluFdcSdR,cplGluFucSuL,cplGluFucSuR,cplGluGluVGL,cplGluGluVGR,cplhhSdcSd,           &
& cplhhSucSu,cplHpmSucSd,cplSdcHpmcSu,cplSdcSdVG,cplSdcSdVP,cplSdcSdVZ,cplSdcSucVWm,     &
& cplSucSdVWm,cplSucSuVG,cplSucSuVP,cplSucSuVZ,cplVGVGVG,ctcplGluFdcSdL,ctcplGluFdcSdR,  &
& ctcplGluFucSuL,ctcplGluFucSuR,ctcplGluGluVGL,ctcplGluGluVGR,ZcplAhSdcSd,               &
& ZcplAhSucSu,ZcplcFdChiSdL,ZcplcFdChiSdR,ZcplcFdFdAhL,ZcplcFdFdAhR,ZcplcFdFdhhL,        &
& ZcplcFdFdhhR,ZcplcFdFdVGL,ZcplcFdFdVGR,ZcplcFdFdVPL,ZcplcFdFdVPR,ZcplcFdFdVZL,         &
& ZcplcFdFdVZR,ZcplcFdGluSdL,ZcplcFdGluSdR,ZcplcFuChiSuL,ZcplcFuChiSuR,ZcplcFuFuAhL,     &
& ZcplcFuFuAhR,ZcplcFuFuhhL,ZcplcFuFuhhR,ZcplcFuFuVGL,ZcplcFuFuVGR,ZcplcFuFuVPL,         &
& ZcplcFuFuVPR,ZcplcFuFuVZL,ZcplcFuFuVZR,ZcplcFuGluSuL,ZcplcFuGluSuR,ZcplChiFdcSdL,      &
& ZcplChiFdcSdR,ZcplChiFucSuL,ZcplChiFucSuR,ZcplGluFdcSdL,ZcplGluFdcSdR,ZcplGluFucSuL,   &
& ZcplGluFucSuR,ZcplGluGluVGL,ZcplGluGluVGR,ZcplhhSdcSd,ZcplhhSucSu,ZcplSdcSdVG,         &
& ZcplSdcSdVP,ZcplSdcSdVZ,ZcplSucSuVG,ZcplSucSuVP,ZcplSucSuVZ,ZcplVGVGVG,ZRUZD,          &
& ZRUZVI,ZRUZVR,ZRUZU,ZRUZE,ZRUZH,ZRUZA,ZRUZP,ZRUZN,ZRUUV,ZRUUM,ZRUUP,ZRUZEL,            &
& ZRUZER,ZRUZDL,ZRUZDR,ZRUZUL,ZRUZUR,MLambda,em,gs,deltaM,kont,gP1LGlu)

Implicit None
Real(dp),Intent(in) :: g1,g2,g3,mHd2,mHu2,ms2

Complex(dp),Intent(in) :: Yd(3,3),Ye(3,3),lam,kap,lamN(3,3),Yu(3,3),Yv(3,3),Td(3,3),Te(3,3),Tlam,               &
& Tk,TLN(3,3),Tu(3,3),Tv(3,3),mq2(3,3),ml2(3,3),md2(3,3),mu2(3,3),me2(3,3),              &
& mv2(3,3),M1,M2,M3

Real(dp),Intent(in) :: vd,vu,vS

Real(dp),Intent(in) :: MAh(3),MAh2(3),MCha(2),MCha2(2),MChi(5),MChi2(5),MFd(3),MFd2(3),MFe(3),               &
& MFe2(3),MFu(3),MFu2(3),MFv(6),MFv2(6),MGlu,MGlu2,Mhh(3),Mhh2(3),MHpm(2),               &
& MHpm2(2),MSd(6),MSd2(6),MSe(6),MSe2(6),MSu(6),MSu2(6),MSvIm(6),MSvIm2(6),              &
& MSvRe(6),MSvRe2(6),MVWm,MVWm2,MVZ,MVZ2,TW,v,ZA(3,3),ZH(3,3),ZP(2,2),ZZ(2,2),betaH

Complex(dp),Intent(in) :: pG,UM(2,2),UP(2,2),UV(6,6),ZD(6,6),ZDL(3,3),ZDR(3,3),ZE(6,6),ZEL(3,3),ZER(3,3),       &
& ZN(5,5),ZU(6,6),ZUL(3,3),ZUR(3,3),ZVI(6,6),ZVR(6,6),ZW(2,2)

Real(dp),Intent(in) :: dg1,dg2,dg3,dmHd2,dmHu2,dms2,dvd,dvu,dvS,dZH(3,3),dZA(3,3),dZP(2,2),dSinTW,           &
& dCosTW,dTanTW

Complex(dp),Intent(in) :: dYd(3,3),dTd(3,3),dYe(3,3),dTe(3,3),dlam,dTlam,dkap,dTk,dlamN(3,3),dTLN(3,3),         &
& dYu(3,3),dTu(3,3),dYv(3,3),dTv(3,3),dmq2(3,3),dml2(3,3),dmd2(3,3),dmu2(3,3),           &
& dme2(3,3),dmv2(3,3),dM1,dM2,dM3,dpG,dZD(6,6),dZVI(6,6),dZVR(6,6),dZU(6,6),             &
& dZE(6,6),dZN(5,5),dUV(6,6),dUM(2,2),dUP(2,2),dZEL(3,3),dZER(3,3),dZDL(3,3),            &
& dZDR(3,3),dZUL(3,3),dZUR(3,3)

Complex(dp),Intent(in) :: cplAhSdcSd(3,6,6),cplAhSucSu(3,6,6),cplcChaFdcSuL(2,3,6),cplcChaFdcSuR(2,3,6),        &
& cplcFdChiSdL(3,5,6),cplcFdChiSdR(3,5,6),cplcFdFdAhL(3,3,3),cplcFdFdAhR(3,3,3),         &
& cplcFdFdhhL(3,3,3),cplcFdFdhhR(3,3,3),cplcFdFdVGL(3,3),cplcFdFdVGR(3,3),               &
& cplcFdFdVPL(3,3),cplcFdFdVPR(3,3),cplcFdFdVZL(3,3),cplcFdFdVZR(3,3),cplcFdFuHpmL(3,3,2),&
& cplcFdFuHpmR(3,3,2),cplcFdFuVWmL(3,3),cplcFdFuVWmR(3,3),cplcFdGluSdL(3,6),             &
& cplcFdGluSdR(3,6),cplcFuChiSuL(3,5,6),cplcFuChiSuR(3,5,6),cplcFuFdcHpmL(3,3,2),        &
& cplcFuFdcHpmR(3,3,2),cplcFuFdcVWmL(3,3),cplcFuFdcVWmR(3,3),cplcFuFuAhL(3,3,3),         &
& cplcFuFuAhR(3,3,3),cplcFuFuhhL(3,3,3),cplcFuFuhhR(3,3,3),cplcFuFuVGL(3,3),             &
& cplcFuFuVGR(3,3),cplcFuFuVPL(3,3),cplcFuFuVPR(3,3),cplcFuFuVZL(3,3),cplcFuFuVZR(3,3),  &
& cplcFuGluSuL(3,6),cplcFuGluSuR(3,6),cplChaFucSdL(2,3,6),cplChaFucSdR(2,3,6),           &
& cplChiFdcSdL(5,3,6),cplChiFdcSdR(5,3,6),cplChiFucSuL(5,3,6),cplChiFucSuR(5,3,6),       &
& cplGluFdcSdL(3,6),cplGluFdcSdR(3,6),cplGluFucSuL(3,6),cplGluFucSuR(3,6),               &
& cplGluGluVGL,cplGluGluVGR,cplhhSdcSd(3,6,6),cplhhSucSu(3,6,6),cplHpmSucSd(2,6,6),      &
& cplSdcHpmcSu(6,2,6),cplSdcSdVG(6,6),cplSdcSdVP(6,6),cplSdcSdVZ(6,6),cplSdcSucVWm(6,6), &
& cplSucSdVWm(6,6),cplSucSuVG(6,6),cplSucSuVP(6,6),cplSucSuVZ(6,6),cplVGVGVG,            &
& ctcplGluFdcSdL(3,6),ctcplGluFdcSdR(3,6),ctcplGluFucSuL(3,6),ctcplGluFucSuR(3,6),       &
& ctcplGluGluVGL,ctcplGluGluVGR,ZcplAhSdcSd(3,6,6),ZcplAhSucSu(3,6,6),ZcplcFdChiSdL(3,5,6),&
& ZcplcFdChiSdR(3,5,6),ZcplcFdFdAhL(3,3,3),ZcplcFdFdAhR(3,3,3),ZcplcFdFdhhL(3,3,3),      &
& ZcplcFdFdhhR(3,3,3),ZcplcFdFdVGL(3,3),ZcplcFdFdVGR(3,3),ZcplcFdFdVPL(3,3),             &
& ZcplcFdFdVPR(3,3),ZcplcFdFdVZL(3,3),ZcplcFdFdVZR(3,3),ZcplcFdGluSdL(3,6),              &
& ZcplcFdGluSdR(3,6),ZcplcFuChiSuL(3,5,6),ZcplcFuChiSuR(3,5,6),ZcplcFuFuAhL(3,3,3),      &
& ZcplcFuFuAhR(3,3,3),ZcplcFuFuhhL(3,3,3),ZcplcFuFuhhR(3,3,3),ZcplcFuFuVGL(3,3),         &
& ZcplcFuFuVGR(3,3),ZcplcFuFuVPL(3,3),ZcplcFuFuVPR(3,3),ZcplcFuFuVZL(3,3),               &
& ZcplcFuFuVZR(3,3),ZcplcFuGluSuL(3,6),ZcplcFuGluSuR(3,6),ZcplChiFdcSdL(5,3,6),          &
& ZcplChiFdcSdR(5,3,6),ZcplChiFucSuL(5,3,6),ZcplChiFucSuR(5,3,6),ZcplGluFdcSdL(3,6),     &
& ZcplGluFdcSdR(3,6),ZcplGluFucSuL(3,6),ZcplGluFucSuR(3,6),ZcplGluGluVGL,ZcplGluGluVGR,  &
& ZcplhhSdcSd(3,6,6),ZcplhhSucSu(3,6,6),ZcplSdcSdVG(6,6),ZcplSdcSdVP(6,6),               &
& ZcplSdcSdVZ(6,6),ZcplSucSuVG(6,6),ZcplSucSuVP(6,6),ZcplSucSuVZ(6,6),ZcplVGVGVG

Real(dp), Intent(in) :: em, gs
Complex(dp),Intent(in) :: ZfVG,ZffG,ZfVP,ZfVZ,ZfVWm,ZfSd(6,6),ZfSvIm(6,6),ZfSvRe(6,6),ZfSu(6,6),ZfSe(6,6),      &
& Zfhh(3,3),ZfAh(3,3),ZfHpm(2,2),ZfL0(5,5),ZfFV(6,6),ZfLm(2,2),ZfLp(2,2),ZfFEL(3,3),     &
& ZfFER(3,3),ZfFDL(3,3),ZfFDR(3,3),ZfFUL(3,3),ZfFUR(3,3),ZfVPVZ,ZfVZVP

Real(dp),Intent(in) :: MSdOS(6),MSd2OS(6),MSvImOS(6),MSvIm2OS(6),MSvReOS(6),MSvRe2OS(6),MSuOS(6),            &
& MSu2OS(6),MSeOS(6),MSe2OS(6),MhhOS(3),Mhh2OS(3),MAhOS(3),MAh2OS(3),MHpmOS(2),          &
& MHpm2OS(2),MChiOS(5),MChi2OS(5),MFvOS(6),MFv2OS(6),MChaOS(2),MCha2OS(2),               &
& MFeOS(3),MFe2OS(3),MFdOS(3),MFd2OS(3),MFuOS(3),MFu2OS(3),MGluOS,MGlu2OS,               &
& MVZOS,MVZ2OS,MVWmOS,MVWm2OS,ZHOS(3,3),ZAOS(3,3),ZPOS(2,2)

Complex(dp),Intent(in) :: ZDOS(6,6),ZVIOS(6,6),ZVROS(6,6),ZUOS(6,6),ZEOS(6,6),ZNOS(5,5),UVOS(6,6),              &
& UMOS(2,2),UPOS(2,2),ZELOS(3,3),ZEROS(3,3),ZDLOS(3,3),ZDROS(3,3),ZULOS(3,3),            &
& ZUROS(3,3)

Complex(dp),Intent(in) :: ZRUZD(6,6),ZRUZVI(6,6),ZRUZVR(6,6),ZRUZU(6,6),ZRUZE(6,6),ZRUZH(3,3),ZRUZA(3,3),       &
& ZRUZP(2,2),ZRUZN(5,5),ZRUUV(6,6),ZRUUM(2,2),ZRUUP(2,2),ZRUZEL(3,3),ZRUZER(3,3),        &
& ZRUZDL(3,3),ZRUZDR(3,3),ZRUZUL(3,3),ZRUZUR(3,3)

Real(dp), Intent(in) :: MLambda, deltaM
Real(dp), Intent(out) :: gP1LGlu(1,49)
Integer, Intent(out) :: kont
Real(dp) :: MVG,MVP,MVG2,MVP2, helfactor, phasespacefactor
Integer :: i1,i2,i3,i4, isave, gt1, gt2, gt3

Complex(dp) :: ZRUZDc(6, 6)
Complex(dp) :: ZRUZVIc(6, 6)
Complex(dp) :: ZRUZVRc(6, 6)
Complex(dp) :: ZRUZUc(6, 6)
Complex(dp) :: ZRUZEc(6, 6)
Complex(dp) :: ZRUZHc(3, 3)
Complex(dp) :: ZRUZAc(3, 3)
Complex(dp) :: ZRUZPc(2, 2)
Complex(dp) :: ZRUZNc(5, 5)
Complex(dp) :: ZRUUVc(6, 6)
Complex(dp) :: ZRUUMc(2, 2)
Complex(dp) :: ZRUUPc(2, 2)
Complex(dp) :: ZRUZELc(3, 3)
Complex(dp) :: ZRUZERc(3, 3)
Complex(dp) :: ZRUZDLc(3, 3)
Complex(dp) :: ZRUZDRc(3, 3)
Complex(dp) :: ZRUZULc(3, 3)
Complex(dp) :: ZRUZURc(3, 3)
Real(dp) :: MRPGluToFdcSd(3,6),MRGGluToFdcSd(3,6), MRPZGluToFdcSd(3,6),MRGZGluToFdcSd(3,6)
Real(dp) :: MVPGluToFdcSd(3,6)
Real(dp) :: RMsqTreeGluToFdcSd(3,6),RMsqWaveGluToFdcSd(3,6),RMsqVertexGluToFdcSd(3,6)
Complex(dp) :: AmpTreeGluToFdcSd(2,3,6),AmpWaveGluToFdcSd(2,3,6)=(0._dp,0._dp),AmpVertexGluToFdcSd(2,3,6)&
 & ,AmpVertexIRosGluToFdcSd(2,3,6),AmpVertexIRdrGluToFdcSd(2,3,6), AmpSumGluToFdcSd(2,3,6), AmpSum2GluToFdcSd(2,3,6)
Complex(dp) :: AmpTreeZGluToFdcSd(2,3,6),AmpWaveZGluToFdcSd(2,3,6),AmpVertexZGluToFdcSd(2,3,6)
Real(dp) :: AmpSqGluToFdcSd(3,6),  AmpSqTreeGluToFdcSd(3,6)
Real(dp) :: MRPGluToFucSu(3,6),MRGGluToFucSu(3,6), MRPZGluToFucSu(3,6),MRGZGluToFucSu(3,6)
Real(dp) :: MVPGluToFucSu(3,6)
Real(dp) :: RMsqTreeGluToFucSu(3,6),RMsqWaveGluToFucSu(3,6),RMsqVertexGluToFucSu(3,6)
Complex(dp) :: AmpTreeGluToFucSu(2,3,6),AmpWaveGluToFucSu(2,3,6)=(0._dp,0._dp),AmpVertexGluToFucSu(2,3,6)&
 & ,AmpVertexIRosGluToFucSu(2,3,6),AmpVertexIRdrGluToFucSu(2,3,6), AmpSumGluToFucSu(2,3,6), AmpSum2GluToFucSu(2,3,6)
Complex(dp) :: AmpTreeZGluToFucSu(2,3,6),AmpWaveZGluToFucSu(2,3,6),AmpVertexZGluToFucSu(2,3,6)
Real(dp) :: AmpSqGluToFucSu(3,6),  AmpSqTreeGluToFucSu(3,6)
Real(dp) :: MRPGluToChiVG(5),MRGGluToChiVG(5), MRPZGluToChiVG(5),MRGZGluToChiVG(5)
Real(dp) :: MVPGluToChiVG(5)
Real(dp) :: RMsqTreeGluToChiVG(5),RMsqWaveGluToChiVG(5),RMsqVertexGluToChiVG(5)
Complex(dp) :: AmpTreeGluToChiVG(4,5),AmpWaveGluToChiVG(4,5)=(0._dp,0._dp),AmpVertexGluToChiVG(4,5)&
 & ,AmpVertexIRosGluToChiVG(4,5),AmpVertexIRdrGluToChiVG(4,5), AmpSumGluToChiVG(4,5), AmpSum2GluToChiVG(4,5)
Complex(dp) :: AmpTreeZGluToChiVG(4,5),AmpWaveZGluToChiVG(4,5),AmpVertexZGluToChiVG(4,5)
Real(dp) :: AmpSqGluToChiVG(5),  AmpSqTreeGluToChiVG(5)
Real(dp) :: MRPGluToGluAh(3),MRGGluToGluAh(3), MRPZGluToGluAh(3),MRGZGluToGluAh(3)
Real(dp) :: MVPGluToGluAh(3)
Real(dp) :: RMsqTreeGluToGluAh(3),RMsqWaveGluToGluAh(3),RMsqVertexGluToGluAh(3)
Complex(dp) :: AmpTreeGluToGluAh(2,3),AmpWaveGluToGluAh(2,3)=(0._dp,0._dp),AmpVertexGluToGluAh(2,3)&
 & ,AmpVertexIRosGluToGluAh(2,3),AmpVertexIRdrGluToGluAh(2,3), AmpSumGluToGluAh(2,3), AmpSum2GluToGluAh(2,3)
Complex(dp) :: AmpTreeZGluToGluAh(2,3),AmpWaveZGluToGluAh(2,3),AmpVertexZGluToGluAh(2,3)
Real(dp) :: AmpSqGluToGluAh(3),  AmpSqTreeGluToGluAh(3)
Real(dp) :: MRPGluToGluhh(3),MRGGluToGluhh(3), MRPZGluToGluhh(3),MRGZGluToGluhh(3)
Real(dp) :: MVPGluToGluhh(3)
Real(dp) :: RMsqTreeGluToGluhh(3),RMsqWaveGluToGluhh(3),RMsqVertexGluToGluhh(3)
Complex(dp) :: AmpTreeGluToGluhh(2,3),AmpWaveGluToGluhh(2,3)=(0._dp,0._dp),AmpVertexGluToGluhh(2,3)&
 & ,AmpVertexIRosGluToGluhh(2,3),AmpVertexIRdrGluToGluhh(2,3), AmpSumGluToGluhh(2,3), AmpSum2GluToGluhh(2,3)
Complex(dp) :: AmpTreeZGluToGluhh(2,3),AmpWaveZGluToGluhh(2,3),AmpVertexZGluToGluhh(2,3)
Real(dp) :: AmpSqGluToGluhh(3),  AmpSqTreeGluToGluhh(3)
Real(dp) :: MRPGluToGluVG,MRGGluToGluVG, MRPZGluToGluVG,MRGZGluToGluVG
Real(dp) :: MVPGluToGluVG
Real(dp) :: RMsqTreeGluToGluVG,RMsqWaveGluToGluVG,RMsqVertexGluToGluVG
Complex(dp) :: AmpTreeGluToGluVG(4),AmpWaveGluToGluVG(4)=(0._dp,0._dp),AmpVertexGluToGluVG(4)&
 & ,AmpVertexIRosGluToGluVG(4),AmpVertexIRdrGluToGluVG(4), AmpSumGluToGluVG(4), AmpSum2GluToGluVG(4)
Complex(dp) :: AmpTreeZGluToGluVG(4),AmpWaveZGluToGluVG(4),AmpVertexZGluToGluVG(4)
Real(dp) :: AmpSqGluToGluVG,  AmpSqTreeGluToGluVG
Real(dp) :: MRPGluToGluVP,MRGGluToGluVP, MRPZGluToGluVP,MRGZGluToGluVP
Real(dp) :: MVPGluToGluVP
Real(dp) :: RMsqTreeGluToGluVP,RMsqWaveGluToGluVP,RMsqVertexGluToGluVP
Complex(dp) :: AmpTreeGluToGluVP(4),AmpWaveGluToGluVP(4)=(0._dp,0._dp),AmpVertexGluToGluVP(4)&
 & ,AmpVertexIRosGluToGluVP(4),AmpVertexIRdrGluToGluVP(4), AmpSumGluToGluVP(4), AmpSum2GluToGluVP(4)
Complex(dp) :: AmpTreeZGluToGluVP(4),AmpWaveZGluToGluVP(4),AmpVertexZGluToGluVP(4)
Real(dp) :: AmpSqGluToGluVP,  AmpSqTreeGluToGluVP
Real(dp) :: MRPGluToGluVZ,MRGGluToGluVZ, MRPZGluToGluVZ,MRGZGluToGluVZ
Real(dp) :: MVPGluToGluVZ
Real(dp) :: RMsqTreeGluToGluVZ,RMsqWaveGluToGluVZ,RMsqVertexGluToGluVZ
Complex(dp) :: AmpTreeGluToGluVZ(4),AmpWaveGluToGluVZ(4)=(0._dp,0._dp),AmpVertexGluToGluVZ(4)&
 & ,AmpVertexIRosGluToGluVZ(4),AmpVertexIRdrGluToGluVZ(4), AmpSumGluToGluVZ(4), AmpSum2GluToGluVZ(4)
Complex(dp) :: AmpTreeZGluToGluVZ(4),AmpWaveZGluToGluVZ(4),AmpVertexZGluToGluVZ(4)
Real(dp) :: AmpSqGluToGluVZ,  AmpSqTreeGluToGluVZ
Write(*,*) "Calculating one-loop decays of Glu "
kont = 0
MVG = MLambda
MVP = MLambda
MVG2 = MLambda**2
MVP2 = MLambda**2

ZRUZDc = Conjg(ZRUZD)
ZRUZVIc = Conjg(ZRUZVI)
ZRUZVRc = Conjg(ZRUZVR)
ZRUZUc = Conjg(ZRUZU)
ZRUZEc = Conjg(ZRUZE)
ZRUZHc = Conjg(ZRUZH)
ZRUZAc = Conjg(ZRUZA)
ZRUZPc = Conjg(ZRUZP)
ZRUZNc = Conjg(ZRUZN)
ZRUUVc = Conjg(ZRUUV)
ZRUUMc = Conjg(ZRUUM)
ZRUUPc = Conjg(ZRUUP)
ZRUZELc = Conjg(ZRUZEL)
ZRUZERc = Conjg(ZRUZER)
ZRUZDLc = Conjg(ZRUZDL)
ZRUZDRc = Conjg(ZRUZDR)
ZRUZULc = Conjg(ZRUZUL)
ZRUZURc = Conjg(ZRUZUR)

 ! Counter
isave = 1

If (.not.CalcLoopDecay_LoopInducedOnly) Then
!----------------
! Fd Conjg(Sd)
!----------------

!Tree Level
  If (.not.ExternalZfactors) Then
Call Amplitude_Tree_seesaw1nmssm_GluToFdcSd(cplGluFdcSdL,cplGluFdcSdR,MFd,            &
& MGlu,MSd,MFd2,MGlu2,MSd2,AmpTreeGluToFdcSd)

  Else
Call Amplitude_Tree_seesaw1nmssm_GluToFdcSd(ZcplGluFdcSdL,ZcplGluFdcSdR,              &
& MFd,MGlu,MSd,MFd2,MGlu2,MSd2,AmpTreeGluToFdcSd)

  End if

!Real Corrections
If (OSkinematics) Then
  If (.not.ExternalZfactors) Then
 ! OS and no Z-factors
Call Gamma_Real_seesaw1nmssm_GluToFdcSd(MLambda,em,gs,cplGluFdcSdL,cplGluFdcSdR,      &
& MFdOS,MGluOS,MSdOS,MRPGluToFdcSd,MRGGluToFdcSd)

  Else
 ! OS and Z-factors
Call Gamma_Real_seesaw1nmssm_GluToFdcSd(MLambda,em,gs,ZcplGluFdcSdL,ZcplGluFdcSdR,    &
& MFdOS,MGluOS,MSdOS,MRPGluToFdcSd,MRGGluToFdcSd)

  End if
Else
 ! DR and no Z-factors
  If (.not.ExternalZfactors) Then
Call Gamma_Real_seesaw1nmssm_GluToFdcSd(MLambda,em,gs,cplGluFdcSdL,cplGluFdcSdR,      &
& MFd,MGlu,MSd,MRPGluToFdcSd,MRGGluToFdcSd)

  Else
 ! DR and Z-factors
Call Gamma_Real_seesaw1nmssm_GluToFdcSd(MLambda,em,gs,ZcplGluFdcSdL,ZcplGluFdcSdR,    &
& MFd,MGlu,MSd,MRPGluToFdcSd,MRGGluToFdcSd)

  End if
End if

!Self-energy Corrections
Call Amplitude_WAVE_seesaw1nmssm_GluToFdcSd(cplGluFdcSdL,cplGluFdcSdR,ctcplGluFdcSdL, &
& ctcplGluFdcSdR,MFd,MFd2,MGlu,MGlu2,MSd,MSd2,ZfFDL,ZfFDR,ZffG,ZfSd,AmpWaveGluToFdcSd)

!Vertex Corrections
Call Amplitude_VERTEX_seesaw1nmssm_GluToFdcSd(MAh,MCha,MChi,MFd,MFu,MGlu,             &
& Mhh,MHpm,MSd,MSu,MVG,MVP,MVWm,MVZ,MAh2,MCha2,MChi2,MFd2,MFu2,MGlu2,Mhh2,               &
& MHpm2,MSd2,MSu2,MVG2,MVP2,MVWm2,MVZ2,cplcFdFdAhL,cplcFdFdAhR,cplAhSdcSd,               &
& cplChaFucSdL,cplChaFucSdR,cplChiFdcSdL,cplChiFdcSdR,cplGluFdcSdL,cplGluFdcSdR,         &
& cplcFdFdhhL,cplcFdFdhhR,cplcFdFdVGL,cplcFdFdVGR,cplcFdFdVPL,cplcFdFdVPR,               &
& cplcFdFdVZL,cplcFdFdVZR,cplcChaFdcSuL,cplcChaFdcSuR,cplcFuFdcHpmL,cplcFuFdcHpmR,       &
& cplcFuFdcVWmL,cplcFuFdcVWmR,cplGluFucSuL,cplGluFucSuR,cplGluGluVGL,cplGluGluVGR,       &
& cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,cplhhSdcSd,cplHpmSucSd,            &
& cplSdcSdVG,cplSdcSdVP,cplSdcSdVZ,cplSucSdVWm,AmpVertexGluToFdcSd)

If (ShiftIRdiv) Then
Call Amplitude_IR_VERTEX_seesaw1nmssm_GluToFdcSd(MAh,MCha,MChi,MFd,MFu,               &
& MGlu,Mhh,MHpm,MSd,MSu,MVG,MVP,MVWm,MVZ,MAh2,MCha2,MChi2,MFd2,MFu2,MGlu2,               &
& Mhh2,MHpm2,MSd2,MSu2,MVG2,MVP2,MVWm2,MVZ2,cplcFdFdAhL,cplcFdFdAhR,cplAhSdcSd,          &
& cplChaFucSdL,cplChaFucSdR,cplChiFdcSdL,cplChiFdcSdR,cplGluFdcSdL,cplGluFdcSdR,         &
& cplcFdFdhhL,cplcFdFdhhR,cplcFdFdVGL,cplcFdFdVGR,cplcFdFdVPL,cplcFdFdVPR,               &
& cplcFdFdVZL,cplcFdFdVZR,cplcChaFdcSuL,cplcChaFdcSuR,cplcFuFdcHpmL,cplcFuFdcHpmR,       &
& cplcFuFdcVWmL,cplcFuFdcVWmR,cplGluFucSuL,cplGluFucSuR,cplGluGluVGL,cplGluGluVGR,       &
& cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,cplhhSdcSd,cplHpmSucSd,            &
& cplSdcSdVG,cplSdcSdVP,cplSdcSdVZ,cplSucSdVWm,AmpVertexIRdrGluToFdcSd)

 If (ExternalZfactors) Then
  If (OSkinematics) Then
 ! OS and Z-factors
Call Amplitude_IR_VERTEX_seesaw1nmssm_GluToFdcSd(MAhOS,MChaOS,MChiOS,MFdOS,           &
& MFuOS,MGluOS,MhhOS,MHpmOS,MSdOS,MSuOS,MVG,MVP,MVWmOS,MVZOS,MAh2OS,MCha2OS,             &
& MChi2OS,MFd2OS,MFu2OS,MGlu2OS,Mhh2OS,MHpm2OS,MSd2OS,MSu2OS,MVG2,MVP2,MVWm2OS,          &
& MVZ2OS,cplcFdFdAhL,cplcFdFdAhR,cplAhSdcSd,cplChaFucSdL,cplChaFucSdR,cplChiFdcSdL,      &
& cplChiFdcSdR,ZcplGluFdcSdL,ZcplGluFdcSdR,cplcFdFdhhL,cplcFdFdhhR,cplcFdFdVGL,          &
& cplcFdFdVGR,cplcFdFdVPL,cplcFdFdVPR,cplcFdFdVZL,cplcFdFdVZR,cplcChaFdcSuL,             &
& cplcChaFdcSuR,cplcFuFdcHpmL,cplcFuFdcHpmR,cplcFuFdcVWmL,cplcFuFdcVWmR,cplGluFucSuL,    &
& cplGluFucSuR,cplGluGluVGL,cplGluGluVGR,cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,         &
& cplcFuGluSuR,cplhhSdcSd,cplHpmSucSd,cplSdcSdVG,cplSdcSdVP,cplSdcSdVZ,cplSucSdVWm,      &
& AmpVertexIRosGluToFdcSd)

   Else
 ! DR and Z-factors
Call Amplitude_IR_VERTEX_seesaw1nmssm_GluToFdcSd(MAh,MCha,MChi,MFd,MFu,               &
& MGlu,Mhh,MHpm,MSd,MSu,MVG,MVP,MVWm,MVZ,MAh2,MCha2,MChi2,MFd2,MFu2,MGlu2,               &
& Mhh2,MHpm2,MSd2,MSu2,MVG2,MVP2,MVWm2,MVZ2,cplcFdFdAhL,cplcFdFdAhR,cplAhSdcSd,          &
& cplChaFucSdL,cplChaFucSdR,cplChiFdcSdL,cplChiFdcSdR,ZcplGluFdcSdL,ZcplGluFdcSdR,       &
& cplcFdFdhhL,cplcFdFdhhR,cplcFdFdVGL,cplcFdFdVGR,cplcFdFdVPL,cplcFdFdVPR,               &
& cplcFdFdVZL,cplcFdFdVZR,cplcChaFdcSuL,cplcChaFdcSuR,cplcFuFdcHpmL,cplcFuFdcHpmR,       &
& cplcFuFdcVWmL,cplcFuFdcVWmR,cplGluFucSuL,cplGluFucSuR,cplGluGluVGL,cplGluGluVGR,       &
& cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,cplhhSdcSd,cplHpmSucSd,            &
& cplSdcSdVG,cplSdcSdVP,cplSdcSdVZ,cplSucSdVWm,AmpVertexIRosGluToFdcSd)

 End if
 Else
  If (OSkinematics) Then
 ! OS and no Z-factors
Call Amplitude_IR_VERTEX_seesaw1nmssm_GluToFdcSd(MAhOS,MChaOS,MChiOS,MFdOS,           &
& MFuOS,MGluOS,MhhOS,MHpmOS,MSdOS,MSuOS,MVG,MVP,MVWmOS,MVZOS,MAh2OS,MCha2OS,             &
& MChi2OS,MFd2OS,MFu2OS,MGlu2OS,Mhh2OS,MHpm2OS,MSd2OS,MSu2OS,MVG2,MVP2,MVWm2OS,          &
& MVZ2OS,cplcFdFdAhL,cplcFdFdAhR,cplAhSdcSd,cplChaFucSdL,cplChaFucSdR,cplChiFdcSdL,      &
& cplChiFdcSdR,cplGluFdcSdL,cplGluFdcSdR,cplcFdFdhhL,cplcFdFdhhR,cplcFdFdVGL,            &
& cplcFdFdVGR,cplcFdFdVPL,cplcFdFdVPR,cplcFdFdVZL,cplcFdFdVZR,cplcChaFdcSuL,             &
& cplcChaFdcSuR,cplcFuFdcHpmL,cplcFuFdcHpmR,cplcFuFdcVWmL,cplcFuFdcVWmR,cplGluFucSuL,    &
& cplGluFucSuR,cplGluGluVGL,cplGluGluVGR,cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,         &
& cplcFuGluSuR,cplhhSdcSd,cplHpmSucSd,cplSdcSdVG,cplSdcSdVP,cplSdcSdVZ,cplSucSdVWm,      &
& AmpVertexIRosGluToFdcSd)

   Else
 ! DR and no Z-factors
Call Amplitude_IR_VERTEX_seesaw1nmssm_GluToFdcSd(MAh,MCha,MChi,MFd,MFu,               &
& MGlu,Mhh,MHpm,MSd,MSu,MVG,MVP,MVWm,MVZ,MAh2,MCha2,MChi2,MFd2,MFu2,MGlu2,               &
& Mhh2,MHpm2,MSd2,MSu2,MVG2,MVP2,MVWm2,MVZ2,cplcFdFdAhL,cplcFdFdAhR,cplAhSdcSd,          &
& cplChaFucSdL,cplChaFucSdR,cplChiFdcSdL,cplChiFdcSdR,cplGluFdcSdL,cplGluFdcSdR,         &
& cplcFdFdhhL,cplcFdFdhhR,cplcFdFdVGL,cplcFdFdVGR,cplcFdFdVPL,cplcFdFdVPR,               &
& cplcFdFdVZL,cplcFdFdVZR,cplcChaFdcSuL,cplcChaFdcSuR,cplcFuFdcHpmL,cplcFuFdcHpmR,       &
& cplcFuFdcVWmL,cplcFuFdcVWmR,cplGluFucSuL,cplGluFucSuR,cplGluGluVGL,cplGluGluVGR,       &
& cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,cplhhSdcSd,cplHpmSucSd,            &
& cplSdcSdVG,cplSdcSdVP,cplSdcSdVZ,cplSucSdVWm,AmpVertexIRosGluToFdcSd)

 End if
 End if
AmpVertexGluToFdcSd = AmpVertexGluToFdcSd -  AmpVertexIRdrGluToFdcSd! +  AmpVertexIRosGluToFdcSd ! Shift added later
End if

 ! Add Z-factors to have external fields on-shell
 If (ExternalZfactors) Then
! Final State 1
AmpWaveZGluToFdcSd=0._dp
AmpVertexZGluToFdcSd=0._dp
Do gt1=1,3
  Do gt2=1,3
AmpWaveZGluToFdcSd(1,gt2,:) = AmpWaveZGluToFdcSd(1,gt2,:)+ZRUZDL(gt2,gt1)*AmpWaveGluToFdcSd(1,gt1,:)
AmpVertexZGluToFdcSd(1,gt2,:)= AmpVertexZGluToFdcSd(1,gt2,:)+ZRUZDL(gt2,gt1)*AmpVertexGluToFdcSd(1,gt1,:)
AmpWaveZGluToFdcSd(2,gt2,:) = AmpWaveZGluToFdcSd(2,gt2,:)+ZRUZDRc(gt2,gt1)*AmpWaveGluToFdcSd(2,gt1,:)
AmpVertexZGluToFdcSd(2,gt2,:)= AmpVertexZGluToFdcSd(2,gt2,:)+ZRUZDRc(gt2,gt1)*AmpVertexGluToFdcSd(2,gt1,:)
 End Do
End Do
AmpWaveGluToFdcSd=AmpWaveZGluToFdcSd
AmpVertexGluToFdcSd= AmpVertexZGluToFdcSd
! Final State 2
AmpWaveZGluToFdcSd=0._dp
AmpVertexZGluToFdcSd=0._dp
Do gt1=1,6
  Do gt2=1,6
AmpWaveZGluToFdcSd(:,:,gt2) = AmpWaveZGluToFdcSd(:,:,gt2)+ZRUZDc(gt2,gt1)*AmpWaveGluToFdcSd(:,:,gt1)
AmpVertexZGluToFdcSd(:,:,gt2)= AmpVertexZGluToFdcSd(:,:,gt2)+ZRUZDc(gt2,gt1)*AmpVertexGluToFdcSd(:,:,gt1)
 End Do
End Do
AmpWaveGluToFdcSd=AmpWaveZGluToFdcSd
AmpVertexGluToFdcSd= AmpVertexZGluToFdcSd
End if
If (ShiftIRdiv) Then
AmpVertexGluToFdcSd = AmpVertexGluToFdcSd  +  AmpVertexIRosGluToFdcSd
End if

 !Square the amplitude
If (DebugLoopDecays) Then
Write(*,*) "------------------ Glu->Fd conj[Sd] -----------------------"
End if
If (.not.SquareFullAmplitudeDecays) Then
 AmpSumGluToFdcSd = AmpTreeGluToFdcSd
 AmpSum2GluToFdcSd = AmpTreeGluToFdcSd + 2._dp*AmpWaveGluToFdcSd + 2._dp*AmpVertexGluToFdcSd
Else
 AmpSumGluToFdcSd = AmpTreeGluToFdcSd + AmpWaveGluToFdcSd + AmpVertexGluToFdcSd
 AmpSum2GluToFdcSd = AmpTreeGluToFdcSd + AmpWaveGluToFdcSd + AmpVertexGluToFdcSd
End If
If (OnlyTreeLevelContributions) Then
 AmpSumGluToFdcSd = AmpTreeGluToFdcSd
 AmpSum2GluToFdcSd = AmpTreeGluToFdcSd
End if
gt1=1
i4 = isave
  Do gt2=1,3
    Do gt3=1,6
If (((OSkinematics).and.(MGluOS.gt.(MFdOS(gt2)+MSdOS(gt3)))).or.((.not.OSkinematics).and.(MGlu.gt.(MFd(gt2)+MSd(gt3))))) Then
 If (DebugLoopDecays) Then
  Write(*,*) gt2, gt3
  AmpSum2GluToFdcSd = AmpTreeGluToFdcSd
If (OSkinematics) Then
  Call SquareAmp_FtoFS(MGluOS,MFdOS(gt2),MSdOS(gt3),AmpSumGluToFdcSd(:,gt2, gt3),AmpSum2GluToFdcSd(:,gt2, gt3),AmpSqGluToFdcSd(gt2, gt3))
Else
  Call SquareAmp_FtoFS(MGlu,MFd(gt2),MSd(gt3),AmpSumGluToFdcSd(:,gt2, gt3),AmpSum2GluToFdcSd(:,gt2, gt3),AmpSqGluToFdcSd(gt2, gt3))
End if
  Write(*,*) "TREE x TREE: ",AmpSqGluToFdcSd(gt2, gt3)
  AmpSum2GluToFdcSd = 2._dp*AmpWaveGluToFdcSd
If (OSkinematics) Then
  Call SquareAmp_FtoFS(MGluOS,MFdOS(gt2),MSdOS(gt3),AmpSumGluToFdcSd(:,gt2, gt3),AmpSum2GluToFdcSd(:,gt2, gt3),AmpSqGluToFdcSd(gt2, gt3))
Else
  Call SquareAmp_FtoFS(MGlu,MFd(gt2),MSd(gt3),AmpSumGluToFdcSd(:,gt2, gt3),AmpSum2GluToFdcSd(:,gt2, gt3),AmpSqGluToFdcSd(gt2, gt3))
End if
  Write(*,*) "TREE x WAVE: ",AmpSqGluToFdcSd(gt2, gt3)
  AmpSum2GluToFdcSd = 2._dp*AmpVertexGluToFdcSd
If (OSkinematics) Then
  Call SquareAmp_FtoFS(MGluOS,MFdOS(gt2),MSdOS(gt3),AmpSumGluToFdcSd(:,gt2, gt3),AmpSum2GluToFdcSd(:,gt2, gt3),AmpSqGluToFdcSd(gt2, gt3))
Else
  Call SquareAmp_FtoFS(MGlu,MFd(gt2),MSd(gt3),AmpSumGluToFdcSd(:,gt2, gt3),AmpSum2GluToFdcSd(:,gt2, gt3),AmpSqGluToFdcSd(gt2, gt3))
End if
  Write(*,*) "TREE x VERTEX: ",AmpSqGluToFdcSd(gt2, gt3)
  AmpSum2GluToFdcSd = AmpTreeGluToFdcSd + 2._dp*AmpWaveGluToFdcSd + 2._dp*AmpVertexGluToFdcSd
If (OSkinematics) Then
  Call SquareAmp_FtoFS(MGluOS,MFdOS(gt2),MSdOS(gt3),AmpSumGluToFdcSd(:,gt2, gt3),AmpSum2GluToFdcSd(:,gt2, gt3),AmpSqGluToFdcSd(gt2, gt3))
Else
  Call SquareAmp_FtoFS(MGlu,MFd(gt2),MSd(gt3),AmpSumGluToFdcSd(:,gt2, gt3),AmpSum2GluToFdcSd(:,gt2, gt3),AmpSqGluToFdcSd(gt2, gt3))
End if
  Write(*,*) "TREE x (TREE+WAVE+VERTEX): ",AmpSqGluToFdcSd(gt2, gt3)
 End if
If (OSkinematics) Then
  AmpSum2GluToFdcSd = AmpTreeGluToFdcSd
  Call SquareAmp_FtoFS(MGluOS,MFdOS(gt2),MSdOS(gt3),AmpSumGluToFdcSd(:,gt2, gt3),AmpSum2GluToFdcSd(:,gt2, gt3),AmpSqGluToFdcSd(gt2, gt3))
  AmpSqTreeGluToFdcSd(gt2, gt3) = AmpSqGluToFdcSd(gt2, gt3)
  AmpSum2GluToFdcSd = + 2._dp*AmpWaveGluToFdcSd + 2._dp*AmpVertexGluToFdcSd
  Call SquareAmp_FtoFS(MGluOS,MFdOS(gt2),MSdOS(gt3),AmpSumGluToFdcSd(:,gt2, gt3),AmpSum2GluToFdcSd(:,gt2, gt3),AmpSqGluToFdcSd(gt2, gt3))
  AmpSqGluToFdcSd(gt2, gt3) = AmpSqGluToFdcSd(gt2, gt3) + AmpSqTreeGluToFdcSd(gt2, gt3)
Else
  AmpSum2GluToFdcSd = AmpTreeGluToFdcSd
  Call SquareAmp_FtoFS(MGlu,MFd(gt2),MSd(gt3),AmpSumGluToFdcSd(:,gt2, gt3),AmpSum2GluToFdcSd(:,gt2, gt3),AmpSqGluToFdcSd(gt2, gt3))
  AmpSqTreeGluToFdcSd(gt2, gt3) = AmpSqGluToFdcSd(gt2, gt3)
  AmpSum2GluToFdcSd = + 2._dp*AmpWaveGluToFdcSd + 2._dp*AmpVertexGluToFdcSd
  Call SquareAmp_FtoFS(MGlu,MFd(gt2),MSd(gt3),AmpSumGluToFdcSd(:,gt2, gt3),AmpSum2GluToFdcSd(:,gt2, gt3),AmpSqGluToFdcSd(gt2, gt3))
  AmpSqGluToFdcSd(gt2, gt3) = AmpSqGluToFdcSd(gt2, gt3) + AmpSqTreeGluToFdcSd(gt2, gt3)
End if
Else
  AmpSqGluToFdcSd(gt2, gt3) = 0._dp
End if

! Calculate Partial widths
helfactor = 2._dp
If (AmpSqGluToFdcSd(gt2, gt3).le.0._dp) Then
  gP1LGlu(gt1,i4) = 0._dp
Else
If (OSkinematics) Then
  gP1LGlu(gt1,i4) = 1._dp*GammaTPS(MGluOS,MFdOS(gt2),MSdOS(gt3),helfactor*AmpSqGluToFdcSd(gt2, gt3))
Else
  gP1LGlu(gt1,i4) = 1._dp*GammaTPS(MGlu,MFd(gt2),MSd(gt3),helfactor*AmpSqGluToFdcSd(gt2, gt3))
End if
If ((Abs(MRPGluToFdcSd(gt2, gt3)).gt.1.0E-20_dp).or.(Abs(MRGGluToFdcSd(gt2, gt3)).gt.1.0E-20_dp)) Then
  phasespacefactor = 1._dp
End if
 If (DebugLoopDecays) Write(*,*) "virtual", gP1LGlu(gt1,i4)
 ! Adding real corrections
If ((Abs(MRPGluToFdcSd(gt2, gt3)).gt.1.0E-20_dp).or.(Abs(MRGGluToFdcSd(gt2, gt3)).gt.1.0E-20_dp)) Then
 If (.not.OnlyTreeLevelContributions) Then
   If (DebugLoopDecays) Write(*,*) "real", phasespacefactor*1._dp*helfactor*(MRPGluToFdcSd(gt2, gt3) + MRGGluToFdcSd(gt2, gt3))
  gP1LGlu(gt1,i4) = gP1LGlu(gt1,i4) + phasespacefactor*1._dp*helfactor*(MRPGluToFdcSd(gt2, gt3) + MRGGluToFdcSd(gt2, gt3))
   If (DebugLoopDecays) Write(*,*) "sum",  gP1LGlu(gt1,i4)
  End if
End if
End if
i4=i4+1

    End do
  End do
isave = i4
End If
If (.not.CalcLoopDecay_LoopInducedOnly) Then
!----------------
! Fu Conjg(Su)
!----------------

!Tree Level
  If (.not.ExternalZfactors) Then
Call Amplitude_Tree_seesaw1nmssm_GluToFucSu(cplGluFucSuL,cplGluFucSuR,MFu,            &
& MGlu,MSu,MFu2,MGlu2,MSu2,AmpTreeGluToFucSu)

  Else
Call Amplitude_Tree_seesaw1nmssm_GluToFucSu(ZcplGluFucSuL,ZcplGluFucSuR,              &
& MFu,MGlu,MSu,MFu2,MGlu2,MSu2,AmpTreeGluToFucSu)

  End if

!Real Corrections
If (OSkinematics) Then
  If (.not.ExternalZfactors) Then
 ! OS and no Z-factors
Call Gamma_Real_seesaw1nmssm_GluToFucSu(MLambda,em,gs,cplGluFucSuL,cplGluFucSuR,      &
& MFuOS,MGluOS,MSuOS,MRPGluToFucSu,MRGGluToFucSu)

  Else
 ! OS and Z-factors
Call Gamma_Real_seesaw1nmssm_GluToFucSu(MLambda,em,gs,ZcplGluFucSuL,ZcplGluFucSuR,    &
& MFuOS,MGluOS,MSuOS,MRPGluToFucSu,MRGGluToFucSu)

  End if
Else
 ! DR and no Z-factors
  If (.not.ExternalZfactors) Then
Call Gamma_Real_seesaw1nmssm_GluToFucSu(MLambda,em,gs,cplGluFucSuL,cplGluFucSuR,      &
& MFu,MGlu,MSu,MRPGluToFucSu,MRGGluToFucSu)

  Else
 ! DR and Z-factors
Call Gamma_Real_seesaw1nmssm_GluToFucSu(MLambda,em,gs,ZcplGluFucSuL,ZcplGluFucSuR,    &
& MFu,MGlu,MSu,MRPGluToFucSu,MRGGluToFucSu)

  End if
End if

!Self-energy Corrections
Call Amplitude_WAVE_seesaw1nmssm_GluToFucSu(cplGluFucSuL,cplGluFucSuR,ctcplGluFucSuL, &
& ctcplGluFucSuR,MFu,MFu2,MGlu,MGlu2,MSu,MSu2,ZffG,ZfFUL,ZfFUR,ZfSu,AmpWaveGluToFucSu)

!Vertex Corrections
Call Amplitude_VERTEX_seesaw1nmssm_GluToFucSu(MAh,MCha,MChi,MFd,MFu,MGlu,             &
& Mhh,MHpm,MSd,MSu,MVG,MVP,MVWm,MVZ,MAh2,MCha2,MChi2,MFd2,MFu2,MGlu2,Mhh2,               &
& MHpm2,MSd2,MSu2,MVG2,MVP2,MVWm2,MVZ2,cplcFuFuAhL,cplcFuFuAhR,cplAhSucSu,               &
& cplChaFucSdL,cplChaFucSdR,cplChiFucSuL,cplChiFucSuR,cplGluFdcSdL,cplGluFdcSdR,         &
& cplcChaFdcSuL,cplcChaFdcSuR,cplGluFucSuL,cplGluFucSuR,cplcFuFuhhL,cplcFuFuhhR,         &
& cplcFdFuHpmL,cplcFdFuHpmR,cplcFuFuVGL,cplcFuFuVGR,cplcFuFuVPL,cplcFuFuVPR,             &
& cplcFdFuVWmL,cplcFdFuVWmR,cplcFuFuVZL,cplcFuFuVZR,cplGluGluVGL,cplGluGluVGR,           &
& cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,cplhhSucSu,cplSdcHpmcSu,           &
& cplSdcSucVWm,cplSucSuVG,cplSucSuVP,cplSucSuVZ,AmpVertexGluToFucSu)

If (ShiftIRdiv) Then
Call Amplitude_IR_VERTEX_seesaw1nmssm_GluToFucSu(MAh,MCha,MChi,MFd,MFu,               &
& MGlu,Mhh,MHpm,MSd,MSu,MVG,MVP,MVWm,MVZ,MAh2,MCha2,MChi2,MFd2,MFu2,MGlu2,               &
& Mhh2,MHpm2,MSd2,MSu2,MVG2,MVP2,MVWm2,MVZ2,cplcFuFuAhL,cplcFuFuAhR,cplAhSucSu,          &
& cplChaFucSdL,cplChaFucSdR,cplChiFucSuL,cplChiFucSuR,cplGluFdcSdL,cplGluFdcSdR,         &
& cplcChaFdcSuL,cplcChaFdcSuR,cplGluFucSuL,cplGluFucSuR,cplcFuFuhhL,cplcFuFuhhR,         &
& cplcFdFuHpmL,cplcFdFuHpmR,cplcFuFuVGL,cplcFuFuVGR,cplcFuFuVPL,cplcFuFuVPR,             &
& cplcFdFuVWmL,cplcFdFuVWmR,cplcFuFuVZL,cplcFuFuVZR,cplGluGluVGL,cplGluGluVGR,           &
& cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,cplhhSucSu,cplSdcHpmcSu,           &
& cplSdcSucVWm,cplSucSuVG,cplSucSuVP,cplSucSuVZ,AmpVertexIRdrGluToFucSu)

 If (ExternalZfactors) Then
  If (OSkinematics) Then
 ! OS and Z-factors
Call Amplitude_IR_VERTEX_seesaw1nmssm_GluToFucSu(MAhOS,MChaOS,MChiOS,MFdOS,           &
& MFuOS,MGluOS,MhhOS,MHpmOS,MSdOS,MSuOS,MVG,MVP,MVWmOS,MVZOS,MAh2OS,MCha2OS,             &
& MChi2OS,MFd2OS,MFu2OS,MGlu2OS,Mhh2OS,MHpm2OS,MSd2OS,MSu2OS,MVG2,MVP2,MVWm2OS,          &
& MVZ2OS,cplcFuFuAhL,cplcFuFuAhR,cplAhSucSu,cplChaFucSdL,cplChaFucSdR,cplChiFucSuL,      &
& cplChiFucSuR,cplGluFdcSdL,cplGluFdcSdR,cplcChaFdcSuL,cplcChaFdcSuR,ZcplGluFucSuL,      &
& ZcplGluFucSuR,cplcFuFuhhL,cplcFuFuhhR,cplcFdFuHpmL,cplcFdFuHpmR,cplcFuFuVGL,           &
& cplcFuFuVGR,cplcFuFuVPL,cplcFuFuVPR,cplcFdFuVWmL,cplcFdFuVWmR,cplcFuFuVZL,             &
& cplcFuFuVZR,cplGluGluVGL,cplGluGluVGR,cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,          &
& cplcFuGluSuR,cplhhSucSu,cplSdcHpmcSu,cplSdcSucVWm,cplSucSuVG,cplSucSuVP,               &
& cplSucSuVZ,AmpVertexIRosGluToFucSu)

   Else
 ! DR and Z-factors
Call Amplitude_IR_VERTEX_seesaw1nmssm_GluToFucSu(MAh,MCha,MChi,MFd,MFu,               &
& MGlu,Mhh,MHpm,MSd,MSu,MVG,MVP,MVWm,MVZ,MAh2,MCha2,MChi2,MFd2,MFu2,MGlu2,               &
& Mhh2,MHpm2,MSd2,MSu2,MVG2,MVP2,MVWm2,MVZ2,cplcFuFuAhL,cplcFuFuAhR,cplAhSucSu,          &
& cplChaFucSdL,cplChaFucSdR,cplChiFucSuL,cplChiFucSuR,cplGluFdcSdL,cplGluFdcSdR,         &
& cplcChaFdcSuL,cplcChaFdcSuR,ZcplGluFucSuL,ZcplGluFucSuR,cplcFuFuhhL,cplcFuFuhhR,       &
& cplcFdFuHpmL,cplcFdFuHpmR,cplcFuFuVGL,cplcFuFuVGR,cplcFuFuVPL,cplcFuFuVPR,             &
& cplcFdFuVWmL,cplcFdFuVWmR,cplcFuFuVZL,cplcFuFuVZR,cplGluGluVGL,cplGluGluVGR,           &
& cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,cplhhSucSu,cplSdcHpmcSu,           &
& cplSdcSucVWm,cplSucSuVG,cplSucSuVP,cplSucSuVZ,AmpVertexIRosGluToFucSu)

 End if
 Else
  If (OSkinematics) Then
 ! OS and no Z-factors
Call Amplitude_IR_VERTEX_seesaw1nmssm_GluToFucSu(MAhOS,MChaOS,MChiOS,MFdOS,           &
& MFuOS,MGluOS,MhhOS,MHpmOS,MSdOS,MSuOS,MVG,MVP,MVWmOS,MVZOS,MAh2OS,MCha2OS,             &
& MChi2OS,MFd2OS,MFu2OS,MGlu2OS,Mhh2OS,MHpm2OS,MSd2OS,MSu2OS,MVG2,MVP2,MVWm2OS,          &
& MVZ2OS,cplcFuFuAhL,cplcFuFuAhR,cplAhSucSu,cplChaFucSdL,cplChaFucSdR,cplChiFucSuL,      &
& cplChiFucSuR,cplGluFdcSdL,cplGluFdcSdR,cplcChaFdcSuL,cplcChaFdcSuR,cplGluFucSuL,       &
& cplGluFucSuR,cplcFuFuhhL,cplcFuFuhhR,cplcFdFuHpmL,cplcFdFuHpmR,cplcFuFuVGL,            &
& cplcFuFuVGR,cplcFuFuVPL,cplcFuFuVPR,cplcFdFuVWmL,cplcFdFuVWmR,cplcFuFuVZL,             &
& cplcFuFuVZR,cplGluGluVGL,cplGluGluVGR,cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,          &
& cplcFuGluSuR,cplhhSucSu,cplSdcHpmcSu,cplSdcSucVWm,cplSucSuVG,cplSucSuVP,               &
& cplSucSuVZ,AmpVertexIRosGluToFucSu)

   Else
 ! DR and no Z-factors
Call Amplitude_IR_VERTEX_seesaw1nmssm_GluToFucSu(MAh,MCha,MChi,MFd,MFu,               &
& MGlu,Mhh,MHpm,MSd,MSu,MVG,MVP,MVWm,MVZ,MAh2,MCha2,MChi2,MFd2,MFu2,MGlu2,               &
& Mhh2,MHpm2,MSd2,MSu2,MVG2,MVP2,MVWm2,MVZ2,cplcFuFuAhL,cplcFuFuAhR,cplAhSucSu,          &
& cplChaFucSdL,cplChaFucSdR,cplChiFucSuL,cplChiFucSuR,cplGluFdcSdL,cplGluFdcSdR,         &
& cplcChaFdcSuL,cplcChaFdcSuR,cplGluFucSuL,cplGluFucSuR,cplcFuFuhhL,cplcFuFuhhR,         &
& cplcFdFuHpmL,cplcFdFuHpmR,cplcFuFuVGL,cplcFuFuVGR,cplcFuFuVPL,cplcFuFuVPR,             &
& cplcFdFuVWmL,cplcFdFuVWmR,cplcFuFuVZL,cplcFuFuVZR,cplGluGluVGL,cplGluGluVGR,           &
& cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,cplhhSucSu,cplSdcHpmcSu,           &
& cplSdcSucVWm,cplSucSuVG,cplSucSuVP,cplSucSuVZ,AmpVertexIRosGluToFucSu)

 End if
 End if
AmpVertexGluToFucSu = AmpVertexGluToFucSu -  AmpVertexIRdrGluToFucSu! +  AmpVertexIRosGluToFucSu ! Shift added later
End if

 ! Add Z-factors to have external fields on-shell
 If (ExternalZfactors) Then
! Final State 1
AmpWaveZGluToFucSu=0._dp
AmpVertexZGluToFucSu=0._dp
Do gt1=1,3
  Do gt2=1,3
AmpWaveZGluToFucSu(1,gt2,:) = AmpWaveZGluToFucSu(1,gt2,:)+ZRUZUL(gt2,gt1)*AmpWaveGluToFucSu(1,gt1,:)
AmpVertexZGluToFucSu(1,gt2,:)= AmpVertexZGluToFucSu(1,gt2,:)+ZRUZUL(gt2,gt1)*AmpVertexGluToFucSu(1,gt1,:)
AmpWaveZGluToFucSu(2,gt2,:) = AmpWaveZGluToFucSu(2,gt2,:)+ZRUZURc(gt2,gt1)*AmpWaveGluToFucSu(2,gt1,:)
AmpVertexZGluToFucSu(2,gt2,:)= AmpVertexZGluToFucSu(2,gt2,:)+ZRUZURc(gt2,gt1)*AmpVertexGluToFucSu(2,gt1,:)
 End Do
End Do
AmpWaveGluToFucSu=AmpWaveZGluToFucSu
AmpVertexGluToFucSu= AmpVertexZGluToFucSu
! Final State 2
AmpWaveZGluToFucSu=0._dp
AmpVertexZGluToFucSu=0._dp
Do gt1=1,6
  Do gt2=1,6
AmpWaveZGluToFucSu(:,:,gt2) = AmpWaveZGluToFucSu(:,:,gt2)+ZRUZUc(gt2,gt1)*AmpWaveGluToFucSu(:,:,gt1)
AmpVertexZGluToFucSu(:,:,gt2)= AmpVertexZGluToFucSu(:,:,gt2)+ZRUZUc(gt2,gt1)*AmpVertexGluToFucSu(:,:,gt1)
 End Do
End Do
AmpWaveGluToFucSu=AmpWaveZGluToFucSu
AmpVertexGluToFucSu= AmpVertexZGluToFucSu
End if
If (ShiftIRdiv) Then
AmpVertexGluToFucSu = AmpVertexGluToFucSu  +  AmpVertexIRosGluToFucSu
End if

 !Square the amplitude
If (DebugLoopDecays) Then
Write(*,*) "------------------ Glu->Fu conj[Su] -----------------------"
End if
If (.not.SquareFullAmplitudeDecays) Then
 AmpSumGluToFucSu = AmpTreeGluToFucSu
 AmpSum2GluToFucSu = AmpTreeGluToFucSu + 2._dp*AmpWaveGluToFucSu + 2._dp*AmpVertexGluToFucSu
Else
 AmpSumGluToFucSu = AmpTreeGluToFucSu + AmpWaveGluToFucSu + AmpVertexGluToFucSu
 AmpSum2GluToFucSu = AmpTreeGluToFucSu + AmpWaveGluToFucSu + AmpVertexGluToFucSu
End If
If (OnlyTreeLevelContributions) Then
 AmpSumGluToFucSu = AmpTreeGluToFucSu
 AmpSum2GluToFucSu = AmpTreeGluToFucSu
End if
gt1=1
i4 = isave
  Do gt2=1,3
    Do gt3=1,6
If (((OSkinematics).and.(MGluOS.gt.(MFuOS(gt2)+MSuOS(gt3)))).or.((.not.OSkinematics).and.(MGlu.gt.(MFu(gt2)+MSu(gt3))))) Then
 If (DebugLoopDecays) Then
  Write(*,*) gt2, gt3
  AmpSum2GluToFucSu = AmpTreeGluToFucSu
If (OSkinematics) Then
  Call SquareAmp_FtoFS(MGluOS,MFuOS(gt2),MSuOS(gt3),AmpSumGluToFucSu(:,gt2, gt3),AmpSum2GluToFucSu(:,gt2, gt3),AmpSqGluToFucSu(gt2, gt3))
Else
  Call SquareAmp_FtoFS(MGlu,MFu(gt2),MSu(gt3),AmpSumGluToFucSu(:,gt2, gt3),AmpSum2GluToFucSu(:,gt2, gt3),AmpSqGluToFucSu(gt2, gt3))
End if
  Write(*,*) "TREE x TREE: ",AmpSqGluToFucSu(gt2, gt3)
  AmpSum2GluToFucSu = 2._dp*AmpWaveGluToFucSu
If (OSkinematics) Then
  Call SquareAmp_FtoFS(MGluOS,MFuOS(gt2),MSuOS(gt3),AmpSumGluToFucSu(:,gt2, gt3),AmpSum2GluToFucSu(:,gt2, gt3),AmpSqGluToFucSu(gt2, gt3))
Else
  Call SquareAmp_FtoFS(MGlu,MFu(gt2),MSu(gt3),AmpSumGluToFucSu(:,gt2, gt3),AmpSum2GluToFucSu(:,gt2, gt3),AmpSqGluToFucSu(gt2, gt3))
End if
  Write(*,*) "TREE x WAVE: ",AmpSqGluToFucSu(gt2, gt3)
  AmpSum2GluToFucSu = 2._dp*AmpVertexGluToFucSu
If (OSkinematics) Then
  Call SquareAmp_FtoFS(MGluOS,MFuOS(gt2),MSuOS(gt3),AmpSumGluToFucSu(:,gt2, gt3),AmpSum2GluToFucSu(:,gt2, gt3),AmpSqGluToFucSu(gt2, gt3))
Else
  Call SquareAmp_FtoFS(MGlu,MFu(gt2),MSu(gt3),AmpSumGluToFucSu(:,gt2, gt3),AmpSum2GluToFucSu(:,gt2, gt3),AmpSqGluToFucSu(gt2, gt3))
End if
  Write(*,*) "TREE x VERTEX: ",AmpSqGluToFucSu(gt2, gt3)
  AmpSum2GluToFucSu = AmpTreeGluToFucSu + 2._dp*AmpWaveGluToFucSu + 2._dp*AmpVertexGluToFucSu
If (OSkinematics) Then
  Call SquareAmp_FtoFS(MGluOS,MFuOS(gt2),MSuOS(gt3),AmpSumGluToFucSu(:,gt2, gt3),AmpSum2GluToFucSu(:,gt2, gt3),AmpSqGluToFucSu(gt2, gt3))
Else
  Call SquareAmp_FtoFS(MGlu,MFu(gt2),MSu(gt3),AmpSumGluToFucSu(:,gt2, gt3),AmpSum2GluToFucSu(:,gt2, gt3),AmpSqGluToFucSu(gt2, gt3))
End if
  Write(*,*) "TREE x (TREE+WAVE+VERTEX): ",AmpSqGluToFucSu(gt2, gt3)
 End if
If (OSkinematics) Then
  AmpSum2GluToFucSu = AmpTreeGluToFucSu
  Call SquareAmp_FtoFS(MGluOS,MFuOS(gt2),MSuOS(gt3),AmpSumGluToFucSu(:,gt2, gt3),AmpSum2GluToFucSu(:,gt2, gt3),AmpSqGluToFucSu(gt2, gt3))
  AmpSqTreeGluToFucSu(gt2, gt3) = AmpSqGluToFucSu(gt2, gt3)
  AmpSum2GluToFucSu = + 2._dp*AmpWaveGluToFucSu + 2._dp*AmpVertexGluToFucSu
  Call SquareAmp_FtoFS(MGluOS,MFuOS(gt2),MSuOS(gt3),AmpSumGluToFucSu(:,gt2, gt3),AmpSum2GluToFucSu(:,gt2, gt3),AmpSqGluToFucSu(gt2, gt3))
  AmpSqGluToFucSu(gt2, gt3) = AmpSqGluToFucSu(gt2, gt3) + AmpSqTreeGluToFucSu(gt2, gt3)
Else
  AmpSum2GluToFucSu = AmpTreeGluToFucSu
  Call SquareAmp_FtoFS(MGlu,MFu(gt2),MSu(gt3),AmpSumGluToFucSu(:,gt2, gt3),AmpSum2GluToFucSu(:,gt2, gt3),AmpSqGluToFucSu(gt2, gt3))
  AmpSqTreeGluToFucSu(gt2, gt3) = AmpSqGluToFucSu(gt2, gt3)
  AmpSum2GluToFucSu = + 2._dp*AmpWaveGluToFucSu + 2._dp*AmpVertexGluToFucSu
  Call SquareAmp_FtoFS(MGlu,MFu(gt2),MSu(gt3),AmpSumGluToFucSu(:,gt2, gt3),AmpSum2GluToFucSu(:,gt2, gt3),AmpSqGluToFucSu(gt2, gt3))
  AmpSqGluToFucSu(gt2, gt3) = AmpSqGluToFucSu(gt2, gt3) + AmpSqTreeGluToFucSu(gt2, gt3)
End if
Else
  AmpSqGluToFucSu(gt2, gt3) = 0._dp
End if

! Calculate Partial widths
helfactor = 2._dp
If (AmpSqGluToFucSu(gt2, gt3).le.0._dp) Then
  gP1LGlu(gt1,i4) = 0._dp
Else
If (OSkinematics) Then
  gP1LGlu(gt1,i4) = 1._dp*GammaTPS(MGluOS,MFuOS(gt2),MSuOS(gt3),helfactor*AmpSqGluToFucSu(gt2, gt3))
Else
  gP1LGlu(gt1,i4) = 1._dp*GammaTPS(MGlu,MFu(gt2),MSu(gt3),helfactor*AmpSqGluToFucSu(gt2, gt3))
End if
If ((Abs(MRPGluToFucSu(gt2, gt3)).gt.1.0E-20_dp).or.(Abs(MRGGluToFucSu(gt2, gt3)).gt.1.0E-20_dp)) Then
  phasespacefactor = 1._dp
End if
 If (DebugLoopDecays) Write(*,*) "virtual", gP1LGlu(gt1,i4)
 ! Adding real corrections
If ((Abs(MRPGluToFucSu(gt2, gt3)).gt.1.0E-20_dp).or.(Abs(MRGGluToFucSu(gt2, gt3)).gt.1.0E-20_dp)) Then
 If (.not.OnlyTreeLevelContributions) Then
   If (DebugLoopDecays) Write(*,*) "real", phasespacefactor*1._dp*helfactor*(MRPGluToFucSu(gt2, gt3) + MRGGluToFucSu(gt2, gt3))
  gP1LGlu(gt1,i4) = gP1LGlu(gt1,i4) + phasespacefactor*1._dp*helfactor*(MRPGluToFucSu(gt2, gt3) + MRGGluToFucSu(gt2, gt3))
   If (DebugLoopDecays) Write(*,*) "sum",  gP1LGlu(gt1,i4)
  End if
End if
End if
i4=i4+1

    End do
  End do
isave = i4
End If
!----------------
! Chi VG
!----------------

If (LoopInducedDecaysOS) Then

!Self-energy Corrections

!Vertex Corrections
 If (ExternalZfactors) Then
Call Amplitude_VERTEX_seesaw1nmssm_GluToChiVG(MChiOS,MFdOS,MFuOS,MGluOS,              &
& MSdOS,MSuOS,MVG,MChi2OS,MFd2OS,MFu2OS,MGlu2OS,MSd2OS,MSu2OS,MVG2,ZcplChiFdcSdL,        &
& ZcplChiFdcSdR,ZcplChiFucSuL,ZcplChiFucSuR,ZcplcFdChiSdL,ZcplcFdChiSdR,ZcplcFuChiSuL,   &
& ZcplcFuChiSuR,ZcplGluFdcSdL,ZcplGluFdcSdR,ZcplcFdFdVGL,ZcplcFdFdVGR,ZcplGluFucSuL,     &
& ZcplGluFucSuR,ZcplcFuFuVGL,ZcplcFuFuVGR,ZcplcFdGluSdL,ZcplcFdGluSdR,ZcplcFuGluSuL,     &
& ZcplcFuGluSuR,ZcplSdcSdVG,ZcplSucSuVG,AmpVertexGluToChiVG)

 Else
Call Amplitude_VERTEX_seesaw1nmssm_GluToChiVG(MChiOS,MFdOS,MFuOS,MGluOS,              &
& MSdOS,MSuOS,MVG,MChi2OS,MFd2OS,MFu2OS,MGlu2OS,MSd2OS,MSu2OS,MVG2,cplChiFdcSdL,         &
& cplChiFdcSdR,cplChiFucSuL,cplChiFucSuR,cplcFdChiSdL,cplcFdChiSdR,cplcFuChiSuL,         &
& cplcFuChiSuR,cplGluFdcSdL,cplGluFdcSdR,cplcFdFdVGL,cplcFdFdVGR,cplGluFucSuL,           &
& cplGluFucSuR,cplcFuFuVGL,cplcFuFuVGR,cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,           &
& cplcFuGluSuR,cplSdcSdVG,cplSucSuVG,AmpVertexGluToChiVG)

 End if
Else

!Self-energy Corrections

!Vertex Corrections
Call Amplitude_VERTEX_seesaw1nmssm_GluToChiVG(MChi,MFd,MFu,MGlu,MSd,MSu,              &
& MVG,MChi2,MFd2,MFu2,MGlu2,MSd2,MSu2,MVG2,cplChiFdcSdL,cplChiFdcSdR,cplChiFucSuL,       &
& cplChiFucSuR,cplcFdChiSdL,cplcFdChiSdR,cplcFuChiSuL,cplcFuChiSuR,cplGluFdcSdL,         &
& cplGluFdcSdR,cplcFdFdVGL,cplcFdFdVGR,cplGluFucSuL,cplGluFucSuR,cplcFuFuVGL,            &
& cplcFuFuVGR,cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,cplSdcSdVG,            &
& cplSucSuVG,AmpVertexGluToChiVG)

End if

 !Square the amplitude
If (DebugLoopDecays) Then
Write(*,*) "------------------ Glu->Chi VG -----------------------"
End if
If (OnlyTreeLevelContributions) Then
 AmpSumGluToChiVG = 0._dp
 AmpSum2GluToChiVG = 0._dp
Else
 AmpSumGluToChiVG = AmpVertexGluToChiVG + AmpWaveGluToChiVG
 AmpSum2GluToChiVG = AmpVertexGluToChiVG + AmpWaveGluToChiVG
End If
gt1=1
i4 = isave
  Do gt2=1,5
If (((OSkinematics).and.(MGluOS.gt.(MChiOS(gt2)+0.))).or.((.not.OSkinematics).and.(MGlu.gt.(MChi(gt2)+MVG)))) Then
If (OSkinematics) Then
  Call SquareAmp_FtoFV(MGluOS,MChiOS(gt2),0._dp,AmpSumGluToChiVG(:,gt2),AmpSum2GluToChiVG(:,gt2),AmpSqGluToChiVG(gt2))
Else
  Call SquareAmp_FtoFV(MGlu,MChi(gt2),MVG,AmpSumGluToChiVG(:,gt2),AmpSum2GluToChiVG(:,gt2),AmpSqGluToChiVG(gt2))
End if
Else
  AmpSqGluToChiVG(gt2) = 0._dp
End if

! Calculate Partial widths
helfactor = 2._dp
If (AmpSqGluToChiVG(gt2).le.0._dp) Then
  gP1LGlu(gt1,i4) = 0._dp
Else
If (OSkinematics) Then
  gP1LGlu(gt1,i4) = 1._dp*GammaTPS(MGluOS,MChiOS(gt2),0._dp,helfactor*AmpSqGluToChiVG(gt2))
Else
  gP1LGlu(gt1,i4) = 1._dp*GammaTPS(MGlu,MChi(gt2),MVG,helfactor*AmpSqGluToChiVG(gt2))
End if
If ((Abs(MRPGluToChiVG(gt2)).gt.1.0E-20_dp).or.(Abs(MRGGluToChiVG(gt2)).gt.1.0E-20_dp)) Then
  phasespacefactor = 1._dp
End if
 If (DebugLoopDecays) Write(*,*) "virtual", gP1LGlu(gt1,i4)
End if
i4=i4+1

  End do
isave = i4
!----------------
! Glu Ah
!----------------

If (LoopInducedDecaysOS) Then

!Self-energy Corrections

!Vertex Corrections
 If (ExternalZfactors) Then
Call Amplitude_VERTEX_seesaw1nmssm_GluToGluAh(MAhOS,MFdOS,MFuOS,MGluOS,               &
& MSdOS,MSuOS,MAh2OS,MFd2OS,MFu2OS,MGlu2OS,MSd2OS,MSu2OS,ZcplcFdFdAhL,ZcplcFdFdAhR,      &
& ZcplcFuFuAhL,ZcplcFuFuAhR,ZcplAhSdcSd,ZcplAhSucSu,ZcplGluFdcSdL,ZcplGluFdcSdR,         &
& ZcplGluFucSuL,ZcplGluFucSuR,ZcplcFdGluSdL,ZcplcFdGluSdR,ZcplcFuGluSuL,ZcplcFuGluSuR,   &
& AmpVertexGluToGluAh)

 Else
Call Amplitude_VERTEX_seesaw1nmssm_GluToGluAh(MAhOS,MFdOS,MFuOS,MGluOS,               &
& MSdOS,MSuOS,MAh2OS,MFd2OS,MFu2OS,MGlu2OS,MSd2OS,MSu2OS,cplcFdFdAhL,cplcFdFdAhR,        &
& cplcFuFuAhL,cplcFuFuAhR,cplAhSdcSd,cplAhSucSu,cplGluFdcSdL,cplGluFdcSdR,               &
& cplGluFucSuL,cplGluFucSuR,cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,         &
& AmpVertexGluToGluAh)

 End if
Else

!Self-energy Corrections

!Vertex Corrections
Call Amplitude_VERTEX_seesaw1nmssm_GluToGluAh(MAh,MFd,MFu,MGlu,MSd,MSu,               &
& MAh2,MFd2,MFu2,MGlu2,MSd2,MSu2,cplcFdFdAhL,cplcFdFdAhR,cplcFuFuAhL,cplcFuFuAhR,        &
& cplAhSdcSd,cplAhSucSu,cplGluFdcSdL,cplGluFdcSdR,cplGluFucSuL,cplGluFucSuR,             &
& cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,AmpVertexGluToGluAh)

End if

 !Square the amplitude
If (DebugLoopDecays) Then
Write(*,*) "------------------ Glu->Glu Ah -----------------------"
End if
If (OnlyTreeLevelContributions) Then
 AmpSumGluToGluAh = 0._dp
 AmpSum2GluToGluAh = 0._dp
Else
 AmpSumGluToGluAh = AmpVertexGluToGluAh + AmpWaveGluToGluAh
 AmpSum2GluToGluAh = AmpVertexGluToGluAh + AmpWaveGluToGluAh
End If
gt1=1
i4 = isave
    Do gt3=2,3
If (((OSkinematics).and.(MGluOS.gt.(MGluOS+MAhOS(gt3)))).or.((.not.OSkinematics).and.(MGlu.gt.(MGlu+MAh(gt3))))) Then
If (OSkinematics) Then
  Call SquareAmp_FtoFS(MGluOS,MGluOS,MAhOS(gt3),AmpSumGluToGluAh(:,gt3),AmpSum2GluToGluAh(:,gt3),AmpSqGluToGluAh(gt3))
Else
  Call SquareAmp_FtoFS(MGlu,MGlu,MAh(gt3),AmpSumGluToGluAh(:,gt3),AmpSum2GluToGluAh(:,gt3),AmpSqGluToGluAh(gt3))
End if
Else
  AmpSqGluToGluAh(gt3) = 0._dp
End if

! Calculate Partial widths
helfactor = 2._dp
If (AmpSqGluToGluAh(gt3).le.0._dp) Then
  gP1LGlu(gt1,i4) = 0._dp
Else
If (OSkinematics) Then
  gP1LGlu(gt1,i4) = 1._dp*GammaTPS(MGluOS,MGluOS,MAhOS(gt3),helfactor*AmpSqGluToGluAh(gt3))
Else
  gP1LGlu(gt1,i4) = 1._dp*GammaTPS(MGlu,MGlu,MAh(gt3),helfactor*AmpSqGluToGluAh(gt3))
End if
If ((Abs(MRPGluToGluAh(gt3)).gt.1.0E-20_dp).or.(Abs(MRGGluToGluAh(gt3)).gt.1.0E-20_dp)) Then
  phasespacefactor = 1._dp
End if
 If (DebugLoopDecays) Write(*,*) "virtual", gP1LGlu(gt1,i4)
End if
i4=i4+1

    End do
isave = i4
!----------------
! Glu hh
!----------------

If (LoopInducedDecaysOS) Then

!Self-energy Corrections

!Vertex Corrections
 If (ExternalZfactors) Then
Call Amplitude_VERTEX_seesaw1nmssm_GluToGluhh(MFdOS,MFuOS,MGluOS,MhhOS,               &
& MSdOS,MSuOS,MFd2OS,MFu2OS,MGlu2OS,Mhh2OS,MSd2OS,MSu2OS,ZcplGluFdcSdL,ZcplGluFdcSdR,    &
& ZcplcFdFdhhL,ZcplcFdFdhhR,ZcplGluFucSuL,ZcplGluFucSuR,ZcplcFuFuhhL,ZcplcFuFuhhR,       &
& ZcplcFdGluSdL,ZcplcFdGluSdR,ZcplcFuGluSuL,ZcplcFuGluSuR,ZcplhhSdcSd,ZcplhhSucSu,       &
& AmpVertexGluToGluhh)

 Else
Call Amplitude_VERTEX_seesaw1nmssm_GluToGluhh(MFdOS,MFuOS,MGluOS,MhhOS,               &
& MSdOS,MSuOS,MFd2OS,MFu2OS,MGlu2OS,Mhh2OS,MSd2OS,MSu2OS,cplGluFdcSdL,cplGluFdcSdR,      &
& cplcFdFdhhL,cplcFdFdhhR,cplGluFucSuL,cplGluFucSuR,cplcFuFuhhL,cplcFuFuhhR,             &
& cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,cplhhSdcSd,cplhhSucSu,             &
& AmpVertexGluToGluhh)

 End if
Else

!Self-energy Corrections

!Vertex Corrections
Call Amplitude_VERTEX_seesaw1nmssm_GluToGluhh(MFd,MFu,MGlu,Mhh,MSd,MSu,               &
& MFd2,MFu2,MGlu2,Mhh2,MSd2,MSu2,cplGluFdcSdL,cplGluFdcSdR,cplcFdFdhhL,cplcFdFdhhR,      &
& cplGluFucSuL,cplGluFucSuR,cplcFuFuhhL,cplcFuFuhhR,cplcFdGluSdL,cplcFdGluSdR,           &
& cplcFuGluSuL,cplcFuGluSuR,cplhhSdcSd,cplhhSucSu,AmpVertexGluToGluhh)

End if

 !Square the amplitude
If (DebugLoopDecays) Then
Write(*,*) "------------------ Glu->Glu hh -----------------------"
End if
If (OnlyTreeLevelContributions) Then
 AmpSumGluToGluhh = 0._dp
 AmpSum2GluToGluhh = 0._dp
Else
 AmpSumGluToGluhh = AmpVertexGluToGluhh + AmpWaveGluToGluhh
 AmpSum2GluToGluhh = AmpVertexGluToGluhh + AmpWaveGluToGluhh
End If
gt1=1
i4 = isave
    Do gt3=1,3
If (((OSkinematics).and.(MGluOS.gt.(MGluOS+MhhOS(gt3)))).or.((.not.OSkinematics).and.(MGlu.gt.(MGlu+Mhh(gt3))))) Then
If (OSkinematics) Then
  Call SquareAmp_FtoFS(MGluOS,MGluOS,MhhOS(gt3),AmpSumGluToGluhh(:,gt3),AmpSum2GluToGluhh(:,gt3),AmpSqGluToGluhh(gt3))
Else
  Call SquareAmp_FtoFS(MGlu,MGlu,Mhh(gt3),AmpSumGluToGluhh(:,gt3),AmpSum2GluToGluhh(:,gt3),AmpSqGluToGluhh(gt3))
End if
Else
  AmpSqGluToGluhh(gt3) = 0._dp
End if

! Calculate Partial widths
helfactor = 2._dp
If (AmpSqGluToGluhh(gt3).le.0._dp) Then
  gP1LGlu(gt1,i4) = 0._dp
Else
If (OSkinematics) Then
  gP1LGlu(gt1,i4) = 1._dp*GammaTPS(MGluOS,MGluOS,MhhOS(gt3),helfactor*AmpSqGluToGluhh(gt3))
Else
  gP1LGlu(gt1,i4) = 1._dp*GammaTPS(MGlu,MGlu,Mhh(gt3),helfactor*AmpSqGluToGluhh(gt3))
End if
If ((Abs(MRPGluToGluhh(gt3)).gt.1.0E-20_dp).or.(Abs(MRGGluToGluhh(gt3)).gt.1.0E-20_dp)) Then
  phasespacefactor = 1._dp
End if
 If (DebugLoopDecays) Write(*,*) "virtual", gP1LGlu(gt1,i4)
End if
i4=i4+1

    End do
isave = i4
!----------------
! Glu VG
!----------------

If (LoopInducedDecaysOS) Then

!Self-energy Corrections
 If (ExternalZfactors) Then
Call Amplitude_WAVE_seesaw1nmssm_GluToGluVG(ZcplGluGluVGL,ZcplGluGluVGR,              &
& ctcplGluGluVGL,ctcplGluGluVGR,MGluOS,MGlu2OS,MVG,MVG2,ZffG,ZfVG,AmpWaveGluToGluVG)

 Else
Call Amplitude_WAVE_seesaw1nmssm_GluToGluVG(cplGluGluVGL,cplGluGluVGR,ctcplGluGluVGL, &
& ctcplGluGluVGR,MGluOS,MGlu2OS,MVG,MVG2,ZffG,ZfVG,AmpWaveGluToGluVG)

 End if

!Vertex Corrections
 If (ExternalZfactors) Then
Call Amplitude_VERTEX_seesaw1nmssm_GluToGluVG(MFdOS,MFuOS,MGluOS,MSdOS,               &
& MSuOS,MVG,MFd2OS,MFu2OS,MGlu2OS,MSd2OS,MSu2OS,MVG2,ZcplGluFdcSdL,ZcplGluFdcSdR,        &
& ZcplcFdFdVGL,ZcplcFdFdVGR,ZcplGluFucSuL,ZcplGluFucSuR,ZcplcFuFuVGL,ZcplcFuFuVGR,       &
& ZcplGluGluVGL,ZcplGluGluVGR,ZcplcFdGluSdL,ZcplcFdGluSdR,ZcplcFuGluSuL,ZcplcFuGluSuR,   &
& ZcplSdcSdVG,ZcplSucSuVG,ZcplVGVGVG,AmpVertexGluToGluVG)

 Else
Call Amplitude_VERTEX_seesaw1nmssm_GluToGluVG(MFdOS,MFuOS,MGluOS,MSdOS,               &
& MSuOS,MVG,MFd2OS,MFu2OS,MGlu2OS,MSd2OS,MSu2OS,MVG2,cplGluFdcSdL,cplGluFdcSdR,          &
& cplcFdFdVGL,cplcFdFdVGR,cplGluFucSuL,cplGluFucSuR,cplcFuFuVGL,cplcFuFuVGR,             &
& cplGluGluVGL,cplGluGluVGR,cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,         &
& cplSdcSdVG,cplSucSuVG,cplVGVGVG,AmpVertexGluToGluVG)

 End if
Else

!Self-energy Corrections
Call Amplitude_WAVE_seesaw1nmssm_GluToGluVG(cplGluGluVGL,cplGluGluVGR,ctcplGluGluVGL, &
& ctcplGluGluVGR,MGlu,MGlu2,MVG,MVG2,ZffG,ZfVG,AmpWaveGluToGluVG)

!Vertex Corrections
Call Amplitude_VERTEX_seesaw1nmssm_GluToGluVG(MFd,MFu,MGlu,MSd,MSu,MVG,               &
& MFd2,MFu2,MGlu2,MSd2,MSu2,MVG2,cplGluFdcSdL,cplGluFdcSdR,cplcFdFdVGL,cplcFdFdVGR,      &
& cplGluFucSuL,cplGluFucSuR,cplcFuFuVGL,cplcFuFuVGR,cplGluGluVGL,cplGluGluVGR,           &
& cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,cplSdcSdVG,cplSucSuVG,             &
& cplVGVGVG,AmpVertexGluToGluVG)

End if

 !Square the amplitude
If (DebugLoopDecays) Then
Write(*,*) "------------------ Glu->Glu VG -----------------------"
End if
If (OnlyTreeLevelContributions) Then
 AmpSumGluToGluVG = 0._dp
 AmpSum2GluToGluVG = 0._dp
Else
 AmpSumGluToGluVG = AmpVertexGluToGluVG + AmpWaveGluToGluVG
 AmpSum2GluToGluVG = AmpVertexGluToGluVG + AmpWaveGluToGluVG
End If
gt1=1
i4 = isave
If (((OSkinematics).and.(MGluOS.gt.(MGluOS+0.))).or.((.not.OSkinematics).and.(MGlu.gt.(MGlu+MVG)))) Then
If (OSkinematics) Then
  Call SquareAmp_FtoFV(MGluOS,MGluOS,0._dp,AmpSumGluToGluVG(:),AmpSum2GluToGluVG(:),AmpSqGluToGluVG)
Else
  Call SquareAmp_FtoFV(MGlu,MGlu,MVG,AmpSumGluToGluVG(:),AmpSum2GluToGluVG(:),AmpSqGluToGluVG)
End if
Else
  AmpSqGluToGluVG = 0._dp
End if

! Calculate Partial widths
helfactor = 2._dp
If (AmpSqGluToGluVG.le.0._dp) Then
  gP1LGlu(gt1,i4) = 0._dp
Else
If (OSkinematics) Then
  gP1LGlu(gt1,i4) = 3._dp*GammaTPS(MGluOS,MGluOS,0._dp,helfactor*AmpSqGluToGluVG)
Else
  gP1LGlu(gt1,i4) = 3._dp*GammaTPS(MGlu,MGlu,MVG,helfactor*AmpSqGluToGluVG)
End if
If ((Abs(MRPGluToGluVG).gt.1.0E-20_dp).or.(Abs(MRGGluToGluVG).gt.1.0E-20_dp)) Then
  phasespacefactor = 1._dp
End if
 If (DebugLoopDecays) Write(*,*) "virtual", gP1LGlu(gt1,i4)
End if
i4=i4+1

isave = i4
!----------------
! Glu VP
!----------------

If (LoopInducedDecaysOS) Then

!Self-energy Corrections

!Vertex Corrections
 If (ExternalZfactors) Then
Call Amplitude_VERTEX_seesaw1nmssm_GluToGluVP(MFdOS,MFuOS,MGluOS,MSdOS,               &
& MSuOS,MVP,MFd2OS,MFu2OS,MGlu2OS,MSd2OS,MSu2OS,MVP2,ZcplGluFdcSdL,ZcplGluFdcSdR,        &
& ZcplcFdFdVPL,ZcplcFdFdVPR,ZcplGluFucSuL,ZcplGluFucSuR,ZcplcFuFuVPL,ZcplcFuFuVPR,       &
& ZcplcFdGluSdL,ZcplcFdGluSdR,ZcplcFuGluSuL,ZcplcFuGluSuR,ZcplSdcSdVP,ZcplSucSuVP,       &
& AmpVertexGluToGluVP)

 Else
Call Amplitude_VERTEX_seesaw1nmssm_GluToGluVP(MFdOS,MFuOS,MGluOS,MSdOS,               &
& MSuOS,MVP,MFd2OS,MFu2OS,MGlu2OS,MSd2OS,MSu2OS,MVP2,cplGluFdcSdL,cplGluFdcSdR,          &
& cplcFdFdVPL,cplcFdFdVPR,cplGluFucSuL,cplGluFucSuR,cplcFuFuVPL,cplcFuFuVPR,             &
& cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,cplSdcSdVP,cplSucSuVP,             &
& AmpVertexGluToGluVP)

 End if
Else

!Self-energy Corrections

!Vertex Corrections
Call Amplitude_VERTEX_seesaw1nmssm_GluToGluVP(MFd,MFu,MGlu,MSd,MSu,MVP,               &
& MFd2,MFu2,MGlu2,MSd2,MSu2,MVP2,cplGluFdcSdL,cplGluFdcSdR,cplcFdFdVPL,cplcFdFdVPR,      &
& cplGluFucSuL,cplGluFucSuR,cplcFuFuVPL,cplcFuFuVPR,cplcFdGluSdL,cplcFdGluSdR,           &
& cplcFuGluSuL,cplcFuGluSuR,cplSdcSdVP,cplSucSuVP,AmpVertexGluToGluVP)

End if

 !Square the amplitude
If (DebugLoopDecays) Then
Write(*,*) "------------------ Glu->Glu VP -----------------------"
End if
If (OnlyTreeLevelContributions) Then
 AmpSumGluToGluVP = 0._dp
 AmpSum2GluToGluVP = 0._dp
Else
 AmpSumGluToGluVP = AmpVertexGluToGluVP + AmpWaveGluToGluVP
 AmpSum2GluToGluVP = AmpVertexGluToGluVP + AmpWaveGluToGluVP
End If
gt1=1
i4 = isave
If (((OSkinematics).and.(MGluOS.gt.(MGluOS+0.))).or.((.not.OSkinematics).and.(MGlu.gt.(MGlu+MVP)))) Then
If (OSkinematics) Then
  Call SquareAmp_FtoFV(MGluOS,MGluOS,0._dp,AmpSumGluToGluVP(:),AmpSum2GluToGluVP(:),AmpSqGluToGluVP)
Else
  Call SquareAmp_FtoFV(MGlu,MGlu,MVP,AmpSumGluToGluVP(:),AmpSum2GluToGluVP(:),AmpSqGluToGluVP)
End if
Else
  AmpSqGluToGluVP = 0._dp
End if

! Calculate Partial widths
helfactor = 2._dp
If (AmpSqGluToGluVP.le.0._dp) Then
  gP1LGlu(gt1,i4) = 0._dp
Else
If (OSkinematics) Then
  gP1LGlu(gt1,i4) = 1._dp*GammaTPS(MGluOS,MGluOS,0._dp,helfactor*AmpSqGluToGluVP)
Else
  gP1LGlu(gt1,i4) = 1._dp*GammaTPS(MGlu,MGlu,MVP,helfactor*AmpSqGluToGluVP)
End if
If ((Abs(MRPGluToGluVP).gt.1.0E-20_dp).or.(Abs(MRGGluToGluVP).gt.1.0E-20_dp)) Then
  phasespacefactor = 1._dp
End if
 If (DebugLoopDecays) Write(*,*) "virtual", gP1LGlu(gt1,i4)
End if
i4=i4+1

isave = i4
!----------------
! Glu VZ
!----------------

If (LoopInducedDecaysOS) Then

!Self-energy Corrections

!Vertex Corrections
 If (ExternalZfactors) Then
Call Amplitude_VERTEX_seesaw1nmssm_GluToGluVZ(MFdOS,MFuOS,MGluOS,MSdOS,               &
& MSuOS,MVZOS,MFd2OS,MFu2OS,MGlu2OS,MSd2OS,MSu2OS,MVZ2OS,ZcplGluFdcSdL,ZcplGluFdcSdR,    &
& ZcplcFdFdVZL,ZcplcFdFdVZR,ZcplGluFucSuL,ZcplGluFucSuR,ZcplcFuFuVZL,ZcplcFuFuVZR,       &
& ZcplcFdGluSdL,ZcplcFdGluSdR,ZcplcFuGluSuL,ZcplcFuGluSuR,ZcplSdcSdVZ,ZcplSucSuVZ,       &
& AmpVertexGluToGluVZ)

 Else
Call Amplitude_VERTEX_seesaw1nmssm_GluToGluVZ(MFdOS,MFuOS,MGluOS,MSdOS,               &
& MSuOS,MVZOS,MFd2OS,MFu2OS,MGlu2OS,MSd2OS,MSu2OS,MVZ2OS,cplGluFdcSdL,cplGluFdcSdR,      &
& cplcFdFdVZL,cplcFdFdVZR,cplGluFucSuL,cplGluFucSuR,cplcFuFuVZL,cplcFuFuVZR,             &
& cplcFdGluSdL,cplcFdGluSdR,cplcFuGluSuL,cplcFuGluSuR,cplSdcSdVZ,cplSucSuVZ,             &
& AmpVertexGluToGluVZ)

 End if
Else

!Self-energy Corrections

!Vertex Corrections
Call Amplitude_VERTEX_seesaw1nmssm_GluToGluVZ(MFd,MFu,MGlu,MSd,MSu,MVZ,               &
& MFd2,MFu2,MGlu2,MSd2,MSu2,MVZ2,cplGluFdcSdL,cplGluFdcSdR,cplcFdFdVZL,cplcFdFdVZR,      &
& cplGluFucSuL,cplGluFucSuR,cplcFuFuVZL,cplcFuFuVZR,cplcFdGluSdL,cplcFdGluSdR,           &
& cplcFuGluSuL,cplcFuGluSuR,cplSdcSdVZ,cplSucSuVZ,AmpVertexGluToGluVZ)

End if

 !Square the amplitude
If (DebugLoopDecays) Then
Write(*,*) "------------------ Glu->Glu VZ -----------------------"
End if
If (OnlyTreeLevelContributions) Then
 AmpSumGluToGluVZ = 0._dp
 AmpSum2GluToGluVZ = 0._dp
Else
 AmpSumGluToGluVZ = AmpVertexGluToGluVZ + AmpWaveGluToGluVZ
 AmpSum2GluToGluVZ = AmpVertexGluToGluVZ + AmpWaveGluToGluVZ
End If
gt1=1
i4 = isave
If (((OSkinematics).and.(MGluOS.gt.(MGluOS+MVZOS))).or.((.not.OSkinematics).and.(MGlu.gt.(MGlu+MVZ)))) Then
If (OSkinematics) Then
  Call SquareAmp_FtoFV(MGluOS,MGluOS,MVZOS,AmpSumGluToGluVZ(:),AmpSum2GluToGluVZ(:),AmpSqGluToGluVZ)
Else
  Call SquareAmp_FtoFV(MGlu,MGlu,MVZ,AmpSumGluToGluVZ(:),AmpSum2GluToGluVZ(:),AmpSqGluToGluVZ)
End if
Else
  AmpSqGluToGluVZ = 0._dp
End if

! Calculate Partial widths
helfactor = 2._dp
If (AmpSqGluToGluVZ.le.0._dp) Then
  gP1LGlu(gt1,i4) = 0._dp
Else
If (OSkinematics) Then
  gP1LGlu(gt1,i4) = 1._dp*GammaTPS(MGluOS,MGluOS,MVZOS,helfactor*AmpSqGluToGluVZ)
Else
  gP1LGlu(gt1,i4) = 1._dp*GammaTPS(MGlu,MGlu,MVZ,helfactor*AmpSqGluToGluVZ)
End if
If ((Abs(MRPGluToGluVZ).gt.1.0E-20_dp).or.(Abs(MRGGluToGluVZ).gt.1.0E-20_dp)) Then
  phasespacefactor = 1._dp
End if
 If (DebugLoopDecays) Write(*,*) "virtual", gP1LGlu(gt1,i4)
End if
i4=i4+1

isave = i4
End Subroutine OneLoopDecay_Glu

End Module Wrapper_OneLoopDecay_Glu_seesaw1nmssm


text:
      MODULE ti1tv_I
      INTERFACE
!...Translated by Charlotte Froese Fischer
!                       Gediminas Gaigalas  10/05/17
      SUBROUTINE TI1TV(CIIN,NCSF,NCIV,I,L,T,NSHL,CIOUT,NTESTG)
      USE vast_kind_param, ONLY:  DOUBLE
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: ncsf, nciv, nshl, i, L, ntestg
      REAL(DOUBLE), DIMENSION(ncsf, nciv), INTENT(IN) :: ciin
      REAL(DOUBLE), DIMENSION(nshl), INTENT(IN) :: t
      REAL(DOUBLE), DIMENSION(ncsf, nciv), INTENT(OUT) :: ciout
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: sylas/grasp-continuum

max_issues_repo_name: sylas/grasp-continuum

max_forks_repo_name: sylas/grasp-continuum


! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download

!> Prints out a quick table of obs types and counts, overall start and
!> stop times, and metadata strings and counts.
!>
!> You can also get more/different info by running the obs_diag program.
!>
!> This program works on any obs_seq file (e.g. input to perfect_model_obs,
!> input to filter, output from filter).  It can print out the first/last
!> obs times in the file, and the count of each type of obs.
!>
!> If the DART QC data is available it can print out the total number
!> of obs for each QC category, and it can print out the QC counts per
!> obs type.  This program is intended to be a quick look at an unknown
!> obs_seq file to see what types and how many obs it contains. if it
!> has been through filter, it prints what happened to the obs (assimilated,
!> failed forward operator, all rejected by outlier threshold, etc).
!>
!> It has an option to open a separate output file and write single line
!> summaries of the QC values.  This output is intended to be ingested into
!> matlab or other plotting programs. (pro tip: see the matlab 'cell'
!> documentation for how to read in data that contains strings as well as
!> numbers.)
!>
!> A missing piece for this purpose is that obs_seq.final files do not
!> contain the analysis time of the model itself.  This may be what you
!> want to use for the X axis on a plot.  Tim suggested some use of the
!> schedule module might make it possible to generate periodic analysis
!> times although they'd have to be sync'd with the input filenames.
!> This is still TBD.

program obs_assim_count

use        types_mod, only : r8, missing_r8, metadatalength, obstypelength
use    utilities_mod, only : initialize_utilities,            &
                             find_namelist_in_file, check_namelist_read,       &
                             error_handler, E_ERR, E_MSG, nmlfileunit,         &
                             do_nml_file, do_nml_term, get_next_filename,      &
                             open_file, close_file, finalize_utilities,        &
                             file_exist
use     location_mod, only : location_type, get_location, set_location,        &
                             LocationName, read_location, operator(/=),        &
                             write_location
use      obs_def_mod, only : obs_def_type, get_obs_def_time, get_obs_def_type_of_obs,&
                             get_obs_def_location, read_obs_def,               &
                             set_obs_def_time
use     obs_kind_mod, only : max_defined_types_of_obs, get_name_for_type_of_obs, &
                             get_index_for_type_of_obs, read_type_of_obs_table
use time_manager_mod, only : time_type, operator(>), print_time, set_time,     &
                             print_date, set_calendar_type,                    &
                             operator(/=), get_calendar_type, NO_CALENDAR,     &
                             operator(-)
use obs_sequence_mod, only : obs_sequence_type, obs_type, write_obs_seq,       &
                             init_obs, assignment(=), get_obs_def,             &
                             init_obs_sequence, static_init_obs_sequence,      &
                             read_obs_seq_header, read_obs_seq, get_num_obs,   &
                             get_first_obs, get_last_obs, get_next_obs,        &
                             insert_obs_in_seq, get_num_copies, get_num_qc,    &
                             get_copy_meta_data, get_qc_meta_data,             &
                             set_copy_meta_data, set_qc_meta_data,             &
                             destroy_obs, destroy_obs_sequence,                &
                             delete_seq_head, delete_seq_tail,                 &
                             get_num_key_range, get_obs_key, get_qc,           &
                             copy_partial_obs, get_next_obs_from_key,          &
                             get_obs_def, set_obs_def

implicit none

character(len=*), parameter :: source   = 'obs_assim_count.f90'

type(obs_sequence_type) :: seq_in
integer                 :: size_seq_in
integer                 :: num_copies_in, num_qc_in
integer                 :: iunit, io, ifile
integer                 :: max_num_obs, file_id
character(len = 32)     :: read_format
character(len = 256)    :: filename_in
logical                 :: pre_I_format, cal
character(len = 512)    :: msgstring, msgstring1, msgstring2

! could go into namelist if you wanted more control
integer, parameter      :: print_every = 5000

! lazy, pick big number.  make it bigger if too small.
integer, parameter :: max_obs_input_types = 800

! 8 is now failed vert convert
integer, parameter :: MAX_DART_QC = 8

integer, parameter :: QC_STR_LEN = 20

character(len=QC_STR_LEN) :: qc_meaning(0:MAX_DART_QC) = &
 (/ "Assimilated_OK      ", &
    "Evaluated_OK        ", &
    "Prior_only_Assim_OK ", &
    "Prior_only_Eval_OK  ", &
    "Forward_Op_failed   ", &
    "Namelist_Excluded   ", &
    "Bad_incoming_QC     ", &
    "Outlier_Failed      ", &
    "Bad_Vert_Convert    " /)

! if the user gives us an output file and asks for counts_only,
! print some stuff to the log & stdout too, but open this file
! and print the counts there as well.
integer :: output_unit = -1

! first time through a loop
logical :: first_time

!----------------------------------------------------------------
! Namelist input with default values

character(len = 256) :: obs_sequence_name = ''
character(len = 256) :: obs_sequence_list = ''
character(len = 256) :: output_file = ''

logical :: stats_by_obs_type = .false.
logical :: counts_only       = .false.
logical :: print_metadata    = .false.
logical :: print_time_info   = .true.

character(len=32) :: calendar    = 'Gregorian'

namelist /obs_assim_count_nml/ &
         obs_sequence_name, obs_sequence_list, output_file, &
         calendar, stats_by_obs_type, counts_only, print_metadata, &
         print_time_info

!----------------------------------------------------------------
! Start of the program:
!
! Process each input observation sequence file in turn, collecting
! and printing stats from each input file.
!----------------------------------------------------------------

call setup()

! Read the namelist entry
call find_namelist_in_file("input.nml", "obs_assim_count_nml", iunit)
read(iunit, nml = obs_assim_count_nml, iostat = io)
call check_namelist_read(iunit, io, "obs_assim_count_nml")

! Record the namelist values used for the run ...
if (do_nml_file()) write(nmlfileunit, nml=obs_assim_count_nml)
if (do_nml_term()) write(     *     , nml=obs_assim_count_nml)

! the default is a gregorian calendar.  if you are using a different type
! set it in the namelist.  this only controls how it prints out the first
! and last timestamps in the obs_seq files.
call set_calendar_type(calendar)

! set a logical to see if we have a calendar or not
cal = (get_calendar_type() /= NO_CALENDAR)

! Check the user input for sanity
if ((obs_sequence_name /= '') .and. (obs_sequence_list /= '')) then
   write(msgstring1,*)'specify "obs_sequence_name" or "obs_sequence_list"'
   write(msgstring2,*)'set other to an empty string ... i.e. ""'
   call error_handler(E_ERR, 'obs_assim_count', msgstring1, &
                   source, text2=msgstring2)
endif

if (counts_only) &
   call error_handler(E_MSG, 'obs_assim_count', 'printing summary counts only')

if (output_file /= '') &
   output_unit = open_file(output_file, action='write', form='formatted')

! if you add anything to the namelist, you can process it here.

! end of namelist processing and setup

ifile = 0
ObsFileLoop : do      ! until we run out of names

   ifile = ifile + 1

   if (obs_sequence_list == '') then
      if (ifile > 1) exit ObsFileLoop
      filename_in = obs_sequence_name
   else
      filename_in = get_next_filename(obs_sequence_list,ifile)
      if (filename_in == '') exit ObsFileLoop
   endif

   if ( .not. file_exist(filename_in) ) then
      write(msgstring1,*)'cannot open ', trim(filename_in)
      call error_handler(E_ERR,'obs_assim_count:',msgstring1,source)
   endif

   ! Read in information about observation sequence so we can allocate
   ! observations. We need info about how many copies, qc values, etc.

   call read_obs_seq_header(filename_in, num_copies_in, num_qc_in, &
      size_seq_in, max_num_obs, file_id, read_format, pre_I_format, &
      close_the_file = .true.)

   if (max_num_obs == 0) then
      write(msgstring,*) 'No obs in input sequence file ', trim(filename_in)
      call error_handler(E_MSG,'obs_assim_count',msgstring)
      cycle ObsFileLoop
   endif

   write(msgstring, *) 'Starting to process input sequence file: '
   write(msgstring1,*)  trim(filename_in)
   call error_handler(E_MSG,'obs_assim_count',msgstring, &
                      text2=msgstring1)

   call read_obs_seq(filename_in, 0, 0, 0, seq_in)

   ! sanity check - ensure the linked list times are in increasing time order
   call validate_obs_seq_time(seq_in, filename_in)

   ! the counting up is done here now.
   call compute_and_print_obs_seq_info(seq_in, filename_in, output_unit)

   ! clean up
   call destroy_obs_sequence(seq_in)

enddo ObsFileLoop

! see comment where the output_unit variable is declared
if (output_file /= '') &
   call close_file(output_unit)

call shutdown()

!---------------------------------------------------------------------
! end of main program.
!---------------------------------------------------------------------

contains

!---------------------------------------------------------------------
subroutine setup()

! Initialize modules used that require it
call initialize_utilities('obs_assim_count')
call static_init_obs_sequence()

end subroutine setup

!---------------------------------------------------------------------
subroutine shutdown()

call finalize_utilities('obs_assim_count')

end subroutine shutdown

!---------------------------------------------------------------------
! you can get more info by running the obs_diag program, but this
! prints out a quick table of obs types and counts, overall start and
! stop times, and metadata strings and counts.

! and this one counts up, if there is a 'posterior ensemble mean' copy,
! how many are missing_r8 and how many are not.   it could also count
! up the 'DART quality control' settings?  start with the latter for now.

subroutine compute_and_print_obs_seq_info(seq_in, filename, ounit)

type(obs_sequence_type), intent(in) :: seq_in
character(len=*),        intent(in) :: filename
integer,                 intent(in) :: ounit

type(obs_type)     :: obs, next_obs
type(obs_def_type) :: this_obs_def
logical            :: is_there_one, is_this_last
integer            :: size_seq_in
integer            :: i, qc_int, obt
integer            :: this_obs_type
integer            :: type_count(0:max_defined_types_of_obs), identity_count, qc_count(0:MAX_DART_QC), qcindex
integer            :: qc_count_by_type(0:MAX_DART_QC, 0:max_defined_types_of_obs)
real(r8)           :: qcval(1)
character(len=32)  :: this_obs_name

! Initialize counters
type_count(:) = 0
identity_count = 0
qc_count(:) = 0
qc_count_by_type(:,:) = 0

size_seq_in = get_num_obs(seq_in)
if (size_seq_in == 0) then
   msgstring = 'Obs_seq file '//trim(filename)//' is empty, skipping.'
   call error_handler(E_MSG,'obs_assim_count',msgstring)
   return
endif

! Initialize individual observation variables
call init_obs(     obs, get_num_copies(seq_in), get_num_qc(seq_in))
call init_obs(next_obs, get_num_copies(seq_in), get_num_qc(seq_in))

! find the dart qc copy, if there is one
qcindex = get_dartqc_index(seq_in, filename)

! blank line
call error_handler(E_MSG,'',' ')

write(msgstring,*) 'Processing sequence file: ', trim(filename)
call error_handler(E_MSG,'',msgstring)

if (print_metadata) &
   call print_seq_metadata(seq_in, filename)

!-------------------------------------------------------------
! Start to process obs from seq_in
!--------------------------------------------------------------
is_there_one = get_first_obs(seq_in, obs)

if ( .not. is_there_one )  then
   write(msgstring,*)'no first observation in ',trim(filename)
   call error_handler(E_MSG,'obs_assim_count', msgstring)
endif

! process it here
is_this_last = .false.

! blank line
call error_handler(E_MSG, '', ' ')

call get_obs_def(obs, this_obs_def)
if (print_time_info) then
   call print_time(get_obs_def_time(this_obs_def), ' First obs time: ')
   ! should, but does not work with NO_CALENDAR
   if (cal) call print_date(get_obs_def_time(this_obs_def), ' First obs date: ')
endif

! collect data counts
ObsLoop : do while ( .not. is_this_last)

   call get_obs_def(obs, this_obs_def)
   this_obs_type = get_obs_def_type_of_obs(this_obs_def)
   if (this_obs_type < 0) then
      identity_count = identity_count + 1
   else
      type_count(this_obs_type) = type_count(this_obs_type) + 1
   endif
   if (qcindex > 0) then
      call get_qc(obs, qcval, qcindex)
      qc_int = nint(qcval(1))
      qc_count(qc_int) = qc_count(qc_int) + 1
      if (this_obs_type >= 0) &
         qc_count_by_type(qc_int,this_obs_type) = qc_count_by_type(qc_int,this_obs_type) + 1
   endif

   call get_next_obs(seq_in, obs, next_obs, is_this_last)
   if (.not. is_this_last) then
      obs = next_obs
   else
      if (print_time_info) then
         call print_time(get_obs_def_time(this_obs_def), '  Last obs time: ')
         if (cal) call print_date(get_obs_def_time(this_obs_def), '  Last obs date: ')
         call error_handler(E_MSG, '', ' ')
      endif
   endif

enddo ObsLoop

! how to output - depends on 'counts_only', 'stats_by_obs_type' settings and
! whether input file has DART QCs or not.

! FIXME: this is very messy.  get it working and then refactor so there
! isn't such a deep set of nested if/else/thens.  the cases are:
!  - are you asking for counts only or a full summary?
!  - are you asking for stats by obs type or only totals?
!  - does this file contain DART QC info?
! which makes 6 different code blocks below, each with slightly
! different output strings and formats.

if (counts_only) then
   if (qcindex > 0) then
      if (stats_by_obs_type) then
         first_time = .true.
         do obt = 0, max_defined_types_of_obs
            if (sum(qc_count_by_type(:,obt)) <= 0) cycle
            this_obs_name = get_name_for_type_of_obs(obt)

            ! print a header, and then try to make the rest of the lines have
            ! all the info on a single line.  (easier to grep, and to read into
            ! matlab for stats, i hope.)
            if (first_time) call print_header('DART QC results:', &
               'Obs_type, QC value, QC meaning, count_of_that_QC_value, total_obs_of_this_type', &
               first_time)

            do i=0, MAX_DART_QC
               if (qc_count_by_type(i,obt) > 0) then
                  write(msgstring, '(a34,a4,i1,a22,2(i12))') trim(this_obs_name), ' QC ', i, qc_meaning(i), &
                                                          qc_count_by_type(i,obt), sum(qc_count_by_type(:,obt))
                  call error_handler(E_MSG, '', msgstring)
                  if (ounit >= 0) write(ounit, '(A)') trim(msgstring)
               endif
            enddo
         enddo
      else
         first_time = .true.

         ! print a header, and then try to make the rest of the lines have
         ! all the info on a single line.  (easier to grep, and to read into
         ! matlab for stats, i hope.)
         if (first_time) call print_header('DART QC results:', &
            'QC value, count_of_that_QC_value, total_obs', &
            first_time)

         do i=0, MAX_DART_QC
            if (qc_count(i) > 0) then
               write(msgstring, '(a22,a4,i1,2(i12))') qc_meaning(i), ' QC ', i, qc_count(i), sum(qc_count(:))
               call error_handler(E_MSG, '', msgstring)
               if (ounit >= 0) write(ounit, '(A)') trim(msgstring)
            endif
         enddo
      endif
   else
      if (stats_by_obs_type) then
         first_time = .true.
         do obt = 0, max_defined_types_of_obs
            if (type_count(obt) <= 0) cycle
            this_obs_name = get_name_for_type_of_obs(obt)

            ! print a header, and then try to make the rest of the lines have
            ! all the info on a single line.  (easier to grep, and to read into
            ! matlab for stats, i hope.)
            if (first_time) call print_header('DART Obs results:', &
               'Obs_type, total_obs_of_this_type', &
               first_time)

            write(msgstring, '(a34,1(i12))') trim(this_obs_name), type_count(obt)
            call error_handler(E_MSG, '', msgstring)
            if (ounit >= 0) write(ounit, '(A)') trim(msgstring)
         enddo
      else
         first_time = .true.
         call print_header('DART QC results:', &
            'QC value, count_of_that_QC_value, total_obs', &
            first_time)

         do i=0, MAX_DART_QC
            if (qc_count(i) > 0) then
               write(msgstring, '(3(i12))') i, qc_count(i), sum(qc_count(:))
               call error_handler(E_MSG, '', msgstring)
               if (ounit >= 0) write(ounit, '(A)') trim(msgstring)
            endif
         enddo
      endif
   endif
else

   write(msgstring, *) 'Number of obs processed  :          ', size_seq_in
   call error_handler(E_MSG, '', msgstring)
   write(msgstring, *) '---------------------------------------------------------'
   call error_handler(E_MSG, '', msgstring)
   do i = 1, max_defined_types_of_obs
      if (type_count(i) > 0) then
         write(msgstring, '(a32,i12,a)') trim(get_name_for_type_of_obs(i)), &
                                        type_count(i), ' obs'
         call error_handler(E_MSG, '', msgstring)
      endif
   enddo
   if (identity_count > 0) then
      write(msgstring, '(a32,i12,a)') 'Identity observations', &
                                     identity_count, ' obs'
      call error_handler(E_MSG, '', msgstring)
   endif
   if (qcindex > 0) then
      call error_handler(E_MSG, '', ' ')
      write(msgstring, *) 'DART QC results for Total obs: ', sum(qc_count(:))
      call error_handler(E_MSG, '', msgstring)
      do i=0, MAX_DART_QC
         if (qc_count(i) > 0) then
            write(msgstring, '(a22,a4,i1,i12)') qc_meaning(i), ' QC ', i, qc_count(i)
            call error_handler(E_MSG, '', msgstring)
         endif
      enddo

      if (stats_by_obs_type) then
         do obt = 0, max_defined_types_of_obs
            if (sum(qc_count_by_type(:,obt)) <= 0) cycle
            this_obs_name = get_name_for_type_of_obs(obt)

            call error_handler(E_MSG, '', ' ')
            write(msgstring, *) 'Total '//trim(this_obs_name)//' obs: ', sum(qc_count_by_type(:,obt))
            call error_handler(E_MSG, '', msgstring)
            do i=0, MAX_DART_QC
               if (qc_count_by_type(i,obt) > 0) then
                  write(msgstring, '(a22,a4,i1,i12)') qc_meaning(i), ' QC ', i, &
                                                      qc_count_by_type(i,obt)
                  call error_handler(E_MSG, '', msgstring)
               endif
            enddo
         enddo
      endif
   endif

   ! another blank line
   call error_handler(E_MSG, '', ' ')
endif

! Time to clean up

call destroy_obs(     obs)
call destroy_obs(next_obs)

end subroutine compute_and_print_obs_seq_info

!---------------------------------------------------------------------

! something like this eventually belongs in the obs_seq_mod code, but
! for now try it out here.  we just fixed a hole in the interactive create
! routine which would silently let you create out-of-time-order
! linked lists, which gave no errors but didn't assimilate the
! right obs at the right time when running filter.   this runs
! through the times in the entire sequence, ensuring they are
! monotonically increasing in time.  this should help catch any
! bad files which were created with older versions of code, or
! created by users with python...

subroutine validate_obs_seq_time(seq, filename)

type(obs_sequence_type), intent(in) :: seq
character(len=*),        intent(in) :: filename

type(obs_type)          :: obs, next_obs
type(obs_def_type)      :: this_obs_def
logical                 :: is_there_one, is_this_last
integer                 :: size_seq, obs_count
integer                 :: key
type(time_type)         :: last_time, this_time

! make sure there are obs left to process before going on.
size_seq = get_num_obs(seq)
if (size_seq == 0) then
   msgstring = 'Obs_seq file '//trim(filename)//' is empty.'
   call error_handler(E_MSG,'obs_assim_count:validate',msgstring)
   return
endif

! Initialize individual observation variables
call init_obs(     obs, get_num_copies(seq), get_num_qc(seq))
call init_obs(next_obs, get_num_copies(seq), get_num_qc(seq))

obs_count = 0

!-------------------------------------------------------------
! Start to process obs from seq
!--------------------------------------------------------------
is_there_one = get_first_obs(seq, obs)

! we already tested for 0 obs above, so there should be a first obs here.
if ( .not. is_there_one )  then
   write(msgstring,*)'no first obs in sequence ' // trim(filename)
   call error_handler(E_ERR,'obs_assim_count:validate', msgstring, source)
   return
endif

is_this_last = .false.
last_time = set_time(0, 0)
ObsLoop : do while ( .not. is_this_last)

   call get_obs_def(obs, this_obs_def)
   this_time = get_obs_def_time(this_obs_def)

   if (last_time > this_time) then
      ! bad time order of observations in linked list
      call print_time(last_time, ' previous timestamp: ')
      if (cal) call print_date(last_time, '      calendar date: ')
      call print_time(this_time, '     next timestamp: ')
      if (cal) call print_date(this_time, '      calendar date: ')

      key = get_obs_key(obs)
      write(msgstring1,*)'obs number ', key, ' has earlier time than previous obs'
      write(msgstring2,*)'observations must be in increasing time order, file ' // trim(filename)
      call error_handler(E_ERR,'obs_assim_count:validate', msgstring2, &
                         source, text2=msgstring1)
   endif

   last_time = this_time
   obs_count = obs_count + 1

   call get_next_obs(seq, obs, next_obs, is_this_last)
   if (.not. is_this_last) obs = next_obs

enddo ObsLoop

! clean up
call destroy_obs(     obs)
call destroy_obs(next_obs)

! technically not a time validation, but easy to check.  obs_count should never
! be larger than size_seq - that's a fatal error.  obs_count < size_seq would
! suggest there are obs in the file that aren't part of the linked list.
! this does not necessarily indicate a fatal error but it's not a common
! situation and might indicate someone should check on the file.
if (obs_count /= size_seq) then
   write(msgstring,*) 'input sequence ', trim(filename)
   call error_handler(E_MSG,'obs_assim_count:validate', msgstring)

   write(msgstring,*) 'total obs in file: ', size_seq, '  obs in linked list: ', obs_count
   if (obs_count > size_seq) then
      ! this is a fatal error
      write(msgstring1,*) 'linked list obs_count > total size_seq, should not happen'
      call error_handler(E_ERR,'obs_assim_count:validate', msgstring, &
                         source, text2=msgstring1)
   else
      ! just warning msg
      write(msgstring1,*) 'only observations in linked list will be processed'
      call error_handler(E_MSG,'obs_assim_count:validate', msgstring, &
                         source, text2=msgstring1)
   endif
endif

end subroutine validate_obs_seq_time

!---------------------------------------------------------------------
! print out the metadata strings, trimmed

subroutine print_seq_metadata(seq, fname)

type(obs_sequence_type), intent(in) :: seq
character(len=*) :: fname

integer :: num_copies , num_qc, i
character(len=metadatalength) :: str
character(len=255) :: msgstring3

num_copies = get_num_copies(seq)
num_qc     = get_num_qc(    seq)

if ( num_copies < 0 .or. num_qc < 0 ) then
   write(msgstring3,*)' illegal copy or obs count in file '//trim(fname)
   call error_handler(E_ERR, 'obs_assim_count', msgstring3, source)
endif

MetaDataLoop : do i=1, num_copies
   str = get_copy_meta_data(seq,i)

   write(msgstring,*)'Data Metadata: ',trim(str)
   call error_handler(E_MSG, '', msgstring)

enddo MetaDataLoop

QCMetaData : do i=1, num_qc
   str = get_qc_meta_data(seq,i)

   write(msgstring,*)'  QC Metadata: ', trim(str)
   call error_handler(E_MSG, '', msgstring)

enddo QCMetaData

end subroutine print_seq_metadata

!---------------------------------------------------------------------
! return the index number of the dart qc copy (-1 if none)

function get_dartqc_index(seq, fname)

type(obs_sequence_type), intent(in) :: seq
character(len=*) :: fname
integer :: get_dartqc_index

integer :: num_qc, i
character(len=metadatalength) :: str
character(len=255) :: msgstring3

num_qc = get_num_qc(seq)

if ( num_qc < 0 ) then
   write(msgstring3,*)' illegal qc metadata count in file '//trim(fname)
   call error_handler(E_ERR, 'obs_assim_count', msgstring3, source)
endif

QCMetaData : do i=1, num_qc
   str = get_qc_meta_data(seq,i)

   if (str == 'DART quality control') then
      get_dartqc_index = i
      return
   endif

enddo QCMetaData

get_dartqc_index = -1

end function get_dartqc_index

!---------------------------------------------------------------------

subroutine print_header(string1, string2, first)

character(len=*), intent(in)    :: string1
character(len=*), intent(in)    :: string2
logical,          intent(inout) :: first

if (.not. first) return

call error_handler(E_MSG, '', ' ')
call error_handler(E_MSG, '', string1)
call error_handler(E_MSG, '', string2)
call error_handler(E_MSG, '', ' ')

first = .false.

end subroutine print_header

!---------------------------------------------------------------------

end program obs_assim_count


text:
      subroutine nplotter_ttbar(p,wt,wt2,switch)
c--- Variable passed in to this routine:
c
c---      p:  4-momenta of particles in the format p(i,4)
c---          with the particles numbered according to the input file
c---          and components labelled by (px,py,pz,E)
c
c---     wt:  weight of this event
c
c---    wt2:  weight^2 of this event
c
c--- switch:  an integer equal to 0 or 1, depending on the type of event
c---                0  --> lowest order, virtual or real radiation
c---                1  --> counterterm for real radiation
      implicit none
      include 'vegas_common.f'
      include 'constants.f'
      include 'histo.f'
      include 'masses.f'
      include 'plabel.f'
      include 'process.f'
      include 'outputflags.f'
      double precision p(mxpart,4),wt,wt2,yrap,pt,ptjet(mxpart),
     & pttwo,mll,ptl,yl,ptttb,yttb,dot,tiny,
     & ptt,pttb,yt,ytb,mttb,mwp,mwm,mlb
      double precision plep(4),ptop(4),pleprest(4),ptoprest(4),mj1j2
      integer switch,n,nplotmax,j,
     & jetindex(mxpart),iorder(mxpart),ijet
      character*4 tag
      logical dilepton,failed
      parameter(tiny=1d-8)
      logical, save::first=.true.
      common/nplotmax/nplotmax
ccccc!$omp threadprivate(first,/nplotmax/,mj1j2)

************************************************************************
*                                                                      *
*     INITIAL BOOKKEEPING                                              *
*                                                                      *
************************************************************************

c--- determine whether or not process represents dilepton channel
      if ((case .eq. 'tt_bbl') .or. (case .eq. 'tt_ldk')
     &.or.(case .eq. 'tt_bbu') .or. (case .eq. 'tt_udk') ) then
        dilepton=.true.
      else
        dilepton=.false.
      endif

      if (first) then
c--- Initialize histograms, without computing any quantities; instead
c--- set them to dummy values
        tag='book'
        mj1j2=0d0
        goto 99
      else
c--- Add event in histograms
        tag='plot'
      endif

************************************************************************
*                                                                      *
*     DEFINITIONS OF QUANTITIES TO PLOT                                *
*                                                                      *
************************************************************************

       ijet=0
       do j=3,9
         if ((plabel(j) .eq. 'pp') .or. (plabel(j) .eq. 'bq')
     &   .or.(plabel(j) .eq. 'ba')) then
           if (p(j,4) .gt. tiny) then
             ijet=ijet+1
             jetindex(ijet)=j
             ptjet(ijet)=pt(j,p)
           endif
         endif
       enddo
       call arraysort(ijet,ptjet,iorder)

       if (ijet .ge. 2) then
c--- Two hardest jets are now indexed by
c--- jetindex(iorder(1)) and jetindex(iorder(2))
       mj1j2=((p(jetindex(iorder(1)),4)+p(jetindex(iorder(2)),4))**2
     &       -(p(jetindex(iorder(1)),1)+p(jetindex(iorder(2)),1))**2
     &       -(p(jetindex(iorder(1)),2)+p(jetindex(iorder(2)),2))**2
     &       -(p(jetindex(iorder(1)),3)+p(jetindex(iorder(2)),3))**2)
       mj1j2=dsqrt(max(mj1j2,0d0))
       else
       mj1j2=-1d0
       endif

************************************************************************
*                                                                      *
*     FILL HISTOGRAMS                                                  *
*                                                                      *
************************************************************************

c--- Call histogram routines
   99 continue

c--- Book and fill ntuple if that option is set, remembering to divide
c--- by # of iterations now that is handled at end for regular histograms
      if (creatent .eqv. .true.) then
        call bookfill(tag,p,wt/dfloat(itmx))
      endif

c--- "n" will count the number of histograms
      n=nextnplot

c--- Syntax of "bookplot" routine is:
c
c---   call bookplot(n,tag,titlex,var,wt,wt2,xmin,xmax,dx,llplot)
c
c---        n:  internal number of histogram
c---      tag:  "book" to initialize histogram, "plot" to fill
c---   titlex:  title of histogram
c---      var:  value of quantity being plotted
c---       wt:  weight of this event (passed in)
c---      wt2:  weight of this event (passed in)
c---     xmin:  lowest value to bin
c---     xmax:  highest value to bin
c---       dx:  bin width
c---   llplot:  equal to "lin"/"log" for linear/log scale

c--- This section tailored for comparison with CMS PAS TOP-11-013

      if (dilepton) then
c--- dilepton=specific plots

c--- pt(l+ and l-)
      call bookplot(n,tag,'pt(l+ and l-)',pt(4,p),wt,wt2,
     & 0d0,400d0,10d0,'log')
      call bookplot(n,tag,'pt(l+ and l-)',pt(7,p),wt,wt2,
     & 0d0,400d0,10d0,'log')
      n=n+1

c--- eta(l+ and l-)
      call bookplot(n,tag,'eta(l+ and l-)',yrap(4,p),wt,wt2,
     & -2.4d0,2.4d0,0.2d0,'lin')
      call bookplot(n,tag,'eta(l+ and l-)',yrap(7,p),wt,wt2,
     & -2.4d0,2.4d0,0.2d0,'lin')
      n=n+1

c--- pt(l+,l-)
      call bookplot(n,tag,'pt(l+,l-)',pttwo(4,7,p),wt,wt2,
     & 0d0,400d0,10d0,'log')
      n=n+1

c--- m(l+,l-)
      mll=dsqrt(max(2d0*dot(p,4,7),0d0))
      call bookplot(n,tag,'m(l+,l-)',mll,wt,wt2,
     & 0d0,400d0,10d0,'log')
      n=n+1

      else

c--- lepton+jets-specific plots
      if (plabel(4) .eq. 'ea') then
        ptl=pt(4,p)
        yl=yrap(4,p)
      else
        ptl=pt(7,p)
        yl=yrap(7,p)
      endif

c--- pt(lepton)
      call bookplot(n,tag,'pt(lepton)',ptl,wt,wt2,
     & 0d0,200d0,10d0,'log')
      n=n+1

c--- eta(lepton)
      call bookplot(n,tag,'eta(lepton)',yl,wt,wt2,
     & -2.2d0,2.2d0,0.2d0,'lin')
      n=n+1

      endif

c--- these plots are common to dilepton and lepton+jets channels
      call topreconstruct(p,failed,
     & ptt,yt,pttb,ytb,ptttb,yttb,mttb,mwp,mwm)

      if ((failed) .and. (first .eqv. .false.)) then
c--- make sure to increment by the number of plots in the 'else' section
      n=n+15

      else

c--- pt(t)
      call bookplot(n,tag,'pt(t)',ptt,wt,wt2,
     & 0d0,400d0,10d0,'log')
      n=n+1

c--- pt(tbar)
      call bookplot(n,tag,'pt(tbar)',pttb,wt,wt2,
     & 0d0,400d0,10d0,'log')
      n=n+1

c--- pt(t and tbar)
      call bookplot(n,tag,'pt(t and tbar)',ptt,wt,wt2,
     & 0d0,400d0,10d0,'log')
      call bookplot(n,tag,'pt(t and tbar)',pttb,wt,wt2,
     & 0d0,400d0,10d0,'log')
      n=n+1

c--- y(t)
      call bookplot(n,tag,'y(t)',yt,wt,wt2,
     & -2.6d0,2.6d0,0.2d0,'lin')
      n=n+1

c--- y(tbar)
      call bookplot(n,tag,'y(tbar)',ytb,wt,wt2,
     & -2.6d0,2.6d0,0.2d0,'lin')
      n=n+1

c--- y(t and tbar)
      call bookplot(n,tag,'y(t and tbar)',yt,wt,wt2,
     & -2.6d0,2.6d0,0.2d0,'lin')
      call bookplot(n,tag,'y(t and tbar)',ytb,wt,wt2,
     & -2.6d0,2.6d0,0.2d0,'lin')
      n=n+1

c--- pt(t,tbar)
      call bookplot(n,tag,'pt(t,tbar)',ptttb,wt,wt2,
     & -0.01d0,399.99d0,10d0,'log')
      n=n+1

c--- y(t,tbar)
      call bookplot(n,tag,'y(t,tbar)',yttb,wt,wt2,
     & -2.6d0,2.6d0,0.2d0,'lin')
      n=n+1

c--- m(t,tbar)
      call bookplot(n,tag,'m(t,tbar)',mttb,wt,wt2,
     & 0d0,1400d0,25d0,'log')
      n=n+1

c--- m(W+ candidate)
      call bookplot(n,tag,'m(W+ candidate)',mwp,wt,wt2,
     & 0d0,200d0,5d0,'log')
      n=n+1

c--- m(W- candidate)
      call bookplot(n,tag,'m(W- candidate)',mwm,wt,wt2,
     & 0d0,200d0,5d0,'log')
      n=n+1
      endif

c--- This section tailored for comparison with CMS PAS TOP-11-013

c--- plots for comparison with arXiv:0907.3090
      if (dilepton) then

c--- pt(l+)
      call bookplot(n,tag,'pt(l+)',pt(4,p),wt,wt2,
     & 0d0,1000d0,50d0,'log')
      n=n+1
c--- pt(l-)
      call bookplot(n,tag,'pt(l-)',pt(7,p),wt,wt2,
     & 0d0,1000d0,50d0,'log')
      n=n+1

c--- eta(l+)
      call bookplot(n,tag,'eta(l+)',yrap(4,p),wt,wt2,
     & -4d0,4d0,0.2d0,'lin')
      n=n+1
c--- eta(l-)
      call bookplot(n,tag,'eta(l-)',yrap(7,p),wt,wt2,
     & -4d0,4d0,0.2d0,'lin')
      n=n+1
      endif

c--- m(l+,b)
      mlb=dsqrt((p(4,4)+p(5,4))**2-(p(4,1)+p(5,1))**2
     &         -(p(4,2)+p(5,2))**2-(p(4,3)+p(5,3))**2)
      call bookplot(n,tag,'m(l+,b)',mlb,wt,wt2,
     & 0d0,200d0,5d0,'log')
      n=n+1
c--- m(l-,bb)
      mlb=dsqrt((p(6,4)+p(7,4))**2-(p(6,1)+p(7,1))**2
     &         -(p(6,2)+p(7,2))**2-(p(6,3)+p(7,3))**2)
      call bookplot(n,tag,'m(l-,bb)',mlb,wt,wt2,
     & 0d0,200d0,5d0,'log')
      n=n+1
c--- plots for comparison with arXiv:0907.3090

c--- plots for computing ttbar asymmetry

      if ((failed) .and. (first .eqv. .false.)) then
c--- make sure to increment by the number of plots in the 'else' section
      n=n+3

      else

c--- inclusive case
      call bookplot(n,tag,'y(t)-y(tbar)',yt-ytb,wt,wt2,
     & -4d0,4d0,4d0,'lin')
      n=n+1

c--- for bins of |y(t)-y(tbar)| of (0,0.5), (0.5,1), (1,1.5)
      call bookplot(n,tag,'y(t)-y(tbar)',yt-ytb,wt,wt2,
     & -1.5d0,1.5d0,0.5d0,'lin')
      n=n+1

c--- for bins of |y(t)-y(tbar)| of (1.5,4.5)
      call bookplot(n,tag,'y(t)-y(tbar)',yt-ytb,wt,wt2,
     & -4.5d0,4.5d0,3d0,'lin')
      n=n+1

      endif

c-- lepton asymmetry
      if (dilepton .eqv. .false.) then
        if (plabel(4) .eq. 'ea') then
        call bookplot(n,tag,'Q(lep)*y(lep)',yrap(4,p),wt,wt2,
     &   -4d0,4d0,4d0,'lin')
        n=n+1
        endif
        if (plabel(7) .eq. 'el') then
        call bookplot(n,tag,'Q(lep)*y(lep)',-yrap(7,p),wt,wt2,
     &   -4d0,4d0,4d0,'lin')
        n=n+1
        endif
      endif

c--- plots for computing ttbar asymmetry

c--- invariant mass of two leading jets
      call bookplot(n,tag,'m(j1,j2)',mj1j2,wt,wt2, 0d0,300d0,5d0,'lin')
      n=n+1

c--- compute energy of particle 4 in rest frame of top
      do j=1,4
      plep(j)=p(4,j)
      ptop(j)=p(3,j)+p(4,j)+p(5,j)
      ptoprest(j)=0d0
      enddo
      ptoprest(4)=mt
      call boostx(plep,ptop,ptoprest,pleprest)
      call bookplot(n,tag,'2*E(4)/mt (top rest frame)',
     & 2d0*pleprest(4)/mt,wt,wt2,0d0,1d0,0.01d0,'lin')
      n=n+1

c--- compute energy of particle 7 in rest frame of top
      do j=1,4
      plep(j)=p(7,j)
      ptop(j)=p(6,j)+p(7,j)+p(8,j)
      ptoprest(j)=0d0
      enddo
      ptoprest(4)=mt
      call boostx(plep,ptop,ptoprest,pleprest)
      call bookplot(n,tag,'2*E(7)/mt (top rest frame)',
     & 2d0*pleprest(4)/mt,wt,wt2,0d0,1d0,0.01d0,'lin')
      n=n+1

c--- single-particle plots
      do j=3,8
        call genplot1(p,j,tag,wt,wt2,n)
      enddo
c--- two-particle plots
      call genplot2(p,3,4,tag,wt,wt2,n)
      call genplot2(p,7,8,tag,wt,wt2,n)
      call genplot2(p,5,6,tag,wt,wt2,n)
c--- three-particle plots
      call genplot3(p,3,4,5,tag,wt,wt2,n)
      call genplot3(p,6,7,8,tag,wt,wt2,n)

c--- additional plots that may be present at NLO
      if (abs(p(9,4)) .gt. 1d-8) then
        call genplot1(p,9,tag,wt,wt2,n)
        call genplot3(p,7,8,9,tag,wt,wt2,n)
      else
        n=n+2
      endif

************************************************************************
*                                                                      *
*     FINAL BOOKKEEPING                                                *
*                                                                      *
************************************************************************

c--- We have over-counted the number of histograms by 1 at this point
      n=n-1

c--- Ensure the built-in maximum number of histograms is not exceeded
      call checkmaxhisto(n)

c--- Set the maximum number of plots, on the first call
      if (first) then
        first=.false.
        nplotmax=n
      endif

      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
subroutine getfilename(pid,onum,vnum,gnum,spseq,anumb,enum,enumos,detectorname,prodtype,fileout)
use precision
implicit none
!import vars
integer :: pid,onum,vnum,gnum,spseq,anumb,enum,enumos
character(8) :: detectorname,prodtype
character(200), dimension(3) :: fileout !name of FITS files
!local vars
character(200) :: basename

!There are three FITS files made.
!1) unconvolved Image
!2) convolved Image
!3) convolved Image with oversampled resolution.

!create basename
write(basename,'(A2,I0.5,I0.3,I0.3,A1,I0.2,I1,I0.2,I0.5,A1,A,A1,A)') 'jw', &
	pid,onum,vnum,'_',gnum,spseq,anumb,enum,'_',detectorname,'_',prodtype
fileout(1) = trim(basename)//'.fits'
fileout(2) = trim(basename)//'_c.fits'

!deal with oversampled case.
write(basename,'(A2,I0.5,I0.3,I0.3,A1,I0.2,I1,I0.2,I0.5,A1,A,A1,A)') 'jw', &
	pid,onum,vnum,'_',gnum,spseq,anumb,enumos,'_',detectorname,'_',prodtype
fileout(3) = trim(basename)//'_ovc.fits'

return
end subroutine getfilename

max_stars_repo_name: AntoineDarveau/jwst-mtl

max_issues_repo_name: AntoineDarveau/jwst-mtl

max_forks_repo_name: AntoineDarveau/jwst-mtl


text:
program hello
real :: x

print *, "Hello world!"
x = sqrt(25.0)
print *, x

end program hello

max_stars_repo_name: jlant/playground

max_issues_repo_name: jlant/playground

max_forks_repo_name: jlant/playground


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details

      SUBROUTINE INFACT (ML, MS, IFOUND, FACT, IIN, N19, N20, FACTOR,
     &   NLPS, IFLINE, ILLIST, LINKL, LINKS, ADDLNK)
C***********************************************************************

C  SUBROUTINE INFACT = ENTERS FACTORS INTO THE DATABASE

C***********************************************************************

      DIMENSION FACTOR (ML), IIN (IFOUND)
      DIMENSION NLPS (MS), IFLINE (MS), ILLIST (MS * 3)
      DIMENSION LINKL (2, ML), LINKS (2, MS)

      LOGICAL ADDLNK

      IF (FACT .LE. 0.)FACT = 1.0
      DO 110 I = 1, IFOUND
         JJ = IIN (I)

C  APPLY INTERVALS TO A SIDE

         IF (JJ .LT. 0) THEN
            JJ = - JJ
            CALL LTSORT (MS, LINKS, JJ, IPNTR, ADDLNK)
            IF ( (JJ .GT. N20) .OR. (IPNTR .LE. 0)) THEN
               WRITE ( * , 10000)JJ
            ELSE
               DO 100 JJJ = IFLINE (IPNTR),
     &            IFLINE (IPNTR) + NLPS (IPNTR) - 1
                  CALL LTSORT (ML, LINKL, ILLIST (JJJ), JPNTR, ADDLNK)
                  IF ( (ILLIST (JJJ) .GT. N19) .OR. (IPNTR .LE. 0))
     &               THEN
                     WRITE ( * , 10010)ILLIST (JJJ)
                  ELSE
                     FACTOR (JPNTR) = FACT
                  ENDIF
  100          CONTINUE
            ENDIF
         ELSE

C  INPUT A FACTOR FOR A LINE

            CALL LTSORT (ML, LINKL, JJ, IPNTR, ADDLNK)
            IF ( (JJ .GT. N19) .OR. (IPNTR .LE. 0)) THEN
               WRITE ( * , 10010)JJ
            ELSE
               FACTOR (IPNTR) = FACT
            ENDIF
         ENDIF
  110 CONTINUE

      RETURN

10000 FORMAT (' SIDE NO:', I5, ' IS NOT IN THE DATABASE',  / ,
     &   ' THUS NO INTERVALS CAN BE ENTERED')
10010 FORMAT (' LINE NO:', I5, ' IS NOT IN THE DATABASE',  / ,
     &   ' THUS NO INTERVAL CAN BE ENTERED')

      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas