DISPATCH
cooling_Gnedin.f90
1 ! cooling_Gnedin.f90 $Id: 87e885e42c6cff1bc93bae7e5bf8084992e568ec $
2 !-------------------------------------------------------------------------------
3 FUNCTION cooling_gnedin (self, temp)
4  implicit none
5  class(cooling_t):: self
6  real(kind=dp) cooling_gnedin, temp, p
7  integer i
8  real, parameter:: dlgt=0.1, lgt0=1.0
9  integer, parameter:: ntab=81
10  real(kind=dp), parameter, dimension(ntab):: ln25 = (/ &
11  -8.6450,-8.1238,-7.6490,-7.2158,-6.8207,-6.4650,-6.1428,-5.8379, &
12  -5.5115,-5.1155,-4.6465,-4.1402,-3.6593,-3.2304,-2.8680,-2.5690, &
13  -2.3228,-2.1203,-1.9505,-1.8042,-1.6724,-1.5460,-1.2870,-1.0269, &
14  -0.7525,-0.4544,-0.1468, 0.1931, 0.5939, 1.1394, 2.0933, 3.9022, &
15  6.1748, 6.5501, 6.2719, 6.2451, 6.4187, 6.7583, 7.2167, 7.6392, &
16  7.7969, 7.7480, 7.7677, 7.8474, 7.8755, 7.4928, 7.2349, 7.1854, &
17  6.9207, 6.6918, 6.6560, 6.5481, 6.3418, 5.9320, 5.4885, 5.1625, &
18  4.9726, 4.9265, 4.9445, 4.9822, 4.9843, 4.9280, 4.8032, 4.6959, &
19  4.6625, 4.6812, 4.7414, 4.8195, 4.9030, 5.0019, 5.0907, 5.1745, &
20  5.2423, 5.2943, 5.3318, 5.3552, 5.3790, 5.3932, 5.4009, 5.4076, &
21  5.4205 /)
22 !...............................................................................
23  if (temp < t_mc) then
24  cooling_gnedin=0.0
25  else
26  p=1.+(log10(temp)-lgt0)/dlgt
27  i=min(max(int(p),1),ntab-1)
28  p=p-i
29  cooling_gnedin=1e-25*exp((1-p)*ln25(i)+p*ln25(i+1))
30  endif
31 END