double precision function integrand(lnM,z)
  USE cosmo
  USE sigma
  USE growth
  USE angular_distance
  USE nfw
  USE multipole
! The SZ power spectrum is given by C_l = int dz dV/dz int dlnM dn/dlnM*Tsz^2.
! This routine returns the integrand, dV/dz*dn/dlnM*Tsz^2. Note that Tsz is 
! in units of uK, Rayleigh-Jeans limit.
! December 13, 2010: E.Komatsu
  IMPLICIT none
  double precision, intent(IN) :: lnM,z
  double precision :: xin=1d-5,xout=3d0 ! pressure profile is cut at r=xout*rvir
  double precision :: omega,tol=1d-6,zin=30d0,pi=3.1415926535d0
  double precision :: mvir,rvir,ls,rs,delc,eta0,rhogas0,kTgas0
  double precision :: rhoc,Ez,da,dvdz,transform2d,tsz,ygas,ypgas,rombint1
  double precision :: deltac=1.6865d0,mf,dndlnRh,dndlnMh,lnnu,dlnnudlnRh,Rh
  double precision :: g,scalefactor,m180d,m200d
  real :: chebev,lnsigma2,lnRh,dlnsigma2dlnRh
  integer :: i
  external ygas,ypgas,rombint1
  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
  ls=da(z)/rs
! calculate the Fourier transform of Tsz
  rhogas0=7.96d13*(obh2/om0)*(mvir/1d15)/rvir**3d0 &
         *cvir**2d0/(1d0+cvir)**2d0/ygas(cvir)/(dlog(1d0+cvir)-cvir/(1d0+cvir))
  eta0=2.235d0+0.202d0*(cvir-5d0)-1.16d-3*(cvir-5d0)**2d0
  kTgas0=8.8d0*eta0*(mvir/1d15)/rvir
  ! See, e.g., Appendix D2 of Komatsu et al., arXiv:1001.4538
  !            Eq.(2) of Komatsu&Seljak (2002)
  transform2d=55d0*(rhogas0/1d14)*(kTgas0/8d0)*(rs/h0) &
       *rombint1(ypgas,cvir*xin,cvir*xout,tol,ell/ls)*(4d0*pi)/ls**2d0
  Tsz=-2d0*283d0*(transform2d/50d0) ! uK, Rayleigh-Jeans limit
! calculate mass function
  ! Sheth&Tormen and Tinker 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^2, to be integrated over lnM and z
  integrand=dvdz*dndlnMh*Tsz**2d0
  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 ygas(x) ! r/rs
  USE nfw
  IMPLICIT none
  double precision :: x
  double precision :: gam,eta0,B
  gam=1.137d0+8.94d-2*dlog(cvir/5d0)-3.68d-3*(cvir-5d0)
  eta0=2.235d0+0.202d0*(cvir-5d0)-1.16d-3*(cvir-5d0)**2d0
  B=3d0/eta0*(gam-1d0)/gam/(dlog(1d0+cvir)/cvir-1d0/(1d0+cvir))
  ygas=(1d0-B*(1d0-dlog(1d0+x)/x))**(1d0/(gam-1d0))
  return
END FUNCTION ygas
!-----------------------------------------------------------------------------
DOUBLE PRECISION FUNCTION ypgas(x,y) ! x=r/rs, y=l/ls, ypgas=ygas^gamma
  USE nfw
  IMPLICIT none
  double precision :: x,y
  double precision :: gam,eta0,B
  gam=1.137d0+8.94d-2*dlog(cvir/5d0)-3.68d-3*(cvir-5d0)
  eta0=2.235d0+0.202d0*(cvir-5d0)-1.16d-3*(cvir-5d0)**2d0
  B=3d0/eta0*(gam-1d0)/gam/(dlog(1d0+cvir)/cvir-1d0/(1d0+cvir))
  ypgas=(1d0-B*(1d0-dlog(1d0+x)/x))**(gam/(gam-1d0)) &
       *x**2d0*dsin(y*x)/(y*x)
  return
END FUNCTION ypgas