c------------------------------------------------------------------------------
c The CDM variance in spheres containing mass m.
c This fits the variance derived from the BBKS CDM power spectrum.
c The units of mass are Msol/h chosen to scale as the mass within a sphere
c of radius r Mpc/h. The scaling of the amplitude with h is such that the
c rms for a sphere of 8 Mpc/h is given by sigma_8 for any h and omega
c In these units if h is varied the characteristic radius goes like
c 1/(omega h) and the characteristic mass like 1/(h^3 omega^2)
c This version also returns the logarithmic slope, alpha, of sigma(m).
c
c If nspec<10 then a powerlaw P(k)=k^nspec is assumed instead of CDM
c
c
      real function sigmacdm(m,alpha)
*************************************variables*********************************
      implicit none
      integer ifirst,igwave,itrans
      real m,ms,a1,a2,a3,p,scla,sclm,m8,sigma,
     x     alpha,fac1,fac2,curl,curlp
      save scla,sclm,ifirst
      parameter( a1=8.6594e-12, a2=3.5 ,a3=1.628e+09, p=0.255) 
      real omega0,lambda0,h0,omegab,nspec,gamma,sigma8
      common /cosmology/ omega0,lambda0,h0,omegab
      common /power_spectrum/ nspec,gamma,sigma8,igwave,itrans
      data ifirst/1/
*******************************************************************************
      if (ifirst.eq.1) then
c       Compute the required scaling factors sclm and scla
        sclm = gamma**3 / omega0
        m8=6.005e+14*omega0      !The mass within an 8Mpc/h sphere
        ms=m8*sclm
        if (nspec.gt.10) then
c         The fit to CDM for Gamma=1
           fac1=a2*ms**-0.1
           fac2=a3*ms**-0.63
           curl=(fac1+fac2)
           curlp=curl**p
           sigma = 1.0 / (a1*ms**0.67 * (1.0+curlp )**(1.0/p) )
        else
           sigma = ms**(-(nspec+3.0)/6.0)
        end if
        scla=sigma8/sigma      !scales sigma_8 to required value
c
        ifirst=0   !indicates first call complete and sclm and scla are set
      end if
c        ----------------------------------------------------
      ms=m*sclm
      if (nspec.gt.10) then
c        The fit to CDM for Gamma=1
         fac1=a2*ms**-0.1
         fac2=a3*ms**-0.63
         curl=(fac1+fac2)
         curlp=curl**p
         sigmacdm=scla  / (a1*ms**0.67 * (1.0+curlp )**(1.0/p) )
         alpha=-0.67 + (0.1*fac1+0.63*fac2)/(curl*(1.0/curlp+1.0))
      else
         alpha= -(nspec+3.0)/6.0
         sigma = scla * ms**alpha
      end if
      return
      end