-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathcooling_h.f90
175 lines (147 loc) · 5.33 KB
/
cooling_h.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
!>
!! \brief This module contains data and subroutines for radiative cooling
!!
!! Module for Capreole / C2-Ray (f90)
!!
!! \b Author: Garrelt Mellema
!!
!! \b Date:
!!
!! \b Version: Non-equilibrium H-only cooling
module radiative_cooling
! This module file contains subroutines for radiative cooling
! - coolin - calculate cooling rate
! - setup_cool - setup cooling table
! Version: H-only cooling
use precision, only: dp
use abundances, only: abu_he
implicit none
integer,parameter,private :: temppoints=801 !< number of points in cooling tables
real(kind=dp),dimension(temppoints),private :: h0_cool !< H0 cooling table
real(kind=dp),dimension(temppoints),private :: h1_cool !< H+ cooling table
real(kind=dp),dimension(temppoints),private :: he0_cool !< He0 cooling table
real(kind=dp),dimension(temppoints),private :: he1_cool !< He1 cooling table
real(kind=dp),dimension(temppoints),private :: he2_cool !< He2 cooling table
real(kind=dp),private :: mintemp !< lowest log10(temperature) in table
real(kind=dp),private :: dtemp !< steps in log10(temperature) in table
real(kind=dp),private :: dumi
contains
!===========================================================================
!> Calculate the cooling rate
function coolin(nucldens,eldens,x_HI, x_HII, x_HeI, x_HeII, x_HeIII, temp0)
real(kind=dp) :: coolin
real(kind=dp),intent(in) :: nucldens !< number density
real(kind=dp),intent(in) :: eldens !< electron density
real(kind=dp),intent(in) :: x_HI
real(kind=dp),intent(in) :: x_HII
real(kind=dp),intent(in) :: x_HeI
real(kind=dp),intent(in) :: x_HeII
real(kind=dp),intent(in) :: x_HeIII
real(kind=dp),intent(in) :: temp0 !< temperature
real(kind=dp) :: tpos, dtpos
integer :: itpos,itpos1
tpos=(log10(temp0)-mintemp)/dtemp+1.0d0
itpos=min(temppoints-1,max(1,int(tpos)))
dtpos=tpos-real(itpos)
itpos1=min(temppoints,itpos+1)
! Cooling curve
coolin=nucldens*eldens*(( &
x_HI*(h0_cool(itpos)+ &
(h0_cool(itpos1)-h0_cool(itpos))*dtpos)+ &
x_HII*(h1_cool(itpos)+ &
(h1_cool(itpos1)-h1_cool(itpos))*dtpos))*(1.0_dp-abu_he)+ &
(x_HeI*(he0_cool(itpos)+ &
(he0_cool(itpos1)-he0_cool(itpos))*dtpos)+&
x_HeII*(he1_cool(itpos)+ &
(he1_cool(itpos1)-he1_cool(itpos))*dtpos)+&
x_HeIII*(he2_cool(itpos)+ &
(he2_cool(itpos1)-he2_cool(itpos))*dtpos))*abu_he)
end function coolin
!===========================================================================
!> Read in and set up the cooling table(s)
subroutine setup_cool ()
real(kind=dp),dimension(temppoints) :: temp
integer :: itemp
integer :: element,ion,nchck
! Open cooling table (H0)
open(unit=22,file='tables/H0-cool.tab',status='old')
! Read the cooling data
read(22,*) element,ion,nchck
if (nchck.eq.0) then
do itemp=1,temppoints
read(22,*) temp(itemp),h0_cool(itemp)
enddo
else
write(*,*) 'Error reading cooling tables'
endif
close(22)
mintemp=temp(1)
dtemp=temp(2)-temp(1)
! not needed: maxtemp=temp(temppoints)
! Open cooling table (H1)
open(unit=22,file='tables/H1-cool-B.tab',status='old')
! Read the cooling data
read(22,*) element,ion,nchck
if (nchck.eq.0) then
do itemp=1,temppoints
read(22,*) temp(itemp),h1_cool(itemp)
enddo
else
write(*,*) 'Error reading cooling tables'
endif
close(22)
! Open cooling table (He0)
!open(unit=22,file='../tables/He0-cool.tab',status='old')
! Since it is not clear what is in the old table, we use now:
! The data from the new table is compiled from:
! collisional ionization from Hui&Gnedin 1997
open(unit=22,file='tables/He0-cool_new.tab',status='old')
! Read the cooling data
read(22,*) element,ion,nchck
! if (nchck.eq.0) then
do itemp=1,temppoints
read(22,*) temp(itemp),he0_cool(itemp)
enddo
! else
! write(*,*) 'Error reading cooling tables'
! endif
close(22)
! Open cooling table (He1)
!open(unit=22,file='../tables/He1-cool.tab',status='old')
! Since it is not clear what is in the old table, we use now:
! The data from the new table is compiled from:
! free-free and recombination B from Hummer&Storey 1998
! collisional excitaion and ionization + dielectronic recomb from Hui&Gnedin 1997
open(unit=22,file='tables/He1-cool_new_nocollion.tab',status='old')
! Read the cooling data
read(22,*) element,ion,nchck
if (nchck.eq.0) then
do itemp=1,temppoints
read(22,*) temp(itemp),he1_cool(itemp)
enddo
else
write(*,*) 'Error reading cooling tables'
endif
close(22)
! Open cooling table (He2)
open(unit=22,file='tables/He2-cool.tab',status='old')
! Read the cooling data
read(22,*) element,ion,nchck
if (nchck.eq.0) then
do itemp=1,temppoints
read(22,*) temp(itemp),he2_cool(itemp)
enddo
else
write(*,*) 'Error reading cooling tables'
endif
close(22)
! Convert cooling to from log to linear
do itemp=1,temppoints
h0_cool(itemp)=10.0d0**h0_cool(itemp)
h1_cool(itemp)=10.0d0**h1_cool(itemp)
he0_cool(itemp)=10.0d0**he0_cool(itemp)
he1_cool(itemp)=10.0d0**he1_cool(itemp)
he2_cool(itemp)=10.0d0**he2_cool(itemp)
enddo
end subroutine setup_cool
end module radiative_cooling