MPI-AMRVAC  2.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_radiative_cooling.t
Go to the documentation of this file.
1 !> module radiative cooling -- add optically thin radiative cooling for HD and MHD
2 !>
3 !> Assumptions: full ionize plasma dominated by H and He, ionization equilibrium
4 !> Formula: Q=-n_H*n_e*f(T), positive f(T) function is pre-computed and tabulated
5 !> Developed by Allard Jan van Marle, Rony Keppens, and Chun Xia
6 !> Cooling tables extend to 1O^9 K, (17.11.2009) AJvM
7 !> Included table by Smith (18.11.2009) AJvM
8 !> Included luminosity output option (18.11.2009) AJvM
9 !> Included two cooling curves from Cloudy code supplied by Wang Ye (12.11.2011) AJvM
10 !> Included a solar coronal cooling curve (JCcorona) supplied by J. Colgan (2008) ApJ
11 !> Modulized and simplified by Chun Xia (2017)
13 !
14 ! these tables contain log_10 temperature values and corresponding
15 ! log_10 luminosity values. The simulation-dependent temperature and luminosity
16 ! scaling parameters are supposed to be provided in the user file.
17 ! All tables have been extended to at least T=10^9 K using a pure Bremmstrahlung
18 ! relationship of Lambda~sqrt(T). This to ensure that a purely explicit calculation
19 ! without timestep check is only used for extremely high temperatures.
20 ! (Except for the SPEX curve, which is more complicated and therefore simply stops
21 ! at the official upper limit of log(T) = 8.16)
22 !
23  use mod_global_parameters, only: std_len
24  use mod_physics
25  implicit none
26  ! parameters used for implicit cooling source calculations
27 
28  !> Coefficent of cooling time step
29  double precision, private :: cfrac
30 
31  !> Lower limit of temperature
32  double precision, private :: tlow
33 
34  !> Helium abundance over Hydrogen
35  double precision, private :: he_abundance
36 
37  !> resolution of temperature in interpolated tables
38  integer, private :: ncool
39 
40  !> Name of cooling curve
41  character(len=std_len), private :: coolcurve
42 
43  !> Name of cooling method
44  character(len=std_len), private :: coolmethod
45 
46  !> Fixed temperature not lower than tlow
47  logical, private :: tfix
48 
49  !> Add cooling source in a split way (.true.) or un-split way (.false.)
50  logical :: rc_split=.false.
51 
52  !> Index of the density (in the w array)
53  integer, private, parameter :: rho_ = 1
54 
55  !> Indices of the momentum density
56  integer, allocatable, private, protected :: mom(:)
57 
58  !> Index of the energy density
59  integer, private, protected :: e_
60 
61  !> The adiabatic index
62  double precision, private :: rc_gamma
63 
64  double precision, allocatable :: tcool(:), lcool(:), dldtcool(:)
65  double precision, allocatable :: yc(:), invyc(:)
66  double precision :: tref, lref, tcoolmin,tcoolmax
67  double precision :: lgtcoolmin, lgtcoolmax, lgstep
68 
69  integer :: n_dm , n_mb , n_mlcosmol &
70  , n_mlwc , n_mlsolar1, n_spex &
72  , n_dm_2
73 
74  double precision :: t_dm(1:71), t_mb(1:51), t_mlcosmol(1:71) &
75  , t_mlwc(1:71), t_mlsolar1(1:71), t_spex(1:110) &
76  , t_jccorona(1:45), t_cl_ism(1:151), t_cl_solar(1:151)&
77  , t_dm_2(1:76)
78 
79  double precision :: l_dm(1:71), l_mb(1:51), l_mlcosmol(1:71) &
80  , l_mlwc(1:71), l_mlsolar1(1:71), l_spex(1:110) &
81  , l_jccorona(1:45), l_cl_ism(1:151), l_cl_solar(1:151)&
82  , l_dm_2(1:76)
83 
84  double precision :: nenh_spex(1:110)
85 
86  data n_jccorona / 45 /
87 
88  data t_jccorona / 4.00000, 4.14230, 4.21995, 4.29761, 4.37528, &
89  4.45294, 4.53061, 4.60827, 4.68593, 4.76359, &
90  4.79705, 4.83049, 4.86394, 4.89739, 4.93084, &
91  4.96428, 4.99773, 5.03117, 5.06461, 5.17574, &
92  5.28684, 5.39796, 5.50907, 5.62018, 5.73129, &
93  5.84240, 5.95351, 6.06461, 6.17574, 6.28684, &
94  6.39796, 6.50907, 6.62018, 6.73129, 6.84240, &
95  6.95351, 7.06461, 7.17574, 7.28684, 7.39796, &
96  7.50907, 7.62018, 7.73129, 7.84240, 7.95351 /
97 
98  data l_jccorona / -200.18883, -100.78630, -30.60384, -22.68481, -21.76445, &
99  -21.67936, -21.54218, -21.37958, -21.25172, -21.17584, &
100  -21.15783, -21.14491, -21.13527, -21.12837, -21.12485, &
101  -21.12439, -21.12642, -21.12802, -21.12548, -21.08965, &
102  -21.08812, -21.19542, -21.34582, -21.34839, -21.31701, &
103  -21.29072, -21.28900, -21.34104, -21.43122, -21.62448, &
104  -21.86694, -22.02897, -22.08051, -22.06057, -22.01973, &
105  -22.00000, -22.05161, -22.22175, -22.41452, -22.52581, &
106  -22.56914, -22.57486, -22.56151, -22.53969, -22.51490 /
107 
108  data n_dm / 71 /
109 
110  data t_dm / 2.0, 2.1, 2.2, 2.3, 2.4, &
111  2.5, 2.6, 2.7, 2.8, 2.9, &
112  3.0, 3.1, 3.2, 3.3, 3.4, &
113  3.5, 3.6, 3.7, 3.8, 3.9, &
114  4.0, 4.1, 4.2, 4.3, 4.4, &
115  4.5, 4.6, 4.7, 4.8, 4.9, &
116  5.0, 5.1, 5.2, 5.3, 5.4, &
117  5.5, 5.6, 5.7, 5.8, 5.9, &
118  6.0, 6.1, 6.2, 6.3, 6.4, &
119  6.5, 6.6, 6.7, 6.8, 6.9, &
120  7.0, 7.1, 7.2, 7.3, 7.4, &
121  7.5, 7.6, 7.7, 7.8, 7.9, &
122  8.0, 8.1, 8.2, 8.3, 8.4, &
123  8.5, 8.6, 8.7, 8.8, 8.9, &
124  9.0 /
125 
126 
127  data l_dm / -26.523, -26.398, -26.301, -26.222, -26.097 &
128  , -26.011, -25.936, -25.866, -25.807, -25.754 &
129  , -25.708, -25.667, -25.630, -25.595, -25.564 &
130  , -25.534, -25.506, -25.479, -25.453, -25.429 &
131  , -25.407, -23.019, -21.762, -21.742, -21.754 &
132  , -21.730, -21.523, -21.455, -21.314, -21.229 &
133  , -21.163, -21.126, -21.092, -21.060, -21.175 &
134  , -21.280, -21.390, -21.547, -21.762, -22.050 &
135  , -22.271, -22.521, -22.646, -22.660, -22.676 &
136  , -22.688, -22.690, -22.662, -22.635, -22.609 &
137  , -22.616, -22.646, -22.697, -22.740, -22.788 &
138  , -22.815, -22.785, -22.754, -22.728, -22.703 &
139  , -22.680, -22.630, -22.580, -22.530, -22.480 &
140  , -22.430, -22.380, -22.330, -22.280, -22.230 &
141  , -22.180 /
142 
143  data n_mb / 51 /
144 
145  data t_mb / 4.0, 4.1, 4.2, 4.3, 4.4, &
146  4.5, 4.6, 4.7, 4.8, 4.9, &
147  5.0, 5.1, 5.2, 5.3, 5.4, &
148  5.5, 5.6, 5.7, 5.8, 5.9, &
149  6.0, 6.1, 6.2, 6.3, 6.4, &
150  6.5, 6.6, 6.7, 6.8, 6.9, &
151  7.0, 7.1, 7.2, 7.3, 7.4, &
152  7.5, 7.6, 7.7, 7.8, 7.9, &
153  8.0, 8.1, 8.2, 8.3, 8.4, &
154  8.5, 8.6, 8.7, 8.8, 8.9, &
155  9.0 /
156 
157  data l_mb / -23.133, -22.895, -22.548, -22.285, -22.099 &
158  , -21.970, -21.918, -21.826, -21.743, -21.638 &
159  , -21.552, -21.447, -21.431, -21.418, -21.461 &
160  , -21.570, -21.743, -21.832, -21.908, -21.981 &
161  , -22.000, -21.998, -21.992, -22.013, -22.095 &
162  , -22.262, -22.397, -22.445, -22.448, -22.446 &
163  , -22.448, -22.465, -22.575, -22.725, -22.749 &
164  , -22.768, -22.753, -22.717, -22.678, -22.637 &
165  , -22.603, -22.553, -22.503, -22.453, -22.403 &
166  , -22.353, -22.303, -22.253, -22.203, -22.153 &
167  , -22.103 /
168 
169  data n_mlcosmol / 71 /
170 
171  data t_mlcosmol / &
172  2.0, 2.1, 2.2, 2.3, 2.4, &
173  2.5, 2.6, 2.7, 2.8, 2.9, &
174  3.0, 3.1, 3.2, 3.3, 3.4, &
175  3.5, 3.6, 3.7, 3.8, 3.9, &
176  4.0, 4.1, 4.2, 4.3, 4.4, &
177  4.5, 4.6, 4.7, 4.8, 4.9, &
178  5.0, 5.1, 5.2, 5.3, 5.4, &
179  5.5, 5.6, 5.7, 5.8, 5.9, &
180  6.0, 6.1, 6.2, 6.3, 6.4, &
181  6.5, 6.6, 6.7, 6.8, 6.9, &
182  7.0, 7.1, 7.2, 7.3, 7.4, &
183  7.5, 7.6, 7.7, 7.8, 7.9, &
184  8.0, 8.1, 8.2, 8.3, 8.4, &
185  8.5, 8.6, 8.7, 8.8, 8.9, &
186  9.0 /
187 
188  data l_mlcosmol / &
189  -99.000, -99.000, -99.000, -99.000, -99.000 &
190  , -99.000, -99.000, -99.000, -99.000, -99.000 &
191  , -99.000, -99.000, -99.000, -99.000, -99.000 &
192  , -99.000, -44.649, -38.362, -33.324, -29.292 &
193  , -26.063, -23.532, -22.192, -22.195, -22.454 &
194  , -22.676, -22.909, -22.925, -22.499, -22.276 &
195  , -22.440, -22.688, -22.917, -23.116, -23.274 &
196  , -23.394, -23.472, -23.516, -23.530, -23.525 &
197  , -23.506, -23.478, -23.444, -23.408, -23.368 &
198  , -23.328, -23.286, -23.244, -23.201, -23.157 &
199  , -23.114, -23.070, -23.026, -22.981, -22.937 &
200  , -22.893, -22.848, -22.803, -22.759, -22.714 &
201  , -22.669, -22.619, -22.569, -22.519, -22.469 &
202  , -22.419, -22.369, -22.319, -22.269, -22.190 &
203  , -22.169 /
204 
205  data n_mlwc / 71 /
206 
207  data t_mlwc / &
208  2.0, 2.1, 2.2, 2.3, 2.4, &
209  2.5, 2.6, 2.7, 2.8, 2.9, &
210  3.0, 3.1, 3.2, 3.3, 3.4, &
211  3.5, 3.6, 3.7, 3.8, 3.9, &
212  4.0, 4.1, 4.2, 4.3, 4.4, &
213  4.5, 4.6, 4.7, 4.8, 4.9, &
214  5.0, 5.1, 5.2, 5.3, 5.4, &
215  5.5, 5.6, 5.7, 5.8, 5.9, &
216  6.0, 6.1, 6.2, 6.3, 6.4, &
217  6.5, 6.6, 6.7, 6.8, 6.9, &
218  7.0, 7.1, 7.2, 7.3, 7.4, &
219  7.5, 7.6, 7.7, 7.8, 7.9, &
220  8.0, 8.1, 8.2, 8.3, 8.4, &
221  8.5, 8.6, 8.7, 8.8, 8.9, &
222  9.0 /
223 
224  data l_mlwc / &
225  -21.192, -21.160, -21.150, -21.150, -21.166 &
226  , -21.191, -21.222, -21.264, -21.308, -21.357 &
227  , -21.408, -21.449, -21.494, -21.544, -21.587 &
228  , -21.638, -21.686, -21.736, -21.780, -21.800 &
229  , -21.744, -21.547, -21.208, -20.849, -20.345 &
230  , -19.771, -19.409, -19.105, -18.827, -18.555 &
231  , -18.460, -18.763, -19.168, -19.334, -19.400 &
232  , -19.701, -20.090, -20.288, -20.337, -20.301 &
233  , -20.233, -20.275, -20.363, -20.508, -20.675 &
234  , -20.856, -21.025, -21.159, -21.256, -21.320 &
235  , -21.354, -21.366, -21.361, -21.343, -21.317 &
236  , -21.285, -21.250, -21.212, -21.172, -21.131 &
237  , -21.089, -21.039, -20.989, -20.939, -20.889 &
238  , -20.839, -20.789, -20.739, -20.689, -20.639 &
239  , -20.589 /
240 
241  data n_mlsolar1 / 71 /
242 
243  data t_mlsolar1 / &
244  2.0, 2.1, 2.2, 2.3, 2.4, &
245  2.5, 2.6, 2.7, 2.8, 2.9, &
246  3.0, 3.1, 3.2, 3.3, 3.4, &
247  3.5, 3.6, 3.7, 3.8, 3.9, &
248  4.0, 4.1, 4.2, 4.3, 4.4, &
249  4.5, 4.6, 4.7, 4.8, 4.9, &
250  5.0, 5.1, 5.2, 5.3, 5.4, &
251  5.5, 5.6, 5.7, 5.8, 5.9, &
252  6.0, 6.1, 6.2, 6.3, 6.4, &
253  6.5, 6.6, 6.7, 6.8, 6.9, &
254  7.0, 7.1, 7.2, 7.3, 7.4, &
255  7.5, 7.6, 7.7, 7.8, 7.9, &
256  8.0, 8.1, 8.2, 8.3, 8.4, &
257  8.5, 8.6, 8.7, 8.8, 8.9, &
258  9.0 /
259 
260  data l_mlsolar1 / &
261  -26.983, -26.951, -26.941, -26.940, -26.956 &
262  , -26.980, -27.011, -27.052, -27.097, -27.145 &
263  , -27.195, -27.235, -27.279, -27.327, -27.368 &
264  , -27.415, -27.456, -27.485, -27.468, -27.223 &
265  , -25.823, -23.501, -22.162, -22.084, -22.157 &
266  , -22.101, -21.974, -21.782, -21.542, -21.335 &
267  , -21.251, -21.275, -21.236, -21.173, -21.167 &
268  , -21.407, -21.670, -21.788, -21.879, -22.008 &
269  , -22.192, -22.912, -22.918, -22.887, -22.929 &
270  , -23.023, -23.094, -23.117, -23.108, -23.083 &
271  , -23.049, -23.011, -22.970, -22.928, -22.885 &
272  , -22.842, -22.798, -22.754, -22.709, -22.665 &
273  , -22.620, -22.570, -22.520, -22.470, -22.420 &
274  , -22.370, -22.320, -22.270, -22.220, -22.170 &
275  , -22.120 /
276 
277 
278  data n_spex / 110 /
279 
280  data t_spex / &
281  3.80, 3.84, 3.88, 3.92, 3.96, &
282  4.00, 4.04, 4.08, 4.12, 4.16, &
283  4.20, 4.24, 4.28, 4.32, 4.36, &
284  4.40, 4.44, 4.48, 4.52, 4.56, &
285  4.60, 4.64, 4.68, 4.72, 4.76, &
286  4.80, 4.84, 4.88, 4.92, 4.96, &
287  5.00, 5.04, 5.08, 5.12, 5.16, &
288  5.20, 5.24, 5.28, 5.32, 5.36, &
289  5.40, 5.44, 5.48, 5.52, 5.56, &
290  5.60, 5.64, 5.68, 5.72, 5.76, &
291  5.80, 5.84, 5.88, 5.92, 5.96, &
292  6.00, 6.04, 6.08, 6.12, 6.16, &
293  6.20, 6.24, 6.28, 6.32, 6.36, &
294  6.40, 6.44, 6.48, 6.52, 6.56, &
295  6.60, 6.64, 6.68, 6.72, 6.76, &
296  6.80, 6.84, 6.88, 6.92, 6.96, &
297  7.00, 7.04, 7.08, 7.12, 7.16, &
298  7.20, 7.24, 7.28, 7.32, 7.36, &
299  7.40, 7.44, 7.48, 7.52, 7.56, &
300  7.60, 7.64, 7.68, 7.72, 7.76, &
301  7.80, 7.84, 7.88, 7.92, 7.96, &
302  8.00, 8.04, 8.08, 8.12, 8.16 /
303 
304  data l_spex / &
305  -25.7331, -25.0383, -24.4059, -23.8288, -23.3027 &
306  , -22.8242, -22.3917, -22.0067, -21.6818, -21.4529 &
307  , -21.3246, -21.3459, -21.4305, -21.5293, -21.6138 &
308  , -21.6615, -21.6551, -21.5919, -21.5092, -21.4124 &
309  , -21.3085, -21.2047, -21.1067, -21.0194, -20.9413 &
310  , -20.8735, -20.8205, -20.7805, -20.7547, -20.7455 &
311  , -20.7565, -20.7820, -20.8008, -20.7994, -20.7847 &
312  , -20.7687, -20.7590, -20.7544, -20.7505, -20.7545 &
313  , -20.7888, -20.8832, -21.0450, -21.2286, -21.3737 &
314  , -21.4573, -21.4935, -21.5098, -21.5345, -21.5863 &
315  , -21.6548, -21.7108, -21.7424, -21.7576, -21.7696 &
316  , -21.7883, -21.8115, -21.8303, -21.8419, -21.8514 &
317  , -21.8690, -21.9057, -21.9690, -22.0554, -22.1488 &
318  , -22.2355, -22.3084, -22.3641, -22.4033, -22.4282 &
319  , -22.4408, -22.4443, -22.4411, -22.4334, -22.4242 &
320  , -22.4164, -22.4134, -22.4168, -22.4267, -22.4418 &
321  , -22.4603, -22.4830, -22.5112, -22.5449, -22.5819 &
322  , -22.6177, -22.6483, -22.6719, -22.6883, -22.6985 &
323  , -22.7032, -22.7037, -22.7008, -22.6950, -22.6869 &
324  , -22.6769, -22.6655, -22.6531, -22.6397, -22.6258 &
325  , -22.6111, -22.5964, -22.5816, -22.5668, -22.5519 &
326  , -22.5367, -22.5216, -22.5062, -22.4912, -22.4753 /
327 
328 
329  data nenh_spex / &
330  0.000013264, 0.000042428, 0.000088276, 0.00017967 &
331  , 0.00084362, 0.0034295, 0.013283, 0.042008 &
332  , 0.12138, 0.30481, 0.53386, 0.76622 &
333  , 0.89459, 0.95414, 0.98342 &
334  , 1.0046, 1.0291, 1.0547 &
335  , 1.0767, 1.0888 &
336  , 1.0945, 1.0972, 1.0988 &
337  , 1.1004, 1.1034 &
338  , 1.1102, 1.1233, 1.1433 &
339  , 1.1638, 1.1791 &
340  , 1.1885, 1.1937, 1.1966 &
341  , 1.1983, 1.1993 &
342  , 1.1999, 1.2004, 1.2008 &
343  , 1.2012, 1.2015 &
344  , 1.2020, 1.2025, 1.2030 &
345  , 1.2035, 1.2037 &
346  , 1.2039, 1.2040, 1.2041 &
347  , 1.2042, 1.2044 &
348  , 1.2045, 1.2046, 1.2047 &
349  , 1.2049, 1.2050 &
350  , 1.2051, 1.2053, 1.2055 &
351  , 1.2056, 1.2058 &
352  , 1.2060, 1.2062, 1.2065 &
353  , 1.2067, 1.2070 &
354  , 1.2072, 1.2075, 1.2077 &
355  , 1.2078, 1.2079 &
356  , 1.2080, 1.2081, 1.2082 &
357  , 1.2083, 1.2083 &
358  , 1.2084, 1.2084, 1.2085 &
359  , 1.2085, 1.2086 &
360  , 1.2086, 1.2087, 1.2087 &
361  , 1.2088, 1.2088 &
362  , 1.2089, 1.2089, 1.2089 &
363  , 1.2089, 1.2089 &
364  , 1.2090, 1.2090, 1.2090 &
365  , 1.2090, 1.2090 &
366  , 1.2090, 1.2090, 1.2090 &
367  , 1.2090, 1.2090 &
368  , 1.2090, 1.2090, 1.2090 &
369  , 1.2090, 1.2090 &
370  , 1.2090, 1.2090, 1.2090 &
371  , 1.2090, 1.2090 /
372  !
373  ! To be used together with the SPEX table for the SPEX_DM option
374  ! Assuming an ionization fraction of 10^-3
375  !
376  data n_dm_2 / 76 /
377 
378  data t_dm_2 / 1.00, 1.04, 1.08, 1.12, 1.16, 1.20 &
379  , 1.24, 1.28, 1.32, 1.36, 1.40 &
380  , 1.44, 1.48, 1.52, 1.56, 1.60 &
381  , 1.64, 1.68, 1.72, 1.76, 1.80 &
382  , 1.84, 1.88, 1.92, 1.96, 2.00 &
383  , 2.04, 2.08, 2.12, 2.16, 2.20 &
384  , 2.24, 2.28, 2.32, 2.36, 2.40 &
385  , 2.44, 2.48, 2.52, 2.56, 2.60 &
386  , 2.64, 2.68, 2.72, 2.76, 2.80 &
387  , 2.84, 2.88, 2.92, 2.96, 3.00 &
388  , 3.04, 3.08, 3.12, 3.16, 3.20 &
389  , 3.24, 3.28, 3.32, 3.36, 3.40 &
390  , 3.44, 3.48, 3.52, 3.56, 3.60 &
391  , 3.64, 3.68, 3.72, 3.76, 3.80 &
392  , 3.84, 3.88, 3.92, 3.96, 4.00 /
393 
394 
395  data l_dm_2 / -30.0377, -29.7062, -29.4055, -29.1331, -28.8864, -28.6631 &
396  , -28.4614, -28.2791, -28.1146, -27.9662, -27.8330 &
397  , -27.7129, -27.6052, -27.5088, -27.4225, -27.3454 &
398  , -27.2767, -27.2153, -27.1605, -27.1111, -27.0664 &
399  , -27.0251, -26.9863, -26.9488, -26.9119, -26.8742 &
400  , -26.8353, -26.7948, -26.7523, -26.7080, -26.6619 &
401  , -26.6146, -26.5666, -26.5183, -26.4702, -26.4229 &
402  , -26.3765, -26.3317, -26.2886, -26.2473, -26.2078 &
403  , -26.1704, -26.1348, -26.1012, -26.0692, -26.0389 &
404  , -26.0101, -25.9825, -25.9566, -25.9318, -25.9083 &
405  , -25.8857, -25.8645, -25.8447, -25.8259, -25.8085 &
406  , -25.7926, -25.7778, -25.7642, -25.7520, -25.7409 &
407  , -25.7310, -25.7222, -25.7142, -25.7071, -25.7005 &
408  , -25.6942, -25.6878, -25.6811, -25.6733, -25.6641 &
409  , -25.6525, -25.6325, -25.6080, -25.5367, -25.4806 /
410 
411  data n_cl_solar / 151 /
412 
413  data t_cl_solar / &
414  1.000 , 1.050 , 1.100 , 1.150 , &
415  1.200 , 1.250 , 1.300 , 1.350 , &
416  1.400 , 1.450 , 1.500 , 1.550 , &
417  1.600 , 1.650 , 1.700 , 1.750 , &
418  1.800 , 1.850 , 1.900 , 1.950 , &
419  2.000 , 2.050 , 2.100 , 2.150 , &
420  2.200 , 2.250 , 2.300 , 2.350 , &
421  2.400 , 2.450 , 2.500 , 2.550 , &
422  2.600 , 2.650 , 2.700 , 2.750 , &
423  2.800 , 2.850 , 2.900 , 2.950 , &
424  3.000 , 3.050 , 3.100 , 3.150 , &
425  3.200 , 3.250 , 3.300 , 3.350 , &
426  3.400 , 3.450 , 3.500 , 3.550 , &
427  3.600 , 3.650 , 3.700 , 3.750 , &
428  3.800 , 3.850 , 3.900 , 3.950 , &
429  4.000 , 4.050 , 4.100 , 4.150 , &
430  4.200 , 4.250 , 4.300 , 4.350 , &
431  4.400 , 4.450 , 4.500 , 4.550 , &
432  4.600 , 4.650 , 4.700 , 4.750 , &
433  4.800 , 4.850 , 4.900 , 4.950 , &
434  5.000 , 5.050 , 5.100 , 5.150 , &
435  5.200 , 5.250 , 5.300 , 5.350 , &
436  5.400 , 5.450 , 5.500 , 5.550 , &
437  5.600 , 5.650 , 5.700 , 5.750 , &
438  5.800 , 5.850 , 5.900 , 5.950 , &
439  6.000 , 6.050 , 6.100 , 6.150 , &
440  6.200 , 6.250 , 6.300 , 6.350 , &
441  6.400 , 6.450 , 6.500 , 6.550 , &
442  6.600 , 6.650 , 6.700 , 6.750 , &
443  6.800 , 6.850 , 6.900 , 6.950 , &
444  7.000 , 7.050 , 7.100 , 7.150 , &
445  7.200 , 7.250 , 7.300 , 7.350 , &
446  7.400 , 7.450 , 7.500 , 7.550 , &
447  7.600 , 7.650 , 7.700 , 7.750 , &
448  7.800 , 7.850 , 7.900 , 7.950 , &
449  8.000 , 8.100 , 8.200 , 8.300 , &
450  8.400 , 8.500 , 8.600 , 8.700 , &
451  8.800 , 8.900 , 9.00 /
452 
453 
454  data l_cl_solar / &
455  -28.375 , -28.251 , -28.137 , -28.029 , &
456  -27.929 , -27.834 , -27.745 , -27.662 , &
457  -27.584 , -27.512 , -27.445 , -27.383 , &
458  -27.326 , -27.273 , -27.223 , -27.175 , &
459  -27.128 , -27.079 , -27.027 , -26.972 , &
460  -26.911 , -26.846 , -26.777 , -26.705 , &
461  -26.632 , -26.554 , -26.479 , -26.407 , &
462  -26.338 , -26.274 , -26.213 , -26.156 , &
463  -26.101 , -26.049 , -25.999 , -25.949 , &
464  -25.901 , -25.852 , -25.803 , -25.754 , &
465  -25.707 , -25.662 , -25.621 , -25.588 , &
466  -25.561 , -25.538 , -25.518 , -25.497 , &
467  -25.475 , -25.452 , -25.426 , -25.400 , &
468  -25.374 , -25.333 , -25.295 , -25.261 , &
469  -25.228 , -25.189 , -25.136 , -25.053 , &
470  -24.888 , -24.454 , -23.480 , -22.562 , &
471  -22.009 , -21.826 , -21.840 , -21.905 , &
472  -21.956 , -21.971 , -21.958 , -21.928 , &
473  -21.879 , -21.810 , -21.724 , -21.623 , &
474  -21.512 , -21.404 , -21.321 , -21.273 , &
475  -21.250 , -21.253 , -21.275 , -21.287 , &
476  -21.282 , -21.275 , -21.272 , -21.267 , &
477  -21.281 , -21.357 , -21.496 , -21.616 , &
478  -21.677 , -21.698 , -21.708 , -21.730 , &
479  -21.767 , -21.793 , -21.794 , -21.787 , &
480  -21.787 , -21.802 , -21.826 , -21.859 , &
481  -21.911 , -21.987 , -22.082 , -22.173 , &
482  -22.253 , -22.325 , -22.392 , -22.448 , &
483  -22.487 , -22.512 , -22.524 , -22.528 , &
484  -22.524 , -22.516 , -22.507 , -22.501 , &
485  -22.502 , -22.511 , -22.533 , -22.565 , &
486  -22.600 , -22.630 , -22.648 , -22.656 , &
487  -22.658 , -22.654 , -22.647 , -22.634 , &
488  -22.619 , -22.602 , -22.585 , -22.566 , &
489  -22.546 , -22.525 , -22.505 , -22.480 , &
490  -22.465 , -22.415 , -22.365 , -22.315 , &
491  -22.265 , -22.215 , -22.165 , -22.115 , &
492  -22.065 , -22.015 , -21.965 /
493 
494 
495  data n_cl_ism / 151 /
496 
497  data t_cl_ism / &
498  1.000 , 1.050 , 1.100 , 1.150 , &
499  1.200 , 1.250 , 1.300 , 1.350 , &
500  1.400 , 1.450 , 1.500 , 1.550 , &
501  1.600 , 1.650 , 1.700 , 1.750 , &
502  1.800 , 1.850 , 1.900 , 1.950 , &
503  2.000 , 2.050 , 2.100 , 2.150 , &
504  2.200 , 2.250 , 2.300 , 2.350 , &
505  2.400 , 2.450 , 2.500 , 2.550 , &
506  2.600 , 2.650 , 2.700 , 2.750 , &
507  2.800 , 2.850 , 2.900 , 2.950 , &
508  3.000 , 3.050 , 3.100 , 3.150 , &
509  3.200 , 3.250 , 3.300 , 3.350 , &
510  3.400 , 3.450 , 3.500 , 3.550 , &
511  3.600 , 3.650 , 3.700 , 3.750 , &
512  3.800 , 3.850 , 3.900 , 3.950 , &
513  4.000 , 4.050 , 4.100 , 4.150 , &
514  4.200 , 4.250 , 4.300 , 4.350 , &
515  4.400 , 4.450 , 4.500 , 4.550 , &
516  4.600 , 4.650 , 4.700 , 4.750 , &
517  4.800 , 4.850 , 4.900 , 4.950 , &
518  5.000 , 5.050 , 5.100 , 5.150 , &
519  5.200 , 5.250 , 5.300 , 5.350 , &
520  5.400 , 5.450 , 5.500 , 5.550 , &
521  5.600 , 5.650 , 5.700 , 5.750 , &
522  5.800 , 5.850 , 5.900 , 5.950 , &
523  6.000 , 6.050 , 6.100 , 6.150 , &
524  6.200 , 6.250 , 6.300 , 6.350 , &
525  6.400 , 6.450 , 6.500 , 6.550 , &
526  6.600 , 6.650 , 6.700 , 6.750 , &
527  6.800 , 6.850 , 6.900 , 6.950 , &
528  7.000 , 7.050 , 7.100 , 7.150 , &
529  7.200 , 7.250 , 7.300 , 7.350 , &
530  7.400 , 7.450 , 7.500 , 7.550 , &
531  7.600 , 7.650 , 7.700 , 7.750 , &
532  7.800 , 7.850 , 7.900 , 7.950 , &
533  8.000 , 8.100 , 8.200 , 8.300 , &
534  8.400 , 8.500 , 8.600 , 8.700 , &
535  8.800 , 8.900 , 9.00 /
536 
537  data l_cl_ism / &
538  -28.365 , -28.242 , -28.127 , -28.020 , &
539  -27.919 , -27.825 , -27.736 , -27.653 , &
540  -27.575 , -27.504 , -27.437 , -27.376 , &
541  -27.319 , -27.267 , -27.220 , -27.176 , &
542  -27.134 , -27.095 , -27.058 , -27.021 , &
543  -26.985 , -26.948 , -26.910 , -26.870 , &
544  -26.827 , -26.775 , -26.721 , -26.664 , &
545  -26.608 , -26.552 , -26.495 , -26.437 , &
546  -26.378 , -26.317 , -26.255 , -26.190 , &
547  -26.123 , -26.053 , -25.984 , -25.913 , &
548  -25.847 , -25.786 , -25.736 , -25.702 , &
549  -25.678 , -25.662 , -25.649 , -25.636 , &
550  -25.621 , -25.604 , -25.587 , -25.571 , &
551  -25.562 , -25.526 , -25.505 , -25.499 , &
552  -25.499 , -25.491 , -25.468 , -25.410 , &
553  -25.268 , -24.888 , -23.702 , -22.624 , &
554  -22.036 , -21.843 , -21.854 , -21.924 , &
555  -21.986 , -22.017 , -22.021 , -22.005 , &
556  -21.964 , -21.896 , -21.806 , -21.699 , &
557  -21.580 , -21.463 , -21.370 , -21.312 , &
558  -21.284 , -21.290 , -21.322 , -21.345 , &
559  -21.354 , -21.366 , -21.385 , -21.396 , &
560  -21.414 , -21.483 , -21.600 , -21.696 , &
561  -21.742 , -21.759 , -21.776 , -21.816 , &
562  -21.885 , -21.939 , -21.946 , -21.918 , &
563  -21.873 , -21.818 , -21.756 , -21.689 , &
564  -21.618 , -21.547 , -21.475 , -21.403 , &
565  -21.331 , -21.260 , -21.188 , -21.114 , &
566  -21.039 , -20.963 , -20.887 , -20.810 , &
567  -20.734 , -20.657 , -20.581 , -20.505 , &
568  -20.429 , -20.352 , -20.276 , -20.200 , &
569  -20.125 , -20.049 , -19.973 , -19.898 , &
570  -19.822 , -19.747 , -19.671 , -19.596 , &
571  -19.520 , -19.445 , -19.370 , -19.295 , &
572  -19.220 , -19.144 , -19.069 , -18.994 , &
573  -18.919 , -18.869 , -18.819 , -18.769 , &
574  -18.719 , -18.669 , -18.619 , -18.569 , &
575  -18.519 , -18.469 , -18.419 /
576 
577  contains
578 
579  !> Read this module"s parameters from a file
580  subroutine rc_params_read(files)
582  character(len=*), intent(in) :: files(:)
583  integer :: n
584 
585  namelist /rc_list/ coolcurve, coolmethod, ncool, cfrac, tlow, tfix
586 
587  do n = 1, size(files)
588  open(unitpar, file=trim(files(n)), status="old")
589  read(unitpar, rc_list, end=111)
590 111 close(unitpar)
591  end do
592 
593  end subroutine rc_params_read
594 
595  !> Radiative cooling initialization
596  subroutine radiative_cooling_init(phys_gamma,He_abund)
598 
599  double precision, intent(in) :: phys_gamma,He_abund
600 
601  double precision, dimension(:), allocatable :: t_table
602  double precision, dimension(:), allocatable :: L_table
603  double precision :: ratt, Lerror
604  double precision ::fact1, fact2, fact3, dL1, dL2
605  double precision :: tstep, Lstep
606  integer :: ntable, i, j, ic, nwx,idir
607  logical :: jump
608 
609  rc_gamma=phys_gamma
610  he_abundance=he_abund
611  ncool=4000
612  coolcurve='JCcorona'
613  coolmethod='exact'
614  cfrac=0.1d0
615  tlow=bigdouble
616  tfix=.false.
617 
619 
620  ! Determine flux variables
621  nwx = 1 ! rho (density)
622 
623  allocate(mom(ndir))
624  do idir = 1, ndir
625  nwx = nwx + 1
626  mom(idir) = nwx ! momentum density
627  end do
628 
629  nwx = nwx + 1
630  e_ = nwx ! energy density
631 
632  allocate(tcool(1:ncool), lcool(1:ncool), dldtcool(1:ncool))
633  allocate(yc(1:ncool), invyc(1:ncool))
634 
635  tcool(1:ncool) = zero
636  lcool(1:ncool) = zero
637  dldtcool(1:ncool) = zero
638 
639  ! Read in the selected cooling curve
640  select case(coolcurve)
641 
642  case('JCcorona')
643  if(mype ==0) &
644  print *,'Use Colgan & Feldman (2008) cooling curve'
645  if(mype ==0) &
646  print *,'This version only till 10000 K, beware for floor T treatment'
647 
648  ntable = n_jccorona
649 
650  allocate(t_table(1:ntable))
651  allocate(l_table(1:ntable))
652  t_table(1:ntable) = t_jccorona(1:n_jccorona)
653  l_table(1:ntable) = l_jccorona(1:n_jccorona)
654 
655  case('DM')
656  if(mype ==0) &
657  print *,'Use Delgano & McCray (1972) cooling curve'
658 
659  ntable = n_dm
660 
661  allocate(t_table(1:ntable))
662  allocate(l_table(1:ntable))
663  t_table(1:ntable) = t_dm(1:n_dm)
664  l_table(1:ntable) = l_dm(1:n_dm)
665 
666  case('MB')
667  if(mype ==0) &
668  write(*,'(3a)') 'Use MacDonald & Bailey (1981) cooling curve '&
669  ,'as implemented in ZEUS-3D, with the values '&
670  ,'from Delgano & McCRay (1972) for low temperatures.'
671 
672  ntable = n_mb + 20
673 
674  allocate(t_table(1:ntable))
675  allocate(l_table(1:ntable))
676 
677  t_table(1:ntable) = t_dm(1:21)
678  l_table(1:ntable) = l_dm(1:21)
679  t_table(22:ntable) = t_mb(2:n_mb)
680  l_table(22:ntable) = l_mb(2:n_mb)
681 
682  case('MLcosmol')
683  if(mype ==0) &
684  print *,'Use Mellema & Lundqvist (2002) cooling curve '&
685  ,'for zero metallicity '
686 
687  ntable = n_mlcosmol
688 
689  allocate(t_table(1:ntable))
690  allocate(l_table(1:ntable))
691  t_table(1:ntable) = t_mlcosmol(1:n_mlcosmol)
692  l_table(1:ntable) = l_mlcosmol(1:n_mlcosmol)
693 
694  case('MLwc')
695  if(mype ==0) &
696  print *,'Use Mellema & Lundqvist (2002) cooling curve '&
697  ,'for WC-star metallicity '
698 
699  ntable = n_mlwc
700 
701  allocate(t_table(1:ntable))
702  allocate(l_table(1:ntable))
703  t_table(1:ntable) = t_mlwc(1:n_mlwc)
704  l_table(1:ntable) = l_mlwc(1:n_mlwc)
705 
706  case('MLsolar1')
707  if(mype ==0) &
708  print *,'Use Mellema & Lundqvist (2002) cooling curve '&
709  ,'for solar metallicity '
710 
711  ntable = n_mlsolar1
712 
713  allocate(t_table(1:ntable))
714  allocate(l_table(1:ntable))
715  t_table(1:ntable) = t_mlsolar1(1:n_mlsolar1)
716  l_table(1:ntable) = l_mlsolar1(1:n_mlsolar1)
717 
718  case('cloudy_ism')
719  if(mype ==0) &
720  print *,'Use Cloudy based cooling curve '&
721  ,'for ism metallicity '
722 
723  ntable = n_cl_ism
724 
725  allocate(t_table(1:ntable))
726  allocate(l_table(1:ntable))
727  t_table(1:ntable) = t_cl_ism(1:n_cl_ism)
728  l_table(1:ntable) = l_cl_ism(1:n_cl_ism)
729 
730  case('cloudy_solar')
731  if(mype ==0) &
732  print *,'Use Cloudy based cooling curve '&
733  ,'for solar metallicity '
734 
735  ntable = n_cl_solar
736 
737  allocate(t_table(1:ntable))
738  allocate(l_table(1:ntable))
739  t_table(1:ntable) = t_cl_solar(1:n_cl_solar)
740  l_table(1:ntable) = l_cl_solar(1:n_cl_solar)
741 
742  case('SPEX')
743  if(mype ==0) &
744  print *,'Use SPEX cooling curve (Schure et al. 2009) '&
745  ,'for solar metallicity '
746 
747  ntable = n_spex
748 
749  allocate(t_table(1:ntable))
750  allocate(l_table(1:ntable))
751  t_table(1:ntable) = t_spex(1:n_spex)
752  l_table(1:ntable) = l_spex(1:n_spex) + log10(nenh_spex(1:n_spex))
753 
754  case('SPEX_DM')
755  if(mype ==0) then
756  print *, 'Use SPEX cooling curve for solar metallicity above 10^4 K. '
757  print *, 'At lower temperatures,use Dalgarno & McCray (1972), '
758  print *, 'with a pre-set ionization fraction of 10^-3. '
759  print *, 'as described by Schure et al. (2009). '
760  endif
761 
762  ntable = n_spex + n_dm_2 - 6
763 
764  allocate(t_table(1:ntable))
765  allocate(l_table(1:ntable))
766  t_table(1:n_dm_2-1) = t_dm_2(1:n_dm_2-1)
767  l_table(1:n_dm_2-1) = l_dm_2(1:n_dm_2-1)
768  t_table(n_dm_2:ntable) = t_spex(6:n_spex)
769  l_table(n_dm_2:ntable) = l_spex(6:n_spex) + log10(nenh_spex(6:n_spex))
770  case default
771  call mpistop("This coolingcurve is unknown")
772  end select
773 
774  ! create cooling table(s) for use in amrvac
775  tcoolmax = t_table(ntable)
776  tcoolmin = t_table(1)
777  ratt = (tcoolmax-tcoolmin)/( dble(ncool-1) + smalldouble)
778 
779  tcool(1) = tcoolmin
780  lcool(1) = l_table(1)
781 
782  tcool(ncool) = tcoolmax
783  lcool(ncool) = l_table(ntable)
784 
785  do i=2,ncool ! loop to create one table
786  tcool(i) = tcool(i-1)+ratt
787  do j=1,ntable-1 ! loop to create one spot on a table
788  ! Second order polynomial interpolation, except at the outer edge,
789  ! or in case of a large jump.
790  if(tcool(i) < t_table(j+1)) then
791  if(j.eq. ntable-1 )then
792  fact1 = (tcool(i)-t_table(j+1)) &
793  /(t_table(j)-t_table(j+1))
794 
795  fact2 = (tcool(i)-t_table(j)) &
796  /(t_table(j+1)-t_table(j))
797 
798  lcool(i) = l_table(j)*fact1 + l_table(j+1)*fact2
799  exit
800  else
801  dl1 = l_table(j+1)-l_table(j)
802  dl2 = l_table(j+2)-l_table(j+1)
803  jump =(max(dabs(dl1),dabs(dl2)) > 2*min(dabs(dl1),dabs(dl2)))
804  endif
805 
806  if( jump ) then
807  fact1 = (tcool(i)-t_table(j+1)) &
808  /(t_table(j)-t_table(j+1))
809 
810  fact2 = (tcool(i)-t_table(j)) &
811  /(t_table(j+1)-t_table(j))
812 
813  lcool(i) = l_table(j)*fact1 + l_table(j+1)*fact2
814  exit
815  else
816  fact1 = ((tcool(i)-t_table(j+1)) &
817  * (tcool(i)-t_table(j+2))) &
818  / ((t_table(j)-t_table(j+1)) &
819  * (t_table(j)-t_table(j+2)))
820 
821  fact2 = ((tcool(i)-t_table(j)) &
822  * (tcool(i)-t_table(j+2))) &
823  / ((t_table(j+1)-t_table(j)) &
824  * (t_table(j+1)-t_table(j+2)))
825 
826  fact3 = ((tcool(i)-t_table(j)) &
827  * (tcool(i)-t_table(j+1))) &
828  / ((t_table(j+2)-t_table(j)) &
829  * (t_table(j+2)-t_table(j+1)))
830 
831  lcool(i) = l_table(j)*fact1 + l_table(j+1)*fact2 &
832  + l_table(j+2)*fact3
833  exit
834  endif
835  endif
836  enddo ! end loop to find create one spot on a table
837  enddo ! end loop to create one table
838 
839  ! Go from logarithmic to actual values.
840  tcool(1:ncool) = 10.0d0**tcool(1:ncool)
841  lcool(1:ncool) = 10.0d0**lcool(1:ncool)
842 
843  ! Scale both T and Lambda
844  tcool(1:ncool) = tcool(1:ncool) / unit_temperature
845  lcool(1:ncool) = lcool(1:ncool) * unit_numberdensity**2 * unit_time / unit_pressure * (1.d0+2.d0*he_abundance)
846  tcoolmin = tcool(1)+smalldouble ! avoid pointless interpolation
847  ! smaller value for lowest temperatures from cooling table and user's choice
848  if (tlow==bigdouble) tlow=tcoolmin
849  tcoolmax = tcool(ncool)
850 
851  lgtcoolmin = dlog10(tcoolmin)
852  lgtcoolmax = dlog10(tcoolmax)
853  lgstep = (lgtcoolmax-lgtcoolmin) * 1.d0 / (ncool-1)
854 
855  dldtcool(1) = (lcool(2)-lcool(1))/(tcool(2)-tcool(1))
856  dldtcool(ncool) = (lcool(ncool)-lcool(ncool-1))/(tcool(ncool)-tcool(ncool-1))
857 
858  do i=2,ncool-1
859  dldtcool(i) = (lcool(i+1)-lcool(i-1))/(tcool(i+1)-tcool(i-1))
860  enddo
861 
862  deallocate(t_table)
863  deallocate(l_table)
864 
865  if( coolmethod == 'exact' ) then
866  tref = tcoolmax
867  lref = lcool(ncool)
868  yc(ncool) = zero
869  do i=ncool-1, 1, -1
870  yc(i) = yc(i+1)
871  do j=1,100
872  tstep = 1.0d-2*(tcool(i+1)-tcool(i))
873  call findl(tcool(i+1)-j*tstep, lstep)
874  yc(i) = yc(i) + lref/tref*tstep/lstep
875  enddo
876  enddo
877  endif
878 
879  end subroutine radiative_cooling_init
880 
881  subroutine cooling_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
883 
884  integer, intent(in) :: ixI^L, ixO^L
885  double precision, intent(in) :: dx^D, x(ixi^s,1:ndim), w(ixi^s,1:nw)
886  double precision, intent(inout) :: dtnew
887 
888  double precision :: etherm(ixi^s)
889  double precision :: L1,Tlocal1, ptherm(ixi^s), lum(ixi^s)
890  double precision :: plocal, rholocal
891  double precision :: Lmax
892  integer :: ix^D
893  !
894  ! Limit timestep to avoid cooling problems when using explicit cooling
895  !
896 
897  if(coolmethod == 'explicit1') then
898  call phys_get_pthermal(w,x,ixi^l,ixo^l,ptherm)
899  {do ix^db = ixo^lim^db\}
900  plocal = ptherm(ix^d)
901  rholocal = w(ix^d,rho_)
902  ! Tlocal = P/rho
903  tlocal1 = max(plocal/(rholocal),smalldouble)
904  ! Determine explicit cooling
905  ! If temperature is below floor level, no cooling.
906  ! Stop wasting time and go to next gridpoint.
907  ! If the temperature is higher than the maximum,
908  ! assume Bremmstrahlung
909  if( tlocal1<=tcoolmin ) then
910  l1 = zero
911  else if( tlocal1>=tcoolmax )then
912  l1 = lcool(ncool)*sqrt(tlocal1/tcoolmax)
913  l1 = l1*(rholocal**2)
914  else
915  call findl(tlocal1,l1)
916  l1 = l1*(rholocal**2)
917  endif
918  lum(ix^d) = l1
919  {enddo^d&\}
920  etherm(ixo^s)=ptherm(ixo^s)/(rc_gamma-1.d0)
921  dtnew = cfrac*minval(etherm(ixo^s)/max(lum(ixo^s),smalldouble))
922  endif
923 
924  end subroutine cooling_get_dt
925 
926  subroutine getvar_cooling(ixI^L,ixO^L,w,x,coolrate)
927  !
928  ! Create extra variable to show cooling rate in the output
929  ! Uses a simple explicit scheme.
930  ! N.B. Since there is no knowledge of the timestep size,
931  ! there is no upper limit for the cooling rate.
933 
934  integer, intent(in) :: ixI^L,ixO^L
935  double precision, intent(in) :: x(ixi^s,1:ndim)
936  double precision :: w(ixi^s,1:nw)
937  double precision, intent(out):: coolrate(ixi^s)
938 
939  double precision :: etherm(ixi^s)
940  double precision :: L1,Tlocal1, ptherm(ixi^s)
941  double precision :: plocal, rholocal
942  double precision :: emin
943 
944  integer :: ix^D
945 
946  call phys_get_pthermal(w,x,ixi^l,ixo^l,ptherm)
947 
948  {do ix^db = ixo^lim^db\}
949  plocal = ptherm(ix^d)
950  rholocal = w(ix^d,rho_)
951  ! Tlocal = P/rho
952  tlocal1 = max(plocal/(rholocal),smalldouble)
953  !
954  ! Determine explicit cooling
955  !
956  if( tlocal1<=tcoolmin ) then
957  l1 = zero
958  else if( tlocal1>=tcoolmax )then
959  l1 = lcool(ncool)*sqrt(tlocal1/tcoolmax)
960  l1 = l1*(rholocal**2)
961  else
962  call findl(tlocal1,l1)
963  l1 = l1*(rholocal**2)
964  endif
965  coolrate(ix^d) = l1
966  {enddo^d&\}
967 
968  end subroutine getvar_cooling
969 
970  subroutine getvar_cooling_exact(qdt, ixI^L, ixO^L, wCT, w, x, coolrate)
971  !
972  ! Calculates cooling rate using the exact cooling method,
973  ! for usage in eg. source_terms subroutine.
974  ! The TEF must be known, so this routine can only be used
975  ! together with the "exact" cooling method.
977 
978  integer, intent(in) :: ixI^L, ixO^L
979  double precision, intent(in) :: qdt, x(ixi^s, 1:ndim), wCT(ixi^s, 1:nw)
980  double precision :: w(ixi^s, 1:nw)
981  double precision, intent(out) :: coolrate(ixi^s)
982 
983  double precision :: y1, y2, l1, l2, tc
984  double precision :: plocal, rholocal, tlocal1, tlocal2, invgam
985  double precision :: ptherm(ixi^s), pnew(ixi^s)
986  double precision :: emin, Lmax, fact
987 
988  integer :: ix^D
989 
990  ! Check cooling method
991  if( coolmethod /= 'exact') then
992  call mpistop("Subroutine getvar_cooling_exact needs the exact cooling method")
993  endif
994 
995  call phys_get_pthermal(wct, x, ixi^l, ixo^l, ptherm)
996  call phys_get_pthermal(w, x, ixi^l, ixo^l, pnew)
997 
998  fact = lref*qdt/tref
999  invgam = 1.d0/(rc_gamma-1.d0)
1000 
1001  {do ix^db = ixo^lim^db\}
1002  plocal = ptherm(ix^d)
1003  rholocal = wct(ix^d, rho_)
1004  tlocal1 = max(plocal/rholocal, smalldouble)
1005 
1006  emin = w(ix^d, rho_) * tlow * invgam
1007  lmax = max(zero, ( pnew(ix^d)*invgam - emin ) / qdt)
1008 
1009  ! No cooling if temperature is below floor level.
1010  ! Assuming Bremmstrahlung if temperature is higher than maximum.
1011  if( tlocal1 <= tcoolmin) then
1012  l1 = zero
1013  else if( tlocal1 >= tcoolmax ) then
1014  l1 = lcool(ncool) * sqrt(tlocal1 / tcoolmax)
1015  l1 = l1 * (rholocal**2)
1016  l1 = min(l1, lmax)
1017  else
1018  call findl(tlocal1, l1)
1019  call findy(tlocal1, y1)
1020  tc = tlocal1 * invgam / (rholocal * l1)
1021  y2 = y1 + (tlocal1 * fact) / (l1 * tc)
1022  call findt(tlocal2, y2)
1023 
1024  if( tlocal2 <= tcoolmin ) then
1025  l1 = lmax
1026  else
1027  l1 = (tlocal1 - tlocal2) * invgam / (rholocal * qdt)
1028  endif
1029 
1030  l1 = l1 * (rholocal**2)
1031  l1 = min(l1, lmax)
1032  endif
1033  coolrate(ix^d) = l1
1034  {enddo^d&\}
1035 
1036  end subroutine getvar_cooling_exact
1037 
1038  subroutine radiative_cooling_add_source(qdt,ixI^L,ixO^L,wCT,w,x,&
1039  qsourcesplit,active)
1041  ! w[iw]=w[iw]+qdt*S[wCT,x] where S is the source based on wCT within ixO
1043 
1044  integer, intent(in) :: ixI^L, ixO^L
1045  double precision, intent(in) :: qdt, x(ixi^s,1:ndim), wCT(ixi^s,1:nw)
1046  double precision, intent(inout) :: w(ixi^s,1:nw)
1047  logical, intent(in) :: qsourcesplit
1048  logical, intent(inout) :: active
1049 
1050  if(qsourcesplit .eqv. rc_split) then
1051  active = .true.
1052  select case(coolmethod)
1053  case ('explicit1')
1054  if(mype==0)then
1055  if(it==1) then
1056  write(*,*)'Fully explicit cooling is not completely safe in this version'
1057  write(*,*)'PROCEED WITH CAUTION!'
1058  endif
1059  endif
1060  call cool_explicit1(qdt,ixi^l,ixo^l,wct,w,x)
1061  case ('explicit2')
1062  call cool_explicit2(qdt,ixi^l,ixo^l,wct,w,x)
1063  case ('semiimplicit')
1064  call cool_semiimplicit(qdt,ixi^l,ixo^l,wct,w,x)
1065  case ('implicit')
1066  call cool_implicit(qdt,ixi^l,ixo^l,wct,w,x)
1067  case ('exact')
1068  call cool_exact(qdt,ixi^l,ixo^l,wct,w,x)
1069  case default
1070  call mpistop("This cooling method is unknown")
1071  end select
1072  if( tfix ) call floortemperature(qdt,ixi^l,ixo^l,wct,w,x)
1073  end if
1074 
1075  end subroutine radiative_cooling_add_source
1076 
1077  subroutine floortemperature(qdt,ixI^L,ixO^L,wCT,w,x)
1078  ! Force minimum temperature to a fixed temperature
1080 
1081  integer, intent(in) :: ixI^L, ixO^L
1082  double precision, intent(in) :: qdt, x(ixi^s,1:ndim), wCT(ixi^s,1:nw)
1083  double precision, intent(inout) :: w(ixi^s,1:nw)
1084 
1085  double precision :: etherm(ixi^s), emin, tfloor
1086  integer :: ix^D
1087 
1088  tfloor = tlow
1089 
1090  call phys_get_pthermal(w,x,ixi^l,ixo^l,etherm)
1091 
1092  {do ix^db = ixo^lim^db\}
1093  emin = w(ix^d,rho_)*tfloor/(rc_gamma-1.d0)
1094  etherm(ix^d) = etherm(ix^d)/(rc_gamma-1.d0)
1095  if( etherm(ix^d) < emin ) w(ix^d,e_)=w(ix^d,e_)-etherm(ix^d)+emin
1096  {enddo^d&\}
1097 
1098  end subroutine floortemperature
1099 
1100  subroutine cool_explicit1(qdt,ixI^L,ixO^L,wCT,w,x)
1101  ! explicit cooling routine that depends on getdt to
1102  ! adjust the timestep. Accurate but incredibly slow
1104 
1105  integer, intent(in) :: ixI^L, ixO^L
1106  double precision, intent(in) :: qdt, x(ixi^s,1:ndim), wCT(ixi^s,1:nw)
1107  double precision, intent(inout) :: w(ixi^s,1:nw)
1108 
1109  double precision :: L1,Tlocal1, ptherm(ixi^s),pnew(ixi^s)
1110  double precision :: plocal, rholocal
1111  double precision :: emin, Lmax
1112 
1113  integer :: ix^D
1114  integer :: icool
1115 
1116  call phys_get_pthermal(wct,x,ixi^l,ixo^l,ptherm)
1117  call phys_get_pthermal(w,x,ixi^l,ixo^l,pnew)
1118 
1119  {do ix^db = ixo^lim^db\}
1120  plocal = ptherm(ix^d)
1121  rholocal = wct(ix^d,rho_)
1122  emin = rholocal*tlow/(rc_gamma-1.d0)
1123  lmax = max(zero,pnew(ix^d)/(rc_gamma-1.d0)-emin)/qdt
1124  ! Tlocal = P/rho
1125  tlocal1 = max(plocal/(rholocal),smalldouble)
1126  !
1127  ! Determine explicit cooling
1128  !
1129  ! If temperature is below floor level, no cooling.
1130  ! Stop wasting time and go to next gridpoint.
1131  ! If the temperature is higher than the maximum,
1132  ! assume Bremmstrahlung
1133  if( tlocal1<=tcoolmin ) then
1134  l1 = zero
1135  else if( tlocal1>=tcoolmax )then
1136  l1 = lcool(ncool)*sqrt(tlocal1/tcoolmax)
1137  l1 = l1*(rholocal**2)
1138  l1 = min(l1,lmax)
1139  w(ix^d,e_) = w(ix^d,e_)-l1*qdt
1140  else
1141  call findl(tlocal1,l1)
1142  l1 = l1*(rholocal**2)
1143  l1 = min(l1,lmax)
1144  w(ix^d,e_) = w(ix^d,e_)-l1*qdt
1145 
1146  endif
1147  {enddo^d&\}
1148 
1149  end subroutine cool_explicit1
1150 
1151  subroutine cool_explicit2(qdt,ixI^L,ixO^L,wCT,w,x)
1152  !
1153  ! explicit cooling routine that does a series
1154  ! of small forward integration steps, to make
1155  ! sure the amount of cooling remains correct
1156  !
1157  ! Not as accurate as 'explicit1', but a lot faster
1158  ! tends to overestimate cooling
1159 
1161 
1162  integer, intent(in) :: ixI^L, ixO^L
1163  double precision, intent(in) :: qdt, x(ixi^s,1:ndim), wCT(ixi^s,1:nw)
1164  double precision, intent(inout) :: w(ixi^s,1:nw)
1165 
1166  double precision :: Ltest, etherm, de
1167  double precision :: dtmax, dtstep
1168 
1169  integer :: idt,ndtstep
1170 
1171  double precision :: L1,Tlocal1,ptherm(ixi^s),pnew(ixi^s)
1172  double precision :: plocal, rholocal
1173  double precision :: emin, Lmax
1174 
1175  integer :: ix^D
1176  integer :: icool
1177 
1178  call phys_get_pthermal(wct,x,ixi^l,ixo^l,ptherm)
1179  call phys_get_pthermal(w,x,ixi^l,ixo^l,pnew )
1180 
1181 
1182  {do ix^db = ixo^lim^db\}
1183  ! Calculate explicit cooling value
1184  dtmax = qdt
1185  plocal = ptherm(ix^d)
1186  etherm = plocal/(rc_gamma-1.d0)
1187 
1188  rholocal = wct(ix^d,rho_)
1189  emin = rholocal*tlow/(rc_gamma-1.d0)
1190  lmax = max(zero,(pnew(ix^d)/(rc_gamma-1.d0))-emin)/qdt
1191  ! Tlocal = P/rho
1192  tlocal1 = plocal/(rholocal)
1193  !
1194  ! Determine explicit cooling
1195  !
1196  ! If temperature is below floor level, no cooling.
1197  ! Stop wasting time and go to next gridpoint.
1198  ! If the temperature is higher than the maximum,
1199  ! assume Bremmstrahlung
1200  if( tlocal1<=tcoolmin ) then
1201  ltest = zero
1202  else if( tlocal1>=tcoolmax )then
1203  ltest = lcool(ncool)*sqrt(tlocal1/tcoolmax)
1204  ltest = l1*(rholocal**2)
1205  ltest = min(l1,lmax)
1206  if( dtmax>cfrac*etherm/ltest) dtmax = cfrac*etherm/ltest
1207  else
1208  call findl(tlocal1,ltest)
1209  ltest = ltest*(rholocal**2)
1210  ltest = min(ltest,lmax)
1211  if( dtmax>cfrac*etherm/ltest) dtmax = cfrac*etherm/ltest
1212  endif
1213  !
1214  ! Calculate number of steps for cooling
1215  !
1216  ndtstep = max(nint(qdt/dtmax),1)+1
1217  dtstep = qdt/ndtstep
1218  !
1219  ! Use explicit cooling value for first step
1220  !
1221  de = ltest*dtstep
1222  etherm = etherm - de
1223 
1224  do idt=2,ndtstep
1225  plocal = etherm*(rc_gamma-1.d0)
1226  lmax = max(zero,etherm-emin)/dtstep
1227  ! Tlocal = P/rho
1228  tlocal1 = plocal/(rholocal)
1229  if( tlocal1<=tcoolmin ) then
1230  l1 = zero
1231  exit
1232  else if( tlocal1>=tcoolmax )then
1233  l1 = lcool(ncool)*sqrt(tlocal1/tcoolmax)
1234  l1 = l1*(rholocal**2)
1235  l1 = min(l1,lmax)
1236  else
1237  call findl(tlocal1,l1)
1238  l1 = l1*(rholocal**2)
1239  l1 = min(l1,lmax)
1240  endif
1241 
1242  de = de + l1*dtstep
1243  etherm = etherm - l1*dtstep
1244  enddo
1245  w(ix^d,e_) = w(ix^d,e_) -de
1246  {enddo^d&\}
1247 
1248  end subroutine cool_explicit2
1249 
1250  subroutine cool_semiimplicit(qdt,ixI^L,ixO^L,wCT,w,x)
1251  !
1252  ! Semi-implicit cooling method based on a two point
1253  ! average
1254  !
1255  ! Fast, but tends to underestimate cooling
1257 
1258  integer, intent(in) :: ixI^L, ixO^L
1259  double precision, intent(in) :: qdt, x(ixi^s,1:ndim), wCT(ixi^s,1:nw)
1260  double precision, intent(inout) :: w(ixi^s,1:nw)
1261 
1262  double precision :: L1,L2,Tlocal1, Tlocal2
1263  double precision :: etemp
1264 
1265  double precision :: plocal, rholocal
1266  double precision :: emin, Lmax
1267 
1268  double precision :: ptherm(ixi^s),pnew(ixi^s)
1269  integer :: ix^D
1270 
1271  call phys_get_pthermal(wct,x,ixi^l,ixo^l,ptherm)
1272  call phys_get_pthermal(w,x,ixi^l,ixo^l,pnew )
1273 
1274  {do ix^db = ixo^lim^db\}
1275  plocal = ptherm(ix^d)
1276  rholocal = wct(ix^d,rho_)
1277  emin = rholocal*tlow/(rc_gamma-1.d0)
1278  lmax = max(zero,pnew(ix^d)/(rc_gamma-1.d0)-emin)/qdt
1279  ! Tlocal = P/rho
1280  tlocal1 = max(plocal/(rholocal),smalldouble)
1281  !
1282  ! Determine explicit cooling at present temperature
1283  !
1284  ! If temperature is below floor level, no cooling.
1285  ! Stop wasting time and go to next gridpoint.
1286  ! If the temperature is higher than the maximum,
1287  ! assume Bremmstrahlung
1288  if( tlocal1<=tcoolmin ) then
1289  l1 = zero
1290  l2 = zero
1291  else
1292  if( tlocal1>=tcoolmax )then
1293  l1 = lcool(ncool)*sqrt(tlocal1/tcoolmax)
1294  else
1295  call findl(tlocal1,l1)
1296  end if
1297  l1 = l1*(rholocal**2)
1298  etemp = plocal/(rc_gamma-1.d0) - l1*qdt
1299  tlocal2 = etemp*(rc_gamma-1.d0)/(rholocal)
1300  !
1301  ! Determine explicit cooling at new temperature
1302  !
1303  if( tlocal2<=tcoolmin ) then
1304  l2 = zero
1305  else if( tlocal2>=tcoolmax )then
1306  l2 = lcool(ncool)*sqrt(tlocal2/tcoolmax)
1307  else
1308  call findl(tlocal2,l2)
1309  end if
1310  l2 = l2*(rholocal**2)
1311  w(ix^d,e_) = w(ix^d,e_) - min(half*(l1+l2),lmax)*qdt
1312  endif
1313  {enddo^d&\}
1314 
1315  end subroutine cool_semiimplicit
1316 
1317  subroutine cool_implicit(qdt,ixI^L,ixO^L,wCT,w,x)
1318  !
1319  ! Implicit cooling method based on a half-step
1320  ! refinement algorithm
1321  !
1323 
1324  integer, intent(in) :: ixI^L, ixO^L
1325  double precision, intent(in) :: qdt, x(ixi^s,1:ndim), wCT(ixi^s,1:nw)
1326  double precision, intent(inout) :: w(ixi^s,1:nw)
1327 
1328  double precision :: Ltemp,Tlocal1,Tnew,f1,f2,ptherm(ixi^s), pnew(ixi^s)
1329 
1330  double precision :: plocal, rholocal, elocal
1331  double precision :: emin, Lmax, eold, enew, estep
1332  integer, parameter :: maxiter = 100
1333  double precision, parameter :: e_error = 1.0d-6
1334 
1335  integer :: ix^D, j
1336 
1337  call phys_get_pthermal(wct,x,ixi^l,ixo^l,ptherm)
1338  call phys_get_pthermal(w,x,ixi^l,ixo^l,pnew )
1339 
1340 
1341  {do ix^db = ixo^lim^db\}
1342  plocal = ptherm(ix^d)
1343  elocal = plocal/(rc_gamma-1.d0)
1344  rholocal = wct(ix^d,rho_)
1345  emin = rholocal*tlow/(rc_gamma-1.d0)
1346  lmax = max(zero,pnew(ix^d)/(rc_gamma-1.d0)-emin)/qdt
1347  ! Tlocal = P/rho
1348  tlocal1 = max(plocal/(rholocal),smalldouble)
1349  !
1350  ! Determine explicit cooling at present temperature
1351  !
1352  ! If temperature is below floor level, no cooling.
1353  ! Stop wasting time and go to next gridpoint.
1354  ! If the temperature is higher than the maximum,
1355  ! assume Bremmstrahlung
1356  if( tlocal1<=tcoolmin ) then
1357  ltemp = zero
1358  else
1359  eold = elocal
1360  enew = elocal
1361  estep = -(smalldouble)
1362  f2 = 1.d0
1363  do j=1,maxiter+1
1364  if( j>maxiter ) call mpistop("Implicit cooling exceeds maximum iterations")
1365  tnew = enew*(rc_gamma-1.d0)/(rholocal)
1366  if( tnew<=tcoolmin ) then
1367  ltemp = lmax
1368  exit
1369  else if( tnew>=tcoolmax )then
1370  ltemp = lcool(ncool)*sqrt(tnew/tcoolmax)
1371  else
1372  call findl(tnew,ltemp)
1373  end if
1374  ltemp = ltemp*(rholocal**2)
1375  eold = enew + ltemp*qdt
1376 
1377  f1 = elocal -eold
1378  if( abs(half*f1/(elocal+eold)) < e_error ) exit
1379 
1380  if(j==1) estep = max((elocal-emin)*half,smalldouble)
1381  if(f1*f2 < zero ) estep = -half*estep
1382  f2 = f1
1383  enew = enew +estep
1384  enddo
1385  endif
1386  w(ix^d,e_) = w(ix^d,e_) - min(ltemp,lmax)*qdt
1387  {enddo^d&\}
1388 
1389  end subroutine cool_implicit
1390 
1391  subroutine cool_exact(qdt,ixI^L,ixO^L,wCT,w,x)
1392  !
1393  ! Cooling routine using exact integration method from Townsend 2009
1394  !
1396 
1397  integer, intent(in) :: ixI^L, ixO^L
1398  double precision, intent(in) :: qdt, x(ixi^s,1:ndim), wCT(ixi^s,1:nw)
1399  double precision, intent(inout) :: w(ixi^s,1:nw)
1400 
1401  double precision :: Y1, tc, Y2
1402  double precision :: L1,Tlocal1, ptherm(ixi^s), Tlocal2, pnew(ixi^s)
1403  double precision :: plocal, rholocal, invgam
1404  double precision :: emin, Lmax, fact
1405 
1406  integer :: ix^D
1407  integer :: icool
1408 
1409  call phys_get_pthermal(wct,x,ixi^l,ixo^l,ptherm)
1410  call phys_get_pthermal(w,x,ixi^l,ixo^l,pnew )
1411 
1412  fact = lref*qdt/tref
1413 
1414  invgam=1.d0/(rc_gamma-1.d0)
1415  {do ix^db = ixo^lim^db\}
1416  plocal = ptherm(ix^d)
1417  rholocal = wct(ix^d,rho_)
1418 
1419  emin = w(ix^d,rho_)*tlow*invgam
1420  lmax = max(zero,(pnew(ix^d)*invgam-emin)/qdt)
1421 
1422  ! Tlocal = P/rho
1423  tlocal1 = max(plocal/rholocal,smalldouble)
1424  !
1425  ! Determine explicit cooling
1426  !
1427  ! If temperature is below floor level, no cooling.
1428  ! Stop wasting time and go to next gridpoint.
1429  ! If the temperature is higher than the maximum,
1430  ! assume Bremmstrahlung
1431  if( tlocal1<=tcoolmin ) then
1432  l1 = zero
1433  else if( tlocal1>=tcoolmax )then
1434  l1 = lcool(ncool)*sqrt(tlocal1/tcoolmax)
1435  l1 = l1*(rholocal**2)
1436  l1 = min(l1,lmax)
1437  w(ix^d,e_) = w(ix^d,e_)-l1*qdt
1438  else
1439  call findl(tlocal1,l1)
1440  call findy(tlocal1,y1)
1441  tc = tlocal1*invgam/(rholocal*l1)
1442  y2 = y1+(tlocal1*fact)/(l1*tc)
1443  call findt(tlocal2,y2)
1444 
1445  if(tlocal2<=tcoolmin) then
1446  l1 = lmax
1447  else
1448  l1 = (tlocal1-tlocal2)*invgam/(rholocal*qdt)
1449  endif
1450 
1451  l1 = l1*(rholocal**2)
1452  l1 = min(l1,lmax)
1453  w(ix^d,e_) = w(ix^d,e_)-l1*qdt
1454  endif
1455  {enddo^d&\}
1456 
1457  end subroutine cool_exact
1458 
1459  subroutine findl (tpoint,Lpoint)
1460  ! Fast search option to find correct point
1461  ! in cooling curve
1463 
1464  double precision,intent(IN) :: tpoint
1465  double precision, intent(OUT) :: Lpoint
1466  integer :: ipoint
1467  integer :: jl,jc,jh
1468  double precision :: lgtp
1469 
1470  lgtp = dlog10(tpoint)
1471  jl = int((lgtp - lgtcoolmin) / lgstep) + 1
1472  lpoint = lcool(jl)+ (tpoint-tcool(jl)) &
1473  * (lcool(jl+1)-lcool(jl)) &
1474  / (tcool(jl+1)-tcool(jl))
1475 
1476 ! if (tpoint == tcoolmin) then
1477 ! Lpoint = Lcool(1)
1478 ! else if (tpoint == tcoolmax) then
1479 ! Lpoint = Lcool(ncool)
1480 ! else
1481 ! jl=0
1482 ! jh=ncool+1
1483 ! do
1484 ! if (jh-jl <= 1) exit
1485 ! jc=(jh+jl)/2
1486 ! if (tpoint >= tcool(jc)) then
1487 ! jl=jc
1488 ! else
1489 ! jh=jc
1490 ! end if
1491 ! end do
1492 ! ! Linear interpolation to obtain correct cooling
1493 ! Lpoint = Lcool(jl)+ (tpoint-tcool(jl)) &
1494 ! * (Lcool(jl+1)-Lcool(jl)) &
1495 ! / (tcool(jl+1)-tcool(jl))
1496 ! end if
1497 
1498  end subroutine findl
1499 
1500  subroutine findy (tpoint,Ypoint)
1502  ! Fast search option to find correct point in cooling time
1503 
1505 
1506  double precision,intent(IN) :: tpoint
1507  double precision, intent(OUT) :: Ypoint
1508  integer :: ipoint
1509  integer :: jl,jc,jh
1510  double precision :: lgtp
1511 
1512  lgtp = dlog10(tpoint)
1513  jl = int((lgtp - lgtcoolmin) / lgstep) + 1
1514  ypoint = yc(jl)+ (tpoint-tcool(jl)) &
1515  * (yc(jl+1)-yc(jl)) &
1516  / (tcool(jl+1)-tcool(jl))
1517 
1518  ! integer i
1519  !
1520  ! if (tpoint == tcoolmin) then
1521  ! Ypoint = Yc(1)
1522  ! else if (tpoint == tcoolmax) then
1523  ! Ypoint = Yc(ncool)
1524  ! else
1525  ! jl=0
1526  ! jh=ncool+1
1527  ! do
1528  ! if (jh-jl <= 1) exit
1529  ! jc=(jh+jl)/2
1530  ! if (tpoint >= tcool(jc)) then
1531  ! jl=jc
1532  ! else
1533  ! jh=jc
1534  ! end if
1535  ! end do
1536  ! ! Linear interpolation to obtain correct value
1537  ! Ypoint = Yc(jl)+ (tpoint-tcool(jl)) &
1538  ! * (Yc(jl+1)-Yc(jl)) &
1539  ! / (tcool(jl+1)-tcool(jl))
1540  ! end if
1541 
1542  end subroutine findy
1543 
1544  subroutine findt (tpoint,Ypoint)
1546  ! Fast search option to find correct temperature
1547  ! from cooling time. Only possible this way because T is a monotonously
1548  ! decreasing function
1550 
1551  double precision,intent(OUT) :: tpoint
1552  double precision, intent(IN) :: Ypoint
1553  integer :: ipoint
1554  integer :: jl,jc,jh
1555  integer i
1556 
1557  if (ypoint >= yc(1)) then
1558  tpoint = tcoolmin
1559  else if (ypoint == yc(ncool)) then
1560  tpoint = tcoolmax
1561  else
1562  jl=0
1563  jh=ncool+1
1564  do
1565  if (jh-jl <= 1) exit
1566  jc=(jh+jl)/2
1567  if (ypoint <= yc(jc)) then
1568  jl=jc
1569  else
1570  jh=jc
1571  end if
1572  end do
1573  ! Linear interpolation to obtain correct temperature
1574  tpoint = tcool(jl)+ (ypoint-yc(jl)) &
1575  * (tcool(jl+1)-tcool(jl)) &
1576  / (yc(jl+1)-yc(jl))
1577  end if
1578  end subroutine findt
1579 
1580  subroutine finddldt (tpoint,dLpoint)
1582  ! Fast search option to find correct point
1583  ! in derivative of cooling curve
1585 
1586  double precision,intent(IN) :: tpoint
1587  double precision, intent(OUT) :: dLpoint
1588  integer :: ipoint
1589  integer :: jl,jc,jh
1590  double precision :: lgtp
1591 
1592  lgtp = dlog10(tpoint)
1593  jl = int((lgtp - lgtcoolmin) / lgstep) + 1
1594  dlpoint = dldtcool(jl)+ (tpoint-tcool(jl)) &
1595  * (dldtcool(jl+1)-dldtcool(jl)) &
1596  / (tcool(jl+1)-tcool(jl))
1597 
1598 ! if (tpoint == tcoolmin) then
1599 ! dLpoint = dLdtcool(1)
1600 ! else if (tpoint == tcoolmax) then
1601 ! dLpoint = dLdtcool(ncool)
1602 ! else
1603 ! jl=0
1604 ! jh=ncool+1
1605 ! do
1606 ! if (jh-jl <= 1) exit
1607 ! jc=(jh+jl)/2
1608 ! if (tpoint >= tcool(jc)) then
1609 ! jl=jc
1610 ! else
1611 ! jh=jc
1612 ! end if
1613 ! end do
1614 ! ! Linear interpolation to obtain correct cooling derivative
1615 ! dLpoint = dLdtcool(jl)+ (tpoint-tcool(jl)) &
1616 ! * (dLdtcool(jl+1)-dLdtcool(jl)) &
1617 ! / (tcool(jl+1)-tcool(jl))
1618 ! end if
1619  end subroutine finddldt
1620 
1621 end module mod_radiative_cooling
This module contains definitions of global parameters and variables and some generic functions/subrou...
double precision, dimension(1:110) l_spex
double precision, dimension(1:71) l_mlcosmol
double precision, dimension(1:151) t_cl_solar
double precision, dimension(:), allocatable lcool
subroutine getvar_cooling_exact(qdt, ixIL, ixOL, wCT, w, x, coolrate)
subroutine radiative_cooling_add_source(qdt, ixIL, ixOL, wCT, w, x, qsourcesplit, active)
integer ndir
Number of spatial dimensions (components) for vector variables.
subroutine cool_semiimplicit(qdt, ixIL, ixOL, wCT, w, x)
subroutine finddldt(tpoint, dLpoint)
double precision, dimension(1:71) t_mlwc
subroutine cool_exact(qdt, ixIL, ixOL, wCT, w, x)
subroutine cooling_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
double precision, dimension(1:71) l_dm
subroutine floortemperature(qdt, ixIL, ixOL, wCT, w, x)
double precision, dimension(1:45) t_jccorona
subroutine cool_explicit1(qdt, ixIL, ixOL, wCT, w, x)
double precision, dimension(1:71) l_mlsolar1
double precision, dimension(1:71) l_mlwc
subroutine rc_params_read(files)
Read this module"s parameters from a file.
subroutine findy(tpoint, Ypoint)
subroutine radiative_cooling_init(phys_gamma, He_abund)
Radiative cooling initialization.
double precision, dimension(1:71) t_mlcosmol
double precision, dimension(1:76) t_dm_2
subroutine findt(tpoint, Ypoint)
double precision, dimension(1:110) nenh_spex
subroutine cool_implicit(qdt, ixIL, ixOL, wCT, w, x)
subroutine cool_explicit2(qdt, ixIL, ixOL, wCT, w, x)
double precision, dimension(:), allocatable tcool
integer, parameter unitpar
double precision, dimension(1:110) t_spex
integer, dimension(:), allocatable, parameter d
integer mype
The rank of the current MPI task.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
subroutine findl(tpoint, Lpoint)
double precision unit_pressure
Physical scaling factor for pressure.
double precision unit_temperature
Physical scaling factor for temperature.
double precision, dimension(1:51) t_mb
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: comm_lib.t:198
double precision, dimension(1:151) l_cl_solar
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision, dimension(:), allocatable invyc
subroutine getvar_cooling(ixIL, ixOL, w, x, coolrate)
double precision, dimension(1:45) l_jccorona
double precision, dimension(:), allocatable dldtcool
double precision, dimension(1:76) l_dm_2
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
logical rc_split
Add cooling source in a split way (.true.) or un-split way (.false.)
module radiative cooling – add optically thin radiative cooling for HD and MHD
double precision unit_time
Physical scaling factor for time.
integer it
Number of time steps taken.
double precision, dimension(1:51) l_mb
double precision, dimension(1:151) l_cl_ism
double precision, dimension(1:151) t_cl_ism
double precision unit_numberdensity
Physical scaling factor for number density.
procedure(sub_get_pthermal), pointer phys_get_pthermal
Definition: mod_physics.t:59
double precision, dimension(1:71) t_dm
double precision, dimension(:), allocatable yc
double precision, dimension(1:71) t_mlsolar1