def METHOD_NAME(cp,karr,z,integrate=True,fname_out=None) :
    cosmo=ccl.Cosmology(Omega_c=cp['Om_m']-cp['Om_b'],Omega_b=cp['Om_b'],h=cp['h'],
                        sigma8=cp['sig8'],n_s=cp['n'],
                        transfer_function='eisenstein_hu',matter_power_spectrum='linear')
    lmarr=LMMIN+(LMMAX-LMMIN)*(np.arange(NM)+0.5)/NM #log(M*h/M_sum)
    dlm=(LMMAX-LMMIN)/NM
    lmarrb=np.zeros(NM+2); lmarrb[0]=lmarr[0]-dlm; lmarrb[1:-1]=lmarr; lmarrb[-1]=lmarr[-1]+dlm
    marr=10.**lmarr #M (in Msun/h) M_h Ms/h = M Ms
    marrb=10.**lmarrb #M (in Msun/h) M_h Ms/h = M Ms
    sigmarrb=ccl.sigmaM(cosmo,marrb/cp['h'],1./(1+z))
    sigmarr=sigmarrb[1:-1]
    dlsMdl10M=np.fabs((sigmarrb[2:]-sigmarrb[:-2])/(2*dlm))/sigmarr
    omega_z=cp['Om_m']*(1+z)**3/(cp['Om_m']*(1+z)**3+1-cp['Om_m'])
    delta_c=params['DELTA_C']*(1+0.012299*np.log10(omega_z))
    Delta_v=(18*np.pi**2+82*(omega_z-1)-39*(omega_z-1)**2)/omega_z
    fM=st_gs(sigmarr,delta_c)*dlsMdl10M
    bM=st_b1(sigmarr,delta_c)
    cM=duffy_c(marr,z)
    cMf=interp1d(lmarr,cM,kind='linear',bounds_error=False,fill_value=cM[0])
    rhoM=ccl.rho_x(cosmo,1.,'matter')/cp['h']**2
    #Compute mass function normalization
    fM0=1-np.sum(fM)*dlm
    fbM0=1-np.sum(fM*bM)*dlm
    rvM=(3*marr/(4*np.pi*rhoM*Delta_v))**0.333333333333
    rsM=rvM/cM
    rsM0=(3*10.**LM0/(4*np.pi*rhoM*Delta_v))**0.333333333/duffy_c(10.**LM0,z)
    rsMf=interp1d(lmarr,rsM,kind='linear',bounds_error=False,fill_value=rsM0)
    f1hf=interp1d(lmarr,marr*fM/rhoM,kind='linear',bounds_error=False,fill_value=0)
    f1h0=fM0*10.**LM0/rhoM
    f2hf=interp1d(lmarr,fM*bM,kind='linear',bounds_error=False,fill_value=0)
    f2h0=fbM0
    #NFW profile
    def u_nfw(lm,k) :
        x=k*rsMf(lm) #k*r_s
        c=cMf(lm)
        sic,cic=sici(x*(1+c))
        six,cix=sici(x)
        fsin=np.sin(x)*(sic-six)
        fcos=np.cos(x)*(cic-cix)
        fsic=np.sin(c*x)/(x*(1+c))
        fcon=np.log(1.+c)-c/(1+c)
        return (fsin+fcos-fsic)/fcon
    def p1h(k) :
        def integ(lm) :
            u=u_nfw(lm,k)
            f1h=f1hf(lm)
            return u*u*f1h
        if integrate :
            return quad(integ,lmarr[0],lmarr[-1])[0]+f1h0*(u_nfw(LM0,k))**2
        else :
            return np.sum(integ(lmarr))*dlm+f1h0*(u_nfw(LM0,k))**2
    def b2h(k) :
        def integ(lm) :
            u=u_nfw(lm,k)
            f2h=f2hf(lm)
            return u*f2h
        if integrate :
            return quad(integ,lmarr[0],lmarr[-1])[0]+f2h0*u_nfw(LM0,k)
        else :
            return np.sum(integ(lmarr))*dlm+f2h0*u_nfw(LM0,k)
    p1harr=np.array([p1h(kk) for kk in karr])
    b2harr=np.array([b2h(kk) for kk in karr])
    pklin=ccl.linear_matter_power(cosmo,karr*cp['h'],1./(1+z))*cp['h']**3
    p2harr=pklin*b2harr**2
    ptarr=p1harr+p2harr
    np.savetxt(fname_out,np.transpose([karr,pklin,p2harr,p1harr,ptarr]))
    return pklin,pklin*b2harr**2,p1harr,pklin*b2harr**2+p1harr