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