diff -r -c -b -B -N CAMB/inidriver.F90 CAMB_addon/inidriver.F90
*** CAMB/inidriver.F90	2005-03-31 00:37:13.000000000 +0100
--- CAMB_addon/inidriver.F90	2006-02-20 20:32:38.000000000 +0000
***************
*** 171,193 ****
               stop 'Unknown recombination'
             end if
  
             P%InitPower%nn = Ini_Read_Int('initial_power_num')
             if (P%InitPower%nn>nnmax) stop 'Too many initial power spectra - increase nnmax in InitialPower'
             P%InitPower%rat(:) = 1
             do i=1, P%InitPower%nn
                write (numstr,*) i 
                numstr=adjustl(numstr)
                P%InitPower%an(i) = &
                     Ini_Read_Double('scalar_spectral_index('//trim(numstr)//')')
  
                P%InitPower%n_run(i) = &
                     Ini_Read_Double('scalar_nrun('//trim(numstr)//')',0._dl)
      
                if (P%WantTensors) then
                   P%InitPower%ant(i) = Ini_Read_Double('tensor_spectral_index('//trim(numstr)//')')
                   P%InitPower%rat(i) = Ini_Read_Double('initial_ratio('//trim(numstr)//')')
                end if              
! 
                P%InitPower%ScalarPowerAmp(i) = Ini_Read_Double('scalar_amp('//trim(numstr)//')',1.d0) 
                !Always need this as may want to set tensor amplitude even if scalars not computed
             end do
--- 171,214 ----
               stop 'Unknown recombination'
             end if
  
+ !addon
+             P%InitPower%WantExternalScalar = Ini_Read_Logical('use_external_scalar',.false.)
+             P%InitPower%WantExternalTensor = Ini_Read_Logical('use_external_tensor',.false.)
+ !end addon
+ 
             P%InitPower%nn = Ini_Read_Int('initial_power_num')
             if (P%InitPower%nn>nnmax) stop 'Too many initial power spectra - increase nnmax in InitialPower'
             P%InitPower%rat(:) = 1
             do i=1, P%InitPower%nn
                write (numstr,*) i 
                numstr=adjustl(numstr)
+ !addon
+               if (P%WantScalars .or. P%WantTransfer) then
+                  if (P%InitPower%WantExternalScalar) then
+                     P%InitPower%ScalarPowerFileName(i) & 
+                          = Ini_Read_String('scalar_input_file('//trim(numstr)//')')
+                  else
+ 
                      P%InitPower%an(i) = &
                           Ini_Read_Double('scalar_spectral_index('//trim(numstr)//')')
  
                      P%InitPower%n_run(i) = &
                           Ini_Read_Double('scalar_nrun('//trim(numstr)//')',0._dl)
  
+                  endif
+               endif
+ 
                if (P%WantTensors) then
+ 
+                  if (P%InitPower%WantExternalTensor) then
+                     P%InitPower%TensorPowerFileName(i) & 
+                          = Ini_Read_String('tensor_input_file('//trim(numstr)//')')
+                  else
                      P%InitPower%ant(i) = Ini_Read_Double('tensor_spectral_index('//trim(numstr)//')')
                      P%InitPower%rat(i) = Ini_Read_Double('initial_ratio('//trim(numstr)//')')
+                  endif
                end if              
! !end addon
                P%InitPower%ScalarPowerAmp(i) = Ini_Read_Double('scalar_amp('//trim(numstr)//')',1.d0) 
                !Always need this as may want to set tensor amplitude even if scalars not computed
             end do
***************
*** 247,252 ****
--- 268,282 ----
         AccuracyBoost  = Ini_Read_Double('accuracy_boost',1.d0)
         lAccuracyBoost = Ini_Read_Double('l_accuracy_boost',1.d0)
         lSampleBoost   = Ini_Read_Double('l_sample_boost',1.d0)
+ !addon     
+        DoClsInterpol  = Ini_Read_Logical('do_cls_interpol',.false.)
+        ComputeEveryl  = Ini_Read_Logical('compute_every_l',.false.)
+   
+        if (FeedbackLevel > 0) then
+           write (*,*)'do_cls_interpol is',DoClsInterpol
+           write (*,*)'compute_every_l is',ComputeEveryl          
+        endif
+ !end addon
  
         if (outroot /= '') then
           call Ini_SaveReadValues(trim(outroot) //'params.ini',1)
***************
*** 290,296 ****
  #endif
          end if
  
!              
    
          end program driver
  
--- 320,329 ----
  #endif
          end if
  
! !addon
!         if (P%InitPower%WantExternalScalar.or.P%InitPower%WantExternalTensor) &
!              call FreeExternalData(P%InitPower)
! !end addon
    
          end program driver
  
diff -r -c -b -B -N CAMB/Makefile CAMB_addon/Makefile
*** CAMB/Makefile	2006-03-14 11:07:24.000000000 +0000
--- CAMB_addon/Makefile	2006-03-26 19:13:14.000000000 +0100
***************
*** 1,57 ****
! #CAMB Makefile
! #Edit for your compiler
! #Note there are many ifc versions, some of which behave oddly
! 
! #Intel 8 , -openmp toggles mutli-processor:
! F90C     = ifort
! FFLAGS = -openmp -ip -O2 -Vaxlib -W0 -WB -fpp2
! 
! # Intel 9 on IA-64 (eg. COSMOS)
! # (do "module load icomp90" before compiling)
! #F90C = ifort
! #FFLAGS = -openmp -fpp2 -w -O3 -ip -mP2OPT_hlo_prefetch=F
! 
! #Intel ifc, add -openmp for multi-processor (some have bugs):
! #F90C     = ifc
! #FFLAGS = -O2 -Vaxlib -ip -W0 -WB -quiet -fpp2
! #some systems can can also add e.g. -tpp7 -xW
! 
! #SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems.
! #F90C     = f90
! #FFLAGS  = -Ofast -mp
! 
! #Digital/Compaq fortran, -omp toggles multi-processor
! #F90C    = f90
! #FFLAGS  = -omp -O4 -arch host -math_library fast -tune host -fpe1
! 
! #Absoft ProFortran, single processor:
! #F90C     = f95
! #FFLAGS = -O2 -cpu:athlon -s -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_"
! 
! #NAGF95, single processor:
! #F90C     = f95
! #FFLAGS = -DNAGF95 -O3
! 
! #G95 compiler
! #F90C   = g95
! #FFLAGS = -O2
! 
! #Sun, single processor:
! #F90C     = f90
! #FFLAGS = -O5
! 
! #Sun V880
! #F90C = mpf90
! #FFLAGS =  -O4 -openmp -ftrap=%none -dalign -DMPI
! 
! #Sun parallel enterprise:
! #F90C     = f95
! #FFLAGS =  -O2 -xarch=native64 -openmp -ftrap=%none
! #try removing -openmp if get bus errors. -03, -04 etc are dodgy.
! 
! #IBM XL Fortran, multi-processor (run gmake)
! #F90C     = xlf90_r
! #FFLAGS  = -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
  
  
  #Files containing evolution equations initial power spectrum module
--- 1,33 ----
! # >>> DESIGNED FOR GMAKE <<<
! 
! #CAMB System unified Makefile
! 
! ext=$(shell uname | cut -c1-3)
! 
! ifeq ($(ext),IRI)
! F90C=f90
! FFLAGS = -O3 -mp -64
! endif
! 
! ifeq ($(ext),Lin)
! F90C=ifort
! FFLAGS= -O -openmp  -fpp
! endif 
! 
! ifeq ($(ext),OSF)
! F90C=f90
! FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
! endif
! 
! ifeq ($(ext),Sun)
! F90C=f90
! FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
! endif
! 
! ifeq ($(ext),AIX)
! F90C=xlf90_r
! FFLAGS= -O -q64 -qsmp=omp -qthreaded  -qmaxmem=-1 -qfree=f90 -qsuffix=f=f90:cpp=F90
! endif
  
  
  #Files containing evolution equations initial power spectrum module
***************
*** 83,99 ****
  CAMBOBJ       = subroutines.o inifile.o $(POWERSPECTRUM).o recfast.o modules.o \
  	bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o cmbmain.o camb.o
  
! default: camb
  
! all: camb $(CAMBLIB)
  
! camb: $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) -o $@
  
  $(CAMBLIB): $(CAMBOBJ)
  	ar -r $@ $?
  
! camb_fits: writefits.f90 $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) -DWRITE_FITS -o $@
  
  %.o: %.f90
--- 59,75 ----
  CAMBOBJ       = subroutines.o inifile.o $(POWERSPECTRUM).o recfast.o modules.o \
  	bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o cmbmain.o camb.o
  
! default: camb.$(ext)
  
! all: camb.$(ext) $(CAMBLIB)
  
! camb.$(ext): $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) -o $@
  
  $(CAMBLIB): $(CAMBOBJ)
  	ar -r $@ $?
  
! camb_fits.$(ext): writefits.f90 $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) -DWRITE_FITS -o $@
  
  %.o: %.f90
diff -r -c -b -B -N CAMB/modules.f90 CAMB_addon/modules.f90
*** CAMB/modules.f90	2006-03-16 11:17:06.000000000 +0000
--- CAMB_addon/modules.f90	2006-03-17 22:37:18.000000000 +0000
***************
*** 186,192 ****
        integer, parameter:: l0max=4000
  
  !     lmax is max possible number of l's evaluated
!       integer, parameter :: lmax_arr = 100+l0max/7
   
          contains
        
--- 186,197 ----
        integer, parameter:: l0max=4000
  
  !     lmax is max possible number of l's evaluated
! !addon
! !      integer, parameter :: lmax_arr = 100+l0max/7
!       integer, parameter :: lmax_arr = 100 + l0max
!       logical :: DoClsInterpol=.true.
!       logical :: ComputeEveryl=.true.
! !end addon
   
          contains
        
***************
*** 539,544 ****
--- 544,567 ----
          integer lind, lvar, step,top,bot,ls(lmax_arr)
          real(dl) AScale
      
+ !addon    
+         if (ComputeEveryl) then
+            write(*,*) 'Warning: computing every l'           
+            if (max_l.gt.lmax_arr) then
+               write(*,*)'lmax_arr maxl',lmax_arr,max_l
+               stop
+            endif
+            lind=0
+            do lvar=2, max_l
+               lind=lind+1
+               ls(lind)=lvar 
+            end do
+ 
+            lSet%l0=lind
+            lSet%l = ls(1:lind)
+ 
+         else
+ !end addon
          Ascale=scale/lSampleBoost       
        
          lind=0
***************
*** 641,646 ****
--- 664,673 ----
          lSet%l0=lind
          lSet%l(1:lind) = ls(1:lind)
          
+ !addon
+         endif
+ !end addon
+         
        end subroutine initlval
  
        subroutine InterpolateClArr(lSet,iCl, all_Cl, max_ind)
***************
*** 660,665 ****
--- 687,694 ----
        xl = real(lSet%l(1:lSet%l0),dl)
        call spline(xl,iCl(1),max_ind,cllo,clhi,ddCl(1))
       
+ !addon
+       if (DoClsInterpol) then
              llo=1
              do il=lmin,lSet%l(max_ind)
                 xi=il
***************
*** 676,681 ****
--- 705,727 ----
                
              end do
  
+       else
+ 
+          llo=1
+          do il=2,lSet%l(max_ind)
+             xi=il
+             if ((xi > lSet%l(llo+1)).and.(llo < max_ind)) then
+                llo=llo+1
+             end if
+ 
+             all_Cl(il) = iCl(llo+1)
+ 
+          end do
+          all_Cl(2)=iCl(1)
+       endif
+ !end addon
+ 
+ 
        end subroutine InterpolateClArr
  
   
diff -r -c -b -B -N CAMB/params.ini CAMB_addon/params.ini
*** CAMB/params.ini	2005-03-30 16:05:28.000000000 +0100
--- CAMB_addon/params.ini	2006-02-20 17:31:02.000000000 +0000
***************
*** 64,69 ****
--- 64,76 ----
  
  #Initial power spectrum, amplitude, spectral index and running
  initial_power_num         = 1
+ 
+ #read the power spectra form ascii files in the form (Mpc^-1): k P(k)
+ use_external_scalar       = F
+ use_external_tensor       = F
+ scalar_input_file(1)      = scalar_power.dat
+ tensor_input_file(1)      = tensor_power.dat
+ 
  scalar_amp(1)             = 2.3e-9
  scalar_spectral_index(1)  = 1
  scalar_nrun(1)            = 0
***************
*** 100,105 ****
--- 107,118 ----
  #Matter power spectrum output against k/h in units of h^3 Mpc^{-3}
  transfer_matterpower(1) = matterpower.dat
  
+ #Use spline interpolation for the Cls
+ do_cls_interpol = T
+ 
+ #Remove undersampling of the Cls: if true lSampleBoost has no effect
+ compute_every_l = F
+ 
  
  #Output files not produced if blank. make camb_fits to use use the FITS setting.
  scalar_output_file = scalCls.dat
diff -r -c -b -B -N CAMB/power_tilt.f90 CAMB_addon/power_tilt.f90
*** CAMB/power_tilt.f90	2005-03-31 00:37:13.000000000 +0100
--- CAMB_addon/power_tilt.f90	2006-02-20 20:35:53.000000000 +0000
***************
*** 51,60 ****
         real(dl) k_0_scalar, k_0_tensor
         real(dl) ScalarPowerAmp(nnmax)
   
        end Type InitialPowerParams
  
!       real(dl) curv  !Curvature contant, set in InitializePowers     
  
        
        Type(InitialPowerParams) :: P
    
--- 51,80 ----
         real(dl) k_0_scalar, k_0_tensor
         real(dl) ScalarPowerAmp(nnmax)
   
+ !addon 
+        logical :: WantExternalScalar
+        logical :: WantExternalTensor
+        character(80) :: ScalarPowerFileName(nnmax),TensorPowerFileName(nnmax)
+ !end addon      
+ 
        end Type InitialPowerParams
  
! !addon
!       type ExternalData
! !only pointers may be defered shape in derived type
!          real(dl) :: ScalarMomentumMin, ScalarMomentumMax
!          real(dl), dimension(:), pointer :: PtrScalarPower, PtrScalarLnMomentum
!          real(dl), dimension(:), pointer :: PtrScalarPower2
!          real(dl) :: TensorMomentumMin, TensorMomentumMax
!          real(dl), dimension(:), pointer :: PtrTensorPower, PtrTensorLnMomentum
!          real(dl), dimension(:), pointer :: PtrTensorPower2
!       end type ExternalData
  
+       type(ExternalData), dimension(:), pointer  :: PtrExternalData => null()
+ !end addon      
+ 
+ 
+       real(dl) curv  !Curvature contant, set in InitializePowers     
        
        Type(InitialPowerParams) :: P
    
***************
*** 63,68 ****
--- 83,91 ----
        public  InitializePowers, ScalarPower, TensorPower
        public nnmax,InitialPowerParams,Power_Descript, Power_Name, SetDefPowerParams
  
+ !addon
+       public FreeExternalData
+ !end addon
      contains
         
    
***************
*** 77,82 ****
--- 100,110 ----
           AP%k_0_scalar = 0.05
           AP%k_0_tensor = 0.05
           AP%ScalarPowerAmp = 1
+ !addon
+ !         nullify(PtrExternalData)
+          AP%WantExternalScalar = .false.
+          AP%WantExternalTensor = .false.
+ !end addon        
  
         end subroutine SetDefPowerParams
  
***************
*** 87,92 ****
--- 115,128 ----
  
          real(dl) acurv
  
+ !addon
+         integer :: i,ioflag,j
+         integer :: nPowerUnit=111
+         integer :: nData
+         real(dl), save :: sbuffer1,sbuffer2
+         real(dl), parameter :: spl_large=1.e40_dl     
+ !end addon
+ 
          P = AParamSet
          if (P%nn > nnmax) then
             write (*,*) 'To use ',P%nn,'power spectra you need to increase'
***************
*** 97,105 ****
--- 133,295 ----
  
  !Write implementation specific code here...
  
+ !addon
+ 
+ !create and point to P%nn ExternalData objects (pointers and reals)
+         if ((P%WantExternalScalar).or.(P%WantExternalTensor)) then
+            if (.not.associated(PtrExternalData)) then
+               allocate(PtrExternalData(P%nn)) 
+            else
+               write(*,*)'InitializePowers: PtrExternalData already Initialized!'
+               return
+            endif
+         endif
+          
+         if (P%WantExternalScalar) then
+                         
+            do i=1,P%nn
+                
+               write(*,*)'using primordial scalar power spectrum: ',P%ScalarPowerFileName(i)
+                
+               open(unit=nPowerUnit,file=P%ScalarPowerFileName(i),action='read' &
+                    ,form='formatted',status='old')             
+ !count the number of records
+               nData=0
+               do                                    
+                  read(nPowerUnit,*,iostat=ioflag)                  
+                  if (ioflag.lt.0) exit
+                  nData=nData + 1                                    
+               enddo
+ !create and point to
+               allocate(PtrExternalData(i)%PtrScalarPower(nData) & 
+                    ,PtrExternalData(i)%PtrScalarLnMomentum(nData))
+               
+               rewind(nPowerUnit)
+               nData=1
+               do                                   
+                  read(nPowerUnit,*,iostat=ioflag) sbuffer1,sbuffer2
+                  if (ioflag.lt.0) exit
+                  PtrExternalData(i)%PtrScalarLnMomentum(nData) = log(sbuffer1)
+                  PtrExternalData(i)%PtrScalarPower(nData) = sbuffer2      
+                  nData=nData+1
+               enddo
+               nData=nData - 1 
+               close(nPowerUnit,iostat=ioflag)
+             
+ !keep trace of the accessible range              
+               PtrExternalData(i)%ScalarMomentumMin = exp(PtrExternalData(i)%PtrScalarLnMomentum(1))
+               PtrExternalData(i)%ScalarMomentumMax = exp(PtrExternalData(i)%PtrScalarLnMomentum(nData))
+               
+ 
+ !2nd order derivatives allocation
+               allocate(PtrExternalData(i)%PtrScalarPower2(nData))
+ !splining in log k
+               call spline(PtrExternalData(i)%PtrScalarLnMomentum(:) &
+                    ,PtrExternalData(i)%PtrScalarPower(:),nData, spl_large &
+                    , spl_large,PtrExternalData(i)%PtrScalarPower2(:))
+               
+            enddo
+ 
+         endif
+ 
+ 
+ !the same for tensors
+ 
+         if (P%WantExternalTensor) then
+             
+            do i=1,P%nn
+               
+               write(*,*)'using primordial tensor power spectrum: ',P%TensorPowerFileName(i)
+ 
+               open(unit=nPowerUnit,file=P%TensorPowerFileName(i) &
+                    ,action='read',form='formatted',status='old')     
+ !count the number of records
+               nData=0
+               do                                    
+                  read(nPowerUnit,*,iostat=ioflag)                  
+                  if (ioflag.lt.0) exit
+                  nData=nData + 1                                    
+               enddo
+ !store the data
+               allocate(PtrExternalData(i)%PtrTensorPower(nData) & 
+                    ,PtrExternalData(i)%PtrTensorLnMomentum(nData))
+                
+               rewind(nPowerUnit)
+               nData=1
+               do                                   
+                  read(nPowerUnit,*,iostat=ioflag) sbuffer1,sbuffer2
+                  if (ioflag.lt.0) exit
+                  PtrExternalData(i)%PtrTensorLnMomentum(nData) = log(sbuffer1)
+                  PtrExternalData(i)%PtrTensorPower(nData) = sbuffer2      
+                  nData=nData+1
+               enddo
+               nData=nData - 1 
+               close(nPowerUnit)
+ 
+ !keep trace of the accessible range              
+               PtrExternalData(i)%TensorMomentumMin = exp(PtrExternalData(i)%PtrTensorLnMomentum(1))
+               PtrExternalData(i)%TensorMomentumMax = exp(PtrExternalData(i)%PtrTensorLnMomentum(nData))
+               
+ 
+ !2nd order derivatives
+               allocate(PtrExternalData(i)%PtrTensorPower2(nData))
+ !splining in log k
+               call spline(PtrExternalData(i)%PtrTensorLnMomentum(:) &
+                    ,PtrExternalData(i)%PtrTensorPower(:),nData, spl_large &
+                    , spl_large,PtrExternalData(i)%PtrTensorPower2(:))
+               
+            enddo
+ 
+             
+         endif
+          
         end subroutine InitializePowers
         
  
+ 
+        subroutine FreeExternalData(AP)
+          type (InitialPowerParams), intent(in) :: AP
+          integer :: i
+ 
+          if (AP%nn.ne.size(PtrExternalData,1)) then
+             write(*,*)'InitialPower.FreeExternalData : wrong initial power params'
+             stop
+          endif
+ 
+          if (AP%WantExternalScalar) then
+             do i=1, AP%nn               
+ !free heap arrays and disassociate ExternalData pointers            
+                deallocate(PtrExternalData(i)%PtrScalarPower)
+                deallocate(PtrExternalData(i)%PtrScalarPower2)
+                deallocate(PtrExternalData(i)%PtrScalarLnMomentum)
+ !should still point to ExternalData objects
+                if (.not.associated(PtrExternalData)) stop 'InitialPower.FreeExternalData : ooops'
+             enddo
+          endif
+          if (AP%WantExternalTensor) then
+             do i=1, AP%nn               
+                deallocate(PtrExternalData(i)%PtrTensorPower)
+                deallocate(PtrExternalData(i)%PtrTensorPower2)
+                deallocate(PtrExternalData(i)%PtrTensorLnMomentum)
+                if (.not.associated(PtrExternalData)) stop 'InitialPower.FreeExternalData : ooops'
+             enddo
+          endif
+ !Free remaining objects and disassociate pointer to ExternalData
+          if (associated(PtrExternalData)) then
+             deallocate(PtrExternalData)
+             PtrExternalData => null()
+          else
+             if (AP%WantExternalScalar.or.AP%WantExternalTensor) then
+ !should never happen               
+                stop 'InitialPower.FreeExternalData : ooops'
+             endif
+          endif
+ 
+        end subroutine FreeExternalData
+ 
+ !end addon
+ 
+ 
        function ScalarPower(k,in)
  
         !"in" gives the index of the power to return for this k
***************
*** 120,128 ****
--- 310,340 ----
          real(dl) ScalarPower,k, lnrat
          integer in
  
+ !addon
+         integer :: nData       
+ 
+         if (P%WantExternalScalar) then
+            nData=size(PtrExternalData(in)%PtrScalarLnMomentum,1)
+ 
+            if ((k.ge.PtrExternalData(in)%ScalarMomentumMin) &
+                 .and.(k.le.PtrExternalData(in)%ScalarMomentumMax)) then
+               call splinterpol(PtrExternalData(in)%PtrScalarLnMomentum &
+                    ,PtrExternalData(in)%PtrScalarPower &
+                    ,PtrExternalData(in)%PtrScalarPower2,nData,log(k),ScalarPower)
+ !confusing              ScalarPower = P%ScalarPowerAmp(in)*ScalarPower
+            else
+               write(*,*)'InitPower.ScalarPower() : out of range: k= ',k
+               ScalarPower=0_dl              
+            endif
+ 
+         else
+ 
            lnrat = log(k/P%k_0_scalar)
            ScalarPower=P%ScalarPowerAmp(in)*exp((P%an(in)-1)*lnrat + P%n_run(in)/2*lnrat**2)   
   
+         endif
+ !end addon
+ 
        end function ScalarPower
  
        
***************
*** 143,149 ****
--- 355,384 ----
     
          integer in
  
+ !addon
+         integer :: nData
+         
+ 
+         if (P%WantExternalTensor) then
+            nData=size(PtrExternalData(in)%PtrTensorLnMomentum,1)
+            if ((k.ge.PtrExternalData(in)%TensorMomentumMin) &
+                 .and.(k.le.PtrExternalData(in)%TensorMomentumMax)) then
+               call splinterpol(PtrExternalData(in)%PtrTensorLnMomentum &
+                    ,PtrExternalData(in)%PtrTensorPower &
+                    ,PtrExternalData(in)%PtrTensorPower2,nData,log(k),TensorPower)
+ 
+ !confusing              TensorPower = P%rat(in)*P%ScalarPowerAmp(in)*TensorPower
+            else
+               write(*,*)'InitPower.TensorPower() : out of range: k= ',k
+               TensorPower=0_dl               
+            endif
+         else
+ 
          TensorPower=P%rat(in)*P%ScalarPowerAmp(in)*exp(P%ant(in)*log(k/P%k_0_tensor))
+ 
+         endif
+ !end addon
+         
          if (curv < 0) TensorPower=TensorPower*tanh(PiByTwo*sqrt(-k**2/curv-3)) 
  
        end function TensorPower
diff -r -c -b -B -N CAMB/subroutines.f90 CAMB_addon/subroutines.f90
*** CAMB/subroutines.f90	2005-03-31 00:37:13.000000000 +0100
--- CAMB_addon/subroutines.f90	2006-02-20 17:04:32.000000000 +0000
***************
*** 962,964 ****
--- 962,997 ----
  !  end abort action
  !
        end subroutine dverk
+ 
+ 
+ !addon
+ !the numrec splint renamed
+ 
+ SUBROUTINE splinterpol(xa,ya,y2a,n,x,y)
+   use Precision
+   INTEGER n
+   REAL(dl) x,y,xa(n),y2a(n),ya(n)
+   INTEGER k,khi,klo
+   REAL(dl) a,b,h
+   klo=1
+   khi=n
+ 1 if (khi-klo.gt.1) then
+      k=(khi+klo)/2
+      if(xa(k).gt.x)then
+         khi=k
+      else
+         klo=k
+      endif
+      goto 1
+   endif
+   h=xa(khi)-xa(klo)
+   if (h.eq.0.) stop 'bad xa input in splint'
+   a=(xa(khi)-x)/h
+   b=(x-xa(klo))/h
+   y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))&
+        *(h**2)/6.
+   return
+ END SUBROUTINE splinterpol
+ 
+ 
+ !end addon
