MPI-AMRVAC 3.1
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 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 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 ! The TEF must be known, so this routine can only be used
1317 ! together with the "exact" cooling method.
1318 use mod_global_parameters
1319
1320 integer, intent(in) :: ixI^L, ixO^L
1321 double precision, intent(in) :: qdt, x(ixI^S, 1:ndim), wCT(ixI^S, 1:nw)
1322 double precision :: w(ixI^S, 1:nw)
1323 double precision, intent(out) :: coolrate(ixI^S)
1324 type(rc_fluid), intent(in) :: fl
1325 double precision :: y1, y2, l1, tlocal2
1326 double precision :: Te(ixI^S), pnew(ixI^S), rho(ixI^S), rhonew(ixI^S)
1327 double precision :: emin, Lmax, fact, Rfactor(ixI^S), pth(ixI^S)
1328 integer :: ix^D
1329
1330 ! Check cooling method
1331 if( fl%coolmethod /= 'exact') then
1332 call mpistop("Subroutine getvar_cooling_exact needs the exact cooling method")
1333 end if
1334
1335 call fl%get_pthermal(wct, x, ixi^l, ixo^l, pth)
1336 call fl%get_rho(wct, x, ixi^l, ixo^l, rho)
1337 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1338 te(ixo^s)=pth(ixo^s)/(rho(ixo^s)*rfactor(ixo^s))
1339
1340 call fl%get_pthermal(w, x, ixi^l, ixo^l, pnew)
1341 call fl%get_rho(w, x, ixi^l, ixo^l, rhonew)
1342
1343 fact=fl%lref*qdt/fl%tref
1344
1345 {do ix^db = ixo^lim^db\}
1346 emin = rhonew(ix^d) * fl%tlow * rfactor(ix^d) * invgam
1347 lmax = max(zero, ( pnew(ix^d)*invgam - emin ) / qdt)
1348
1349 ! No cooling if temperature is below floor level.
1350 ! Assuming Bremsstrahlung if temperature is higher than maximum.
1351 if( te(ix^d)<= fl%tcoolmin) then
1352 l1 = zero
1353 else if( te(ix^d)>= fl%tcoolmax ) then
1354 call calc_l_extended(te(ix^d), l1, fl)
1355 l1 = l1 * rho(ix^d)**2
1356 l1 = min(l1, lmax)
1357 else
1358 call findy(te(ix^d), y1, fl)
1359 y2 = y1 + fact * rho(ix^d)*rc_gamma_1
1360 call findt(tlocal2, y2, fl)
1361 if( tlocal2 <= fl%tcoolmin ) then
1362 l1 = lmax
1363 else
1364 l1 = (te(ix^d)- tlocal2)*rho(ix^d)*rfactor(ix^d)*invgam/qdt
1365 end if
1366 l1 = min(l1, lmax)
1367 end if
1368 if(slab_uniform .and. fl%rad_cut .and. x(ix^d,ndim) .le. fl%rad_cut_hgt) then
1369 l1 = l1*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1370 end if
1371 coolrate(ix^d) = l1
1372 {end do\}
1373 end subroutine getvar_cooling_exact
1374
1375 subroutine radiative_cooling_add_source(qdt,ixI^L,ixO^L,wCT,wCTprim,w,x,&
1376 qsourcesplit,active,fl)
1377 ! w[iw]=w[iw]+qdt*S[wCT,x] where S is the source based on wCT within ixO
1378 use mod_global_parameters
1379 integer, intent(in) :: ixI^L, ixO^L
1380 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw), wCTprim(ixI^S,1:nw)
1381 double precision, intent(inout) :: w(ixI^S,1:nw)
1382 logical, intent(in) :: qsourcesplit
1383 logical, intent(inout) :: active
1384 type(rc_fluid), intent(in) :: fl
1385 double precision, allocatable, dimension(:^D&) :: Lequi
1386
1387 if(qsourcesplit .eqv.fl%rc_split) then
1388 active = .true.
1389 select case(fl%coolmethod)
1390 case ('explicit1')
1391 if(mype==0)then
1392 if(it==1) then
1393 write(*,*)'Fully explicit cooling is not completely safe in this version'
1394 write(*,*)'PROCEED WITH CAUTION!'
1395 endif
1396 endif
1397 call cool_explicit1(qdt,ixi^l,ixo^l,wct,w,x,fl)
1398 case ('explicit2')
1399 call cool_explicit2(qdt,ixi^l,ixo^l,wct,w,x,fl)
1400 case ('semiimplicit')
1401 call cool_semiimplicit(qdt,ixi^l,ixo^l,wct,w,x,fl)
1402 case ('implicit')
1403 call cool_implicit(qdt,ixi^l,ixo^l,wct,w,x,fl)
1404 case ('exact')
1405 call cool_exact(qdt,ixi^l,ixo^l,wct,wctprim,w,x,fl)
1406 case default
1407 call mpistop("This cooling method is unknown")
1408 end select
1409 if(fl%has_equi) then
1410 allocate(lequi(ixi^s))
1411 call get_cool_equi(qdt,ixi^l,ixo^l,wct,w,x,fl,lequi)
1412 w(ixo^s,fl%e_) = w(ixo^s,fl%e_)+lequi(ixo^s)
1413 deallocate(lequi)
1414 endif
1415 if( fl%Tfix ) call floortemperature(qdt,ixi^l,ixo^l,wct,w,x,fl)
1416 end if
1417 end subroutine radiative_cooling_add_source
1418
1419 subroutine floortemperature(qdt,ixI^L,ixO^L,wCT,w,x,fl)
1420 ! Force minimum temperature to a fixed temperature
1421 use mod_global_parameters
1422 integer, intent(in) :: ixI^L, ixO^L
1423 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw)
1424 double precision, intent(inout) :: w(ixI^S,1:nw)
1425 type(rc_fluid), intent(in) :: fl
1426 double precision :: etherm(ixI^S), rho(ixI^S), Rfactor(ixI^S),emin
1427 integer :: ix^D
1428
1429 call fl%get_pthermal(w,x,ixi^l,ixo^l,etherm)
1430 call fl%get_rho(w,x,ixi^l,ixo^l,rho)
1431 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1432 {do ix^db = ixo^lim^db\}
1433 emin = rho(ix^d)*fl%tlow*rfactor(ix^d)
1434 if(etherm(ix^d) < emin) then
1435 w(ix^d,fl%e_)=w(ix^d,fl%e_)+(emin-etherm(ix^d))*invgam
1436 end if
1437 {end do\}
1438 end subroutine floortemperature
1439
1440 subroutine get_cool_equi(qdt,ixI^L,ixO^L,wCT,w,x,fl,res)
1441 ! explicit cooling routine that depends on getdt to
1442 ! adjust the timestep. Accurate but incredibly slow
1443 use mod_global_parameters
1444
1445 integer, intent(in) :: ixI^L, ixO^L
1446 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw)
1447 double precision, intent(inout) :: w(ixI^S,1:nw)
1448 type(rc_fluid), intent(in) :: fl
1449 double precision, intent(out) :: res(ixI^S)
1450
1451 double precision :: pth(ixI^S),rho(ixI^S),Rfactor(ixI^S),L1,Tlocal2
1452 double precision :: Te(ixI^S)
1453 double precision :: emin, Lmax
1454 double precision :: Y1, Y2
1455 double precision :: de, emax,fact
1456 integer :: ix^D
1457
1458 call fl%get_pthermal_equi(wct,x,ixi^l,ixo^l,pth)
1459 call fl%get_rho_equi(wct,x,ixi^l,ixo^l,rho)
1460 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1461 te(ixo^s)=pth(ixo^s)/(rho(ixo^s)*rfactor(ixo^s))
1462
1463 res=0d0
1464
1465 if(fl%coolmethod == 'exact') then
1466
1467 fact = fl%lref*qdt/fl%tref
1468 {do ix^db = ixo^lim^db\}
1469 emin = rho(ix^d)*fl%tlow*rfactor(ix^d)*invgam
1470 lmax = max(zero,(pth(ix^d)*invgam-emin)/qdt)
1471 emax = max(zero, pth(ix^d)*invgam-emin)
1472 ! Determine explicit cooling
1473 ! If temperature is below floor level, no cooling.
1474 ! Stop wasting time and go to next gridpoint.
1475 ! If the temperature is higher than the maximum,
1476 ! assume Bremsstrahlung
1477 if( te(ix^d)<=fl%tcoolmin ) then
1478 l1 = zero
1479 else if( te(ix^d)>=fl%tcoolmax )then
1480 call calc_l_extended(te(ix^d), l1,fl)
1481 l1 = l1*rho(ix^d)**2
1482 if(phys_trac) then
1483 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1484 l1=l1*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1485 end if
1486 end if
1487 l1 = min(l1,lmax)
1488 res(ix^d) = l1*qdt
1489 else
1490 call findy(te(ix^d),y1,fl)
1491 y2 = y1 + fact * rho(ix^d)*rc_gamma_1
1492 call findt(tlocal2,y2,fl)
1493 if(tlocal2<=fl%tcoolmin) then
1494 de = emax
1495 else
1496 de = (te(ix^d)-tlocal2)*rho(ix^d)*rfactor(ix^d)*invgam
1497 end if
1498 if(phys_trac) then
1499 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1500 de=de*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1501 end if
1502 end if
1503 de = min(de,emax)
1504 res(ix^d) = de
1505 end if
1506 {end do\}
1507 else
1508 {do ix^db = ixo^lim^db\}
1509 emin = rho(ix^d)*fl%tlow*rfactor(ix^d)*invgam
1510 lmax = max(zero,pth(ix^d)*invgam-emin)/qdt
1511 ! Determine explicit cooling
1512 ! If temperature is below floor level, no cooling.
1513 ! Stop wasting time and go to next gridpoint.
1514 ! If the temperature is higher than the maximum,
1515 ! assume Bremsstrahlung
1516 if( te(ix^d)<=fl%tcoolmin ) then
1517 l1 = zero
1518 else if( te(ix^d)>=fl%tcoolmax )then
1519 call calc_l_extended(te(ix^d), l1,fl)
1520 else
1521 call findl(te(ix^d),l1,fl)
1522 end if
1523 l1 = l1*rho(ix^d)**2
1524 if(phys_trac) then
1525 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1526 l1=l1*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1527 end if
1528 end if
1529 l1 = min(l1,lmax)
1530 res(ix^d) =l1*qdt
1531 {end do\}
1532 end if
1533 end subroutine get_cool_equi
1534
1535 subroutine cool_explicit1(qdt,ixI^L,ixO^L,wCT,w,x,fl)
1536 ! explicit cooling routine that depends on getdt to
1537 ! adjust the timestep. Accurate but incredibly slow
1538 use mod_global_parameters
1539
1540 integer, intent(in) :: ixI^L, ixO^L
1541 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw)
1542 double precision, intent(inout) :: w(ixI^S,1:nw)
1543 type(rc_fluid), intent(in) :: fl
1544
1545 double precision :: L1,pth(ixI^S),pnew(ixI^S),rho(ixI^S),Rfactor(ixI^S)
1546 double precision :: Te(ixI^S)
1547 double precision :: emin, Lmax
1548 integer :: ix^D
1549
1550 call fl%get_pthermal(wct,x,ixi^l,ixo^l,pth)
1551 call fl%get_pthermal(w,x,ixi^l,ixo^l,pnew)
1552 call fl%get_rho(wct,x,ixi^l,ixo^l,rho)
1553 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1554 te(ixo^s)=pth(ixo^s)/(rho(ixo^s)*rfactor(ixo^s))
1555
1556 {do ix^db = ixo^lim^db\}
1557 emin = rho(ix^d)*fl%tlow*rfactor(ix^d)*invgam
1558 lmax = max(zero,pnew(ix^d)*invgam-emin)/qdt
1559 ! Determine explicit cooling
1560 ! If temperature is below floor level, no cooling.
1561 ! Stop wasting time and go to next gridpoint.
1562 ! If the temperature is higher than the maximum,
1563 ! assume Bremsstrahlung
1564 if( te(ix^d)<=fl%tcoolmin ) then
1565 l1 = zero
1566 else if( te(ix^d)>=fl%tcoolmax )then
1567 call calc_l_extended(te(ix^d), l1,fl)
1568 l1 = l1*rho(ix^d)**2
1569 if(phys_trac) then
1570 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1571 l1=l1*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1572 end if
1573 end if
1574 l1 = min(l1,lmax)
1575 else
1576 call findl(te(ix^d),l1,fl)
1577 l1 = l1*rho(ix^d)**2
1578 if(phys_trac) then
1579 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1580 l1=l1*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1581 end if
1582 end if
1583 l1 = min(l1,lmax)
1584 end if
1585 if(slab_uniform .and. fl%rad_cut .and. x(ix^d,ndim) .le. fl%rad_cut_hgt) then
1586 l1 = l1*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1587 end if
1588 w(ix^d,fl%e_) = w(ix^d,fl%e_)-l1*qdt
1589 {end do\}
1590 end subroutine cool_explicit1
1591
1592 subroutine cool_explicit2(qdt,ixI^L,ixO^L,wCT,w,x,fl)
1593 ! explicit cooling routine that does a series
1594 ! of small forward integration steps, to make
1595 ! sure the amount of cooling remains correct
1596 ! Not as accurate as 'explicit1', but a lot faster
1597 ! tends to overestimate cooling
1598 use mod_global_parameters
1599 integer, intent(in) :: ixI^L, ixO^L
1600 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw)
1601 double precision, intent(inout) :: w(ixI^S,1:nw)
1602 type(rc_fluid), intent(in) :: fl
1603 double precision :: Ltest, etherm, de
1604 double precision :: dtmax, dtstep
1605 double precision :: L1,pth(ixI^S),pnew(ixI^S),rho(ixI^S),Rfactor(ixI^S)
1606 double precision :: Tlocal1,plocal,Te(ixI^S)
1607 double precision :: emin, Lmax
1608 integer :: idt,ndtstep
1609 integer :: ix^D
1610
1611 call fl%get_pthermal(wct,x,ixi^l,ixo^l,pth)
1612 call fl%get_pthermal(w,x,ixi^l,ixo^l,pnew)
1613 call fl%get_rho(wct,x,ixi^l,ixo^l,rho)
1614 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1615 te(ixo^s)=pth(ixo^s)/(rho(ixo^s)*rfactor(ixo^s))
1616
1617 {do ix^db = ixo^lim^db\}
1618 ! Calculate explicit cooling value
1619 dtmax = qdt
1620 etherm = pth(ix^d)*invgam
1621 emin = rho(ix^d)*fl%tlow*rfactor(ix^d)*invgam
1622 lmax = max(zero,pnew(ix^d)*invgam-emin)/qdt
1623 ! Determine explicit cooling
1624 ! If temperature is below floor level, no cooling.
1625 ! Stop wasting time and go to next gridpoint.
1626 ! If the temperature is higher than the maximum,
1627 ! assume Bremmstrahlung
1628 if( te(ix^d)<=fl%tcoolmin ) then
1629 ltest = zero
1630 else if( te(ix^d)>=fl%tcoolmax )then
1631 call calc_l_extended(te(ix^d), ltest,fl)
1632 ltest = l1*rho(ix^d)**2
1633 if(phys_trac) then
1634 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1635 ltest=ltest*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1636 end if
1637 end if
1638 ltest = min(l1,lmax)
1639 if( dtmax>fl%cfrac*etherm/ltest) dtmax = fl%cfrac*etherm/ltest
1640 else
1641 call findl(te(ix^d),ltest,fl)
1642 ltest = ltest*rho(ix^d)**2
1643 if(phys_trac) then
1644 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1645 ltest=ltest*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1646 end if
1647 end if
1648 ltest = min(ltest,lmax)
1649 if( dtmax>fl%cfrac*etherm/ltest) dtmax = fl%cfrac*etherm/ltest
1650 end if
1651 ! Calculate number of steps for cooling
1652 ndtstep = max(nint(qdt/dtmax),1)+1
1653 dtstep = qdt/ndtstep
1654 ! Use explicit cooling value for first step
1655 de = ltest*dtstep
1656 etherm = etherm - de
1657
1658 do idt=2,ndtstep
1659 plocal = etherm*rc_gamma_1
1660 lmax = max(zero,etherm-emin)/dtstep
1661 ! Tlocal = P/(rho*R)
1662 tlocal1 = plocal/(rho(ix^d)*rfactor(ix^d))
1663 if( tlocal1<=fl%tcoolmin ) then
1664 l1 = zero
1665 exit
1666 else if( tlocal1>=fl%tcoolmax )then
1667 call calc_l_extended(tlocal1, l1,fl)
1668 l1 = l1*rho(ix^d)**2
1669 if(phys_trac) then
1670 if(tlocal1<block%wextra(ix^d,fl%Tcoff_)) then
1671 l1=l1*sqrt((tlocal1/block%wextra(ix^d,fl%Tcoff_))**5)
1672 end if
1673 end if
1674 l1 = min(l1,lmax)
1675 else
1676 call findl(tlocal1,l1,fl)
1677 l1 = l1*rho(ix^d)**2
1678 if(phys_trac) then
1679 if(tlocal1<block%wextra(ix^d,fl%Tcoff_)) then
1680 l1=l1*sqrt((tlocal1/block%wextra(ix^d,fl%Tcoff_))**5)
1681 end if
1682 end if
1683 l1 = min(l1,lmax)
1684 end if
1685 de = de + l1*dtstep
1686 etherm = etherm - l1*dtstep
1687 end do
1688 if(slab_uniform .and. fl%rad_cut .and. x(ix^d,ndim) .le. fl%rad_cut_hgt) then
1689 de = de*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1690 end if
1691 w(ix^d,fl%e_) = w(ix^d,fl%e_) -de
1692 {end do\}
1693 end subroutine cool_explicit2
1694
1695 subroutine cool_semiimplicit(qdt,ixI^L,ixO^L,wCT,w,x,fl)
1696 ! Semi-implicit cooling method based on a two point average
1697 ! Fast, but tends to underestimate cooling
1698 use mod_global_parameters
1699 integer, intent(in) :: ixI^L, ixO^L
1700 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw)
1701 double precision, intent(inout) :: w(ixI^S,1:nw)
1702 type(rc_fluid), intent(in) :: fl
1703 double precision :: L1,L2,Tlocal2
1704 double precision :: etemp
1705 double precision :: emin, Lmax
1706 double precision :: pth(ixI^S),pnew(ixI^S),rho(ixI^S),Rfactor(ixI^S),Te(ixI^S)
1707 integer :: ix^D
1708
1709 call fl%get_pthermal(wct,x,ixi^l,ixo^l,pth)
1710 call fl%get_pthermal(w,x,ixi^l,ixo^l,pnew)
1711 call fl%get_rho(wct,x,ixi^l,ixo^l,rho)
1712 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1713 te(ixo^s)=pth(ixo^s)/(rho(ixo^s)*rfactor(ixo^s))
1714
1715 {do ix^db = ixo^lim^db\}
1716 emin = rho(ix^d)*fl%tlow*rfactor(ix^d)*invgam
1717 lmax = max(zero,pnew(ix^d)*invgam-emin)/qdt
1718 ! Determine explicit cooling at present temperature
1719 !
1720 ! If temperature is below floor level, no cooling.
1721 ! Stop wasting time and go to next gridpoint.
1722 ! If the temperature is higher than the maximum,
1723 ! assume Bremsstrahlung
1724 if( te(ix^d)<=fl%tcoolmin ) then
1725 l1 = zero
1726 l2 = zero
1727 else
1728 if( te(ix^d)>=fl%tcoolmax ) then
1729 call calc_l_extended(te(ix^d), l1,fl)
1730 else
1731 call findl(te(ix^d),l1,fl)
1732 end if
1733 l1 = l1*rho(ix^d)**2
1734 if(phys_trac) then
1735 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1736 l1=l1*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1737 end if
1738 end if
1739 etemp = pth(ix^d)*invgam - l1*qdt
1740 tlocal2 = etemp*rc_gamma_1/(rho(ix^d)*rfactor(ix^d))
1741 ! Determine explicit cooling at new temperature
1742 if( tlocal2<=fl%tcoolmin ) then
1743 l2 = zero
1744 else if( tlocal2>=fl%tcoolmax )then
1745 call calc_l_extended(tlocal2, l2,fl)
1746 else
1747 call findl(tlocal2,l2,fl)
1748 end if
1749 l2 = l2*rho(ix^d)**2
1750 if(phys_trac) then
1751 if(tlocal2<block%wextra(ix^d,fl%Tcoff_)) then
1752 l2=l2*sqrt((tlocal2/block%wextra(ix^d,fl%Tcoff_))**5)
1753 end if
1754 end if
1755 if(slab_uniform .and. fl%rad_cut .and. x(ix^d,ndim) .le. fl%rad_cut_hgt) then
1756 l1 = l1*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1757 l2 = l2*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1758 end if
1759 w(ix^d,fl%e_) = w(ix^d,fl%e_) - min(half*(l1+l2),lmax)*qdt
1760 end if
1761 {end do\}
1762 end subroutine cool_semiimplicit
1763
1764 subroutine cool_implicit(qdt,ixI^L,ixO^L,wCT,w,x,fl)
1765 ! Implicit cooling method based on a half-step
1766 ! refinement algorithm
1767 use mod_global_parameters
1768 integer, intent(in) :: ixI^L, ixO^L
1769 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw)
1770 double precision, intent(inout) :: w(ixI^S,1:nw)
1771 type(rc_fluid), intent(in) :: fl
1772 double precision :: Ltemp,Tnew,f1,f2,pth(ixI^S), pnew(ixI^S), rho(ixI^S), Rfactor(ixI^S)
1773 double precision :: elocal, Te(ixI^S)
1774 double precision :: emin, Lmax, eold, enew, estep
1775 double precision, parameter :: e_error = 1.0d-6
1776 integer, parameter :: maxiter = 100
1777 integer :: ix^D, j
1778
1779 call fl%get_pthermal(wct,x,ixi^l,ixo^l,pth)
1780 call fl%get_pthermal(w,x,ixi^l,ixo^l,pnew)
1781 call fl%get_rho(wct,x,ixi^l,ixo^l,rho)
1782 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1783 te(ixo^s)=pth(ixo^s)/(rho(ixo^s)*rfactor(ixo^s))
1784
1785 {do ix^db = ixo^lim^db\}
1786 elocal = pth(ix^d)*invgam
1787 emin = rho(ix^d)*fl%tlow*rfactor(ix^d)*invgam
1788 lmax = max(zero,pnew(ix^d)*invgam-emin)/qdt
1789 ! Determine explicit cooling at present temperature
1790 ! If temperature is below floor level, no cooling.
1791 ! Stop wasting time and go to next gridpoint.
1792 ! If the temperature is higher than the maximum,
1793 ! assume Bremsstrahlung
1794 if( te(ix^d)<=fl%tcoolmin ) then
1795 ltemp = zero
1796 else
1797 eold = elocal
1798 enew = elocal
1799 estep = -(smalldouble)
1800 f2 = 1.d0
1801 do j=1,maxiter+1
1802 if( j>maxiter ) call mpistop("Implicit cooling exceeds maximum iterations")
1803 tnew = enew*rc_gamma_1/(rho(ix^d)*rfactor(ix^d))
1804 if( tnew<=fl%tcoolmin ) then
1805 ltemp = lmax
1806 exit
1807 else if( tnew>=fl%tcoolmax )then
1808 call calc_l_extended(tnew, ltemp,fl)
1809 else
1810 call findl(tnew,ltemp,fl)
1811 end if
1812 ltemp = ltemp*rho(ix^d)**2
1813 eold = enew + ltemp*qdt
1814 f1 = elocal -eold
1815 if(abs(half*f1/(elocal+eold)) < e_error) exit
1816 if(phys_trac) then
1817 if(tnew<block%wextra(ix^d,fl%Tcoff_)) then
1818 ltemp=ltemp*sqrt((tnew/block%wextra(ix^d,fl%Tcoff_))**5)
1819 end if
1820 end if
1821 if(j==1) estep = max((elocal-emin)*half,smalldouble)
1822 if(f1*f2 < zero) estep = -half*estep
1823 f2 = f1
1824 enew = enew +estep
1825 end do
1826 end if
1827 if(slab_uniform .and. fl%rad_cut .and. x(ix^d,ndim) .le. fl%rad_cut_hgt) then
1828 ltemp = ltemp*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1829 end if
1830 w(ix^d,fl%e_) = w(ix^d,fl%e_) - min(ltemp,lmax)*qdt
1831 {end do\}
1832 end subroutine cool_implicit
1833
1834 subroutine cool_exact(qdt,ixI^L,ixO^L,wCT,wCTprim,w,x,fl)
1835 ! Cooling routine using exact integration method from Townsend 2009
1836 use mod_global_parameters
1837 integer, intent(in) :: ixI^L, ixO^L
1838 double precision, intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw), wCTprim(ixI^S,1:nw)
1839 double precision, intent(inout) :: w(ixI^S,1:nw)
1840 type(rc_fluid), intent(in) :: fl
1841 double precision :: Y1, Y2
1842 double precision :: L1, pth(ixI^S), Tlocal2, pnew(ixI^S)
1843 double precision :: rho(ixI^S), Te(ixI^S), rhonew(ixI^S), Rfactor(ixI^S)
1844 double precision :: emin, Lmax, fact
1845 double precision :: de, emax
1846 integer :: ix^D
1847
1848 call fl%get_rho(wct,x,ixi^l,ixo^l,rho)
1849 call fl%get_var_Rfactor(wct,x,ixi^l,ixo^l,rfactor)
1850 if(phys_equi_pe) then
1851 ! need pressure splitting
1852 call fl%get_pthermal(wct,x,ixi^l,ixo^l,te)
1853 te(ixo^s)=te(ixo^s)/(rho(ixo^s)*rfactor(ixo^s))
1854 else
1855 te(ixo^s)=wctprim(ixo^s,iw_e)/(rho(ixo^s)*rfactor(ixo^s))
1856 end if
1857 call fl%get_pthermal(w,x,ixi^l,ixo^l,pnew)
1858 call fl%get_rho(w,x,ixi^l,ixo^l,rhonew)
1859
1860 fact = fl%lref*qdt/fl%tref
1861
1862 {do ix^db = ixo^lim^db\}
1863 emin = rhonew(ix^d)*fl%tlow*rfactor(ix^d)*invgam
1864 lmax = max(zero,pnew(ix^d)*invgam-emin)/qdt
1865 emax = max(zero,pnew(ix^d)*invgam-emin)
1866 ! Determine explicit cooling
1867 ! If temperature is below floor level, no cooling.
1868 ! Stop wasting time and go to next gridpoint.
1869 ! If the temperature is higher than the maximum,
1870 ! assume Bremsstrahlung
1871 if( te(ix^d)<=fl%tcoolmin ) then
1872 l1 = zero
1873 else if( te(ix^d)>=fl%tcoolmax )then
1874 call calc_l_extended(te(ix^d), l1,fl)
1875 l1 = l1*rho(ix^d)**2
1876 if(phys_trac) then
1877 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1878 l1=l1*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1879 end if
1880 end if
1881 l1 = min(l1,lmax)
1882 if(slab_uniform .and. fl%rad_cut .and. x(ix^d,ndim) .le. fl%rad_cut_hgt) then
1883 l1 = l1*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1884 end if
1885 w(ix^d,fl%e_) = w(ix^d,fl%e_)-l1*qdt
1886 else
1887 call findy(te(ix^d),y1,fl)
1888 y2 = y1 + fact*rho(ix^d)*rc_gamma_1
1889 call findt(tlocal2,y2,fl)
1890 if(tlocal2<=fl%tcoolmin) then
1891 de = emax
1892 else
1893 de = (te(ix^d)-tlocal2)*rho(ix^d)*rfactor(ix^d)*invgam
1894 end if
1895 if(phys_trac) then
1896 if(te(ix^d)<block%wextra(ix^d,fl%Tcoff_)) then
1897 de=de*sqrt((te(ix^d)/block%wextra(ix^d,fl%Tcoff_))**5)
1898 end if
1899 end if
1900 de = min(de,emax)
1901 if(slab_uniform .and. fl%rad_cut .and. x(ix^d,ndim) .le. fl%rad_cut_hgt) then
1902 de = de*exp(-(x(ix^d,ndim)-fl%rad_cut_hgt)**2/fl%rad_cut_dey**2)
1903 end if
1904 w(ix^d,fl%e_) = w(ix^d,fl%e_)-de
1905 end if
1906 {end do\}
1907 end subroutine cool_exact
1908
1909 subroutine calc_l_extended (tpoint, lpoint,fl)
1910 ! Calculate l for t beyond tcoolmax
1911 ! Assumes Bremsstrahlung for the interpolated tables
1912 ! Uses the power law for piecewise power laws
1913 double precision, intent(IN) :: tpoint
1914 double precision, intent(OUT) :: lpoint
1915 type(rc_fluid), intent(in) :: fl
1916
1917 if(fl%isPPL) then
1918 lpoint =fl%l_PPL(fl%n_PPL) * ( tpoint / fl%t_PPL(fl%n_PPL) )**fl%a_PPL(fl%n_PPL)
1919 else
1920 lpoint = fl%Lcool(fl%ncool) * sqrt( tpoint / fl%tcoolmax)
1921 end if
1922 end subroutine calc_l_extended
1923
1924 subroutine findl (tpoint,Lpoint,fl)
1925 ! Fast search option to find correct point
1926 ! in cooling curve
1927 use mod_global_parameters
1928
1929 double precision,intent(IN) :: tpoint
1930 double precision, intent(OUT) :: Lpoint
1931 type(rc_fluid), intent(in) :: fl
1932
1933 double precision :: lgtp
1934 integer :: jl,jc,jh,i
1935
1936 if(fl%isPPL) then
1937 i = maxloc(fl%t_PPL, dim=1, mask=fl%t_PPL<tpoint)
1938 lpoint = fl%l_PPL(i) * (tpoint / fl%t_PPL(i))**fl%a_PPL(i)
1939 else
1940 lgtp = dlog10(tpoint)
1941 jl = int((lgtp - fl%lgtcoolmin) /fl%lgstep) + 1
1942 lpoint = fl%Lcool(jl)+ (tpoint-fl%tcool(jl)) &
1943 * (fl%Lcool(jl+1)-fl%Lcool(jl)) &
1944 / (fl%tcool(jl+1)-fl%tcool(jl))
1945 end if
1946
1947! if (tpoint == fl%tcoolmin) then
1948! Lpoint = fl%Lcool(1)
1949! else if (tpoint == fl%tcoolmax) then
1950! Lpoint = fl%Lcool(fl%ncool)
1951! else
1952! jl=0
1953! jh=fl%ncool+1
1954! do
1955! if (jh-jl <= 1) exit
1956! jc=(jh+jl)/2
1957! if (tpoint >= fl%tcool(jc)) then
1958! jl=jc
1959! else
1960! jh=jc
1961! end if
1962! end do
1963! ! Linear interpolation to obtain correct cooling
1964! Lpoint = fl%Lcool(jl)+ (tpoint-fl%tcool(jl)) &
1965! * (fl%Lcool(jl+1)-fl%Lcool(jl)) &
1966! / (fl%tcool(jl+1)-fl%tcool(jl))
1967! end if
1968 end subroutine findl
1969
1970 subroutine findy (tpoint,Ypoint,fl)
1971 ! Fast search option to find correct point in cooling time (TEF)
1972 use mod_global_parameters
1973
1974 double precision,intent(IN) :: tpoint
1975 double precision, intent(OUT) :: Ypoint
1976 type(rc_fluid), intent(in) :: fl
1977
1978 double precision :: lgtp
1979 double precision :: y_extra,factor
1980 integer :: jl,jc,jh,i
1981
1982 if(fl%isPPL) then
1983 i = maxloc(fl%t_PPL, dim=1, mask=fl%t_PPL<tpoint)
1984 factor = fl%l_PPL(fl%n_PPL+1) * fl%t_PPL(i) / (fl%l_PPL(i) * fl%t_PPL(fl%n_PPL+1))
1985 if(fl%a_PPL(i)==1.d0) then
1986 y_extra = log( fl%t_PPL(i) / tpoint )
1987 else
1988 y_extra = 1 / (1 - fl%a_PPL(i)) * (1 - ( fl%t_PPL(i) / tpoint )**(fl%a_PPL(i)-1) )
1989 end if
1990 ypoint = fl%y_PPL(i) + factor*y_extra
1991 else
1992 lgtp = dlog10(tpoint)
1993 jl = int((lgtp - fl%lgtcoolmin) / fl%lgstep) + 1
1994 ypoint = fl%Yc(jl)+ (tpoint-fl%tcool(jl)) &
1995 * (fl%Yc(jl+1)-fl%Yc(jl)) &
1996 / (fl%tcool(jl+1)-fl%tcool(jl))
1997 end if
1998
1999 ! integer i
2000 !
2001 ! if (tpoint == tcoolmin) then
2002 ! Ypoint = Yc(1)
2003 ! else if (tpoint == tcoolmax) then
2004 ! Ypoint = Yc(ncool)
2005 ! else
2006 ! jl=0
2007 ! jh=ncool+1
2008 ! do
2009 ! if (jh-jl <= 1) exit
2010 ! jc=(jh+jl)/2
2011 ! if (tpoint >= tcool(jc)) then
2012 ! jl=jc
2013 ! else
2014 ! jh=jc
2015 ! end if
2016 ! end do
2017 ! ! Linear interpolation to obtain correct value
2018 ! Ypoint = Yc(jl)+ (tpoint-tcool(jl)) &
2019 ! * (Yc(jl+1)-Yc(jl)) &
2020 ! / (tcool(jl+1)-tcool(jl))
2021 ! end if
2022 end subroutine findy
2023
2024 subroutine findt (tpoint,Ypoint,fl)
2025 ! Fast search option to find correct temperature
2026 ! from temporal evolution function. Only possible this way because T is a monotonously
2027 ! decreasing function for the interpolated tables
2028 ! Uses eq. A7 from Townsend 2009 for piecewise power laws
2029 use mod_global_parameters
2030
2031 double precision,intent(OUT) :: tpoint
2032 double precision, intent(IN) :: Ypoint
2033 type(rc_fluid), intent(in) :: fl
2034
2035 double precision :: factor
2036 integer :: jl,jc,jh,i
2037
2038 if(fl%isPPL) then
2039 i = minloc(fl%y_PPL, dim=1, mask=fl%y_PPL>ypoint)
2040 factor = fl%l_PPL(i) * fl%t_PPL(fl%n_PPL+1) / (fl%l_PPL(fl%n_PPL+1) * fl%t_PPL(i))
2041 if(fl%a_PPL(i)==1.d0) then
2042 tpoint = fl%t_PPL(i) * exp( -1.d0 * factor * ( ypoint - fl%y_PPL(i)))
2043 else
2044 tpoint = fl%t_PPL(i) * (1 - (1 - fl%a_PPL(i)) * factor * (ypoint - fl%y_PPL(i)))**(1 / (1 - fl%a_PPL(i)))
2045 end if
2046 else
2047 if(ypoint >= fl%Yc(1)) then
2048 tpoint = fl%tcoolmin
2049 else if (ypoint == fl%Yc(fl%ncool)) then
2050 tpoint = fl%tcoolmax
2051 else
2052 jl=0
2053 jh=fl%ncool+1
2054 do
2055 if(jh-jl <= 1) exit
2056 jc=(jh+jl)/2
2057 if(ypoint <= fl%Yc(jc)) then
2058 jl=jc
2059 else
2060 jh=jc
2061 end if
2062 end do
2063 ! Linear interpolation to obtain correct temperature
2064 tpoint = fl%tcool(jl)+ (ypoint-fl%Yc(jl)) &
2065 * (fl%tcool(jl+1)-fl%tcool(jl)) &
2066 / (fl%Yc(jl+1)-fl%Yc(jl))
2067 end if
2068 end if
2069 end subroutine findt
2070
2071 subroutine finddldt (tpoint,dLpoint,fl)
2072 ! Fast search option to find correct point
2073 ! in derivative of cooling curve
2074 ! Does not work for the piecewise power laws
2075 use mod_global_parameters
2076
2077 double precision,intent(IN) :: tpoint
2078 double precision, intent(OUT) :: dLpoint
2079 type(rc_fluid), intent(in) :: fl
2080
2081 double precision :: lgtp
2082 integer :: jl,jc,jh
2083
2084 lgtp = dlog10(tpoint)
2085 jl = int((lgtp -fl%lgtcoolmin) / fl%lgstep) + 1
2086 dlpoint = fl%dLdtcool(jl)+ (tpoint-fl%tcool(jl)) &
2087 * (fl%dLdtcool(jl+1)-fl%dLdtcool(jl)) &
2088 / (fl%tcool(jl+1)-fl%tcool(jl))
2089
2090! if (tpoint == tcoolmin) then
2091! dLpoint = dLdtcool(1)
2092! else if (tpoint == tcoolmax) then
2093! dLpoint = dLdtcool(ncool)
2094! else
2095! jl=0
2096! jh=ncool+1
2097! do
2098! if (jh-jl <= 1) exit
2099! jc=(jh+jl)/2
2100! if (tpoint >= tcool(jc)) then
2101! jl=jc
2102! else
2103! jh=jc
2104! end if
2105! end do
2106! ! Linear interpolation to obtain correct cooling derivative
2107! dLpoint = dLdtcool(jl)+ (tpoint-tcool(jl)) &
2108! * (dLdtcool(jl+1)-dLdtcool(jl)) &
2109! / (tcool(jl+1)-tcool(jl))
2110! end if
2111 end subroutine finddldt
2112end 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:774
mod_global_parameters::unitpar
integer, parameter unitpar
file handle for IO
Definition mod_global_parameters.t:372
mod_global_parameters::it
integer it
Number of time steps taken.
Definition mod_global_parameters.t:438
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:768
mod_global_parameters::mype
integer mype
The rank of the current MPI task.
Definition mod_global_parameters.t:223
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:664
mod_global_parameters::slab_uniform
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
Definition mod_global_parameters.t:568
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:1765
mod_radiative_cooling::cool_explicit1
subroutine cool_explicit1(qdt, ixil, ixol, wct, w, x, fl)
Definition mod_radiative_cooling.t:1536
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:1925
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:1441
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::finddldt
subroutine finddldt(tpoint, dlpoint, fl)
Definition mod_radiative_cooling.t:2072
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:1377
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:1910
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:1593
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:1835
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:2025
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:1420
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:1696
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:1971
mod_radiative_cooling::rc_fluid
Definition mod_radiative_cooling.t:49
Generated by doxygen 1.9.8