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