Add: After use Errors
! EFTCAMB MOD START
use EFTinitialization
! EFTCAMB MOD END
Replace:
! EFTCAMB MOD START:
! To debug the code it should be launched on just one thread by
! changing the appropriate entry in the parameter file.
real(dl) :: fixq = 0._dl !Debug output of one q
! EFTCAMB MOD END.
#… Original:
real(dl) :: fixq = 0._dl !Debug output of one q
::::: subroutine cmbmain :::::
Add: After integer q_ix
! EFTCAMB MOD START
logical EFTsuccess
! EFTCAMB MOD END
Add: After WantLateTime = CP%DoLensing .or. num_redshiftwindows > 0
! EFTCAMB MOD START
if (EFTCAMBuseinCOSMOMC==1.and.CP%EFTflag/=0) then
call EFTCAMB_initialization(EFTsuccess)
if (.not.EFTsuccess) stop
end if
! EFTCAMB MOD END
::::: subroutine CalcScalarSources(EV,taustart) :::::
Replace:
! EFTCAMB MOD START: It's always nice to plot the evolution of perturbations in DE/MG models!
! In order to print the user needs to choose the scale by means of fixq and then turno on the
! flag that prints quantities from the output routine.
! It is suggested to run the code on just one core.
if (fixq/=0._dl) then
tol1=tol/exp(AccuracyBoost-1)
write(*,*) 'EFTCAMB: start printing.'
call CreateTxtFile('Results/Debug_Evolution/Files/1_FRW.dat',1)
call CreateTxtFile('Results/Debug_Evolution/Files/2_EFTfunctions.dat',2)
call CreateTxtFile('Results/Debug_Evolution/Files/3_EFTBackground.dat',3)
call CreateTxtFile('Results/Debug_Evolution/Files/4_PiFieldSolution.dat',4)
call CreateTxtFile('Results/Debug_Evolution/Files/5_MetricSolution.dat',5)
call CreateTxtFile('Results/Debug_Evolution/Files/6_DensitySolution.dat',6)
call CreateTxtFile('Results/Debug_Evolution/Files/7_EinsteinEqFactors.dat',7)
call CreateTxtFile('Results/Debug_Evolution/Files/8_PiEqFactors.dat',8)
call CreateTxtFile('Results/Debug_Evolution/Files/9_Sources.dat',9)
call CreateTxtFile('Results/Debug_Evolution/Files/10_MuGamma.dat',10)
call CreateTxtFile('Results/Debug_File/Debug.dat',11)
do j=1,10000
tauend = taustart +REAL(j-1)*(CP%tau0-taustart)/REAL(10000-1)
call GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,c,w)
yprime = 0
call derivs(EV,EV%ScalEqsToPropagate,tau,y,yprime)
call output(EV,y,j,tau,sources)
end do
close(1)
close(2)
close(3)
close(4)
close(5)
close(6)
close(7)
close(8)
close(9)
close(10)
close(11)
write(*,*) 'Stop printing'
stop
end if
! Original code:
!!Example code for plotting out variable evolution
!!if (fixq/=0._dl) then
!! tol1=tol/exp(AccuracyBoost-1)
!! call CreateTxtFile('evolve_q005.txt',1)
!! do j=1,1000
!! tauend = taustart+(j-1)*(CP%tau0-taustart)/1000
!! call GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,c,w)
!! yprime = 0
!! call derivs(EV,EV%ScalEqsToPropagate,tau,y,yprime)
!! adotoa = 1/(y(1)*dtauda(y(1)))
!! ddelta= (yprime(3)*grhoc+yprime(4)*grhob)/(grhob+grhoc)
!! delta=(grhoc*y(3)+grhob*y(4))/(grhob+grhoc)
!! growth= ddelta/delta/adotoa
!! write (1,'(7E15.5)') tau, delta, growth, y(3), y(4), y(EV%g_ix), y(1)
!! end do
!! close(1)
!! stop
!!end if
! EFTCAMB MOD END
::::: subroutine GetTransfer(EV,tau) :::::
Replace:
! EFTCAMB MOD START: this seems to create some problems. It should be safe to turn it off.
if (CP%EFTflag/=0) then
if (CP%Transfer%high_precision) atol=atol
else
if (CP%Transfer%high_precision) atol=atol/10000
end if
! Original code:
! if (CP%Transfer%high_precision) atol=atol/10000
! EFTCAMB MOD END.
----------------------------------------------------------------------------------------------------------------------------
Modify:
modules.f90
::::: module ModelParams :::::
Add: After integer :: Nu_mass_numbers(max_nu) !physical number per eigenstate
! EFTCAMB MOD START
! 1) Definition of flags:
integer :: EFTflag
integer :: EFTwDE
integer :: PureEFTmodelOmega
integer :: PureEFTmodelAlpha1, PureEFTmodelAlpha2, PureEFTmodelAlpha3
integer :: PureEFTmodelAlpha4, PureEFTmodelAlpha5, PureEFTmodelAlpha6
integer :: DesignerEFTmodel
integer :: MatchingEFTmodel
! 2) Definition of model parameters:
real(dl) :: EFTw0, EFTwa, EFTwn, EFTwat, EFtw2, EFTw3
real(dl) :: EFTOmega0, EFTOmegaExp
real(dl) :: EFTAlpha10, EFTAlpha1Exp, EFTAlpha20, EFTAlpha2Exp
real(dl) :: EFTAlpha30, EFTAlpha3Exp, EFTAlpha40, EFTAlpha4Exp
real(dl) :: EFTAlpha50, EFTAlpha5Exp, EFTAlpha60, EFTAlpha6Exp
real(dl) :: EFTB0
! EFTCAMB MOD END
::::: module MassiveNu :::::
Replace:
! EFTCAMB MOD START: compatibility with massive neutrinos
real(dl), dimension(:), allocatable :: r1,p1,dr1,dp1,ddr1, ddp1, dddp1
! Original code:
! real(dl), dimension(:), allocatable :: r1,p1,dr1,dp1,ddr1
! EFTCAMB MOD END.
Replace:
! EFTCAMB MOD START: compatibility with massive neutrinos
public const,Nu_Init,Nu_background, Nu_rho, Nu_drho, nqmax0, nqmax, &
nu_int_kernel, nu_q, Nu_pidot, Nu_pidotdot
! Original code:
! public const,Nu_Init,Nu_background, Nu_rho, Nu_drho, nqmax0, nqmax, &
! nu_int_kernel, nu_q
! EFTCAMB MOD END.
::::: subroutine Nu_init :::::
Replace:
! EFTCAMB MOD START: compatibility with massive neutrinos
allocate(r1(nrhopn),p1(nrhopn),dr1(nrhopn),dp1(nrhopn),ddr1(nrhopn),ddp1(nrhopn),dddp1(nrhopn))
! Original code:
! allocate(r1(nrhopn),p1(nrhopn),dr1(nrhopn),dp1(nrhopn),ddr1(nrhopn))
! EFTCAMB MOD END.
Add: After call splder(dr1,ddr1,nrhopn,spline_data)
! EFTCAMB MOD START: compatibility with massive neutrinos
call splder(dp1,ddp1,nrhopn,spline_data)
call splder(ddp1,dddp1,nrhopn,spline_data)
! EFTCAMB MOD END.
Add: After end function Nu_drho
! EFTCAMB MOD START: compatibility with massive neutrinos
function Nu_pidot(am,adotoa,presnu) result (presnudot)
use precision
use ModelParams
real(dl) adotoa,presnu,presnudot
real(dl) d
real(dl), intent(IN) :: am
integer i
if (am< am_minp) then
presnudot = -2*const2*am**2*adotoa/3._dl
else if (am>am_maxp) then
presnudot = -((15._dl*(4._dl*am**2*zeta5 -189._dl*Zeta7))/(8._dl*am**3*const))*adotoa
else
d=log(am/am_min)/dlnam+1._dl
i=int(d)
d=d-i
presnudot = dp1(i)+d*(ddp1(i)+d*(3._dl*(dp1(i+1)-dp1(i))-2._dl*ddp1(i) &
-ddp1(i+1)+d*(ddp1(i)+ddp1(i+1)+2._dl*(dp1(i)-dp1(i+1)))))
presnudot=presnu*adotoa*presnudot/dlnam
end if
end function Nu_pidot
function Nu_pidotdot(am,adotoa,Hdot,presnu,presnudot) result (presnudotdot)
use precision
use ModelParams
real(dl) adotoa,Hdot,presnu,presnudot,presnudotdot
real(dl) d
real(dl), intent(in) :: am
integer i
if (am< am_minp) then
presnudotdot = presnudot*(adotoa +Hdot/adotoa) +am**2*adotoa**2*(-2._dl*const2/3._dl)
else if (am>am_maxp) then
presnudotdot = presnudot*(adotoa +Hdot/adotoa) +am**2*adotoa**2*(&
&-((15._dl*zeta5)/(am**3*const)) + (15._dl*(4._dl*am**2*zeta5 -189._dl*Zeta7))/(2._dl*am**5*const))
else
d=log(am/am_min)/dlnam+1._dl
i=int(d)
d=d-i
presnudotdot = ddp1(i)+d*(dddp1(i)+d*(3._dl*(ddp1(i+1)-ddp1(i))-2._dl*dddp1(i) &
-dddp1(i+1)+d*(dddp1(i)+dddp1(i+1)+2._dl*(ddp1(i)-ddp1(i+1)))))
presnudotdot = +adotoa**2*presnu*presnudotdot/dlnam +Hdot/adotoa*presnudot +presnudot**2/presnu
end if
end function Nu_pidotdot
! EFTCAMB MOD END.
----------------------------------------------------------------------------------------------------------------------------
Modify:
inidriver.F90
::::: program driver :::::
Add: After
#ifdef NAGF95
use F90_UNIX
#endif
! EFTCAMB MOD START
use EFTdef
! EFTCAMB MOD END
Add: After P%tcmb = Ini_Read_Double('temp_cmb',COBE_CMBTemp)
! EFTCAMB MOD START
! 1) Initialization of EFTCAMB flags.
P%EFTflag = Ini_Read_Int('EFTflag',0)
P%EFTwDE = Ini_Read_Int('EFTwDE',0)
P%PureEFTmodelOmega = Ini_Read_Int('PureEFTmodelOmega',0)
P%PureEFTmodelAlpha1 = Ini_Read_Int('PureEFTmodelAlpha1',0)
P%PureEFTmodelAlpha2 = Ini_Read_Int('PureEFTmodelAlpha2',0)
P%PureEFTmodelAlpha3 = Ini_Read_Int('PureEFTmodelAlpha3',0)
P%PureEFTmodelAlpha4 = Ini_Read_Int('PureEFTmodelAlpha4',0)
P%PureEFTmodelAlpha5 = Ini_Read_Int('PureEFTmodelAlpha5',0)
P%PureEFTmodelAlpha6 = Ini_Read_Int('PureEFTmodelAlpha6',0)
P%DesignerEFTmodel = Ini_Read_Int('DesignerEFTmodel',1)
P%MatchingEFTmodel = Ini_Read_Int('MatchingEFTmodel',1)
! 2) Initialization of EFTCAMB model parameters.
P%EFTw0 = Ini_Read_Double('EFTw0',-1._dl)
P%EFTwa = Ini_Read_Double('EFTwa',0._dl)
P%EFTwn = Ini_Read_Double('EFTwn',2._dl)
P%EFTwat = Ini_Read_Double('EFTwat',1._dl)
P%EFtw2 = Ini_Read_Double('EFtw2',0._dl)
P%EFTw3 = Ini_Read_Double('EFTw3',0._dl)
P%EFTOmega0 = Ini_Read_Double('EFTOmega0')
P%EFTOmegaExp = Ini_Read_Double('EFTOmegaExp')
P%EFTAlpha10 = Ini_Read_Double('EFTAlpha10')
P%EFTAlpha1Exp = Ini_Read_Double('EFTAlpha1Exp')
P%EFTAlpha20 = Ini_Read_Double('EFTAlpha20')
P%EFTAlpha2Exp = Ini_Read_Double('EFTAlpha2Exp')
P%EFTAlpha30 = Ini_Read_Double('EFTAlpha30')
P%EFTAlpha3Exp = Ini_Read_Double('EFTAlpha3Exp')
P%EFTAlpha40 = Ini_Read_Double('EFTAlpha40')
P%EFTAlpha4Exp = Ini_Read_Double('EFTAlpha4Exp')
P%EFTAlpha50 = Ini_Read_Double('EFTAlpha50')
P%EFTAlpha5Exp = Ini_Read_Double('EFTAlpha5Exp')
P%EFTAlpha60 = Ini_Read_Double('EFTAlpha60')
P%EFTAlpha6Exp = Ini_Read_Double('EFTAlpha6Exp')
P%EFTB0 = Ini_Read_Double('EFTB0')
! EFTCAMB MOD END
----------------------------------------------------------------------------------------------------------------------------
Modify:
equations.f90 => equations_EFT.f90
::::: function dtauda(a) :::::
Add: After use LambdaGeneral
! EFTCAMB MOD START
use EFTfunctions
! EFTCAMB MOD END
Replace:
! EFTCAMB MOD START
if (CP%EFTflag==0) then ! Original CAMB code.
if (w_lam == -1._dl) then
grhoa2=grhoa2+grhov*a2**2
else
grhoa2=grhoa2+grhov*a**(1-3*w_lam)
end if
else if (CP%EFTflag/=0) then
grhoa2=grhoa2+grhov*EFTw(a,3)*a2
end if
! EFTCAMB MOD END
::::: module GaugeInterface :::::
Add: After use Errors
! EFTCAMB MOD START
use EFTfunctions
! EFTCAMB MOD END
Add: After real(dl) denlkt(4,max_l_evolve),Kft(max_l_evolve)
! EFTCAMB MOD START
logical EFTCAMBactive
real(dl) EFTturnOnTime
! EFTCAMB MOD END
::::: recursive subroutine GaugeInterface_EvolveScal(EV,tau,y,tauend,tol1,ind,c,w) :::::
Add: After tau_switch_nu_massive= noSwitch
! EFTCAMB MOD START: add the dark energy switch.
if (CP%EFTflag==0) EV%EFTturnOnTime = noSwitch
! EFTCAMB MOD END
Replace:
! EFTCAMB MOD START: add the dark energy switch.
next_switch = min(tau_switch_ktau, tau_switch_nu_massless,EV%TightSwitchoffTime, tau_switch_nu_massive, &
tau_switch_no_nu_multpoles, tau_switch_no_phot_multpoles, tau_switch_nu_nonrel,noSwitch, EV%EFTturnOnTime)
! Original code:
! next_switch = min(tau_switch_ktau, tau_switch_nu_massless,EV%TightSwitchoffTime, tau_switch_nu_massive, &
! tau_switch_no_nu_multpoles, tau_switch_no_phot_multpoles, tau_switch_nu_nonrel,noSwitch)
! EFTCAMB MOD END
Add: Inside else if (next_switch==tau_switch_no_phot_multpoles) then After EV=EVout
! EFTCAMB MOD START: effect of the DE switch.
else if (next_switch==EV%EFTturnOnTime) then
ind=1
EVout%EFTCAMBactive = .true.
EVout%EFTturnOnTime = noSwitch
call SetupScalarArrayIndices(EVout)
call CopyScalarVariableArray(y,yout, EV, EVout)
y=yout
EV=EVout
call EFTpiInitialConditions(y,EV,tau)
! EFTCAMB MOD STOP
::::: subroutine SetupScalarArrayIndices(EV, max_num_eqns) :::::
Replace:
! EFTCAMB MOD START
if ((w_lam /= -1 .and. w_Perturb).or.CP%EFTflag/=0) then
EV%w_ix = neq+1
neq=neq+2
maxeq=maxeq+2
else
EV%w_ix=0
end if
! Original code:
! !Dark energy
! if (w_lam /= -1 .and. w_Perturb) then
! EV%w_ix = neq+1
! neq=neq+2
! maxeq=maxeq+2
! else
! EV%w_ix=0
! end if
! EFTCAMB MOD END
::::: subroutine CopyScalarVariableArray(y,yout, EV, EVout) :::::
Replace:
! EFTCAMB MOD START
if ((w_lam /= -1.and.w_Perturb).or.&
&(CP%EFTflag/=0.and.EVout%EFTCAMBactive).or.&
&(CP%EFTflag/=0.and.EV%EFTCAMBactive)) then
yout(EVout%w_ix)=y(EV%w_ix)
yout(EVout%w_ix+1)=y(EV%w_ix+1)
end if
! Original code:
! if (w_lam /= -1 .and. w_Perturb) then
! yout(EVout%w_ix)=y(EV%w_ix)
! yout(EVout%w_ix+1)=y(EV%w_ix+1)
! end if
! EFTCAMB MOD END
::::: subroutine GetNumEqns(EV) :::::
Add: After EV%MassiveNuApprox=.false.
! EFTCAMB MOD START:
EV%EFTCAMBactive = .false.
! EFTCAMB MOD END
::::: subroutine MassiveNuVars(EV,y,a,grho,gpres,dgrho,dgq, wnu_arr) :::::
Replace:
! EFTCAMB MOD START: compatibility with massive neutrinos
subroutine MassiveNuVars(EV,y,a,grho,gpres,dgrho,dgq, wnu_arr, dgp)
! Original code:
! subroutine MassiveNuVars(EV,y,a,grho,gpres,dgrho,dgq, wnu_arr)
! EFTCAMB MOD END.
implicit none
type(EvolutionVars) EV
real(dl) :: y(EV%nvar), a, grho,gpres,dgrho,dgq
real(dl), intent(out), optional :: wnu_arr(max_nu)
! EFTCAMB MOD START: compatibility with massive neutrinos
real(dl), intent(out), optional :: dgp
! EFTCAMB MOD END.
Add: After dgq = dgq + grhonu_t*qnu
! EFTCAMB MOD START: compatibility with massive neutrinos
if (present(dgp)) then
dgp = dgp + gpnu_t*clxnu
end if
! EFTCAMB MOD END.
-1-