-
Notifications
You must be signed in to change notification settings - Fork 10
/
Copy pathplanck.F90
327 lines (260 loc) · 10.9 KB
/
planck.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
! Include shortname defintions, so that the F77 code does not have to be modified to
! reference the CARMA structure.
#include "carma_globaer.h"
module planck
contains
!! This routine calculates the planck intensity.
!!
!! This algorithm is based upon eqn 1.2.4 from Liou[2002].
!!
!! @author Chuck Bardeen
!! @version Jan-2010
function planckIntensity(wvl, temp)
! types
use carma_precision_mod
use carma_enums_mod
use carma_constants_mod
use carma_types_mod
use carmastate_mod
use carma_mod
implicit none
real(kind=f), intent(in) :: wvl !! wavelength (cm)
real(kind=f), intent(in) :: temp !! temperature (K)
real(kind=f) :: planckIntensity !! Planck intensity (erg/s/cm2/sr/cm)
! Local declarations
real(kind=f), parameter :: C = 2.9979e10_f ! Speed of light [cm/s]
real(kind=f), parameter :: H = 6.62608e-27_f ! Planck constant [erg s]
! Calculate the planck intensity.
planckIntensity = 2._f * H * C**2 / ((wvl**5) * (exp(H * C / (BK * wvl * temp)) - 1._f))
! Return the planck intensity to the caller.
return
end function
!! This routine calculates the total planck intensity from the specified
!! wavelength to a wavelength of 0.
!!
!! This algorithm is based upon Widger and Woodall, BAMS, 1976 as
!! indicated at http://www.spectralcalc.com/blackbody/appendixA.html.
!!
!! @author Chuck Bardeen
!! @version Aug-2011
function planckIntensityWidger1976(wvl, temp, miniter)
! types
use carma_precision_mod
use carma_enums_mod
use carma_constants_mod
use carma_types_mod
use carmastate_mod
use carma_mod
implicit none
real(kind=f), intent(in) :: wvl !! band center wavelength (cm)
real(kind=f), intent(in) :: temp !! temperature (K)
integer, intent(in) :: miniter !! minimum iterations
real(kind=f) :: planckIntensityWidger1976 !! Planck intensity (erg/s/cm2/sr/cm)
! Local Variables
real(kind=f), parameter :: C = 299792458.0_f ! Speed of light [m/s]
real(kind=f), parameter :: H = 6.6260693e-34_f ! Planck constant [J s]
real(kind=f), parameter :: BZ = 1.380658e-23_f ! Boltzman constant
real(kind=f) :: c1, x, x2, x3, sumJ, dn, sigma
integer :: iter, n
sigma = 1._f / wvl
c1 = H * C / BZ
x = c1 * 100._f * sigma / temp
x2 = x * x
x3 = x2 * x
! Use fewer iterations, since speed is more important than accuracy for
! the particle heating code, and even with fewer iterations the results
! with CAM bands still show good accuracy.
! iter = min(512, int(2._f + 20._f / x))
iter = min(miniter, int(2._f + 20._f / x))
sumJ = 0._f
do n = 1, iter
dn = 1._f / n
sumJ = sumJ + exp(-n*x) * (x3 + (3.0_f * x2 + 6.0_f * (x + dn) * dn) * dn) * dn
end do
! Convert results from W/m2/sr to erg/cm2/s/sr
planckIntensityWidger1976 = 2.0_f * H * (C**2) * ((temp / c1) ** 4) * sumJ * 1e7_f / 1e4_f
return
end function
!! This routine calculates the average planck intensity in the wavelength
!! band defined by wvl and dwvl.
!!
!! This algorithm is based upon Widger and Woodall, BAMS, 1976 as
!! indicated at http://www.spectralcalc.com/blackbody/appendixA.html.
!!
!! @author Chuck Bardeen
!! @version Aug-2011
function planckBandIntensityWidger1976(wvl, dwvl, temp, miniter)
! types
use carma_precision_mod
use carma_enums_mod
use carma_constants_mod
use carma_types_mod
use carmastate_mod
use carma_mod
implicit none
real(kind=f), intent(in) :: wvl !! band center wavelength (cm)
real(kind=f), intent(in) :: dwvl !! band width (cm)
real(kind=f), intent(in) :: temp !! temperature (K)
integer, intent(in) :: miniter !! minimum iterations
real(kind=f) :: planckBandIntensityWidger1976 !! Planck intensity (erg/s/cm2/sr/cm)
! Calculate the integral from the edges to 0 and subtract.
planckBandIntensityWidger1976 = &
(planckIntensityWidger1976(wvl + (dwvl / 2._f), temp, miniter) &
- planckIntensityWidger1976(wvl - (dwvl / 2._f), temp, miniter)) / dwvl
return
end function
!! This routine calculates the average planck intensity in the wavelength
!! band defined by wvl and dwvl.
!!
!! This algorithm does a brute force integral by dividing the band into
!! small sub-bands. This routine can be slow.
!!
!! @author Chuck Bardeen
!! @version Aug-2011
function planckBandIntensity(wvl, dwvl, temp, iter)
! types
use carma_precision_mod
use carma_enums_mod
use carma_constants_mod
use carma_types_mod
use carmastate_mod
use carma_mod
implicit none
real(kind=f), intent(in) :: wvl !! band center wavelength (cm)
real(kind=f), intent(in) :: dwvl !! band width (cm)
real(kind=f), intent(in) :: temp !! temperature (K)
integer, intent(in) :: iter !! number of iterations
real(kind=f) :: planckBandIntensity !! Planck intensity (erg/s/cm2/sr/cm)
! Local Variables
real(kind=f) :: wstart ! Starting wavelength (cm)
real(kind=f) :: ddwave ! sub-band width (cm)
integer :: i
wstart = wvl - (dwvl / 2._f)
ddwave = dwvl / iter
planckBandIntensity = 0._f
do i = 1, iter
planckBandIntensity = planckBandIntensity + planckIntensity(wstart + (i - 0.5_f) * ddwave, temp) * ddwave
end do
planckBandIntensity = planckBandIntensity / dwvl
return
end function
!! This routine calculates the average planck intensity in the wavelength
!! band defined by wvl and dwvl.
!!
!! error computed on full spectrum compared to planck function. Band-levels may be different
!! 8.9% error with 5 quadrature points in [100 micrometer, 1 millimeter]
!! 1.7% error with 10 quadrature points in [100 micrometer, 1 millimeter]
!! 0.001% error with 100 quadrature points in [100 micrometer, 1 millimeter]
!!
!! NOTE: This code was design to work with the CAM RRTMG band structure, it may not work as
!! well with arbitrary bands.
!!
!! NOTE: For most RRTMG bands, 3 quadrature points are probably sufficient, but testing is
!! left to the reader.
!!
!! @author Andrew Conley, Chuck Bardeen
!! @version Aug-2011
function planckBandIntensityConley2011(wvl, dwvl, temp, iter)
! types
use carma_precision_mod
use carma_enums_mod
use carma_constants_mod
use carma_types_mod
use carmastate_mod
use carma_mod
implicit none
real(kind=f), intent(in) :: wvl !! band center wavelength (cm)
real(kind=f), intent(in) :: dwvl !! band width (cm)
real(kind=f), intent(in) :: temp !! temperature (K)
integer, intent(in) :: iter !! number of iterations
real(kind=f) :: planckBandIntensityConley2011 !! Planck intensity (erg/s/cm2/sr/cm)
real(kind=f) :: half = 0.5_f
real(kind=f) :: third= 1._f / 3._f
real(kind=f) :: sixth= 1._f / 6._f
real(kind=f) :: tfth = 1._f /24._f
real(kind=f) :: k = 1.3806488e-23_f ! boltzmann J/K
real(kind=f) :: c = 2.99792458e8_f ! light m/s
real(kind=f) :: h = 6.62606957e-34_f ! planck J s
real(kind=f) :: sigma = 5.670373e-8_f ! stef-bolt W/m/m/k/k/k/k
real(kind=f) :: lambda1 ! wavelength m (lower bound)
real(kind=f) :: lambda2 ! wavelength m (upper bound)
! quadrature iteration
integer :: i,inumber
! internal temporary variables
real(kind=f) :: fr1, fr2 ! frequency bounds of partition
real(kind=f) :: kt ! k_boltzmann * temperature
real(kind=f) :: l1,l2 ! lower and upper bounds of (wavelength)
real(kind=f) :: dellam ! fraction multiplier for next lambda interval
real(kind=f) :: t1,t3 ! 2nd and 4th order terms
real(kind=f) :: total, total2 ! 2nd and 4th order cumulative partial integral
real(kind=f) :: e,d,em1i,di,ci ! exponential terms appearing in integral
real(kind=f) :: dfr,m,a,o,tt,mi ! terms appearing in integral
real(kind=f) :: argexp ! argument to exponent
real(kind=f) :: coeff ! front coefficient of integral
real(kind=f) :: planck ! planck function
inumber = iter ! number of partitions
!initialize
total = 0._f ! partial (cumulative) integral (4th order)
! total2 = 0._f ! partial (cumulative) integral (2nd order)
kt = k*temp
lambda1 = (wvl - (dwvl / 2._f)) * 1e-2_f
lambda2 = (wvl + (dwvl / 2._f)) * 1e-2_f
ci = 1._f/c
if (inumber .gt. 1) then
l1 = lambda1
dellam = exp(log(lambda1/lambda2)/inumber)
l2 = l1/dellam
fr1 = c/l2
fr2 = c/l1
else
dellam = 1._f ! meaningless
l1 = lambda1
l2 = lambda2
fr1 = c/l2
fr2 = c/l1
endif
! accumulate integral by stepping (backwards) through partions of frequency
do i = 1,inumber
! constants
dfr = half * (fr2-fr1) ! half-range freq interval
m = half * (fr1+fr2) ! mean freq
mi = 1._f/m
a = h/kt ! alpha
argexp = a*m
if (argexp .lt. 0.5_f) then
e = 1._f + &
argexp + &
(argexp*argexp)*half + &
(argexp*argexp*argexp)*sixth + &
(argexp*argexp*argexp*argexp)*tfth
em1i = 1._f/(e - 1._f )
di = e*em1i
else if (argexp .lt. 20.0_f) then
e = exp(argexp)
em1i = 1._f/(e - 1._f )
di = e*em1i
else
e = 1.e+20_f ! exp(20) is large. Use this for frequency >> Temperature
em1i = 1.e-20_f
di = 1._f
endif
! frontpiece
coeff = 2._f*h*m*m*m*ci*ci*em1i
! integrals
o = fr2-fr1 ! int 1 deps
tt = 2._f*(dfr*dfr*dfr)*third ! int eps^2 deps
! term and 4th order correction
t1 = 1._f
t3 = 3._f*mi*mi - 3._f*a*di*mi + a*a*di*di - half*a*a*di
! t3 could be made more stable by placing (-) terms in denominator of pade approx.
! sum it up. Total is 4th order, total2 is 2nd order
total = total + coeff*(o*t1+tt*t3)
! total2 = total2 + coeff*o*t1
fr2 = fr1
fr1 = fr1 * dellam
enddo
! Convert to erg/cm2/s/sr/cm
planckBandIntensityConley2011 = total * 1e7_f / 1e4_f / dwvl
return
end function
end module planck