MPI-AMRVAC 3.2
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
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 ionized 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 or a piecewise power law
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)
12!> Included piecewise power law functionality by Joris Hermans (01.2020)
13module mod_radiative_cooling
14! For the interpolatable tables: 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 Bremsstrahlung
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 use mod_global_parameters, only: std_len
23 use mod_physics
24 use mod_comm_lib, only: mpistop
25 implicit none
26
27 !> Helium abundance over Hydrogen
28 double precision, private :: He_abundance
29
30 !> The adiabatic index
31 double precision, private :: rc_gamma
32
33 !> The adiabatic index minus 1
34 double precision, private :: rc_gamma_1
35
36 !> inverse of the adiabatic index minus 1
37 double precision, private :: invgam
38
39 abstract interface
40 subroutine get_subr1(w,x,ixI^L,ixO^L,res)
41 use mod_global_parameters
42 integer, intent(in) :: ixI^L, ixO^L
43 double precision, intent(in) :: w(ixI^S,nw)
44 double precision, intent(in) :: x(ixI^S,1:ndim)
45 double precision, intent(out):: res(ixI^S)
46 end subroutine get_subr1
47 end interface
48
49 type rc_fluid
50
51 double precision :: rad_cut_hgt
52 double precision :: rad_cut_dey
53
54 ! these are set in init method
55 double precision, allocatable :: tcool(:), lcool(:), dldtcool(:)
56 double precision, allocatable :: yc(:), invyc(:)
57 double precision :: tref, lref, tcoolmin,tcoolmax
58 double precision :: lgtcoolmin, lgtcoolmax, lgstep
59
60 ! The piecewise powerlaw (PPL) tabels and variabels
61 ! x_* en t_* are given as log_10
62 double precision, allocatable :: y_ppl(:), t_ppl(:), l_ppl(:), a_ppl(:)
63
64 !> Lower limit of temperature
65 double precision :: tlow
66
67 !> Coefficent of cooling time step
68 double precision :: cfrac
69
70 !> Index of the energy density
71 integer :: e_
72 !> Index of cut off temperature for TRAC
73 integer :: tcoff_
74
75 ! these are set as parameters
76 !> Resolution of temperature in interpolated tables
77 integer :: ncool
78
79 integer :: n_ppl
80
81 !> Fixed temperature not lower than tlow
82 logical :: tfix
83
84 !> Add cooling source in a split way (.true.) or un-split way (.false.)
85 logical :: rc_split
86
87 logical :: isppl = .false.
88
89 !> cutoff radiative cooling below rad_cut_hgt
90 logical :: rad_cut
91 ! these are to be set directly
92 logical :: has_equi = .false.
93
94 !> Name of cooling curve
95 character(len=std_len) :: coolcurve
96
97 !> Name of cooling method
98 character(len=std_len) :: coolmethod
99
100 procedure(get_subr1), pointer, nopass :: get_rho => null()
101 procedure(get_subr1), pointer, nopass :: get_rho_equi => null()
102 procedure(get_subr1), pointer, nopass :: get_pthermal => null()
103 procedure(get_subr1), pointer, nopass :: get_pthermal_equi => null()
104 procedure(get_subr1), pointer, nopass :: get_var_rfactor => null()
105
106 end type rc_fluid
107
108 double precision :: t_hildner(1:6), t_fm(1:5), t_rosner(1:10), t_klimchuk(1:8), &
109 t_spex_dm_rough(1:8), t_spex_dm_fine(1:15)
110
111 double precision :: x_hildner(1:5), x_fm(1:4), x_rosner(1:9), x_klimchuk(1:7), &
112 x_spex_dm_rough(1:7), x_spex_dm_fine(1:14)
113
114 double precision :: a_hildner(1:5), a_fm(1:4), a_rosner(1:9), a_klimchuk(1:7), &
115 a_spex_dm_rough(1:7), a_spex_dm_fine(1:14)
116
117 integer :: n_hildner, n_fm, n_rosner, n_klimchuk, n_spex_dm_rough, n_spex_dm_fine
118
119 data n_hildner / 5 /
120
121 data t_hildner / 3.00000, 4.17609, 4.90309, 5.47712, 5.90309, 10.00000 /
122
123 data x_hildner / -53.30803, -29.92082, -21.09691, -7.40450, -16.25885 /
124
125 data a_hildner / 7.4, 1.8, 0.0, -2.5, -1.0 /
126
127 data n_fm / 4 /
128
129 data t_fm / 3.00000, 4.30103, 5.60206, 7.00000, 10.00000 /
130
131 data x_fm / -31.15813, -27.50227, -18.25885, -35.75893 /
132
133 data a_fm / 2.0, 1.15, -0.5, 2.0 /
134
135 data n_rosner / 9 /
136
137 data t_rosner / 3.00000, 3.89063, 4.30195, 4.57500, 4.90000, &
138 5.40000, 5.77000, 6.31500, 7.60457, 10.00000 /
139
140 data x_rosner / -69.900, -48.307, -21.850, -31.000, -21.200, &
141 -10.400, -21.940, -17.730, -26.602 /
142
143 data a_rosner / 11.70, 6.15, 0.00, 2.00, 0.00, &
144 -2.00, 0.00, -0.67, 0.50 /
145
146 data n_klimchuk / 7 /
147
148 data t_klimchuk / 3.00, 4.97, 5.67, 6.18, 6.55, 6.90, 7.63, 10.00 /
149
150 data x_klimchuk / -30.96257, -16.05208, -21.72125, -12.45223, &
151 -24.46092, -15.26043, -26.70774 /
152
153 data a_klimchuk / 2.00, -1.00, 0.00, -1.50, 0.33, -1.00, 0.50 /
154
155 data n_spex_dm_rough / 7 /
156
157 data t_spex_dm_rough / 1.000, 1.572, 3.992, 4.165, 5.221, 5.751, 7.295, 8.160 /
158
159 data x_spex_dm_rough / -34.286, -28.282, -108.273, -26.662, -9.729, -17.550, -24.767 /
160
161 data a_spex_dm_rough / 4.560, 0.740, 20.777, 1.182, -2.061, -0.701, 0.288 /
162
163 data n_spex_dm_fine / 14 /
164
165 data t_spex_dm_fine / 1.000, 1.422, 2.806, 3.980, 4.177, 4.443, 4.832, 5.397, &
166 5.570, 5.890, 6.232, 6.505, 6.941, 7.385, 8.160 /
167
168 data x_spex_dm_fine / -35.314, -29.195, -26.912, -108.273, -18.971, -32.195, -21.217, &
169 -0.247, -15.415, -19.275, -9.387, -22.476, -17.437, -25.026 /
170
171 data a_spex_dm_fine / 5.452, 1.150, 0.337, 20.777, -0.602, 2.374, 0.102, &
172 -3.784, -1.061, -0.406, -1.992, 0.020, -0.706, 0.321 /
173
174 ! Interpolatable tables
175
176 double precision :: t_dm(1:71), t_mb(1:51), t_mlcosmol(1:71) &
177 , t_mlwc(1:71), t_mlsolar1(1:71), t_spex(1:110) &
178 , t_jccorona(1:45), t_cl_ism(1:151), t_cl_solar(1:151)&
179 , t_dm_2(1:76), t_dere(1:101), t_colgan(1:55)
180
181 double precision :: l_dm(1:71), l_mb(1:51), l_mlcosmol(1:71) &
182 , l_mlwc(1:71), l_mlsolar1(1:71), l_spex(1:110) &
183 , l_jccorona(1:45), l_cl_ism(1:151), l_cl_solar(1:151) &
184 , l_dm_2(1:76), l_dere_corona(1:101), l_dere_photo(1:101)&
185 , l_colgan(1:55)
186
187 double precision :: nenh_spex(1:110)
188
189 integer :: n_dm , n_mb , n_mlcosmol &
190 , n_mlwc , n_mlsolar1, n_spex &
191 , n_jccorona, n_cl_ism , n_cl_solar &
192 , n_dm_2 , n_dere , n_colgan
193
194 data n_jccorona / 45 /
195
196 data t_jccorona / 4.00000, 4.14230, 4.21995, 4.29761, 4.37528, &
197 4.45294, 4.53061, 4.60827, 4.68593, 4.76359, &
198 4.79705, 4.83049, 4.86394, 4.89739, 4.93084, &
199 4.96428, 4.99773, 5.03117, 5.06461, 5.17574, &
200 5.28684, 5.39796, 5.50907, 5.62018, 5.73129, &
201 5.84240, 5.95351, 6.06461, 6.17574, 6.28684, &
202 6.39796, 6.50907, 6.62018, 6.73129, 6.84240, &
203 6.95351, 7.06461, 7.17574, 7.28684, 7.39796, &
204 7.50907, 7.62018, 7.73129, 7.84240, 7.95351 /
205
206 data l_jccorona / -200.18883, -100.78630, -30.60384, -22.68481, -21.76445, &
207 -21.67936, -21.54218, -21.37958, -21.25172, -21.17584, &
208 -21.15783, -21.14491, -21.13527, -21.12837, -21.12485, &
209 -21.12439, -21.12642, -21.12802, -21.12548, -21.08965, &
210 -21.08812, -21.19542, -21.34582, -21.34839, -21.31701, &
211 -21.29072, -21.28900, -21.34104, -21.43122, -21.62448, &
212 -21.86694, -22.02897, -22.08051, -22.06057, -22.01973, &
213 -22.00000, -22.05161, -22.22175, -22.41452, -22.52581, &
214 -22.56914, -22.57486, -22.56151, -22.53969, -22.51490 /
215
216 data n_dm / 71 /
217
218 data t_dm / 2.0, 2.1, 2.2, 2.3, 2.4, &
219 2.5, 2.6, 2.7, 2.8, 2.9, &
220 3.0, 3.1, 3.2, 3.3, 3.4, &
221 3.5, 3.6, 3.7, 3.8, 3.9, &
222 4.0, 4.1, 4.2, 4.3, 4.4, &
223 4.5, 4.6, 4.7, 4.8, 4.9, &
224 5.0, 5.1, 5.2, 5.3, 5.4, &
225 5.5, 5.6, 5.7, 5.8, 5.9, &
226 6.0, 6.1, 6.2, 6.3, 6.4, &
227 6.5, 6.6, 6.7, 6.8, 6.9, &
228 7.0, 7.1, 7.2, 7.3, 7.4, &
229 7.5, 7.6, 7.7, 7.8, 7.9, &
230 8.0, 8.1, 8.2, 8.3, 8.4, &
231 8.5, 8.6, 8.7, 8.8, 8.9, &
232 9.0 /
233
234 data l_dm / -26.523, -26.398, -26.301, -26.222, -26.097 &
235 , -26.011, -25.936, -25.866, -25.807, -25.754 &
236 , -25.708, -25.667, -25.630, -25.595, -25.564 &
237 , -25.534, -25.506, -25.479, -25.453, -25.429 &
238 , -25.407, -23.019, -21.762, -21.742, -21.754 &
239 , -21.730, -21.523, -21.455, -21.314, -21.229 &
240 , -21.163, -21.126, -21.092, -21.060, -21.175 &
241 , -21.280, -21.390, -21.547, -21.762, -22.050 &
242 , -22.271, -22.521, -22.646, -22.660, -22.676 &
243 , -22.688, -22.690, -22.662, -22.635, -22.609 &
244 , -22.616, -22.646, -22.697, -22.740, -22.788 &
245 , -22.815, -22.785, -22.754, -22.728, -22.703 &
246 , -22.680, -22.630, -22.580, -22.530, -22.480 &
247 , -22.430, -22.380, -22.330, -22.280, -22.230 &
248 , -22.180 /
249
250 data n_mb / 51 /
251
252 data t_mb / 4.0, 4.1, 4.2, 4.3, 4.4, &
253 4.5, 4.6, 4.7, 4.8, 4.9, &
254 5.0, 5.1, 5.2, 5.3, 5.4, &
255 5.5, 5.6, 5.7, 5.8, 5.9, &
256 6.0, 6.1, 6.2, 6.3, 6.4, &
257 6.5, 6.6, 6.7, 6.8, 6.9, &
258 7.0, 7.1, 7.2, 7.3, 7.4, &
259 7.5, 7.6, 7.7, 7.8, 7.9, &
260 8.0, 8.1, 8.2, 8.3, 8.4, &
261 8.5, 8.6, 8.7, 8.8, 8.9, &
262 9.0 /
263
264 data l_mb / -23.133, -22.895, -22.548, -22.285, -22.099 &
265 , -21.970, -21.918, -21.826, -21.743, -21.638 &
266 , -21.552, -21.447, -21.431, -21.418, -21.461 &
267 , -21.570, -21.743, -21.832, -21.908, -21.981 &
268 , -22.000, -21.998, -21.992, -22.013, -22.095 &
269 , -22.262, -22.397, -22.445, -22.448, -22.446 &
270 , -22.448, -22.465, -22.575, -22.725, -22.749 &
271 , -22.768, -22.753, -22.717, -22.678, -22.637 &
272 , -22.603, -22.553, -22.503, -22.453, -22.403 &
273 , -22.353, -22.303, -22.253, -22.203, -22.153 &
274 , -22.103 /
275
276 data n_mlcosmol / 71 /
277
278 data t_mlcosmol / &
279 2.0, 2.1, 2.2, 2.3, 2.4, &
280 2.5, 2.6, 2.7, 2.8, 2.9, &
281 3.0, 3.1, 3.2, 3.3, 3.4, &
282 3.5, 3.6, 3.7, 3.8, 3.9, &
283 4.0, 4.1, 4.2, 4.3, 4.4, &
284 4.5, 4.6, 4.7, 4.8, 4.9, &
285 5.0, 5.1, 5.2, 5.3, 5.4, &
286 5.5, 5.6, 5.7, 5.8, 5.9, &
287 6.0, 6.1, 6.2, 6.3, 6.4, &
288 6.5, 6.6, 6.7, 6.8, 6.9, &
289 7.0, 7.1, 7.2, 7.3, 7.4, &
290 7.5, 7.6, 7.7, 7.8, 7.9, &
291 8.0, 8.1, 8.2, 8.3, 8.4, &
292 8.5, 8.6, 8.7, 8.8, 8.9, &
293 9.0 /
294
295 data l_mlcosmol / &
296 -99.000, -99.000, -99.000, -99.000, -99.000 &
297 , -99.000, -99.000, -99.000, -99.000, -99.000 &
298 , -99.000, -99.000, -99.000, -99.000, -99.000 &
299 , -99.000, -44.649, -38.362, -33.324, -29.292 &
300 , -26.063, -23.532, -22.192, -22.195, -22.454 &
301 , -22.676, -22.909, -22.925, -22.499, -22.276 &
302 , -22.440, -22.688, -22.917, -23.116, -23.274 &
303 , -23.394, -23.472, -23.516, -23.530, -23.525 &
304 , -23.506, -23.478, -23.444, -23.408, -23.368 &
305 , -23.328, -23.286, -23.244, -23.201, -23.157 &
306 , -23.114, -23.070, -23.026, -22.981, -22.937 &
307 , -22.893, -22.848, -22.803, -22.759, -22.714 &
308 , -22.669, -22.619, -22.569, -22.519, -22.469 &
309 , -22.419, -22.369, -22.319, -22.269, -22.190 &
310 , -22.169 /
311
312 data n_mlwc / 71 /
313
314 data t_mlwc / &
315 2.0, 2.1, 2.2, 2.3, 2.4, &
316 2.5, 2.6, 2.7, 2.8, 2.9, &
317 3.0, 3.1, 3.2, 3.3, 3.4, &
318 3.5, 3.6, 3.7, 3.8, 3.9, &
319 4.0, 4.1, 4.2, 4.3, 4.4, &
320 4.5, 4.6, 4.7, 4.8, 4.9, &
321 5.0, 5.1, 5.2, 5.3, 5.4, &
322 5.5, 5.6, 5.7, 5.8, 5.9, &
323 6.0, 6.1, 6.2, 6.3, 6.4, &
324 6.5, 6.6, 6.7, 6.8, 6.9, &
325 7.0, 7.1, 7.2, 7.3, 7.4, &
326 7.5, 7.6, 7.7, 7.8, 7.9, &
327 8.0, 8.1, 8.2, 8.3, 8.4, &
328 8.5, 8.6, 8.7, 8.8, 8.9, &
329 9.0 /
330
331 data l_mlwc / &
332 -21.192, -21.160, -21.150, -21.150, -21.166 &
333 , -21.191, -21.222, -21.264, -21.308, -21.357 &
334 , -21.408, -21.449, -21.494, -21.544, -21.587 &
335 , -21.638, -21.686, -21.736, -21.780, -21.800 &
336 , -21.744, -21.547, -21.208, -20.849, -20.345 &
337 , -19.771, -19.409, -19.105, -18.827, -18.555 &
338 , -18.460, -18.763, -19.168, -19.334, -19.400 &
339 , -19.701, -20.090, -20.288, -20.337, -20.301 &
340 , -20.233, -20.275, -20.363, -20.508, -20.675 &
341 , -20.856, -21.025, -21.159, -21.256, -21.320 &
342 , -21.354, -21.366, -21.361, -21.343, -21.317 &
343 , -21.285, -21.250, -21.212, -21.172, -21.131 &
344 , -21.089, -21.039, -20.989, -20.939, -20.889 &
345 , -20.839, -20.789, -20.739, -20.689, -20.639 &
346 , -20.589 /
347
348 data n_mlsolar1 / 71 /
349
350 data t_mlsolar1 / &
351 2.0, 2.1, 2.2, 2.3, 2.4, &
352 2.5, 2.6, 2.7, 2.8, 2.9, &
353 3.0, 3.1, 3.2, 3.3, 3.4, &
354 3.5, 3.6, 3.7, 3.8, 3.9, &
355 4.0, 4.1, 4.2, 4.3, 4.4, &
356 4.5, 4.6, 4.7, 4.8, 4.9, &
357 5.0, 5.1, 5.2, 5.3, 5.4, &
358 5.5, 5.6, 5.7, 5.8, 5.9, &
359 6.0, 6.1, 6.2, 6.3, 6.4, &
360 6.5, 6.6, 6.7, 6.8, 6.9, &
361 7.0, 7.1, 7.2, 7.3, 7.4, &
362 7.5, 7.6, 7.7, 7.8, 7.9, &
363 8.0, 8.1, 8.2, 8.3, 8.4, &
364 8.5, 8.6, 8.7, 8.8, 8.9, &
365 9.0 /
366
367 data l_mlsolar1 / &
368 -26.983, -26.951, -26.941, -26.940, -26.956 &
369 , -26.980, -27.011, -27.052, -27.097, -27.145 &
370 , -27.195, -27.235, -27.279, -27.327, -27.368 &
371 , -27.415, -27.456, -27.485, -27.468, -27.223 &
372 , -25.823, -23.501, -22.162, -22.084, -22.157 &
373 , -22.101, -21.974, -21.782, -21.542, -21.335 &
374 , -21.251, -21.275, -21.236, -21.173, -21.167 &
375 , -21.407, -21.670, -21.788, -21.879, -22.008 &
376 , -22.192, -22.912, -22.918, -22.887, -22.929 &
377 , -23.023, -23.094, -23.117, -23.108, -23.083 &
378 , -23.049, -23.011, -22.970, -22.928, -22.885 &
379 , -22.842, -22.798, -22.754, -22.709, -22.665 &
380 , -22.620, -22.570, -22.520, -22.470, -22.420 &
381 , -22.370, -22.320, -22.270, -22.220, -22.170 &
382 , -22.120 /
383
384 data n_spex / 110 /
385
386 data t_spex / &
387 3.80, 3.84, 3.88, 3.92, 3.96, &
388 4.00, 4.04, 4.08, 4.12, 4.16, &
389 4.20, 4.24, 4.28, 4.32, 4.36, &
390 4.40, 4.44, 4.48, 4.52, 4.56, &
391 4.60, 4.64, 4.68, 4.72, 4.76, &
392 4.80, 4.84, 4.88, 4.92, 4.96, &
393 5.00, 5.04, 5.08, 5.12, 5.16, &
394 5.20, 5.24, 5.28, 5.32, 5.36, &
395 5.40, 5.44, 5.48, 5.52, 5.56, &
396 5.60, 5.64, 5.68, 5.72, 5.76, &
397 5.80, 5.84, 5.88, 5.92, 5.96, &
398 6.00, 6.04, 6.08, 6.12, 6.16, &
399 6.20, 6.24, 6.28, 6.32, 6.36, &
400 6.40, 6.44, 6.48, 6.52, 6.56, &
401 6.60, 6.64, 6.68, 6.72, 6.76, &
402 6.80, 6.84, 6.88, 6.92, 6.96, &
403 7.00, 7.04, 7.08, 7.12, 7.16, &
404 7.20, 7.24, 7.28, 7.32, 7.36, &
405 7.40, 7.44, 7.48, 7.52, 7.56, &
406 7.60, 7.64, 7.68, 7.72, 7.76, &
407 7.80, 7.84, 7.88, 7.92, 7.96, &
408 8.00, 8.04, 8.08, 8.12, 8.16 /
409
410 data l_spex / &
411 -25.7331, -25.0383, -24.4059, -23.8288, -23.3027 &
412 , -22.8242, -22.3917, -22.0067, -21.6818, -21.4529 &
413 , -21.3246, -21.3459, -21.4305, -21.5293, -21.6138 &
414 , -21.6615, -21.6551, -21.5919, -21.5092, -21.4124 &
415 , -21.3085, -21.2047, -21.1067, -21.0194, -20.9413 &
416 , -20.8735, -20.8205, -20.7805, -20.7547, -20.7455 &
417 , -20.7565, -20.7820, -20.8008, -20.7994, -20.7847 &
418 , -20.7687, -20.7590, -20.7544, -20.7505, -20.7545 &
419 , -20.7888, -20.8832, -21.0450, -21.2286, -21.3737 &
420 , -21.4573, -21.4935, -21.5098, -21.5345, -21.5863 &
421 , -21.6548, -21.7108, -21.7424, -21.7576, -21.7696 &
422 , -21.7883, -21.8115, -21.8303, -21.8419, -21.8514 &
423 , -21.8690, -21.9057, -21.9690, -22.0554, -22.1488 &
424 , -22.2355, -22.3084, -22.3641, -22.4033, -22.4282 &
425 , -22.4408, -22.4443, -22.4411, -22.4334, -22.4242 &
426 , -22.4164, -22.4134, -22.4168, -22.4267, -22.4418 &
427 , -22.4603, -22.4830, -22.5112, -22.5449, -22.5819 &
428 , -22.6177, -22.6483, -22.6719, -22.6883, -22.6985 &
429 , -22.7032, -22.7037, -22.7008, -22.6950, -22.6869 &
430 , -22.6769, -22.6655, -22.6531, -22.6397, -22.6258 &
431 , -22.6111, -22.5964, -22.5816, -22.5668, -22.5519 &
432 , -22.5367, -22.5216, -22.5062, -22.4912, -22.4753 /
433
434 data nenh_spex / &
435 0.000013264, 0.000042428, 0.000088276, 0.00017967 &
436 , 0.00084362, 0.0034295, 0.013283, 0.042008 &
437 , 0.12138, 0.30481, 0.53386, 0.76622 &
438 , 0.89459, 0.95414, 0.98342 &
439 , 1.0046, 1.0291, 1.0547 &
440 , 1.0767, 1.0888 &
441 , 1.0945, 1.0972, 1.0988 &
442 , 1.1004, 1.1034 &
443 , 1.1102, 1.1233, 1.1433 &
444 , 1.1638, 1.1791 &
445 , 1.1885, 1.1937, 1.1966 &
446 , 1.1983, 1.1993 &
447 , 1.1999, 1.2004, 1.2008 &
448 , 1.2012, 1.2015 &
449 , 1.2020, 1.2025, 1.2030 &
450 , 1.2035, 1.2037 &
451 , 1.2039, 1.2040, 1.2041 &
452 , 1.2042, 1.2044 &
453 , 1.2045, 1.2046, 1.2047 &
454 , 1.2049, 1.2050 &
455 , 1.2051, 1.2053, 1.2055 &
456 , 1.2056, 1.2058 &
457 , 1.2060, 1.2062, 1.2065 &
458 , 1.2067, 1.2070 &
459 , 1.2072, 1.2075, 1.2077 &
460 , 1.2078, 1.2079 &
461 , 1.2080, 1.2081, 1.2082 &
462 , 1.2083, 1.2083 &
463 , 1.2084, 1.2084, 1.2085 &
464 , 1.2085, 1.2086 &
465 , 1.2086, 1.2087, 1.2087 &
466 , 1.2088, 1.2088 &
467 , 1.2089, 1.2089, 1.2089 &
468 , 1.2089, 1.2089 &
469 , 1.2090, 1.2090, 1.2090 &
470 , 1.2090, 1.2090 &
471 , 1.2090, 1.2090, 1.2090 &
472 , 1.2090, 1.2090 &
473 , 1.2090, 1.2090, 1.2090 &
474 , 1.2090, 1.2090 &
475 , 1.2090, 1.2090, 1.2090 &
476 , 1.2090, 1.2090 /
477 !
478 ! To be used together with the SPEX table for the SPEX_DM option
479 ! Assuming an ionization fraction of 10^-3
480 !
481 data n_dm_2 / 76 /
482
483 data t_dm_2 / 1.00, 1.04, 1.08, 1.12, 1.16, 1.20 &
484 , 1.24, 1.28, 1.32, 1.36, 1.40 &
485 , 1.44, 1.48, 1.52, 1.56, 1.60 &
486 , 1.64, 1.68, 1.72, 1.76, 1.80 &
487 , 1.84, 1.88, 1.92, 1.96, 2.00 &
488 , 2.04, 2.08, 2.12, 2.16, 2.20 &
489 , 2.24, 2.28, 2.32, 2.36, 2.40 &
490 , 2.44, 2.48, 2.52, 2.56, 2.60 &
491 , 2.64, 2.68, 2.72, 2.76, 2.80 &
492 , 2.84, 2.88, 2.92, 2.96, 3.00 &
493 , 3.04, 3.08, 3.12, 3.16, 3.20 &
494 , 3.24, 3.28, 3.32, 3.36, 3.40 &
495 , 3.44, 3.48, 3.52, 3.56, 3.60 &
496 , 3.64, 3.68, 3.72, 3.76, 3.80 &
497 , 3.84, 3.88, 3.92, 3.96, 4.00 /
498
499 data l_dm_2 / -30.0377, -29.7062, -29.4055, -29.1331, -28.8864, -28.6631 &
500 , -28.4614, -28.2791, -28.1146, -27.9662, -27.8330 &
501 , -27.7129, -27.6052, -27.5088, -27.4225, -27.3454 &
502 , -27.2767, -27.2153, -27.1605, -27.1111, -27.0664 &
503 , -27.0251, -26.9863, -26.9488, -26.9119, -26.8742 &
504 , -26.8353, -26.7948, -26.7523, -26.7080, -26.6619 &
505 , -26.6146, -26.5666, -26.5183, -26.4702, -26.4229 &
506 , -26.3765, -26.3317, -26.2886, -26.2473, -26.2078 &
507 , -26.1704, -26.1348, -26.1012, -26.0692, -26.0389 &
508 , -26.0101, -25.9825, -25.9566, -25.9318, -25.9083 &
509 , -25.8857, -25.8645, -25.8447, -25.8259, -25.8085 &
510 , -25.7926, -25.7778, -25.7642, -25.7520, -25.7409 &
511 , -25.7310, -25.7222, -25.7142, -25.7071, -25.7005 &
512 , -25.6942, -25.6878, -25.6811, -25.6733, -25.6641 &
513 , -25.6525, -25.6325, -25.6080, -25.5367, -25.4806 /
514
515 data n_cl_solar / 151 /
516
517 data t_cl_solar / &
518 1.000 , 1.050 , 1.100 , 1.150 , &
519 1.200 , 1.250 , 1.300 , 1.350 , &
520 1.400 , 1.450 , 1.500 , 1.550 , &
521 1.600 , 1.650 , 1.700 , 1.750 , &
522 1.800 , 1.850 , 1.900 , 1.950 , &
523 2.000 , 2.050 , 2.100 , 2.150 , &
524 2.200 , 2.250 , 2.300 , 2.350 , &
525 2.400 , 2.450 , 2.500 , 2.550 , &
526 2.600 , 2.650 , 2.700 , 2.750 , &
527 2.800 , 2.850 , 2.900 , 2.950 , &
528 3.000 , 3.050 , 3.100 , 3.150 , &
529 3.200 , 3.250 , 3.300 , 3.350 , &
530 3.400 , 3.450 , 3.500 , 3.550 , &
531 3.600 , 3.650 , 3.700 , 3.750 , &
532 3.800 , 3.850 , 3.900 , 3.950 , &
533 4.000 , 4.050 , 4.100 , 4.150 , &
534 4.200 , 4.250 , 4.300 , 4.350 , &
535 4.400 , 4.450 , 4.500 , 4.550 , &
536 4.600 , 4.650 , 4.700 , 4.750 , &
537 4.800 , 4.850 , 4.900 , 4.950 , &
538 5.000 , 5.050 , 5.100 , 5.150 , &
539 5.200 , 5.250 , 5.300 , 5.350 , &
540 5.400 , 5.450 , 5.500 , 5.550 , &
541 5.600 , 5.650 , 5.700 , 5.750 , &
542 5.800 , 5.850 , 5.900 , 5.950 , &
543 6.000 , 6.050 , 6.100 , 6.150 , &
544 6.200 , 6.250 , 6.300 , 6.350 , &
545 6.400 , 6.450 , 6.500 , 6.550 , &
546 6.600 , 6.650 , 6.700 , 6.750 , &
547 6.800 , 6.850 , 6.900 , 6.950 , &
548 7.000 , 7.050 , 7.100 , 7.150 , &
549 7.200 , 7.250 , 7.300 , 7.350 , &
550 7.400 , 7.450 , 7.500 , 7.550 , &
551 7.600 , 7.650 , 7.700 , 7.750 , &
552 7.800 , 7.850 , 7.900 , 7.950 , &
553 8.000 , 8.100 , 8.200 , 8.300 , &
554 8.400 , 8.500 , 8.600 , 8.700 , &
555 8.800 , 8.900 , 9.00 /
556
557 data l_cl_solar / &
558 -28.375 , -28.251 , -28.137 , -28.029 , &
559 -27.929 , -27.834 , -27.745 , -27.662 , &
560 -27.584 , -27.512 , -27.445 , -27.383 , &
561 -27.326 , -27.273 , -27.223 , -27.175 , &
562 -27.128 , -27.079 , -27.027 , -26.972 , &
563 -26.911 , -26.846 , -26.777 , -26.705 , &
564 -26.632 , -26.554 , -26.479 , -26.407 , &
565 -26.338 , -26.274 , -26.213 , -26.156 , &
566 -26.101 , -26.049 , -25.999 , -25.949 , &
567 -25.901 , -25.852 , -25.803 , -25.754 , &
568 -25.707 , -25.662 , -25.621 , -25.588 , &
569 -25.561 , -25.538 , -25.518 , -25.497 , &
570 -25.475 , -25.452 , -25.426 , -25.400 , &
571 -25.374 , -25.333 , -25.295 , -25.261 , &
572 -25.228 , -25.189 , -25.136 , -25.053 , &
573 -24.888 , -24.454 , -23.480 , -22.562 , &
574 -22.009 , -21.826 , -21.840 , -21.905 , &
575 -21.956 , -21.971 , -21.958 , -21.928 , &
576 -21.879 , -21.810 , -21.724 , -21.623 , &
577 -21.512 , -21.404 , -21.321 , -21.273 , &
578 -21.250 , -21.253 , -21.275 , -21.287 , &
579 -21.282 , -21.275 , -21.272 , -21.267 , &
580 -21.281 , -21.357 , -21.496 , -21.616 , &
581 -21.677 , -21.698 , -21.708 , -21.730 , &
582 -21.767 , -21.793 , -21.794 , -21.787 , &
583 -21.787 , -21.802 , -21.826 , -21.859 , &
584 -21.911 , -21.987 , -22.082 , -22.173 , &
585 -22.253 , -22.325 , -22.392 , -22.448 , &
586 -22.487 , -22.512 , -22.524 , -22.528 , &
587 -22.524 , -22.516 , -22.507 , -22.501 , &
588 -22.502 , -22.511 , -22.533 , -22.565 , &
589 -22.600 , -22.630 , -22.648 , -22.656 , &
590 -22.658 , -22.654 , -22.647 , -22.634 , &
591 -22.619 , -22.602 , -22.585 , -22.566 , &
592 -22.546 , -22.525 , -22.505 , -22.480 , &
593 -22.465 , -22.415 , -22.365 , -22.315 , &
594 -22.265 , -22.215 , -22.165 , -22.115 , &
595 -22.065 , -22.015 , -21.965 /
596
597 data n_cl_ism / 151 /
598
599 data t_cl_ism / &
600 1.000 , 1.050 , 1.100 , 1.150 , &
601 1.200 , 1.250 , 1.300 , 1.350 , &
602 1.400 , 1.450 , 1.500 , 1.550 , &
603 1.600 , 1.650 , 1.700 , 1.750 , &
604 1.800 , 1.850 , 1.900 , 1.950 , &
605 2.000 , 2.050 , 2.100 , 2.150 , &
606 2.200 , 2.250 , 2.300 , 2.350 , &
607 2.400 , 2.450 , 2.500 , 2.550 , &
608 2.600 , 2.650 , 2.700 , 2.750 , &
609 2.800 , 2.850 , 2.900 , 2.950 , &
610 3.000 , 3.050 , 3.100 , 3.150 , &
611 3.200 , 3.250 , 3.300 , 3.350 , &
612 3.400 , 3.450 , 3.500 , 3.550 , &
613 3.600 , 3.650 , 3.700 , 3.750 , &
614 3.800 , 3.850 , 3.900 , 3.950 , &
615 4.000 , 4.050 , 4.100 , 4.150 , &
616 4.200 , 4.250 , 4.300 , 4.350 , &
617 4.400 , 4.450 , 4.500 , 4.550 , &
618 4.600 , 4.650 , 4.700 , 4.750 , &
619 4.800 , 4.850 , 4.900 , 4.950 , &
620 5.000 , 5.050 , 5.100 , 5.150 , &
621 5.200 , 5.250 , 5.300 , 5.350 , &
622 5.400 , 5.450 , 5.500 , 5.550 , &
623 5.600 , 5.650 , 5.700 , 5.750 , &
624 5.800 , 5.850 , 5.900 , 5.950 , &
625 6.000 , 6.050 , 6.100 , 6.150 , &
626 6.200 , 6.250 , 6.300 , 6.350 , &
627 6.400 , 6.450 , 6.500 , 6.550 , &
628 6.600 , 6.650 , 6.700 , 6.750 , &
629 6.800 , 6.850 , 6.900 , 6.950 , &
630 7.000 , 7.050 , 7.100 , 7.150 , &
631 7.200 , 7.250 , 7.300 , 7.350 , &
632 7.400 , 7.450 , 7.500 , 7.550 , &
633 7.600 , 7.650 , 7.700 , 7.750 , &
634 7.800 , 7.850 , 7.900 , 7.950 , &
635 8.000 , 8.100 , 8.200 , 8.300 , &
636 8.400 , 8.500 , 8.600 , 8.700 , &
637 8.800 , 8.900 , 9.00 /
638
639 data l_cl_ism / &
640 -28.365 , -28.242 , -28.127 , -28.020 , &
641 -27.919 , -27.825 , -27.736 , -27.653 , &
642 -27.575 , -27.504 , -27.437 , -27.376 , &
643 -27.319 , -27.267 , -27.220 , -27.176 , &
644 -27.134 , -27.095 , -27.058 , -27.021 , &
645 -26.985 , -26.948 , -26.910 , -26.870 , &
646 -26.827 , -26.775 , -26.721 , -26.664 , &
647 -26.608 , -26.552 , -26.495 , -26.437 , &
648 -26.378 , -26.317 , -26.255 , -26.190 , &
649 -26.123 , -26.053 , -25.984 , -25.913 , &
650 -25.847 , -25.786 , -25.736 , -25.702 , &
651 -25.678 , -25.662 , -25.649 , -25.636 , &
652 -25.621 , -25.604 , -25.587 , -25.571 , &
653 -25.562 , -25.526 , -25.505 , -25.499 , &
654 -25.499 , -25.491 , -25.468 , -25.410 , &
655 -25.268 , -24.888 , -23.702 , -22.624 , &
656 -22.036 , -21.843 , -21.854 , -21.924 , &
657 -21.986 , -22.017 , -22.021 , -22.005 , &
658 -21.964 , -21.896 , -21.806 , -21.699 , &
659 -21.580 , -21.463 , -21.370 , -21.312 , &
660 -21.284 , -21.290 , -21.322 , -21.345 , &
661 -21.354 , -21.366 , -21.385 , -21.396 , &
662 -21.414 , -21.483 , -21.600 , -21.696 , &
663 -21.742 , -21.759 , -21.776 , -21.816 , &
664 -21.885 , -21.939 , -21.946 , -21.918 , &
665 -21.873 , -21.818 , -21.756 , -21.689 , &
666 -21.618 , -21.547 , -21.475 , -21.403 , &
667 -21.331 , -21.260 , -21.188 , -21.114 , &
668 -21.039 , -20.963 , -20.887 , -20.810 , &
669 -20.734 , -20.657 , -20.581 , -20.505 , &
670 -20.429 , -20.352 , -20.276 , -20.200 , &
671 -20.125 , -20.049 , -19.973 , -19.898 , &
672 -19.822 , -19.747 , -19.671 , -19.596 , &
673 -19.520 , -19.445 , -19.370 , -19.295 , &
674 -19.220 , -19.144 , -19.069 , -18.994 , &
675 -18.919 , -18.869 , -18.819 , -18.769 , &
676 -18.719 , -18.669 , -18.619 , -18.569 , &
677 -18.519 , -18.469 , -18.419 /
678
679 data n_dere / 101 /
680
681 data t_dere / 4.00, 4.05, 4.10, 4.15, 4.20, 4.25, 4.30, 4.35, &
682 4.40, 4.45, 4.50, 4.55, 4.60, 4.65, 4.70, 4.75, &
683 4.80, 4.85, 4.90, 4.95, 5.00, 5.05, 5.10, 5.15, &
684 5.20, 5.25, 5.30, 5.35, 5.40, 5.45, 5.50, 5.55, &
685 5.60, 5.65, 5.70, 5.75, 5.80, 5.85, 5.90, 5.95, &
686 6.00, 6.05, 6.10, 6.15, 6.20, 6.25, 6.30, 6.35, &
687 6.40, 6.45, 6.50, 6.55, 6.60, 6.65, 6.70, 6.75, &
688 6.80, 6.85, 6.90, 6.95, 7.00, 7.05, 7.10, 7.15, &
689 7.20, 7.25, 7.30, 7.35, 7.40, 7.45, 7.50, 7.55, &
690 7.60, 7.65, 7.70, 7.75, 7.80, 7.85, 7.90, 7.95, &
691 8.00, 8.05, 8.10, 8.15, 8.20, 8.25, 8.30, 8.35, &
692 8.40, 8.45, 8.50, 8.55, 8.60, 8.65, 8.70, 8.75, &
693 8.80, 8.85, 8.90, 8.95, 9.00 /
694
695 data l_dere_corona / &
696 -23.00744648, -22.55439580, -22.15614458, -21.83268267, -21.64589156, &
697 -21.61618463, -21.68402965, -21.79048499, -21.87614836, -21.91009489, &
698 -21.89962945, -21.86012091, -21.79588002, -21.71669877, -21.62342304, &
699 -21.52143350, -21.41793664, -21.33068312, -21.27736608, -21.25181197, &
700 -21.24184538, -21.25806092, -21.27901426, -21.27164622, -21.24412514, &
701 -21.21467016, -21.19586057, -21.18309616, -21.18708664, -21.24108811, &
702 -21.35163999, -21.45099674, -21.48678240, -21.47625353, -21.45222529, &
703 -21.43297363, -21.42596873, -21.42021640, -21.41005040, -21.40120949, &
704 -21.40450378, -21.42250820, -21.44977165, -21.47755577, -21.51144928, &
705 -21.55909092, -21.63451202, -21.73754891, -21.85078089, -21.95467702, &
706 -22.03526908, -22.08990945, -22.11804503, -22.12436006, -22.11633856, &
707 -22.10072681, -22.08301995, -22.06701918, -22.05650548, -22.05551733, &
708 -22.06803389, -22.10072681, -22.15926677, -22.24033216, -22.32882716, &
709 -22.40782324, -22.47108330, -22.51855737, -22.54975089, -22.57024772, &
710 -22.58004425, -22.58335949, -22.58169871, -22.57511836, -22.56543110, &
711 -22.55439580, -22.54060751, -22.52724355, -22.51427857, -22.49894074, &
712 -22.48545225, -22.47108330, -22.45593196, -22.44009337, -22.42365865, &
713 -22.40671393, -22.38933984, -22.37059040, -22.35163999, -22.33161408, &
714 -22.31158018, -22.29073004, -22.26921772, -22.24795155, -22.22621356, &
715 -22.20411998, -22.18111459, -22.15864053, -22.13608262, -22.11294562, &
716 -22.08937560 /
717
718 data l_dere_photo / &
719 -23.25649024, -22.74232143, -22.28988263, -21.93554201, -21.74232143, &
720 -21.72353820, -21.80410035, -21.91009489, -22.00000000, -22.04143612, &
721 -22.04575749, -22.02502801, -21.97469413, -21.89962945, -21.80410035, &
722 -21.69250396, -21.57511836, -21.46852108, -21.38827669, -21.33629907, &
723 -21.31069114, -21.31966449, -21.33724217, -21.32882716, -21.30189945, &
724 -21.27408837, -21.25649024, -21.24718357, -21.25649024, -21.32148162, &
725 -21.46092390, -21.61083392, -21.70774393, -21.74958000, -21.76955108, &
726 -21.79588002, -21.84466396, -21.88941029, -21.91009489, -21.91721463, &
727 -21.92811799, -21.94692156, -21.96657624, -21.99139983, -22.01412464, &
728 -22.04963515, -22.10402527, -22.17848647, -22.26201267, -22.34390180, &
729 -22.41680123, -22.47237010, -22.50723961, -22.52578374, -22.53313238, &
730 -22.53165267, -22.52578374, -22.51855737, -22.51286162, -22.51286162, &
731 -22.52143350, -22.54211810, -22.57839607, -22.62708800, -22.67366414, &
732 -22.70996539, -22.73282827, -22.74714697, -22.75202673, -22.74958000, &
733 -22.74232143, -22.73282827, -22.72124640, -22.70553377, -22.69036983, &
734 -22.67162040, -22.65364703, -22.63638802, -22.61618463, -22.59687948, &
735 -22.57675413, -22.55752023, -22.53610701, -22.51570016, -22.49485002, &
736 -22.47366072, -22.45222529, -22.42945706, -22.40782324, -22.38510278, &
737 -22.36251027, -22.33913452, -22.31605287, -22.29242982, -22.26921772, &
738 -22.24565166, -22.22184875, -22.19859629, -22.17457388, -22.15058059, &
739 -22.12609840 /
740
741 data n_colgan / 55 /
742
743 data t_colgan / 4.06460772, 4.14229559, 4.21995109, 4.29760733, 4.37527944, 4.45293587, &
744 4.53060946, 4.60826923, 4.68592974, 4.76359269, 4.79704583, 4.83049243, &
745 4.86394114, 4.89738514, 4.93083701, 4.96428321, 4.99773141, 5.03116600, &
746 5.06460772, 5.17574368, 5.28683805, 5.39795738, 5.50906805, 5.62017771, &
747 5.73129054, 5.84240328, 5.95351325, 6.06460772, 6.17574368, 6.28683805, &
748 6.39795738, 6.50906805, 6.62017771, 6.73129054, 6.84240328, 6.95351325, &
749 7.06460772, 7.17574368, 7.28683805, 7.39795738, 7.50906805, 7.62017771, &
750 7.73129054, 7.84240328, 7.95351325, 8.06460772, 8.17574368, 8.28683805, &
751 8.39795738, 8.50906805, 8.62017771, 8.73129054, 8.84240328, 8.95351325, &
752 9.06460772 /
753
754 data l_colgan / -22.18883401, -21.78629635, -21.60383554, -21.68480662, -21.76444630, &
755 -21.67935529, -21.54217864, -21.37958284, -21.25171892, -21.17584161, &
756 -21.15783402, -21.14491111, -21.13526945, -21.12837453, -21.12485189, &
757 -21.12438898, -21.12641785, -21.12802448, -21.12547760, -21.08964778, &
758 -21.08812360, -21.19542445, -21.34582346, -21.34839251, -21.31700703, &
759 -21.29072156, -21.28900309, -21.34104468, -21.43122351, -21.62448270, &
760 -21.86694036, -22.02897478, -22.08050874, -22.06057061, -22.01973295, &
761 -22.00000434, -22.05161149, -22.22175466, -22.41451671, -22.52581288, &
762 -22.56913516, -22.57485721, -22.56150512, -22.53968863, -22.51490350, &
763 -22.48895932, -22.46071057, -22.42908363, -22.39358639, -22.35456791, &
764 -22.31261375, -22.26827428, -22.22203698, -22.17422996, -22.12514145 /
765
766 contains
767
768 !> Radiative cooling initialization
769 subroutine radiative_cooling_init_params(phys_gamma,He_abund)
770 use mod_global_parameters
771 double precision, intent(in) :: phys_gamma,He_abund
772
773 rc_gamma=phys_gamma
774 he_abundance=he_abund
775 end subroutine radiative_cooling_init_params
776
777 subroutine radiative_cooling_init(fl,read_params)
778 use mod_global_parameters
779 interface
780 subroutine read_params(fl)
781 use mod_global_parameters, only: unitpar,par_files
782 import rc_fluid
783 type(rc_fluid), intent(inout) :: fl
784
785 end subroutine read_params
786 end interface
787
788 type(rc_fluid), intent(inout) :: fl
789
790 double precision, dimension(:), allocatable :: t_table
791 double precision, dimension(:), allocatable :: L_table
792 double precision :: ratt, fact1, fact2, fact3, dL1, dL2
793 double precision :: tstep, Lstep
794 integer :: ntable, i, j
795 logical :: jump
796 Character(len=65) :: PPL_curves(1:6)
797
798 fl%ncool=4000
799 fl%coolcurve='JCcorona'
800 fl%coolmethod='exact'
801 fl%tlow=bigdouble
802 fl%Tfix=.false.
803 fl%rc_split=.false.
804 fl%rad_cut=.false.
805 fl%rad_cut_hgt=0.5d0
806 fl%rad_cut_dey=0.15d0
807 call read_params(fl)
808
809 if(fl%rc_split) any_source_split=.true.
810
811 ! Checks if coolcurve is a piecewise power law (PPL)
812 ppl_curves = [Character(len=65) :: 'Hildner','FM', 'Rosner', 'Klimchuk','SPEX_DM_rough','SPEX_DM_fine']
813 do i=1,size(ppl_curves)
814 if (ppl_curves(i)==fl%coolcurve) then
815 fl%isPPL = .true.
816 end if
817 end do
818
819 ! Init for PPL
820 if (fl%isPPL) then
821 ! Read in tables and create t_PPL, l_PPL, a_PPL
822 select case(fl%coolcurve)
823
824 case('Hildner')
825 if(mype ==0) &
826 print *,'Use Hildner (1974) piecewise power law'
827 fl%n_PPL = n_hildner
828 allocate(fl%t_PPL(1:fl%n_PPL+1), fl%l_PPL(1:fl%n_PPL+1))
829 allocate(fl%a_PPL(1:fl%n_PPL))
830 fl%t_PPL(1:fl%n_PPL+1) = t_hildner(1:n_hildner+1)
831 fl%a_PPL(1:fl%n_PPL) = a_hildner(1:n_hildner)
832 fl%l_PPL(1:fl%n_PPL) = 10.d0**x_hildner(1:n_hildner) * (10.d0**fl%t_PPL(1:fl%n_PPL))**fl%a_PPL(1:fl%n_PPL)
833
834 case('FM')
835 if(mype==0) &
836 print *,'Use Forbes and Malherbe (1991)-like piecewise power law'
837 fl%n_PPL = n_fm
838 allocate(fl%t_PPL(1:fl%n_PPL+1), fl%l_PPL(1:fl%n_PPL+1))
839 allocate(fl%a_PPL(1:fl%n_PPL))
840 fl%t_PPL(1:fl%n_PPL+1) = t_fm(1:n_fm+1)
841 fl%a_PPL(1:fl%n_PPL) = a_fm(1:n_fm)
842 fl%l_PPL(1:fl%n_PPL) = 10.d0**x_fm(1:n_fm) * (10.d0**fl%t_PPL(1:fl%n_PPL))**fl%a_PPL(1:fl%n_PPL)
843
844 case('Rosner')
845 if(mype==0) &
846 print *,'Use piecewise power law according to Rosner (1978)'
847 if(mype ==0) &
848 print *,'and extended by Priest (1982) from Van Der Linden (1991)'
849 fl%n_PPL = n_rosner
850 allocate(fl%t_PPL(1:fl%n_PPL+1), fl%l_PPL(1:fl%n_PPL+1))
851 allocate(fl%a_PPL(1:fl%n_PPL))
852 fl%t_PPL(1:fl%n_PPL+1) = t_rosner(1:n_rosner+1)
853 fl%a_PPL(1:fl%n_PPL) = a_rosner(1:n_rosner)
854 fl%l_PPL(1:fl%n_PPL) = 10.d0**x_rosner(1:n_rosner) * (10.d0**fl%t_PPL(1:fl%n_PPL))**fl%a_PPL(1:fl%n_PPL)
855
856 case('Klimchuk')
857 if(mype==0) &
858 print *,'Use Klimchuk (2008) piecewise power law'
859 fl%n_PPL = n_klimchuk
860 allocate(fl%t_PPL(1:fl%n_PPL+1), fl%l_PPL(1:fl%n_PPL+1))
861 allocate(fl%a_PPL(1:fl%n_PPL))
862 fl%t_PPL(1:fl%n_PPL+1) = t_klimchuk(1:n_klimchuk+1)
863 fl%a_PPL(1:fl%n_PPL) = a_klimchuk(1:n_klimchuk)
864 fl%l_PPL(1:fl%n_PPL) = 10.d0**x_klimchuk(1:n_klimchuk) * (10.d0**fl%t_PPL(1:fl%n_PPL))**fl%a_PPL(1:fl%n_PPL)
865
866 case('SPEX_DM_rough')
867 if(mype==0) &
868 print *,'Use the rough piece wise power law fit to the SPEX_DM curve (2009)'
869 fl%n_PPL = n_spex_dm_rough
870 allocate(fl%t_PPL(1:fl%n_PPL+1), fl%l_PPL(1:fl%n_PPL+1))
871 allocate(fl%a_PPL(1:fl%n_PPL))
872 fl%t_PPL(1:fl%n_PPL+1) = t_spex_dm_rough(1:n_spex_dm_rough+1)
873 fl%a_PPL(1:fl%n_PPL) = a_spex_dm_rough(1:n_spex_dm_rough)
874 fl%l_PPL(1:fl%n_PPL) = 10.d0**x_spex_dm_rough(1:n_spex_dm_rough) * (10.d0**fl%t_PPL(1:fl%n_PPL))**fl%a_PPL(1:fl%n_PPL)
875
876 case('SPEX_DM_fine')
877 if(mype==0) &
878 print *,'Use the fine, detailed piece wise power law fit to the SPEX_DM curve (2009)'
879 fl%n_PPL = n_spex_dm_fine
880 allocate(fl%t_PPL(1:fl%n_PPL+1), fl%l_PPL(1:fl%n_PPL+1))
881 allocate(fl%a_PPL(1:fl%n_PPL))
882 fl%t_PPL(1:fl%n_PPL+1) = t_spex_dm_fine(1:n_spex_dm_fine+1)
883 fl%a_PPL(1:fl%n_PPL) = a_spex_dm_fine(1:n_spex_dm_fine)
884 fl%l_PPL(1:fl%n_PPL) = 10.d0**x_spex_dm_fine(1:n_spex_dm_fine) * (10.d0**fl%t_PPL(1:fl%n_PPL))**fl%a_PPL(1:fl%n_PPL)
885
886 case default
887 call mpistop("This piecewise power law is unknown")
888 end select
889
890 ! Go from logarithmic to actual values.
891 fl%t_PPL(1:fl%n_PPL+1) = 10.d0**fl%t_PPL(1:fl%n_PPL+1)
892 ! Change unit of table if SI is used instead of cgs
893 if (si_unit) fl%l_PPL(1:fl%n_PPL) = fl%l_PPL(1:fl%n_PPL) * 10.0d0**(-13)
894
895 ! Make dimensionless
896 fl%t_PPL(1:fl%n_PPL+1) = fl%t_PPL(1:fl%n_PPL+1) / unit_temperature
897 fl%l_PPL(1:fl%n_PPL) = fl%l_PPL(1:fl%n_PPL) * unit_numberdensity**2 * unit_time / unit_pressure * (1.d0+2.d0*he_abundance)
898
899 ! Set tref en lref
900 fl%l_PPL(fl%n_PPL+1) = fl%l_PPL(fl%n_PPL) * ( fl%t_PPL(fl%n_PPL+1) / fl%t_PPL(fl%n_PPL) )**fl%a_PPL(fl%n_PPL)
901 fl%lref = fl%l_PPL(fl%n_PPL+1)
902 fl%tref = fl%t_PPL(fl%n_PPL+1)
903
904 ! Set tcoolmin and tcoolmax
905 fl%tcoolmin = fl%t_PPL(1)
906 fl%tcoolmax = fl%t_PPL(fl%n_PPL+1)
907 ! smaller value for lowest temperatures from cooling table and user's choice
908 if (fl%tlow==bigdouble) fl%tlow=fl%tcoolmin
909 !create y_PPL
910 call create_y_ppl(fl)
911
912 else
913
914 ! Init for interpolatable tables
915 allocate(fl%tcool(1:fl%ncool), fl%Lcool(1:fl%ncool), fl%dLdtcool(1:fl%ncool))
916 allocate(fl%Yc(1:fl%ncool), fl%invYc(1:fl%ncool))
917
918 fl%tcool(1:fl%ncool) = zero
919 fl%Lcool(1:fl%ncool) = zero
920 fl%dLdtcool(1:fl%ncool) = zero
921
922 ! Read in the selected cooling curve
923 select case(fl%coolcurve)
924
925 case('JCcorona')
926 if(mype ==0) &
927 print *,'Use Colgan & Feldman (2008) cooling curve'
928 if(mype ==0) &
929 print *,'This version only till 10000 K, beware for floor T treatment'
930 ntable = n_jccorona
931 allocate(t_table(1:ntable))
932 allocate(l_table(1:ntable))
933 t_table(1:ntable) = t_jccorona(1:n_jccorona)
934 l_table(1:ntable) = l_jccorona(1:n_jccorona)
935
936 case('DM')
937 if(mype ==0) &
938 print *,'Use Dalgarno & McCray (1972) cooling curve'
939 ntable = n_dm
940 allocate(t_table(1:ntable))
941 allocate(l_table(1:ntable))
942 t_table(1:ntable) = t_dm(1:n_dm)
943 l_table(1:ntable) = l_dm(1:n_dm)
944
945 case('MB')
946 if(mype ==0) &
947 write(*,'(3a)') 'Use MacDonald & Bailey (1981) cooling curve '&
948 ,'as implemented in ZEUS-3D, with the values '&
949 ,'from Dalgarno & McCRay (1972) for low temperatures.'
950 ntable = n_mb + 20
951 allocate(t_table(1:ntable))
952 allocate(l_table(1:ntable))
953 t_table(1:ntable) = t_dm(1:21)
954 l_table(1:ntable) = l_dm(1:21)
955 t_table(22:ntable) = t_mb(2:n_mb)
956 l_table(22:ntable) = l_mb(2:n_mb)
957
958 case('MLcosmol')
959 if(mype ==0) &
960 print *,'Use Mellema & Lundqvist (2002) cooling curve '&
961 ,'for zero metallicity '
962 ntable = n_mlcosmol
963 allocate(t_table(1:ntable))
964 allocate(l_table(1:ntable))
965 t_table(1:ntable) = t_mlcosmol(1:n_mlcosmol)
966 l_table(1:ntable) = l_mlcosmol(1:n_mlcosmol)
967
968 case('MLwc')
969 if(mype ==0) &
970 print *,'Use Mellema & Lundqvist (2002) cooling curve '&
971 ,'for WC-star metallicity '
972 ntable = n_mlwc
973 allocate(t_table(1:ntable))
974 allocate(l_table(1:ntable))
975 t_table(1:ntable) = t_mlwc(1:n_mlwc)
976 l_table(1:ntable) = l_mlwc(1:n_mlwc)
977
978 case('MLsolar1')
979 if(mype ==0) &
980 print *,'Use Mellema & Lundqvist (2002) cooling curve '&
981 ,'for solar metallicity '
982 ntable = n_mlsolar1
983 allocate(t_table(1:ntable))
984 allocate(l_table(1:ntable))
985 t_table(1:ntable) = t_mlsolar1(1:n_mlsolar1)
986 l_table(1:ntable) = l_mlsolar1(1:n_mlsolar1)
987
988 case('cloudy_ism')
989 if(mype ==0) &
990 print *,'Use Cloudy based cooling curve '&
991 ,'for ism metallicity '
992 ntable = n_cl_ism
993 allocate(t_table(1:ntable))
994 allocate(l_table(1:ntable))
995 t_table(1:ntable) = t_cl_ism(1:n_cl_ism)
996 l_table(1:ntable) = l_cl_ism(1:n_cl_ism)
997
998 case('cloudy_solar')
999 if(mype ==0) &
1000 print *,'Use Cloudy based cooling curve '&
1001 ,'for solar metallicity '
1002 ntable = n_cl_solar
1003 allocate(t_table(1:ntable))
1004 allocate(l_table(1:ntable))
1005 t_table(1:ntable) = t_cl_solar(1:n_cl_solar)
1006 l_table(1:ntable) = l_cl_solar(1:n_cl_solar)
1007
1008 case('SPEX')
1009 if(mype ==0) &
1010 print *,'Use SPEX cooling curve (Schure et al. 2009) '&
1011 ,'for solar metallicity '
1012 ntable = n_spex
1013 allocate(t_table(1:ntable))
1014 allocate(l_table(1:ntable))
1015 t_table(1:ntable) = t_spex(1:n_spex)
1016 l_table(1:ntable) = l_spex(1:n_spex) + log10(nenh_spex(1:n_spex))
1017
1018 case('SPEX_DM')
1019 if(mype ==0) then
1020 print *, 'Use SPEX cooling curve for solar metallicity above 10^4 K. '
1021 print *, 'At lower temperatures,use Dalgarno & McCray (1972), '
1022 print *, 'with a pre-set ionization fraction of 10^-3. '
1023 print *, 'as described by Schure et al. (2009). '
1024 endif
1025 ntable = n_spex + n_dm_2 - 6
1026 allocate(t_table(1:ntable))
1027 allocate(l_table(1:ntable))
1028 t_table(1:n_dm_2-1) = t_dm_2(1:n_dm_2-1)
1029 l_table(1:n_dm_2-1) = l_dm_2(1:n_dm_2-1)
1030 t_table(n_dm_2:ntable) = t_spex(6:n_spex)
1031 l_table(n_dm_2:ntable) = l_spex(6:n_spex) + log10(nenh_spex(6:n_spex))
1032
1033 case('Dere_corona')
1034 if(mype ==0) &
1035 print *,'Use Dere (2009) cooling curve for solar corona'
1036 ntable = n_dere
1037 allocate(t_table(1:ntable))
1038 allocate(l_table(1:ntable))
1039 t_table(1:ntable) = t_dere(1:n_dere)
1040 l_table(1:ntable) = l_dere_corona(1:n_dere)
1041
1042 case('Dere_corona_DM')
1043 if(mype==0)&
1044 print *, 'Combination of Dere_corona (2009) for high temperatures and'
1045 if(mype==0)&
1046 print *, 'Dalgarno & McCray (1972), DM2, for low temperatures'
1047 ntable = n_dere + n_dm_2 - 1
1048 allocate(t_table(1:ntable))
1049 allocate(l_table(1:ntable))
1050 t_table(1:n_dm_2-1) = t_dm_2(1:n_dm_2-1)
1051 l_table(1:n_dm_2-1) = l_dm_2(1:n_dm_2-1)
1052 t_table(n_dm_2:ntable) = t_dere(1:n_dere)
1053 l_table(n_dm_2:ntable) = l_dere_corona(1:n_dere)
1054
1055 case('Dere_photo')
1056 if(mype ==0) &
1057 print *,'Use Dere (2009) cooling curve for solar photophere'
1058 ntable = n_dere
1059 allocate(t_table(1:ntable))
1060 allocate(l_table(1:ntable))
1061 t_table(1:ntable) = t_dere(1:n_dere)
1062 l_table(1:ntable) = l_dere_photo(1:n_dere)
1063
1064 case('Dere_photo_DM')
1065 if(mype==0)&
1066 print *, 'Combination of Dere_photo (2009) for high temperatures and'
1067 if(mype==0)&
1068 print *, 'Dalgarno & McCray (1972), DM2, for low temperatures'
1069 ntable = n_dere + n_dm_2 - 1
1070 allocate(t_table(1:ntable))
1071 allocate(l_table(1:ntable))
1072 t_table(1:n_dm_2-1) = t_dm_2(1:n_dm_2-1)
1073 l_table(1:n_dm_2-1) = l_dm_2(1:n_dm_2-1)
1074 t_table(n_dm_2:ntable) = t_dere(1:n_dere)
1075 l_table(n_dm_2:ntable) = l_dere_photo(1:n_dere)
1076
1077 case('Colgan')
1078 if(mype==0) &
1079 print *, 'Use Colgan (2008) cooling curve'
1080 ntable = n_colgan
1081 allocate(t_table(1:ntable))
1082 allocate(l_table(1:ntable))
1083 t_table(1:ntable) = t_colgan(1:n_colgan)
1084 l_table(1:ntable) = l_colgan(1:n_colgan)
1085
1086 case('Colgan_DM')
1087 if(mype==0)&
1088 print *, 'Combination of Colgan (2008) for high temperatures and'
1089 if(mype==0)&
1090 print *, 'Dalgarno & McCray (1972), DM2, for low temperatures'
1091 ntable = n_colgan + n_dm_2
1092 allocate(t_table(1:ntable))
1093 allocate(l_table(1:ntable))
1094 t_table(1:n_dm_2) = t_dm_2(1:n_dm_2)
1095 l_table(1:n_dm_2) = l_dm_2(1:n_dm_2)
1096 t_table(n_dm_2+1:ntable) = t_colgan(1:n_colgan)
1097 l_table(n_dm_2+1:ntable) = l_colgan(1:n_colgan)
1098
1099 case default
1100 call mpistop("This coolingcurve is unknown")
1101 end select
1102
1103 ! create cooling table(s) for use in amrvac
1104 fl%tcoolmax = t_table(ntable)
1105 fl%tcoolmin = t_table(1)
1106 ratt = (fl%tcoolmax-fl%tcoolmin)/( dble(fl%ncool-1) + smalldouble)
1107
1108 fl%tcool(1) = fl%tcoolmin
1109 fl%Lcool(1) = l_table(1)
1110
1111 fl%tcool(fl%ncool) = fl%tcoolmax
1112 fl%Lcool(fl%ncool) = l_table(ntable)
1113
1114 do i=2,fl%ncool ! loop to create one table
1115 fl%tcool(i) = fl%tcool(i-1)+ratt
1116 do j=1,ntable-1 ! loop to create one spot on a table
1117 ! Second order polynomial interpolation, except at the outer edge,
1118 ! or in case of a large jump.
1119 if(fl%tcool(i) < t_table(j+1)) then
1120 if(j.eq. ntable-1 )then
1121 fact1 = (fl%tcool(i)-t_table(j+1)) &
1122 /(t_table(j)-t_table(j+1))
1123 fact2 = (fl%tcool(i)-t_table(j)) &
1124 /(t_table(j+1)-t_table(j))
1125 fl%Lcool(i) = l_table(j)*fact1 + l_table(j+1)*fact2
1126 exit
1127 else
1128 dl1 = l_table(j+1)-l_table(j)
1129 dl2 = l_table(j+2)-l_table(j+1)
1130 jump =(max(dabs(dl1),dabs(dl2)) > 2*min(dabs(dl1),dabs(dl2)))
1131 end if
1132 if( jump ) then
1133 fact1 = (fl%tcool(i)-t_table(j+1)) &
1134 /(t_table(j)-t_table(j+1))
1135 fact2 = (fl%tcool(i)-t_table(j)) &
1136 /(t_table(j+1)-t_table(j))
1137 fl%Lcool(i) = l_table(j)*fact1 + l_table(j+1)*fact2
1138 exit
1139 else
1140 fact1 = ((fl%tcool(i)-t_table(j+1)) &
1141 * (fl%tcool(i)-t_table(j+2))) &
1142 / ((t_table(j)-t_table(j+1)) &
1143 * (t_table(j)-t_table(j+2)))
1144 fact2 = ((fl%tcool(i)-t_table(j)) &
1145 * (fl%tcool(i)-t_table(j+2))) &
1146 / ((t_table(j+1)-t_table(j)) &
1147 * (t_table(j+1)-t_table(j+2)))
1148 fact3 = ((fl%tcool(i)-t_table(j)) &
1149 * (fl%tcool(i)-t_table(j+1))) &
1150 / ((t_table(j+2)-t_table(j)) &
1151 * (t_table(j+2)-t_table(j+1)))
1152 fl%Lcool(i) = l_table(j)*fact1 + l_table(j+1)*fact2 &
1153 + l_table(j+2)*fact3
1154 exit
1155 end if
1156 end if
1157 end do ! end loop to find create one spot on a table
1158 end do ! end loop to create one table
1159
1160 ! Go from logarithmic to actual values.
1161 fl%tcool(1:fl%ncool) = 10.0d0**fl%tcool(1:fl%ncool)
1162 fl%Lcool(1:fl%ncool) = 10.0d0**fl%Lcool(1:fl%ncool)
1163
1164 ! Change unit of table if SI is used instead of cgs
1165 if (si_unit) fl%Lcool(1:fl%ncool) = fl%Lcool(1:fl%ncool) * 10.0d0**(-13)
1166
1167 ! Scale both T and Lambda
1168 fl%tcool(1:fl%ncool) = fl%tcool(1:fl%ncool) / unit_temperature
1169 fl%Lcool(1:fl%ncool) = fl%Lcool(1:fl%ncool) * unit_numberdensity**2 * unit_time / unit_pressure * (1.d0+2.d0*he_abundance)
1170
1171 fl%tcoolmin = fl%tcool(1)+smalldouble ! avoid pointless interpolation
1172 ! smaller value for lowest temperatures from cooling table and user's choice
1173 if (fl%tlow==bigdouble) fl%tlow=fl%tcoolmin
1174 fl%tcoolmax = fl%tcool(fl%ncool)
1175 fl%lgtcoolmin = dlog10(fl%tcoolmin)
1176 fl%lgtcoolmax = dlog10(fl%tcoolmax)
1177 fl%lgstep = (fl%lgtcoolmax-fl%lgtcoolmin) * 1.d0 / (fl%ncool-1)
1178 fl%dLdtcool(1) = (fl%Lcool(2)-fl%Lcool(1))/(fl%tcool(2)-fl%tcool(1))
1179 fl%dLdtcool(fl%ncool) = (fl%Lcool(fl%ncool)-fl%Lcool(fl%ncool-1))/(fl%tcool(fl%ncool)-fl%tcool(fl%ncool-1))
1180
1181 do i=2,fl%ncool-1
1182 fl%dLdtcool(i) = (fl%Lcool(i+1)-fl%Lcool(i-1))/(fl%tcool(i+1)-fl%tcool(i-1))
1183 end do
1184
1185 deallocate(t_table)
1186 deallocate(l_table)
1187
1188 if( fl%coolmethod == 'exact' ) then
1189 fl%tref = fl%tcoolmax
1190 fl%lref = fl%Lcool(fl%ncool)
1191 fl%Yc(fl%ncool) = zero
1192 do i=fl%ncool-1, 1, -1
1193 fl%Yc(i) = fl%Yc(i+1)
1194 do j=1,100
1195 tstep = 1.0d-2*(fl%tcool(i+1)-fl%tcool(i))
1196 call findl(fl%tcool(i+1)-j*tstep, lstep, fl)
1197 fl%Yc(i) = fl%Yc(i) + fl%lref/fl%tref*tstep/lstep
1198 end do
1199 end do
1200 end if
1201 end if
1202
1203 rc_gamma_1=rc_gamma-1.d0
1204 invgam = 1.d0/rc_gamma_1
1205
1206 end subroutine radiative_cooling_init
1207
1208 subroutine create_y_ppl(fl)
1209 ! creates the constants of integration needed for solving
1210 ! the cooling law exact for a piecewise power law
1211 ! In correspondence with eq. A6 of Townsend (2009)
1212 use mod_global_parameters
1213 type(rc_fluid) :: fl
1214 double precision :: y_extra, factor
1215 integer :: i
1216
1217 allocate(fl%y_PPL(1:fl%n_PPL+1))
1218
1219 fl%y_PPL(1:fl%n_PPL+1) = zero
1220
1221 do i=fl%n_PPL, 1, -1
1222 factor = fl%l_PPL(fl%n_PPL+1) * fl%t_PPL(i) / (fl%l_PPL(i) * fl%t_PPL(fl%n_PPL+1))
1223 if (fl%a_PPL(i) == 1.d0) then
1224 y_extra = log( fl%t_PPL(i) / fl%t_PPL(i+1) )
1225 else
1226 y_extra = 1 / (1 - fl%a_PPL(i)) * (1 - ( fl%t_PPL(i) / fl%t_PPL(i+1) )**(fl%a_PPL(i)-1) )
1227 end if
1228 fl%y_PPL(i) = fl%y_PPL(i+1) - factor*y_extra
1229 end do
1230 end subroutine create_y_ppl
1231
1232 subroutine cooling_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x,fl)
1233 use mod_global_parameters
1234
1235 integer, intent(in) :: ixI^L, ixO^L
1236 double precision, intent(in) :: dx^D, x(ixI^S,1:ndim), w(ixI^S,1:nw)
1237 type(rc_fluid), intent(in) :: fl
1238 double precision, intent(inout) :: dtnew
1239
1240 double precision :: etherm(ixI^S), rho(ixI^S), Rfactor(ixI^S)
1241 double precision :: L1,Te(ixI^S), pth(ixI^S), lum(ixI^S)
1242 integer :: ix^D
1243 !
1244 ! Limit timestep to avoid cooling problems when using explicit cooling
1245 !
1246 if(fl%coolmethod == 'explicit1') then
1247 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
1248 call fl%get_rho(w,x,ixi^l,ixo^l,rho)
1249 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,rfactor)
1250 te(ixo^s)=pth(ixo^s)/(rho(ixo^s)*rfactor(ixo^s))
1251 {do ix^db = ixo^lim^db\}
1252 ! Determine explicit cooling
1253 ! If temperature is below floor level, no cooling.
1254 ! Stop wasting time and go to next gridpoint.
1255 ! If the temperature is higher than the maximum,
1256 ! assume Bremsstrahlung
1257 if( te(ix^d)<=fl%tcoolmin ) then
1258 l1 = zero
1259 else if( te(ix^d)>=fl%tcoolmax )then
1260 call calc_l_extended(te(ix^d), l1, fl)
1261 l1 = l1*rho(ix^d)**2
1262 else
1263 call findl(te(ix^d),l1,fl)
1264 l1 = l1*rho(ix^d)**2
1265 end if
1266 lum(ix^d) = l1
1267 {end do\}
1268 etherm(ixo^s)=pth(ixo^s)*invgam
1269 dtnew =fl%cfrac*minval(etherm(ixo^s)/max(lum(ixo^s),smalldouble))
1270 end if
1271 end subroutine cooling_get_dt
1272
1273 subroutine getvar_cooling(ixI^L,ixO^L,w,x,coolrate,fl)
1274 ! Create extra variable to show cooling rate in the output
1275 ! Uses a simple explicit scheme.
1276 ! N.B. Since there is no knowledge of the timestep size,
1277 ! there is no upper limit for the cooling rate.
1278 use mod_global_parameters
1279
1280 integer, intent(in) :: ixI^L,ixO^L
1281 double precision, intent(in) :: x(ixI^S,1:ndim)
1282 double precision :: w(ixI^S,1:nw)
1283 double precision, intent(out):: coolrate(ixI^S)
1284 type(rc_fluid), intent(in) :: fl
1285
1286 double precision :: pth(ixI^S),rho(ixI^S)
1287 double precision :: L1,Te(ixI^S),Rfactor(ixI^S)
1288 integer :: ix^D
1289
1290 call fl%get_pthermal(w,x,ixi^l,ixo^l,pth)
1291 call fl%get_rho(w,x,ixi^l,ixo^l,rho)
1292 call fl%get_var_Rfactor(w,x,ixi^l,ixo^l,rfactor)
1293 te(ixo^s) = pth(ixo^s) / (rho(ixo^s)*rfactor(ixo^s))
1294
1295 {do ix^db = ixo^lim^db\}
1296 ! Determine explicit cooling
1297 if(te(ix^d) <= fl%tcoolmin) then
1298 l1 = zero
1299 else if(te(ix^d) >= fl%tcoolmax)then
1300 call calc_l_extended(te(ix^d),l1,fl)
1301 l1 = l1*rho(ix^d)**2
1302 else
1303 call findl(te(ix^d),l1,fl)
1304 l1 = l1*rho(ix^d)**2
1305 end if
1306 if(slab_uniform .and. fl%rad_cut .and. x(ix^d,ndim) .le. fl%rad_cut_hgt) then
1307 l1 = l1*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1308 end if
1309 coolrate(ix^d) = l1
1310 {end do\}
1311 end subroutine getvar_cooling
1312
1313 subroutine getvar_cooling_exact(qdt, ixI^L, ixO^L, wCT, w, x, coolrate, fl)
1314 ! Calculates cooling rate using the exact cooling method,
1315 ! for usage in eg. source_terms subroutine.
1316 ! this routine can only be used together with the "exact" cooling method.
1317 use mod_global_parameters
1318
1319 integer, intent(in) :: ixI^L, ixO^L
1320 double precision, intent(in) :: qdt, x(ixI^S, 1:ndim), wCT(ixI^S, 1:nw)
1321 double precision :: w(ixI^S, 1:nw)
1322 double precision, intent(out) :: coolrate(ixI^S)
1323 type(rc_fluid), intent(in) :: fl
1324 double precision :: y1, y2, l1, tlocal2
1325 double precision :: Te(ixI^S), pnew(ixI^S), rho(ixI^S), rhonew(ixI^S)
1326 double precision :: emin, Lmax, fact, Rfactor(ixI^S), pth(ixI^S)
1327 integer :: ix^D
1328
1329 ! Check cooling method
1330 if( fl%coolmethod /= 'exact') then
1331 call mpistop("Subroutine getvar_cooling_exact needs the exact cooling method")
1332 end if
1333
1334 call fl%get_pthermal(wct, x, ixi^l, ixo^l, pth)
1335 call fl%get_rho(wct, x, ixi^l, ixo^l, rho)
1336 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1337 te(ixo^s)=pth(ixo^s)/(rho(ixo^s)*rfactor(ixo^s))
1338
1339 call fl%get_pthermal(w, x, ixi^l, ixo^l, pnew)
1340 call fl%get_rho(w, x, ixi^l, ixo^l, rhonew)
1341
1342 fact=fl%lref*qdt/fl%tref
1343
1344 {do ix^db = ixo^lim^db\}
1345 emin = rhonew(ix^d) * fl%tlow * rfactor(ix^d) * invgam
1346 lmax = max(zero, ( pnew(ix^d)*invgam - emin ) / qdt)
1347
1348 ! No cooling if temperature is below floor level.
1349 ! Assuming Bremsstrahlung if temperature is higher than maximum.
1350 if( te(ix^d)<= fl%tcoolmin) then
1351 l1 = zero
1352 else if( te(ix^d)>= fl%tcoolmax ) then
1353 call calc_l_extended(te(ix^d), l1, fl)
1354 l1 = l1 * rho(ix^d)**2
1355 l1 = min(l1, lmax)
1356 else
1357 call findy(te(ix^d), y1, fl)
1358 y2 = y1 + fact * rho(ix^d)*rc_gamma_1
1359 call findt(tlocal2, y2, fl)
1360 if( tlocal2 <= fl%tcoolmin ) then
1361 l1 = lmax
1362 else
1363 l1 = (te(ix^d)- tlocal2)*rho(ix^d)*rfactor(ix^d)*invgam/qdt
1364 end if
1365 l1 = min(l1, lmax)
1366 end if
1367 if(slab_uniform .and. fl%rad_cut .and. x(ix^d,ndim) .le. fl%rad_cut_hgt) then
1368 l1 = l1*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1369 end if
1370 coolrate(ix^d) = l1
1371 {end do\}
1372 end subroutine getvar_cooling_exact
1373
1374 subroutine radiative_cooling_add_source(qdt,ixI^L,ixO^L,wCT,wCTprim,w,x,&
1375 qsourcesplit,active,fl)
1376 ! w[iw]=w[iw]+qdt*S[wCT,x] where S is the source based on wCT within ixO
1377 use mod_global_parameters
1378 integer, intent(in) :: ixI^L, ixO^L
1379 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw), wCTprim(ixI^S,1:nw)
1380 double precision, intent(inout) :: w(ixI^S,1:nw)
1381 logical, intent(in) :: qsourcesplit
1382 logical, intent(inout) :: active
1383 type(rc_fluid), intent(in) :: fl
1384 double precision, allocatable, dimension(:^D&) :: Lequi
1385
1386 if(qsourcesplit .eqv.fl%rc_split) then
1387 active = .true.
1388 select case(fl%coolmethod)
1389 case ('explicit1')
1390 if(mype==0)then
1391 if(it==1) then
1392 write(*,*)'Fully explicit cooling is not completely safe in this version'
1393 write(*,*)'PROCEED WITH CAUTION!'
1394 endif
1395 endif
1396 call cool_explicit1(qdt,ixi^l,ixo^l,wct,w,x,fl)
1397 case ('explicit2')
1398 call cool_explicit2(qdt,ixi^l,ixo^l,wct,w,x,fl)
1399 case ('semiimplicit')
1400 call cool_semiimplicit(qdt,ixi^l,ixo^l,wct,w,x,fl)
1401 case ('implicit')
1402 call cool_implicit(qdt,ixi^l,ixo^l,wct,w,x,fl)
1403 case ('exact')
1404 call cool_exact(qdt,ixi^l,ixo^l,wct,wctprim,w,x,fl)
1405 case default
1406 call mpistop("This cooling method is unknown")
1407 end select
1408 if(fl%has_equi) then
1409 allocate(lequi(ixi^s))
1410 call get_cool_equi(qdt,ixi^l,ixo^l,wct,w,x,fl,lequi)
1411 w(ixo^s,fl%e_) = w(ixo^s,fl%e_)+lequi(ixo^s)
1412 deallocate(lequi)
1413 endif
1414 if( fl%Tfix ) call floortemperature(qdt,ixi^l,ixo^l,wct,w,x,fl)
1415 end if
1416 end subroutine radiative_cooling_add_source
1417
1418 subroutine floortemperature(qdt,ixI^L,ixO^L,wCT,w,x,fl)
1419 ! Force minimum temperature to a fixed temperature
1420 use mod_global_parameters
1421 integer, intent(in) :: ixI^L, ixO^L
1422 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw)
1423 double precision, intent(inout) :: w(ixI^S,1:nw)
1424 type(rc_fluid), intent(in) :: fl
1425 double precision :: etherm(ixI^S), rho(ixI^S), Rfactor(ixI^S),emin
1426 integer :: ix^D
1427
1428 call fl%get_pthermal(w,x,ixi^l,ixo^l,etherm)
1429 call fl%get_rho(w,x,ixi^l,ixo^l,rho)
1430 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1431 {do ix^db = ixo^lim^db\}
1432 emin = rho(ix^d)*fl%tlow*rfactor(ix^d)
1433 if(etherm(ix^d) < emin) then
1434 w(ix^d,fl%e_)=w(ix^d,fl%e_)+(emin-etherm(ix^d))*invgam
1435 end if
1436 {end do\}
1437 end subroutine floortemperature
1438
1439 subroutine get_cool_equi(qdt,ixI^L,ixO^L,wCT,w,x,fl,res)
1440 ! explicit cooling routine that depends on getdt to
1441 ! adjust the timestep. Accurate but incredibly slow
1442 use mod_global_parameters
1443
1444 integer, intent(in) :: ixI^L, ixO^L
1445 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw)
1446 double precision, intent(inout) :: w(ixI^S,1:nw)
1447 type(rc_fluid), intent(in) :: fl
1448 double precision, intent(out) :: res(ixI^S)
1449
1450 double precision :: pth(ixI^S),rho(ixI^S),Rfactor(ixI^S),L1,Tlocal2
1451 double precision :: Te(ixI^S)
1452 double precision :: emin, Lmax
1453 double precision :: Y1, Y2
1454 double precision :: de, emax,fact
1455 integer :: ix^D
1456
1457 call fl%get_pthermal_equi(wct,x,ixi^l,ixo^l,pth)
1458 call fl%get_rho_equi(wct,x,ixi^l,ixo^l,rho)
1459 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1460 te(ixo^s)=pth(ixo^s)/(rho(ixo^s)*rfactor(ixo^s))
1461
1462 res=0d0
1463
1464 if(fl%coolmethod == 'exact') then
1465
1466 fact = fl%lref*qdt/fl%tref
1467 {do ix^db = ixo^lim^db\}
1468 emin = rho(ix^d)*fl%tlow*rfactor(ix^d)*invgam
1469 lmax = max(zero,(pth(ix^d)*invgam-emin)/qdt)
1470 emax = max(zero, pth(ix^d)*invgam-emin)
1471 ! Determine explicit cooling
1472 ! If temperature is below floor level, no cooling.
1473 ! Stop wasting time and go to next gridpoint.
1474 ! If the temperature is higher than the maximum,
1475 ! assume Bremsstrahlung
1476 if( te(ix^d)<=fl%tcoolmin ) then
1477 l1 = zero
1478 else if( te(ix^d)>=fl%tcoolmax )then
1479 call calc_l_extended(te(ix^d), l1,fl)
1480 l1 = l1*rho(ix^d)**2
1481 if(phys_trac) then
1482 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1483 l1=l1*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1484 end if
1485 end if
1486 l1 = min(l1,lmax)
1487 res(ix^d) = l1*qdt
1488 else
1489 call findy(te(ix^d),y1,fl)
1490 y2 = y1 + fact * rho(ix^d)*rc_gamma_1
1491 call findt(tlocal2,y2,fl)
1492 if(tlocal2<=fl%tcoolmin) then
1493 de = emax
1494 else
1495 de = (te(ix^d)-tlocal2)*rho(ix^d)*rfactor(ix^d)*invgam
1496 end if
1497 if(phys_trac) then
1498 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1499 de=de*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1500 end if
1501 end if
1502 de = min(de,emax)
1503 res(ix^d) = de
1504 end if
1505 {end do\}
1506 else
1507 {do ix^db = ixo^lim^db\}
1508 emin = rho(ix^d)*fl%tlow*rfactor(ix^d)*invgam
1509 lmax = max(zero,pth(ix^d)*invgam-emin)/qdt
1510 ! Determine explicit cooling
1511 ! If temperature is below floor level, no cooling.
1512 ! Stop wasting time and go to next gridpoint.
1513 ! If the temperature is higher than the maximum,
1514 ! assume Bremsstrahlung
1515 if( te(ix^d)<=fl%tcoolmin ) then
1516 l1 = zero
1517 else if( te(ix^d)>=fl%tcoolmax )then
1518 call calc_l_extended(te(ix^d), l1,fl)
1519 else
1520 call findl(te(ix^d),l1,fl)
1521 end if
1522 l1 = l1*rho(ix^d)**2
1523 if(phys_trac) then
1524 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1525 l1=l1*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1526 end if
1527 end if
1528 l1 = min(l1,lmax)
1529 res(ix^d) =l1*qdt
1530 {end do\}
1531 end if
1532 end subroutine get_cool_equi
1533
1534 subroutine cool_explicit1(qdt,ixI^L,ixO^L,wCT,w,x,fl)
1535 ! explicit cooling routine that depends on getdt to
1536 ! adjust the timestep. Accurate but incredibly slow
1537 use mod_global_parameters
1538
1539 integer, intent(in) :: ixI^L, ixO^L
1540 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw)
1541 double precision, intent(inout) :: w(ixI^S,1:nw)
1542 type(rc_fluid), intent(in) :: fl
1543
1544 double precision :: L1,pth(ixI^S),pnew(ixI^S),rho(ixI^S),Rfactor(ixI^S)
1545 double precision :: Te(ixI^S)
1546 double precision :: emin, Lmax
1547 integer :: ix^D
1548
1549 call fl%get_pthermal(wct,x,ixi^l,ixo^l,pth)
1550 call fl%get_pthermal(w,x,ixi^l,ixo^l,pnew)
1551 call fl%get_rho(wct,x,ixi^l,ixo^l,rho)
1552 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1553 te(ixo^s)=pth(ixo^s)/(rho(ixo^s)*rfactor(ixo^s))
1554
1555 {do ix^db = ixo^lim^db\}
1556 emin = rho(ix^d)*fl%tlow*rfactor(ix^d)*invgam
1557 lmax = max(zero,pnew(ix^d)*invgam-emin)/qdt
1558 ! Determine explicit cooling
1559 ! If temperature is below floor level, no cooling.
1560 ! Stop wasting time and go to next gridpoint.
1561 ! If the temperature is higher than the maximum,
1562 ! assume Bremsstrahlung
1563 if( te(ix^d)<=fl%tcoolmin ) then
1564 l1 = zero
1565 else if( te(ix^d)>=fl%tcoolmax )then
1566 call calc_l_extended(te(ix^d), l1,fl)
1567 l1 = l1*rho(ix^d)**2
1568 if(phys_trac) then
1569 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1570 l1=l1*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1571 end if
1572 end if
1573 l1 = min(l1,lmax)
1574 else
1575 call findl(te(ix^d),l1,fl)
1576 l1 = l1*rho(ix^d)**2
1577 if(phys_trac) then
1578 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1579 l1=l1*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1580 end if
1581 end if
1582 l1 = min(l1,lmax)
1583 end if
1584 if(slab_uniform .and. fl%rad_cut .and. x(ix^d,ndim) .le. fl%rad_cut_hgt) then
1585 l1 = l1*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1586 end if
1587 w(ix^d,fl%e_) = w(ix^d,fl%e_)-l1*qdt
1588 {end do\}
1589 end subroutine cool_explicit1
1590
1591 subroutine cool_explicit2(qdt,ixI^L,ixO^L,wCT,w,x,fl)
1592 ! explicit cooling routine that does a series
1593 ! of small forward integration steps, to make
1594 ! sure the amount of cooling remains correct
1595 ! Not as accurate as 'explicit1', but a lot faster
1596 ! tends to overestimate cooling
1597 use mod_global_parameters
1598 integer, intent(in) :: ixI^L, ixO^L
1599 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw)
1600 double precision, intent(inout) :: w(ixI^S,1:nw)
1601 type(rc_fluid), intent(in) :: fl
1602 double precision :: Ltest, etherm, de
1603 double precision :: dtmax, dtstep
1604 double precision :: L1,pth(ixI^S),pnew(ixI^S),rho(ixI^S),Rfactor(ixI^S)
1605 double precision :: Tlocal1,plocal,Te(ixI^S)
1606 double precision :: emin, Lmax
1607 integer :: idt,ndtstep
1608 integer :: ix^D
1609
1610 call fl%get_pthermal(wct,x,ixi^l,ixo^l,pth)
1611 call fl%get_pthermal(w,x,ixi^l,ixo^l,pnew)
1612 call fl%get_rho(wct,x,ixi^l,ixo^l,rho)
1613 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1614 te(ixo^s)=pth(ixo^s)/(rho(ixo^s)*rfactor(ixo^s))
1615
1616 {do ix^db = ixo^lim^db\}
1617 ! Calculate explicit cooling value
1618 dtmax = qdt
1619 etherm = pth(ix^d)*invgam
1620 emin = rho(ix^d)*fl%tlow*rfactor(ix^d)*invgam
1621 lmax = max(zero,pnew(ix^d)*invgam-emin)/qdt
1622 ! Determine explicit cooling
1623 ! If temperature is below floor level, no cooling.
1624 ! Stop wasting time and go to next gridpoint.
1625 ! If the temperature is higher than the maximum,
1626 ! assume Bremmstrahlung
1627 if( te(ix^d)<=fl%tcoolmin ) then
1628 ltest = zero
1629 else if( te(ix^d)>=fl%tcoolmax )then
1630 call calc_l_extended(te(ix^d), ltest,fl)
1631 ltest = l1*rho(ix^d)**2
1632 if(phys_trac) then
1633 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1634 ltest=ltest*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1635 end if
1636 end if
1637 ltest = min(l1,lmax)
1638 if( dtmax>fl%cfrac*etherm/ltest) dtmax = fl%cfrac*etherm/ltest
1639 else
1640 call findl(te(ix^d),ltest,fl)
1641 ltest = ltest*rho(ix^d)**2
1642 if(phys_trac) then
1643 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1644 ltest=ltest*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1645 end if
1646 end if
1647 ltest = min(ltest,lmax)
1648 if( dtmax>fl%cfrac*etherm/ltest) dtmax = fl%cfrac*etherm/ltest
1649 end if
1650 ! Calculate number of steps for cooling
1651 ndtstep = max(nint(qdt/dtmax),1)+1
1652 dtstep = qdt/ndtstep
1653 ! Use explicit cooling value for first step
1654 de = ltest*dtstep
1655 etherm = etherm - de
1656
1657 do idt=2,ndtstep
1658 plocal = etherm*rc_gamma_1
1659 lmax = max(zero,etherm-emin)/dtstep
1660 ! Tlocal = P/(rho*R)
1661 tlocal1 = plocal/(rho(ix^d)*rfactor(ix^d))
1662 if( tlocal1<=fl%tcoolmin ) then
1663 l1 = zero
1664 exit
1665 else if( tlocal1>=fl%tcoolmax )then
1666 call calc_l_extended(tlocal1, l1,fl)
1667 l1 = l1*rho(ix^d)**2
1668 if(phys_trac) then
1669 if(tlocal1<block%wextra(ix^d,fl%Tcoff_)) then
1670 l1=l1*sqrt((tlocal1/block%wextra(ix^d,fl%Tcoff_))**5)
1671 end if
1672 end if
1673 l1 = min(l1,lmax)
1674 else
1675 call findl(tlocal1,l1,fl)
1676 l1 = l1*rho(ix^d)**2
1677 if(phys_trac) then
1678 if(tlocal1<block%wextra(ix^d,fl%Tcoff_)) then
1679 l1=l1*sqrt((tlocal1/block%wextra(ix^d,fl%Tcoff_))**5)
1680 end if
1681 end if
1682 l1 = min(l1,lmax)
1683 end if
1684 de = de + l1*dtstep
1685 etherm = etherm - l1*dtstep
1686 end do
1687 if(slab_uniform .and. fl%rad_cut .and. x(ix^d,ndim) .le. fl%rad_cut_hgt) then
1688 de = de*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1689 end if
1690 w(ix^d,fl%e_) = w(ix^d,fl%e_) -de
1691 {end do\}
1692 end subroutine cool_explicit2
1693
1694 subroutine cool_semiimplicit(qdt,ixI^L,ixO^L,wCT,w,x,fl)
1695 ! Semi-implicit cooling method based on a two point average
1696 ! Fast, but tends to underestimate cooling
1697 use mod_global_parameters
1698 integer, intent(in) :: ixI^L, ixO^L
1699 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw)
1700 double precision, intent(inout) :: w(ixI^S,1:nw)
1701 type(rc_fluid), intent(in) :: fl
1702 double precision :: L1,L2,Tlocal2
1703 double precision :: etemp
1704 double precision :: emin, Lmax
1705 double precision :: pth(ixI^S),pnew(ixI^S),rho(ixI^S),Rfactor(ixI^S),Te(ixI^S)
1706 integer :: ix^D
1707
1708 call fl%get_pthermal(wct,x,ixi^l,ixo^l,pth)
1709 call fl%get_pthermal(w,x,ixi^l,ixo^l,pnew)
1710 call fl%get_rho(wct,x,ixi^l,ixo^l,rho)
1711 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1712 te(ixo^s)=pth(ixo^s)/(rho(ixo^s)*rfactor(ixo^s))
1713
1714 {do ix^db = ixo^lim^db\}
1715 emin = rho(ix^d)*fl%tlow*rfactor(ix^d)*invgam
1716 lmax = max(zero,pnew(ix^d)*invgam-emin)/qdt
1717 ! Determine explicit cooling at present temperature
1718 !
1719 ! If temperature is below floor level, no cooling.
1720 ! Stop wasting time and go to next gridpoint.
1721 ! If the temperature is higher than the maximum,
1722 ! assume Bremsstrahlung
1723 if( te(ix^d)<=fl%tcoolmin ) then
1724 l1 = zero
1725 l2 = zero
1726 else
1727 if( te(ix^d)>=fl%tcoolmax ) then
1728 call calc_l_extended(te(ix^d), l1,fl)
1729 else
1730 call findl(te(ix^d),l1,fl)
1731 end if
1732 l1 = l1*rho(ix^d)**2
1733 if(phys_trac) then
1734 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1735 l1=l1*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1736 end if
1737 end if
1738 etemp = pth(ix^d)*invgam - l1*qdt
1739 tlocal2 = etemp*rc_gamma_1/(rho(ix^d)*rfactor(ix^d))
1740 ! Determine explicit cooling at new temperature
1741 if( tlocal2<=fl%tcoolmin ) then
1742 l2 = zero
1743 else if( tlocal2>=fl%tcoolmax )then
1744 call calc_l_extended(tlocal2, l2,fl)
1745 else
1746 call findl(tlocal2,l2,fl)
1747 end if
1748 l2 = l2*rho(ix^d)**2
1749 if(phys_trac) then
1750 if(tlocal2<block%wextra(ix^d,fl%Tcoff_)) then
1751 l2=l2*sqrt((tlocal2/block%wextra(ix^d,fl%Tcoff_))**5)
1752 end if
1753 end if
1754 if(slab_uniform .and. fl%rad_cut .and. x(ix^d,ndim) .le. fl%rad_cut_hgt) then
1755 l1 = l1*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1756 l2 = l2*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1757 end if
1758 w(ix^d,fl%e_) = w(ix^d,fl%e_) - min(half*(l1+l2),lmax)*qdt
1759 end if
1760 {end do\}
1761 end subroutine cool_semiimplicit
1762
1763 subroutine cool_implicit(qdt,ixI^L,ixO^L,wCT,w,x,fl)
1764 ! Implicit cooling method based on a half-step
1765 ! refinement algorithm
1766 use mod_global_parameters
1767 integer, intent(in) :: ixI^L, ixO^L
1768 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw)
1769 double precision, intent(inout) :: w(ixI^S,1:nw)
1770 type(rc_fluid), intent(in) :: fl
1771 double precision :: Ltemp,Tnew,f1,f2,pth(ixI^S), pnew(ixI^S), rho(ixI^S), Rfactor(ixI^S)
1772 double precision :: elocal, Te(ixI^S)
1773 double precision :: emin, Lmax, eold, enew, estep
1774 double precision, parameter :: e_error = 1.0d-6
1775 integer, parameter :: maxiter = 100
1776 integer :: ix^D, j
1777
1778 call fl%get_pthermal(wct,x,ixi^l,ixo^l,pth)
1779 call fl%get_pthermal(w,x,ixi^l,ixo^l,pnew)
1780 call fl%get_rho(wct,x,ixi^l,ixo^l,rho)
1781 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1782 te(ixo^s)=pth(ixo^s)/(rho(ixo^s)*rfactor(ixo^s))
1783
1784 {do ix^db = ixo^lim^db\}
1785 elocal = pth(ix^d)*invgam
1786 emin = rho(ix^d)*fl%tlow*rfactor(ix^d)*invgam
1787 lmax = max(zero,pnew(ix^d)*invgam-emin)/qdt
1788 ! Determine explicit cooling at present temperature
1789 ! If temperature is below floor level, no cooling.
1790 ! Stop wasting time and go to next gridpoint.
1791 ! If the temperature is higher than the maximum,
1792 ! assume Bremsstrahlung
1793 if( te(ix^d)<=fl%tcoolmin ) then
1794 ltemp = zero
1795 else
1796 eold = elocal
1797 enew = elocal
1798 estep = -(smalldouble)
1799 f2 = 1.d0
1800 do j=1,maxiter+1
1801 if( j>maxiter ) call mpistop("Implicit cooling exceeds maximum iterations")
1802 tnew = enew*rc_gamma_1/(rho(ix^d)*rfactor(ix^d))
1803 if( tnew<=fl%tcoolmin ) then
1804 ltemp = lmax
1805 exit
1806 else if( tnew>=fl%tcoolmax )then
1807 call calc_l_extended(tnew, ltemp,fl)
1808 else
1809 call findl(tnew,ltemp,fl)
1810 end if
1811 ltemp = ltemp*rho(ix^d)**2
1812 eold = enew + ltemp*qdt
1813 f1 = elocal -eold
1814 if(abs(half*f1/(elocal+eold)) < e_error) exit
1815 if(phys_trac) then
1816 if(tnew<block%wextra(ix^d,fl%Tcoff_)) then
1817 ltemp=ltemp*sqrt((tnew/block%wextra(ix^d,fl%Tcoff_))**5)
1818 end if
1819 end if
1820 if(j==1) estep = max((elocal-emin)*half,smalldouble)
1821 if(f1*f2 < zero) estep = -half*estep
1822 f2 = f1
1823 enew = enew +estep
1824 end do
1825 end if
1826 if(slab_uniform .and. fl%rad_cut .and. x(ix^d,ndim) .le. fl%rad_cut_hgt) then
1827 ltemp = ltemp*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1828 end if
1829 w(ix^d,fl%e_) = w(ix^d,fl%e_) - min(ltemp,lmax)*qdt
1830 {end do\}
1831 end subroutine cool_implicit
1832
1833 subroutine cool_exact(qdt,ixI^L,ixO^L,wCT,wCTprim,w,x,fl)
1834 ! Cooling routine using exact integration method from Townsend 2009
1835 use mod_global_parameters
1836 integer, intent(in) :: ixI^L, ixO^L
1837 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw), wCTprim(ixI^S,1:nw)
1838 double precision, intent(inout) :: w(ixI^S,1:nw)
1839 type(rc_fluid), intent(in) :: fl
1840 double precision :: Y1, Y2
1841 double precision :: L1, pth(ixI^S), Tlocal2, pnew(ixI^S)
1842 double precision :: rho(ixI^S), Te(ixI^S), rhonew(ixI^S), Rfactor(ixI^S)
1843 double precision :: emin, Lmax, fact
1844 double precision :: de, emax
1845 integer :: ix^D
1846
1847 call fl%get_rho(wct,x,ixi^l,ixo^l,rho)
1848 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1849 if(phys_equi_pe) then
1850 ! need pressure splitting
1851 call fl%get_pthermal(wct,x,ixi^l,ixo^l,te)
1852 te(ixo^s)=te(ixo^s)/(rho(ixo^s)*rfactor(ixo^s))
1853 else
1854 te(ixo^s)=wctprim(ixo^s,iw_e)/(rho(ixo^s)*rfactor(ixo^s))
1855 end if
1856 call fl%get_pthermal(w,x,ixi^l,ixo^l,pnew)
1857 call fl%get_rho(w,x,ixi^l,ixo^l,rhonew)
1858
1859 fact = fl%lref*qdt/fl%tref
1860
1861 {do ix^db = ixo^lim^db\}
1862 emin = rhonew(ix^d)*fl%tlow*rfactor(ix^d)*invgam
1863 lmax = max(zero,pnew(ix^d)*invgam-emin)/qdt
1864 emax = max(zero,pnew(ix^d)*invgam-emin)
1865 ! Determine explicit cooling
1866 ! If temperature is below floor level, no cooling.
1867 ! Stop wasting time and go to next gridpoint.
1868 ! If the temperature is higher than the maximum,
1869 ! assume Bremsstrahlung
1870 if( te(ix^d)<=fl%tcoolmin ) then
1871 l1 = zero
1872 else if( te(ix^d)>=fl%tcoolmax )then
1873 call calc_l_extended(te(ix^d), l1,fl)
1874 l1 = l1*rho(ix^d)**2
1875 if(phys_trac) then
1876 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1877 l1=l1*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1878 end if
1879 end if
1880 l1 = min(l1,lmax)
1881 if(slab_uniform .and. fl%rad_cut .and. x(ix^d,ndim) .le. fl%rad_cut_hgt) then
1882 l1 = l1*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1883 end if
1884 w(ix^d,fl%e_) = w(ix^d,fl%e_)-l1*qdt
1885 else
1886 call findy(te(ix^d),y1,fl)
1887 y2 = y1 + fact*rho(ix^d)*rc_gamma_1
1888 call findt(tlocal2,y2,fl)
1889 if(tlocal2<=fl%tcoolmin) then
1890 de = emax
1891 else
1892 de = (te(ix^d)-tlocal2)*rho(ix^d)*rfactor(ix^d)*invgam
1893 end if
1894 if(phys_trac) then
1895 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1896 de=de*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1897 end if
1898 end if
1899 de = min(de,emax)
1900 if(slab_uniform .and. fl%rad_cut .and. x(ix^d,ndim) .le. fl%rad_cut_hgt) then
1901 de = de*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1902 end if
1903 w(ix^d,fl%e_) = w(ix^d,fl%e_)-de
1904 end if
1905 {end do\}
1906 end subroutine cool_exact
1907
1908 subroutine calc_l_extended (tpoint, lpoint,fl)
1909 ! Calculate l for t beyond tcoolmax
1910 ! Assumes Bremsstrahlung for the interpolated tables
1911 ! Uses the power law for piecewise power laws
1912 double precision, intent(IN) :: tpoint
1913 double precision, intent(OUT) :: lpoint
1914 type(rc_fluid), intent(in) :: fl
1915
1916 if(fl%isPPL) then
1917 lpoint =fl%l_PPL(fl%n_PPL) * ( tpoint / fl%t_PPL(fl%n_PPL) )**fl%a_PPL(fl%n_PPL)
1918 else
1919 lpoint = fl%Lcool(fl%ncool) * sqrt( tpoint / fl%tcoolmax)
1920 end if
1921 end subroutine calc_l_extended
1922
1923 subroutine findl (tpoint,Lpoint,fl)
1924 ! Fast search option to find correct point
1925 ! in cooling curve
1926 use mod_global_parameters
1927
1928 double precision,intent(IN) :: tpoint
1929 double precision, intent(OUT) :: Lpoint
1930 type(rc_fluid), intent(in) :: fl
1931
1932 double precision :: lgtp
1933 integer :: jl,i
1934
1935 if(fl%isPPL) then
1936 i = maxloc(fl%t_PPL, dim=1, mask=fl%t_PPL<tpoint)
1937 lpoint = fl%l_PPL(i) * (tpoint / fl%t_PPL(i))**fl%a_PPL(i)
1938 else
1939 lgtp = dlog10(tpoint)
1940 jl = int((lgtp - fl%lgtcoolmin) /fl%lgstep) + 1
1941 lpoint = fl%Lcool(jl)+ (tpoint-fl%tcool(jl)) &
1942 * (fl%Lcool(jl+1)-fl%Lcool(jl)) &
1943 / (fl%tcool(jl+1)-fl%tcool(jl))
1944 end if
1945
1946 end subroutine findl
1947
1948 subroutine findy (tpoint,Ypoint,fl)
1949 ! Fast search option to find correct point in cooling time
1950 use mod_global_parameters
1951
1952 double precision,intent(IN) :: tpoint
1953 double precision, intent(OUT) :: Ypoint
1954 type(rc_fluid), intent(in) :: fl
1955
1956 double precision :: lgtp
1957 double precision :: y_extra,factor
1958 integer :: jl,i
1959
1960 if(fl%isPPL) then
1961 i = maxloc(fl%t_PPL, dim=1, mask=fl%t_PPL<tpoint)
1962 factor = fl%l_PPL(fl%n_PPL+1) * fl%t_PPL(i) / (fl%l_PPL(i) * fl%t_PPL(fl%n_PPL+1))
1963 if(fl%a_PPL(i)==1.d0) then
1964 y_extra = log( fl%t_PPL(i) / tpoint )
1965 else
1966 y_extra = 1 / (1 - fl%a_PPL(i)) * (1 - ( fl%t_PPL(i) / tpoint )**(fl%a_PPL(i)-1) )
1967 end if
1968 ypoint = fl%y_PPL(i) + factor*y_extra
1969 else
1970 lgtp = dlog10(tpoint)
1971 jl = int((lgtp - fl%lgtcoolmin) / fl%lgstep) + 1
1972 ypoint = fl%Yc(jl)+ (tpoint-fl%tcool(jl)) &
1973 * (fl%Yc(jl+1)-fl%Yc(jl)) &
1974 / (fl%tcool(jl+1)-fl%tcool(jl))
1975 end if
1976
1977 end subroutine findy
1978
1979 subroutine findt (tpoint,Ypoint,fl)
1980 ! Fast search option to find correct temperature
1981 ! from temporal evolution function. Only possible this way because T is a monotonously
1982 ! decreasing function for the interpolated tables
1983 ! Uses eq. A7 from Townsend 2009 for piecewise power laws
1984 use mod_global_parameters
1985
1986 double precision,intent(OUT) :: tpoint
1987 double precision, intent(IN) :: Ypoint
1988 type(rc_fluid), intent(in) :: fl
1989
1990 double precision :: factor
1991 integer :: jl,jc,jh,i
1992
1993 if(fl%isPPL) then
1994 i = minloc(fl%y_PPL, dim=1, mask=fl%y_PPL>ypoint)
1995 factor = fl%l_PPL(i) * fl%t_PPL(fl%n_PPL+1) / (fl%l_PPL(fl%n_PPL+1) * fl%t_PPL(i))
1996 if(fl%a_PPL(i)==1.d0) then
1997 tpoint = fl%t_PPL(i) * exp( -1.d0 * factor * ( ypoint - fl%y_PPL(i)))
1998 else
1999 tpoint = fl%t_PPL(i) * (1 - (1 - fl%a_PPL(i)) * factor * (ypoint - fl%y_PPL(i)))**(1 / (1 - fl%a_PPL(i)))
2000 end if
2001 else
2002 if(ypoint >= fl%Yc(1)) then
2003 tpoint = fl%tcoolmin
2004 else if (ypoint == fl%Yc(fl%ncool)) then
2005 tpoint = fl%tcoolmax
2006 else
2007 jl=0
2008 jh=fl%ncool+1
2009 do
2010 if(jh-jl <= 1) exit
2011 jc=(jh+jl)/2
2012 if(ypoint <= fl%Yc(jc)) then
2013 jl=jc
2014 else
2015 jh=jc
2016 end if
2017 end do
2018 ! Linear interpolation to obtain correct temperature
2019 tpoint = fl%tcool(jl)+ (ypoint-fl%Yc(jl)) &
2020 * (fl%tcool(jl+1)-fl%tcool(jl)) &
2021 / (fl%Yc(jl+1)-fl%Yc(jl))
2022 end if
2023 end if
2024 end subroutine findt
2025
2026end module mod_radiative_cooling
mod_radiative_cooling::get_subr1
Definition mod_radiative_cooling.t:40
mod_comm_lib
Definition mod_comm_lib.t:1
mod_comm_lib::mpistop
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition mod_comm_lib.t:208
mod_global_parameters
This module contains definitions of global parameters and variables and some generic functions/subrou...
Definition mod_global_parameters.t:4
mod_global_parameters::block
type(state), pointer block
Block pointer for using one block and its previous state.
Definition mod_global_parameters.t:769
mod_global_parameters::unitpar
integer, parameter unitpar
file handle for IO
Definition mod_global_parameters.t:370
mod_global_parameters::it
integer it
Number of time steps taken.
Definition mod_global_parameters.t:436
mod_global_parameters::par_files
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
Definition mod_global_parameters.t:763
mod_global_parameters::mype
integer mype
The rank of the current MPI task.
Definition mod_global_parameters.t:221
mod_global_parameters::d
double precision, dimension(:), allocatable, parameter d
Definition mod_global_parameters.t:23
mod_global_parameters::phys_trac
logical phys_trac
Use TRAC for MHD or 1D HD.
Definition mod_global_parameters.t:662
mod_global_parameters::slab_uniform
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
Definition mod_global_parameters.t:566
mod_physics
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition mod_physics.t:4
mod_radiative_cooling
module radiative cooling – add optically thin radiative cooling for HD and MHD
Definition mod_radiative_cooling.t:13
mod_radiative_cooling::t_fm
double precision, dimension(1:5) t_fm
Definition mod_radiative_cooling.t:108
mod_radiative_cooling::n_mlcosmol
integer n_mlcosmol
Definition mod_radiative_cooling.t:189
mod_radiative_cooling::l_mlcosmol
double precision, dimension(1:71) l_mlcosmol
Definition mod_radiative_cooling.t:181
mod_radiative_cooling::n_jccorona
integer n_jccorona
Definition mod_radiative_cooling.t:189
mod_radiative_cooling::getvar_cooling
subroutine getvar_cooling(ixil, ixol, w, x, coolrate, fl)
Definition mod_radiative_cooling.t:1274
mod_radiative_cooling::l_mb
double precision, dimension(1:51) l_mb
Definition mod_radiative_cooling.t:181
mod_radiative_cooling::cool_implicit
subroutine cool_implicit(qdt, ixil, ixol, wct, w, x, fl)
Definition mod_radiative_cooling.t:1764
mod_radiative_cooling::cool_explicit1
subroutine cool_explicit1(qdt, ixil, ixol, wct, w, x, fl)
Definition mod_radiative_cooling.t:1535
mod_radiative_cooling::t_cl_solar
double precision, dimension(1:151) t_cl_solar
Definition mod_radiative_cooling.t:176
mod_radiative_cooling::l_dm_2
double precision, dimension(1:76) l_dm_2
Definition mod_radiative_cooling.t:181
mod_radiative_cooling::radiative_cooling_init_params
subroutine radiative_cooling_init_params(phys_gamma, he_abund)
Radiative cooling initialization.
Definition mod_radiative_cooling.t:770
mod_radiative_cooling::findl
subroutine findl(tpoint, lpoint, fl)
Definition mod_radiative_cooling.t:1924
mod_radiative_cooling::l_spex
double precision, dimension(1:110) l_spex
Definition mod_radiative_cooling.t:181
mod_radiative_cooling::get_cool_equi
subroutine get_cool_equi(qdt, ixil, ixol, wct, w, x, fl, res)
Definition mod_radiative_cooling.t:1440
mod_radiative_cooling::cooling_get_dt
subroutine cooling_get_dt(w, ixil, ixol, dtnew, dxd, x, fl)
Definition mod_radiative_cooling.t:1233
mod_radiative_cooling::n_fm
integer n_fm
Definition mod_radiative_cooling.t:117
mod_radiative_cooling::a_hildner
double precision, dimension(1:5) a_hildner
Definition mod_radiative_cooling.t:114
mod_radiative_cooling::t_hildner
double precision, dimension(1:6) t_hildner
Definition mod_radiative_cooling.t:108
mod_radiative_cooling::radiative_cooling_init
subroutine radiative_cooling_init(fl, read_params)
Definition mod_radiative_cooling.t:778
mod_radiative_cooling::x_spex_dm_rough
double precision, dimension(1:7) x_spex_dm_rough
Definition mod_radiative_cooling.t:111
mod_radiative_cooling::l_colgan
double precision, dimension(1:55) l_colgan
Definition mod_radiative_cooling.t:181
mod_radiative_cooling::n_spex
integer n_spex
Definition mod_radiative_cooling.t:189
mod_radiative_cooling::t_rosner
double precision, dimension(1:10) t_rosner
Definition mod_radiative_cooling.t:108
mod_radiative_cooling::t_colgan
double precision, dimension(1:55) t_colgan
Definition mod_radiative_cooling.t:176
mod_radiative_cooling::a_fm
double precision, dimension(1:4) a_fm
Definition mod_radiative_cooling.t:114
mod_radiative_cooling::n_cl_solar
integer n_cl_solar
Definition mod_radiative_cooling.t:189
mod_radiative_cooling::n_cl_ism
integer n_cl_ism
Definition mod_radiative_cooling.t:189
mod_radiative_cooling::create_y_ppl
subroutine create_y_ppl(fl)
Definition mod_radiative_cooling.t:1209
mod_radiative_cooling::t_dere
double precision, dimension(1:101) t_dere
Definition mod_radiative_cooling.t:176
mod_radiative_cooling::n_hildner
integer n_hildner
Definition mod_radiative_cooling.t:117
mod_radiative_cooling::n_dere
integer n_dere
Definition mod_radiative_cooling.t:189
mod_radiative_cooling::n_dm
integer n_dm
Definition mod_radiative_cooling.t:189
mod_radiative_cooling::t_mlcosmol
double precision, dimension(1:71) t_mlcosmol
Definition mod_radiative_cooling.t:176
mod_radiative_cooling::t_spex_dm_fine
double precision, dimension(1:15) t_spex_dm_fine
Definition mod_radiative_cooling.t:108
mod_radiative_cooling::radiative_cooling_add_source
subroutine radiative_cooling_add_source(qdt, ixil, ixol, wct, wctprim, w, x, qsourcesplit, active, fl)
Definition mod_radiative_cooling.t:1376
mod_radiative_cooling::x_hildner
double precision, dimension(1:5) x_hildner
Definition mod_radiative_cooling.t:111
mod_radiative_cooling::calc_l_extended
subroutine calc_l_extended(tpoint, lpoint, fl)
Definition mod_radiative_cooling.t:1909
mod_radiative_cooling::t_mb
double precision, dimension(1:51) t_mb
Definition mod_radiative_cooling.t:176
mod_radiative_cooling::a_spex_dm_fine
double precision, dimension(1:14) a_spex_dm_fine
Definition mod_radiative_cooling.t:114
mod_radiative_cooling::cool_explicit2
subroutine cool_explicit2(qdt, ixil, ixol, wct, w, x, fl)
Definition mod_radiative_cooling.t:1592
mod_radiative_cooling::l_dere_photo
double precision, dimension(1:101) l_dere_photo
Definition mod_radiative_cooling.t:181
mod_radiative_cooling::n_spex_dm_fine
integer n_spex_dm_fine
Definition mod_radiative_cooling.t:117
mod_radiative_cooling::n_rosner
integer n_rosner
Definition mod_radiative_cooling.t:117
mod_radiative_cooling::cool_exact
subroutine cool_exact(qdt, ixil, ixol, wct, wctprim, w, x, fl)
Definition mod_radiative_cooling.t:1834
mod_radiative_cooling::l_jccorona
double precision, dimension(1:45) l_jccorona
Definition mod_radiative_cooling.t:181
mod_radiative_cooling::findt
subroutine findt(tpoint, ypoint, fl)
Definition mod_radiative_cooling.t:1980
mod_radiative_cooling::t_dm_2
double precision, dimension(1:76) t_dm_2
Definition mod_radiative_cooling.t:176
mod_radiative_cooling::n_dm_2
integer n_dm_2
Definition mod_radiative_cooling.t:189
mod_radiative_cooling::l_dere_corona
double precision, dimension(1:101) l_dere_corona
Definition mod_radiative_cooling.t:181
mod_radiative_cooling::l_cl_solar
double precision, dimension(1:151) l_cl_solar
Definition mod_radiative_cooling.t:181
mod_radiative_cooling::l_dm
double precision, dimension(1:71) l_dm
Definition mod_radiative_cooling.t:181
mod_radiative_cooling::x_rosner
double precision, dimension(1:9) x_rosner
Definition mod_radiative_cooling.t:111
mod_radiative_cooling::x_klimchuk
double precision, dimension(1:7) x_klimchuk
Definition mod_radiative_cooling.t:111
mod_radiative_cooling::a_spex_dm_rough
double precision, dimension(1:7) a_spex_dm_rough
Definition mod_radiative_cooling.t:114
mod_radiative_cooling::n_mlwc
integer n_mlwc
Definition mod_radiative_cooling.t:189
mod_radiative_cooling::x_fm
double precision, dimension(1:4) x_fm
Definition mod_radiative_cooling.t:111
mod_radiative_cooling::t_jccorona
double precision, dimension(1:45) t_jccorona
Definition mod_radiative_cooling.t:176
mod_radiative_cooling::t_mlwc
double precision, dimension(1:71) t_mlwc
Definition mod_radiative_cooling.t:176
mod_radiative_cooling::floortemperature
subroutine floortemperature(qdt, ixil, ixol, wct, w, x, fl)
Definition mod_radiative_cooling.t:1419
mod_radiative_cooling::t_mlsolar1
double precision, dimension(1:71) t_mlsolar1
Definition mod_radiative_cooling.t:176
mod_radiative_cooling::nenh_spex
double precision, dimension(1:110) nenh_spex
Definition mod_radiative_cooling.t:187
mod_radiative_cooling::n_mb
integer n_mb
Definition mod_radiative_cooling.t:189
mod_radiative_cooling::n_spex_dm_rough
integer n_spex_dm_rough
Definition mod_radiative_cooling.t:117
mod_radiative_cooling::t_cl_ism
double precision, dimension(1:151) t_cl_ism
Definition mod_radiative_cooling.t:176
mod_radiative_cooling::t_klimchuk
double precision, dimension(1:8) t_klimchuk
Definition mod_radiative_cooling.t:108
mod_radiative_cooling::n_colgan
integer n_colgan
Definition mod_radiative_cooling.t:189
mod_radiative_cooling::n_klimchuk
integer n_klimchuk
Definition mod_radiative_cooling.t:117
mod_radiative_cooling::l_cl_ism
double precision, dimension(1:151) l_cl_ism
Definition mod_radiative_cooling.t:181
mod_radiative_cooling::l_mlwc
double precision, dimension(1:71) l_mlwc
Definition mod_radiative_cooling.t:181
mod_radiative_cooling::cool_semiimplicit
subroutine cool_semiimplicit(qdt, ixil, ixol, wct, w, x, fl)
Definition mod_radiative_cooling.t:1695
mod_radiative_cooling::getvar_cooling_exact
subroutine getvar_cooling_exact(qdt, ixil, ixol, wct, w, x, coolrate, fl)
Definition mod_radiative_cooling.t:1314
mod_radiative_cooling::t_spex_dm_rough
double precision, dimension(1:8) t_spex_dm_rough
Definition mod_radiative_cooling.t:108
mod_radiative_cooling::a_klimchuk
double precision, dimension(1:7) a_klimchuk
Definition mod_radiative_cooling.t:114
mod_radiative_cooling::a_rosner
double precision, dimension(1:9) a_rosner
Definition mod_radiative_cooling.t:114
mod_radiative_cooling::l_mlsolar1
double precision, dimension(1:71) l_mlsolar1
Definition mod_radiative_cooling.t:181
mod_radiative_cooling::t_spex
double precision, dimension(1:110) t_spex
Definition mod_radiative_cooling.t:176
mod_radiative_cooling::x_spex_dm_fine
double precision, dimension(1:14) x_spex_dm_fine
Definition mod_radiative_cooling.t:111
mod_radiative_cooling::t_dm
double precision, dimension(1:71) t_dm
Definition mod_radiative_cooling.t:176
mod_radiative_cooling::n_mlsolar1
integer n_mlsolar1
Definition mod_radiative_cooling.t:189
mod_radiative_cooling::findy
subroutine findy(tpoint, ypoint, fl)
Definition mod_radiative_cooling.t:1949
mod_radiative_cooling::rc_fluid
Definition mod_radiative_cooling.t:49
Generated by doxygen 1.9.8