diff -r -c -b -B -N cosmomc/camb/inidriver.F90 cosmomc_addon/camb/inidriver.F90
*** cosmomc/camb/inidriver.F90	2005-03-31 00:37:14.000000000 +0100
--- cosmomc_addon/camb/inidriver.F90	2006-02-20 20:32:38.000000000 +0000
***************
*** 171,193 ****
               stop 'Unknown recombination'
             end if
  
             P%InitPower%nn = Ini_Read_Int('initial_power_num')
             if (P%InitPower%nn>nnmax) stop 'Too many initial power spectra - increase nnmax in InitialPower'
             P%InitPower%rat(:) = 1
             do i=1, P%InitPower%nn
                write (numstr,*) i 
                numstr=adjustl(numstr)
                P%InitPower%an(i) = &
                     Ini_Read_Double('scalar_spectral_index('//trim(numstr)//')')
  
                P%InitPower%n_run(i) = &
                     Ini_Read_Double('scalar_nrun('//trim(numstr)//')',0._dl)
      
                if (P%WantTensors) then
                   P%InitPower%ant(i) = Ini_Read_Double('tensor_spectral_index('//trim(numstr)//')')
                   P%InitPower%rat(i) = Ini_Read_Double('initial_ratio('//trim(numstr)//')')
                end if              
! 
                P%InitPower%ScalarPowerAmp(i) = Ini_Read_Double('scalar_amp('//trim(numstr)//')',1.d0) 
                !Always need this as may want to set tensor amplitude even if scalars not computed
             end do
--- 171,214 ----
               stop 'Unknown recombination'
             end if
  
+ !addon
+             P%InitPower%WantExternalScalar = Ini_Read_Logical('use_external_scalar',.false.)
+             P%InitPower%WantExternalTensor = Ini_Read_Logical('use_external_tensor',.false.)
+ !end addon
+ 
             P%InitPower%nn = Ini_Read_Int('initial_power_num')
             if (P%InitPower%nn>nnmax) stop 'Too many initial power spectra - increase nnmax in InitialPower'
             P%InitPower%rat(:) = 1
             do i=1, P%InitPower%nn
                write (numstr,*) i 
                numstr=adjustl(numstr)
+ !addon
+               if (P%WantScalars .or. P%WantTransfer) then
+                  if (P%InitPower%WantExternalScalar) then
+                     P%InitPower%ScalarPowerFileName(i) & 
+                          = Ini_Read_String('scalar_input_file('//trim(numstr)//')')
+                  else
+ 
                      P%InitPower%an(i) = &
                           Ini_Read_Double('scalar_spectral_index('//trim(numstr)//')')
  
                      P%InitPower%n_run(i) = &
                           Ini_Read_Double('scalar_nrun('//trim(numstr)//')',0._dl)
  
+                  endif
+               endif
+ 
                if (P%WantTensors) then
+ 
+                  if (P%InitPower%WantExternalTensor) then
+                     P%InitPower%TensorPowerFileName(i) & 
+                          = Ini_Read_String('tensor_input_file('//trim(numstr)//')')
+                  else
                      P%InitPower%ant(i) = Ini_Read_Double('tensor_spectral_index('//trim(numstr)//')')
                      P%InitPower%rat(i) = Ini_Read_Double('initial_ratio('//trim(numstr)//')')
+                  endif
                end if              
! !end addon
                P%InitPower%ScalarPowerAmp(i) = Ini_Read_Double('scalar_amp('//trim(numstr)//')',1.d0) 
                !Always need this as may want to set tensor amplitude even if scalars not computed
             end do
***************
*** 247,252 ****
--- 268,282 ----
         AccuracyBoost  = Ini_Read_Double('accuracy_boost',1.d0)
         lAccuracyBoost = Ini_Read_Double('l_accuracy_boost',1.d0)
         lSampleBoost   = Ini_Read_Double('l_sample_boost',1.d0)
+ !addon     
+        DoClsInterpol  = Ini_Read_Logical('do_cls_interpol',.false.)
+        ComputeEveryl  = Ini_Read_Logical('compute_every_l',.false.)
+   
+        if (FeedbackLevel > 0) then
+           write (*,*)'do_cls_interpol is',DoClsInterpol
+           write (*,*)'compute_every_l is',ComputeEveryl          
+        endif
+ !end addon
  
         if (outroot /= '') then
           call Ini_SaveReadValues(trim(outroot) //'params.ini',1)
***************
*** 290,296 ****
  #endif
          end if
  
!              
    
          end program driver
  
--- 320,329 ----
  #endif
          end if
  
! !addon
!         if (P%InitPower%WantExternalScalar.or.P%InitPower%WantExternalTensor) &
!              call FreeExternalData(P%InitPower)
! !end addon
    
          end program driver
  
diff -r -c -b -B -N cosmomc/camb/Makefile cosmomc_addon/camb/Makefile
*** cosmomc/camb/Makefile	2006-03-14 11:07:24.000000000 +0000
--- cosmomc_addon/camb/Makefile	2006-03-26 19:13:14.000000000 +0100
***************
*** 1,57 ****
! #CAMB Makefile
! #Edit for your compiler
! #Note there are many ifc versions, some of which behave oddly
! 
! #Intel 8 , -openmp toggles mutli-processor:
! F90C     = ifort
! FFLAGS = -openmp -ip -O2 -Vaxlib -W0 -WB -fpp2
! 
! # Intel 9 on IA-64 (eg. COSMOS)
! # (do "module load icomp90" before compiling)
! #F90C = ifort
! #FFLAGS = -openmp -fpp2 -w -O3 -ip -mP2OPT_hlo_prefetch=F
! 
! #Intel ifc, add -openmp for multi-processor (some have bugs):
! #F90C     = ifc
! #FFLAGS = -O2 -Vaxlib -ip -W0 -WB -quiet -fpp2
! #some systems can can also add e.g. -tpp7 -xW
! 
! #SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems.
! #F90C     = f90
! #FFLAGS  = -Ofast -mp
! 
! #Digital/Compaq fortran, -omp toggles multi-processor
! #F90C    = f90
! #FFLAGS  = -omp -O4 -arch host -math_library fast -tune host -fpe1
! 
! #Absoft ProFortran, single processor:
! #F90C     = f95
! #FFLAGS = -O2 -cpu:athlon -s -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_"
! 
! #NAGF95, single processor:
! #F90C     = f95
! #FFLAGS = -DNAGF95 -O3
! 
! #G95 compiler
! #F90C   = g95
! #FFLAGS = -O2
! 
! #Sun, single processor:
! #F90C     = f90
! #FFLAGS = -O5
! 
! #Sun V880
! #F90C = mpf90
! #FFLAGS =  -O4 -openmp -ftrap=%none -dalign -DMPI
! 
! #Sun parallel enterprise:
! #F90C     = f95
! #FFLAGS =  -O2 -xarch=native64 -openmp -ftrap=%none
! #try removing -openmp if get bus errors. -03, -04 etc are dodgy.
! 
! #IBM XL Fortran, multi-processor (run gmake)
! #F90C     = xlf90_r
! #FFLAGS  = -qsmp=omp -qsuffix=f=f90:cpp=F90 -O3 -qstrict -qarch=pwr3 -qtune=pwr3
  
  
  #Files containing evolution equations initial power spectrum module
--- 1,33 ----
! # >>> DESIGNED FOR GMAKE <<<
! 
! #CAMB System unified Makefile
! 
! ext=$(shell uname | cut -c1-3)
! 
! ifeq ($(ext),IRI)
! F90C=f90
! FFLAGS = -O3 -mp -64
! endif
! 
! ifeq ($(ext),Lin)
! F90C=ifort
! FFLAGS= -O -openmp  -fpp
! endif 
! 
! ifeq ($(ext),OSF)
! F90C=f90
! FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
! endif
! 
! ifeq ($(ext),Sun)
! F90C=f90
! FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
! endif
! 
! ifeq ($(ext),AIX)
! F90C=xlf90_r
! FFLAGS= -O -q64 -qsmp=omp -qthreaded  -qmaxmem=-1 -qfree=f90 -qsuffix=f=f90:cpp=F90
! endif
  
  
  #Files containing evolution equations initial power spectrum module
***************
*** 83,99 ****
  CAMBOBJ       = subroutines.o inifile.o $(POWERSPECTRUM).o recfast.o modules.o \
  	bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o cmbmain.o camb.o
  
! default: camb
  
! all: camb $(CAMBLIB)
  
! camb: $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) -o $@
  
  $(CAMBLIB): $(CAMBOBJ)
  	ar -r $@ $?
  
! camb_fits: writefits.f90 $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) -DWRITE_FITS -o $@
  
  %.o: %.f90
--- 59,75 ----
  CAMBOBJ       = subroutines.o inifile.o $(POWERSPECTRUM).o recfast.o modules.o \
  	bessels.o $(EQUATIONS).o $(NONLINEAR).o lensing.o cmbmain.o camb.o
  
! default: camb.$(ext)
  
! all: camb.$(ext) $(CAMBLIB)
  
! camb.$(ext): $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) $(CAMBOBJ) $(DRIVER) -o $@
  
  $(CAMBLIB): $(CAMBOBJ)
  	ar -r $@ $?
  
! camb_fits.$(ext): writefits.f90 $(CAMBOBJ) $(DRIVER)
  	$(F90C) $(F90FLAGS) -I$(HEALPIXDIR)/include $(CAMBOBJ) writefits.f90 $(DRIVER) $(HEALPIXLD) -DWRITE_FITS -o $@
  
  %.o: %.f90
diff -r -c -b -B -N cosmomc/camb/modules.f90 cosmomc_addon/camb/modules.f90
*** cosmomc/camb/modules.f90	2006-03-16 11:07:20.000000000 +0000
--- cosmomc_addon/camb/modules.f90	2006-03-17 22:37:18.000000000 +0000
***************
*** 186,192 ****
        integer, parameter:: l0max=4000
  
  !     lmax is max possible number of l's evaluated
!       integer, parameter :: lmax_arr = 100+l0max/7
   
          contains
        
--- 186,197 ----
        integer, parameter:: l0max=4000
  
  !     lmax is max possible number of l's evaluated
! !addon
! !      integer, parameter :: lmax_arr = 100+l0max/7
!       integer, parameter :: lmax_arr = 100 + l0max
!       logical :: DoClsInterpol=.true.
!       logical :: ComputeEveryl=.true.
! !end addon
   
          contains
        
***************
*** 539,544 ****
--- 544,567 ----
          integer lind, lvar, step,top,bot,ls(lmax_arr)
          real(dl) AScale
      
+ !addon    
+         if (ComputeEveryl) then
+            write(*,*) 'Warning: computing every l'           
+            if (max_l.gt.lmax_arr) then
+               write(*,*)'lmax_arr maxl',lmax_arr,max_l
+               stop
+            endif
+            lind=0
+            do lvar=2, max_l
+               lind=lind+1
+               ls(lind)=lvar 
+            end do
+ 
+            lSet%l0=lind
+            lSet%l = ls(1:lind)
+ 
+         else
+ !end addon
          Ascale=scale/lSampleBoost       
        
          lind=0
***************
*** 641,646 ****
--- 664,673 ----
          lSet%l0=lind
          lSet%l(1:lind) = ls(1:lind)
          
+ !addon
+         endif
+ !end addon
+         
        end subroutine initlval
  
        subroutine InterpolateClArr(lSet,iCl, all_Cl, max_ind)
***************
*** 660,665 ****
--- 687,694 ----
        xl = real(lSet%l(1:lSet%l0),dl)
        call spline(xl,iCl(1),max_ind,cllo,clhi,ddCl(1))
       
+ !addon
+       if (DoClsInterpol) then
              llo=1
              do il=lmin,lSet%l(max_ind)
                 xi=il
***************
*** 676,681 ****
--- 705,727 ----
                
              end do
  
+       else
+ 
+          llo=1
+          do il=2,lSet%l(max_ind)
+             xi=il
+             if ((xi > lSet%l(llo+1)).and.(llo < max_ind)) then
+                llo=llo+1
+             end if
+ 
+             all_Cl(il) = iCl(llo+1)
+ 
+          end do
+          all_Cl(2)=iCl(1)
+       endif
+ !end addon
+ 
+ 
        end subroutine InterpolateClArr
  
   
diff -r -c -b -B -N cosmomc/camb/params.ini cosmomc_addon/camb/params.ini
*** cosmomc/camb/params.ini	2005-03-30 16:05:28.000000000 +0100
--- cosmomc_addon/camb/params.ini	2006-02-20 17:31:02.000000000 +0000
***************
*** 64,69 ****
--- 64,76 ----
  
  #Initial power spectrum, amplitude, spectral index and running
  initial_power_num         = 1
+ 
+ #read the power spectra form ascii files in the form (Mpc^-1): k P(k)
+ use_external_scalar       = F
+ use_external_tensor       = F
+ scalar_input_file(1)      = scalar_power.dat
+ tensor_input_file(1)      = tensor_power.dat
+ 
  scalar_amp(1)             = 2.3e-9
  scalar_spectral_index(1)  = 1
  scalar_nrun(1)            = 0
***************
*** 100,105 ****
--- 107,118 ----
  #Matter power spectrum output against k/h in units of h^3 Mpc^{-3}
  transfer_matterpower(1) = matterpower.dat
  
+ #Use spline interpolation for the Cls
+ do_cls_interpol = T
+ 
+ #Remove undersampling of the Cls: if true lSampleBoost has no effect
+ compute_every_l = F
+ 
  
  #Output files not produced if blank. make camb_fits to use use the FITS setting.
  scalar_output_file = scalCls.dat
diff -r -c -b -B -N cosmomc/camb/power_tilt.f90 cosmomc_addon/camb/power_tilt.f90
*** cosmomc/camb/power_tilt.f90	2005-03-31 00:37:14.000000000 +0100
--- cosmomc_addon/camb/power_tilt.f90	2006-02-20 20:35:53.000000000 +0000
***************
*** 51,60 ****
         real(dl) k_0_scalar, k_0_tensor
         real(dl) ScalarPowerAmp(nnmax)
   
        end Type InitialPowerParams
  
!       real(dl) curv  !Curvature contant, set in InitializePowers     
  
        
        Type(InitialPowerParams) :: P
    
--- 51,80 ----
         real(dl) k_0_scalar, k_0_tensor
         real(dl) ScalarPowerAmp(nnmax)
   
+ !addon 
+        logical :: WantExternalScalar
+        logical :: WantExternalTensor
+        character(80) :: ScalarPowerFileName(nnmax),TensorPowerFileName(nnmax)
+ !end addon      
+ 
        end Type InitialPowerParams
  
! !addon
!       type ExternalData
! !only pointers may be defered shape in derived type
!          real(dl) :: ScalarMomentumMin, ScalarMomentumMax
!          real(dl), dimension(:), pointer :: PtrScalarPower, PtrScalarLnMomentum
!          real(dl), dimension(:), pointer :: PtrScalarPower2
!          real(dl) :: TensorMomentumMin, TensorMomentumMax
!          real(dl), dimension(:), pointer :: PtrTensorPower, PtrTensorLnMomentum
!          real(dl), dimension(:), pointer :: PtrTensorPower2
!       end type ExternalData
  
+       type(ExternalData), dimension(:), pointer  :: PtrExternalData => null()
+ !end addon      
+ 
+ 
+       real(dl) curv  !Curvature contant, set in InitializePowers     
        
        Type(InitialPowerParams) :: P
    
***************
*** 63,68 ****
--- 83,91 ----
        public  InitializePowers, ScalarPower, TensorPower
        public nnmax,InitialPowerParams,Power_Descript, Power_Name, SetDefPowerParams
  
+ !addon
+       public FreeExternalData
+ !end addon
      contains
         
    
***************
*** 77,82 ****
--- 100,110 ----
           AP%k_0_scalar = 0.05
           AP%k_0_tensor = 0.05
           AP%ScalarPowerAmp = 1
+ !addon
+ !         nullify(PtrExternalData)
+          AP%WantExternalScalar = .false.
+          AP%WantExternalTensor = .false.
+ !end addon        
  
         end subroutine SetDefPowerParams
  
***************
*** 87,92 ****
--- 115,128 ----
  
          real(dl) acurv
  
+ !addon
+         integer :: i,ioflag,j
+         integer :: nPowerUnit=111
+         integer :: nData
+         real(dl), save :: sbuffer1,sbuffer2
+         real(dl), parameter :: spl_large=1.e40_dl     
+ !end addon
+ 
          P = AParamSet
          if (P%nn > nnmax) then
             write (*,*) 'To use ',P%nn,'power spectra you need to increase'
***************
*** 97,105 ****
--- 133,295 ----
  
  !Write implementation specific code here...
  
+ !addon
+ 
+ !create and point to P%nn ExternalData objects (pointers and reals)
+         if ((P%WantExternalScalar).or.(P%WantExternalTensor)) then
+            if (.not.associated(PtrExternalData)) then
+               allocate(PtrExternalData(P%nn)) 
+            else
+               write(*,*)'InitializePowers: PtrExternalData already Initialized!'
+               return
+            endif
+         endif
+          
+         if (P%WantExternalScalar) then
+                         
+            do i=1,P%nn
+                
+               write(*,*)'using primordial scalar power spectrum: ',P%ScalarPowerFileName(i)
+                
+               open(unit=nPowerUnit,file=P%ScalarPowerFileName(i),action='read' &
+                    ,form='formatted',status='old')             
+ !count the number of records
+               nData=0
+               do                                    
+                  read(nPowerUnit,*,iostat=ioflag)                  
+                  if (ioflag.lt.0) exit
+                  nData=nData + 1                                    
+               enddo
+ !create and point to
+               allocate(PtrExternalData(i)%PtrScalarPower(nData) & 
+                    ,PtrExternalData(i)%PtrScalarLnMomentum(nData))
+               
+               rewind(nPowerUnit)
+               nData=1
+               do                                   
+                  read(nPowerUnit,*,iostat=ioflag) sbuffer1,sbuffer2
+                  if (ioflag.lt.0) exit
+                  PtrExternalData(i)%PtrScalarLnMomentum(nData) = log(sbuffer1)
+                  PtrExternalData(i)%PtrScalarPower(nData) = sbuffer2      
+                  nData=nData+1
+               enddo
+               nData=nData - 1 
+               close(nPowerUnit,iostat=ioflag)
+             
+ !keep trace of the accessible range              
+               PtrExternalData(i)%ScalarMomentumMin = exp(PtrExternalData(i)%PtrScalarLnMomentum(1))
+               PtrExternalData(i)%ScalarMomentumMax = exp(PtrExternalData(i)%PtrScalarLnMomentum(nData))
+               
+ 
+ !2nd order derivatives allocation
+               allocate(PtrExternalData(i)%PtrScalarPower2(nData))
+ !splining in log k
+               call spline(PtrExternalData(i)%PtrScalarLnMomentum(:) &
+                    ,PtrExternalData(i)%PtrScalarPower(:),nData, spl_large &
+                    , spl_large,PtrExternalData(i)%PtrScalarPower2(:))
+               
+            enddo
+ 
+         endif
+ 
+ 
+ !the same for tensors
+ 
+         if (P%WantExternalTensor) then
+             
+            do i=1,P%nn
+               
+               write(*,*)'using primordial tensor power spectrum: ',P%TensorPowerFileName(i)
+ 
+               open(unit=nPowerUnit,file=P%TensorPowerFileName(i) &
+                    ,action='read',form='formatted',status='old')     
+ !count the number of records
+               nData=0
+               do                                    
+                  read(nPowerUnit,*,iostat=ioflag)                  
+                  if (ioflag.lt.0) exit
+                  nData=nData + 1                                    
+               enddo
+ !store the data
+               allocate(PtrExternalData(i)%PtrTensorPower(nData) & 
+                    ,PtrExternalData(i)%PtrTensorLnMomentum(nData))
+                
+               rewind(nPowerUnit)
+               nData=1
+               do                                   
+                  read(nPowerUnit,*,iostat=ioflag) sbuffer1,sbuffer2
+                  if (ioflag.lt.0) exit
+                  PtrExternalData(i)%PtrTensorLnMomentum(nData) = log(sbuffer1)
+                  PtrExternalData(i)%PtrTensorPower(nData) = sbuffer2      
+                  nData=nData+1
+               enddo
+               nData=nData - 1 
+               close(nPowerUnit)
+ 
+ !keep trace of the accessible range              
+               PtrExternalData(i)%TensorMomentumMin = exp(PtrExternalData(i)%PtrTensorLnMomentum(1))
+               PtrExternalData(i)%TensorMomentumMax = exp(PtrExternalData(i)%PtrTensorLnMomentum(nData))
+               
+ 
+ !2nd order derivatives
+               allocate(PtrExternalData(i)%PtrTensorPower2(nData))
+ !splining in log k
+               call spline(PtrExternalData(i)%PtrTensorLnMomentum(:) &
+                    ,PtrExternalData(i)%PtrTensorPower(:),nData, spl_large &
+                    , spl_large,PtrExternalData(i)%PtrTensorPower2(:))
+               
+            enddo
+ 
+             
+         endif
+          
         end subroutine InitializePowers
         
  
+ 
+        subroutine FreeExternalData(AP)
+          type (InitialPowerParams), intent(in) :: AP
+          integer :: i
+ 
+          if (AP%nn.ne.size(PtrExternalData,1)) then
+             write(*,*)'InitialPower.FreeExternalData : wrong initial power params'
+             stop
+          endif
+ 
+          if (AP%WantExternalScalar) then
+             do i=1, AP%nn               
+ !free heap arrays and disassociate ExternalData pointers            
+                deallocate(PtrExternalData(i)%PtrScalarPower)
+                deallocate(PtrExternalData(i)%PtrScalarPower2)
+                deallocate(PtrExternalData(i)%PtrScalarLnMomentum)
+ !should still point to ExternalData objects
+                if (.not.associated(PtrExternalData)) stop 'InitialPower.FreeExternalData : ooops'
+             enddo
+          endif
+          if (AP%WantExternalTensor) then
+             do i=1, AP%nn               
+                deallocate(PtrExternalData(i)%PtrTensorPower)
+                deallocate(PtrExternalData(i)%PtrTensorPower2)
+                deallocate(PtrExternalData(i)%PtrTensorLnMomentum)
+                if (.not.associated(PtrExternalData)) stop 'InitialPower.FreeExternalData : ooops'
+             enddo
+          endif
+ !Free remaining objects and disassociate pointer to ExternalData
+          if (associated(PtrExternalData)) then
+             deallocate(PtrExternalData)
+             PtrExternalData => null()
+          else
+             if (AP%WantExternalScalar.or.AP%WantExternalTensor) then
+ !should never happen               
+                stop 'InitialPower.FreeExternalData : ooops'
+             endif
+          endif
+ 
+        end subroutine FreeExternalData
+ 
+ !end addon
+ 
+ 
        function ScalarPower(k,in)
  
         !"in" gives the index of the power to return for this k
***************
*** 120,128 ****
--- 310,340 ----
          real(dl) ScalarPower,k, lnrat
          integer in
  
+ !addon
+         integer :: nData       
+ 
+         if (P%WantExternalScalar) then
+            nData=size(PtrExternalData(in)%PtrScalarLnMomentum,1)
+ 
+            if ((k.ge.PtrExternalData(in)%ScalarMomentumMin) &
+                 .and.(k.le.PtrExternalData(in)%ScalarMomentumMax)) then
+               call splinterpol(PtrExternalData(in)%PtrScalarLnMomentum &
+                    ,PtrExternalData(in)%PtrScalarPower &
+                    ,PtrExternalData(in)%PtrScalarPower2,nData,log(k),ScalarPower)
+ !confusing              ScalarPower = P%ScalarPowerAmp(in)*ScalarPower
+            else
+               write(*,*)'InitPower.ScalarPower() : out of range: k= ',k
+               ScalarPower=0_dl              
+            endif
+ 
+         else
+ 
            lnrat = log(k/P%k_0_scalar)
            ScalarPower=P%ScalarPowerAmp(in)*exp((P%an(in)-1)*lnrat + P%n_run(in)/2*lnrat**2)   
   
+         endif
+ !end addon
+ 
        end function ScalarPower
  
        
***************
*** 143,149 ****
--- 355,384 ----
     
          integer in
  
+ !addon
+         integer :: nData
+         
+ 
+         if (P%WantExternalTensor) then
+            nData=size(PtrExternalData(in)%PtrTensorLnMomentum,1)
+            if ((k.ge.PtrExternalData(in)%TensorMomentumMin) &
+                 .and.(k.le.PtrExternalData(in)%TensorMomentumMax)) then
+               call splinterpol(PtrExternalData(in)%PtrTensorLnMomentum &
+                    ,PtrExternalData(in)%PtrTensorPower &
+                    ,PtrExternalData(in)%PtrTensorPower2,nData,log(k),TensorPower)
+ 
+ !confusing              TensorPower = P%rat(in)*P%ScalarPowerAmp(in)*TensorPower
+            else
+               write(*,*)'InitPower.TensorPower() : out of range: k= ',k
+               TensorPower=0_dl               
+            endif
+         else
+ 
          TensorPower=P%rat(in)*P%ScalarPowerAmp(in)*exp(P%ant(in)*log(k/P%k_0_tensor))
+ 
+         endif
+ !end addon
+         
          if (curv < 0) TensorPower=TensorPower*tanh(PiByTwo*sqrt(-k**2/curv-3)) 
  
        end function TensorPower
diff -r -c -b -B -N cosmomc/camb/subroutines.f90 cosmomc_addon/camb/subroutines.f90
*** cosmomc/camb/subroutines.f90	2005-03-31 00:37:14.000000000 +0100
--- cosmomc_addon/camb/subroutines.f90	2006-02-20 17:04:32.000000000 +0000
***************
*** 962,964 ****
--- 962,997 ----
  !  end abort action
  !
        end subroutine dverk
+ 
+ 
+ !addon
+ !the numrec splint renamed
+ 
+ SUBROUTINE splinterpol(xa,ya,y2a,n,x,y)
+   use Precision
+   INTEGER n
+   REAL(dl) x,y,xa(n),y2a(n),ya(n)
+   INTEGER k,khi,klo
+   REAL(dl) a,b,h
+   klo=1
+   khi=n
+ 1 if (khi-klo.gt.1) then
+      k=(khi+klo)/2
+      if(xa(k).gt.x)then
+         khi=k
+      else
+         klo=k
+      endif
+      goto 1
+   endif
+   h=xa(khi)-xa(klo)
+   if (h.eq.0.) stop 'bad xa input in splint'
+   a=(xa(khi)-x)/h
+   b=(x-xa(klo))/h
+   y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))&
+        *(h**2)/6.
+   return
+ END SUBROUTINE splinterpol
+ 
+ 
+ !end addon
diff -r -c -b -B -N cosmomc/distparams.ini cosmomc_addon/distparams.ini
*** cosmomc/distparams.ini	2006-03-17 23:39:38.000000000 +0000
--- cosmomc_addon/distparams.ini	2006-03-26 18:33:33.000000000 +0100
***************
*** 32,37 ****
--- 32,38 ----
  
  plot_meanlikes = T
  shade_meanlikes = T
+ td_meanlikes = F
  
  # if non-zero, output _thin file, thinned by thin_factor
  thin_factor = 0
diff -r -c -b -B -N cosmomc/params.ini cosmomc_addon/params.ini
*** cosmomc/params.ini	2006-03-17 20:54:54.000000000 +0000
--- cosmomc_addon/params.ini	2006-03-26 18:44:37.000000000 +0100
***************
*** 116,122 ****
  #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
--- 116,122 ----
  #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
***************
*** 132,137 ****
--- 132,140 ----
  #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 -N cosmomc/source/CMB_Cls_simple.f90 cosmomc_addon/source/CMB_Cls_simple.f90
*** cosmomc/source/CMB_Cls_simple.f90	2006-03-17 16:18:11.000000000 +0000
--- cosmomc_addon/source/CMB_Cls_simple.f90	2006-03-18 14:16:08.000000000 +0000
***************
*** 266,271 ****
--- 266,289 ----
     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)
***************
*** 294,302 ****
              P%Transfer%high_precision=.true.
              P%Transfer%k_per_logint=3*nint(AccuracyLevel)
            end if
!           AccuracyBoost = AccuracyLevel
!           lAccuracyBoost = AccuracyLevel
!           lSampleBoost = AccuracyLevel
            P%AccurateReionization = .true.
          end if
                  
--- 312,323 ----
              P%Transfer%high_precision=.true.
              P%Transfer%k_per_logint=3*nint(AccuracyLevel)
            end if
! !addon
! !moved up for screen output
! !          AccuracyBoost = AccuracyLevel
! !          lAccuracyBoost = AccuracyLevel
! !          lSampleBoost = AccuracyLevel
! !end addon
            P%AccurateReionization = .true.
          end if
                  
diff -r -c -b -B -N cosmomc/source/distparams.ini cosmomc_addon/source/distparams.ini
*** cosmomc/source/distparams.ini	1970-01-01 01:00:00.000000000 +0100
--- cosmomc_addon/source/distparams.ini	2006-03-26 18:33:33.000000000 +0100
***************
*** 0 ****
--- 1,132 ----
+ #Params for "getdist" - for processing .txt chain information
+ #Sampple for using 'Kosowsky' parameters theta and Z
+ 
+ columnnum = 22
+ file_root = chains/test
+ out_root = 
+ 
+ #If generated chain at higher temperature can cool for computing results
+ cool = 1
+ 
+ #If 0 assume 1 and no chain filename prefixes
+ chain_num = 4
+ first_chain = 
+ exclude_chain = 
+ 
+ #width of Gaussian smoothing - Should check plots are robust to changes in 
+ #this parameter. Narrow diagonal distributions need larger number
+ num_bins = 20
+ 
+ #For disgarding burn-in if using raw chains
+ #if < 1 interpreted as a fraction of the total number of rows (0.3 ignores first 30% of lines)
+ ignore_rows = 0.5
+ 
+ #if T produced B&W printer friendly output
+ B&W = F
+ #version of MatLab you are using
+ matlab_version = 7
+ 
+ #Include these in 1D plots for comparison - must have same parameters
+ compare_num = 0
+ compare1 = basic6_cmb
+ 
+ plot_meanlikes = T
+ shade_meanlikes = T
+ td_meanlikes = F
+ 
+ # if non-zero, output _thin file, thinned by thin_factor
+ thin_factor = 0
+ #Do probabilistic importance sampling to single samples
+ make_single_samples = F
+ single_thin = 4
+ 
+ #Do simple importance sampling
+ adjust_priors = F
+ #Map parameters to other derived quantities
+ map_params = F
+ 
+ #Use a Gaussian smoothing with width of the bin size
+ #Otherwise use top hat bins
+ smoothing = T
+ 
+ contour1 = 0.68
+ contour2 = 0.95
+ 
+ #if we only want 2D plots agains a particular variable
+ plot_2D_param = 0
+ 
+ #if above zero, instead plot just these combinations:
+ #if both zero it will plot most correlated variables
+ plot_2D_num = 0
+ plot1 = 8 1
+ plot2 =
+ 
+ #number of sample plots, colored by third parameter
+ #if last parameter is 0 or -1 colored by the parameter most correlated
+ #with one of the eigenvector directions (e.g. parallel or orthogonal to degeneracy)
+ num_3D_plots = 1
+ 3D_plot1 = 1 2 3
+ 
+ #Output 2D plots for param combos with 1D marginalized plots along the diagonal
+ triangle_plot = T
+ 
+ #Number of parameters to get covariance matrix for 
+ #If you are going to use the output as a proposal density make sure
+ #you have map_params = F, and the dimension equal to the number of MCMC parameters
+ cov_matrix_dimension = 13
+ 
+ #e.g. colormap('jet')
+ matlab_colscheme =
+ 
+ #Parameters to use. If zero use all parameters which have lables.
+ plotparams_num = 0
+ plotparams = 1 2 3 11
+ 
+ #labels for the parameters
+ lab1 = \Omega_b h^2
+ lab2 = \Omega_c h^2
+ lab3 = \theta
+ lab4 = \tau
+ #Omk
+ lab5 =
+ #nufrac 
+ lab6 = 
+ #w
+ lab7 = 
+ 
+ lab8 = n_s
+ #n_t
+ lab9 =  
+ lab10 = n_run
+ 
+ lab11 = log[10^{10} A_s]
+ #amp ratio
+ lab12 = 
+ lab13 =  
+ 
+ lab14 = \Omega_\Lambda
+ lab15 = Age/GYr
+ lab16 = \Omega_m
+ lab17 = \sigma_8
+ lab18 = z_{re}
+ #r_{10}
+ lab19 = 
+ lab20 = H_0
+ 
+ #Need to give limits if prior cuts off distribution where not very small
+ limits4 = 0.01 N
+ limits6 = 0 N
+ limits7 = -1 N
+ 
+ #compute two-tail marginalized limits irrespective of limits settings above
+ #(otherwise limits are two-tail only for those parameters without limits)
+ force_twotail = F
+ 
+ #PCA - analysis output in file file_root.PCA
+ #number of parameter to do PCA for
+ PCA_num = 0
+ PCA_normparam = 13
+ #The parameters to use
+ PCA_params = 3 13 16
+ #L for log(x), M for log(-x), N for no log
+ PCA_func   = LLLL
\ No newline at end of file
diff -r -c -b -B -N cosmomc/source/driver.F90 cosmomc_addon/source/driver.F90
*** cosmomc/source/driver.F90	2006-03-17 21:11:02.000000000 +0000
--- cosmomc_addon/source/driver.F90	2006-03-26 18:38:43.000000000 +0100
***************
*** 50,61 ****
  
  #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')
  
--- 50,69 ----
  
  #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')
  
***************
*** 77,82 ****
--- 85,95 ----
  
          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 -N cosmomc/source/GetDist.f90 cosmomc_addon/source/GetDist.f90
*** cosmomc/source/GetDist.f90	2006-03-17 16:18:11.000000000 +0000
--- cosmomc_addon/source/GetDist.f90	2006-03-26 18:40:43.000000000 +0100
***************
*** 66,71 ****
--- 66,74 ----
          logical force_twotail 
          logical smoothing, plot_meanlikes, plot_NDcontours
          logical shade_meanlikes, make_single_samples
+ !addon
+         logical td_meanlikes
+ !end addon
          integer single_thin,cust2DPlots(max_cols**2)
          integer ix_min(max_cols),ix_max(max_cols)
          real center(max_cols),width(max_cols)   
***************
*** 1146,1160 ****
   !             write (aunit,fmt) (/(I, I=ix_min(j), ix_max(j), 1)/)*width(j)+center(j)
                if (DoShade) then
                if (shade_meanlikes) then
                  write (aunit,'(a)') 'load (fullfile(''plot_data'',''' // trim(plotfile) //'_likes''));'
                  write (aunit,*) 'contourf(x1,x2,'//trim(plotfile)//'_likes,64);'
                else     
                  write (aunit,*) 'contourf(x1,x2,pts,64);'
                end if
! 
                write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;'
!               end if
!               
             
                if (num_contours /= 0) then
  
--- 1149,1176 ----
   !             write (aunit,fmt) (/(I, I=ix_min(j), ix_max(j), 1)/)*width(j)+center(j)
                if (DoShade) then
                if (shade_meanlikes) then
+ !addon
                  write (aunit,'(a)') 'load (fullfile(''plot_data'',''' // trim(plotfile) //'_likes''));'
+                  if (td_meanlikes) then
+                     write (aunit,*) 'surfl(x1,x2,'//trim(plotfile)//'_likes);'
+                  else
                  write (aunit,*) 'contourf(x1,x2,'//trim(plotfile)//'_likes,64);'
+                  endif
+ !end addon
                else     
                  write (aunit,*) 'contourf(x1,x2,pts,64);'
                end if
! !addon
! !              write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;'
!               if (td_meanlikes) then
!                  write (aunit,*) 'set(gca,''climmode'',''manual'');shading interp; hold on;'
!                  write (aunit,*) 'grid off;box on;camproj(''perspective'');'
!                  write (aunit,*) 'set(gca,''climmode'',''manual'');alpha(0.9);axis tight'
!               else
                write (aunit,*) 'set(gca,''climmode'',''manual''); shading flat; hold on;'
!               endif
! !end addon
!               endif
             
                if (num_contours /= 0) then
  
***************
*** 1472,1477 ****
--- 1488,1496 ----
            write (*,*) 'WARNING: Mapping params - .covmat file is new params.'
  
      shade_meanlikes = Ini_Read_Logical('shade_meanlikes',.false.)   
+ !addon
+     td_meanlikes = Ini_Read_Logical('td_meanlikes',.false.)  
+ !end addon
  
      matlab_col = Ini_Read_String('matlab_colscheme')
      smoothing = Ini_Read_Logical('smoothing',.true.)
***************
*** 2004,2010 ****
--- 2023,2036 ----
        write (51,*)  'set(gcf, ''PaperUnits'',''inches'');'
       
        if (plot_row < plot_col .and. plot_row*plot_col>9) then
+ !addon
+ !            write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 10]);';
+          if (td_meanlikes) then
+             write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 7 9]);';
+          else
              write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 10]);';
+          endif
+ !end addon           
            else
              write (51,*) 'set(gcf, ''PaperPosition'',[ 0 0 8 8]);';
        end if
diff -r -c -b -B -N cosmomc/source/Makefile cosmomc_addon/source/Makefile
*** cosmomc/source/Makefile	2006-03-17 22:10:50.000000000 +0000
--- cosmomc_addon/source/Makefile	2006-03-26 18:42:02.000000000 +0100
***************
*** 1,160 ****
! #Settings here are for Intel ifort on CITA cluster; uncomment others below
! #You may need to edit the library paths for MKL for Intel
! #Beware of using optmizations which lose accuracy - may give errors when running
! #With Visual Fortan open cosmomc.dsw workspace file instead.
! 
! #Edit for CFITSIO, HEALPIX directories
! cfitsio = /home/pascal/include/cfitsio
! healpix = /home/pascal/include/Healpix_2.00
! WMAP3 = /home/antlewis/WMAP3
! 
! #Intel MPI (assuming mpif77 set to point to ifc)
! #change -lmkl_ia32 to -lmkl_p3 for MKL versions earlier than 6 (6 needed for ifc 8+) 
! F90C     = mpif77
! FFLAGS = -O2 -Vaxlib -ip -W0 -WB -openmp -fpp -DMPI
! LAPACKL = -L/opt/intel/mkl61/lib/32 -lmkl_lapack -lmkl_ia32 -lguide -lpthread
! INCLUDE = -I../camb 
  
! # Use with Intel 9 on the Altix (eg. COSMOS)
! # (do "module load icomp90" before compiling, or add to .bashrc)
! #F90C = ifort
! #FFLAGS = -O3 -w -fpp2 -DMPI 
! #LAPACKL = -lscs -lguide -lmpi -ldl
! #INCLUDE = -I../camb
! 
! #Intel fortran 8, check you have the latest update from the Intel web pages
! #See Makefile_intel for ifc 7.1 or lower (some versions have problems)
! #F90C     = ifort
! #FFLAGS = -O2 -Vaxlib -ip -W0 -WB -openmp -fpp
! #LAPACKL = -L/opt/intel/mkl61/lib/32 -lmkl_lapack -lmkl_ia32 -lguide -lpthread
! #INCLUDE = -I../camb 
! 
! #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
! #INCLUDE = -I../camb
! 
! #SGI, -mp toggles multi-processor. Use -O2 if -Ofast gives problems.
! #Not various versions of the compiler are buggy giving erroneous seg faults with -mp.
! #Version 7.3 is OK, currently version 7.4 is bugged, as are some earlier versions.
! #F90C     = f90
  #LAPACKL = -lcomplib.sgimath
! #INCLUDE = -I../camb 
! #FFLAGS  = -Ofast -mp
  
! #Digital/Compaq fortran, -omp toggles multi-processor
! #F90C    = f90
! #FFLAGS  = -omp -O4 -arch host -math_library fast -tune host -fpe1
! #LAPACKL  = -lcxml
! #INCLUDE = -I../camb 
! 
! #Intel Itanium efc, -openmp toggles multi-processor:
! #F90C     = efc
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp
! #LAPACKL = -L/opt/intel/mkl/lib/64 -lmkl_lapack -lmkl_ipf -lguide -lpthread -static
! #INCLUDE = -I../camb 
! 
! #Intel Itanium efc MPI:
! # mpirun -np 4 ./cosmomc params.ini
! #F90C     = efc
! #FFLAGS = -O2 -Vaxlib -W0 -WB -openmp -fpp -lmpi -lscs  -ldl -DMPI
! #LAPACKL = -L/opt/intel/mkl/lib/64 -lmkl_lapack -lmkl_ipf -limf
! #INCLUDE = -I../camb
! 
! #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
! #INCLUDE = -I../camb
! 
! #Absoft ProFortran, single processor, set -cpu:[type] for your local system 
! #F90C     = f95
! #FFLAGS = -O2 -s -cpu:athlon -lU77 -w -YEXT_NAMES="LCS" -YEXT_SFX="_" 
! #LAPACKL =  -llapack -lblas -lg2c 
! #INCLUDE = -I../camb -p../camb
! 
! #NAGF95, single processor:
! #F90C     = f95
! #FFLAGS = -DNAGF95 -O3
! #LAPACKL = -llapack -lblas -lg2c
! #INCLUDE = -I../camb 
! 
! #Sun, single processor:
! #F90C     = f90
! #FFLAGS = -O5
! #LAPACKL = -lsunperf -lfsu
! #INCLUDE = -I../camb -M../camb
! 
! #Sun MPI
! #F90C = mpf90
! #FFLAGS =  -O4 -openmp -ftrap=%none -dalign -DMPI
! #LAPACKL = -lsunperf -lfsu -lmpi_mt
! #INCLUDE = -I../camb -M../camb
! 
! #Sun parallel enterprise:
! #F90C     = f95
! #FFLAGS =  -O4 -xarch=native64 -openmp -ftrap=%none
! #LAPACKL = -lsunperf -lfsu
! #INCLUDE = -I../camb -M../camb
! 
! #IBM XL Fortran, multi-processor (run "module load lapack" then run "gmake")
! # 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 = -I../camb -lessl
! #LAPACKL =
  
  PROPOSE = propose.o
  CLSFILE = CMB_Cls_simple.o
  
  #Can use params_H if you prefer more generic parameters
  PARAMETERIZATION = params_CMB.o
  
! WMAP3_inc = -I$(cfitsio)/include -I$(healpix)/include -I$(WMAP3)
! WMAP3_lib = -L$(cfitsio)/lib -L$(healpix)/lib -L$(WMAP3) \
! 	-lhealpix -lcfitsio
  
! CLSLIB  = $(WMAP3_lib) -L../camb -lcamb $(LAPACKL) 
  
! F90FLAGS = $(FFLAGS) $(WMAP3_inc) $(INCLUDE)
  
  DISTFILES = utils.o settings.o GetDist.o
  
! WMAPOBJS = $(WMAP3)/read_archive_map.o \
! 	$(WMAP3)/read_fits.o \
! 	$(WMAP3)/WMAP_3yr_options.o \
! 	$(WMAP3)/WMAP_3yr_util.o \
! 	$(WMAP3)/WMAP_3yr_tt_pixlike.o \
! 	$(WMAP3)/WMAP_3yr_tt_beam_and_ptsrc_corr.o \
! 	$(WMAP3)/WMAP_3yr_teeebb_pixlike.o \
! 	$(WMAP3)/WMAP_3yr_likelihood.o
  
  OBJFILES= $(WMAPOBJS) utils.o settings.o cmbtypes.o cmbdata.o WeakLen.o twodf.o mpk.o \
  	supernovae.o lya.o\
  	$(CLSFILE) paramdef.o $(PROPOSE) $(PARAMETERIZATION) calclike.o \
  	conjgrad_wrapper.o EstCovmat.o postprocess.o MCMC.o driver.o
  
! default: cosmomc
  
! all : cosmomc getdist
  
  .f.o:
  	f77 $(F90FLAGS) -c $<
  
  %.o: %.f90
! 	$(F90C) $(F90FLAGS) -c $*.f90
  
  %.o: %.F90
! 	$(F90C) $(F90FLAGS) -c $*.F90
! 
  
- cosmomc: $(OBJFILES) 	
- 	$(F90C) -o ../cosmomc $(OBJFILES) $(CLSLIB) $(F90FLAGS)
  
  
  clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist: $(DISTFILES)
! 	$(F90C) -o ../getdist $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
--- 1,102 ----
! # >>> 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
! LAPACKL = -lcomplib.sgimath_mp
! INCLUDE = -I../camb
! CFITSIODIR = /people/ringeval/usr
! HEALPIXDIR = /people/ringeval/usr
! endif
! 
! ifeq ($(ext),Lin)
! F90C=ifort
! FFLAGS= -O -openmp -fpp
! WMAPFLAGS= -O -save -fpp
! LAPACKL = -lmkl -lmkl_lapack
! INCLUDE = -I../camb 
! CFITSIODIR = /opt/intel/cfitsio-3.006
! HEALPIXDIR = /opt/intel/healpix-2.0.1
! endif 
! 
! ifeq ($(ext),OSF)
! F90C=f90
! FFLAGS= -omp -O -arch host -math_library fast -tune host -fpe1
! WMAPFLAGS= $(FFLAGS)
! LAPACKL  = -lcxml
! INCLUDE = -I../camb 
! endif
  
! ifeq ($(ext),Sun)
! F90C=f90
! FFLAGS= -O4 -xarch=native64 -openmp -ftrap=%none
! WMAPFLAGS= $(FFLAGS)
! LAPACKL = -lsunperf -lfsu
! INCLUDE = -I../camb -M../camb
! endif
! 
! ifeq ($(ext),AIX)
! F90C    = mpxlf90_r
! FFLAGS  = -O4 -WF,-DIBMXL,-DMPI -qstrict -qsmp=omp -qmaxmem=-1 -qsuffix=f=f90:cpp=F90
! WMAPFLAGS= $(FFLAGS)
! INCLUDE = -I../camb -lessl
! LAPACKL =
! endif
  
+ WMAPDIR = ../WMAP
  PROPOSE = propose.o
  CLSFILE = CMB_Cls_simple.o
  
  #Can use params_H if you prefer more generic parameters
  PARAMETERIZATION = params_CMB.o
  
! WMAPINCLUDE =  -I$(CFITSIODIR)/include -I$(HEALPIXDIR)/include -I$(WMAPDIR)
  
! CLSLIB  = -L$(HEALPIXDIR)/lib -L$(CFITSIODIR)/lib \
! 	-L../camb -lcamb -lhealpix -lcfitsio $(LAPACKL)
  
! F90FLAGS = $(FFLAGS) $(INCLUDE)
  
  DISTFILES = utils.o settings.o GetDist.o
  
! WMAPOBJS = read_archive_map.o read_fits.o WMAP_3yr_options.o WMAP_3yr_util.o \
! 	WMAP_3yr_tt_pixlike.o WMAP_3yr_tt_beam_and_ptsrc_corr.o \
! 	WMAP_3yr_teeebb_pixlike.o WMAP_3yr_likelihood.o
  
  OBJFILES= $(WMAPOBJS) utils.o settings.o cmbtypes.o cmbdata.o WeakLen.o twodf.o mpk.o \
  	supernovae.o lya.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)
  
  .f.o:
  	f77 $(F90FLAGS) -c $<
  
+ %.o: $(WMAPDIR)/%.f90
+ 	$(F90C) $(WMAPFLAGS) $(WMAPINCLUDE) -c $<
+ 
  %.o: %.f90
! 	$(F90C) $(F90FLAGS) -c $<
  
  %.o: %.F90
! 	$(F90C) $(F90FLAGS) -c $<
  
  
+ cosmomc.$(ext): $(OBJFILES) 	
+ 	$(F90C) -o ../$@ $(OBJFILES) $(CLSLIB) $(F90FLAGS)
  
  clean:
  	rm -f *.o *.mod *.d *.pc *.obj ../core
  
! getdist.$(ext): $(DISTFILES)
! 	$(F90C) -o ../$@ $(DISTFILES) $(CLSLIB) $(F90FLAGS) 
diff -r -c -b -B -N cosmomc/source/MCMC.f90 cosmomc_addon/source/MCMC.f90
*** cosmomc/source/MCMC.f90	2006-03-17 16:18:11.000000000 +0000
--- cosmomc_addon/source/MCMC.f90	2006-03-26 19:31:09.000000000 +0100
***************
*** 313,320 ****
      end if
  
      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
   
--- 313,322 ----
      end if
  
      Like = GetLogLike(Trial) 
! !addon    if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:', CurLike
!         if (Feedback > 1) write (*,*) 'Likelihood: ', Like, 'Current Like:'&
!              , CurLike, 'Rank', MPIRank
! !end addon
      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 -N cosmomc/source/params.ini cosmomc_addon/source/params.ini
*** cosmomc/source/params.ini	1970-01-01 01:00:00.000000000 +0100
--- cosmomc_addon/source/params.ini	2006-03-26 18:44:37.000000000 +0100
***************
*** 0 ****
--- 1,183 ----
+ #Sample parameters for cosmomc in default parameterization
+ 
+ #Root name for files produced
+ file_root = chains/test
+ 
+ #action = 0:  MCMC, action=1: postprocess .data file, action=2: find best fit point only
+ action = 0
+ 
+ #Maximum number of chain steps
+ samples = 200000
+ 
+ #Feedback level ( 2=lots,1=chatty,0=none)
+ feedback = 1
+ 
+ #Temperature at which to Monte-Carlo
+ temperature = 1
+ 
+ #filenames for CMB datasets
+ #Note you may need to change lmax in cmbtypes.f90 to use small scales (e.g. lmax=2100)
+ cmb_numdatasets = 1
+ cmb_dataset1 = WMAP
+ 
+ cmb_dataset2 = data/ACBAR_lge800_.dataset
+ cmb_dataset3 = data/CBI2.0_lge600.dataset
+ cmb_dataset4 = data/VSA_lge600.dataset
+ 
+ #filenames for matter power spectrum datasets, incl twodf
+ mpk_numdatasets = 1
+ mpk_dataset1 = twodf
+ #mpk_dataset1 = data/sdss.dataset
+ 
+ #if true, use HALOFIT for non-linear corrections (astro-ph/0207664).
+ #note lyman-alpha (lya) code assumes linear spectrum
+ nonlinear_pk = F
+ 
+ use_CMB = T
+ use_HST = T
+ use_mpk = F
+ use_clusters = F
+ use_BBN = F
+ use_Age_Tophat_Prior = T
+ use_SN = F
+ use_lya = F
+ use_min_zre = 0
+ 
+ #Force computation of sigma_8 even if use_mpk = F
+ get_sigma8 = F
+ 
+ #1: Simple Metropolis, 2: slice sampling, 3: slice sampling fast parameters, 4: directional gridding
+ sampling_method = 1
+ 
+ #if sampling_method =4, iterations per gridded direction
+ directional_grid_steps = 20
+ 
+ #use fast-slow parameter distinctions to speed up
+ use_fast_slow = T
+ 
+ #Can use covariance matrix for proposal density, otherwise use settings below
+ #Covariance matrix can be produced using "getdist" program.
+ propose_matrix = params_CMB.covmat
+ 
+ #If propose_matrix is blank (first run), can try to use numerical Hessian to 
+ #estimate a good propose matrix. As a byproduct you also get an approx best fit point
+ estimate_propose_matrix = F
+ 
+ #Tolerance on log likelihood to use when estimating best fit point
+ delta_loglike = 2
+ 
+ #Scale of proposal relative to covariance; 2.4 is recommended by astro-ph/0405462 for Gaussians
+ #If propose_matrix is much broader than the new distribution, make proportionately smaller
+ #Generally make smaller if your acceptance rate is too low
+ propose_scale = 2.4
+ 
+ #Increase to oversample fast parameters more, e.g. if space is odd shape
+ oversample_fast = 1
+ 
+ #if non-zero number of steps between sample info dumped to file file_root.data
+ indep_sample = 0
+ 
+ #number of samples to disgard at start; usually set to zero and remove later
+ burn_in = 0
+ 
+ #If zero set automatically
+ num_threads = 0
+ 
+ #MPI mode multi-chain options (recommended)
+ #MPI_Converge_Stop is a (variance of chain means)/(mean of variances) parameter that can be used to stop the chains
+ #Set to a negative number not to use this feature. Does not guarantee good accuracy of confidence limits.
+ MPI_Converge_Stop = 0.03
+ 
+ #Do initial period of slice sampling; may be good idea if 
+ #cov matrix or widths are likely to be very poor estimates
+ MPI_StartSliceSampling  = F
+ 
+ #Can optionally also check for convergence of confidence limits (after MPI_Converge_Stop reached)
+ #Can be good idea as small value of MPI_Converge_Stop does not (necessarily) imply good exploration of tails
+ MPI_Check_Limit_Converge = F
+ 
+ #if MPI_Check_Limit_Converge = T, give tail fraction to check (checks both tails):
+ MPI_Limit_Converge = 0.025
+ #permitted quantile chain variance in units of the standard deviation (small values v slow):
+ MPI_Limit_Converge_Err = 0.3
+ #which parameter's tails to check. If zero, check all parameters:
+ MPI_Limit_Param = 0
+ 
+ #if MPI_LearnPropose = T, the proposal density is continally updated from the covariance of samples so far (since burn in)
+ MPI_LearnPropose = T
+ #can set a value of converge at which to stop updating covariance (so that it becomes rigorously Markovian)
+ #e.g. MPI_R_StopProposeUpdate = 0.4 will stop updating when (variance of chain means)/(mean of variances) < 0.4
+ MPI_R_StopProposeUpdate = 0
+ 
+ 
+ #if blank this is set from system clock
+ 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 = T
+ 
+ 
+ #CAMB parameters
+ #If we are including tensors
+ compute_tensors = F
+ #Initial power spectrum amplitude point (Mpc^{-1})
+ pivot_k = 0.05
+ #If using tensors, enforce n_T = -A_T/(8A_s)
+ inflation_consistency = F
+ 
+ #Whether the CMB should be lensed (slows a lot unless also computing matter power)
+ CMB_lensing = F
+ #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
+ redo_theory = F
+ redo_cls = F
+ redo_pk = F
+ redo_skip = 0
+ redo_outroot = 
+ redo_thin = 1
+ redo_add = F
+ redo_from_text = F
+ #If large difference in log likelihoods may need to offset to give sensible weights
+ #for exp(difference in likelihoods)
+ redo_likeoffset =  0
+ 
+ #parameter start center, min, max, start width, st. dev. estimate
+ #ombh2
+ param1 = 0.0223 0.005 0.1 0.001 0.001
+ #omdmh2
+ param2 = 0.105 0.01 0.99 0.01 0.01
+ #theta
+ param3 = 1.04 0.5 10 0.002 0.002
+ #tau
+ param4 = 0.09 0.01 0.8 0.03 0.03
+ 
+ #omk
+ param5 = 0 0 0 0 0
+ #nufrac
+ param6 = 0 0 0 0 0
+ #w
+ param7 = -1 -1 -1 0 0
+ 
+ #n_s
+ param8 = 0.95 0.5 1.5 0.03 0.01
+ #n_t
+ param9 = 0 0 0 0 0
+ #n_run
+ param10 = 0 0 0 0 0
+ 
+ #log[10^10 A_s]
+ param11 = 3 2.7 4 0.01 0.01
+ #amp_ratio
+ param12 = 0 0 0 0 0
+ #unused fast parameter
+ param13 = 0 0 0 0 0
+ 
diff -r -c -b -B -N cosmomc/source/settings.f90 cosmomc_addon/source/settings.f90
*** cosmomc/source/settings.f90	2006-03-17 16:18:11.000000000 +0000
--- cosmomc_addon/source/settings.f90	2006-02-20 17:27:32.000000000 +0000
***************
*** 5,10 ****
--- 5,15 ----
    implicit none
  
    real :: AccuracyLevel = 1.
+ !addon
+   real :: lAccuracyLevel = 1.
+   real :: lSampleLevel = 1.
+   logical :: CambComputeEveryl = .false.
+ !end addon
    !Set to >1 to use CAMB etc on higher accuracy settings. 
    !Does not affect MCMC (except making it all slower)
  
***************
*** 43,48 ****
--- 48,56 ----
  
    integer :: num_threads = 0
    integer :: instance = 0
+ !addon
+   integer :: instance_shift = 0
+ !end addon
    integer :: MPIchains = 1, MPIrank = 0
  
    logical :: Use_LSS = .true.
