!Module with subroutines for generating Cls. AL 30/3/01.
!
!To use first call Initialize(name_of_flat_bessels,max_multipole,Max_wavenumber)
!For each model fill in a structure of type CAMBparams and call GetCls(my_CAMBparams)
!The Cls will then be in the ModelData (modules.f90) module, clts for Cl,temp,scalar, etc.

        module CAMBinterface
         use ModelData
         implicit none
!Possible values for initial conditions
         integer, parameter :: initial_adiabatic=1,initial_iso_CDM=2,  &
                               initial_iso_baryon=3,initial_seed=4

        type CAMBparams
    
         logical   :: Out_scalar,Out_tensor

         integer   :: Max_l
         real*8    :: Max_k

         real*8    :: Om_baryon,Om_cdm,Om_lambda,Om_nu
         real*8    :: Hubble_0,T_CMB,Y_He,Num_Nu_massless,Num_Nu_massive

!If you write your own Powers module these will be ignored
         real*8    :: Scalar_spectral_index,Tensor_spectral_index
         real*8    :: Tensor_scalar_ratio !this is meaningless unless flat SCDM model

         integer   :: Scalar_initial_condition 
!must be one of values above
     
         logical   :: RECFAST_recombination

!Reionization settings - used if Reionization=.true.
!if  Reionization_use_optical_depth=.false.  then specify redshift and xe instead
         logical   :: Reionization,  Re_use_optical_depth 
         real*8    :: Re_optical_depth
!redshift, ionization fraction (e.g. 50 0.2)
         real*8    :: Re_redshift,Re_ion_fraction

        end type CAMBparams

!      itflag - 0 scal only, 1 tens+scal, 2 tens only
        integer itflag
        
        integer rcflag !1 use RECFAST
        integer lmaxinitial
        real*8 MaxWavenumberInitial

     contains
     
       subroutine Initialize(BesselsFilename,MaxMultipole,MaxWavenumber)
         use lvalues
         use SpherBessels
         character (LEN=80), intent(in) ::  BesselsFilename
         integer, intent(in)      ::  MaxMultipole
         real*8, intent(in)       ::  MaxWavenumber
       
         lmo = MaxMultipole
         akmax0 = MaxWavenumber 
         
         if (.not.allocated(ajl).or. lmo>lmaxinitial.or. MaxWavenumber>MaxWavenumberInitial)then
         lmaxinitial=lmo
         MaxWavenumberInitial=akmax0
         
         call InitSpherBessels(BesselsFilename) 
         end if
     
       end subroutine Initialize

 
       subroutine GetCls(Params)
        use ModelParams
        use ModelData
        use Reionization
        use Transfer
        use SpherBessels
        use USpherBessels
        use lvalues
        use InitialPower
        implicit none

        type(CAMBparams), intent(in) :: Params

        integer riflag
        
        lmo = Params%Max_l
        akmax0 = Params%Max_k


        WantCls=.true.
        WantTransfer=.false.
        WantVars=.false.
        ictpres=0
        
        ntf=0
      
        omegab=Params%Om_baryon
        omegac=Params%Om_cdm
        omegav=Params%Om_lambda
        omegan=Params%Om_nu

         h0=Params%Hubble_0
        tcmb=Params%T_CMB
        yhe=Params%Y_He
        annur=Params%Num_Nu_massless
        annunr=Params%Num_Nu_massive

        akmaxt=akmaxt*(h0/100.0d0)
        if (h0 < 25.d0.or.h0 > 100.d0) then
          write(*,*) &
            '  Warning: Hubble_ has units of km/s/Mpc.  Your value is weird.'
        end if
        if (tcmb < 2.7d0.or.tcmb > 2.8d0) then
          write(*,*) &
            '  Warning: T_CMB has units of K.  Your value is weird.'
        end if

        if (yhe < 0.2d0.or.yhe > 0.3d0) then
          write(*,*) &
            '  Warning: Y_He is the Helium fraction of baryons.', &
            '  Your value is weird.'
        end if
        if (annunr < 0.or.annunr > 3.1) then
          write(*,*) &
            'Warning: Num_Nu_massive is strange'
          write(*,*) '  Illegal choice. '
      
        end if
        if (annur < 0.or.annur > 3.1) then
          write(*,*) &
            'Warning: Num_Nu_massive is strange'
          write(*,*) '  Illegal choice.'
        end if
        if (annunr < 1.and.omegan > 0.0) then
          write(*,*) &
            'Warning: Num_nu_massive should be 1, 2, or 3', &
            'For non zero omegan'
           write(*,*) '  Illegal choice. '
        end if
       
        UseRECFAST = Params%RECFAST_recombination
        
        riflag=0
        if (Params%Reionization) then
           if (Params%Re_use_optical_depth) then
              riflag = 1
           else
              riflag =2
           end if
        end if

        zri=0.0d0
        rif=0.0d0
        optdlss=0.0d0

        if (riflag == 1) then
           optdlss=  Params%Re_optical_depth
           rif=1.0d0
        end if

        if (riflag == 2) then
           zri=Params%Re_redshift
           rif=Params%Re_ion_fraction
           zristp=0.07d0*zri-1.0d0
           if (zristp < 0.0) zristp=0.0d0
        end if
        
           rat(1)=1.d0
         
           itflag=0
            nn=1
            !just use one initial power spectrum for the moment
           if (Params%Out_scalar) then
              an(1)=Params%Scalar_spectral_index
           end if
           if (Params%Out_tensor) then
              if (Params%Out_scalar) then 
                 itflag=1 
              else 
                 itflag =2
              end if
           
              ant(1)=Params%Tensor_spectral_index
              rat(1)=Params%Tensor_scalar_ratio
           end if

        
          lensflag=0
         
        if (itflag /= 2) then
           initfl=Params%Scalar_initial_condition
        else
           initfl=1
        end if
       
        WantScalars=itflag /= 2
        WantTensors=itflag /= 0
        
        call CalcSecondaryVars !set other derived variables in ModelParams (modules.f90)

        if (flat.and..not.allocated(ajl)) then
           write(*,*) 'For flat models you must call Initialize first.'
           stop
        end if
       
        call initlval !Needed before InitSpherBessels
         
        call cmbmain
    
        if (OutputNormalization == outCOBE) call COBEnormalize
  
     end subroutine GetCls

 end module CAMBinterface