Issue with sobroutine inithermo

Use of Cobaya. camb, CLASS, cosmomc, compilers, etc.
Post Reply
raul gonzalez
Posts: 3
Joined: December 05 2013
Affiliation: research

Issue with sobroutine inithermo

Post by raul gonzalez » April 05 2014

Hello everyone
I have work with a fortran program
one of its subroutine is subroutine inithermo that Compute and save unperturbed baryon temperature and ionization fraction as a function of time.

unfortunately i do not understand how this subroutine work.
Are you Can help me ? are you can describe its line & parameters ?

here is the sobroutine:
implicit double precision (a-h,o-z)
c
common /lingerinc/ omegab,omegac,omegav,omegan,h0,tcmb,yhe,
& nnur,nnunr
common /genparm/ grhom,grhog,grhor,adotrad
c
parameter (barssc0=9.1820d-14)
c
parameter (nthermo=10000)
dimension tb(nthermo),cs2(nthermo),xe(nthermo)
dimension dtb(nthermo),dcs2(nthermo),dxe(nthermo)
common /thermod/ tauminn,dlntau,tb,cs2,xe,dtb,dcs2,dxe
save /thermod/
c
thomc0=5.0577d-8*tcmb**4
tauminn=taumin
dlntau=log(taumax/taumin)/(nthermo-1)
c
c Initial conditions: assume radiation-dominated universe.
tau0=taumin
adot0=adotrad
a0=adotrad*taumin
c Assume that any entropy generation occurs before taumin.
c This gives wrong temperature before pair annihilation, but
c the error is harmless.
tb(1)=tcmb/a0
xe0=1.0d0
x1=0.0d0
x2=1.0d0
xe(1)=xe0+0.25d0*yhe/(1.0d0-yhe)*(x1+2*x2)
barssc=barssc0*(1.d0-0.75d0*yhe+(1.d0-yhe)*xe(1))
cs2(1)=4.0d0/3.0d0*barssc*tb(1)
c
do 10 i=2,nthermo
tau=taumin*exp((i-1)*dlntau)
dtau=tau-tau0
c Integrate Friedmann equation using inverse trapezoidal rule.
a=a0+adot0*dtau
a2=a*a
call nu1(a,rhonu,pnu)
grho=grhom*(omegac+omegab)/a+(grhog+grhor*(nnur+nnunr*rhonu))/a2
2 +grhom*omegav*a2
adot=sqrt(grho/3.0d0)*a
a=a0+2.0d0*dtau/(1.0d0/adot0+1.0d0/adot)
c Baryon temperature evolution: adiabatic except for Thomson cooling.
c Use quadrature solution.
tg0=tcmb/a0
ahalf=0.5d0*(a0+a)
adothalf=0.5d0*(adot0+adot)
c fe=number of free electrons divided by total number of free baryon
c particles (e+p+H+He). Evaluate at timstep i-1 for convenience; if
c more accuracy is required (unlikely) then this can be iterated with
c the solution of the ionization equation.
fe=(1.d0-yhe)*xe(i-1)/(1.d0-0.75d0*yhe+(1.d0-yhe)*xe(i-1))
thomc=thomc0*fe/adothalf/ahalf**3
etc=exp(-thomc*(a-a0))
a2t=a0*a0*(tb(i-1)-tg0)*etc-tcmb/thomc*(1.d0-etc)
tb(i)=tcmb/a+a2t/(a*a)
c Integrate ionization equation.
tbhalf=0.5d0*(tb(i-1)+tb(i))
call ionize(tbhalf,ahalf,adothalf,dtau,xe0)
call ionhe(tb(i),a,xe0,x1,x2)
xe(i)=xe0+0.25d0*yhe/(1.0d0-yhe)*(x1+2*x2)
c Baryon sound speed squared (over c**2).
adotoa=adot/a
dtbdla=-2.0d0*tb(i)-thomc*a2t/a
barssc=barssc0*(1.d0-0.75d0*yhe+(1.d0-yhe)*xe(i))
cs2(i)=barssc*tb(i)*(1-dtbdla/tb(i)/3.0d0)
c
a0=a
tau0=tau
adot0=adot
10 continue
c
call splini
call splder(xe,dxe,nthermo)
call splder(cs2,dcs2,nthermo)
call splder(tb,dtb,nthermo)
c
return
end

Post Reply