double precision function integrand(lnM,z)
  USE battaglia
  USE cosmo
  USE sigma
  USE growth
  USE angular_distance
  USE nfw
! The mean Compton Y is given by <Y> = int dz dV/dz int dlnM dn/dlnM*Y.
! This routine returns the integrand, dV/dz*dn/dlnM*Y. The pressure profile
! is taken from Battaglia et al., ApJ, 758, 75 (2012).
! August 2, 2015: E.Komatsu
  IMPLICIT none
  double precision, intent(IN) :: lnM,z
  double precision :: xin=1d-5,xout=2.5d0 ! pressure profile is cut at r=xout*rvir
  double precision :: omega,tol=1d-6,pi=3.1415926535d0
  double precision :: zin=1d-5 ! redshift of the input matter power spectrum
  double precision :: mvir,rvir,rs,delc,r200,l200,P200
  double precision :: rhoc,Ez,da,dvdz,transform2d,tsz,pgnfw,rombint
  double precision :: deltac=1.6865d0,mf,dndlnRh,dndlnMh,lnnu,dlnnudlnRh,Rh
  double precision :: g,scalefactor,m180d,m200d,fb
  real :: chebev,lnsigma2,lnRh,dlnsigma2dlnRh
  integer :: i
  external pgnfw,rombint
  mvir=dexp(lnM)
! compute omega's, rhoc, and delc
  omega=om0*(1d0+z)**3d0/Ez(z)**2d0 ! Omega(z); E(z)=H(z)/H0
  rhoc=2.775d11*Ez(z)**2d0 ! critical density in units of h^2 M_sun/Mpc^3
  ! Eq.(6) of Bryan & Norman (1998), valid only for a flat universe!
  delc=18d0*pi*pi+82d0*(omega-1d0)-39d0*(omega-1d0)**2d0
  ! One may also use Eq.(C19) of Nakamura & Suto (1997):
  !!$ delc=omega*18d0*pi*pi*(1d0+0.4093d0*(1d0/omega-1d0)**0.9052d0)
! compute virial and NFW parameters
  rvir=(3d0*mvir/4d0/pi/delc/rhoc)**(1d0/3d0) ! virial radius, h^-1 Mpc
  ! Concentration parameter from Duffy et al. (2008)
  cvir=7.85d0*(mvir/2d12)**(-0.081)/(1d0+z)**0.71
  ! One may also use the concentration parameter from Seljak (2000) 
  !!$ cvir=10d0*(mvir/3.42d12)**(-0.2)/(1d0+z)
  rs=rvir/cvir ! NFW scale radius, h^-1 Mpc
! convert to m200 for Battaglia et al.'s profile
  zz=z
  CALL mvir2mdel(mvir,rs,cvir,200d0*rhoc,m200)
  r200=(3d0*m200/4d0/pi/200d0/rhoc)**(1d0/3d0) ! h^-1 Mpc
  l200=da(z)/r200
! calculate the Fourier transform of Tsz
  ! See, e.g., Appendix D1 of Komatsu et al., arXiv:1001.4538
  !            Eq.(2) of Komatsu&Seljak (2002)
  fb=obh2/(om0*h0**2d0)
  P200=0.50439d0*200d0*(h0*Ez(z))**2d0*fb*(m200/1d15)/(2d0*r200) ! eV/cm^3
  transform2d=P200*rombint(pgnfw,xin,xout*rvir/r200,tol)*(4d0*3.14159d0)/l200**2d0*(r200/h0)
  Tsz=283d0*(transform2d/50d0)/2.725d6 ! compton y
  !!$Tsz=-2d0*283d0*(transform2d/50d0) ! uK, Rayleigh-Jeans limit
  !!$Tsz=-0.953652*283d0*(transform2d/50d0) ! uK, 150GHz
! calculate mass function
  ! Sheth&Tormen, Tinker et al., and Bocquet et al.'s mass functions are given for the 
  ! overdensity mass M200d (with respect to the mean mass density rather than 
  ! the critical density).
  CALL mvir2mdel(mvir,rs,cvir,200d0*omega*rhoc,m200d)
  Rh=(3d0*m200d/4d0/pi/om0/2.775d11)**(1d0/3d0) ! h^-1 Mpc
  ! Alternatively, one may wish to use Jenkins et al.'s mass function,
  ! which is given for the overdensity mass M180d (with respect to the mean 
  ! mass density rather than the critical density).
  !!$ CALL mvir2mdel(mvir,rs,cvir,180d0*omega*rhoc,m180d)
  !!$ Rh=(3d0*m180d/4d0/pi/om0/2.775d11)**(1d0/3d0) ! h^-1 Mpc
  lnRh=real(dlog(Rh))
  lnsigma2=CHEBEV(lnR1,lnR2,c,ndim,lnRh)          ! ln(sigma^2)
  dlnsigma2dlnRh=CHEBEV(lnR1,lnR2,cder,ndim,lnRh) ! dln(sigma^2)/dlnRh
  scalefactor=g(z)/g(zin)*(1d0+zin)/(1d0+z)
  lnnu=2d0*dlog(deltac)-dble(lnsigma2)-2d0*dlog(scalefactor) ! ln(nu)
  dlnnudlnRh=-dble(dlnsigma2dlnRh)     ! dln(nu)/dlnRh
  dndlnRh=(3d0/4d0/pi)*dlnnudlnRh*mf(lnnu,z)/Rh**3d0
  !!$ dndlnRh=(3d0/4d0/pi)*dlnnudlnRh*mf(lnnu)/Rh**3d0 ! if using ST or Jenkins
  dndlnMh=dndlnRh/3d0 ! in units of h^3 Mpc^-3
! volume element
  dvdz=(1d0+z)**2d0*da(z)**2d0*2998d0/Ez(z) ! h^-3 Mpc^3
! This is the integrand: dV/dz*dn/dlnM*Tsz, to be integrated over lnM and z
  integrand=dvdz*dndlnMh*Tsz
  return
end function integrand
!-----------------------------------------------------------------------------
DOUBLE PRECISION FUNCTION Ez(z)
  USE cosmo
  ! Ez = H(z)/H0 (dimensionless)
  ! x=a
  IMPLICIT none
  DOUBLE PRECISION, intent(IN) :: z
  double precision :: x,ok0,or0=0d0 ! radiation density is ignored!
  ok0=1d0-om0-ode0
  x=1d0/(1d0+z)
  Ez=dsqrt(om0/x**3d0+ok0/x**2d0+or0/x**4d0+ode0/x**(3d0+3d0*w))
  return
END FUNCTION Ez
!-----------------------------------------------------------------------------
DOUBLE PRECISION FUNCTION pgnfw(x) ! x=r/r200
! Ref: Battaglia et al., ApJ, 758, 75 (2012)
  USE cosmo
  USE battaglia
  IMPLICIT none
  double precision :: x
  double precision :: xc,g=-0.3d0,a=1d0,b,P0
  xc=0.497d0*(m200/h0/1d14)**(-0.00865d0)*(1d0+zz)**(0.731d0)
  P0=18.1d0*(m200/h0/1d14)**(0.154d0)*(1d0+zz)**(-0.758d0)
  b=4.35d0*(m200/h0/1d14)**(0.0393d0)*(1d0+zz)**(0.415d0)
  pgnfw=P0*(x/xc)**g/(1d0+(x/xc)**a)**b*x**2d0
  return
END FUNCTION pgnfw