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