-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathcolh2diss.f
120 lines (120 loc) · 3.27 KB
/
colh2diss.f
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
c=======================================================================
c///////////////////// SUBROUTINE COLH2DISS \\\\\\\\\\\\\\\\\\\\\\\\c
subroutine colh2diss(t, f1, f2, f3, f4, f5, f6, f7)
c
c COMPUTE DENSITY DEPENDENT COLLISIONAL H2 DISSOCIATION RATE
c
c written by: Tom Abel
c date:
c modified1: Feb, 2000 by Greg Bryan; adapted to AMR
c modified2: 2005 by Alexei Razoumov; adapted to FTTE
c
c PURPOSE:
c Computes the 7 temperature dependant functions required to
c generate the density-dependant k13 rate.
c
c compute density dependent collisional H2 dissociation by HI
c data from Martin, Schwartz, Mandy, 1996, ApJ 461, 265
c returns log (base 10) of the rate coefficient in cm^3/s
c of the reaction:
c H2 + H -> 3 H
c Tom Abel 2000
c
c UNITS:
c log10(cgs)
c
c PARAMETERS:
c
c INPUTS:
C T is the gas temperature in Kelvin
c
c OUTPUTS:
c f1-7: rates as given below
c
c-----------------------------------------------------------------------
c
implicit NONE
double precision f1, f2, f3, f4, f5, f6, f7, t
double precision tl, a1, a, b, b1, c, c1, d, ftd
double precision y(21)
c
c Set fitting function values to return very small CID
c (these are returned if T is out of range).
c
f1 = 1.0e-20
f2 = 1.0e-20
f3 = 1.0e-20
f4 = 1.0e-20
f5 = 1.0
f6 = 1.0
f7 = 0.0
c
c do not use 1.0e-20 values for temperatures below 500 K also do not
c return values for temperatures above 1 million Kelvin. Note that
c data and fits are only accurate for t< 5 10^5. However,
c collisional dissociation by electrons will be dominant above this
c temperature, anyhow.
c
if (t .le. 500.) then
c CID = -60.
return
endif
if (t .ge. 1.e6) then
c CID = -60.
return
endif
c
c fitting parameters
c
y(1) = -1.784239D+02
y(2) = -6.842243D+01
y(3) = 4.320243D+01
y(4) = -4.633167D+00
y(5) = 6.970086D+01
y(6) = 4.087038D+04
y(7) = -2.370570D+04
y(8) = 1.288953D+02
y(9) = -5.391334D+01
y(10) = 5.315517D+00
y(11) = -1.973427D+01
y(12) = 1.678095D+04
y(13) = -2.578611D+04
y(14) = 1.482123D+01
y(15) = -4.890915D+00
y(16) = 4.749030D-01
y(17) = -1.338283D+02
y(18) = -1.164408D+00
y(19) = 8.227443D-01
y(20) = 5.864073D-01
y(21) = -2.056313D+00
c
tl=log10(t)
c high density limit
a = y(1)+y(2)*tl+y(3)*tl*tl+y(4)*tl*tl*tl
$ +y(5)*dlog10(1.0+y(6)/t)
a1= y(7)/t
c low density limit
b = y(8)+y(9)*tl+y(10)*tl*tl+y(11)*dlog10(1.0+y(12)/t)
b1= y(13)/t
c critical density
c = y(14)+y(15)*tl+y(16)*tl*tl+y(17)/t
c1 = y(18)+c
d = y(19)+y(20)*exp(-t/1850.0)+y(21)*exp(-t/440.0)
c tabulate the following temperature dependent coefficients:
f1 = a
f2 = (a-b)
f3 = a1
f4 = (a1-b1)
f5 = 10.**c
f6 = 10.**c1
f7 = d
c
c and then get the dissociation rate in cm^3/s with this
c (note: this is in log10)
c
c CID = f1-f2/(1.0+(nh/f5)**f7)
c $ +f3-f4/(1.0+(nh/f6)**f7)
c
c
return
end