Module AddLoopFunctions
Use Control
Use Settings
!Use Mathematics, Only: Li2, CLi2, Delt
!Use LoopFunctions
Use LoopFunctions !, Only: DB0, A0, B0, B1, B00, B11, Cget, Log1minusX, GetRenormalizationScale

! private variables
Real(dp), Private :: mudim2 = 0._dp !, divergence = 0._dp
Real(dp), Private :: xloss = 0.125_dp
! Integer, Private :: divonly = 0 !when set to 1, only UV divergent parts are taken
                                !not for special C routines
Real(dp), Private :: controlmdiff = 1.0e-15_dp !controls the mass difference
                                               !when comparing two masses
Real(dp) :: CCut = 10.e-15_dp

! taken parts from Stefan Liebler

Contains

Logical Function EqualMasses(x,y)
Implicit None
Real(dp), Intent(in) :: x,y
Real(dp) :: epsD
epsD = 1.0E-20_dp
If ((Max(Abs(x),Abs(y)).lt.epsD)) Then
 EqualMasses = .True.
Else
 If ((Abs(x-y)/Max(Abs(x),Abs(y))).le.epsD) Then
  EqualMasses = .True.
!   Write(*,*) x,y
 Else
  EqualMasses = .False.
 End if
End if
End Function EqualMasses

Complex(dp) Function C0_IR(rp1,rp2,rp3,m1,m2,m3)
Implicit None
Real(dp), Intent(in) :: rp1, rp2, rp3, m1, m2, m3
Complex(dp) :: pp1, pp2, pp3
Real(dp) :: epsIR
pp1=cmplx(rp1,0._dp,dp)
pp2=cmplx(rp2,0._dp,dp)
pp3=cmplx(rp3,0._dp,dp)

! Write(*,*) pp1,pp2,pp3

! If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.rp1).and.(m3.eq.rp2)) Then
If ((EqualMasses(m2,Mass_Regulator_PhotonGluon**2)).and.(EqualMasses(m1,rp1)).and.(EqualMasses(m3,rp2))) Then
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2/sqrt(m1*m3))
! C0_IR = zLog((pp1 + pp2 - pp3 + zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3)))/(2.*zSqrt(pp1*pp2)))/(epsIR*zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3)))
C0_IR = getC_IR(pp3,pp1,pp2,Mass_Regulator_PhotonGluon)
Else If ((EqualMasses(m1,Mass_Regulator_PhotonGluon**2)).and.(EqualMasses(m2,rp1)).and.(EqualMasses(m3,rp3))) Then
! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.rp1).and.(m3.eq.rp3)) Then
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2/sqrt(m2*m3))
! C0_IR = zLog((pp1 - pp2 + pp3 + zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3)))/(2.*zSqrt(pp1*pp3)))/(epsIR*zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3)))
C0_IR = getC_IR(pp2,pp1,pp3,Mass_Regulator_PhotonGluon)
! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.rp2).and.(m1.eq.rp3)) Then
Else If ((EqualMasses(m3,Mass_Regulator_PhotonGluon**2)).and.(EqualMasses(m1,rp3)).and.(EqualMasses(m2,rp2))) Then
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2/sqrt(m1*m2))
! C0_IR = zLog((-pp1 + pp2 + pp3 + zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3)))/(2.*zSqrt(pp2*pp3)))/(epsIR*zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3)))
C0_IR = getC_IR(pp1,pp2,pp3,Mass_Regulator_PhotonGluon)
Else
C0_IR = 0._dp
End if

If (Abs(C0_IR).ne.Abs(C0_IR)) Then
Write(*,*) "NaN in C0_IR", rp1,rp2,rp3,m1,m2,m3
C0_IR= 0._dp
End if

End Function C0_IR

Complex(dp) Function getC_IR(pp1,pp2,pp3,reg)
Complex(dp), Intent(in) :: pp1, pp2, pp3
Real(dp), Intent(in) :: reg
Real(dp) :: epsIR
Complex(dp) :: Creg

If ((Abs((pp2/pp1)).lt.1E-8_dp).and.(Abs((pp3/pp1)).lt.1E-8_dp)) Then
  Creg=cmplx(reg,0._dp,dp)
!  getC_IR = ((pp1*(pp2 + pp1) + pp3*(4*pp2 + pp1))*zLog(Creg**2)*zLog(-1._dp/pp1))/pp1**3
 epsIR = zlog(Creg**2/(-pp1))
 getC_IR = epsIR/pp1*zlog(zSqrt(pp2*pp3)/(-pp1))
!  Write(*,*) "A", pp1, pp2, pp3, getC_IR
else if (Abs((pp2/pp1)).lt.1E-8_dp) Then
 Creg=cmplx(reg,0._dp,dp)
 epsIR = zlog(zSqrt(pp2)*(pp3-pp1)/(Creg**2*zSqrt(pp3)))
 getC_IR = 1._dp/(pp1-pp3)*epsIR*zLog((pp3-pp1)/zSqrt(pp2*pp3))
!  Write(*,*) "B", pp1, pp2, pp3, getC_IR
else if (Abs((pp3/pp1)).lt.1E-8_dp) Then
 Creg=cmplx(reg,0._dp,dp)
 epsIR = zlog(zSqrt(pp3)*(pp2-pp1)/(Creg**2*zSqrt(pp2)))
 getC_IR = 1._dp/(pp1-pp2)*epsIR*zLog((pp2-pp1)/zSqrt(pp2*pp3))
!  Write(*,*) "C", pp1, pp2, pp3, getC_IR
else
 epsIR = 1._dp/log(reg**2/sqrt(pp2*pp3))
 getC_IR = zLog((-pp1 + pp2 + pp3 + zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3)))/(2.*zSqrt(pp2*pp3)))/ &
  & (epsIR*zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3)))
!   Write(*,*) "D", pp1, pp2, pp3,getC_IR
end if

If (Abs(getC_IR).ne.Abs(getC_IR)) Then
 getC_IR = ((pp1*(pp2 + pp3) + pp3*(4*pp2 + pp1))*zLog(Creg**2)*zLog(-1._dp/pp1))/pp1**3
End if

! Write(*,*) pp1,pp2,pp3,reg, getC_IR

End Function getC_IR

Complex(dp) Function C1_IR(rp1,rp2,rp3,m1,m2,m3)
Implicit None
Real(dp), Intent(in) :: rp1, rp2, rp3, m1, m2, m3
Complex(dp) :: pp1, pp2, pp3
Real(dp) :: epsIR
pp1=cmplx(rp1,0._dp,dp)
pp2=cmplx(rp2,0._dp,dp)
pp3=cmplx(rp3,0._dp,dp)

! If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.rp1).and.(m3.eq.rp2)) Then
If ((EqualMasses(m2,Mass_Regulator_PhotonGluon**2)).and.(EqualMasses(m1,rp1)).and.(EqualMasses(m3,rp2))) Then
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2/sqrt(m1*m3))
! C1_IR = -(zLog((pp1 + pp2 - pp3 + zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3)))/(2.*zSqrt(pp1*pp2)))/(epsIR*zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3))))
C1_IR = -getC_IR(pp3,pp1,pp2,Mass_Regulator_PhotonGluon)
Else
C1_IR = 0._dp
End if

If (Abs(C1_IR).ne.Abs(C1_IR)) Then
! Write(*,*) "NaN in C1_IR", rp1,rp2,rp3,m1,m2,m3
C1_IR=0._dp
End if

End Function C1_IR

Complex(dp) Function C2_IR(rp1,rp2,rp3,m1,m2,m3)
Implicit None
Real(dp), Intent(in) :: rp1, rp2, rp3, m1, m2, m3
Complex(dp) :: pp1, pp2, pp3
Real(dp) :: epsIR
pp1=cmplx(rp1,0._dp,dp)
pp2=cmplx(rp2,0._dp,dp)
pp3=cmplx(rp3,0._dp,dp)

! If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.rp2).and.(m1.eq.rp3)) Then
If ((EqualMasses(m3,Mass_Regulator_PhotonGluon**2)).and.(EqualMasses(m1,rp3)).and.(EqualMasses(m2,rp2))) Then
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2/sqrt(m1*m2))
! C2_IR = -(zLog((-pp1 + pp2 + pp3 + zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3)))/(2.*zSqrt(pp2*pp3)))/(epsIR*zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3))))
C2_IR = -getC_IR(pp1,pp2,pp3,Mass_Regulator_PhotonGluon)
Else
C2_IR = 0._dp
End if

If (Abs(C2_IR).ne.Abs(C2_IR)) Then
! Write(*,*) "NaN in C2_IR", rp1,rp2,rp3,m1,m2,m3
C2_IR=0._dp
End if

End Function C2_IR

Complex(dp) Function C00_IR(rp1,rp2,rp3,m1,m2,m3)
Implicit None
Real(dp), Intent(in) :: rp1, rp2, rp3, m1, m2, m3
Complex(dp) :: pp1, pp2, pp3
Real(dp) :: epsIR
pp1=cmplx(rp1,0._dp,dp)
pp2=cmplx(rp2,0._dp,dp)
pp3=cmplx(rp3,0._dp,dp)

C00_IR = 0._dp

If (Abs(C00_IR).ne.Abs(C00_IR)) Then
! Write(*,*) "NaN in C00_IR", rp1,rp2,rp3,m1,m2,m3
End if

End Function C00_IR

Complex(dp) Function C11_IR(rp1,rp2,rp3,m1,m2,m3)
Implicit None
Real(dp), Intent(in) :: rp1, rp2, rp3, m1, m2, m3
Complex(dp) :: pp1, pp2, pp3
Real(dp) :: epsIR
pp1=cmplx(rp1,0._dp,dp)
pp2=cmplx(rp2,0._dp,dp)
pp3=cmplx(rp3,0._dp,dp)

! If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.rp1).and.(m3.eq.rp2)) Then
If ((EqualMasses(m2,Mass_Regulator_PhotonGluon**2)).and.(EqualMasses(m1,rp1)).and.(EqualMasses(m3,rp2))) Then
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2/sqrt(m1*m3))
! C11_IR = zLog((pp1 + pp2 - pp3 + zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3)))/(2.*zSqrt(pp1*pp2)))/(epsIR*zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3)))
C11_IR = getC_IR(pp3,pp1,pp2,Mass_Regulator_PhotonGluon)
Else
C11_IR = 0._dp
End if

If (Abs(C11_IR).ne.Abs(C11_IR)) Then
! Write(*,*) "NaN in C11_IR", rp1,rp2,rp3,m1,m2,m3
C11_IR=0._dp
End if

End Function C11_IR

Complex(dp) Function C12_IR(rp1,rp2,rp3,m1,m2,m3)
Implicit None
Real(dp), Intent(in) :: rp1, rp2, rp3, m1, m2, m3
Complex(dp) :: pp1, pp2, pp3
Real(dp) :: epsIR
pp1=cmplx(rp1,0._dp,dp)
pp2=cmplx(rp2,0._dp,dp)
pp3=cmplx(rp3,0._dp,dp)

C12_IR = 0._dp

If (Abs(C12_IR).ne.Abs(C12_IR)) Then
! Write(*,*) "NaN in C12_IR", rp1,rp2,rp3,m1,m2,m3
C12_IR=0._dp
End if

End Function C12_IR

Complex(dp) Function C22_IR(rp1,rp2,rp3,m1,m2,m3)
Implicit None
Real(dp), Intent(in) :: rp1, rp2, rp3, m1, m2, m3
Complex(dp) :: pp1, pp2, pp3
Real(dp) :: epsIR
pp1=cmplx(rp1,0._dp,dp)
pp2=cmplx(rp2,0._dp,dp)
pp3=cmplx(rp3,0._dp,dp)

! If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.rp2).and.(m1.eq.rp3)) Then
If ((EqualMasses(m3,Mass_Regulator_PhotonGluon**2)).and.(EqualMasses(m1,rp3)).and.(EqualMasses(m2,rp2))) Then
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2/sqrt(m1*m2))
! C22_IR = zLog((-pp1 + pp2 + pp3 + zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3)))/(2.*zSqrt(pp2*pp3)))/(epsIR*zSqrt(pp1**2 + (pp2 - pp3)**2 - 2*pp1*(pp2 + pp3)))
C22_IR = getC_IR(pp1,pp2,pp3,Mass_Regulator_PhotonGluon)
Else
C22_IR = 0._dp
End if

If (Abs(C22_IR).ne.Abs(C22_IR)) Then
! Write(*,*) "NaN in C22_IR", rp1,rp2,rp3,m1,m2,m3
C22_IR=0._dp
End if

End Function C22_IR

!
!
! Complex(dp) Function C0_IR(p1,p2,p3,m1,m2,m3)
! Implicit None
! Real(dp), Intent(in) :: p1, p2, p3, m1, m2, m3
! Real(dp) :: epsIR
! ! If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m3.eq.p3)) Then
! ! C0_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m3.eq.p3)) Then
! ! C0_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m2.eq.p2)) Then
! ! C0_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p2).and.(m2.eq.p3)) Then
! ! C0_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m2.eq.p3)) Then
! ! C0_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m3.eq.p2)) Then
! ! C0_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p2).and.(m3.eq.p3)) Then
! ! C0_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m3.eq.p3)) Then
! ! C0_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m1.eq.p2)) Then
! ! C0_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p2).and.(m2.eq.p3)) Then
! ! C0_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m2.eq.p3)) Then
! ! C0_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m1.eq.p2)) Then
! ! C0_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p2).and.(m1.eq.p3)) Then
! ! C0_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m1.eq.p3)) Then
! ! C0_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m3.eq.p2)) Then
! ! C0_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m1.eq.p3)) Then
! If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m1.eq.p3)) Then
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2/sqrt(m1*m2))
!  C0_IR = zLog((-p1 + p2 + p3 + zSqrt(Cmplx(p1**2 + (p2 - p3)**2 - 2._dp*p1*(p2 + p3),0._dp,dp)))/(2._dp*zSqrt(Cmplx(p2*p3,0._dp,dp))))/ &
!  & (epsIR*zSqrt(Cmplx(p1**2 + (p2 - p3)**2 - 2._dp*p1*(p2 + p3),0._dp,dp)))
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m1.eq.p3)) Then
! ! C0_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m2.eq.p2)) Then
! ! C0_IR = 0
! Else
! C0_IR = 0._dp
! End if
!
!
! ! if ((Abs(C0_IR)).ne.(Abs(C0_IR)))  Then
! ! Write(*,*) "NaN in C0"
! ! Write(*,*) "input:", p1,p2,p3,m1,m2,m3,epsIR
! ! End if
!
! ! Write(*,*) "C0", C0_IR
!
! End Function C0_IR

!
!
! Complex(dp) Function C1_IR(p1,p2,p3,m1,m2,m3)
! Implicit None
! Real(dp), Intent(in) :: p1, p2, p3, m1, m2, m3
! Real(dp) :: epsIR
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2)
! ! If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m3.eq.p3)) Then
! ! C1_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m3.eq.p3)) Then
! ! C1_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m2.eq.p2)) Then
! ! C1_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p2).and.(m2.eq.p3)) Then
! ! C1_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m2.eq.p3)) Then
! ! C1_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m3.eq.p2)) Then
! ! C1_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p2).and.(m3.eq.p3)) Then
! ! C1_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m3.eq.p3)) Then
! ! C1_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m1.eq.p2)) Then
! ! C1_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p2).and.(m2.eq.p3)) Then
! ! C1_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m2.eq.p3)) Then
! ! C1_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m1.eq.p2)) Then
! ! C1_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p2).and.(m1.eq.p3)) Then
! ! C1_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m1.eq.p3)) Then
! ! C1_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m3.eq.p2)) Then
! ! C1_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m1.eq.p3)) Then
! ! C1_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m1.eq.p3)) Then
! ! C1_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m2.eq.p2)) Then
! ! C1_IR = 0
! ! Else
! C1_IR = 0._dp
! ! End if
! End Function C1_IR
!
!
! Complex(dp) Function C2_IR(p1,p2,p3,m1,m2,m3)
! Implicit None
! Real(dp), Intent(in) :: p1, p2, p3, m1, m2, m3
! Real(dp) :: epsIR
! Complex(dp) :: arg
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2)
! ! If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m3.eq.p3)) Then
! ! C2_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m3.eq.p3)) Then
! ! C2_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m2.eq.p2)) Then
! ! C2_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p2).and.(m2.eq.p3)) Then
! ! C2_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m2.eq.p3)) Then
! ! C2_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m3.eq.p2)) Then
! ! C2_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p2).and.(m3.eq.p3)) Then
! ! C2_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m3.eq.p3)) Then
! ! C2_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m1.eq.p2)) Then
! ! C2_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p2).and.(m2.eq.p3)) Then
! ! C2_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m2.eq.p3)) Then
! ! C2_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m1.eq.p2)) Then
! ! C2_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p2).and.(m1.eq.p3)) Then
! ! C2_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m1.eq.p3)) Then
! ! C2_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m3.eq.p2)) Then
! ! C2_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m1.eq.p3)) Then
! If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m1.eq.p3)) Then
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2/sqrt(m1*m2))
!  arg = (-p1 + p2 + p3 + zSqrt(CMPLX(p1**2 + (p2 - p3)**2 - 2._dp*p1*(p2 + p3),0._dp,dp)))/(2._dp*zSqrt(CMPLX(p2*p3,0._dp,dp)))
!  C2_IR = -(zLog(arg)/(epsIR*zSqrt(CMPLX(p1**2 + (p2 - p3)**2 - 2._dp*p1*(p2 + p3),0._dp,dp))))
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m1.eq.p3)) Then
! ! C2_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m2.eq.p2)) Then
! ! C2_IR = 0
! Else
! C2_IR = 0._dp
! End if
!
! ! if ((Abs(C2_IR)).ne.(Abs(C2_IR)))  Then
! ! Write(*,*) "NaN in C2"
! ! Write(*,*) "input:", p1,p2,p3,m1,m2,m3,epsIR
! ! Write(*,*) "arg, Log(arg)", arg, zLog(arg)
! ! Write(*,*) "part I, II", (-p1 + p2 + p3 + zSqrt(CMPLX(p1**2 + (p2 - p3)**2 - 2._dp*p1*(p2 + p3),0._dp,dp))), (2._dp*zSqrt(CMPLX(p2*p3,0._dp,dp)))
! ! End if
!
! End Function C2_IR
!
!
! Complex(dp) Function C00_IR(p1,p2,p3,m1,m2,m3)
! Implicit None
! Real(dp), Intent(in) :: p1, p2, p3, m1, m2, m3
! Real(dp) :: epsIR
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2)
! ! If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m3.eq.p3)) Then
! ! C00_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m3.eq.p3)) Then
! ! C00_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m2.eq.p2)) Then
! ! C00_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p2).and.(m2.eq.p3)) Then
! ! C00_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m2.eq.p3)) Then
! ! C00_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m3.eq.p2)) Then
! ! C00_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p2).and.(m3.eq.p3)) Then
! ! C00_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m3.eq.p3)) Then
! ! C00_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m1.eq.p2)) Then
! ! C00_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p2).and.(m2.eq.p3)) Then
! ! C00_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m2.eq.p3)) Then
! ! C00_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m1.eq.p2)) Then
! ! C00_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p2).and.(m1.eq.p3)) Then
! ! C00_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m1.eq.p3)) Then
! ! C00_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m3.eq.p2)) Then
! ! C00_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m1.eq.p3)) Then
! ! C00_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m1.eq.p3)) Then
! ! C00_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m2.eq.p2)) Then
! ! C00_IR = 0
! ! Else
! C00_IR = 0._dp
! ! End if
! End Function C00_IR
!
!
! Complex(dp) Function C11_IR(p1,p2,p3,m1,m2,m3)
! Implicit None
! Real(dp), Intent(in) :: p1, p2, p3, m1, m2, m3
! Real(dp) :: epsIR
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2)
! ! If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m3.eq.p3)) Then
! ! C11_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m3.eq.p3)) Then
! ! C11_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m2.eq.p2)) Then
! ! C11_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p2).and.(m2.eq.p3)) Then
! ! C11_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m2.eq.p3)) Then
! ! C11_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m3.eq.p2)) Then
! ! C11_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p2).and.(m3.eq.p3)) Then
! ! C11_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m3.eq.p3)) Then
! ! C11_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m1.eq.p2)) Then
! ! C11_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p2).and.(m2.eq.p3)) Then
! ! C11_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m2.eq.p3)) Then
! ! C11_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m1.eq.p2)) Then
! ! C11_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p2).and.(m1.eq.p3)) Then
! ! C11_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m1.eq.p3)) Then
! ! C11_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m3.eq.p2)) Then
! ! C11_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m1.eq.p3)) Then
! ! C11_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m1.eq.p3)) Then
! ! C11_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m2.eq.p2)) Then
! ! C11_IR = 0
! ! Else
! C11_IR = 0._dp
! ! End if
! End Function C11_IR
!
!
! Complex(dp) Function C12_IR(p1,p2,p3,m1,m2,m3)
! Implicit None
! Real(dp), Intent(in) :: p1, p2, p3, m1, m2, m3
! Real(dp) :: epsIR
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2)
! ! If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m3.eq.p3)) Then
! ! C12_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m3.eq.p3)) Then
! ! C12_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m2.eq.p2)) Then
! ! C12_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p2).and.(m2.eq.p3)) Then
! ! C12_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m2.eq.p3)) Then
! ! C12_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m3.eq.p2)) Then
! ! C12_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p2).and.(m3.eq.p3)) Then
! ! C12_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m3.eq.p3)) Then
! ! C12_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m1.eq.p2)) Then
! ! C12_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p2).and.(m2.eq.p3)) Then
! ! C12_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m2.eq.p3)) Then
! ! C12_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m1.eq.p2)) Then
! ! C12_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p2).and.(m1.eq.p3)) Then
! ! C12_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m1.eq.p3)) Then
! ! C12_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m3.eq.p2)) Then
! ! C12_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m1.eq.p3)) Then
! ! C12_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m1.eq.p3)) Then
! ! C12_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m2.eq.p2)) Then
! ! C12_IR = 0
! ! Else
! C12_IR = 0._dp
! ! End if
! End Function C12_IR
!
!
! Complex(dp) Function C22_IR(p1,p2,p3,m1,m2,m3)
! Implicit None
! Real(dp), Intent(in) :: p1, p2, p3, m1, m2, m3
! Real(dp) :: epsIR
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2)
! ! If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m3.eq.p3)) Then
! ! C22_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m3.eq.p3)) Then
! ! C22_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m2.eq.p2)) Then
! ! C22_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p2).and.(m2.eq.p3)) Then
! ! C22_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m2.eq.p3)) Then
! ! C22_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p1).and.(m3.eq.p2)) Then
! ! C22_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p2).and.(m3.eq.p3)) Then
! ! C22_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m3.eq.p3)) Then
! ! C22_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m1.eq.p2)) Then
! ! C22_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m1.eq.p2).and.(m2.eq.p3)) Then
! ! C22_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m2.eq.p3)) Then
! ! C22_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m1.eq.p2)) Then
! ! C22_IR = 0
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p2).and.(m1.eq.p3)) Then
! ! C22_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m1.eq.p3)) Then
! ! C22_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p1).and.(m3.eq.p2)) Then
! ! C22_IR = 0
! ! Else If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m1.eq.p3)) Then
! If ((m3.eq.Mass_Regulator_PhotonGluon**2).and.(m2.eq.p2).and.(m1.eq.p3)) Then
! epsIR = 1._dp/log(Mass_Regulator_PhotonGluon**2/sqrt(m1*m2))
!  C22_IR = zLog(((-p1 + p2 + p3 + zSqrt(Cmplx(p1**2 + (p2 - p3)**2 - 2*p1*(p2 + p3),0._dp,dp)))/ &
!  & (2._dp*zSqrt(Cmplx(p2*p3,0._dp,dp)))))/(epsIR*zSqrt(Cmplx(p1**2 + (p2 - p3)**2 - 2*p1*(p2 + p3),0._dp,dp)))
! ! Else If ((m2.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m1.eq.p3)) Then
! ! C22_IR = 0
! ! Else If ((m1.eq.Mass_Regulator_PhotonGluon**2).and.(m3.eq.p1).and.(m2.eq.p2)) Then
! ! C22_IR = 0
! Else
! C22_IR = 0._dp
! End if
!
!
! ! if ((Abs(C22_IR)).ne.(Abs(C22_IR))) Then
! ! Write(*,*) "NaN in C22"
! ! Write(*,*) "input:", p1,p2,p3,m1,m2,m3,epsIR
! ! End if
!
! End Function C22_IR

Complex(dp) Function SA_B0(xp, xm1, xm2)
 Implicit None
  Real(dp), Intent(in) :: xp, xm1, xm2

! FS: changed 15/07/2017
    If ((xm1.eq.0._dp).and.(xm2.eq.0._dp).and.(xp.eq.0._dp)) Then
      SA_B0=0._dp ! to prevent numerical problems.
      Return
    End if
!

  If (divonly.Eq.1) Then
   If (IRdivonly) Then
    SA_B0 = 0._dp
   Else
    SA_B0 = divergence
   End if
   Return
  End If

  If (IRdivonly) Then
   SA_B0 = 0._dp
   Return
  End if

  SA_B0 = B0(xp,xm1,xm2)

End Function SA_B0

Complex(dp) Function SA_A0(m)
 Implicit None
  Real(dp), Intent(in) :: m

  If (divonly.Eq.1) Then
   If (IRdivonly) Then
    SA_A0 = 0._dp
   Else
    SA_A0 = m*divergence
   End if
   Return
  End If

  If (IRdivonly) Then
   SA_A0 = 0._dp
   Return
  End if

  SA_A0 = A0(m)

End Function SA_A0

Complex(dp) Function SA_B1(xp, xm1, xm2)
 Implicit None
  Real(dp), Intent(in) :: xp, xm1, xm2

  If ((xm1.eq.0._dp).and.(xm2.eq.0._dp)) Then
    SA_B1=0._dp
    Return
  End if

  If (divonly.Eq.1) Then
   If (IRdivonly) Then
    SA_B1 = 0._dp
   Else
    SA_B1 = - 0.5_dp * divergence
   End if
   Return
  End If

  If (IRdivonly) Then
   SA_B1 = 0._dp
   Return
  End if

  SA_B1 = B1(xp,xm1,xm2)

End Function SA_B1

Complex(dp) Function SA_B00(xp, xm1, xm2)
 Implicit None
  Real(dp), Intent(in) :: xp, xm1, xm2

  If (divonly.Eq.1) Then
   If (IRdivonly) Then
    SA_B00 = 0._dp
   Else
    SA_B00 = divergence * (3._dp * xm1 + 3._dp * xm2 - xp)/12._dp
   End if
   Return
  End If

  SA_B00 = B00(xp,xm1,xm2)

End Function SA_B00

Complex(dp) Function SA_Cget(name, p1, p2, p1p2, m1, m2, m3)
 Implicit None
  Real(dp), Intent(in) :: p1, p2, p1p2, m1, m2, m3
  Character(len=4), Intent(in) :: name

  Complex(dp) :: C_0, C1, C2, C00, C11, C12, C22, C001, C002, C111, C112  &
    & , C122, C222

  Iname = Iname + 1
  NameOfUnit(Iname) = "Cget"

 If (divonly.Eq.1) Then
  if (IRdivonly) Then
    SA_Cget = 0._dp
  Else
   If ((name.Eq."C0").Or.(name.Eq."C1").Or.(name.Eq."C2")) Then
    SA_Cget = ZeroC
   Else If (name.Eq."C00") Then
    SA_Cget = 0.25_dp * divergence
   Else If ((name.Eq."C11").Or.(name.Eq."C12").Or.(name.Eq."C22")) Then
    SA_Cget = ZeroC
   Else If ((name.Eq."C001").Or.(name.Eq."C002")) Then
    SA_Cget = - divergence / 12._dp
   Else If ((name.Eq."C111").Or.(name.Eq."C112").Or.(name.Eq."C122").Or.(name.Eq."C222")) Then
    SA_Cget = ZeroC
   Else
    Write (ErrCan,*) "Problem in function Cget, function "//name
    Write (ErrCan,*) "is not defined"
    If (Errorlevel.Ge.0) Call TerminateProgram
   End If
  End if
!   Return

 Else if (IRdivonly) Then

  If (name.Eq."C0") Then
   SA_Cget = C0_IR(p1,p2,p1p2, m1, m2, m3)
  Else If (name.Eq."C1") Then
   SA_Cget = C1_IR(p1,p2,p1p2, m1, m2, m3)
  Else If (name.Eq."C2") Then
   SA_Cget = C2_IR(p1,p2,p1p2, m1, m2, m3)
  Else If (name.Eq."C11") Then
   SA_Cget = C11_IR(p1,p2,p1p2, m1, m2, m3)
  Else If (name.Eq."C12") Then
   SA_Cget = C12_IR(p1,p2,p1p2, m1, m2, m3)
  Else If (name.Eq."C22") Then
   SA_Cget = C22_IR(p1,p2,p1p2, m1, m2, m3)
  Else If (name.Eq."C00") Then
   SA_Cget = C00_IR(p1,p2,p1p2, m1, m2, m3)
  Else
   Write(*,*) "unknown IR divergence: ", name
  End If

 Else

  SA_Cget = Cget(name, p1, p2, p1p2, m1, m2, m3)

 End if
Iname = Iname - 1

End Function SA_Cget

Complex(dp) Function FloopRXi(p2,m12,m22)
   Implicit None
   Real(dp),Intent(in)::p2,m12,m22

    If (RXi.eq.1._dp) Then
!        FloopRXi=Floop(p2,m12,m22)
!             Write(*,*) "frxi A",FloopRXi
     FloopRXi = SA_A0(m12) - 2._dp * SA_A0(m22) &
       & -  (2._dp * p2 + 2._dp * m12 - m22) * SA_B0(p2,m12,m22)
!             Write(*,*) "frxi B",FloopRXi
    Else
       If ((m12.gt.1.0E-10_dp).and.(m22.gt.1.0E-10_dp)) Then
          FloopRXi=A0(m12)-A0(m22)+((m12-p2)*A0(m22))/m22-&
               & ((m12-p2+m22*RXi)*A0(m22*RXi))/m22+(-m12+m22+p2)*b0(p2,m12,m22)-&
               & (m12-(m12-p2)**2/m22+3._dp*p2)*b0(p2,m12,m22)-((m12-p2)**2*b0(p2,m12,m22*RXi))/m22
       Else
          If (m12.gt.1.0E-10_dp) Then
              FloopRXi=A0(m12)-2._dp*(m12+p2)*B0(p2,0._dp,m12)
          Else If (m22.gt.1.0E-10_dp) Then
              FloopRXi=(-((m22+p2)*A0(m22))+(p2-m22*RXi)*A0(m22*RXi)+(m22-p2)**2*B0(p2,0._dp,m22)-p2**2*B0(p2,0._dp,m22*RXi))/m22
          Else
              FloopRXi=0._dp
          End if
       End if
    End if
     If (Real(FloopRXi,dp).ne.Real(FloopRXi,dp)) Write(*,*) "NaN in FloopRXi", p2,m12,m22

End Function FloopRXi

Complex(dp) Function SVVloop(p2,m12,m22)
    Implicit None

    Real(dp),Intent(in)::p2,m12,m22

    If ((m12.eq.0._dp).and.(m22.eq.0._dp)) Then
      SVVloop = 0._dp
      Return
    End if

    If (RXi.eq.1._dp) Then

        SVVloop = 4._dp*Real(SA_B0(p2,m12,m22)-0.5_dp*rMS,dp)
    Else If (RXi.lt.0.1_dp) Then
	If (m12.gt.1.0E-10_dp) Then
	    If (p2.gt.0.1_dp) Then
		SVVloop = -2._dp*(rMS-A0(m12)/(8._dp*m12)+(RXi*A0(m12))/(8._dp*m12)-A0(m22)/(8._dp*m22)+(RXi*A0(m22))/(8._dp*m22)+A0(m12*RXi)/(8._dp*m12)&
               &-(RXi*A0(m12*RXi))/(8._dp*m12)+A0(m22*RXi)/(8._dp*m22)-&
               &(RXi*A0(m22*RXi))/(8._dp*m22)-(5._dp*B0(p2,m12,m22))/4._dp-(m12*B0(p2,m12,m22))/(8._dp*m22)&
               &-(m22*B0(p2,m12,m22))/(8._dp*m12)+(p2*B0(p2,m12,m22))/(4._dp*m12)+(p2*B0(p2,m12,m22))/(4._dp*m22)-&
               &(p2**2*B0(p2,m12,m22))/(8._dp*m12*m22)+(m12*B0(p2,m12,m22*RXi))/(8._dp*m22)-(p2*B0(p2,m12,m22*RXi))/(4._dp*m22)&
               &+(p2**2*B0(p2,m12,m22*RXi))/(8._dp*m12*m22)-(RXi*B0(p2,m12,m22*RXi))/4._dp-&
               &(p2*RXi*B0(p2,m12,m22*RXi))/(4._dp*m12)+(m22*RXi**2*B0(p2,m12,m22*RXi))/(8._dp*m12)&
               &+(m22*B0(p2,m22,m12*RXi))/(8._dp*m12)-(p2*B0(p2,m22,m12*RXi))/(4._dp*m12)+&
               &(p2**2*B0(p2,m22,m12*RXi))/(8._dp*m12*m22)-(RXi*B0(p2,m22,m12*RXi))/4._dp-(p2*RXi*B0(p2,m22,m12*RXi))/(4._dp*m22)&
               &+(m12*RXi**2*B0(p2,m22,m12*RXi))/(8._dp*m22)-&
               &(p2**2*B0(p2,m12*RXi,m22*RXi))/(8._dp*m12*m22)+(p2*RXi*B0(p2,m12*RXi,m22*RXi))/(4._dp*m12)&
               &+(p2*RXi*B0(p2,m12*RXi,m22*RXi))/(4._dp*m22)-(RXi**2*B0(p2,m12*RXi,m22*RXi))/4._dp-&
               &(m12*RXi**2*B0(p2,m12*RXi,m22*RXi))/(8._dp*m22)-(m22*RXi**2*B0(p2,m12*RXi,m22*RXi))/(8._dp*m12))
!                Write(*,*) SVVloop
	    Else
	     If (m12.eq.m22) Then
	        SVVloop = 0._dp
	     Else
		SVVloop = (2._dp*rMS*(-m12+m22)+3._dp*A0(m12)-3._dp*A0(m22))/(m12-m22)
             End if
	    End if
        Else
	     If (p2.gt.0.1) Then
		SVVloop = (-8._dp*rMS*m22+A0(m22)+2._dp*p2*B0(p2,0._dp,0._dp)+2._dp*(5._dp*m22-p2)*B0(p2,0._dp,m22))/(4._dp*m22)
	    Else
		SVVloop = -2._dp*rMS + 3._dp*B0(0._dp,0._dp,m22)
	    End if
	End if
    Else

       If ((m12.gt.1.0E-10_dp).and.(m22.gt.1.0E-10_dp)) Then
	  If (p2.lt.0.1) Then
	      SVVloop = (-2._dp*rMS*m12+2._dp*rMS*m22+3._dp*A0(m12)-3._dp*A0(m22)+RXi*A0(m12*RXi)-RXi*A0(m22*RXi))/(m12-m22)
	  Else

	   If (m12-m22.lt.1_dp) Then
	     SVVloop=  (p2**2*(B0(p2,m12,m12) - 2*B0(p2,m12,m12*RXi) + B0(p2,m12*RXi,m12*RXi)) + &
     &  4*m12*p2*(-B0(p2,m12,m12) + (1 + RXi)*B0(p2,m12,m12*RXi) - RXi*B0(p2,m12*RXi,m12*RXi)) + &
     &  2*m12*((1 - RXi)*A0(m12) + (-1 + RXi)*A0(m12*RXi) - m12*(4*rMS - 6*B0(p2,m12,m12) + (-1 + RXi)**2*B0(p2,m12,m12*RXi) - 2*RXi**2*B0(p2,m12*RXi,m12*RXi))))/(4*m12**2)
	   Else
	      SVVloop = -2._dp*(rMS-A0(m12)/(8._dp*m12)+(RXi*A0(m12))/(8._dp*m12)-A0(m22)/(8._dp*m22)+(RXi*A0(m22))/(8._dp*m22)+A0(m12*RXi)/(8._dp*m12)&
               &-(RXi*A0(m12*RXi))/(8._dp*m12)+A0(m22*RXi)/(8._dp*m22)-&
               &(RXi*A0(m22*RXi))/(8._dp*m22)-(5._dp*B0(p2,m12,m22))/4._dp-(m12*B0(p2,m12,m22))/(8._dp*m22)&
               &-(m22*B0(p2,m12,m22))/(8._dp*m12)+(p2*B0(p2,m12,m22))/(4._dp*m12)+(p2*B0(p2,m12,m22))/(4._dp*m22)-&
               &(p2**2*B0(p2,m12,m22))/(8._dp*m12*m22)+(m12*B0(p2,m12,m22*RXi))/(8._dp*m22)-(p2*B0(p2,m12,m22*RXi))/(4._dp*m22)&
               &+(p2**2*B0(p2,m12,m22*RXi))/(8._dp*m12*m22)-(RXi*B0(p2,m12,m22*RXi))/4._dp-&
               &(p2*RXi*B0(p2,m12,m22*RXi))/(4._dp*m12)+(m22*RXi**2*B0(p2,m12,m22*RXi))/(8._dp*m12)&
               &+(m22*B0(p2,m22,m12*RXi))/(8._dp*m12)-(p2*B0(p2,m22,m12*RXi))/(4._dp*m12)+&
               &(p2**2*B0(p2,m22,m12*RXi))/(8._dp*m12*m22)-(RXi*B0(p2,m22,m12*RXi))/4._dp-(p2*RXi*B0(p2,m22,m12*RXi))/(4._dp*m22)&
               &+(m12*RXi**2*B0(p2,m22,m12*RXi))/(8._dp*m22)-&
               &(p2**2*B0(p2,m12*RXi,m22*RXi))/(8._dp*m12*m22)+(p2*RXi*B0(p2,m12*RXi,m22*RXi))/(4._dp*m12)&
               &+(p2*RXi*B0(p2,m12*RXi,m22*RXi))/(4._dp*m22)-(RXi**2*B0(p2,m12*RXi,m22*RXi))/4._dp-&
               &(m12*RXi**2*B0(p2,m12*RXi,m22*RXi))/(8._dp*m22)-(m22*RXi**2*B0(p2,m12*RXi,m22*RXi))/(8._dp*m12))
!                Write(*,*) SVVloop
          End if
          End if
       Else If ((m12.gt.1.0E-10_dp).and.(m22.lt.1.0E-10_dp)) Then
	  If (p2.lt.0.1) Then
	      SVVloop = (-8._dp*rMS*m22+(11._dp+RXi)*A0(m22)+(1._dp+3._dp*RXi)*A0(m22*RXi))/(4._dp*m22)
	  Else
	      SVVloop =  -2._dp*rMS+A0(m12)/(4._dp*m12)-(RXi*A0(m12))/(4._dp*m12)-A0(m12*RXi)/(4._dp*m12)+(RXi*A0(m12*RXi))/(4._dp*m12)&
              &+(5._dp*B0(p2,0._dp,m12))/2._dp-(p2*B0(p2,0._dp,m12))/(2._dp*m12)+(RXi*B0(p2,0._dp,m12))/2._dp+(p2*RXi*B0(p2,0._dp,m12))/(2._dp*m12)&
              &+(p2*B0(p2,0._dp,m12*RXi))/(2._dp*m12)+(RXi*B0(p2,0._dp,m12*RXi))/2._dp-(p2*RXi*B0(p2,0._dp,m12*RXi))/(2._dp*m12) + (RXi**2*B0(p2,0._dp,m12*RXi))/2._dp
          End if

       Else If ((m22.gt.1.0E-10_dp).and.(m12.lt.1.0E-10_dp)) Then
	  If (p2.lt.0.1) Then
	      SVVloop = -2._dp*rMS + (1.0_dp*RXi)*B0(0._dp,0._dp,m22*RXi) + (3._dp)*B0(0._dp,0._dp,m22) !+ A0(m22)/m22
	  Else
	       SVVloop =  -2._dp*rMS+A0(m22)/(4._dp*m22)-(RXi*A0(m22))/(4._dp*m22)-A0(m22*RXi)/(4._dp*m22)+(RXi*A0(m22*RXi))/(4._dp*m22)&
              &+(5._dp*B0(p2,0._dp,m22))/2._dp-(p2*B0(p2,0._dp,m22))/(2._dp*m22)+(RXi*B0(p2,0._dp,m22))/2._dp+(p2*RXi*B0(p2,0._dp,m22))/(2._dp*m22)&
              &+(p2*B0(p2,0._dp,m22*RXi))/(2._dp*m22)+(RXi*B0(p2,0._dp,m22*RXi))/2._dp-(p2*RXi*B0(p2,0._dp,m22*RXi))/(2._dp*m22) + (RXi**2*B0(p2,0._dp,m22*RXi))/2._dp
          End if

       End if
    End if
         If (Real(SVVloop,dp).ne.Real(SVVloop,dp)) Write(*,*) "NaN in SVVloop", p2,m12,m22
  End Function SVVloop

Complex(dp) Function VSSloop(p2,m12,m22)
Implicit None
Real(dp),Intent(in)::p2,m12,m22

VSSloop=  -4._dp*SA_B00(p2,m12,m22)

End Function VSSloop

Complex(dp) Function VVSloop(p2,m12,m22)
Implicit None
Real(dp),Intent(in)::p2,m12,m22

If (RXi.eq.1._dp) Then
  VVSloop=  SA_B0(p2,m12,m22)
Else
   VVSloop = B0(p2,m12,m22) + (-B00(p2,m12,m22) + B00(p2,m22,rxi*m12))/m12 !! FeynArts
End if

End Function VVSloop

Complex(dp) Function VVVloop(p2,m12,m22)
Implicit None
Real(dp),Intent(in)::p2,m12,m22
Real(dp) :: b001

If (RXi.eq.1._dp) Then
  VVVloop=   -2._dp*rMS*(m12 + m22 - p2/3._dp) + 10._dp*SA_B00(p2,m12,m22) &
     & + (m12+m22+4._dp*p2)*SA_B0(p2,m12,m22)+  SA_A0(m12) + SA_A0(m22)
Else
  If ((p2.gt.0).and.(m12.gt.0_dp).and.(m22.gt.0._dp)) Then
      VVVloop= ((m12*m22*p2*(-3*m12 - 3*m22 - 40*p2 - 12*(m12 + m22 - 4*p2)*RXi - 9*(m12 + m22)*RXi**2) +  &
     &      6*m22*(-m22**2 - 9*m22*p2 + m12*(-9*m22 + p2*(13 - 2*RXi)) + m12**2*(9 + RXi) + p2**2*(9 + RXi))*A0(m12) -  &
     &      6*m12*(m12**2 + 9*m12*(m22 + p2) - m22**2*(9 + RXi) - p2**2*(9 + RXi) + m22*p2*(-13 + 2*RXi))*A0(m22) +  &
     &      6*(m22*(m22**2 + m22*(9*p2 - m12*RXi) - p2*(-11*m12*RXi + p2*(9 + RXi)))*A0(m12*RXi) + m12*(m12**2 + 9*m12*p2 - 9*p2**2 - (m12*m22 - 11*m22*p2 + p2**2)*RXi)*A0(m22*RXi) -  &
     &         m12**4*B0(p2,m12,m22) - 8*m12**3*m22*B0(p2,m12,m22) + 18*m12**2*m22**2*B0(p2,m12,m22) - 8*m12*m22**3*B0(p2,m12,m22) - m22**4*B0(p2,m12,m22) - 8*m12**3*p2*B0(p2,m12,m22) + &
     &         32*m12**2*m22*p2*B0(p2,m12,m22) + 32*m12*m22**2*p2*B0(p2,m12,m22) - 8*m22**3*p2*B0(p2,m12,m22) + 18*m12**2*p2**2*B0(p2,m12,m22) + 32*m12*m22*p2**2*B0(p2,m12,m22) +  &
     &         18*m22**2*p2**2*B0(p2,m12,m22) - 8*m12*p2**3*B0(p2,m12,m22) - 8*m22*p2**3*B0(p2,m12,m22) - p2**4*B0(p2,m12,m22) + m12**4*B0(p2,m12,m22*RXi) + 8*m12**3*p2*B0(p2,m12,m22*RXi) -  &
     &         18*m12**2*p2**2*B0(p2,m12,m22*RXi) + 8*m12*p2**3*B0(p2,m12,m22*RXi) + p2**4*B0(p2,m12,m22*RXi) - 2*m12**3*m22*RXi*B0(p2,m12,m22*RXi) + 2*m12**2*m22*p2*RXi*B0(p2,m12,m22*RXi) +  &
     &         2*m12*m22*p2**2*RXi*B0(p2,m12,m22*RXi) - 2*m22*p2**3*RXi*B0(p2,m12,m22*RXi) + m12**2*m22**2*RXi**2*B0(p2,m12,m22*RXi) - 2*m12*m22**2*p2*RXi**2*B0(p2,m12,m22*RXi) +  &
     &         m22**2*p2**2*RXi**2*B0(p2,m12,m22*RXi) + m22**4*B0(p2,m22,m12*RXi) + 8*m22**3*p2*B0(p2,m22,m12*RXi) - 18*m22**2*p2**2*B0(p2,m22,m12*RXi) + 8*m22*p2**3*B0(p2,m22,m12*RXi) +  &
     &         p2**4*B0(p2,m22,m12*RXi) - 2*m12*m22**3*RXi*B0(p2,m22,m12*RXi) + 2*m12*m22**2*p2*RXi*B0(p2,m22,m12*RXi) + 2*m12*m22*p2**2*RXi*B0(p2,m22,m12*RXi) - 2*m12*p2**3*RXi*B0(p2,m22,m12*RXi) +  &
     &         m12**2*m22**2*RXi**2*B0(p2,m22,m12*RXi) - 2*m12**2*m22*p2*RXi**2*B0(p2,m22,m12*RXi) + m12**2*p2**2*RXi**2*B0(p2,m22,m12*RXi) -  &
     &         p2**2*(p2**2 - 2*(m12 + m22)*p2*RXi + (m12 - m22)**2*RXi**2)*B0(p2,m12*RXi,m22*RXi))))/(72.*m12*m22*p2) -2._dp*(rMS-1._dp)*(m12 + m22 - p2/3._dp)
  Else ! to be added
     Write(*,*) "Warning from VVVloop", p2,m12,m22
     VVVloop= 0._dp
   End if
! End if
End if
End Function VVVloop

 Complex(dp) Function SA_Gloop(p2,m12,m22)
 !-----------------------------------------------------------------------
 ! calculates the function G as defined in J. Bagger at al, Nucl.Phys.B
 ! written by Werner Porod, 4.8.1999
 ! 18.05.2001: porting to f90
 !-----------------------------------------------------------------------
 Implicit None

  Real(dp), Intent(in) :: p2,m12,m22

  If ((m12.eq.0._dp).and.(m22.eq.0._dp)) then
   SA_Gloop = p2 * SA_B0(p2,m12,m22)

  Else If (m12.eq.0._dp) then
   SA_Gloop = - SA_A0(m22) + (p2 - m22) * SA_B0(p2,m12,m22)

  Else If (m22.eq.0._dp) then
   SA_Gloop = - SA_A0(m12) + (p2 - m12) * SA_B0(p2,m12,m22)

  Else If (m12.eq.m22) then
   SA_Gloop =  - 2._dp * SA_A0(m12) + (p2 - m12 - m22) * SA_B0(p2,m12,m22)

  Else
   SA_Gloop =  - SA_A0(m12) - SA_A0(m22) + (p2 - m12 - m22) * SA_B0(p2,m12,m22)

  End If

 End Function SA_Gloop

 Complex(dp) Function SA_Hloop(p2,m12,m22)
 !-----------------------------------------------------------------------
 ! calculates the function H as defined in J. Bagger at al, Nucl.Phys.B
 ! written by Werner Porod, 4.8.1999
 ! 18.05.2001: porting to f90
 !-----------------------------------------------------------------------
 Implicit None

  Real(dp), Intent(in) :: p2,m12,m22

  SA_Hloop = 4._dp * SA_B00(p2,m12,m22) + SA_GLoop(p2,m12,m22)

 End Function SA_Hloop

Complex(dp) Function VGGloop(p2,m12,m22)
Implicit None
Real(dp),Intent(in)::p2,m12,m22

If (RXi.eq.1._dp) Then
  VGGloop=  SA_B00(p2,m12,m22)
Else
  VGGloop=  SA_B00(p2,m12*RXi,m22*RXi)

End if
End Function VGGloop

!definition of A00, B001, B111

 Complex(dp) Function A00(xm)
 Implicit None
  Real(dp), Intent(in) :: xm

  If (divonly.Eq.1) Then
   A00 = 0.25_dp * divergence * xm * xm
   Return
  End If

  A00 = (1._dp/4._dp)*A0(xm)*xm + (1._dp/8._dp)*xm*xm

 End Function A00

 Complex(dp) Function B001(xp, xm1, xm2)
 Implicit None
  Real(dp), Intent(in) :: xp, xm1, xm2
  Real(dp) :: cf1
  Complex(dp) :: ca0i(2)

  If (divonly.Eq.1) Then
   B001 = divergence * (-2._dp * xm1 - 4._dp * xm2 + xp) / 24._dp
   Return
  End If

  cf1 = xp - xm2 + xm1

  ca0i(1) = A0(xm1)
  ca0i(2) = A0(xm2)

  B001 = 2._dp * xm1 * B1(xp,xm1,xm2) &
    & - ca0i(2) + cf1 * B11(xp,xm1,xm2) &
    & - (1._dp)/(6._dp) &
    & * (2._dp * xm1 + 4._dp * xm2 - xp)

  B001 = (1._dp/8._dp) * B001

 End Function B001

 Complex(dp) Function B111(xp, xm1, xm2)
 Implicit None
  Real(dp), Intent(in) :: xp, xm1, xm2
  Real(dp) :: cf1
  Complex(dp) :: ca0i(2)

  If (divonly.Eq.1) Then
   B111 = - 0.25_dp * divergence
   Return
  End If

  cf1 = xp - xm2 + xm1

  ca0i(1) = A0(xm1)
  ca0i(2) = A0(xm2)

  B111 = ca0i(2) + 3._dp * cf1 * B11(xp,xm1,xm2) &
    & + 2._dp * xm1 * B1(xp,xm1,xm2) &
    & - (1._dp)/(6._dp) &
    & * (2._dp * xm1 + 4._dp * xm2 - xp)

  B111 = - (1._dp/4._dp) * B111
  B111 = B111/xp

 End Function B111

!derivative of B0
 Complex(dp) Function SA_DerB0(xp, xm1, xm2)
 Implicit None
  Real(dp), Intent(in) :: xp, xm1, xm2

  If ((xm1.eq.0._dp).and.(xm2.eq.0._dp)) Then
    SA_DerB0=0._dp
    Return
  End if

  If (divonly.Eq.1) Then
   SA_DerB0 = ZeroC
   Return
  End If

  If (IRdivonly) Then
!    If ((xm1.eq.Mass_Regulator_PhotonGluon**2).and.(xp.eq.xm2)) Then
   If ((EqualMasses(xm1,Mass_Regulator_PhotonGluon**2)).and.(EqualMasses(xm2,xp))) Then
    SA_DerB0 = 1._dp/(2._dp*xm2)*(0._dp - log(Mass_Regulator_PhotonGluon**2/xm2))
!    Else If ((xm2.eq.Mass_Regulator_PhotonGluon**2).and.(xp.eq.xm1)) Then
    Else If ((EqualMasses(xm2,Mass_Regulator_PhotonGluon**2)).and.(EqualMasses(xm1,xp))) Then
    SA_DerB0 = 1._dp/(2._dp*xm1)*(0._dp - log(Mass_Regulator_PhotonGluon**2/xm1))
   Else
    SA_DerB0 = 0._dp
   End if

  Else

  SA_DerB0 = DB0(xp,xm1,xm2)

  End If

 End Function SA_DerB0

 !derivative of B0
 Complex(dp) Function SA_DerB1(xp, xm1, xm2)
 Implicit None
  Real(dp), Intent(in) :: xp, xm1, xm2

!
   If ((xm1.eq.0._dp).and.(xm2.eq.0._dp)) Then
    SA_DerB1=0._dp
    Return
  End if

    If (divonly.Eq.1) Then
     SA_DerB1 = ZeroC
     Return
    End If

  If (IRdivonly) Then
!    If (xm1.eq.Mass_Regulator_PhotonGluon**2) Then
   If (EqualMasses(xm1,Mass_Regulator_PhotonGluon**2)) Then
    SA_DerB1 = ZeroC
   Else If ((EqualMasses(xm2,Mass_Regulator_PhotonGluon**2)).and.(EqualMasses(xm1,xp))) Then
!    Else If ((xm2.eq.Mass_Regulator_PhotonGluon**2).and.(xm1.eq.xp)) Then
    SA_DerB1 = -1._dp/(2._dp*xm1)*(0._dp - log(Mass_Regulator_PhotonGluon**2/xm1))
   Else
    SA_DerB1 = 0._dp
   End if

  Else
  SA_DerB1 = DerB1(xp,xm1,xm2)
  End if
 End Function SA_DerB1

!derivatives of other loop functions
Complex(dp) Function SA_DerA0(m1)
Implicit None
Real(dp), Intent(in) :: m1
SA_DerA0 = 0._dp
End Function SA_DerA0

Complex(dp) Function DerB22(p2,m12,m22)
Implicit None
Real(dp), Intent(in) :: p2,m12,m22
DerB22 = SA_DerB00(p2,m12,m22)
End Function DerB22

Complex(dp) Function SA_DerHloop(p2,m12,m22)
Implicit None
Real(dp), Intent(in) :: p2,m12,m22
SA_DerHloop =   4._dp * SA_DerB00(p2,m12,m22) + SA_DerGLoop(p2,m12,m22)
End Function SA_DerHloop

Complex(dp) Function SA_DerGloop(p2,m12,m22)
Implicit None
Real(dp), Intent(in) :: p2,m12,m22
! SA_DerGloop =  - 2._dp*(m12* SA_DerB0(p2,m12,m22) + p2*SA_DerB1(p2,m12,m22) + SA_B1(p2,m12,m22))
SA_DerGloop =  + (p2 - m12 - m22) * SA_DerB0(p2,m12,m22) + SA_B0(p2,m12,m22)
End Function SA_DerGloop

Complex(dp) Function DerFloopRXi(p2,m12,m22)
Implicit None
Real(dp), Intent(in) :: p2,m12,m22
DerFloopRXi = (-(2._dp*p2 + 2._dp*m12-m22)*SA_DerB0(p2,m12,m22)-2._dp*SA_B0(p2,m12,m22))
End Function DerFloopRXi

Complex(dp) Function DerVSSloop(p2,m12,m22)
Implicit None
Real(dp), Intent(in) :: p2,m12,m22
DerVSSloop = -4._dp*SA_DerB00(p2,m12,m22)
End Function DerVSSloop

Complex(dp) Function DerVVSloop(p2,m12,m22)
Implicit None
Real(dp), Intent(in) :: p2,m12,m22
DerVVSloop = SA_DerB0(p2,m12,m22)
End Function DerVVSloop

Complex(dp) Function DerSVVloop(p2,m12,m22)
Implicit None
Real(dp), Intent(in) :: p2,m12,m22
DerSVVloop = 4._dp*Real(SA_DerB0(p2,m12,m22),dp)
End Function DerSVVloop

Complex(dp) Function  DerVGGloop(p2,m12,m22)
Implicit None
Real(dp), Intent(in) :: p2,m12,m22
DerVGGloop = SA_DerB00(p2,m12,m22)
End Function  DerVGGloop

Complex(dp) Function DerVVVloop(p2,m12,m22)
Implicit None
Real(dp), Intent(in) :: p2,m12,m22
DerVVVloop = 10._dp*SA_DerB00(p2,m12,m22) &
     & + (m12+m22+4._dp*p2)*SA_DerB0(p2,m12,m22)+ 4._dp*SA_B0(p2,m12,m22)

End Function DerVVVloop

!derivatives of B1, B00, B11, B001, B111

 Complex(dp) Function DerB1(xp, xm1, xm2)
 Implicit None
  Real(dp), Intent(in) :: xp, xm1, xm2
  Complex(dp) :: LogA, LogB, LogC, LogD, LogE, sxm2, sxp
  Real(dp) :: mdiff

  If (divonly.Eq.1) Then
   DerB1 = ZeroC
   Return
  End If

  If (xp.Eq.0._dp) Then
    If (xm1.Ne.0._dp) Then
    mdiff = Abs((xm1-xm2)/xm1) !to avoid problems when comparing xm1 with xm2
    Else If (xm2.Ne.0._dp) Then
    mdiff = Abs((xm1-xm2)/xm2)
    Else
    mdiff = 0._dp
    End If
   If ((mdiff.Le.controlmdiff).And.(xm1.Eq.0._dp)) Then !xp=xm1=xm2=0
    write(ErrCan,*) "DerB1(0,0,0) diverges!"
   Else If (mdiff.Le.controlmdiff) Then !xp=0, xm1=xm2
    DerB1 = - 1._dp / (12._dp * xm1)
   Else If (xm1.Eq.0._dp) Then !xp=xm1=0
    DerB1 = - 1._dp / (6._dp * xm2)
   Else If (xm2.Eq.0._dp) Then !xp=xm2=0
    DerB1 = 1._dp / (3._dp * xm1)
   Else !xp=0

    mudim2=GetRenormalizationScale()

    LogA = Log(xm2/xm1)
    LogB = Log(Sqrt(xm1/xm2))
    LogC = Log(Sqrt(xm1*xm2)/mudim2)
    LogD = Log(mudim2/xm1)
    LogE = Log(mudim2/xm2)

    DerB1 = LogA * xm1**3 - 3._dp * LogA * xm1**2 * xm2 &
    & + 3._dp * LogA * xm1 * xm2**2 - LogA * xm2**3 &
    & + 8._dp * LogB * xm1**2 * xm2 + 8._dp * LogB * xm1 * xm2**2 &
    & - 4._dp * LogB * xm2**3 - 2._dp * LogC * xm1**3 &
    & + 6._dp * LogC * xm1**2 * xm2 - 6._dp * LogC * xm1 * xm2**2 &
    & + 2._dp * LogC * xm2**3 - 2._dp * LogD * xm1**3 &
    & + 4._dp * LogD * xm1**2 * xm2 - 2._dp * LogD * xm1 * xm2**2 &
    & + 2._dp * LogE * xm1**2 * xm2 - 4._dp * LogE * xm1 * xm2**2 &
    & + 2._dp * LogE * xm2**3 - 2._dp * xm1**3 - 3._dp * xm1**2 * xm2 &
    & + 6._dp * xm1 * xm2**2 - xm2**3
    DerB1 = DerB1 / (6._dp * (xm1 - xm2)**4)
   End If

  Else If (EqualMasses(xp,xm1).And.(Abs(xm2/xp).Le.(0.0001_dp))) Then
  !relevant in particular for IR divergent cases with small regulator mass xm2
  !Taylor Series of DerB1(xp,xp,xm2) up to xm2^4

   sxm2 = dsqrt(xm2)
   sxp = dsqrt(xp)

   DerB1 = 768._dp * xm2**4 + 7168._dp * xm2**3 * xp &
   & - 10584._dp * Pi * sxm2**5 * sxp**3 - 89600._dp * xm2**2 * xp**2 &
   & + 94080._dp * Pi * sxm2**3 * sxp**5 + 21504._dp * xm2 * xp**3 &
   & - 21._dp * (-6144._dp * xp**4 + 33._dp * Pi * sxm2**7 * sxp &
   & + 5120._dp * Pi * sxm2 * sxp**7) &
   & + 43008._dp * xp**2 * (xm2**2 - 3._dp * xm2 * xp + xp**2) * DLog(xm2/xp)

   DerB1 = DerB1 /(86016._dp * xp**5)

   !up to xm2^6 - makes essentially no difference
!   DerB1 = DerB1 + (131072._dp * xm2**6 - 96525._dp * Pi * sxm2**11 * sxp &
!   & + 720896._dp * xm2**5 * xp - 566280._dp * Pi * sxm2**9 * sxp**3 &
!   & - 4181760._dp * Pi * sxm2**7 * sxp**5 &
!   & - 648806400._dp * Pi * sxm2 * sxp**11 &
!   & + 126720._dp * Pi * xp**2 * (33._dp * sxm2**7 * sxp &
!   & + 5120._dp * sxm2 * sxp**7))/(519045120._dp * xp**7)

  Else !generell case

  DerB1 = (xm1 - xm2) * B0(xp,xm1,xm2) &
   & + (xm2 - xm1) * B0(0._dp,xm1,xm2) &
   & - xp * (xm1 - xm2 + xp) * SA_DerB0(xp,xm1,xm2)
  DerB1 = DerB1 / (2._dp * xp * xp)

  End If

 End Function DerB1

 Complex(dp) Function SA_DerB00(xp, xm1, xm2)
 Implicit None
  Real(dp), Intent(in) :: xp, xm1, xm2
  Complex(dp) :: VB0, VDerB0
  Complex(dp) :: LogA, LogB, LogC, LogD, LogE
  Real(dp) :: mdiff

  mudim2=GetRenormalizationScale()

  If (divonly.Eq.1) Then
   If (IRdivonly) Then
    SA_DerB00 = 0._dp
   Else
    SA_DerB00 = - divergence / 12._dp
    End if
   Return
  End If

  If (xp.Eq.0._dp) Then
    If (xm1.Ne.0._dp) Then
    mdiff = Abs((xm1-xm2)/xm1) !to avoid problems when comparing xm1 with xm2
    Else If (xm2.Ne.0._dp) Then
    mdiff = Abs((xm1-xm2)/xm2)
    Else
    mdiff = 0._dp
    End If
   If ((mdiff.Le.controlmdiff).And.(xm1.Eq.0._dp)) Then !xp=xm1=xm2=0
    write(ErrCan,*) "DerB00(0,0,0) diverges!"
   Else If (mdiff.Le.controlmdiff) Then !xp=0, xm1=xm2
    LogD = Log(mudim2/xm1)
    SA_DerB00 = -(divergence + LogD) / 12._dp
   Else If (xm1.Eq.0._dp) Then !xp=xm1=0
    LogE = Log(mudim2/xm2)
    SA_DerB00 = - (6._dp * divergence + 6._dp * LogE + 5._dp) / 72._dp
   Else If (xm2.Eq.0._dp) Then !xp=xm2=0
    LogD = Log(mudim2/xm1)
    SA_DerB00 = - (6._dp * divergence + 6._dp * LogD + 5._dp) / 72._dp
   Else !xp=0
    LogA = Log(xm2/xm1)
    LogB = Log(Sqrt(xm1/xm2))
    LogC = Log(Sqrt(xm1*xm2)/mudim2)
    LogD = Log(mudim2/xm1)
    LogE = Log(mudim2/xm2)

    SA_DerB00 = - 6._dp * divergence * xm1**3 + 18._dp * divergence * xm1**2 * xm2 &
    & -18._dp * divergence * xm1 * xm2**2 + 6._dp * divergence * xm2**3 &
    & + 6._dp * LogB * xm1**3 - 18._dp * LogB * xm1**2 * xm2 &
    & - 18._dp * LogB * xm1 * xm2**2 + 6._dp * LogB * xm2**3 &
    & + 6._dp * LogC * xm1**3 - 18._dp * LogC * xm1**2 * xm2 &
    & + 18._dp * LogC * xm1 * xm2**2 - 6._dp * LogC * xm2**3 &
    & - 5._dp * xm1**3 + 27._dp * xm1**2 * xm2 - 27._dp * xm1 * xm2**2 &
    & + 5._dp * xm2**3
    SA_DerB00 = SA_DerB00 / (72._dp * (xm1 - xm2)**3)
   End If

  Else !generell case

  VB0 = B0(xp,xm1,xm2)
  VDerB0 = SA_DerB0(xp,xm1,xm2)

  SA_DerB00 = 3._dp * (xm1 * xm1 - 2._dp * xm1 * xm2 &
   & + xm2 * xm2 - xp * xp) * VB0 &
   & - 3._dp * (xm1 - xm2) * (xm1 - xm2) * B0(0._dp,xm1,xm2) &
   & - 3._dp * xp * (xm1 * xm1 + xm2 * xm2 + xp * xp)* VDerB0 &
   & + 6._dp * (xm1 + xm2) * xp * xp * VDerB0 &
   & + 6._dp * xm1 * xm2 * xp * VDerB0 - 2._dp * xp * xp
  SA_DerB00 = SA_DerB00 / (36._dp * xp * xp)

  End If

 End Function SA_DerB00

 Complex(dp) Function DerB11(xp, xm1, xm2)
 Implicit None
  Real(dp), Intent(in) :: xp, xm1, xm2
  Complex(dp) :: VB0, VDerB0
  Complex(dp) :: LogA, LogB, LogC, LogD, LogE
  Real(dp) :: mdiff

  mudim2=GetRenormalizationScale()

  If (divonly.Eq.1) Then
   DerB11 = ZeroC
   Return
  End If

  If (xp.Eq.0._dp) Then
    If (xm1.Ne.0._dp) Then
    mdiff = Abs((xm1-xm2)/xm1) !to avoid problems when comparing xm1 with xm2
    Else If (xm2.Ne.0._dp) Then
    mdiff = Abs((xm1-xm2)/xm2)
    Else
    mdiff = 0._dp
    End If
   If ((mdiff.Le.controlmdiff).And.(xm1.Eq.0._dp)) Then !xp=xm1=xm2=0
    write(ErrCan,*) "DerB11(0,0,0) diverges!"
   Else If (mdiff.Le.controlmdiff) Then !xp=0, xm1=xm2
    DerB11 = 1._dp / (20._dp * xm1)
   Else If (xm1.Eq.0._dp) Then !xp=xm1=0
    DerB11 = 1._dp / (12._dp * xm2)
   Else If (xm2.Eq.0._dp) Then !xp=xm2=0
    DerB11 = 1._dp / (4._dp * xm1)
   Else !xp=0
    LogA = Log(xm2/xm1)
    LogB = Log(Sqrt(xm1/xm2))
    LogC = Log(Sqrt(xm1*xm2)/mudim2)
    LogD = Log(mudim2/xm1)
    LogE = Log(mudim2/xm2)

    DerB11 = -3._dp * LogA * xm1**4 + 12._dp * LogA * xm1**3 * xm2 &
    & - 18._dp * LogA * xm1**2 * xm2**2 + 12._dp * LogA * xm1 * xm2**3 &
    & - 3._dp * LogA * xm2**4 - 30._dp * LogB * xm1**3 * xm2 &
    & - 54._dp * LogB * xm1**2 * xm2**2 + 54._dp * LogB * xm1 * xm2**3 &
    & - 18._dp * LogB * xm2**4 + 6._dp * LogC * xm1**4 &
    & - 30._dp * LogC * xm1**3 * xm2 + 54._dp * LogC * xm1**2 * xm2**2 &
    & - 42._dp * LogC * xm1 * xm2**3 + 12._dp * LogC * xm2**4 &
    & + 6._dp * LogD * xm1**4 - 18._dp * LogD * xm1**3 * xm2 &
    & + 18._dp * LogD * xm1**2 * xm2**2 - 6._dp * LogD * xm1 * xm2**3 &
    & - 12._dp * LogE * xm1**3 * xm2 + 36._dp * LogE * xm1**2 * xm2**2 &
    & - 36._dp * LogE * xm1 * xm2**3 + 12._dp * LogE * xm2**4 &
    & + 6._dp * xm1**4 + 20._dp * xm1**3 * xm2 - 36._dp * xm1**2 * xm2**2 &
    & + 12._dp * xm1 * xm2**3 - 2._dp * xm2**4
    DerB11 = DerB11 / (24._dp * (xm1 - xm2)**5)
   End If

  Else !generell case

  VB0 = B0(xp,xm1,xm2)
  VDerB0 = SA_DerB0(xp,xm1,xm2)

   DerB11 = 4._dp * xm1**2 * SA_B0(0._dp,xm1,xm2) &
    & + 4._dp * xm2**2 * B0(0._dp,xm1,xm2) + 2._dp * xm1 * xp * B0(0._dp,xm1,xm2) &
    & - 2._dp * xm2 * xp * B0(0._dp,xm1,xm2) - 8._dp * xm1 * xm2 * B0(0._dp,xm1,xm2) &
    & - 2._dp * xp * A0(xm2) - 4._dp * VB0 * xm1**2 + 8._dp * VB0 * xm1 * xm2 &
    & - 2._dp * VB0 * xm1 * xp - 4._dp * VB0 * xm2**2 + 4._dp * VB0 * xm2 * xp &
    & + 2._dp * VDerB0 * xm1**2 * xp - 4._dp * VDerB0 * xm1 * xm2 * xp &
    & + 2._dp * VDerB0 * xm1 * xp**2 + 2._dp * VDerB0 * xm2**2 * xp &
    & - 4._dp * VDerB0 * xm2 * xp**2 + 2._dp * VDerB0 * xp**3 + xm1 * xp + xm2 * xp
   DerB11 = DerB11 / (6._dp * xp * xp * xp)

  End If

 End Function DerB11
!
!  Complex(dp) Function DerB001(xp, xm1, xm2)
!  Implicit None
!   Real(dp), Intent(in) :: xp, xm1, xm2
!   Complex(dp) :: VB0, VDerB0
!   Complex(dp) :: LogA, LogB, LogC, LogD, LogE
!   Real(dp) :: mdiff
!
!   If (divonly.Eq.1) Then
!    DerB001 = divergence / 24._dp
!    Return
!   End If
!
!   mudim2=GetRenormalizationScale()
!
!   If (xp.Eq.0._dp) Then
!     If (xm1.Ne.0._dp) Then
!     mdiff = Abs((xm1-xm2)/xm1) !to avoid problems when comparing xm1 with xm2
!     Else If (xm2.Ne.0._dp) Then
!     mdiff = Abs((xm1-xm2)/xm2)
!     Else
!     mdiff = 0._dp
!     End If
!    If ((mdiff.Le.controlmdiff).And.(xm1.Eq.0._dp)) Then !xp=xm1=xm2=0
!     write(ErrCan,*) "DerB001(0,0,0) diverges!"
!    Else If (mdiff.Le.controlmdiff) Then !xp=0, xm1=xm2
!     LogD = Log(mudim2/xm1)
!     DerB001 = (divergence + LogD)/ 24._dp
!    Else If (xm1.Eq.0._dp) Then !xp=xm1=0
!     LogE = Log(mudim2/xm2)
!     DerB001 = (12._dp * divergence + 12._dp * LogE + 7._dp) / 288._dp
!    Else If (xm2.Eq.0._dp) Then !xp=xm2=0
!     LogD = Log(mudim2/xm1)
!     DerB001 = (12._dp * divergence + 12._dp * LogD + 13._dp) / 288._dp
!    Else !xp=0
!     LogA = Log(xm2/xm1)
!     LogB = Log(Sqrt(xm1/xm2))
!     LogC = Log(Sqrt(xm1*xm2)/mudim2)
!     LogD = Log(mudim2/xm1)
!     LogE = Log(mudim2/xm2)
!
!     DerB001 = 12._dp * divergence * xm1**4 - 48._dp * divergence * xm1**3 * xm2 &
!     & + 72._dp * divergence * xm1**2 * xm2**2 - 48._dp * divergence * xm1 * xm2**3 &
!     & + 12._dp * divergence * xm2**4 + 2._dp * LogA * xm1**4 - 12._dp * LogA * xm1**2 * xm2**2 &
!     & + 16._dp * LogA * xm1 * xm2**3 - 6._dp * LogA * xm2**4 - 12._dp * LogB * xm1**4 &
!     & + 48._dp * LogB * xm1**3 * xm2 + 48._dp * LogB * xm1**2 * xm2**2 &
!     & - 12._dp * LogB * xm2**4 - 16._dp * LogC * xm1**4 + 48._dp * LogC * xm1**3 * xm2 &
!     & - 48._dp * LogC * xm1**2 * xm2**2 + 16._dp * LogC * xm1 * xm2**3 &
!     & - 4._dp * LogD * xm1**4 + 12._dp * LogD * xm1**2 * xm2**2 &
!     & - 8._dp * LogD * xm1 * xm2**3 + 12._dp * LogE * xm1**2 * xm2**2 &
!     & - 24._dp * LogE * xm1 * xm2**3 + 12._dp * LogE * xm2**4 &
!     & + 13._dp * xm1**4 - 88._dp * xm1**3 * xm2 + 108._dp * xm1**2 * xm2**2 &
!     & - 40._dp * xm1 * xm2**3 + 7._dp * xm2**4
!     DerB001 = DerB001 / (288._dp * (xm1 - xm2)**4)
!    End If
!
!   Else !generell case
!
!   VB0 = B0(xp,xm1,xm2)
!   VDerB0 = DerB0(xp,xm1,xm2)
!
!   DerB001 = 6._dp * (xm1 - xm2) * (2._dp * xm1**2 - xm1 * (4._dp * xm2 + xp) &
!    &         + 2._dp * xm2 * (xm2-xp)) * B0(0._dp,xm1,xm2) &
!    & + 6._dp * xp * (xm2 - xm1) * A0(xm2) - 12._dp * VB0 * xm1**3 &
!    & + 36._dp * VB0 * xm1**2 * xm2 + 6._dp * VB0 * xm1**2 * xp &
!    & - 36._dp * VB0 * xm1 * xm2**2 + 12._dp * VB0 * xm1 * xm2 * xp &
!    & + 12._dp * VB0 * xm2**3 - 18._dp * VB0 * xm2**2 * xp &
!    & + 6._dp * VB0 * xp**3 + 6._dp * VDerB0 * xm1**3 * xp &
!    & - 18._dp * VDerB0 * xm1**2 * xm2 * xp - 6._dp * VDerB0 * xm1**2 * xp**2 &
!    & + 18._dp * VDerB0 * xm1 * xm2**2 * xp - 12._dp * VDerB0 * xm1 * xm2 * xp**2 &
!    & - 6._dp * VDerB0 * xm1 * xp**3 - 6._dp * VDerB0 * xm2**3 * xp &
!    & + 18._dp * VDerB0 * xm2**2 * xp**2 - 18._dp * VDerB0 * xm2 * xp**3 &
!    & + 6._dp * VDerB0 * xp**4 + 3._dp * xm1**2 * xp - 3._dp * xm2**2 * xp + 4._dp * xp**3
!   DerB001 = DerB001 / (144._dp * xp * xp * xp)
!
!   End If
!
!  End Function DerB001
!
!  Complex(dp) Function DerB111(xp, xm1, xm2)
!  Implicit None
!   Real(dp), Intent(in) :: xp, xm1, xm2
!   Complex(dp) :: VB0, VDerB0
!   Complex(dp) :: LogA, LogB, LogC, LogD, LogE
!   Real(dp) :: mdiff
!
!   If (divonly.Eq.1) Then
!    DerB111 = ZeroC
!    Return
!   End If
!
!   If (xp.Eq.0._dp) Then
!     If (xm1.Ne.0._dp) Then
!     mdiff = Abs((xm1-xm2)/xm1) !to avoid problems when comparing xm1 with xm2
!     Else If (xm2.Ne.0._dp) Then
!     mdiff = Abs((xm1-xm2)/xm2)
!     Else
!     mdiff = 0._dp
!     End If
!    If ((mdiff.Le.controlmdiff).And.(xm1.Eq.0._dp)) Then !xp=xm1=xm2=0
!     write(ErrCan,*) "DerB111(0,0,0) diverges!"
!    Else If (mdiff.Le.controlmdiff) Then !xp=0, xm1=xm2
!     DerB111 = - 1._dp / (30._dp * xm1)
!    Else If (xm1.Eq.0._dp) Then !xp=xm1=0
!     DerB111 = - 1._dp / (20._dp * xm2)
!    Else If (xm2.Eq.0._dp) Then !xp=xm2=0
!     DerB111 = - 1._dp / (5._dp * xm1)
!    Else !xp=0
!     LogA = Log(xm2/xm1)
!     LogB = Log(Sqrt(xm1/xm2))
!     LogC = Log(Sqrt(xm1*xm2)/mudim2)
!     LogD = Log(mudim2/xm1)
!     LogE = Log(mudim2/xm2)
!
!     DerB111 = 15._dp * LogA  * xm1**5 + 21._dp * LogA * xm1**4 * xm2 &
!     & - 9._dp * LogA * xm1**3 * xm2**2 - 249._dp * LogA * xm1**2 * xm2**3 &
!     & + 366._dp * LogA * xm1 * xm2**4 - 144._dp * LogA * xm2**5 &
!     & + 120._dp * LogB * xm1**4 * xm2 - 360._dp * LogB * xm1**2 * xm2**3 &
!     & + 1200._dp * LogB * xm1 * xm2**4 - 720._dp * LogB * xm2**5 &
!     & - 30._dp * LogC * xm1**5 + 102._dp * LogC * xm1**4 * xm2 &
!     & - 138._dp * LogC * xm1**3 * xm2**2 + 522._dp * LogC * xm1**2 * xm2**3 &
!     & - 888._dp * LogC * xm1 * xm2**4 + 432._dp * LogC * xm2**5 &
!     & - 30._dp * LogD * xm1**5 - 30._dp * LogD * xm1**4 * xm2 &
!     & - 60._dp * LogD * xm1**3 * xm2**2 + 330._dp * LogD * xm1**2 * xm2**3 &
!     & - 210._dp * LogD * xm1 * xm2**4 + 132._dp * LogE * xm1**4 * xm2 &
!     & - 78._dp * LogE * xm1**3 * xm2**2 + 192._dp * LogE * xm1**2 * xm2**3 &
!     & - 678._dp * LogE * xm1 * xm2**4 + 432._dp * LogE * xm2**5 &
!     & - 24._dp * xm1**5 - 130._dp * xm1**4 * xm2 + 240._dp * xm1**3 * xm2**2 &
!     & - 120._dp * xm1**2 * xm2**3 + 40._dp * xm1 * xm2**4 - 6._dp * xm2**5
!     DerB111 = DerB111 / (120._dp * (xm1 - xm2)**6)
!    End If
!
!   Else !generell case
!
!   VB0 = B0(xp,xm1,xm2)
!   VDerB0 = DerB0(xp,xm1,xm2)
!
!   DerB111 =3._dp * (xm1-xm2) * (3._dp * xm1**2 - 6._dp * xm1 * xm2 + 2._dp * xm1 * xp &
!    & + 3._dp * xm2**2 - 4._dp * xm2 * xp + xp**2)*B0(0._dp,xm1,xm2) &
!    & - 6._dp * xp *A0(xm2)*(xm1-xm2+xp) - 9._dp * VB0 * xm1**3 &
!    & + 27._dp * VB0 * xm1**2 * xm2 - 6._dp * VB0 * xm1**2 * xp &
!    & - 27._dp * VB0 * xm1 * xm2**2 + 24._dp * VB0 * xm1 * xm2 * xp &
!    & - 3._dp * VB0 * xm1 * xp**2 + 9._dp * VB0 * xm2**3 &
!    & - 18._dp * VB0 * xm2**2 * xp + 9._dp * VB0 * xm2 * xp**2 &
!    & + 3._dp * VDerB0 * xm1**3 * xp - 9._dp * VDerB0 * xm1**2 * xm2 * xp &
!    & + 3._dp * VDerB0 * xm1**2 * xp**2 + 9._dp * VDerB0 * xm1 * xm2**2 * xp &
!    & - 12._dp * VDerB0 * xm1 * xm2 * xp**2 + 3._dp * VDerB0 * xm1 * xp**3 &
!    & - 3._dp * VDerB0 * xm2**3 * xp + 9._dp * VDerB0 * xm2**2 * xp**2 &
!    & - 9._dp * VDerB0 * xm2 * xp**3 + 3._dp * VDerB0 * xp**4 &
!    & + 3._dp * xm1**2 * xp + 2._dp * xm1 * xp**2 - 3._dp * xm2**2 * xp &
!    & + 4._dp * xm2 * xp**2
!   DerB111 = - DerB111 / (12._dp * xp * xp * xp * xp)
!
!   End If
!
!  End Function DerB111

!special passarino veltman integrals with zero masses:
!IR145 = argument 1,4 and 5 are zero
!example: C0_IR145(xp1,xp2,xm) = C0(0,xp1,xp2,0,0,xm)
!notation equivalent to looptools = ff package

 Complex(dp) Function C0_IR145(xp1, xp2, xm)
 Implicit None
  Real(dp), Intent(in) :: xp1, xp2, xm
!  Real(dp) :: s1
  Complex(dp) :: logspecial

  If (xp1.Eq.xp2) Then !xp1 = xp2

   If (xp1.Eq.xm) Then !xp1 = xp2 = xm

    !Call WriteLFerror(3) !divergent
    write(ErrCan,*) "Divergent Passarino Veltman, C0_IR145!"
    C0_IR145 = Huge(1._dp)

   Else !xp1 = xp2 unequal xm

!    calculating log(1-xp2/xm)
!     s1 = 1._dp - xm/xp2
!    If ( Abs(s1) .Lt. xloss ) Then
!     logspecial = Log1minusX(s1)
!    Else
    logspecial = Log(Abs(1._dp - xm/xp2)) + (xm/xp2) * Log(Abs(1._dp/(xp2/xm - 1._dp)))
    If (xp2.Gt.xm) logspecial = logspecial + 2._dp * Cmplx(0._dp,pi,dp) * (1._dp - xm/xp2)

    C0_IR145 = logspecial / (xp2 - xm)

   End If

  Else !the two momenta are not equal

  C0_IR145 = Cget("C0  ",0._dp,xp1,xp2,0._dp,0._dp,xm)

  End If

 End Function C0_IR145

 Complex(dp) Function C2_IR145(xp1, xp2, xm)
 Implicit None
  Real(dp), Intent(in) :: xp1, xp2, xm
  Real(dp) :: s1
  Complex(dp) :: logspecial

  If (xp1.Eq.xp2) Then !xp1 = xp2

   If (xp1.Eq.xm) Then !xp1 = xp2 = xm

    !Call WriteLFerror(3) !divergent
    write(ErrCan,*) "Divergent Passarino Veltman, C2_IR145!"
    C2_IR145 = Huge(1._dp)

   Else !xp1 = xp2 unequal xm

    !calculating log(1-xp2/xm)
     s1 = xp2/xm
    If ( Abs(s1) .Lt. xloss ) Then
     logspecial = Log1minusX(s1)
    Else
     logspecial = Log(Abs(1._dp - s1))
    End If
    If (xp2.Gt.xm) logspecial = logspecial + Cmplx(0._dp,pi,dp)

    C2_IR145 = - xm * logspecial - xp2
    C2_IR145 = C2_IR145 / (xp2 * xp2)

   End If

  Else !the two momenta are not equal

  C2_IR145 = B0(xp1,0._dp,xm) - B0(xp2,0._dp,xm)
  C2_IR145 = C2_IR145 / (xp1 - xp2)

  End If

 End Function C2_IR145

 Real(dp) Function SetRenormalizationScaleAdd(mu2_in)
 Implicit None
  Real(dp), Intent(in) :: mu2_in

  SetRenormalizationScaleAdd = mudim2
  mudim2 = mu2_in

 End Function SetRenormalizationScaleAdd

 Real(dp) Function SetDivergenceAdd(diver_in)
 !-----------------------------------------------------------------------
 ! resets divergence
 !-----------------------------------------------------------------------
 Implicit None
  Real(dp), Intent(in) :: diver_in

  SetDivergenceAdd = divergence
  divergence = diver_in

 End Function SetDivergenceAdd

 Integer Function SetDivonlyAdd(div_in)
 !-----------------------------------------------------------------------
 ! take only divergent parts
 !-----------------------------------------------------------------------
 Implicit None
  Integer, Intent(in) :: div_in

  SetDivonlyAdd = divonly
  divonly = div_in

  If (divonly.eq.1) rMS=0._dp ! remove constant parts in MS scheme

 End Function SetDivonlyAdd

!

End Module AddLoopFunctions


       subroutine Iteration(NLEV,NCONS,AZ_P,AN_P,conv)
       use VAPHFB_Par
       implicit none
       integer III,ii,jj,imm,mm
       integer NLEV,NCONS
       real*8 H_0,AZ_0,AN_0,AZ_P,AN_P,graddd
       real*8 H20_0(1:HO%NLEV*HO%NLEV)
       real*8, DIMENSION(:), allocatable :: ACONS_20 !(NCONS*NLEV*NLEV)
       real*8 bet2t,gam2t,P00t
       real*8 conv,zg_0,EPNP
!      ...........................................................
       print *, ' -> Iteration '
       print *, '    ......... '
       write(*,'(A4,100A12)') 'It','Grad','E_HFB', 'Prot',&
       &'Neut','bet2','gam2','P00','cr_freq'

        if(NCONS.gt.NCONSMAX) &
       & stop 'Error: NCONS is larger than NCONSMAX!'

!      ......... allocate memory
        if(.NOT. ALLOCATED(HFB%RO_0))      &
       &  ALLOCATE(HFB%RO_0(1:HO%NLEV,1:HO%NLEV))
        if(.NOT. ALLOCATED(HFB%Akapa10_0)) &
       & ALLOCATE(HFB%Akapa10_0(1:HO%NLEV,1:HO%NLEV))
        if(.NOT. ALLOCATED(HFB%Akapa01_0)) &
       & ALLOCATE(HFB%Akapa01_0(1:HO%NLEV,1:HO%NLEV))
        if(.NOT. ALLOCATED(ACONS_20))    &
       & ALLOCATE(ACONS_20(1:HO%NLEV*HO%NLEV*NCONS))

!        write(*,'(a,f6.3,a)')' Iteration: allocated ACONS_20 &
!     &  with memory:',sizeof(ACONS_20)/(1024*1024*1024.0),'G byte ..'

!       ........................ start the iteration
        DO III=1,Input%itermax
!       .........................Initialitation
        H_0 =zero
        AN_0=zero
        AZ_0=zero

        do jj=1,NLEV**2
           H20_0(jj)=zero
        end do

!        write(*,*) ' ..... Rho and kappa ...'
        call UV2Density(HFB%U0,HFB%V0,HFB%RO_0,HFB%Akapa10_0,HFB%Akapa01_0,HO%NLEV)

        call N_0(HFB%V0,AN_0,AZ_0,HO%NLEV)

        call betgamm_0(HFB%V0,bet2t,gam2t,HO%NLEV)

        call Qpair_0(cME1B%P_00_10_me,HFB%Akapa10_0,HFB%Akapa01_0,P00t,HO%NLEV)

       do mm=1,NCONS-ipair_const
           imm=(mm-1)*HO%NLEV**2+1
           call F_20(HFB%U0,HFB%V0,ACONS_ME(imm),ACONS_20(imm),NLEV)
        end do

        do mm=NCONS-ipair_const+1,NCONS
           imm=(mm-1)*HO%NLEV**2+1
           call G_20(HFB%U0,HFB%V0,ACONS_ME(imm),ACONS_20(imm),NLEV)
        end do

!        write(*,*) '.... NZPROJ_VAP ........'
!       ........................................
         call NZPROJ_VAP(AZ_P,AN_P,EPNP,H20_0,PNP%NFOM,NLEV)

        if (III.eq.1) then
           do mm=1,NCONS
              alag_0(mm)=0.d0
           end do
        end if

!      .....

!        write(*,*) '.... LAGRANGE_MULT ........'
        call LAGRANGE_MULT(NCONS,H20_0,ACONS_20,alag_0,alag,NLEV)

        HFB%EFermi(0) =  alag_0(1)
        HFB%EFermi(1) =  alag_0(2)
!.......Gradient method.......

        do jj=1,NLEV**2
          graddd=0.d0
          do mm=1,NCONS
           imm=(mm-1)*NLEV**2
           graddd=graddd+alag(mm)*ACONS_20(imm+jj)
          end do

!          HFB%Thoul_0(jj)=  -Input%eta*(H20_0(jj)-graddd) !aqui

!    ... momentum
!    https://distill.pub/2017/momentum/
!          HFB%Thoul_0(jj)=  -Input%eta*(HFB%alpha*HFB%Thoul_old(jj)+H20_0(jj)-graddd) !aqui
          HFB%Thoul_0(jj)=  -Input%eta*(H20_0(jj)-graddd) !aqui

        end do

        call DGEMM ('t','n',NLEV,NLEV,NLEV,1.d0,HFB%Thoul_0,NLEV,HFB%Thoul_0,NLEV,0.d0,HFB%GRAD_0,NLEV)

        zg_0=zero
        do ii=1,NLEV
           zg_0=zg_0+HFB%GRAD_0(ii,ii)
        end do

        conv=zg_0/(Input%eta**2)

        write(*,1001),III,conv,EPNP+H%E0,AZ_P,AN_P,bet2t,gam2t,P00t,alag_0(9) !shared%cf
        if (conv.le.Input%tolgrad) then
           if(Input%bkick.lt.CHOP) then
             print*,'CALCULATION CONVERGED'
             print*,' '
             exit
!      ..........................................
            elseif(Input%bkick.gt.CHOP .and. HFB%iconv.gt.0) then  ! with a nonzero kick stepsize
              print*,'CALCULATION CONVERGED'
              print*,' '
              exit
!     .................... after the first convergence
            elseif(Input%bkick.gt.CHOP .and. HFB%iconv.eq.0) then  ! with a nonzero kick stepsize
              HFB%Et = EPNP+H%E0    ! record the energy
              print*,' kick the state with a Delta_bet2=',Input%bkick
!     .................... record the U,V
              do ii=1,NLEV
              do jj=1,NLEV
                 HFB%Ut(jj,ii)=HFB%U0(jj,ii)
                 HFB%Vt(jj,ii)=HFB%V0(jj,ii)
              end do
              end do
! add a kick onto Q20 after the first convergence
              if(Const%iQB.eq.1) acons_mv(3) = qtc20+Input%bkick/Const%Q2BA
              if(Const%iQB.ne.1) acons_mv(3) = qtc20+Input%bkick
              HFB%iconv   = HFB%iconv + 1
           endif
        endif
1001    format(I4,100F12.4)

        call adjust_cons(NCONS,HFB%U0,HFB%V0,HFB%Thoul_0,HFB%U1,HFB%V1,HFB%ALcholes,NLEV)

!      update U,V
        do ii=1,NLEV
         do jj=1,NLEV
            HFB%U0(jj,ii)=HFB%U1(jj,ii)
            HFB%V0(jj,ii)=HFB%V1(jj,ii)
         end do
        end do

        do jj=1,NLEV**2
            HFB%Thoul_old(jj)=HFB%Thoul_0(jj)
            HFB%Thoul_0(jj)  =zero
        end do
!      ........................ print out intermediate wave functions
!        if(mod(III,20).eq.0) then
        !writing the wave function of check point
!         print *, 'write check point file ..'
!         do ii=1,NLEV
!          do jj=1,NLEV
!             write(12,*) HFB%U0(ii,jj),HFB%V0(ii,jj)
!          end do
!         end do
!         close(12)
!        endif

        if(III.eq.Input%itermax) then
         print*,'MAXIMUM NUMBER OF ITERATIONS REACHED'
         print*,' '
        end if

!       .................................................
        END DO
        !END OF THE GRADIENT ITERATIONS
!       .................................................
        if(Input%bkick.gt.CHOP .and. HFB%Et .lt. EPNP+H%E0) then ! keep the U,V of the first convergence
        print *,' recover the wave function from the first convergence.'
        do ii=1,NLEV
         do jj=1,NLEV
            HFB%U0(jj,ii)=HFB%Ut(jj,ii)
            HFB%V0(jj,ii)=HFB%Vt(jj,ii)
         end do
        end do
        endif

       return
       end


! { dg-do compile }
! PR38709 - ICE-on-invalid on zero-sized array in init-expr.

  INTEGER, PARAMETER :: a(1) = (/ 1 /)
  INTEGER, PARAMETER :: i = a(shape(1))   ! { dg-error "Incompatible ranks" }
END


      subroutine ch_rtpath

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    this subroutine routes bacteria through the stream network

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~
!!    name             |units       |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    hrchwtr(:)       |m^3 H2O     |water stored in reach at beginning of hour
!!    rch_bactlp(:)    |# cfu/100ml |less persistent bacteria stored in reach
!!    rch_bactp(:)     |# cfu/100ml |persistent bacteria stored in reach
!!    rchwtr           |m^3 H2O     |water stored in reach at beginning of day
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    hbactlp(:)   |# cfu/100mL  |less persistent bacteria in reach/outflow
!!                               |during hour
!!    hbactp(:)    |# cfu/100mL  |persistent bacteria in reach/outflow during
!!                               |hour
!!    rch_bactlp(:)|# cfu/100ml  |less persistent bacteria in reach/outflow
!!                               |at end of day
!!    rch_bactp(:) |# cfu/100ml  |persistent bacteria in reach/outflow at end
!!                               |of day
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ii          |none          |counter
!!    initlp      |# cfu/100mL   |bacteria concentration in reach at beginning
!!                               |of hour (less persistent)
!!    initp       |# cfu/100mL   |bacteria concentration in reach at beginning
!!                               |of hour (persistent)
!!    jrch        |none          |reach number
!!    netwtr      |m^3 H2O       |net amount of water in reach during time step
!!    tday        |day           |routing time for the reach
!!    totbactlp   |10^4 cfu      |mass less persistent bacteria
!!    totbactp    |10^4 cfu      |mass persistent bacteria
!!    wtmp        |deg C         |temperature of water in reach
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    Intrinsic: Exp, Max
!!    SWAT: Theta

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use basin_module
      use time_module
      use pathogen_data_module
      use channel_module
      use hydrograph_module, only : ob
      use climate_module
      use constituent_mass_module

      implicit none

      real, external :: Theta  !units         |description
      real :: path_tot         !10^4 cfu      |mass persistent bacteria
      real :: totbactlp        !10^4 cfu      |mass less persistent bacteria
      real :: netwtr           !m^3 H2O       |net amount of water in reach during time step
      real :: initlp           !# cfu/100mL   |bacteria concentration in reach at beginning
      real :: initp            !# cfu/100mL   |bacteria concentration in reach at beginning
      real :: tday             !day           |routing time for the reach
      real :: wtmp             !deg C         |temperature of water in reach
      real :: rchwtr           !m^3 H2O       |water stored in reach at beginning of day
      integer :: iwst          !units         |description
      integer :: ipath         !none          |pathogen counter
      integer :: ipath_db      !none          |pathogen number from data file
      integer :: isp_ini       !none          |soil-plant initialization number from data fil
      integer :: jrch          !none          |reach number
      integer :: iob           !none          |
      integer :: icmd          !none          |

      if (rtwtr > 0. .and. rchdep > 0.) then

        wtmp = 5.0 + 0.75 * wst(iwst)%weat%tave
        if (wtmp <= 0.) wtmp = 0.1

        do ipath = 1, cs_db%num_paths
          !isp_ini = hru(ihru)%dbs%soil_plant_init
          !ipath_db = sol_plt_ini(isp_ini)%path

          !! total pathogen mass in reach
          path_tot = obcs(iob)%hd(1)%path(ipath) * ob(icmd)%hin%flo + ch(jrch)%bactp * rchwtr

          !! compute pathogen die-off
          tday = rttime / 24.0            !! calculate flow duration
          if (tday > 1.0) tday = 1.0
          path_tot = path_tot * Exp(-Theta(path_db(ipath)%do_stream, path_db(ipath)%t_adj,wtmp) * tday)
          path_tot = Max(0., path_tot)

          !! new concentration
          netwtr = ob(icmd)%hin%flo  + rchwtr

  	      !! change made by CS while running region 4; date 2 jan 2006
	      if (path_tot < 1.e-6) path_tot = 0.0
          if (netwtr >= 1.) then
            ch_water(jrch)%path(ipath) = path_tot / netwtr
          else
            ch_water(jrch)%path(ipath) = 0.
          end if

        end do
      end if

      return
      end subroutine ch_rtpath


text:
module mo_test_loop

  use mo_kind, only: i4, i8, sp, dp

  implicit none

  public :: n_save_state
  public :: loop

  interface loop
     module procedure loops, loopd
  end interface loop

  private

  integer(i4), parameter :: n_save_state = 132_i4

CONTAINS

  subroutine loops(seed, SingleIntegerRN, save_state)

    implicit none

    integer(i4),                                    intent(in)    :: seed
    integer(i4),                                    intent(out)   :: SingleIntegerRN
    integer(i4), optional, dimension(n_save_state), intent(inout) :: save_state

    integer(i4)        :: wlen, r, s, a, b, c, d

    integer(i4), save  :: w
    integer(i4), save  :: x(0:127)                 ! x(0) ... x(r-1)
    integer(i4)        :: weyl = 1640531527_i4     !Z'61C88647'       ! Hexadecimal notation
    ! integer(i4)        :: weyl = 1_i4     !Z'61C88647'       ! Hexadecimal notation
    integer(i4)        :: t, v
    integer(i4), save  :: i = -1                   ! i<0 indicates first call
    integer(i4)        :: k

    !$omp   threadprivate(x,i,w)

    wlen = 32
    r = 128
    s = 95
    a = 17
    b = 12
    c = 13
    d = 15

    if ( present(save_state) .and. (seed .eq. 0) ) then
       x(0:r-1) = save_state(1:r)
       i        = save_state(r+1)
       w        = save_state(r+2)
    end if

    If ((i .lt. 0) .or. (seed .ne. 0)) then     ! Initialization necessary
       If (seed .ne. 0) then                   ! v must be nonzero
          v = seed
       else
          v = NOT(seed)
       end if

       do k=wlen, 1, -1                          ! Avoid correlations for close seeds
          ! This recurrence has period of 2^32-1
          v = IEOR(v, ISHFT(v,  13))
          v = IEOR(v, ISHFT(v, -17))
          v = IEOR(v, ISHFT(v,   5))
       end do

       ! Initialize circular array
       w = v
       do k=0,r-1
          w = w + weyl
          v = IEOR(v,ISHFT(v,13))
          v = IEOR(v,ISHFT(v,-17))
          v = IEOR(v,ISHFT(v, 5))
          x(k) = v + w
       end do

       ! Discard first 4*r results (Gimeno)
       i = r-1
       do k = 4*r,1,-1
          i = IAND(i+1,r-1)
          t = x(i)
          v = x(IAND(i+(r-s),r-1))
          t = IEOR(t,ISHFT(t,a))
          t = IEOR(t,ISHFT(t,-b))
          v = IEOR(v,ISHFT(v,c))
          v = IEOR(v,IEOR(t,ISHFT(v,-d)))
          x(i) = v
       end do
    end if ! end of initialization

    ! Apart from initialization (above), this is the generator
    i = IAND(i+1,r-1)
    t = x(i)
    v = x(IAND(i+(r-s),r-1))
    t = IEOR(t,ISHFT(t,a))
    t = IEOR(t,ISHFT(t,-b))
    v = IEOR(v,ISHFT(v,c))
    v = IEOR(v,IEOR(t,ISHFT(v,-d)))
    x(i) = v

    w = w + weyl

    SingleIntegerRN = v+w

    if( present(save_state) ) then
       save_state(1:r)   = x(0:r-1)
       save_state(r+1)   = i
       save_state(r+2)   = w
       if ((r+3) <= n_save_state) save_state(r+3:n_save_state) = 0
    end if

  end subroutine loops

  subroutine loopd(seed,DoubleRealRN,save_state)

    implicit none

    integer(i8),                                    intent(in)    :: seed
    real(DP),                                       intent(out)   :: DoubleRealRN
    integer(i8), optional, dimension(n_save_state), intent(inout) :: save_state

    integer(i8)        :: wlen, r, s, a, b, c, d

    integer(i8), save  :: w
    integer(i8), save  :: x(0:63)                  ! x(0) ... x(r-1)
    integer(i8)        :: weyl = 7046029254386353131_i8
    integer(i8)        :: t,v, tmp
    integer(i8), save  :: i = -1                   ! i<0 indicates first call
    integer(i8)        :: k

    real(DP)            :: t53 = 1.0_DP/9007199254740992.0_DP                     ! = 0.5^53 = 1/2^53

    !$omp   threadprivate(x,i,w)

    ! produces a 53bit Integer Random Number (0...9 007 199 254 740 992) and
    ! scales it afterwards to (0.0,1.0)

    wlen = 64_i8
    r = 64_i8
    s = 53_i8
    a = 33_i8
    b = 26_i8
    c = 27_i8
    d = 29_i8

    if ( present(save_state) .and. (seed .eq. 0_i8) ) then
       x(0:r-1)  = save_state(1:r)
       i         = save_state(r+1)
       w         = save_state(r+2)
    end if

    If ((i .lt. 0) .or. (seed .ne. 0)) then     ! Initialization necessary
       If (seed .ne. 0) then                   ! v must be nonzero
          v = seed
       else
          v = NOT(seed)
       end if

       do k=wlen,1,-1                          ! Avoid correlations for close seeds
          ! This recurrence has period of 2^64-1
          v = IEOR(v,ISHFT(v,7))
          v = IEOR(v,ISHFT(v,-9))
       end do

       ! Initialize circular array
       w = v
       do k=0,r-1
          ! w = w + weyl
          if (w < 0_i8) then
             w = w + weyl
          else if ((huge(w) - w) > weyl) then
             w = w + weyl
          else
             tmp = -(huge(w) - w - weyl)
             w =  tmp - huge(w) - 2_i8
          endif
          v = IEOR(v,ISHFT(v,7))
          v = IEOR(v,ISHFT(v,-9))
          x(k) = v + w
       end do

       ! Discard first 4*r results (Gimeno)
       i = r-1
       do k = 4*r,1,-1
          i = IAND(i+1, r-1)
          t = x(i)
          v = x(IAND(i+(r-s),r-1))
          t = IEOR(t,ISHFT(t,a))
          t = IEOR(t,ISHFT(t,-b))
          v = IEOR(v,ISHFT(v,c))
          v = IEOR(v,IEOR(t,ISHFT(v,-d)))
          x(i) = v
       end do
    end if ! end of initialization

    ! Apart from initialization (above), this is the generator
    v = 0_i8
    Do While (v .eq. 0_i8)
       i = IAND(i+1,r-1)
       t = x(i)
       v = x(IAND(i+(r-s),r-1))
       t = IEOR(t,ISHFT(t,a))
       t = IEOR(t,ISHFT(t,-b))
       v = IEOR(v,ISHFT(v,c))
       v = IEOR(v,IEOR(t,ISHFT(v,-d)))
       x(i) = v
       w = w + weyl
       v = v + w
       v = ISHFT(v,-11)
    End Do

    DoubleRealRN = t53*v

    if( present(save_state) ) then
       save_state(1:r)   = x(0:r-1)
       save_state(r+1)   = i
       save_state(r+2)   = w
       if ((r+3) <= n_save_state) save_state(r+3:n_save_state) = 0
    end if

  end subroutine loopd

end module mo_test_loop

max_stars_repo_licenses: MIT Intel Unlicense NetCDF

max_issues_repo_licenses: MIT Intel Unlicense NetCDF

max_forks_repo_licenses: MIT Intel Unlicense NetCDF


text:
      MODULE m_vacudz
      CONTAINS
      SUBROUTINE vacudz(
     >     e,vz,vz0,nmz,dz,
     <     udz,dudz,ddnv,ue,
     >     duz,u)
c*********************************************************************
c     integrates the vacuum energy derivative wavefunction ue for the
c     energy e>0 using the schrodinger equation. the wavefunction u
c     is necessary (obtained from a previous call to vacuz). the
c     conventions are the same as in vacuz. ue will be orthogonal to
c     u. udz and dudz are the value and normal derivative at the
c     vacuum boundary.
c     based on code by m. weinert
c*********************************************************************
      use m_juDFT
      USE m_constants
      USE m_intgr, ONLY : intgz0
      IMPLICIT NONE
C     ..
C     .. Scalar Arguments ..
      INTEGER, INTENT (IN) :: nmz
      REAL,    INTENT (IN) :: duz,dz,e,vz0
      REAL,    INTENT (OUT):: ddnv,dudz,udz
C     ..
C     .. Array Arguments ..
      REAL, INTENT (IN) ::  u(nmz),vz(nmz)
      REAL, INTENT (OUT):: ue(nmz)
C     ..
C     .. Local Scalars ..
      REAL a,eru,erw,h,sk1,sk2,sk3,sk4,sl1,sl2,sl3,sl4,tnorm,u1,u10,
     +     u1p,vme,w1,w10,w1p,yn,zn,znm
      INTEGER i,it,n,n1
      LOGICAL tail
C     ..
C     .. Local Arrays ..
      REAL we(nmz),wp(nmz)
C     ..
C     .. Data statements ..
      REAL,PARAMETER:: eps=1.e-06
C     ..
      IF (e.GE.vz0)  CALL
     +     juDFT_error("e >= vz0",calledby ="vacudz",hint =
     +     'Vacuum energy-parameter too high')
c---  >    set up initial conditions
      znm = (nmz-1)*dz
      a = sqrt(2.* (vz0-e))
      ue(nmz) = (a*znm-0.5)*u(nmz)/ (a*a)
      we(nmz) = (a*znm-1.5)*u(nmz)/a
      wp(nmz) = 2.* (vz(nmz)-e)*ue(nmz) - 2.*u(nmz)
c---  >    use 4nd order runge-kutta to first few mesh points
      h = dz
      DO  i = 1,4
         n = nmz + 1 - i
         yn = ue(n)
         zn = we(n)
         sk1 = h*zn
         sl1 = h*wp(n)
         sk2 = h* (zn+0.5*sl1)
         sl2 = h* ((vz(n)+vz(n-1)-e)* (yn+0.5*sk1)- (u(n)+u(n-1)))
         sk3 = h* (zn+0.5*sl2)
         sl3 = h* ((vz(n)+vz(n-1)-e)* (yn+0.5*sk2)- (u(n)+u(n-1)))
         sk4 = h* (zn+sl3)
         sl4 = h* (2.* (vz(n-1)-e)* (yn+sk3)-2.*u(n-1))
         ue(n-1) = yn + (sk1+2.*sk2+2.*sk3+sk4)/6.
         we(n-1) = zn + (sl1+2.*sl2+2.*sl3+sl4)/6.
         wp(n-1) = 2.* (vz(n-1)-e)*ue(n-1) - 2.*u(n-1)
      ENDDO
c---  >    use adams-bashforth-moulton predictor-corrector for rest
      DO i = 5,nmz - 1
c---  >    adams-bashforth predictor (order h**5)
         n = nmz + 1 - i
         n1 = n - 1
         u10 = ue(n) + h* (55.*we(n)-59.*we(n+1)+37.*we(n+2)-
     +        9.*we(n+3))/24.
         w10 = we(n) + h* (55.*wp(n)-59.*wp(n+1)+37.*wp(n+2)-
     +        9.*wp(n+3))/24.
c---  >    evaluate derivative at next point( n-1 corresponds to m+1)
         vme = 2.* (vz(n-1)-e)
         u1p = w10
         w1p = vme*u10 - 2.*u(n-1)
c---  >    adams-moulton correctors: (order h**6)
         DO  it = 1,10
            u1 = ue(n) + h* (251.*u1p+646.*we(n)-264.*we(n+1)+
     +           106.*we(n+2)-19.*we(n+3))/720.
            w1 = we(n) + h* (251.*w1p+646.*wp(n)-264.*wp(n+1)+
     +           106.*wp(n+2)-19.*wp(n+3))/720.
c---  >    final evaluation
            u1p = w1
            w1p = vme*u1 - 2.*u(n-1)
c---  >    test quality of corrector and iterate if necessary
            eru = abs(u1-u10)/ (abs(u1)+abs(h*u1p))
            erw = abs(w1-w10)/ (abs(w1)+abs(h*w1p))
            IF (eru.LT.eps .AND. erw.LT.eps) EXIT
            u10 = u1
            w10 = w1
         ENDDO
         IF (it>10) THEN
            WRITE (oUnit,FMT=8000) n1,eru,erw
 8000       FORMAT (' ***vacudz - step may be too big - mesh point',i5,
     +           ', eru,erw=',1p,2e16.7)
         ENDIF
c---  >    store values
         ue(n1) = u1
         we(n1) = w1
         wp(n1) = w1p
      ENDDO
c---  >    ensure orthogonality
      tail = .true.
      DO  i = 1,nmz
         wp(i) = ue(nmz+1-i)*u(nmz+1-i)
      ENDDO
      CALL intgz0(wp,dz,nmz,tnorm,tail)
      ue(:) = ue(:) - tnorm*u(:)
      udz = ue(1)
      dudz = we(1) - tnorm*duz
c---  >    determine normalization
      DO  i = 1,nmz
         wp(i) = ue(nmz+1-i)*ue(nmz+1-i)
      enddo
      CALL intgz0(wp,dz,nmz,ddnv,tail)
      RETURN
      END subroutine
      END

ext: f

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
#include "scalar.fpp"

!=====================================================================!
! Module that defines an interface for SENSITIVITY ANALYSIS of
! physical systems. This interface is a 'decorator' on top of analysis
! interface.
!
! Main steps in sensitivty analysis are:
!
! (a) solving the physical system for state variables
! (b) solving the adjoint/direct sensitivity variables
!
! (i)  evaluating the functions of interest
! (ii) evaluating the derivatives of functions of interest
!
! Author: Komahan Boopathy (komahan@gatech.edu)
!=====================================================================!

module sensitivity_analysis_interface

  use analysis_interface, only : ianalysis => analysis
  use function_interface, only : ifunction => function

  implicit none

  private

  public :: sensitivity_analysis

  !-------------------------------------------------------------------!
  ! Type definition for analysis
  !-------------------------------------------------------------------!

  type, extends(analysis), abstract :: sensitivity_analysis

     ! Sentivity analysis needs the behavior of 'analysis' class
     class(ianalysis), allocatable :: analysis ! used to evaluate the states
     class(ifunction), allocatable :: function ! used to evaluate functions of interest

   contains

     ! Implement methods to evalaute the sensitivity variables

  end type sensitivity_analysis

  ! Interfaces to deferred procedures
  abstract interface
  end interface

contains

end module sensitivity_analysis_interface

max_stars_repo_name: komahanb/time-marching-adjoint

max_issues_repo_name: komahanb/diffusion

max_forks_repo_name: komahanb/diffusion


text:
! { dg-do compile }
! { dg-additional-options "-fopenmp" }
      SUBROUTINE ICHI
      INTEGER :: ARGC
      ARGC = COMMAND_ARGUMENT_COUNT ()

!$OMP PARALLEL
!$ACC PARALLEL                                                          &
!$ACC& COPYIN(ARGC) ! { dg-error "directive cannot be specified within" }
      IF (ARGC .NE. 0) THEN
         STOP 1
      END IF
!$ACC END PARALLEL
!$OMP END PARALLEL

      END SUBROUTINE ICHI

      SUBROUTINE NI
      IMPLICIT NONE
      INTEGER :: I

!$ACC PARALLEL                                                          &
!$OMP& DO ! { dg-error "Wrong OpenACC continuation" }
      DO I = 1, 10
      ENDDO

!$OMP PARALLEL                                                          &
!$ACC& LOOP ! { dg-error "Wrong OpenMP continuation" }
      DO I = 1, 10
      ENDDO
      END SUBROUTINE NI


SUBROUTINE p(k,x,val)

!-------------------------------------------------------------
!
!  Compute the values of P from the Legendre polynomials
!
!-------------------------------------------------------------

USE solvar
IMPLICIT NONE
INTEGER, INTENT(IN) :: k
REAL*8, INTENT(IN) :: x
REAL*8, INTENT(OUT) :: val

IF (k == 0) THEN
   val = 1.0
ELSE IF (k == 1) THEN
   val = x
ELSE IF (k == 2) THEN
   val = 1.5*(x**2) - 0.5
ELSE IF (k == 3) THEN
   val = 2.5*(x**3) - 1.5*x
ELSE IF (k == 4) THEN
   val = (35.0*(x**4) - 30.0*(x**2) + 3.0)/8.0
ELSE IF (k == 5) THEN
   val = (63.0*(x**5) - 70.0*(x**3) + 15.0*x)/8.0
ELSE IF (k == 6) THEN
   val = (231.0*(x**6) - 315.0*(x**4) + 105.0*(x**2) - 5.0)/16.0
ELSE IF (k == 7) THEN
   val = (429.0*(x**7) - 693.0*(x**5) + 315.0*(x**3) - 35.0*x)/16.0
ELSE IF (k == 8) THEN
   val = x*(((26.8125*x*x-43.3125)*x*x+19.6875)*x*x-2.1875)
ELSE
   warn = warn + 1
   WRITE (8,'(//,1X,A)') "WARNING: outside of range that Legendre is currently programmed in this p.f90 version"
   val = 0.0
END IF

RETURN
END SUBROUTINE p


!##############################################################################
! PROGRAM Arrays
!
! ## Use do-loops to fill up arrays with values
!
! This code is published under the GNU General Public License v3
!                         (https://www.gnu.org/licenses/gpl-3.0.en.html)
!
! Authors: Hans Fehr and Fabian Kindermann
!          contact@ce-fortran.com
!
! #VC# VERSION: 1.0  (23 January 2018)
!
!##############################################################################
program Arrays

    implicit none
    real*8 :: x(0:10), y(0:10)
    integer :: j

    ! initialize x and calculate y
    do j = 0, 10
        x(j) = 1d0/10d0*dble(j)
        y(j) = exp(x(j))
    enddo

    ! output table of values
    write(*,'(a)')'         X         Y'
    do j = 0, 10
        write(*,'(2f10.3)')x(j), y(j)
    enddo

end program


text:
http://en.wikipedia.org/wiki/Rosanna_(song) Rosanna is a song by the American rock band Toto, the opening track from its 1982 album Toto IV. This song won the Record of the Year Grammy Award in the 1983 presentations. Rosanna was also nominated for the Song of the Year award. In musician circles, the song is known for its highly influential halftime shuffle, as well as a blowing ending guitar solo played by guitarist Steve Lukather. The music video features Cynthia Rhodes as the lead dancer and Patrick Swayze as one of the male dancers. They would reunite as members of the cast of the 1987 movie Dirty Dancing. —Wikipedia

Also a http://daviswiki.org/Design UCD Design / http://daviswiki.org/KDVS KDVS 90.3 FM / http://daviswiki.org/Creative_Media Creative Media alumni who sporadically reminisces about the days of yore.

http://rosannayau.com rosannayau.com

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


!***********************************/
!	Name:計算設定を読み込むプログラム
!	Alias:ReadConfigulation
!	Description:専用のCalcConfigファイルが必要
!	Type:Configulation
!	Input:CalcConfig(外部入力)
!	Output:Configulation
!	Note:CalcConfigはプログラム本体と同じディレクトリに置くこと
!	Author:Akitaka Toyota
!	Date:2017.10.26
!	Update:2017.11.07
!	Other:
!***********************************/
subroutine ReadConfigulation(Conf, my_rank)
    use StructVar_Mod
    use ConstantVar_Mod
    implicit none
    type(Configulation), intent(inout) :: Conf
    integer :: my_rank
    character :: cAnnotate
    integer :: tolerance_exp
    logical :: debug = .false.

    call make_logfile_name(Conf%cLogName)

    open(unit=my_rank+100, file='CalcConfig',status='unknown')
        read(my_rank+100,*) cAnnotate
        read(my_rank+100,*) Conf%SwitchProgram, cAnnotate
        read(my_rank+100,*) Conf%UseReadRegion,cAnnotate
        read(my_rank+100,*) Conf%UseResume,cAnnotate
        read(my_rank+100,*) Conf%ResumeFileName
        read(my_rank+100,*) Conf%ResumeFileFormat
        read(my_rank+100,*) Conf%ResumeNumber
        read(my_rank+100,*) Conf%cCaseName
        read(my_rank+100,*) Conf%UseReadBoundaryCond,cAnnotate
        read(my_rank+100,*) Conf%UseVariableTime,cAnnotate
        read(my_rank+100,*) Conf%UseSteadyCalc,cAnnotate
        read(my_rank+100,*) Conf%UseLocalTimeStep,cAnnotate
        read(my_rank+100,*) Conf%UseRRK2,cAnnotate
        read(my_rank+100,*) Conf%UseMUSCL,cAnnotate
        read(my_rank+100,*) Conf%KindLimiter,cAnnotate
        read(my_rank+100,*) Conf%UseOverSet, cAnnotate
        read(my_rank+100,*) FixedTimeStep,cAnnotate
        read(my_rank+100,*) OutputInterval, cAnnotate
        read(my_rank+100,*) IterationNumber,cAnnotate
        read(my_rank+100,*) GridNumber, cAnnotate
        read(my_rank+100,*) SigmoidGain, cAnnotate
        read(my_rank+100,*) EntropyCorrection, cAnnotate
        read(my_rank+100,*) VenkatakrishnanParameterK, cAnnotate
        read(my_rank+100,*) CourantFriedrichsLewyCondition, cAnnotate
        read(my_rank+100,*) Conf%TurbulenceModel, cAnnotate
        read(my_rank+100,*) ReynoldsNumber, cAnnotate
        read(my_rank+100,*) Conf%UseFluxMethod, cAnnotate
        read(my_rank+100,*) Converge_Method
        read(my_rank+100,*) Converge_tolerance, cAnnotate
    close(my_rank+100)

    if(trim(adjustl(Conf%cCaseName)) == "null") Conf%cCaseName = ""
    if(Conf%TurbulenceModel /= 0) invicid = .false.
    if(Conf%UseRRK2 == 1) IterationNumber = 2*IterationNumber
    DefaultTimeStep = FixedTimeStep

    if(Conf%UseJobParallel == 0) then
        Conf%my_rank = 0
    end if

    call Show_Configulation(debug)

    return

contains
    subroutine Show_Configulation(debug)
        logical, intent(in) :: debug

        if (debug == .false.) then
            ! no action
        else
            write(6,*) "show_configulation == TRUE"
            write(6,*) "Conf%SwitchProgram", Conf%SwitchProgram
            write(6,*) "Conf%UseReadRegion", Conf%UseReadRegion
            write(6,*) "Conf%UseResume", Conf%UseResume
            write(6,*) "Conf%ResumeFileName", Conf%ResumeFileName
            write(6,*) "Conf%ResumeFileFormat", Conf%ResumeFileFormat
            write(6,*) "Conf%ResumeNumber", Conf%ResumeNumber
            write(6,*) "Conf%UseReadBoundaryCond", Conf%UseReadBoundaryCond
            write(6,*) "Conf%UseVariableTime", Conf%UseVariableTime
            write(6,*) "Conf%UseSteadyCalc", Conf%UseSteadyCalc
            write(6,*) "Conf%UseLocalTimeStep", Conf%UseLocalTimeStep
            write(6,*) "Conf%UseRRK2", Conf%UseRRK2
            write(6,*) "Conf%UseMUSCL", Conf%UseMUSCL
            write(6,*) "Conf%KindLimiter", Conf%KindLimiter
            write(6,*) "Conf%UseOverSet", Conf%UseOverSet
            write(6,*) "FixedTimeStep", FixedTimeStep
            write(6,*) "OutputInterval", OutputInterval
            write(6,*) "IterationNumber", IterationNumber
            write(6,*) "GridNumber", GridNumber
            write(6,*) "SigmoidGain", SigmoidGain
            write(6,*) "EntropyCorrection", EntropyCorrection
            write(6,*) "VenkatakrishnanParameterK", VenkatakrishnanParameterK
            write(6,*) "CourantFriedrichsLewyCondition", CourantFriedrichsLewyCondition
            write(6,*) "Conf%TurbulenceModel", Conf%TurbulenceModel
            write(6,*) "ReynoldsNumber", ReynoldsNumber
            write(6,*) "UseFluxMethod", Conf%UseFluxMethod
        end if

        return
    end subroutine Show_Configulation

    subroutine make_logfile_name(fname)
        character(len=256), intent(inout) :: fname
        integer :: datetime(8)
        character(len=10) :: sys_time(3)

        call date_and_time(sys_time(1), sys_time(2), sys_time(3), datetime)

        fname = "CalcLog_"//trim(adjustl(sys_time(1)))//"_"//trim(adjustl(sys_time(2)))//".txt"
        return
    end subroutine make_logfile_name

end subroutine ReadConfigulation


text:
When the urge to Barbeque hits and you dont own your own propane fired grill it is time to go looking for a public grill.

 City Parks
  Chestnut Park
  Civic Center Park
  La Playa Park
  Mace Ranch Park Covered Picnic Area
  N Street Park
  Northstar Pocket Park
  Oxford Circle Park
  Putah Creek Park
  Redwood Park
  Sandy Motley Park
  Westwood Park
  Whaleback Park

 Campus
  Putah Creek Lodge
  Rec Pool Lodge

 Apartment Complexes
  La Salle
  Temescal Apartments


text:
program biot
!Pierre Navaro
!navaro@unistra.fr
!IRMA http://www-irma.u-strasbg.fr

use sorties, only: gnuplot_output
use initial, only: lecture

implicit none

integer, parameter :: idm = 100000
integer :: nbpart, i, j, k, n, istep, imov, nstep
real :: xp(idm), yp(idm), op(idm), up(idm), vp(idm)
real, dimension(64,64) :: w
real :: time

real(8) :: t0, t1, tcpu

real :: delta
real :: dt

call cpu_time(tcpu)
t0 = getRealTimer()
call lecture( nstep, imov, xp, yp, op, delta, idm, dt, nbpart )
call calcul_w(xp, yp, op, w, nbpart)

open(10, file="initial_w.dat")
do i = 1, 64
do j = 1, 64
   write(10,*) i, j, w(i,j)
end do
write(10,*)
end do
close(10)

time = 0.0

do istep = 1, nstep       !loop over time

   call vitesse(nbpart, xp, yp, op, up, vp, delta)
   call deplace(nbpart, xp, yp, up, vp, dt)
!   call calcul_w(xp, yp, op, w, nbpart)

   time = time + dt

   if ( mod(istep, imov) == 0) &
   call gnuplot_output(istep, time, xp, yp, op, nbpart, nstep)
   write(*,"(i5,1x,1g12.3)")istep, time

end do      !next time step

call cpu_time(tcpu)
t1 = getRealTimer()
write(*,"(5x,' CPU time = ', G15.3)") tcpu
write(*,"(5x,' CALLSITE time = ', G15.3)") t1-t0

contains

subroutine vitesse (nbpart, xp, yp, op, up, vp, delta)
implicit none
integer, intent(in)  :: nbpart
real, intent(in)     :: xp(nbpart), yp(nbpart), op(nbpart)
real, intent(inout)  :: up(nbpart), vp(nbpart)
real, intent(in)     :: delta
real :: xo, yo, dx, dy, usum, vsum, dpi, d2
integer :: j, k

real(8) :: r2, r22, r2a1, r2a13, xm
real(8) :: a1, a12, a122, fi,fj

dpi   = 8.0 * atan( 1.0 )

a1    = 0.01
a12   = a1*a1
a122  = a12*a12

up = 0.0
vp = 0.0

do k = 1, nbpart-1

   xo = xp(k); yo = yp(k)

   do j = k+1 , nbpart
      dx    = xp( j ) - xo
      dy    = yp( j ) - yo
      r2    = dx * dx + dy * dy
      r22   = r2 * r2
      r2a1  = r2 + a12
      r2a13 = r2a1 * r2a1 * r2a1
      xm    = (r22+3.0*a12*r2+4.0*a122) / r2a13
      fi    = dy *  xm
      fj    = - dx * xm
      up(k)  = up(k) + fi * op(j)
      vp(k)  = vp(k) + fj * op(j)
      up(j)  = up(j) - fi * op(k)
      vp(j)  = vp(j) - fj * op(k)
   end do

end do

up = up / dpi
vp = vp / dpi

end subroutine vitesse

subroutine deplace (nbpart, xp, yp, up, vp, dt)
implicit none
integer, intent(in)  :: nbpart
real, intent(inout)  :: xp(nbpart), yp(nbpart)
real, intent(in)     :: up(nbpart), vp(nbpart)
real, intent(in)     :: dt
integer :: k

do k = 1, nbpart
   xp(k) = xp(k) + dt * up(k)
   yp(k) = yp(k) + dt * vp(k)
end do

end subroutine deplace

subroutine calcul_w(xp, yp, op, w, nbpart)
implicit none
integer, intent(in) :: nbpart
real, dimension(nbpart), intent(in) :: xp, yp, op
integer :: ip, jp, kp, i, j
real, parameter :: xmin    = - 2., xmax = 2.
real, parameter :: ymin    = - 2., ymax = 2.
real :: a1, a2, a3, a4, xt, yt
logical :: lflag
integer, parameter :: nxw=64, nyw=64
real, dimension(nxw,nyw), intent(inout) :: w

lflag   = .false.

do j = 1, nyw
do i = 1, nxw
w(i,j) = 0.0
end do
end do

do kp = 1,nbpart

   xt = (xp(kp)-xmin)/(xmax-xmin)*nxw
   yt = (yp(kp)-ymin)/(ymax-ymin)*nyw

   ip = floor(xt)
   jp = floor(yt)

   if (      (ip <= 0 .or. ip >= nxw) &
        .or. (jp <= 0 .or. jp >= nyw)) then

      lflag = .true.

   else

      a1 = (ip+1 - xt) * (jp+1 - yt)
      a2 = (xt   - ip) * (jp+1 - yt)
      a3 = (ip+1 - xt) * (yt   - jp)
      a4 = (xt   - ip) * (yt   - jp)

      w(ip  ,jp  )=w(ip  ,jp  )+a1*op(kp)
      w(ip+1,jp  )=w(ip+1,jp  )+a2*op(kp)
      w(ip  ,jp+1)=w(ip  ,jp+1)+a3*op(kp)
      w(ip+1,jp+1)=w(ip+1,jp+1)+a4*op(kp)

   end if

end do

!if (lflag) then
!
!   write(*,*) " xmin, xmax = ", xmin, xmax
!   write(*,*) " ymin, ymax = ", ymin, ymax
!   write(*,*) " xpmin, xpmax = ", minval(xp(:)), maxval(yp(:))
!   write(*,*) " ypmin, ypmax = ", minval(xp(:)), maxval(yp(:))
!   stop 'Calcul de la w. Revoir les parametres'
!
!end if

end subroutine calcul_w

function getRealTimer()
implicit none
real(8) :: out, getRealTimer
integer(8) :: count, count_rate, max
call system_clock(count, count_rate)
count = count - 1254348000*count_rate
out = count
out = out / count_rate
getRealTimer = out
end function getRealTimer

end program biot

max_stars_repo_name: pnavaro/vortex-pp

max_issues_repo_name: pnavaro/vortex-pp

max_forks_repo_name: pnavaro/vortex-pp


text:
Buildings Breeding is a Music Scene local band that debuted at Delta of Venus in October 2006. They are now based in Los Angeles. They are signed with http://www.devilinthewoodsrecords.com/index_2.htm Devil In The Woods Records. The group is composed of:

  Users/ChristapherLarsen Christapher Larsen  vocals, guitars, bass, beer
  Users/MelanieGlover Melanie Glover  drums, vocals, estrogen
  Kevin Dockter  vocals, guitar
  Chris Vogel  bass, designated driver
  Justin Titsworth  drums and percussion

Former members:
  Users/EvanHart Evan Hart  lead guitar, atmosphere, luscious locks
  Users/JesseMiller Jesse Miller  synths
  Paul Bae  bass
  Nick lastnamenot available  bass

Their song Stacking Up Reasons, was the second track of the 2010 fundraiser collection Audioversity in the Little Tomato, and can be heard http://soundcloud.com/davisca/02buildingsbreedingstackingupreasons here.

Releases and Photos