source: trunk/LMDZ.GENERIC/utilities/aerosols_utilities/bhmie.f90 @ 3436

Last change on this file since 3436 was 3240, checked in by mturbet, 9 months ago

add optpropgen aerosol optical property code in the svn
File size: 9.5 KB
Line 
1module bhmie
2    use math, only: pi
3
4    implicit none
5
6    integer, parameter :: &
7        mxnang = 1000, &  ! Note: important that MXNANG be consistent with dimension of S1 and S2 in calling routine!
8        s_size = 2 * MXNANG - 1
9
10    save
11
12    contains
13        SUBROUTINE BH_MIE(X, REFREL, NANG, S1, S2, QEXT, QSCA, QBACK, GSCA)
14            implicit none
15            ! Declare parameters:
16
17            INTEGER NMXX
18            !      PARAMETER(MXNANG=1000,NMXX=15000)
19            PARAMETER(NMXX = 150000)
20            ! Arguments:
21            INTEGER NANG
22            doubleprecision GSCA, QBACK, QEXT, QSCA, X
23            COMPLEX(kind=8) REFREL
24            COMPLEX(kind=8) S1(s_size), S2(s_size)
25            ! Local variables:
26            INTEGER J, JJ, N, NSTOP, NMX, NN
27            DOUBLE PRECISION CHI, CHI0, CHI1, DANG, DX, EN, FN, P, PSI, PSI0, PSI1, &
28                    THETA, XSTOP, YMOD
29            DOUBLE PRECISION AMU(MXNANG), pi_(MXNANG), PI0(MXNANG), PI1(MXNANG), &
30                    TAU(MXNANG)
31            COMPLEX(kind=8) AN, AN1, BN, BN1, DREFRL, XI, XI1, Y
32            COMPLEX(kind=8) D(NMXX)
33            !***********************************************************************
34            ! Subroutine BHMIE is the Bohren-Huffman Mie scattering subroutine
35            !    to calculate scattering and absorption by a homogenous isotropic
36            !    sphere.
37            ! Given:
38            !    X = 2*pi*a/lambda
39            !    REFREL = (complex refr. index of sphere)/(real index of medium)
40            !    NANG = number of angles between 0 and 90 degrees
41            !           (will calculate 2*NANG-1 directions from 0 to 180 deg.)
42            !           if called with NANG<2, will set NANG=2 and will compute
43            !           scattering for theta=0,90,180.
44            ! Returns:
45            !    S1(1 - 2*NANG-1) = -i*f_22 (incid. E perp. to scatt. plane,
46            !                                scatt. E perp. to scatt. plane)
47            !    S2(1 - 2*NANG-1) = -i*f_11 (incid. E parr. to scatt. plane,
48            !                                scatt. E parr. to scatt. plane)
49            !    QEXT = C_ext/pi*a**2 = efficiency factor for extinction
50            !    QSCA = C_sca/pi*a**2 = efficiency factor for scattering
51            !    QBACK = (dC_sca/domega)/pi*a**2
52            !          = backscattering efficiency [NB: this is (1/4*pi) smaller
53            !            than the "radar backscattering efficiency"; see Bohren &
54            !            Huffman 1983 pp. 120-123]
55            !    GSCA = <cos(theta)> for scattering
56            !
57            ! Original program taken from Bohren and Huffman (1983), Appendix A
58            ! Modified by B.T.Draine, Princeton Univ. Obs., 90/10/26
59            ! in order to compute <cos(theta)>
60            ! 91/05/07 (BTD): Modified to allow NANG=1
61            ! 91/08/15 (BTD): Corrected error (failure to initialize P)
62            ! 91/08/15 (BTD): Modified to enhance vectorizability.
63            ! 91/08/15 (BTD): Modified to make NANG=2 if called with NANG=1
64            ! 91/08/15 (BTD): Changed definition of QBACK.
65            ! 92/01/08 (BTD): Converted to full double precision and double complex
66            !                 eliminated 2 unneed lines of code
67            !                 eliminated redundant variables (e.g. APSI,APSI0)
68            !                 renamed RN -> EN = double precision N
69            !                 Note that DOUBLE COMPLEX and DCMPLX are not part
70            !                 of f77 standard, so this version may not be fully
71            !                 portable.  In event that portable version is
72            !                 needed, use src/bhmie_f77.f
73            ! 93/06/01 (BTD): Changed AMAX1 to generic function MAX
74            !***********************************************************************
75            !*** Safety checks
76            IF(NANG>MXNANG)STOP '***Error: NANG > MXNANG in bhmie'
77            IF(NANG<2)NANG = 2
78
79            DX = X
80            DREFRL = REFREL
81            Y = X * DREFRL
82            YMOD = ABS(Y)
83            !
84            !*** Series expansion terminated after NSTOP terms
85            !    Logarithmic derivatives calculated from NMX on down
86            XSTOP = X + 4d0 * X ** (1d0 / 3d0) + 2d0
87            NMX = int(MAX(XSTOP, YMOD)) + 15
88            ! BTD experiment 91/1/15: add one more term to series and compare results
89            !      NMX=AMAX1(XSTOP,YMOD)+16
90            ! test: compute 7001 wavelengths between .0001 and 1000 micron
91            ! for a=1.0micron SiC grain.  When NMX increased by 1, only a single
92            ! computed number changed (out of 4*7001) and it only changed by 1/8387
93            ! conclusion: we are indeed retaining enough terms in series!
94            NSTOP = int(XSTOP)
95            !
96            IF(NMX>NMXX)THEN
97                WRITE(0, *)'Error: NMX > NMXX=', NMXX, ' for |m|x=', YMOD
98                STOP
99            ENDIF
100            !*** Require NANG.GE.1 in order to calculate scattering intensities
101            DANG = 0d0
102            IF(NANG>1)DANG = 0.5d0 * pi / DBLE(NANG - 1)
103            DO J = 1, NANG
104                THETA = DBLE(J - 1) * DANG
105                AMU(J) = COS(THETA)
106            end do
107            DO J = 1, NANG
108                PI0(J) = 0d0
109                PI1(J) = 1d0
110            end do
111            NN = 2 * NANG - 1
112            DO J = 1, NN
113                S1(J) = (0d0, 0d0)
114                S2(J) = (0d0, 0d0)
115            end do
116            !
117            !*** Logarithmic derivative D(J) calculated by downward recurrence
118            !    beginning with initial value (0.,0.) at J=NMX
119            !
120            D(NMX) = (0d0, 0d0)
121            NN = NMX - 1
122            DO N = 1, NN
123                EN = NMX - N + 1
124                D(NMX - N) = (EN / Y) - (1. / (D(NMX - N + 1) + EN / Y))
125            end do
126            !
127            !*** Riccati-Bessel functions with real argument X
128            !    calculated by upward recurrence
129            !
130            PSI0 = COS(DX)
131            PSI1 = SIN(DX)
132            CHI0 = -SIN(DX)
133            CHI1 = COS(DX)
134            XI1 = DCMPLX(PSI1, -CHI1)
135            QSCA = 0d0
136            GSCA = 0d0
137            P = -1d0
138            DO N = 1, NSTOP
139                EN = N
140                FN = (2d0 * EN + 1d0) / (EN * (EN + 1d0))
141                ! for given N, PSI  = psi_n        CHI  = chi_n
142                !              PSI1 = psi_{n-1}    CHI1 = chi_{n-1}
143                !              PSI0 = psi_{n-2}    CHI0 = chi_{n-2}
144                ! Calculate psi_n and chi_n
145                PSI = (2d0 * EN - 1d0) * PSI1 / DX - PSI0
146                CHI = (2d0 * EN - 1d0) * CHI1 / DX - CHI0
147                XI = DCMPLX(PSI, -CHI)
148                !
149                !*** Store previous values of AN and BN for use
150                !    in computation of g=<cos(theta)>
151                IF(N>1)THEN
152                    AN1 = AN
153                    BN1 = BN
154                ENDIF
155                !
156                !*** Compute AN and BN:
157                AN = (D(N) / DREFRL + EN / DX) * PSI - PSI1
158                AN = AN / ((D(N) / DREFRL + EN / DX) * XI - XI1)
159                BN = (DREFRL * D(N) + EN / DX) * PSI - PSI1
160                BN = BN / ((DREFRL * D(N) + EN / DX) * XI - XI1)
161                !
162                !*** Augment sums for Qsca and g=<cos(theta)>
163                QSCA = QSCA + (2d0 * EN + 1d0) * (ABS(AN)**2 + ABS(BN)**2)
164                GSCA = GSCA + ((2d0 * EN + 1d0) / (EN * (EN + 1d0))) * &
165                        (dble(AN) * dble(BN) + AIMAG(AN) * AIMAG(BN))
166                IF(N>1)THEN
167                    GSCA = GSCA + ((EN - 1d0) * (EN + 1d0) / EN) * &
168                            (dble(AN1) * dble(AN) + AIMAG(AN1) * AIMAG(AN) + &
169                                    dble(BN1) * dble(BN) + AIMAG(BN1) * AIMAG(BN))
170                ENDIF
171                !
172                !*** Now calculate scattering intensity pattern
173                !    First do angles from 0 to 90
174                DO J = 1, NANG
175                    JJ = 2 * NANG - J
176                    pi_(J) = PI1(J)
177                    TAU(J) = EN * AMU(J) * pi_(J) - (EN + 1d0) * PI0(J)
178                    S1(J) = S1(J) + FN * (AN * pi_(J) + BN * TAU(J))
179                    S2(J) = S2(J) + FN * (AN * TAU(J) + BN * pi_(J))
180                end do
181                !
182                !*** Now do angles greater than 90 using pi_ and TAU from
183                !    angles less than 90.
184                !    P=1 for N=1,3,...; P=-1 for N=2,4,...
185                P = -P
186                DO J = 1, NANG - 1
187                    JJ = 2 * NANG - J
188                    S1(JJ) = S1(JJ) + FN * P * (AN * pi_(J) - BN * TAU(J))
189                    S2(JJ) = S2(JJ) + FN * P * (BN * pi_(J) - AN * TAU(J))
190                end do
191                PSI0 = PSI1
192                PSI1 = PSI
193                CHI0 = CHI1
194                CHI1 = CHI
195                XI1 = DCMPLX(PSI1, -CHI1)
196                !
197                !*** Compute pi_n for next value of n
198                !    For each angle J, compute pi_n+1
199                !    from pi_ = pi_n , PI0 = pi_n-1
200                DO J = 1, NANG
201                    PI1(J) = ((2. * EN + 1.) * AMU(J) * pi_(J) - (EN + 1.) * PI0(J)) / EN
202                    PI0(J) = pi_(J)
203                end do
204            end do
205            !
206            !*** Have summed sufficient terms.
207            !    Now compute QSCA,QEXT,QBACK,and GSCA
208            GSCA = 2d0 * GSCA / QSCA
209            QSCA = (2d0 / (DX * DX)) * QSCA
210            QEXT = (4d0 / (DX * DX)) * dble(S1(1))
211            QBACK = (ABS(S1(2 * NANG - 1)) / DX)**2 / pi
212            RETURN
213        END
214end module bhmie
Note: See TracBrowser for help on using the repository browser.