diff --git a/README b/README
new file mode 100644
index 0000000..ec492a8
--- /dev/null
+++ b/README
@@ -0,0 +1,70 @@
+**********************************************************************
+Add-on switches to CAMB and cosmomc:
+- compute all multipoles
+- remove interpolation
+- new fine tuned accuracy parameters
+- use external files for the primordial spectra
+- second argument in mpirun sets the file suffix
+- mean likelihoods plotted in 3D as translucent surfaces
+**********************************************************************
+
+-These patches modify the cosmomc or CAMB codes and therefore may have
+bugs that are not present in the originals. The authors of the
+original versions should not be bothered with problems resulting in
+the use of these patches. Feel free to insult me, or at least contact
+me, in case of bugs!
+
+-Apply the patch in the cosmomc/ directory:
+patch -p1 < purpose_MMYY.patch
+
+*Modifications are traced by the "addon" in the source files.
+
+*Please note that the genuine Makefiles have been heavily changed and
+may require edition to work on your system. They assume the existence
+of a directory "cosmomc/WMAP" which points to the WMAP team likelihood
+code directory. The WMAP data are also assumed to reside in the
+"cosmomc/data" directory.
+
+
+*Modified files and settings 
+
+- camb/modules.f90 
+New switches "DoClsInterpol" and "ComputeEveryl".
+
+- camb/power_tilt.f90
+Can read from external files the scalar and tensor power spectra. Uses
+a spline interpolation in ln(k).
+
+- camb/inidriver.F90
+New keys for the above material.
+
+- camb/params.ini
+New keys for the above material.
+
+-camb/subroutines.f90
+Added subroutine splinterpol
+
+- cosmomc/CMB_Cls_simple.f90
+Added screen output to check camb accuracy default values. New
+switches lSampleLevel, lAccuracyLevel, CambComputeEveryl to control
+the camb accuracy from cosmomc.
+
+- cosmomc/GetDist.f90
+3d transparent plots of the mean likelihoods. Try colormap('copper').
+
+- cosmomc/driver.F90
+Second instance in mpirun sets the file name suffix.
+
+-cosmomc/MCMC.f90
+Add MpiRank output in the logs
+
+- cosmomc/settings.f90
+Define lAccuracyLevel,lSampleLevel, CambComputeEveryl, instance_shift.
+
+- cosmomc/params.ini
+New keys for the above material.
+
+- cosmomc/distparams.ini
+New keys for the above material.
+
+
diff --git a/camb/Makefile b/camb/Makefile
index 5c09a67..1aa151f 100644
--- a/camb/Makefile
+++ b/camb/Makefile
@@ -1,83 +1,115 @@
-#CAMB Makefile
+# >>> 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
+LFLAGS= 
+endif
+
+ifeq ($(ext),Lin)
+F90C=gfortran
+FFLAGS= -O -fopenmp #-DMPI
+LFLAGS= 
+endif 
 
-#Set FISHER=Y to compile bispectrum fisher matrix code
-FISHER=
+ifeq ($(ext),OSF)
+F90C=f90
+FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
+LFLAGS= 
+endif
 
-#Edit for your compiler
-#Note there are many old ifc versions, some of which behave oddly
 
+ifeq ($(ext),Sun)
+F90C=f90
+FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
+LFLAGS= 
+endif
 
-#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 -fast -W0 -WB -fpp2 -vec_report0
-ifneq ($(FISHER),)
-FFLAGS += -mkl
+ifeq ($(ext),AIX)
+F90C=xlf90_r
+FFLAGS= -O4 -qsmp=omp -qmaxmem=-1 -qstrict -qfree=f90 -qsuffix=f=f90:cpp=F90
+LFLAGS= 
 endif
 
-#Gfortran compiler:
-#The options here work in v4.5, delete from RHS in earlier versions (15% slower)
-#if pre v4.3 add -D__GFORTRAN__
-#With v4.6+ try -Ofast -march=native -fopenmp
-#On my machine v4.5 is about 20% slower than ifort
-#F90C     = gfortran
-#FFLAGS =  -O3 -fopenmp -ffast-math -march=native -funroll-loops
-
-
-#Old 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
-
-#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 -Mpreprocess
-
-#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
-
-#Settings for building camb_fits
-#Location of FITSIO and name of library
-FITSDIR       = /home/cpac/cpac-tools/lib
-FITSLIB       = cfitsio
-#Location of HEALPIX for building camb_fits
-HEALPIXDIR    = /home/cpac/cpac-tools/healpix
-
-ifneq ($(FISHER),)
+
+#Files containing evolution equations initial power spectrum module
+EQUATIONS     = equations
+POWERSPECTRUM = power_tilt
+REIONIZATION = reionization
+RECOMBINATION = recfast
+BISPECTRUM = SeparableBispectrum
+DENABLE_FISHER =
+
+#Module doing non-linear scaling
+NONLINEAR     = halofit
+
+#Driver program
+DRIVER        = inidriver.F90
+#DRIVER        = sigma8.f90
+#DRIVER        = tester.f90
+
+CAMBLIB       = libcamb.a
+
+
+ifneq ($(DENABLE_FISHER),)
 FFLAGS += -DFISHER
+LFLAGS += -llapack -lblas
 EXTCAMBFILES = Matrix_utils.o
 else
 EXTCAMBFILES =
 endif
 
-include ./Makefile_main
+HEALPIXDIR=
+
+#Shouldn't need to change anything else...
+
+F90FLAGS      = $(FFLAGS)
+HEALPIXLD     = -L$(HEALPIXDIR)/lib -lhealpix 
+
+CAMBOBJ       = constants.o utils.o $(EXTCAMBFILES) subroutines.o inifile.o $(POWERSPECTRUM).o $(RECOMBINATION).o $(REIONIZATION).o modules.o \
+        bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o $(BISPECTRUM).o cmbmain.o camb.o
+
+default: camb.$(ext)
+
+all: camb.$(ext) $(CAMBLIB)
+
+subroutines.o: constants.o utils.o
+$(POWERSPECTRUM).o: subroutines.o  inifile.o
+$(RECOMBINATION).o: subroutines.o inifile.o
+$(REIONIZATION).o: constants.o inifile.o
+modules.o: $(REIONIZATION).o $(POWERSPECTRUM).o $(RECOMBINATION).o
+bessels.o: modules.o
+$(EQUATIONS).o: bessels.o
+$(NONLINEAR).o:  modules.o
+lensing.o: bessels.o
+$(BISPECTRUM).o: lensing.o modules.o
+cmbmain.o: lensing.o $(NONLINEAR).o $(EQUATIONS).o $(BISPECTRUM).o
+camb.o: cmbmain.o
+Matrix_utils.o: utils.o
+
+
+camb.$(ext): $(CAMBOBJ) $(DRIVER)
+	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) $(LFLAGS) -o $@
+
+$(CAMBLIB): $(CAMBOBJ)
+	ar -r $@ $?
+
+camb_fits: writefits.f90 $(CAMBOBJ) $(DRIVER)
+	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include tils.o $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) $(LFLAGS) -DWRITE_FITS -o $@
+
+%.o: %.f90
+	$(F90C) $(F90FLAGS) -c $*.f90
+
+utils.o:
+	$(F90C) $(F90FLAGS) -c utils.F90	
+
+$(BISPECTRUM).o:
+	$(F90C) $(F90FLAGS) -c $(BISPECTRUM).F90	
+
+Matrix_utils.o:
+	$(F90C) $(F90FLAGS) -c Matrix_utils.F90	
+
+clean:
+	-rm -f *.o *.a *.d core *.mod
diff --git a/camb/inidriver.F90 b/camb/inidriver.F90
index ca231b5..8720dc2 100644
--- a/camb/inidriver.F90
+++ b/camb/inidriver.F90
@@ -266,6 +266,16 @@
        else
         lSampleBoost   = Ini_Read_Double('l_sample_boost',lSampleBoost)
        end if
+!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)
        end if
@@ -320,6 +330,10 @@
 
         call CAMB_cleanup             
   
+        if (P%InitPower%WantExternalScalar.or.P%InitPower%WantExternalTensor) &
+             call FreeExternalData(P%InitPower)
+
+
         end program driver
 
 
diff --git a/camb/modules.f90 b/camb/modules.f90
index 68c74b7..9baf89d 100644
--- a/camb/modules.f90
+++ b/camb/modules.f90
@@ -212,7 +212,11 @@
 
 !     lmax is max possible number of l's evaluated
       integer, parameter :: lmax_arr = l0max 
- 
+
+!addon
+      logical :: DoClsInterpol=.true.
+      logical :: ComputeEveryl=.true.
+!end addon 
         contains
       
 
@@ -574,6 +578,24 @@
         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       
        
         if (lSampleBoost >=50) then
@@ -714,6 +736,10 @@
         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)
@@ -733,6 +759,8 @@
       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
@@ -749,6 +777,23 @@
               
             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 --git a/camb/params.ini b/camb/params.ini
index 95893cc..a1b1552 100644
--- a/camb/params.ini
+++ b/camb/params.ini
@@ -67,6 +67,13 @@ nu_mass_fractions = 1
 
 #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.1e-9
@@ -216,3 +223,9 @@ l_accuracy_boost        = 1
 #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 --git a/camb/power_tilt.f90 b/camb/power_tilt.f90
index c23155f..53066d5 100644
--- a/camb/power_tilt.f90
+++ b/camb/power_tilt.f90
@@ -50,9 +50,28 @@
        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,6 +80,11 @@
  
       public InitialPowerParams, InitialPower_ReadParams, InitializePowers, ScalarPower, TensorPower
       public nnmax,Power_Descript, Power_Name, SetDefPowerParams
+
+!addon
+      public FreeExternalData
+!end addon
+
 !      public 
     contains
        
@@ -76,6 +100,11 @@
          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,6 +115,14 @@
 
         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,9 +133,163 @@
 
 !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,10 +310,30 @@
         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)   
-     
-!         ScalarPower = ScalarPower * (1 + 0.1*cos( lnrat*30 ) )
+       endif
+!end addon
  
       end function ScalarPower
 
@@ -144,7 +355,30 @@
    
         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)) 
 
        
@@ -195,17 +429,36 @@
           
            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 --git a/camb/subroutines.f90 b/camb/subroutines.f90
index b642f34..e453137 100644
--- a/camb/subroutines.f90
+++ b/camb/subroutines.f90
@@ -1131,3 +1131,36 @@
 !  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 --git a/distparams.ini b/distparams.ini
index 95e67d5..2c852c2 100644
--- a/distparams.ini
+++ b/distparams.ini
@@ -49,6 +49,9 @@ compare1 = basic6_cmb
 
 plot_meanlikes = T
 shade_meanlikes = T
+plot_type = stairs
+td_meanlikes = F
+volume_plots = F
 
 # if non-zero, output _thin file, thinned by thin_factor
 thin_factor = 0
diff --git a/params.ini b/params.ini
index 4d574d6..46bff78 100644
--- a/params.ini
+++ b/params.ini
@@ -140,7 +140,7 @@ rand_seed =
 #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
+checkpoint = T
 
 #whether to stop on CAMB error, or continue ignoring point
 stop_on_error=  T
@@ -166,6 +166,9 @@ high_accuracy_default = F
 #(default is about 0.3%, 0.1% if high_accuracy_default =T)
 #accuracy_level can be increased to check stability/higher accuracy, but inefficient
 accuracy_level = 1
+l_accuracy_level = 1
+l_sample_level = 1
+camb_compute_every_l = F
 
 #If action = 1
 redo_likelihoods = T
diff --git a/source/CMB_Cls_simple.f90 b/source/CMB_Cls_simple.f90
index db9d905..f30cbd3 100644
--- a/source/CMB_Cls_simple.f90
+++ b/source/CMB_Cls_simple.f90
@@ -404,6 +404,25 @@ contains
    integer zix
    real redshifts(matter_power_lnzsteps)
 
+!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)
@@ -436,9 +455,12 @@ contains
             P%Transfer%high_precision=.true.
             P%Transfer%kmax=P%Transfer%kmax + 0.2
            end if
-          AccuracyBoost = AccuracyLevel
-          lAccuracyBoost = AccuracyLevel
-          lSampleBoost = AccuracyLevel
+!addon
+!moved up for screen output
+!          AccuracyBoost = AccuracyLevel
+!          lAccuracyBoost = AccuracyLevel
+!          lSampleBoost = AccuracyLevel
+!end addon
           P%AccurateReionization = .true.
         end if
                 
diff --git a/source/GetDist.f90 b/source/GetDist.f90
index b4094ce..f487c65 100644
--- a/source/GetDist.f90
+++ b/source/GetDist.f90
@@ -94,6 +94,10 @@ module MCSamples
         logical force_twotail 
         logical smoothing, plot_meanlikes, plot_NDcontours
         logical shade_meanlikes, make_single_samples
+!addon
+        logical td_meanlikes,volume_plots
+        character(len=6) :: plottype
+!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)   
@@ -1149,9 +1153,190 @@ contains
          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
          integer i,ix1,ix2,wx,wy
@@ -1365,13 +1550,28 @@ contains
  !             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
               
            
@@ -1474,14 +1674,20 @@ contains
         fname =trim(trim(rootname)//trim(numcat('_p',colix(j)-2)))
         write (aunit,'(a)') "load (fullfile('"//trim(plot_data_dir)//"','" // trim(fname)//  '.dat''));'
         write (aunit,'(a)') 'pts='//trim(fname)//';'
-        write (aunit,*) 'plot(pts(:,1),pts(:,2),lineM{1},''LineWidth'',lw1);'
+!addon        
+        !write (aunit,*) 'plot(pts(:,1),pts(:,2),lineM{1},''LineWidth'',lw1);'
+        write (aunit,*) trim(plottype)//'(pts(:,1),pts(:,2),lineM{1},''LineWidth'',lw1);'
+!end addon
         write (aunit,*) 'axis([-Inf,Inf,0,1.1]);axis manual;'
         write (aunit,*) 'set(gca,''FontSize'',axes_fontsize); hold on;'
 
         if (plot_meanlikes) then
          write (aunit,'(a)') "load (fullfile('"//trim(plot_data_dir)//"','" // trim(fname)//  '.likes''));'
          write (aunit,'(a)') 'pts='//trim(fname)//';'
-         write (aunit,*) 'plot(pts(:,1),pts(:,2),lineL{1},''LineWidth'',lw1);'
+!addon
+!         write (aunit,*) 'plot(pts(:,1),pts(:,2),lineL{1},''LineWidth'',lw1);'
+         write (aunit,*) trim(plottype)//'(pts(:,1),pts(:,2),lineL{1},''LineWidth'',lw1);'
+!end addon
         end if
 
         if (has_markers(colix(j))) then
@@ -1495,14 +1701,18 @@ contains
 
           write (aunit,'(a)') "pts = load (fullfile('"//trim(plot_data_dir)//"','" // trim(fname)// '.dat''));'
           fmt = trim(numcat('lineM{',ix1+1)) // '}'        
-          write (aunit,'(a)') 'plot(pts(:,1),pts(:,2),'//trim(fmt)//',''LineWidth'',lw2);'
-    
+!addon
+!          write (aunit,'(a)') 'plot(pts(:,1),pts(:,2),'//trim(fmt)//',''LineWidth'',lw2);'
+          write (aunit,'(a)') trim(plottype)//'(pts(:,1),pts(:,2),'//trim(fmt)//',''LineWidth'',lw2);'
+!end addon
          if (plot_meanlikes) then
 
             write (aunit,'(a)') "pts=load (fullfile('"//trim(plot_data_dir)//"','" //trim(fname)//'.likes''));'
            fmt = trim(numcat('lineL{',ix1+1)) // '}'        
-            write (aunit,'(a)') 'plot(pts(:,1),pts(:,2),'//trim(fmt)//',''LineWidth'',lw2);'
-
+!addon
+!            write (aunit,'(a)') 'plot(pts(:,1),pts(:,2),'//trim(fmt)//',''LineWidth'',lw2);'
+           write (aunit,'(a)') trim(plottype)//'(pts(:,1),pts(:,2),'//trim(fmt)//',''LineWidth'',lw2);'
+!end addon
          end if
 
          end if
@@ -1607,7 +1817,9 @@ program GetDist
 
         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    
@@ -1825,6 +2037,15 @@ program GetDist
           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.)
+    plottype = Ini_Read_String('plot_type',.false.)
+
+    if (plottype =='') then
+       plottype='plot'
+    endif
+!end addon
 
     matlab_col = Ini_Read_String('matlab_colscheme')
     smoothing = Ini_Read_Logical('smoothing',.true.)
@@ -2457,7 +2678,14 @@ program GetDist
               if (line_labels) call WriteMatlabLineLabels(51)
               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
@@ -2668,6 +2896,35 @@ program GetDist
           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
+                call ParamNames_ReadIndices(NameMapping,plot_3D(i), tmp_params, 3)
+                ix1 = tmp_params(1)
+                ix2 = tmp_params(2)
+                ix3 = tmp_params(3)
+                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 --git a/source/MCMC.f90 b/source/MCMC.f90
index 9c3285d..a657e5b 100644
--- a/source/MCMC.f90
+++ b/source/MCMC.f90
@@ -215,7 +215,7 @@ contains
   
     Like = GetLogLike(grid(r))
    
-    if (Feedback > 1) write (*,*) r, 'Likelihood: ', Like, 'Current Like:', CurLike
+    if (Feedback > 1) write (*,*) r, 'Likelihood: ', Like, 'Current Like:', CurLike, 'Rank:',MPIRank
 
     if ((Like /= logZero) .and. (CurLike > Like .or. randexp1() > Like - CurLike)) then
     !Accept
@@ -501,7 +501,7 @@ function WL_Weight(L) result (W)
            output_lines = output_lines +1
            call WriteParams(CurParams, real(mult), CurLike)
           end if
-          if (Feedback > 1) write (*,*) instance, 'Slicing, Current Like:', CurLike
+          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) 
@@ -536,7 +536,7 @@ function WL_Weight(L) result (W)
 
     Like = GetLogLike(Trial) 
 
-    if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:', CurLike
+    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 --git a/source/Makefile b/source/Makefile
index 3a816a7..1fd7df6 100644
--- a/source/Makefile
+++ b/source/Makefile
@@ -1,122 +1,83 @@
-#You may need to edit the library paths for MKL for Intel
-#Beware of using optmizations that lose accuracy - may give errors when running
-
-##Uncomment the next line to include dr7 LRG
-EXTDATA =
-#EXTDATA = LRG
-
-#set WMAP empty not to compile with WMAP
-WMAP = /home/aml1005/WMAP7/likelihood_v4
-
-#Only needed for WMAP
-cfitsio = /usr/local/cfitsio/intel10/64/3.040
-
-#GSL only needed for DR7 LRG
-GSLPATH = /home/aml1005/libs/gsl
-
-IFLAG = -I
-INCLUDE=
-
-#Intel MPI (assuming mpif77 set to point to ifort)
-#these settings for ifort 11.1 and higher; may need to add explicit link directory otherwise
-#Can add -xHost if your cluster is uniform, or specify specific processor optimizations -x...
-F90C     = mpif90
-FFLAGS = -O3 -W0 -WB -openmp -fpp -DMPI -vec_report0 -mkl=parallel
-LAPACKL = -lmpi
-
-#GFortran: defaults for v4.5; if pre v4.3 add -D__GFORTRAN__
-#F90C     = gfortran
-#in earlier versions use FFLAGS = -O2 -ffree-form -x f95-cpp-input -D__GFORTRAN__
-#FFLAGS =  -O3 -fopenmp -ffree-form -x f95-cpp-input -ffast-math -march=native -funroll-loops
-#LAPACKL = -Wl,-framework -Wl,accelerate
-#commented above is (I think) for Mac; this is standard linux (sudo apt-get install liblapack-dev)
-#LAPACKL = -lblas -llapack
-
-#HPCF settings. Use Intel 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/10.2.2.025/lib/em64t -lmkl_lapack -lmkl -lguide -lpthread
-#GSLPATH = /usr/local/Cluster-Apps/gsl/1.9
-#cfitsio = /usr/local/Cluster-Users/cpac/cmb/2.1.0/cfitsio
-#INCLUDE=
-
-#COSMOS: use "module load cosmolib latest"
-#use "runCosmomc" (globally installed) to run, defining required memory usage
-ifeq ($(COSMOHOST),cosmos)
-F90C = ifort
-FFLAGS = -openmp -fast -w -fpp2 -DMPI
-LAPACKL = -mkl=sequential -lmkl_lapack -lmpi
-cfitsio = $(CFITSIO)
-WMAP = $(COSMOLIB)/WMAP7
-GSLPATH = $(GSL_ROOT)
+# >>> DESIGNED FOR GMAKE <<<
+
+# Unified Systems makefile for COSMOMC
+# Add FLAGS -DMPI for using MPI
+
+ext=$(shell uname | cut -c1-3)
+
+CC=cc
+
+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
+GSLDIR =
+endif
+
+ifeq ($(ext),Lin)
+F90C=gfortran
+FFLAGS= -O -cpp -fopenmp
+WMAPFLAGS= -O
+LAPACKL = -llapack -lblas
+INCLUDE = -I../camb 
+CFITSIODIR = /usr
+GSLDIR = /usr
+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 = 
+GSLDIR =
+endif
+
+ifeq ($(ext),Sun)
+F90C=f90
+FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
+WMAPFLAGS= $(FFLAGS)
+LAPACKL = -lsunperf -lfsu
+INCLUDE = -I../camb -M../camb
+CFITSIODIR = 
+GSLDIR =
+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 = 
+GSLDIR =
 endif
 
-#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
-
-#G95; make sure LAPACK and MPI libs also compiled with g95
-#F90C     = mpif90
-#FFLAGS = -O2 -cpp -DMPI
-#LAPACKL = /LAPACK/lapack_LINUX.a /LAPACK/blas_LINUX.a
-
-#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
-
-#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 -Mpreprocess
-#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 =
+EXTDATA = LRG
+EXTINCLUDE = -I$(GSLDIR)/include
+EXTOBJS = bsplinepk.o
+
+
+WMAPDIR = ../WMAP
+WMAPINCLUDE = -I$(CFITSIODIR)/include
+WMAPOBJS = read_archive_map.o \
+	read_fits.o \
+	healpix_types.o \
+	br_mod_dist.o \
+	WMAP_7yr_options.o \
+	WMAP_7yr_util.o \
+	WMAP_7yr_gibbs.o \
+	WMAP_7yr_tt_pixlike.o \
+	WMAP_7yr_tt_beam_ptsrc_chisq.o \
+	WMAP_7yr_teeebb_pixlike.o \
+	WMAP_7yr_tetbeebbeb_pixlike.o \
+	WMAP_7yr_likelihood.o
+
+
 
 PROPOSE = propose.o
 CLSFILE = CMB_Cls_simple.o
@@ -124,51 +85,42 @@ CLSFILE = CMB_Cls_simple.o
 #Can use params_H if you prefer more generic parameters
 PARAMETERIZATION = params_CMB.o
 
-F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) $(IFLAG)../camb $(INCLUDE)
-LINKFLAGS = -L../camb -lcamb $(LAPACKL)
+LINKFLAGS =  -L../camb -lcamb  $(LAPACKL)
 
-DISTFILES = ParamNames.o Matrix_utils.o settings.o IO.o GetDist.o
+F90FLAGS = -DMATRIX_SINGLE $(FFLAGS) $(INCLUDE)
 
-OBJFILES=  ParamNames.o Matrix_utils.o settings.o IO.o cmbtypes.o Planck_like.o  \
-	cmbdata.o WeakLen.o bbn.o lrggettheory.o mpk.o bao.o supernovae.o HST.o SDSSLy-a-v3.o \
+DISTFILES = utils.o ParamNames.o Matrix_utils.o settings.o IO.o GetDist.o
+
+OBJFILES = utils.o ParamNames.o Matrix_utils.o settings.o IO.o cmbtypes.o Planck_like.o \
+	cmbdata.o WeakLen.o bbn.o lrggettheory.o mpk.o bao.o supernovae.o HST.o SDSSLy-a-v3.o\
 	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
 	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
 
-
-ifeq ($(EXTDATA),LRG)
+ifeq ($(EXTDATA),LRG)	
 F90FLAGS += -DDR71RG
-OBJFILES += bsplinepk.o
-GSLINC = -I$(GSLPATH)/include
-LINKFLAGS += -L$(GSLPATH)/lib -lgsl -lgslcblas
-endif
-
-ifneq ($(WMAP),)
-F90FLAGS += $(IFLAG)$(cfitsio)/include $(IFLAG)$(WMAP)
-LINKFLAGS +=  -L$(cfitsio)/lib -L$(WMAP) -lcfitsio
-
-OBJFILES += $(WMAP)/read_archive_map.o \
-	$(WMAP)/read_fits.o \
-	$(WMAP)/healpix_types.o \
-	$(WMAP)/WMAP_7yr_options.o \
-	$(WMAP)/WMAP_7yr_util.o \
-	$(WMAP)/WMAP_7yr_tt_pixlike.o \
-	$(WMAP)/WMAP_7yr_teeebb_pixlike.o \
-	$(WMAP)/WMAP_7yr_likelihood.o \
-	$(WMAP)/WMAP_7yr_gibbs.o \
-	$(WMAP)/WMAP_7yr_tt_beam_ptsrc_chisq.o \
-	$(WMAP)/br_mod_dist.o
+EXTINCLUDE = -I$(GSLDIR)/include
+LINKFLAGS += -lgsl -lgslcblas
+OBJFILES += $(EXTOBJS)
+endif	
+
+ifneq ($(WMAPDIR),)
+F90FLAGS += $(WMAPINCLUDE)
+LINKFLAGS += -lcfitsio
+OBJFILES += $(WMAPOBJS)
 else
 F90FLAGS += -DNOWMAP
 endif
 
-default: cosmomc
 
-all : cosmomc getdist
 
-settings.o: ../camb/libcamb.a
+default: cosmomc.$(ext)
+
+all : cosmomc.$(ext) getdist.$(ext)
+
+settings.o: utils.o
 cmbtypes.o: settings.o
 Planck_like.o: cmbtypes.o
-cmbdata.o: Planck_like.o  $(WMAPOBJS)
+cmbdata.o: Planck_like.o $(WMAPOBJS)
 WeakLen.o: cmbtypes.o
 bbn.o: settings.o
 mpk.o: cmbtypes.o lrggettheory.o
@@ -187,14 +139,20 @@ postprocess.o: calclike.o
 MCMC.o: calclike.o
 driver.o: MCMC.o
 
-export FFLAGS
-export F90C
-
 .f.o:
 	f77 $(F90FLAGS) -c $<
 
 %.o: %.c
-	$(CC) $(GSLINC) -c $*.c
+	$(CC) $(EXTINCLUDE) -c $*.c
+
+%.o: $(WMAPDIR)/%.f90
+	$(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $<
+
+%.o: $(WMAPDIR)/%.F90
+	$(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $<
+
+utils.o: ../camb/utils.F90
+	$(F90C) $(F90FLAGS) -c $<
 
 %.o: %.f90
 	$(F90C) $(F90FLAGS) -c $*.f90
@@ -203,19 +161,12 @@ export F90C
 	$(F90C) $(F90FLAGS) -c $*.F90
 
 
-cosmomc: camb $(OBJFILES)
-	$(F90C) -o ../cosmomc $(OBJFILES) $(LINKFLAGS) $(F90FLAGS)
+cosmomc.$(ext): $(OBJFILES) ../camb/libcamb.a 	
+	$(F90C) -o ../$@ $(OBJFILES) $(LINKFLAGS) $(F90FLAGS)
 
 
-clean: cleancosmomc
-	rm -f ../camb/*.o ../camb/*.obj ../camb/*.mod
-
-cleancosmomc:
+clean:
 	rm -f *.o *.mod *.d *.pc *.obj ../core
 
-
-getdist: camb $(DISTFILES)
-	$(F90C) -o ../getdist $(DISTFILES) $(LINKFLAGS) $(F90FLAGS)
-
-camb:
-	cd ../camb && $(MAKE) --file=Makefile_main libcamb.a
+getdist.$(ext): $(DISTFILES)
+	$(F90C) -o ../$@ $(DISTFILES) $(LINKFLAGS) $(F90FLAGS) 
diff --git a/source/driver.F90 b/source/driver.F90
index c323451..4460128 100644
--- a/source/driver.F90
+++ b/source/driver.F90
@@ -58,12 +58,19 @@ program SolveCosmology
         end if
 
 #ifdef MPI 
-
-        if (instance /= 0) call DoAbort('With MPI should not have second parameter')
+!addon
+!        if (instance /= 0) call DoAbort('With MPI should not have second parameter')
         call mpi_comm_rank(mpi_comm_world,MPIrank,ierror)
-
+        instance_shift=instance
+        write(*,*)'MPI file names shifted by: ',instance_shift
+!end addon
         instance = MPIrank +1 !start at 1 for chains                                             
-        write (numstr,*) instance
+!addon
+!        write (numstr,*) instance
+        instance_shift = MPIrank + instance_shift
+        write (numstr,*) instance_shift
+!end addon
+
         rand_inst = instance
         if (ierror/=MPI_SUCCESS) call DoAbort('MPI fail')
 
@@ -104,6 +111,11 @@ program SolveCosmology
 
         HighAccuracyDefault = Ini_Read_Logical('high_accuracy_default',.false.)
         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
         
         if (Ini_HasKey('highL_unlensed_cl_template')) &
           highL_unlensed_cl_template=  ReadIniFilename(DefIni,'highL_unlensed_cl_template')
diff --git a/source/settings.f90 b/source/settings.f90
index 90d84d0..f163ff6 100644
--- a/source/settings.f90
+++ b/source/settings.f90
@@ -6,6 +6,11 @@ module settings
   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)
 
@@ -66,6 +71,9 @@ module settings
 
   integer :: num_threads = 0
   integer :: instance = 0
+!addon
+  integer :: instance_shift = 0
+!end addon
   integer :: MPIchains = 1, MPIrank = 0
 
   logical :: Use_LSS = .true.