diff -c -r -b -B cosmomc/camb/Makefile cosmomc_addon/camb/Makefile
*** cosmomc/camb/Makefile	Tue Mar 29 13:40:54 2005
--- cosmomc_addon/camb/Makefile	Mon Feb 20 13:33:12 2006
***************
*** 1,49 ****
! #CAMB Makefile
! #Edit for your compiler
! #Note there are many ifc versions, some of which behave oddly
! 
! #Intel 8 , -openmp toggles mutli-processor:
! #(do "module load dev8" to get v8 compilers on COSMOS)
! F90C     = ifort
! FFLAGS = -openmp -ip -O2 -Vaxlib -W0 -WB -fpp2
! 
! #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 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
***************
*** 54,61 ****
  
  #Driver program
  DRIVER        = inidriver.F90
- #DRIVER        = sigma8.f90
- #DRIVER        = tester.f90
  
  #Settings for building camb_fits
  #Location of FITSIO and name of library
--- 38,43 ----
***************
*** 75,91 ****
  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
--- 57,73 ----
  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 -c -r -b -B cosmomc/camb/inidriver.F90 cosmomc_addon/camb/inidriver.F90
*** cosmomc/camb/inidriver.F90	Wed Mar 30 18:37:13 2005
--- cosmomc_addon/camb/inidriver.F90	Mon Feb 20 15:32:38 2006
***************
*** 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 -c -r -b -B cosmomc/camb/modules.f90 cosmomc_addon/camb/modules.f90
*** cosmomc/camb/modules.f90	Fri Apr  8 09:46:12 2005
--- cosmomc_addon/camb/modules.f90	Mon Feb 20 12:03:38 2006
***************
*** 183,189 ****
        integer, parameter:: l0max=4000
  
  !     lmax is max possible number of l's evaluated
!       integer, parameter :: lmax_arr = 100+l0max/7
   
          contains
        
--- 183,194 ----
        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
        
***************
*** 536,541 ****
--- 541,564 ----
          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
***************
*** 638,643 ****
--- 661,670 ----
          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)
***************
*** 657,662 ****
--- 684,691 ----
        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=2,lSet%l(max_ind)
                 xi=il
***************
*** 673,678 ****
--- 702,724 ----
                
              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 -c -r -b -B cosmomc/camb/params.ini cosmomc_addon/camb/params.ini
*** cosmomc/camb/params.ini	Wed Mar 30 10:05:28 2005
--- cosmomc_addon/camb/params.ini	Mon Feb 20 12:31:02 2006
***************
*** 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 -c -r -b -B cosmomc/camb/power_tilt.f90 cosmomc_addon/camb/power_tilt.f90
*** cosmomc/camb/power_tilt.f90	Wed Mar 30 18:37:13 2005
--- cosmomc_addon/camb/power_tilt.f90	Mon Feb 20 15:35:53 2006
***************
*** 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 -c -r -b -B cosmomc/camb/subroutines.f90 cosmomc_addon/camb/subroutines.f90
*** cosmomc/camb/subroutines.f90	Wed Mar 30 18:37:13 2005
--- cosmomc_addon/camb/subroutines.f90	Mon Feb 20 12:04:32 2006
***************
*** 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
diff -c -r -b -B cosmomc/distparams.ini cosmomc_addon/distparams.ini
*** cosmomc/distparams.ini	Thu Mar 31 10:17:45 2005
--- cosmomc_addon/distparams.ini	Mon Feb 20 12:35:03 2006
***************
*** 32,37 ****
--- 32,38 ----
  
  plot_meanlikes = T
  shade_meanlikes = T
+ td_meanlikes = F
  
  # if non-zero, output _thin file, thinned by thin_factor
  thin_factor = 0
diff -c -r -b -B cosmomc/params.ini cosmomc_addon/params.ini
*** cosmomc/params.ini	Thu Jul 21 21:08:37 2005
--- cosmomc_addon/params.ini	Mon Feb 20 12:29:30 2006
***************
*** 110,116 ****
  #if blank this is set from system clock
  rand_seed = 
  
! #Can re-start from the last line of previous run (.txt file)
  continue_from = 
  
  
--- 110,116 ----
  #if blank this is set from system clock
  rand_seed = 
  
! #Can re-start from the last line of previous run (root file)
  continue_from = 
  
  
***************
*** 122,127 ****
--- 122,130 ----
  #increase accuracy_level to run CAMB on higher accuracy
  #(default is about 0.3%, accuracy_level=2 around 0.1% at high l)
  accuracy_level = 1
+ l_accuracy_level = 1
+ l_sample_level = 1
+ camb_compute_every_l = F
  
  #If action = 1
  redo_likelihoods = T
diff -c -r -b -B cosmomc/source/CMB_Cls_simple.f90 cosmomc_addon/source/CMB_Cls_simple.f90
*** cosmomc/source/CMB_Cls_simple.f90	Thu Jul 21 21:37:41 2005
--- cosmomc_addon/source/CMB_Cls_simple.f90	Mon Feb 20 12:24:35 2006
***************
*** 264,269 ****
--- 264,287 ----
     type(CAMBParams)  P 
     integer zix
  
+ !addon
+    AccuracyBoost = AccuracyLevel
+    lAccuracyBoost = lAccuracyLevel
+    lSampleBoost = lSampleLevel
+    ComputeEveryl = CambComputeEveryl
+ 
+    write(*,*)
+    write(*,*) '>>>>>>>>>>>>>> accuracy settings <<<<<<<<<<<<<<'
+    write(*,*) 'Power_Name     =    ',Power_Name
+    write(*,*) 'write_all_cls  =   ',write_all_cls
+    write(*,*) 'AccuracyBoost  = ',AccuracyBoost
+    write(*,*) 'lAccuracyBoost = ',lAccuracyBoost          
+    write(*,*) 'DoClsInterpol  =   ',DoClsInterpol
+    write(*,*) 'ComputeEveryl  =   ',ComputeEveryl
+    write(*,*) '>>>>>>>>>>>>>> accuracy settings <<<<<<<<<<<<<<'
+    write(*,*)
+ !end addon
+ 
          Threadnum =num_threads
          w_lam = -1
          call CAMB_SetDefParams(P)
***************
*** 291,299 ****
              P%Transfer%high_precision=.true.
              P%Transfer%k_per_logint=3*nint(AccuracyLevel)
            end if
!           AccuracyBoost = AccuracyLevel
!           lAccuracyBoost = AccuracyLevel
!           lSampleBoost = AccuracyLevel
            P%AccurateReionization = .true.
          end if
                  
--- 309,320 ----
              P%Transfer%high_precision=.true.
              P%Transfer%k_per_logint=3*nint(AccuracyLevel)
            end if
! !addon
! !moved up for screen output
! !          AccuracyBoost = AccuracyLevel
! !          lAccuracyBoost = AccuracyLevel
! !          lSampleBoost = AccuracyLevel
! !end addon
            P%AccurateReionization = .true.
          end if
                  
diff -c -r -b -B cosmomc/source/GetDist.f90 cosmomc_addon/source/GetDist.f90
*** cosmomc/source/GetDist.f90	Thu Jul 21 21:40:27 2005
--- cosmomc_addon/source/GetDist.f90	Mon Feb 20 12:33:58 2006
***************
*** 66,71 ****
--- 66,74 ----
          logical force_twotail 
          logical smoothing, plot_meanlikes, plot_NDcontours
          logical shade_meanlikes, make_single_samples
+ !addon
+         logical td_meanlikes
+ !end addon
          integer single_thin,cust2DPlots(max_cols**2)
          integer ix_min(max_cols),ix_max(max_cols)
          real center(max_cols),width(max_cols)   
***************
*** 1146,1158 ****
   !             write (aunit,fmt) (/(I, I=ix_min(j), ix_max(j), 1)/)*width(j)+center(j)
                if (DoShade) then
                if (shade_meanlikes) then
                  write (aunit,'(a)') 'load (fullfile(''plot_data'',''' // trim(plotfile) //'_likes''));'
                  write (aunit,*) 'contourf(x1,x2,'//trim(plotfile)//'_likes,64);'
                else     
                  write (aunit,*) 'contourf(x1,x2,pts,64);'
                end if
! 
                write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;'
                end if
                
             
--- 1149,1175 ----
   !             write (aunit,fmt) (/(I, I=ix_min(j), ix_max(j), 1)/)*width(j)+center(j)
                if (DoShade) then
                if (shade_meanlikes) then
+ !addon
                  write (aunit,'(a)') 'load (fullfile(''plot_data'',''' // trim(plotfile) //'_likes''));'
+                  if (td_meanlikes) then
+                     write (aunit,*) 'surfl(x1,x2,'//trim(plotfile)//'_likes);'
+                  else
                  write (aunit,*) 'contourf(x1,x2,'//trim(plotfile)//'_likes,64);'
+                  endif
+ !end addon
                else     
                  write (aunit,*) 'contourf(x1,x2,pts,64);'
                end if
! !addon
! !              write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;'
!               if (td_meanlikes) then
!                  write (aunit,*) 'set(gca,''climmode'',''manual'');shading interp; hold on;'
!                  write (aunit,*) 'grid off;box on;camproj(''perspective'');'
!                  write (aunit,*) 'set(gca,''climmode'',''manual'');alpha(0.9);axis tight'
!               else
                write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;'
+               endif
+ !end addon
                end if
                
             
***************
*** 1471,1476 ****
--- 1488,1496 ----
            write (*,*) 'WARNING: Mapping params - .covmat file is new params.'
  
      shade_meanlikes = Ini_Read_Logical('shade_meanlikes',.false.)   
+ !addon
+     td_meanlikes = Ini_Read_Logical('td_meanlikes',.false.)  
+ !end addon   
  
      matlab_col = Ini_Read_String('matlab_colscheme')
      smoothing = Ini_Read_Logical('smoothing',.true.)
***************
*** 1971,1977 ****
--- 1991,2004 ----
        write (51,*)  'set(gcf, ''PaperUnits'',''inches'');'
       
        if (plot_row < plot_col .and. plot_row*plot_col>9) then
+ !addon
+ !            write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 10]);';
+          if (td_meanlikes) then
+             write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 7 9]);';
+          else
              write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 10]);';
+          endif
+ !end addon           
            else
              write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 8]);';
        end if
diff -c -r -b -B cosmomc/source/Makefile cosmomc_addon/source/Makefile
*** cosmomc/source/Makefile	Thu Sep 30 15:07:02 2004
--- cosmomc_addon/source/Makefile	Mon Feb 20 13:33:34 2006
***************
*** 1,90 ****
! #Settings here are for Intel ifc, uncomment others below:
! #You may need to edit the library paths for MKL for Intel
! #Beware of using optmizations which lose accuracy - may give errors when running
! #With Visual Fortan compile open cosmomc.dsw workspace file instead.
! 
! #Intel MPI (assuming mpif77 set to point to ifc)
! #change -lmkl_ia32 to -lmkl_p3 for MKL versions earlier than 6 (6 needed for ifc 8+) 
! F90C     = mpif77
! FFLAGS = -O2 -Vaxlib -ip -W0 -WB -openmp -fpp -DMPI -quiet
! LAPACKL = -L/opt/intel/mkl61/lib/32 -lmkl_lapack -lmkl_ia32 -lguide -lpthread -static
! INCLUDE = -I../camb 
  
! #Intel fortran 8, check you have the latest update from the Intel web pages
! #See Makefile_intel for ifc 7.1 or lower (some versions have problems)
! #F90C     = ifort
! #FFLAGS = -O2 -Vaxlib -ip -W0 -WB -openmp -fpp
! #LAPACKL = -L/opt/intel/mkl61/lib/32 -lmkl_lapack -lmkl_ia32 -lguide -lpthread -static
! #INCLUDE = -I../camb 
! 
! #Intel fortran 8 on COSMOS (do "module load dev8" to get v8 compilers)
! #F90C     = ifort
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp2 -lmpi -lscs  -ldl -DMPI
! #LAPACKL = -lmkl_lapack -lmkl_ipf -limf
! #INCLUDE = -I../camb
! 
! #SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems.
! #Not various versions of the compiler are buggy giving erroneous seg faults with -mp.
! #Version 7.3 is OK, currently version 7.4 is bugged, as are some earlier versions.
! #F90C     = f90
  #LAPACKL = -lcomplib.sgimath
! #INCLUDE = -I../camb 
! #FFLAGS  = -Ofast -mp
  
! #Digital/Compaq fortran, -omp toggles multi-processor
! #F90C    = f90
! #FFLAGS  = -omp -O4 -arch host -math_library fast -tune host -fpe1
! #LAPACKL  = -lcxml
! #INCLUDE = -I../camb 
! 
! #Intel Itanium efc, -openmp toggles multi-processor:
! #F90C     = efc
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp
! #LAPACKL = -L/opt/intel/mkl/lib/64 -lmkl_lapack -lmkl_ipf -lguide -lpthread -static
! #INCLUDE = -I../camb 
! 
! #Intel Itanium efc MPI:
! # mpirun -np 4 ./cosmomc params.ini
! #F90C     = efc
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp -lmpi -lscs  -ldl -DMPI
! #LAPACKL = -L/opt/intel/mkl/lib/64 -lmkl_lapack -lmkl_ipf -limf
! #INCLUDE = -I../camb
! 
! #Absoft ProFortran, single processor, set -cpu:[type] for your local system 
! #F90C     = f95
! #FFLAGS = -O2 -s -cpu:athlon -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_" 
! #LAPACKL =  -llapack -lblas -lg2c 
! #INCLUDE = -I../camb -p../camb
! 
! #NAGF95, single processor:
! #F90C     = f95
! #FFLAGS = -DNAGF95 -O3
! #LAPACKL = -llapack -lblas -lg2c
! #INCLUDE = -I../camb 
! 
! #Sun, single processor:
! #F90C     = f90
! #FFLAGS = -O5
! #LAPACKL = -lsunperf -lfsu
! #INCLUDE = -I../camb -M../camb
! 
! #Sun parallel enterprise:
! #F90C     = f95
! #FFLAGS =  -O4 -xarch=native64 -openmp -ftrap=%none
! #LAPACKL = -lsunperf -lfsu
! #INCLUDE = -I../camb -M../camb
! 
! #IBM XL Fortran, multi-processor (run "module load lapack" then run "gmake")
! #F90C     = xlf90_r $(LAPACK)
! #FFLAGS  = -WF,-DIBMXL -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
! #INCLUDE = -I../camb -lessl
! #LAPACKL =
  
  
  PROPOSE = propose.o
  CLSFILE = CMB_Cls_simple.o
  
! #Can use params_H if you prefer more generic parameters
  PARAMETERIZATION = params_CMB.o
  
  CLSLIB  = -L../camb -lcamb $(LAPACKL)
--- 1,50 ----
! # >>> DESIGNED FOR GMAKE <<<
! 
! # Unified Systems makefile for COSMOMC
! # Add FLAGS -DMPI for using MPI
  
! ext=$(shell uname | cut -c1-3)
! 
! ifeq ($(ext),IRI)
! F90C= f90
! FFLAGS= -Ofast -mp -n32 -LANG:recursive=ON -lmpi -DMPI
  #LAPACKL = -lcomplib.sgimath
! LAPACKL = -lcomplib.sgimath_mp
! INCLUDE = -I../camb
! endif
  
! ifeq ($(ext),Lin)
! F90C = ifort
! FFLAGS = -O -openmp -fpp
! LAPACKL = -lmkl_lapack -lmkl
! INCLUDE = -I../camb -L/usr/local/lib 
! endif 
! 
! ifeq ($(ext),OSF)
! F90C=f90
! FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
! LAPACKL  = -lcxml
! INCLUDE = -I../camb 
! endif
  
+ ifeq ($(ext),Sun)
+ F90C=f90
+ FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
+ LAPACKL = -lsunperf -lfsu
+ INCLUDE = -I../camb -M../camb
+ endif
+ 
+ ifeq ($(ext),AIX)
+ F90C     = mpxlf90_r
+ FFLAGS  = -O4 -WF,-DIBMXL,-DMPI -qstrict -qsmp=omp -qmaxmem=-1 -qsuffix=f=f90:cpp=F90
+ INCLUDE = -I../camb -lessl
+ LAPACKL =
+ endif
  
  PROPOSE = propose.o
  CLSFILE = CMB_Cls_simple.o
  
! #PARAMETRIZATION = params_H.o
  PARAMETERIZATION = params_CMB.o
  
  CLSLIB  = -L../camb -lcamb $(LAPACKL)
***************
*** 99,107 ****
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
  
! default: cosmomc
  
! all : cosmomc getdist
  
  .f.o:
  	f77 $(F90FLAGS) -c $<
--- 59,67 ----
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
  
! default: cosmomc.$(ext)
  
! all : cosmomc.$(ext) getdist.$(ext)
  
  .f.o:
  	f77 $(F90FLAGS) -c $<
***************
*** 112,124 ****
  %.o: %.F90
  	$(F90C) $(F90FLAGS) -c $*.F90
  
! 
! cosmomc: $(OBJFILES) 	
! 	$(F90C) -o ../cosmomc $(OBJFILES) $(CLSLIB) $(F90FLAGS)
! 
  
  clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist: $(DISTFILES)
! 	$(F90C) -o ../getdist $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
--- 72,82 ----
  %.o: %.F90
  	$(F90C) $(F90FLAGS) -c $*.F90
  
! cosmomc.$(ext): $(OBJFILES) 	
! 	$(F90C) -o ../$@ $(OBJFILES) $(CLSLIB) $(F90FLAGS)
  
  clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist.$(ext): $(DISTFILES)
! 	$(F90C) -o ../$@ $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
Only in cosmomc_addon/source: distparams.ini
diff -c -r -b -B cosmomc/source/driver.F90 cosmomc_addon/source/driver.F90
*** cosmomc/source/driver.F90	Thu Jul 21 21:39:48 2005
--- cosmomc_addon/source/driver.F90	Tue Feb 21 06:43:07 2006
***************
*** 16,23 ****
          
  ! This is a driving routine that illustrates the use of the program.
          character(LEN=Ini_max_string_len) InputFile, LogFile
! 
! 
          logical bad, est_bfp_before_covmat
          integer numsets, nummpksets, i, numtoget, action
          character(LEN=Ini_max_string_len) rootname, baseroot, filename(100), mpk_filename(100), numstr
--- 16,24 ----
          
  ! This is a driving routine that illustrates the use of the program.
          character(LEN=Ini_max_string_len) InputFile, LogFile
! !addon
!         logical :: keepFile = .false.
! !end addon
          logical bad, est_bfp_before_covmat
          integer numsets, nummpksets, i, numtoget, action
          character(LEN=Ini_max_string_len) rootname, baseroot, filename(100), mpk_filename(100), numstr
***************
*** 52,65 ****
  
          propose_scale = Ini_Read_Real('propose_scale',2.4)
          AccuracyLevel = Ini_Read_Real('accuracy_level',1.)
  
  #ifdef MPI 
          
!         if (instance /= 0) call DoStop('With MPI should not have second parameter')
          call mpi_comm_rank(mpi_comm_world,MPIrank,ierror)
  
          instance = MPIrank +1 !start at 1 for chains
!         write (numstr,*) instance
          rand_inst = instance
          if (ierror/=MPI_SUCCESS) call DoStop('MPI fail')
  
--- 53,85 ----
  
          propose_scale = Ini_Read_Real('propose_scale',2.4)
          AccuracyLevel = Ini_Read_Real('accuracy_level',1.)
+ !addon
+ 	lAccuracyLevel = Ini_Read_Real('l_accuracy_level',1.)
+         lSampleLevel = Ini_Read_Real('l_sample_level',1.)
+ 	CambComputeEveryl = Ini_Read_Logical('camb_compute_every_l',.false.)
+ 
+         if (Ini_Read_String('continue_from') /= '') then
+            keepFile = .true.
+         else
+            keepFile = .false.
+         endif
+ !end addon
  
  #ifdef MPI 
          
! !addon                        
! !        if (instance /= 0) call DoStop('With MPI should not have second parameter')
!         instance_shift=instance
!         write(*,*)'MPI file names shifted by: ',instance_shift
! !end addon
          call mpi_comm_rank(mpi_comm_world,MPIrank,ierror)
  
          instance = MPIrank +1 !start at 1 for chains
! !addon
! !        write (numstr,*) instance
!         instance_shift = MPIrank + instance_shift
!         write (numstr,*) instance_shift        
! !end addon
          rand_inst = instance
          if (ierror/=MPI_SUCCESS) call DoStop('MPI fail')
  
***************
*** 99,130 ****
          FeedBack = Ini_Read_Int('feedback',0)
          FileChangeIni = trim(rootname)//'.read'
  
!         if (action /= 1) then
!         LogFile = trim(rootname)//'.log'
! 
!         if (LogFile /= '') then
!          logfile_unit = 49
!          open(unit=logfile_unit,file=LogFile,form='formatted',status='replace')
!          
!         else
!          logfile_unit = 0
!         end if
! 
!         outfile_unit = 48
!         open(unit=outfile_unit,file=trim(rootname)//'.txt',form='formatted',status='replace')
! 
!         indep_sample = Ini_Read_Int('indep_sample')
!         if (indep_sample /=0) then
!           indepfile_unit = 47
!           open(unit=indepfile_unit,file=trim(rootname)//'.data',form='unformatted',status='replace')
!         end if
!  
          Ini_fail_on_not_found = .false.
!         burn_in = Ini_Read_Int('burn_in',0)     
!         sampling_method = Ini_Read_Int('sampling_method',sampling_metropolis)
!         if (sampling_method > 4 .or. sampling_method<1) call DoStop('Unknown sampling method')
!         if (sampling_method==4) directional_grid_steps = Ini_Read_Int('directional_grid_steps',20)
!         end if
  
          numstr = Ini_Read_String('rand_seed')
          if (numstr /= '') then
--- 119,128 ----
          FeedBack = Ini_Read_Int('feedback',0)
          FileChangeIni = trim(rootname)//'.read'
  
! !addon
! !files opening moved after call Initialize()
          Ini_fail_on_not_found = .false.
! !end addon
  
          numstr = Ini_Read_String('rand_seed')
          if (numstr /= '') then
***************
*** 178,183 ****
--- 176,231 ----
  
          call Initialize(Params)
  
+ !addon
+ !otherwise conflicts between reading old data on the same files
+         if (action /= 1) then
+            LogFile = trim(rootname)//'.log'
+ 
+            if (LogFile /= '') then
+               logfile_unit = 49
+               if (keepFile) then
+                  open(unit=logfile_unit,file=LogFile,form='formatted' &
+                       ,status='old',position='append')
+               else
+                  open(unit=logfile_unit,file=LogFile,form='formatted' &
+                       ,status='replace')
+               endif
+            else
+               logfile_unit = 0
+            end if
+ 
+            outfile_unit = 48
+ 
+            if (keepFile) then 
+               open(unit=outfile_unit,file=trim(rootname)//'.txt',form='formatted' &
+                    ,status='old',position='append')
+            else
+               open(unit=outfile_unit,file=trim(rootname)//'.txt',form='formatted' &
+                    ,status='replace')
+            endif
+ 
+            indep_sample = Ini_Read_Int('indep_sample')
+            if (indep_sample /=0) then
+               indepfile_unit = 47
+               if (keepFile) then
+                  open(unit=indepfile_unit,file=trim(rootname)//'.data',form='unformatted'&
+                       ,status='old',position='append')
+               else
+                  open(unit=indepfile_unit,file=trim(rootname)//'.data',form='unformatted'&
+                       ,status='replace')
+               endif
+            end if
+  
+            Ini_fail_on_not_found = .false.
+            burn_in = Ini_Read_Int('burn_in',0)     
+            sampling_method = Ini_Read_Int('sampling_method',sampling_metropolis)
+            if (sampling_method > 4 .or. sampling_method<1) &
+                 call DoStop('Unknown sampling method')
+            if (sampling_method==4) directional_grid_steps &
+                 = Ini_Read_Int('directional_grid_steps',20)
+         end if
+ !end addon
+         
          numsets = Ini_Read_Int('cmb_numdatasets')
          do i= 1, numsets
           filename(i) = Ini_Read_String(numcat('cmb_dataset',i)) 
diff -c -r -b -B cosmomc/source/paramdef.F90 cosmomc_addon/source/paramdef.F90
*** cosmomc/source/paramdef.F90	Thu Jul 21 21:02:59 2005
--- cosmomc_addon/source/paramdef.F90	Tue Feb 21 06:55:43 2006
***************
*** 104,109 ****
--- 104,112 ----
          integer i
          character(LEN=5000) InLine
          character(LEN=120) prop_mat
+ !addon
+         character(LEN=Ini_max_string_len) numstr
+ !end addon
          real center, wid, mult, like
          logical restart
          real pmat(num_params,num_params)
***************
*** 119,126 ****
          restart = InLine /= ''
     
          if (restart) then
             InLine = LastFileLine(InLine)
-       
             read(InLine, *) mult, like, Params%P(1:num_params)
             StartLike = like
          end if
--- 122,138 ----
          restart = InLine /= ''
          
          if (restart) then
+ !addon
+ #ifdef MPI                      
+            write (numstr,*) instance_shift
+            InLine = trim(InLine)//'_'//trim(adjustl(numstr))//'.txt'           
+            write(*,*)'Processus MPIRank =',MpiRank
+ #else
+            InLine = trim(InLine)//'.txt'
+ #endif           
+            write(*,*)'appends data on file: ',trim(InLine)
+ !end addon
             InLine = LastFileLine(InLine)           
             read(InLine, *) mult, like, Params%P(1:num_params)
             StartLike = like
          end if
Only in cosmomc_addon/source: params.ini
diff -c -r -b -B cosmomc/source/settings.f90 cosmomc_addon/source/settings.f90
*** cosmomc/source/settings.f90	Thu Jul 21 21:03:20 2005
--- cosmomc_addon/source/settings.f90	Mon Feb 20 12:27:32 2006
***************
*** 5,10 ****
--- 5,15 ----
    implicit none
  
    real :: AccuracyLevel = 1.
+ !addon
+   real :: lAccuracyLevel = 1.
+   real :: lSampleLevel = 1.
+   logical :: CambComputeEveryl = .false.
+ !end addon
    !Set to >1 to use CAMB etc on higher accuracy settings. 
    !Does not affect MCMC (except making it all slower)
  
***************
*** 43,48 ****
--- 48,56 ----
  
    integer :: num_threads = 0
    integer :: instance = 0
+ !addon
+   integer :: instance_shift = 0
+ !end addon
    integer :: MPIchains = 1, MPIrank = 0
  
    logical :: Use_LSS = .true.
