I think the callable function at line 513 should be : pspec = lambda k: h**3. * ccl.linear_matter_power(self.cosmo, k.flatten() * h, a=1.).reshape(k.shape) because my pyccl version generates 1D power spectra, while array k is 3D.
However, the reshape at line 565,570,575 fails, maybe it is not needed here?