-
Notifications
You must be signed in to change notification settings - Fork 1
/
cooling_h.f90
108 lines (84 loc) · 2.98 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
!>
!! \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
implicit none
integer,parameter,private :: temppoints=81 !< 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),private :: mintemp !< lowest log10(temperature) in table
real(kind=dp),private :: dtemp !< steps in log10(temperature) in table
contains
!===========================================================================
!> Calculate the cooling rate
function coolin(nucldens,eldens,xh,temp0)
real(kind=dp) :: coolin
real(kind=dp),intent(in) :: nucldens !< number density
real(kind=dp),intent(in) :: eldens !< electron density
real(kind=dp),dimension(0:1),intent(in) :: xh !< H ionization fractions
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*( &
xh(0)*(h0_cool(itpos)+ &
(h0_cool(itpos1)-h0_cool(itpos))*dtpos)+ &
xh(1)*(h1_cool(itpos)+ &
(h1_cool(itpos1)-h1_cool(itpos))*dtpos))
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-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),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)
! 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)
enddo
end subroutine setup_cool
end module radiative_cooling