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