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