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 : !! There are several different algorithms that can be used to solve
6 : !! a mie calculation for the optical properties of particles. This
7 : !! routine provides a generic front end to these different mie
8 : !! routines.
9 : !!
10 : !! Current methods are:
11 : !! miess - Original CARMA code, from Toon and Ackerman, supports core/shell
12 : !! bhmie - Homogeneous sphere, from Bohren and Huffman, handles wider range of parameters
13 : !!
14 : !! @author Chuck Bardeen
15 : !! @version 2011
16 0 : subroutine mie(carma, miertn, radius, wavelength, nmonomer, fractaldim, rmonomer, falpha_in, m, rcore, mcore, lqext, lqsca, lasym, rc)
17 :
18 : ! types
19 : use carma_precision_mod
20 : use carma_enums_mod
21 : use carma_constants_mod
22 : use carma_types_mod
23 : use carma_mod
24 : use fractal_meanfield_mod
25 :
26 : implicit none
27 :
28 : type(carma_type), intent(in) :: carma !! the carma object
29 : integer, intent(in) :: miertn !! mie routine enumeration
30 : real(kind=f), intent(in) :: radius !! radius (cm)
31 : real(kind=f), intent(in) :: wavelength !! wavelength (cm)
32 : real(kind=f), intent(in) :: nmonomer !! number of monomers per aggregate [fractal particles only]
33 : real(kind=f), intent(in) :: fractaldim !! fractal dimension [fractal particles only]
34 : real(kind=f), intent(in) :: rmonomer !! monomer size (units?) [fractal particles only]
35 : real(kind=f), intent(in) :: falpha_in !! packing coefficient [fractal particles only]
36 : complex(kind=f), intent(in) :: m !! refractive index particle
37 : real(kind=f), intent(in) :: rcore !! radius core (cm) [core shell only]
38 : complex(kind=f), intent(in) :: mcore !! refractive index core [core shell only]
39 : real(kind=f), intent(out) :: lqext !! EFFICIENCY FACTOR FOR EXTINCTION
40 : real(kind=f), intent(out) :: lqsca !! EFFICIENCY FACTOR FOR SCATTERING
41 : real(kind=f), intent(out) :: lasym !! asymmetry factor
42 : integer, intent(inout) :: rc !! return code, negative indicates failure
43 :
44 :
45 : integer, parameter :: nang = 10 ! Number of angles
46 :
47 : real(kind=f) :: theta(IT)
48 : real(kind=f) :: wvno
49 : real(kind=f) :: rfr
50 : real(kind=f) :: rfi
51 : real(kind=f) :: rfcr
52 : real(kind=f) :: rfci
53 : real(kind=f) :: x
54 : real(kind=f) :: qback
55 : real(kind=f) :: ctbrqs
56 : complex(kind=f) :: s1(2*nang-1)
57 : complex(kind=f) :: s2(2*nang-1)
58 : real(kind=f) :: rmonomer_out
59 : real(kind=f) :: fractaldim_out
60 :
61 : ! Calculate the wave number.
62 0 : wvno = 2._f * PI / wavelength
63 :
64 : ! Select the appropriate routine.
65 0 : if (miertn == I_MIERTN_TOON1981) then
66 :
67 : ! We only care about the forward direction.
68 0 : theta(:) = 0.0_f
69 :
70 0 : rfr = real(m)
71 0 : rfi = aimag(m)
72 :
73 0 : if (rcore .gt. 0.0_f) then
74 0 : rfcr = real(mcore)
75 0 : rfci = aimag(mcore)
76 : else
77 0 : rfcr = rfr
78 0 : rfci = rfi
79 : end if
80 :
81 : call miess(carma, &
82 : radius, &
83 : rfr, &
84 : rfi, &
85 : theta, &
86 : 1, &
87 : lqext, &
88 : lqsca, &
89 : qback,&
90 : ctbrqs, &
91 : rcore, &
92 : rfcr, &
93 : rfci, &
94 : wvno, &
95 0 : rc)
96 :
97 0 : lasym = ctbrqs / lqsca
98 :
99 0 : else if (miertn == I_MIERTN_BOHREN1983) then
100 :
101 0 : x = radius * wvno
102 :
103 : call bhmie(carma, &
104 : x, &
105 : m, &
106 : nang, &
107 : s1, &
108 : s2, &
109 : lqext, &
110 : lqsca, &
111 : qback, &
112 : lasym, &
113 0 : rc)
114 :
115 0 : else if (miertn == I_MIERTN_BOTET1997) then
116 :
117 0 : rfr = real(m)
118 0 : rfi = aimag(m)
119 :
120 0 : if (radius .le. rmonomer) then
121 0 : rmonomer_out = radius
122 0 : fractaldim_out = 3.0_f
123 : else
124 0 : rmonomer_out = rmonomer
125 0 : fractaldim_out = fractaldim
126 : end if
127 :
128 : call fractal_meanfield(carma, & !! carma object
129 : wavelength*1.0e4_f, & !! lambda in microns
130 : rfi, & !! imaginary index of refraction
131 : rfr, & !! real index of refraction
132 : nmonomer, & !! number of monomers
133 : falpha_in, & !! packing coefficient
134 : fractaldim_out, & !! fractal dimension
135 : rmonomer_out, & !! monomer size
136 : 1.0_f, & !! xv,"set to 1"
137 : 0.0_f, & !! angle, set to 0
138 : lqext, & !! extinction efficiency
139 : lqsca, & !! scattering efficiency
140 : lasym, & !! asymmetry parameter
141 0 : rc)
142 :
143 : else
144 0 : if (do_print) write(LUNOPRT, *) "mie::Unknown Mie routine specified."
145 0 : rc = RC_ERROR
146 : end if
147 :
148 : ! The mie code isn't perfect, so don't let it return values that aren't
149 : ! physical.
150 0 : lqext = max(lqext, 0._f)
151 0 : lqsca = max(0._f, min(lqext, lqsca))
152 0 : lasym = max(-1.0_f, min(1.0_f, lasym))
153 :
154 0 : return
155 0 : end subroutine mie
|
