diff -r -c -b -B cosmomc/camb/inidriver.F90 cosmomc_addon/camb/inidriver.F90
*** cosmomc/camb/inidriver.F90	2008-06-09 19:20:30.000000000 +0200
--- cosmomc_addon/camb/inidriver.F90	2008-07-23 17:40:47.000000000 +0200
***************
*** 244,249 ****
--- 244,258 ----
         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)
***************
*** 293,298 ****
--- 302,311 ----
  
          call CAMB_cleanup             
    
+         if (P%InitPower%WantExternalScalar.or.P%InitPower%WantExternalTensor) &
+              call FreeExternalData(P%InitPower)
+ 
+ 
          end program driver
  
  
diff -r -c -b -B cosmomc/camb/Makefile cosmomc_addon/camb/Makefile
*** cosmomc/camb/Makefile	2008-06-09 19:07:38.000000000 +0200
--- cosmomc_addon/camb/Makefile	2008-05-12 10:12:45.000000000 +0200
***************
*** 1,62 ****
! #CAMB Makefile
! #Edit for your compiler
! #Note there are many ifc versions, some of which behave oddly
! 
! #Intel , -openmp toggles mutli-processor:
! #note version 10.0 gives wrong result for lensed when compiled with -openmp [fixed in 10.1]
! F90C     = ifort
! FFLAGS = -openmp -O2 -ip -W0 -WB -fpp2 -vec_report0
! 
! # 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
! 
! #G95 compiler
! #F90C   = g95
! #FFLAGS = -O2
! 
! #Gfortran compiler: if pre v4.3 add -D__GFORTRAN__
! #F90C     = gfc
! #FFLAGS =  -O2 
! 
! #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
! 
! #PGF90
! #F90C = pgf90
! #FFLAGS = -O2 -DESCAPEBACKSLASH
! 
! #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  = -DESCAPEBACKSLASH -DIBMXL -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
  
  
  #Files containing evolution equations initial power spectrum module
--- 1,31 ----
!  # >>> DESIGNED FOR GMAKE <<<
!  #CAMB System unified Makefile 
! 
! ext=$(shell uname | cut -c1-3)
! ifeq ($(ext),IRI)
! F90C=f90	
! FFLAGS = -Ofast -mp -n32 -LANG:recursive=ON
! endif
! 
! ifeq ($(ext),Lin)
! F90C=ifort
! FFLAGS= -O -openmp #-D__GFORTRAN__
! 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= -O4 -qsmp=omp -qmaxmem=-1 -qstrict -qfree=f90 -qsuffix=f=f90:cpp=F90
! endif
  
  
  #Files containing evolution equations initial power spectrum module
***************
*** 89,97 ****
  CAMBOBJ       = subroutines.o inifile.o $(POWERSPECTRUM).o recfast.o $(REIONIZATION).o modules.o \
  	bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o cmbmain.o camb.o
  
! default: camb
  
! all: camb $(CAMBLIB)
  
  subroutines: utils.o
  $(POWERSPECTRUM): subroutines.o
--- 58,66 ----
  CAMBOBJ       = subroutines.o inifile.o $(POWERSPECTRUM).o recfast.o $(REIONIZATION).o modules.o \
  	bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o cmbmain.o camb.o
  
! default: camb.$(ext)
  
! all: camb.$(ext) $(CAMBLIB)
  
  subroutines: utils.o
  $(POWERSPECTRUM): subroutines.o
***************
*** 106,119 ****
  camb.o: cmbmain.o
  
  
! camb: utils.o $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) utils.o $(CAMBOBJ) $(DRIVER) -o $@
  
  $(CAMBLIB): $(CAMBOBJ)
  	ar -r $@ $?
  
  camb_fits: writefits.f90 utils.o $(CAMBOBJ) $(DRIVER)
! 	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include utils.o $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) -DWRITE_FITS -o $@
  
  %.o: %.f90
  	$(F90C) $(F90FLAGS) -c $*.f90
--- 75,88 ----
  camb.o: cmbmain.o
  
  
! camb.$(ext): utils.o $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) utils.o $(CAMBOBJ) $(DRIVER) -o $@
  
  $(CAMBLIB): $(CAMBOBJ)
  	ar -r $@ $?
  
  camb_fits: writefits.f90 utils.o $(CAMBOBJ) $(DRIVER)
! 	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include tils.o $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) -DWRITE_FITS -o $@
  
  %.o: %.f90
  	$(F90C) $(F90FLAGS) -c $*.f90
diff -r -c -b -B cosmomc/camb/modules.f90 cosmomc_addon/camb/modules.f90
*** cosmomc/camb/modules.f90	2008-06-09 19:20:30.000000000 +0200
--- cosmomc_addon/camb/modules.f90	2008-07-23 17:40:56.000000000 +0200
***************
*** 206,212 ****
        integer, parameter:: l0max=4000
  
  !     lmax is max possible number of l's evaluated
!       integer, parameter :: lmax_arr = 100+l0max/7
   
          contains
        
--- 206,217 ----
        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
        
***************
*** 552,557 ****
--- 557,580 ----
          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
***************
*** 669,674 ****
--- 692,701 ----
          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)
***************
*** 688,693 ****
--- 715,722 ----
        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
***************
*** 704,709 ****
--- 733,755 ----
                
              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 cosmomc/camb/params.ini cosmomc_addon/camb/params.ini
*** cosmomc/camb/params.ini	2008-03-26 15:24:18.000000000 +0100
--- cosmomc_addon/camb/params.ini	2008-07-23 17:41:01.000000000 +0200
***************
*** 65,70 ****
--- 65,77 ----
  
  #Initial power spectrum, amplitude, spectral index and running. Pivot k in Mpc^{-1}.
  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
+ 
  pivot_scalar              = 0.05
  pivot_tensor              = 0.05
  scalar_amp(1)             = 2.3e-9
***************
*** 129,134 ****
--- 136,147 ----
  #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
***************
*** 186,188 ****
--- 199,207 ----
  #interpolation errors may be up to 3%
  #Decrease to speed up non-flat models a bit
  l_sample_boost          = 1
+ 
+ #Use spline interpolation for the Cls
+ do_cls_interpol         = F
+  
+ # Remove undersampling of the Cls: if true lSampleBoost has no effect
+ compute_every_l         = T
diff -r -c -b -B cosmomc/camb/power_tilt.f90 cosmomc_addon/camb/power_tilt.f90
*** cosmomc/camb/power_tilt.f90	2008-03-13 15:40:20.000000000 +0100
--- cosmomc_addon/camb/power_tilt.f90	2008-07-23 17:41:09.000000000 +0200
***************
*** 50,58 ****
--- 50,77 ----
         real(dl) rat(nnmax) !ratio of scalar to tensor initial power spectrum amplitudes
         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
***************
*** 61,66 ****
--- 80,90 ----
   
        public InitialPowerParams, InitialPower_ReadParams, InitializePowers, ScalarPower, TensorPower
        public nnmax,Power_Descript, Power_Name, SetDefPowerParams
+ 
+ !addon
+       public FreeExternalData
+ !end addon
+ 
  !      public 
      contains
         
***************
*** 76,81 ****
--- 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
  
***************
*** 86,91 ****
--- 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
+ 
          if (AParamSet%nn > nnmax) then
             write (*,*) 'To use ',AParamSet%nn,'power spectra you need to increase'
             write (*,*) 'nnmax in power_tilt.f90, currently ',nnmax
***************
*** 96,104 ****
--- 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
***************
*** 119,127 ****
--- 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
  
        
***************
*** 142,148 ****
--- 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)) 
  
         
***************
*** 193,209 ****
            
             InitPower%k_0_scalar = Ini_Read_Double_File(Ini,'pivot_scalar',InitPower%k_0_scalar)
             InitPower%k_0_tensor = Ini_Read_Double_File(Ini,'pivot_tensor',InitPower%k_0_tensor) 
             InitPower%nn = Ini_Read_Int('initial_power_num')
             if (InitPower%nn>nnmax) stop 'Too many initial power spectra - increase nnmax in InitialPower'
             InitPower%rat(:) = 1
             do i=1, InitPower%nn
! 
                InitPower%an(i) = Ini_Read_Double_Array_File(Ini,'scalar_spectral_index', i)
                InitPower%n_run(i) = Ini_Read_Double_Array_File(Ini,'scalar_nrun',i,0._dl) 
                
                if (WantTensors) then
                   InitPower%ant(i) = Ini_Read_Double_Array_File(Ini,'tensor_spectral_index',i)
                   InitPower%rat(i) = Ini_Read_Double_Array_File(Ini,'initial_ratio',i)
                end if              
  
                InitPower%ScalarPowerAmp(i) = Ini_Read_Double_Array_File(Ini,'scalar_amp',i,1._dl) 
--- 429,464 ----
            
             InitPower%k_0_scalar = Ini_Read_Double_File(Ini,'pivot_scalar',InitPower%k_0_scalar)
             InitPower%k_0_tensor = Ini_Read_Double_File(Ini,'pivot_tensor',InitPower%k_0_tensor) 
+ 
+ !addon
+             InitPower%WantExternalScalar = Ini_Read_Logical('use_external_scalar',.false.)
+             InitPower%WantExternalTensor = Ini_Read_Logical('use_external_tensor',.false.)
+ !end addon
             InitPower%nn = Ini_Read_Int('initial_power_num')
             if (InitPower%nn>nnmax) stop 'Too many initial power spectra - increase nnmax in InitialPower'
             InitPower%rat(:) = 1
             do i=1, InitPower%nn
! !addon
!               if (InitPower%WantExternalScalar) then
!                  InitPower%ScalarPowerFileName(i) & 
!                       = Ini_Read_String_Array_File(Ini,'scalar_input_file',i)
!               else
! !end addon
                InitPower%an(i) = Ini_Read_Double_Array_File(Ini,'scalar_spectral_index', i)
                InitPower%n_run(i) = Ini_Read_Double_Array_File(Ini,'scalar_nrun',i,0._dl) 
+               endif
                
                if (WantTensors) then
+ !addon
+                  if (InitPower%WantExternalTensor) then
+                     InitPower%TensorPowerFileName(i) & 
+                          = Ini_Read_String_Array_File(Ini,'tensor_input_file',i)
+                  else
+ !end addon
                   InitPower%ant(i) = Ini_Read_Double_Array_File(Ini,'tensor_spectral_index',i)
                   InitPower%rat(i) = Ini_Read_Double_Array_File(Ini,'initial_ratio',i)
+                  endif
+ 
                end if              
  
                InitPower%ScalarPowerAmp(i) = Ini_Read_Double_Array_File(Ini,'scalar_amp',i,1._dl) 
diff -r -c -b -B cosmomc/camb/subroutines.f90 cosmomc_addon/camb/subroutines.f90
*** cosmomc/camb/subroutines.f90	2008-06-09 19:20:30.000000000 +0200
--- cosmomc_addon/camb/subroutines.f90	2008-07-23 17:41:19.000000000 +0200
***************
*** 1145,1147 ****
--- 1145,1180 ----
  !  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 -r -c -b -B cosmomc/distparams.ini cosmomc_addon/distparams.ini
*** cosmomc/distparams.ini	2008-05-02 15:11:47.000000000 +0200
--- cosmomc_addon/distparams.ini	2008-07-23 17:41:32.000000000 +0200
***************
*** 42,47 ****
--- 42,48 ----
  
  plot_meanlikes = T
  shade_meanlikes = T
+ td_meanlikes = F
  
  # if non-zero, output _thin file, thinned by thin_factor
  thin_factor = 0
diff -r -c -b -B cosmomc/params.ini cosmomc_addon/params.ini
*** cosmomc/params.ini	2008-06-07 16:37:44.000000000 +0200
--- cosmomc_addon/params.ini	2008-07-23 17:41:32.000000000 +0200
***************
*** 127,133 ****
  #If true, generate checkpoint files and terminated runs can be restarted using exactly the same command
  #and chains continued from where they stopped
  #With checkpoint=T note you must delete all chains/file_root.* files if you want new chains with an old file_root
! checkpoint = F
  
  
  #CAMB parameters
--- 127,133 ----
  #If true, generate checkpoint files and terminated runs can be restarted using exactly the same command
  #and chains continued from where they stopped
  #With checkpoint=T note you must delete all chains/file_root.* files if you want new chains with an old file_root
! checkpoint = T
  
  
  #CAMB parameters
***************
*** 143,148 ****
--- 143,151 ----
  #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 -r -c -b -B cosmomc/source/CMB_Cls_simple.f90 cosmomc_addon/source/CMB_Cls_simple.f90
*** cosmomc/source/CMB_Cls_simple.f90	2008-03-13 20:19:20.000000000 +0100
--- cosmomc_addon/source/CMB_Cls_simple.f90	2008-07-23 17:41:38.000000000 +0200
***************
*** 288,293 ****
--- 288,311 ----
     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)
***************
*** 319,327 ****
              P%Transfer%high_precision=.true.
              P%Transfer%kmax=P%Transfer%kmax + 0.2
             end if
!           AccuracyBoost = AccuracyLevel
!           lAccuracyBoost = AccuracyLevel
!           lSampleBoost = AccuracyLevel
            P%AccurateReionization = .true.
          end if
                  
--- 337,348 ----
              P%Transfer%high_precision=.true.
              P%Transfer%kmax=P%Transfer%kmax + 0.2
             end if
! !addon
! !moved up for screen output
! !          AccuracyBoost = AccuracyLevel
! !          lAccuracyBoost = AccuracyLevel
! !          lSampleBoost = AccuracyLevel
! !end addon
            P%AccurateReionization = .true.
          end if
                  
diff -r -c -b -B cosmomc/source/driver.F90 cosmomc_addon/source/driver.F90
*** cosmomc/source/driver.F90	2008-06-10 11:13:47.000000000 +0200
--- cosmomc_addon/source/driver.F90	2008-07-23 17:41:44.000000000 +0200
***************
*** 58,69 ****
  
  #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')
  
--- 58,77 ----
  
  #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')
  
***************
*** 85,90 ****
--- 93,103 ----
  
          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.)
+ !end addon
  
          checkpoint = Ini_Read_Logical('checkpoint',.false.)
          if (checkpoint) flush_write = .true.
diff -r -c -b -B cosmomc/source/GetDist.f90 cosmomc_addon/source/GetDist.f90
*** cosmomc/source/GetDist.f90	2008-06-10 11:13:47.000000000 +0200
--- cosmomc_addon/source/GetDist.f90	2008-07-23 17:42:04.000000000 +0200
***************
*** 85,90 ****
--- 85,93 ----
          logical force_twotail 
          logical smoothing, plot_meanlikes, plot_NDcontours
          logical shade_meanlikes, make_single_samples
+ !addon
+         logical td_meanlikes,volume_plots
+ !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)   
***************
*** 1120,1127 ****
--- 1123,1311 ----
           close(40)
  
           end subroutine DoConvergeTests
+ !addon
+         subroutine Get4DPlotData(j,j2,j3)
+          integer, intent(in) :: j,j2,j3
+          integer i,ix1,ix2,ix3,wx,wy,wz, ix1min,ix1max,ix2min,ix2max,ix3min,ix3max
+          real binweight, dist1, dist2, dist3, norm, maxbin
+          real try_b, try_t,try_sum, try_last
+          character(LEN=120) :: plotfile, filename,numstr
+          real, dimension(:,:,:), allocatable :: TheBins3, bins3D, bin3Dlikes, bin3Dmax
+ 
+                 
+          ix1min = 0
+          ix1max = 0
+          ix2min = 0
+          ix2max = 0
+          ix3min = 0
+          ix3max = 0
+ 
+ 
+          do i = 0,nrows-1
+             ix1=nint((coldata(colix(j),i)-center(j))/width(j))
+             ix2=nint((coldata(colix(j2),i)-center(j2))/width(j2))
+             ix3=nint((coldata(colix(j3),i)-center(j3))/width(j3))
+             ix1min = min(ix1min,min(ix_min(j),ix1-2))
+             ix1max = max(ix1max,max(ix_max(j),ix1+2))
+             ix2min = min(ix2min,min(ix_min(j2),ix2-2))
+             ix2max = max(ix2max,max(ix_max(j2),ix2+2))
+             ix3min = min(ix3min,min(ix_min(j3),ix3-2))
+             ix3max = max(ix3max,max(ix_max(j3),ix3+2))           
+          enddo
+ 
+ !         print *,'bounds 1',ix1min,ix1max
+ !         print *,'bounds 2',ix2min,ix2max
+ !         print *,'bounds 3',ix3min,ix3max
+ !         read(*,*)
+          allocate(bins3D(ix1min:ix1max,ix2min:ix2max,ix3min:ix3max))
+          allocate(bin3Dlikes(ix1min:ix1max,ix2min:ix2max,ix3min:ix3max))
+          allocate(bin3Dmax(ix1min:ix1max,ix2min:ix2max,ix3min:ix3max))
+ 
+          bins3D = 0.
+ 
+ 
+          bin3Dlikes = 0.
+          bin3Dmax = 1e30         
+          
+ !         print *,'jj2j3',j,j2,j3
+ !         print *,'center',center(j),center(j2),center(j3)
+ !         print *,'width',width(j),width(j2),width(j3)
+ !         print *,'ixmin',ix_min(j),ix_min(j2),ix_min(j3)
+ !         print *,'ixmax',ix_max(j),ix_max(j2),ix_max(j3)
+                 
+ 
+          do i = 0, nrows-1
+             ix1=nint((coldata(colix(j),i)-center(j))/width(j))
+             ix2=nint((coldata(colix(j2),i)-center(j2))/width(j2))
+             ix3=nint((coldata(colix(j3),i)-center(j3))/width(j3))
           
  
+             if (smoothing) then
+                do wx = max(ix_min(j),ix1-2),min(ix_max(j),ix1+2)
+                   dist1 = (coldata(colix(j),i) - (center(j) + real(wx)*width(j)))**2/width(j)**2
+ 
+                   do wy = max(ix_min(j2),ix2-2),min(ix_max(j2),ix2+2)
+                      dist2 = (coldata(colix(j2),i) - (center(j2) + real(wy)*width(j2)))**2/width(j2)**2
+ 
+                      do wz =  max(ix_min(j3),ix3-2),min(ix_max(j3),ix3+2)
+                         dist3 = (coldata(colix(j3),i) - (center(j3) + real(wz)*width(j3)))**2/width(j3)**2 
+ 
+                         binweight = coldata(1,i)*exp( - (dist1 + dist2 + dist3)/2.)
+                         bins3D(wx,wy,wz) = bins3D(wx,wy,wz) + binweight
+                         bin3Dlikes(wx,wy,wz) = bin3Dlikes(wx,wy,wz) + binweight*exp(meanlike - coldata(2,i))
+ 
+                      end do
+                   end do
+                enddo
+                
+             else
+                bins3D(ix1,ix2,ix3) = bins3D(ix1,ix2,ix3) + coldata(1,i) 
+                bin3Dlikes(ix1,ix2,ix3) = bin3Dlikes(ix1,ix2,ix3) + coldata(1,i)*exp(meanlike-coldata(2,i))
+             end if
+             
+             bin3Dmax(ix1,ix2,ix3) = min(bin3Dmax(ix1,ix2,ix3),real(coldata(2,i)))
+             
+          end do
+          
+          
+          do ix1=ix_min(j),ix_max(j)
+             do ix2 =ix_min(j2),ix_max(j2)
+                do ix3 = ix_min(j3),ix_max(j3)
+                   if (bins3D(ix1,ix2,ix3) >0) bin3Dlikes(ix1,ix2,ix3) = bin3Dlikes(ix1,ix2,ix3)/bins3D(ix1,ix2,ix3)
+                   if (plot_NDcontours) then
+                   !Map maximum likelihood in each bin into a significance value from the full N-D distribution
+                      if (bin3DMax(ix1,ix2,ix3) == 1e30) then
+                         bin3DMax(ix1,ix2,ix3) = 0
+                      else
+                         bin3DMax(ix1,ix2,ix3) = &
+                              sum(coldata(1,0:nrows-1), mask = coldata(2,0:nrows-1) > bin3DMax(ix1,ix2,ix3))/numsamp
+                      end if
+                   end if
+                end do
+             end do
+          enddo
+          
+          if (has_limits(colix(j)) .or. has_limits(colix(j2))) then
+ 
+         !Fix up underweighting near edges. Note this makes edge pixels noisier.
+             do ix1 = ix_min(j), ix_max(j)
+                  do ix2 = ix_min(j2),ix_max(j2)
+                     do ix3 = ix_min(j3),ix_max(j3)                       
+                     if (ix1 ==ix_min(j) .and. has_limits_bot(colix(j))) bins3D(ix1,ix2,ix3) = bins3D(ix1,ix2,ix3)*2
+                     if (ix1 ==ix_max(j) .and. has_limits_top(colix(j))) bins3D(ix1,ix2,ix3) = bins3D(ix1,ix2,ix3)*2
+                     if (ix2 ==ix_min(j2).and. has_limits_bot(colix(j2))) bins3D(ix1,ix2,ix3) = bins3D(ix1,ix2,ix3)*2
+                     if (ix2 ==ix_max(j2).and. has_limits_top(colix(j2))) bins3D(ix1,ix2,ix3) = bins3D(ix1,ix2,ix3)*2
+                     if (ix3 ==ix_min(j3).and. has_limits_bot(colix(j3))) bins3D(ix1,ix2,ix3) = bins3D(ix1,ix2,ix3)*2
+                     if (ix3 ==ix_max(j3).and. has_limits_top(colix(j3))) bins3D(ix1,ix2,ix3) = bins3D(ix1,ix2,ix3)*2
+ 
+                     if (ix1 ==ix_min(j)+1.and. has_limits_bot(colix(j))) bins3D(ix1,ix2,ix3) = bins3D(ix1,ix2,ix3)/0.84
+                     if (ix1 ==ix_max(j)-1.and. has_limits_top(colix(j))) bins3D(ix1,ix2,ix3) = bins3D(ix1,ix2,ix3)/0.84
+                     if (ix2 ==ix_min(j2)+1.and. has_limits_bot(colix(j2))) bins3D(ix1,ix2,ix3) = bins3D(ix1,ix2,ix3)/0.84
+                     if (ix2 ==ix_max(j2)-1.and. has_limits_top(colix(j2))) bins3D(ix1,ix2,ix3) = bins3D(ix1,ix2,ix3)/0.84
+                     if (ix3 ==ix_min(j3)+1.and. has_limits_bot(colix(j3))) bins3D(ix1,ix2,ix3) = bins3D(ix1,ix2,ix3)/0.84
+                     if (ix3 ==ix_max(j3)-1.and. has_limits_top(colix(j3))) bins3D(ix1,ix2,ix3) = bins3D(ix1,ix2,ix3)/0.84
+                  end do
+               end do
+            enddo
+         end if
+         
+ ! Get contour containing contours(:) of the probability
+ 
+              allocate(TheBins3(ix_max(j)-ix_min(j)+1,ix_max(j2)-ix_min(j2)+1,ix_max(j3)-ix_min(j3)+1))
+              TheBins3 = bins3D(ix_min(j):ix_max(j),ix_min(j2):ix_max(j2),ix_min(j3):ix_max(j3))
+                 
+              norm = sum(TheBins3)
+              
+              do ix1 = 1, num_contours
+              try_t = maxval(TheBins3)
+              try_b = 0
+              try_last = -1
+              do
+                try_sum = sum(TheBins3,mask = TheBins3 < (try_b + try_t)/2) 
+                if (try_sum > (1-contours(ix1))*norm) then
+                   try_t = (try_b+try_t)/2
+                else
+                   try_b = (try_b+try_t)/2
+                end if
+                if (try_sum == try_last) exit
+                try_last = try_sum
+          
+              end do
+              contour_levels(ix1) = (try_b+try_t)/2
+ 
+              end do
+              deallocate(TheBins3)
+             
+              
+              plotfile = numcat(trim(numcat(trim(numcat(trim(rootname)//'_4D_',colix(j)-2))//'_',colix(j2)-2))//'_',colix(j3)-2)
+              filename = 'plot_data/'//trim(plotfile)
+              open(unit=49,file=filename,form='formatted',status='replace')
+              maxbin = maxval(bin3Dlikes(ix_min(j):ix_max(j),ix_min(j2):ix_max(j2),ix_min(j3):ix_max(j3)))
+              do ix1 = ix_min(j), ix_max(j)
+                 do ix2 = ix_min(j2), ix_max(j2)
+                    do ix3 = ix_min(j3), ix_max(j3)
+                       write (49,'(6(E16.7))') center(j) + ix1*width(j), center(j2) + ix2*width(j2) &
+                            , center(j3) + ix3*width(j3) , bins3D(ix1,ix2,ix3),bin3Dmax(ix1,ix2,ix3) &
+                            , bin3Dlikes(ix1,ix2,ix3)/maxbin
+                    end do
+                 enddo
+              enddo
+              
+              close(49)
+ 
+ 
+ 
+              open(unit=49,file=trim(filename)//'_cont',form='formatted',status='replace')
+             write(numstr,*) contour_levels(1:num_contours)
+             if (num_contours==1) numstr = trim(numstr)//' '//trim(numstr)
+             write (49,*) trim(numstr)
+             close(49)
+ 
+ 
+              deallocate(bins3D,bin3Dmax,bin3Dlikes)
+      
+         end subroutine Get4DPlotData
+ !end addon
  
          subroutine Get2DPlotData(j,j2)
           integer, intent(in) :: j,j2
***************
*** 1336,1348 ****
--- 1520,1547 ----
   !             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('"//trim(plot_data_dir)//"','" // 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
                
             
***************
*** 1534,1540 ****
  
          real thin_cool
          logical no_tests, auto_label
! 
          logical :: single_column_chain_files, samples_are_chains
          integer sz
          !for single_colum_chain_files    
--- 1733,1741 ----
  
          real thin_cool
          logical no_tests, auto_label
! !addon
!         integer :: w,j3,i2
! !end addon
          logical :: single_column_chain_files, samples_are_chains
          integer sz
          !for single_colum_chain_files    
***************
*** 1697,1702 ****
--- 1898,1907 ----
            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.)  
+     volume_plots = Ini_Read_Logical('volume_plots',.false.)
+ !end addon
  
      matlab_col = Ini_Read_String('matlab_colscheme')
      smoothing = Ini_Read_Logical('smoothing',.true.)
***************
*** 2342,2348 ****
--- 2547,2560 ----
                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
***************
*** 2540,2545 ****
--- 2752,2784 ----
            end if
      
     
+ !addon
+           if (num_3D_plots /=0 .and. .not. no_plots .and. volume_plots) then
+              write (*,*) 'Producing ',num_3D_plots,' 3D volume plots'
+              
+              do i=1, num_3D_plots
+                 read(plot_3D(i),*) ix1,ix2,ix3
+ !                print *,'p1 p2 p3',ix1, ix2, ix3
+                 j=0
+                 j2=0
+                 j3=0
+                 do w = 1,num_vars
+ !                   print *,'colix',colix(w)
+                    if (colix(w)-2.eq.ix1) j = w
+                    if (colix(w)-2.eq.ix2) j2 = w
+                    if (colix(w)-2.eq.ix3) j3 = w
+                 enddo
+ !                print *,'jj2j3',j,j2,j3
+                 if (j*j2*j3.eq.0) then
+                    print *
+                    print *,ix1,ix2,ix3,' not in 1D: skipped'
+                    print *
+                 else
+                    call Get4DPlotData(j,j2,j3)  
+                 endif
+              enddo
+           endif
+ !end addon
  
  !write out stats
  !Marginalized
diff -r -c -b -B cosmomc/source/Makefile cosmomc_addon/source/Makefile
*** cosmomc/source/Makefile	2008-05-06 22:43:49.000000000 +0200
--- cosmomc_addon/source/Makefile	2008-05-12 10:14:19.000000000 +0200
***************
*** 1,109 ****
! #You may need to edit the library paths for MKL for Intel
! #Beware of using optmizations that lose accuracy - may give errors when running
  
! #Edit for CFITSIO directories
! cfitsio = /usr/local/cfitsio/intel10/64/3.040
! #/usr/local/Cluster-Users/cpac/cmb/2.1.0/cfitsio
! 
! WMAP = /home/aml1005/WMAP5/likelihood_v3
! 
! IFLAG = -I
! INCLUDE=
! 
! #Intel MPI (assuming mpif77 set to point to ifort)
! #change -lmkl_ia32 to -lmkl_p3 for MKL versions earlier than 6 (6 needed for ifc 8+) 
! F90C     = mpif90
! FFLAGS = -O2 -ip -W0 -WB -openmp -fpp -DMPI -vec_report0
! LAPACKL = -L/usr/local/intel/mkl/9.1/lib/em64t -lmkl_lapack -lmkl_em64t -lguide -lpthread
! 
! #HPCF settings. Use Inteal 9 or 10.1, not 10.0
! #F90C     = mpif90
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp -DMPI -vec_report0
! #LAPACKL = -L/usr/local/Cluster-Apps/intel/mkl/9.0.018/lib/em64t -lmkl_lapack -lmkl -lguide -lpthread
! #INCLUDE=
! 
! # Use with Intel 9 on the Altix (eg. COSMOS) [this is old]
! # (do "module load icomp90" before compiling, or add to .bashrc)
! #F90C = ifort
! #FFLAGS = -O3 -w -fpp2 -DMPI 
! #LAPACKL = -lscs -lguide -lmpi -ldl
! 
! #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
! 
! #Intel fortran 8 on COSMOS
! #F90C     = ifort
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp2 -lmpi -lscs  -ldl -DMPI
! #LAPACKL = -L/opt/intel/mkl80/lib/64 -lmkl_lapack -lmkl_ipf -limf
! 
! #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
! #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
! 
! #G95; make sure LAPACK and MPI libs also compiled with g95
! #F90C     = mpif90
! #FFLAGS = -O2 -DMPI
! #LAPACKL = /LAPACK/lapack_LINUX.a /LAPACK/blas_LINUX.a
! 
! #GFortran: if pre v4.3 add -D__GFORTRAN__
! #F90C     = gfc
! #FFLAGS = -O2 -ffree-form 
! #LAPACKL = -Wl,-framework -Wl,accelerate
! 
! #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 
! #IFLAG = -p
! 
! #NAGF95, single processor:
! #F90C     = f95
! #FFLAGS = -DNAGF95 -O3
! #LAPACKL = -llapack -lblas -lg2c
! 
! #PGF90
! #F90C = pgf90
! #FFLAGS = -O2 -DESCAPEBACKSLASH
! #LAPACKL = -llapack -lblas
! 
! 
! #Sun, single processor:
! #F90C     = f90
! #FFLAGS = -fast -ftrap=%none
! #LAPACKL = -lsunperf -lfsu
! #LAPACKL = -lsunperf -lfsu -lsocket -lm
! #IFLAG = -M
! 
! #Sun MPI
! #F90C = mpf90
! #FFLAGS =  -O4 -openmp -ftrap=%none -dalign -DMPI
! #LAPACKL = -lsunperf -lfsu -lmpi_mt
! #IFLAG = -M
! 
! #Sun parallel enterprise:
! #F90C     = f95
! #FFLAGS =  -O4 -xarch=native64 -openmp -ftrap=%none
! #LAPACKL = -lsunperf -lfsu
! #IFLAG = -M
! 
! 
! #IBM XL Fortran, multi-processor (run "module load lapack" then run "gmake")
! # See also http://cosmocoffee.info/viewtopic.php?t=326
! #F90C     = xlf90_r $(LAPACK)
! #FFLAGS  = -WF,-DIBMXL -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
! #INCLUDE = -lessl
! #LAPACKL =
  
  PROPOSE = propose.o
  CLSFILE = CMB_Cls_simple.o
--- 1,57 ----
! # >>> 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
! WMAPFLAGS= $(FFLAGS)
! LAPACKL = -lcomplib.sgimath_mp
! INCLUDE = -I../camb
! CFITSIODIR = /people/ringeval/usr
! endif
! 
! ifeq ($(ext),Lin)
! F90C=ifort
! FFLAGS= -O -openmp -fpp
! WMAPFLAGS= -O -fpp
! LAPACKL = -L/opt/intel/mkl/9.1/lib/64 -lmkl -lmkl_lapack
! INCLUDE = -I../camb 
! CFITSIODIR = /opt/intel/cfitsio
! endif 
! 
! ifeq ($(ext),OSF)
! F90C=f90
! FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
! WMAPFLAGS= $(FFLAGS)
! LAPACKL  = -lcxml
! INCLUDE = -I../camb 
! CFITSIODIR = 
! endif
! 
! ifeq ($(ext),Sun)
! F90C=f90
! FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
! WMAPFLAGS= $(FFLAGS)
! LAPACKL = -lsunperf -lfsu
! INCLUDE = -I../camb -M../camb
! CFITSIODIR = 
! endif
! 
! ifeq ($(ext),AIX)
! F90C    = mpxlf90_r
! FFLAGS  = -O4 -WF,-DIBMXL,-DMPI -qstrict -qsmp=omp -qmaxmem=-1 -qsuffix=f=f90:cpp=F90
! WMAPFLAGS= $(FFLAGS)
! LAPACKL = -lessl
! INCLUDE = -I../camb
! CFITSIODIR = 
! endif
! 
! 
! WMAPDIR = ../WMAP
  
  PROPOSE = propose.o
  CLSFILE = CMB_Cls_simple.o
***************
*** 111,144 ****
  #Can use params_H if you prefer more generic parameters
  PARAMETERIZATION = params_CMB.o
  
! WMAP_inc = $(IFLAG)$(cfitsio)/include $(IFLAG)$(WMAP)
! WMAP_lib = -L$(cfitsio)/lib -L$(WMAP) -lcfitsio
  
! CLSLIB  = $(WMAP_lib) -L../camb -lcamb $(LAPACKL) 
  
! F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) $(WMAP_inc) $(IFLAG)../camb $(INCLUDE)
  
  DISTFILES = utils.o Matrix_utils.o settings.o GetDist.o
  
! WMAPOBJS = $(WMAP)/read_archive_map.o \
! 	$(WMAP)/read_fits.o \
! 	$(WMAP)/WMAP_5yr_options.o \
! 	$(WMAP)/WMAP_5yr_util.o \
! 	$(WMAP)/WMAP_5yr_tt_pixlike.o \
! 	$(WMAP)/WMAP_5yr_teeebb_pixlike.o \
! 	$(WMAP)/WMAP_5yr_likelihood.o \
! 	$(WMAP)/WMAP_5yr_gibbs.o \
! 	$(WMAP)/WMAP_5yr_tt_beam_ptsrc_chisq.o \
! 	$(WMAP)/br_mod_dist.o 
  
  OBJFILES= $(WMAPOBJS) utils.o Matrix_utils.o settings.o cmbtypes.o Planck_like.o \
  	cmbdata.o WeakLen.o mpk.o supernovae.o SDSSLy-a-v3.o\
  	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
! default: cosmomc
  
! all : cosmomc getdist
  
  settings.o: utils.o
  cmbtypes.o: settings.o
--- 59,94 ----
  #Can use params_H if you prefer more generic parameters
  PARAMETERIZATION = params_CMB.o
  
! WMAPINCLUDE = -I$(CFITSIODIR)/include
  
! CLSLIB = -L$(CFITSIODIR)/lib -L../camb -lcamb -lcfitsio $(LAPACKL)
  
! F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) $(INCLUDE)
  
  DISTFILES = utils.o Matrix_utils.o settings.o GetDist.o
  
! WMAPOBJS = read_archive_map.o \
! 	read_fits.o \
! 	healpix_types.o \
! 	br_mod_dist.o \
! 	WMAP_5yr_options.o \
! 	WMAP_5yr_util.o \
! 	WMAP_5yr_gibbs.o \
! 	WMAP_5yr_tt_pixlike.o \
! 	WMAP_5yr_tt_beam_ptsrc_chisq.o \
! 	WMAP_5yr_teeebb_pixlike.o \
! 	WMAP_5yr_tetbeebbeb_pixlike.o \
! 	WMAP_5yr_likelihood.o
! 
  
  OBJFILES= $(WMAPOBJS) utils.o Matrix_utils.o settings.o cmbtypes.o Planck_like.o \
  	cmbdata.o WeakLen.o mpk.o supernovae.o SDSSLy-a-v3.o\
  	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
! default: cosmomc.$(ext)
  
! all : cosmomc.$(ext) getdist.$(ext)
  
  settings.o: utils.o
  cmbtypes.o: settings.o
***************
*** 162,167 ****
--- 112,124 ----
  .f.o:
  	f77 $(F90FLAGS) -c $<
  
+ %.o: $(WMAPDIR)/%.f90
+ 	$(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $<
+ 
+ %.o: $(WMAPDIR)/%.F90
+ 	$(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $<
+ 
+ 
  %.o: %.f90
  	$(F90C) $(F90FLAGS) -c $*.f90
  
***************
*** 169,180 ****
  	$(F90C) $(F90FLAGS) -c $*.F90
  
  
! cosmomc: $(OBJFILES) ../camb/libcamb.a 	
! 	$(F90C) -o ../cosmomc $(OBJFILES) $(CLSLIB) $(F90FLAGS)
  
  
  clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist: $(DISTFILES)
! 	$(F90C) -o ../getdist $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
--- 126,137 ----
  	$(F90C) $(F90FLAGS) -c $*.F90
  
  
! cosmomc.$(ext): $(OBJFILES) ../camb/libcamb.a 	
! 	$(F90C) -o ../$@ $(OBJFILES) $(CLSLIB) $(F90FLAGS)
  
  
  clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist.$(ext): $(DISTFILES)
! 	$(F90C) -o ../$@ $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
diff -r -c -b -B cosmomc/source/MCMC.f90 cosmomc_addon/source/MCMC.f90
*** cosmomc/source/MCMC.f90	2008-03-05 11:39:29.000000000 +0100
--- cosmomc_addon/source/MCMC.f90	2008-07-23 17:42:09.000000000 +0200
***************
*** 214,220 ****
    
      Like = GetLogLike(grid(r))
     
!     if (Feedback > 1) write (*,*) r, 'Likelihood: ', Like, 'Current Like:', CurLike
  
      if ((Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike)) then
      !Accept
--- 214,220 ----
    
      Like = GetLogLike(grid(r))
     
!     if (Feedback > 1) write (*,*) r, 'Likelihood: ', Like, 'Current Like:', CurLike, 'Rank:',MPIRank
  
      if ((Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike)) then
      !Accept
***************
*** 499,505 ****
             output_lines = output_lines +1
             call WriteParams(CurParams, real(mult), CurLike)
            end if
!           if (Feedback > 1) write (*,*) instance, 'Slicing, Current Like:', CurLike
            mult = 1
            if (num_slow /=0) call SliceSampleSlowParam(CurParams, CurLike) 
            if (num_fast /=0) call SliceSampleFastParams(CurParams, CurLike) 
--- 499,505 ----
             output_lines = output_lines +1
             call WriteParams(CurParams, real(mult), CurLike)
            end if
!           if (Feedback > 1) write (*,*) instance, 'Slicing, Current Like:', CurLike, 'Rank:',MPIRank
            mult = 1
            if (num_slow /=0) call SliceSampleSlowParam(CurParams, CurLike) 
            if (num_fast /=0) call SliceSampleFastParams(CurParams, CurLike) 
***************
*** 534,540 ****
  
      Like = GetLogLike(Trial) 
  
!     if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:', CurLike
  
  !!!    accpt = (Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike) 
               !Include the min() so that compilers not doing optimal compilation don't complain
--- 534,540 ----
  
      Like = GetLogLike(Trial) 
  
!     if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:', CurLike, 'Rank:',MPIRank
  
  !!!    accpt = (Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike) 
               !Include the min() so that compilers not doing optimal compilation don't complain
diff -r -c -b -B cosmomc/source/settings.f90 cosmomc_addon/source/settings.f90
*** cosmomc/source/settings.f90	2008-06-10 11:13:47.000000000 +0200
--- cosmomc_addon/source/settings.f90	2008-07-23 17:42:14.000000000 +0200
***************
*** 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)
  
***************
*** 52,57 ****
--- 57,65 ----
  
    integer :: num_threads = 0
    integer :: instance = 0
+ !addon
+   integer :: instance_shift = 0
+ !end addon
    integer :: MPIchains = 1, MPIrank = 0
  
    logical :: Use_LSS = .true.
