module_fastj_mie.F
References to this file elsewhere.
1 !**********************************************************************************
2 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
3 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
4 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
5 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
6 !
7 ! Module to Compute Aerosol Optical Properties
8 ! * Primary investigator: James C. Barnard
9 ! * Co-investigators: Rahul A. Zaveri, Richard C. Easter, William I. Gustafson Jr.
10 ! Last update: September 2005
11 !
12 ! Contact:
13 ! Jerome D. Fast, PhD
14 ! Staff Scientist
15 ! Pacific Northwest National Laboratory
16 ! P.O. Box 999, MSIN K9-30
17 ! Richland, WA, 99352
18 ! Phone: (509) 372-6116
19 ! Email: Jerome.Fast@pnl.gov
20 !
21 ! Please report any bugs or problems to Jerome Fast, the WRF-chem implmentation
22 ! team leader for PNNL
23 !
24 ! Terms of Use:
25 ! 1) This module may not be included in commerical package or used for any
26 ! commercial applications without the consent of the PNNL contact.
27 ! 2) This module is provided to the WRF modeling community; however, no portion
28 ! of it can be used separately or in another code without the consent of the
29 ! PNNL contact.
30 ! 3) This module may be used for research, educational, and non-profit purposes
31 ! only. Any other usage must be first approved by the PNNL contact.
32 ! 4) Publications resulting from the usage of this module must use one the
33 ! reference below for proper acknowledgment.
34 !
35 ! Note that the aerosol optical properites are currently tied to the use of Fast-J
36 ! and MOSAIC. Future modifications will make the calculations more generic so the
37 ! code is not tied to the photolysis scheme and the code can be used for both
38 ! modal and sectional treatments.
39 !
40 ! References:
41 ! * Fast, J.D., W.I. Gustafson Jr., R.C. Easter, R.A. Zaveri, J.C. Barnard, E.G.
42 ! Chapman, G.A. Grell, and S.E. Peckham (2005), Evolution of ozone, particulates,
43 ! and aerosol direct radiative forcing in the vicinity of Houston using a fully-
44 ! coupled meteorology-chemistry-aerosol model, Submitted to J. Geophys. Res.
45 !
46 ! Contact Jerome Fast for updates on the status of manuscripts under review.
47 !
48 ! Additional information:
49 ! * www.pnl.gov/atmos_sciences/Jdf/wrfchem.html
50 !
51 ! Support:
52 ! Funding for adapting Fast-J was provided by the U.S. Department of Energy
53 ! under the auspices of Atmospheric Science Program of the Office of Biological
54 ! and Environmental Research the PNNL Laboratory Research and Directed Research
55 ! and Development program.
56 !**********************************************************************************
57 module module_fastj_mie
58
59
60 ! rce 2005-apr-22 - want lunerr = -1 for wrf-chem 3d;
61 ! also, define lunerr here, and it's available everywhere
62 integer, parameter :: lunerr = -1
63
64
65 contains
66
67 !***********************************************************************
68 ! <1.> subr mieaer
69 ! Purpose: calculate aerosol optical depth, single scattering albedo,
70 ! asymmetry factor, extinction, Legendre coefficients, and average aerosol
71 ! size. parameterizes aerosol coefficients using chebychev polynomials
72 ! requires double precision on 32-bit machines
73 ! uses Wiscombe's (1979) mie scattering code
74 ! INPUT
75 ! iclm, jclm -- i,j of grid column being processed
76 ! nbin_a -- number of bins
77 ! number_bin(nbin_a,kmaxd) -- number density in layer, #/cm^3
78 ! radius_wet(nbin_a,kmaxd) -- wet radius, cm
79 ! refindx(nbin_a,kmaxd) --volume averaged complex index of refraction
80 ! dz -- depth of individual cells in column, m
81 ! isecfrm0 - time from start of run, sec
82 ! lpar -- number of grid cells in vertical (via module_fastj_cmnh)
83 ! kmaxd -- predefined maximum allowed levels from module_data_mosaic_other
84 ! passed here via module_fastj_cmnh
85 ! OUTPUT: saved in module_fastj_cmnmie
86 ! real tauaer ! aerosol optical depth
87 ! waer ! aerosol single scattering albedo
88 ! gaer ! aerosol asymmetery factor
89 ! extaer ! aerosol extinction
90 ! l2,l3,l4,l5,l6,l7 ! Legendre coefficients, numbered 0,1,2,......
91 ! sizeaer ! average wet radius
92 !----------------------------------------------------------------------
93 subroutine mieaer( &
94 iclm, jclm, nbin_a, &
95 number_bin, radius_wet, refindx, &
96 dz, isecfrm0, lpar, &
97 sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
98
99 USE module_data_mosaic_other, only : kmaxd
100 USE module_data_mosaic_therm, ONLY: nbin_a_maxd
101 USE module_peg_util, only: peg_message
102
103 IMPLICIT NONE
104 ! subr arguments
105 !jdf
106 integer,parameter :: nspint = 4 ! Num of spectral intervals across
107 ! solar spectrum for FAST-J
108 integer, intent(in) :: lpar
109 real, dimension (nspint, kmaxd+1),intent(out) :: sizeaer,extaer,waer,gaer,tauaer
110 real, dimension (nspint, kmaxd+1),intent(out) :: l2,l3,l4,l5,l6,l7
111 real, dimension (nspint),save :: wavmid !cm
112 data wavmid &
113 / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
114 !jdf
115 integer, intent(in) :: iclm, jclm, nbin_a, isecfrm0
116 real, intent(in), dimension(nbin_a_maxd, kmaxd) :: number_bin
117 real, intent(inout), dimension(nbin_a_maxd, kmaxd) :: radius_wet
118 complex, intent(in) :: refindx(nbin_a_maxd, kmaxd)
119 real, intent(in) :: dz(lpar)
120
121 !local variables
122 real weighte, weights
123 ! various bookeeping variables
124 integer ltype ! total number of indicies of refraction
125 parameter (ltype = 1) ! bracket refractive indices based on information from Rahul, 2002/11/07
126 real x
127 real thesum ! for normalizing things
128 real sizem ! size in microns
129 integer kcallmieaer
130 !
131 integer m, j, nc, klevel
132 real pext ! parameterized specific extinction (cm2/g)
133 real pasm ! parameterized asymmetry factor
134 real ppmom2 ! 2 Lengendre expansion coefficient (numbered 0,1,2,...)
135 real ppmom3 ! 3 Lengendre expansion coefficient (numbered 0,1,2,...)
136 real ppmom4 ! 4 Lengendre expansion coefficient (numbered 0,1,2,...)
137 real ppmom5 ! 5 Lengendre expansion coefficient (numbered 0,1,2,...)
138 real ppmom6 ! 6 Lengendre expansion coefficient (numbered 0,1,2,...)
139 real ppmom7 ! 7 Lengendre expansion coefficient (numbered 0,1,2,...)
140
141 integer ns ! Spectral loop index
142 integer i ! Longitude loop index
143 integer k ! Level loop index
144 real pscat !scattering cross section
145
146 integer prefr,prefi,nrefr,nrefi,nr,ni
147 save nrefr,nrefi
148 parameter (prefr=7,prefi=7)
149
150 complex*16 sforw,sback,tforw(2),tback(2)
151 integer numang,nmom,ipolzn,momdim
152 real*8 pmom(0:7,1)
153 logical perfct,anyang,prnt(2)
154 logical first
155 integer, parameter :: nsiz=200,nlog=30,ncoef=50
156 !
157 real p2(nsiz),p3(nsiz),p4(nsiz),p5(nsiz)
158 real p6(nsiz),p7(nsiz)
159 !
160 real*8 xmu(1)
161 data xmu/1./,anyang/.false./
162 data numang/0/
163 complex*16 s1(1),s2(1)
164 data first/.true./
165 save first
166 real*8 mimcut
167 data perfct/.false./,mimcut/0.0/
168 data nmom/7/,ipolzn/0/,momdim/7/
169 data prnt/.false.,.false./
170 ! coefficients for parameterizing aerosol radiative properties
171 ! in terms of refractive index and wet radius
172 real extp(ncoef,prefr,prefi,nspint) ! specific extinction
173 real albp(ncoef,prefr,prefi,nspint) ! single scat alb
174 real asmp(ncoef,prefr,prefi,nspint) ! asymmetry factor
175 real ascat(ncoef,prefr,prefi,nspint) ! scattering efficiency, JCB 2004/02/09
176 real pmom2(ncoef,prefr,prefi,nspint) ! phase function expansion, #2
177 real pmom3(ncoef,prefr,prefi,nspint) ! phase function expansion, #3
178 real pmom4(ncoef,prefr,prefi,nspint) ! phase function expansion, #4
179 real pmom5(ncoef,prefr,prefi,nspint) ! phase function expansion, #5
180 real pmom6(ncoef,prefr,prefi,nspint) ! phase function expansion, #6
181 real pmom7(ncoef,prefr,prefi,nspint) ! phase function expansion, #7
182
183 save :: extp,albp,asmp,ascat,pmom2,pmom3,pmom4,pmom5,pmom6,pmom7
184 !--------------
185 real cext(ncoef),casm(ncoef),cpmom2(ncoef)
186 real cscat(ncoef) ! JCB 2004/02/09
187 real cpmom3(ncoef),cpmom4(ncoef),cpmom5(ncoef)
188 real cpmom6(ncoef),cpmom7(ncoef)
189 integer itab,jtab
190 real ttab,utab
191
192 ! nsiz = number of wet particle sizes
193 ! crefin = complex refractive index
194 integer n
195 real*8 thesize ! 2 pi radpart / waveleng = size parameter
196 real*8 qext(nsiz) ! array of extinction efficiencies
197 real*8 qsca(nsiz) ! array of scattering efficiencies
198 real*8 gqsc(nsiz) ! array of asymmetry factor * scattering efficiency
199 real asymm(nsiz) ! array of asymmetry factor
200 real scat(nsiz) ! JCB 2004/02/09
201 ! specabs = absorption coeff / unit dry mass
202 ! specscat = scattering coeff / unit dry mass
203 complex*16 crefin,crefd,crefw
204 save crefw
205 real, save :: rmin,rmax ! min, max aerosol size bin
206 data rmin /0.005e-4/ ! rmin in cm. 5e-3 microns min allowable size
207 data rmax /50.e-4 / ! rmax in cm. 50 microns, big particle, max allowable size
208
209 real bma,bpa
210
211 real, save :: xrmin,xrmax,xr
212 real rs(nsiz) ! surface mode radius (cm)
213 real xrad ! normalized aerosol radius
214 real ch(ncoef) ! chebychev polynomial
215
216 real, save :: rhoh2o ! density of liquid water (g/cm3)
217 data rhoh2o/1./
218
219 real refr ! real part of refractive index
220 real refi ! imaginary part of refractive index
221 real qextr4(nsiz) ! extinction, real*4
222
223 real refrmin ! minimum of real part of refractive index
224 real refrmax ! maximum of real part of refractive index
225 real refimin ! minimum of imag part of refractive index
226 real refimax ! maximum of imag part of refractive index
227 real drefr ! increment in real part of refractive index
228 real drefi ! increment in imag part of refractive index
229 real, save :: refrtab(prefr) ! table of real refractive indices for aerosols
230 real, save :: refitab(prefi) ! table of imag refractive indices for aerosols
231 complex specrefndx(ltype) ! refractivr indices
232 real pie,third
233
234 integer irams, jrams
235 ! diagnostic declarations
236 integer kcallmieaer2
237 integer ibin
238 character*80 msg
239
240 #if (defined(CHEM_DBG_I) && defined(CHEM_DBG_J) && defined(CHEM_DBG_K))
241 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
242 !ec diagnostics
243 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
244 !ec run_out.25 has aerosol physical parameter info for bins 1-8
245 !ec and vertical cells 1 to kmaxd.
246 ! ilaporte = 33
247 ! jlaporte = 34
248 if (iclm .eq. CHEM_DBG_I) then
249 if (jclm .eq. CHEM_DBG_J) then
250 ! initial entry
251 if (kcallmieaer2 .eq. 0) then
252 write(*,9099)iclm, jclm
253 9099 format('for cell i = ', i3, 2x, 'j = ', i3)
254 write(*,9100)
255 9100 format( &
256 'isecfrm0', 3x, 'i', 3x, 'j', 3x,'k', 3x, &
257 'ibin', 3x, &
258 'refindx(ibin,k)', 3x, &
259 'radius_wet(ibin,k)', 3x, &
260 'number_bin(ibin,k)' &
261 )
262 end if
263 !ec output for run_out.25
264 do k = 1, lpar
265 do ibin = 1, nbin_a
266 write(*, 9120) &
267 isecfrm0,iclm, jclm, k, ibin, &
268 refindx(ibin,k), &
269 radius_wet(ibin,k), &
270 number_bin(ibin,k)
271 9120 format( i7,3(2x,i4),2x,i4, 4x, 4(e14.6,2x))
272 end do
273 end do
274 kcallmieaer2 = kcallmieaer2 + 1
275 end if
276 end if
277 !ec end print of aerosol physical parameter diagnostics
278 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
279 #endif
280 !
281 ! assign fast-J wavelength, these values are in cm
282 ! wavmid(1)=0.30e-4
283 ! wavmid(2)=0.40e-4
284 ! wavmid(3)=0.60e-4
285 ! wavmid(4)=0.999e-04
286 !
287 pie=4.*atan(1.)
288 third=1./3.
289 if(first)then
290 first=.false.
291
292 ! parameterize aerosol radiative properties in terms of
293 ! relative humidity, surface mode wet radius, aerosol species,
294 ! and wavelength
295
296 ! first find min,max of real and imaginary parts of refractive index
297
298 ! real and imaginary parts of water refractive index
299
300 crefw=cmplx(1.33,0.0)
301 refrmin=real(crefw)
302 refrmax=real(crefw)
303 ! change Rahul's imaginary part of the refractive index from positive to negative
304 refimin=-imag(crefw)
305 refimax=-imag(crefw)
306
307 ! aerosol mode loop
308 specrefndx(1)=cmplx(1.82,-0.74) ! max values from Rahul, 7 Nov 2002
309 !
310
311 do i=1,ltype ! loop over all possible refractive indices
312
313 ! real and imaginary parts of aerosol refractive index
314
315 refrmin=amin1(refrmin,real(specrefndx(ltype)))
316 refrmax=amax1(refrmax,real(specrefndx(ltype)))
317 refimin=amin1(refimin,aimag(specrefndx(ltype)))
318 refimax=amax1(refimax,aimag(specrefndx(ltype)))
319
320 enddo
321
322 drefr=(refrmax-refrmin)
323 if(drefr.gt.1.e-4)then
324 nrefr=prefr
325 drefr=drefr/(nrefr-1)
326 else
327 nrefr=1
328 endif
329
330 drefi=(refimax-refimin)
331 if(drefi.gt.1.e-4)then
332 nrefi=prefi
333 drefi=drefi/(nrefi-1)
334 else
335 nrefi=1
336 endif
337
338
339 !
340 bma=0.5*alog(rmax/rmin) ! JCB
341 bpa=0.5*alog(rmax*rmin) ! JCB
342 xrmin=alog(rmin)
343 xrmax=alog(rmax)
344
345 ! wavelength loop
346
347 do 200 ns=1,nspint
348
349 ! calibrate parameterization with range of refractive indices
350
351 do 120 nr=1,nrefr
352 do 120 ni=1,nrefi
353
354 refrtab(nr)=refrmin+(nr-1)*drefr
355 refitab(ni)=refimin+(ni-1)*drefi
356 crefd=cmplx(refrtab(nr),refitab(ni))
357
358 ! mie calculations of optical efficiencies
359
360 do n=1,nsiz
361 xr=cos(pie*(float(n)-0.5)/float(nsiz))
362 rs(n)=exp(xr*bma+bpa)
363
364
365 ! size parameter and weighted refractive index
366
367 thesize=2.*pie*rs(n)/wavmid(ns)
368 thesize=min(thesize,10000.d0)
369
370 call miev0(thesize,crefd,perfct,mimcut,anyang, &
371 numang,xmu,nmom,ipolzn,momdim,prnt, &
372 qext(n),qsca(n),gqsc(n),pmom,sforw,sback,s1, &
373 s2,tforw,tback )
374 qextr4(n)=qext(n)
375 ! qabs(n)=qext(n)-qsca(n) ! not necessary anymore JCB 2004/02/09
376 scat(n)=qsca(n) ! JCB 2004/02/09
377 asymm(n)=gqsc(n)/qsca(n) ! assume always greater than zero
378 ! coefficients of phase function expansion; note modification by JCB of miev0 coefficients
379 p2(n)=pmom(2,1)/pmom(0,1)*5.0
380 p3(n)=pmom(3,1)/pmom(0,1)*7.0
381 p4(n)=pmom(4,1)/pmom(0,1)*9.0
382 p5(n)=pmom(5,1)/pmom(0,1)*11.0
383 p6(n)=pmom(6,1)/pmom(0,1)*13.0
384 p7(n)=pmom(7,1)/pmom(0,1)*15.0
385 enddo
386 100 continue
387 !
388 call fitcurv(rs,qextr4,extp(1,nr,ni,ns),ncoef,nsiz)
389 call fitcurv(rs,scat,ascat(1,nr,ni,ns),ncoef,nsiz) ! JCB 2004/02/07 - scattering efficiency
390 call fitcurv(rs,asymm,asmp(1,nr,ni,ns),ncoef,nsiz)
391 call fitcurv_nolog(rs,p2,pmom2(1,nr,ni,ns),ncoef,nsiz)
392 call fitcurv_nolog(rs,p3,pmom3(1,nr,ni,ns),ncoef,nsiz)
393 call fitcurv_nolog(rs,p4,pmom4(1,nr,ni,ns),ncoef,nsiz)
394 call fitcurv_nolog(rs,p5,pmom5(1,nr,ni,ns),ncoef,nsiz)
395 call fitcurv_nolog(rs,p6,pmom6(1,nr,ni,ns),ncoef,nsiz)
396 call fitcurv_nolog(rs,p7,pmom7(1,nr,ni,ns),ncoef,nsiz)
397
398 120 continue
399 200 continue
400
401
402 endif
403 ! begin level loop
404 do 2000 klevel=1,lpar
405 ! sum densities for normalization
406 thesum=0.0
407 do m=1,nbin_a
408 thesum=thesum+number_bin(m,klevel)
409 enddo
410 ! Begin spectral loop
411 do 1000 ns=1,nspint
412
413 ! aerosol optical properties
414
415 tauaer(ns,klevel)=0.
416 waer(ns,klevel)=0.
417 gaer(ns,klevel)=0.
418 sizeaer(ns,klevel)=0.0
419 extaer(ns,klevel)=0.0
420 l2(ns,klevel)=0.0
421 l3(ns,klevel)=0.0
422 l4(ns,klevel)=0.0
423 l5(ns,klevel)=0.0
424 l6(ns,klevel)=0.0
425 l7(ns,klevel)=0.0
426 if(thesum.le.0.0)goto 1000 ! set everything = 0 if no aerosol
427
428 ! loop over the bins
429 do m=1,nbin_a ! nbin_a is number of bins
430 ! check to see if there's any aerosol
431 if(number_bin(m,klevel).le.0.0)goto 70 ! no aerosol
432 ! here's the size
433 sizem=radius_wet(m,klevel) ! radius in cm
434 ! check limits of particle size
435 ! rce 2004-dec-07 - use klevel in write statements
436 if(radius_wet(m,klevel).le.rmin)then
437 radius_wet(m,klevel)=rmin
438 write( msg, '(a, 4i4,1x, e11.4)' ) &
439 'FASTJ mie: radius_wet set to rmin,' // &
440 'i,j,k,m,rm(m,k)', iclm, jclm, klevel, m, radius_wet(m,klevel)
441 call peg_message( lunerr, msg )
442 ! write(6,'('' particle size too small '')')
443 endif
444 !
445 if(radius_wet(m,klevel).gt.rmax)then
446 radius_wet(m,klevel)=rmax
447 write( msg, '(a, 4i4,1x, e11.4)' ) &
448 'FASTJ mie: radius_wet set to rmax,' // &
449 'i,j,k,m,rm(m,k)', iclm, jclm, klevel, m, radius_wet(m,klevel)
450 call peg_message( lunerr, msg )
451 ! write(6,'('' particle size too large '')')
452 endif
453 !
454 x=alog(radius_wet(m,klevel)) ! radius in cm
455 !
456 crefin=refindx(m,klevel)
457 refr=real(crefin)
458 ! change Rahul's imaginary part of the index of refraction from positive to negative
459 refi=-imag(crefin)
460 !
461 xrad=x
462
463 thesize=2.0*pie*exp(x)/wavmid(ns)
464 ! normalize size parameter
465 xrad=(2*xrad-xrmax-xrmin)/(xrmax-xrmin)
466 ! retain this diagnostic code
467 if(abs(refr).gt.10.0.or.abs(refr).le.0.001)then
468 write ( msg, '(a,1x, e14.5)' ) &
469 'FASTJ mie /refr/ outside range 1e-3 - 10 ' // &
470 'refr= ', refr
471 call peg_message( lunerr, msg )
472 ! print *,'refr=',refr
473 call exit(1)
474 endif
475 if(abs(refi).gt.10.)then
476 write ( msg, '(a,1x, e14.5)' ) &
477 'FASTJ mie /refi/ >10 ' // &
478 'refi', refi
479 call peg_message( lunerr, msg )
480 ! print *,'refi=',refi
481 call exit(1)
482 endif
483
484 ! interpolate coefficients linear in refractive index
485 ! first call calcs itab,jtab,ttab,utab
486 itab=0
487 call binterp(extp(1,1,1,ns),ncoef,nrefr,nrefi, &
488 refr,refi,refrtab,refitab,itab,jtab, &
489 ttab,utab,cext)
490
491 ! JCB 2004/02/09 -- new code for scattering cross section
492 call binterp(ascat(1,1,1,ns),ncoef,nrefr,nrefi, &
493 refr,refi,refrtab,refitab,itab,jtab, &
494 ttab,utab,cscat)
495 call binterp(asmp(1,1,1,ns),ncoef,nrefr,nrefi, &
496 refr,refi,refrtab,refitab,itab,jtab, &
497 ttab,utab,casm)
498 call binterp(pmom2(1,1,1,ns),ncoef,nrefr,nrefi, &
499 refr,refi,refrtab,refitab,itab,jtab, &
500 ttab,utab,cpmom2)
501 call binterp(pmom3(1,1,1,ns),ncoef,nrefr,nrefi, &
502 refr,refi,refrtab,refitab,itab,jtab, &
503 ttab,utab,cpmom3)
504 call binterp(pmom4(1,1,1,ns),ncoef,nrefr,nrefi, &
505 refr,refi,refrtab,refitab,itab,jtab, &
506 ttab,utab,cpmom4)
507 call binterp(pmom5(1,1,1,ns),ncoef,nrefr,nrefi, &
508 refr,refi,refrtab,refitab,itab,jtab, &
509 ttab,utab,cpmom5)
510 call binterp(pmom6(1,1,1,ns),ncoef,nrefr,nrefi, &
511 refr,refi,refrtab,refitab,itab,jtab, &
512 ttab,utab,cpmom6)
513 call binterp(pmom7(1,1,1,ns),ncoef,nrefr,nrefi, &
514 refr,refi,refrtab,refitab,itab,jtab, &
515 ttab,utab,cpmom7)
516
517 ! chebyshev polynomials
518
519 ch(1)=1.
520 ch(2)=xrad
521 do nc=3,ncoef
522 ch(nc)=2.*xrad*ch(nc-1)-ch(nc-2)
523 enddo
524 ! parameterized optical properties
525
526 pext=0.5*cext(1)
527 do nc=2,ncoef
528 pext=pext+ch(nc)*cext(nc)
529 enddo
530 pext=exp(pext)
531
532 ! JCB 2004/02/09 -- for scattering efficiency
533 pscat=0.5*cscat(1)
534 do nc=2,ncoef
535 pscat=pscat+ch(nc)*cscat(nc)
536 enddo
537 pscat=exp(pscat)
538 !
539 pasm=0.5*casm(1)
540 do nc=2,ncoef
541 pasm=pasm+ch(nc)*casm(nc)
542 enddo
543 pasm=exp(pasm)
544 !
545 ppmom2=0.5*cpmom2(1)
546 do nc=2,ncoef
547 ppmom2=ppmom2+ch(nc)*cpmom2(nc)
548 enddo
549 if(ppmom2.le.0.0)ppmom2=0.0
550 !
551 ppmom3=0.5*cpmom3(1)
552 do nc=2,ncoef
553 ppmom3=ppmom3+ch(nc)*cpmom3(nc)
554 enddo
555 ! ppmom3=exp(ppmom3) ! no exponentiation required
556 if(ppmom3.le.0.0)ppmom3=0.0
557 !
558 ppmom4=0.5*cpmom4(1)
559 do nc=2,ncoef
560 ppmom4=ppmom4+ch(nc)*cpmom4(nc)
561 enddo
562 if(ppmom4.le.0.0.or.sizem.le.0.03e-04)ppmom4=0.0
563 !
564 ppmom5=0.5*cpmom5(1)
565 do nc=2,ncoef
566 ppmom5=ppmom5+ch(nc)*cpmom5(nc)
567 enddo
568 if(ppmom5.le.0.0.or.sizem.le.0.03e-04)ppmom5=0.0
569 !
570 ppmom6=0.5*cpmom6(1)
571 do nc=2,ncoef
572 ppmom6=ppmom6+ch(nc)*cpmom6(nc)
573 enddo
574 if(ppmom6.le.0.0.or.sizem.le.0.03e-04)ppmom6=0.0
575 !
576 ppmom7=0.5*cpmom7(1)
577 do nc=2,ncoef
578 ppmom7=ppmom7+ch(nc)*cpmom7(nc)
579 enddo
580 if(ppmom7.le.0.0.or.sizem.le.0.03e-04)ppmom7=0.0
581 !
582 ! weights:
583 weighte=pext*pie*exp(x)**2 ! JCB, extinction cross section
584 weights=pscat*pie*exp(x)**2 ! JCB, scattering cross section
585 tauaer(ns,klevel)=tauaer(ns,klevel)+weighte*number_bin(m,klevel) ! must be multiplied by deltaZ
586 extaer(ns,klevel)=extaer(ns,klevel)+pext*number_bin(m,klevel)
587 sizeaer(ns,klevel)=sizeaer(ns,klevel)+exp(x)*10000.0* &
588 number_bin(m,klevel)
589 waer(ns,klevel)=waer(ns,klevel)+weights*number_bin(m,klevel) !JCB
590 gaer(ns,klevel)=gaer(ns,klevel)+pasm*weights* &
591 number_bin(m,klevel) !JCB
592 ! need weighting by scattering cross section ? JCB 2004/02/09
593 l2(ns,klevel)=l2(ns,klevel)+weights*ppmom2*number_bin(m,klevel)
594 l3(ns,klevel)=l3(ns,klevel)+weights*ppmom3*number_bin(m,klevel)
595 l4(ns,klevel)=l4(ns,klevel)+weights*ppmom4*number_bin(m,klevel)
596 l5(ns,klevel)=l5(ns,klevel)+weights*ppmom5*number_bin(m,klevel)
597 l6(ns,klevel)=l6(ns,klevel)+weights*ppmom6*number_bin(m,klevel)
598 l7(ns,klevel)=l7(ns,klevel)+weights*ppmom7*number_bin(m,klevel)
599
600 end do ! end of nbin_a loop
601 ! take averages - weighted by cross section - new code JCB 2004/02/09
602 extaer(ns,klevel)=extaer(ns,klevel)/thesum
603 sizeaer(ns,klevel)=sizeaer(ns,klevel)/thesum
604 gaer(ns,klevel)=gaer(ns,klevel)/waer(ns,klevel) ! JCB removed *3 factor 2/9/2004
605 ! because factor is applied in subroutine opmie, file zz01fastj_mod.f
606 l2(ns,klevel)=l2(ns,klevel)/waer(ns,klevel)
607 l3(ns,klevel)=l3(ns,klevel)/waer(ns,klevel)
608 l4(ns,klevel)=l4(ns,klevel)/waer(ns,klevel)
609 l5(ns,klevel)=l5(ns,klevel)/waer(ns,klevel)
610 l6(ns,klevel)=l6(ns,klevel)/waer(ns,klevel)
611 l7(ns,klevel)=l7(ns,klevel)/waer(ns,klevel)
612 ! this must be last!! JCB 2007/02/09
613 waer(ns,klevel)=waer(ns,klevel)/tauaer(ns,klevel) ! JCB
614
615 70 continue ! bail out if no aerosol;go on to next wavelength bin
616
617 1000 continue ! end of wavelength loop
618
619
620 2000 continue ! end of klevel loop
621 !
622 ! before returning, multiply tauaer by depth of individual cells.
623 ! tauaer is in cm, dz in m; multiply dz by 100 to convert from m to cm.
624 do ns = 1, nspint
625 do klevel = 1, lpar
626 tauaer(ns,klevel) = tauaer(ns,klevel) * dz(klevel)* 100.
627 end do
628 end do
629
630 #if (defined(CHEM_DBG_I) && defined(CHEM_DBG_J) && defined(CHEM_DBG_K))
631 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
632 !ec fastj diagnostics
633 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
634 !ec run_out.30 has aerosol optical info for cells 1 to kmaxd.
635 ! ilaporte = 33
636 ! jlaporte = 34
637 if (iclm .eq. CHEM_DBG_I) then
638 if (jclm .eq. CHEM_DBG_J) then
639 ! initial entry
640 if (kcallmieaer .eq. 0) then
641 write(*,909) CHEM_DBG_I, CHEM_DBG_J
642 909 format( ' for cell i = ', i3, ' j = ', i3)
643 write(*,910)
644 910 format( &
645 'isecfrm0', 3x, 'i', 3x, 'j', 3x,'k', 3x, &
646 'dzmfastj', 8x, &
647 'tauaer(1,k)',1x, 'tauaer(2,k)',1x,'tauaer(3,k)',3x, &
648 'tauaer(4,k)',5x, &
649 'waer(1,k)', 7x, 'waer(2,k)', 7x,'waer(3,k)', 7x, &
650 'waer(4,k)', 7x, &
651 'gaer(1,k)', 7x, 'gaer(2,k)', 7x,'gaer(3,k)', 7x, &
652 'gaer(4,k)', 7x, &
653 'extaer(1,k)',5x, 'extaer(2,k)',5x,'extaer(3,k)',5x, &
654 'extaer(4,k)',5x, &
655 'sizeaer(1,k)',4x, 'sizeaer(2,k)',4x,'sizeaer(3,k)',4x, &
656 'sizeaer(4,k)' )
657 end if
658 !ec output for run_out.30
659 do k = 1, lpar
660 write(*, 912) &
661 isecfrm0,iclm, jclm, k, &
662 dz(k) , &
663 (tauaer(n,k), n=1,4), &
664 (waer(n,k), n=1,4), &
665 (gaer(n,k), n=1,4), &
666 (extaer(n,k), n=1,4), &
667 (sizeaer(n,k), n=1,4)
668 912 format( i7,3(2x,i4),2x,21(e14.6,2x))
669 end do
670 kcallmieaer = kcallmieaer + 1
671 end if
672 end if
673 !ec end print of fastj diagnostics
674 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
675 #endif
676
677 return
678 end subroutine mieaer
679 !****************************************************************
680 subroutine fitcurv(rs,yin,coef,ncoef,maxm)
681
682 ! fit y(x) using Chebychev polynomials
683 ! wig 7-Sep-2004: Removed dependency on pre-determined maximum
684 ! array size and replaced with f90 array info.
685
686 USE module_peg_util, only: peg_message
687
688 IMPLICIT NONE
689 ! integer nmodes, nrows, maxm, ncoef
690 ! parameter (nmodes=500,nrows=8)
691 integer, intent(in) :: maxm, ncoef
692
693 ! real rs(nmodes),yin(nmodes),coef(ncoef)
694 ! real x(nmodes),y(nmodes)
695 real, dimension(ncoef) :: coef
696 real, dimension(:) :: rs, yin
697 real x(size(rs)),y(size(yin))
698
699 integer m
700 real xmin, xmax
701 character*80 msg
702
703 !!$ if(maxm.gt.nmodes)then
704 !!$ write ( msg, '(a, 1x,i6)' ) &
705 !!$ 'FASTJ mie nmodes too small in fitcurv, ' // &
706 !!$ 'maxm ', maxm
707 !!$ call peg_message( lunerr, msg )
708 !!$! write(*,*)'nmodes too small in fitcurv',maxm
709 !!$ call exit(1)
710 !!$ endif
711
712 do 100 m=1,maxm
713 x(m)=alog(rs(m))
714 y(m)=alog(yin(m))
715 100 continue
716
717 xmin=x(1)
718 xmax=x(maxm)
719 do 110 m=1,maxm
720 x(m)=(2*x(m)-xmax-xmin)/(xmax-xmin)
721 110 continue
722
723 call chebft(coef,ncoef,maxm,y)
724
725 return
726 end subroutine fitcurv
727 !**************************************************************
728 subroutine fitcurv_nolog(rs,yin,coef,ncoef,maxm)
729
730 ! fit y(x) using Chebychev polynomials
731 ! wig 7-Sep-2004: Removed dependency on pre-determined maximum
732 ! array size and replaced with f90 array info.
733
734 USE module_peg_util, only: peg_message
735 IMPLICIT NONE
736
737 ! integer nmodes, nrows, maxm, ncoef
738 ! parameter (nmodes=500,nrows=8)
739 integer, intent(in) :: maxm, ncoef
740
741 ! real rs(nmodes),yin(nmodes),coef(ncoef)
742 real, dimension(:) :: rs, yin
743 real, dimension(ncoef) :: coef(ncoef)
744 real x(size(rs)),y(size(yin))
745
746 integer m
747 real xmin, xmax
748 character*80 msg
749
750 !!$ if(maxm.gt.nmodes)then
751 !!$ write ( msg, '(a,1x, i6)' ) &
752 !!$ 'FASTJ mie nmodes too small in fitcurv ' // &
753 !!$ 'maxm ', maxm
754 !!$ call peg_message( lunerr, msg )
755 !!$! write(*,*)'nmodes too small in fitcurv',maxm
756 !!$ call exit(1)
757 !!$ endif
758
759 do 100 m=1,maxm
760 x(m)=alog(rs(m))
761 y(m)=yin(m) ! note, no "alog" here
762 100 continue
763
764 xmin=x(1)
765 xmax=x(maxm)
766 do 110 m=1,maxm
767 x(m)=(2*x(m)-xmax-xmin)/(xmax-xmin)
768 110 continue
769
770 call chebft(coef,ncoef,maxm,y)
771
772 return
773 end subroutine fitcurv_nolog
774 !************************************************************************
775 subroutine chebft(c,ncoef,n,f)
776 ! given a function f with values at zeroes x_k of Chebychef polynomial
777 ! T_n(x), calculate coefficients c_j such that
778 ! f(x)=sum(k=1,n) c_k t_(k-1)(y) - 0.5*c_1
779 ! where y=(x-0.5*(xmax+xmin))/(0.5*(xmax-xmin))
780 ! See Numerical Recipes, pp. 148-150.
781
782 IMPLICIT NONE
783 real pi
784 integer ncoef, n
785 parameter (pi=3.14159265)
786 real c(ncoef),f(n)
787
788 ! local variables
789 real fac, thesum
790 integer j, k
791
792 fac=2./n
793 do j=1,ncoef
794 thesum=0
795 do k=1,n
796 thesum=thesum+f(k)*cos((pi*(j-1))*((k-0.5)/n))
797 enddo
798 c(j)=fac*thesum
799 enddo
800 return
801 end subroutine chebft
802 !*************************************************************************
803 subroutine binterp(table,km,im,jm,x,y,xtab,ytab,ix,jy,t,u,out)
804
805 ! bilinear interpolation of table
806 !
807 implicit none
808 integer im,jm,km
809 real table(km,im,jm),xtab(im),ytab(jm),out(km)
810 integer i,ix,ip1,j,jy,jp1,k
811 real x,dx,t,y,dy,u,tu, tuc,tcu,tcuc
812
813 if(ix.gt.0)go to 30
814 if(im.gt.1)then
815 do i=1,im
816 if(x.lt.xtab(i))go to 10
817 enddo
818 10 ix=max0(i-1,1)
819 ip1=min0(ix+1,im)
820 dx=(xtab(ip1)-xtab(ix))
821 if(abs(dx).gt.1.e-20)then
822 t=(x-xtab(ix))/(xtab(ix+1)-xtab(ix))
823 else
824 t=0
825 endif
826 else
827 ix=1
828 ip1=1
829 t=0
830 endif
831 if(jm.gt.1)then
832 do j=1,jm
833 if(y.lt.ytab(j))go to 20
834 enddo
835 20 jy=max0(j-1,1)
836 jp1=min0(jy+1,jm)
837 dy=(ytab(jp1)-ytab(jy))
838 if(abs(dy).gt.1.e-20)then
839 u=(y-ytab(jy))/dy
840 else
841 u=0
842 endif
843 else
844 jy=1
845 jp1=1
846 u=0
847 endif
848 30 continue
849 jp1=min(jy+1,jm)
850 ip1=min(ix+1,im)
851 tu=t*u
852 tuc=t-tu
853 tcuc=1-tuc-u
854 tcu=u-tu
855 do k=1,km
856 out(k)=tcuc*table(k,ix,jy)+tuc*table(k,ip1,jy) &
857 +tu*table(k,ip1,jp1)+tcu*table(k,ix,jp1)
858 enddo
859 return
860 end subroutine binterp
861 !***************************************************************
862 subroutine miev0 ( xx, crefin, perfct, mimcut, anyang, &
863 numang, xmu, nmom, ipolzn, momdim, prnt, &
864 qext, qsca, gqsc, pmom, sforw, sback, s1, &
865 s2, tforw, tback )
866 !
867 ! computes mie scattering and extinction efficiencies; asymmetry
868 ! factor; forward- and backscatter amplitude; scattering
869 ! amplitudes for incident polarization parallel and perpendicular
870 ! to the plane of scattering, as functions of scattering angle;
871 ! coefficients in the legendre polynomial expansions of either the
872 ! unpolarized phase function or the polarized phase matrix;
873 ! and some quantities needed in polarized radiative transfer.
874 !
875 ! calls : biga, ckinmi, small1, small2, testmi, miprnt,
876 ! lpcoef, errmsg
877 !
878 ! i n t e r n a l v a r i a b l e s
879 ! -----------------------------------
880 !
881 ! an,bn mie coefficients little-a-sub-n, little-b-sub-n
882 ! ( ref. 1, eq. 16 )
883 ! anm1,bnm1 mie coefficients little-a-sub-(n-1),
884 ! little-b-sub-(n-1); used in -gqsc- sum
885 ! anp coeffs. in s+ expansion ( ref. 2, p. 1507 )
886 ! bnp coeffs. in s- expansion ( ref. 2, p. 1507 )
887 ! anpm coeffs. in s+ expansion ( ref. 2, p. 1507 )
888 ! when mu is replaced by - mu
889 ! bnpm coeffs. in s- expansion ( ref. 2, p. 1507 )
890 ! when mu is replaced by - mu
891 ! calcmo(k) true, calculate moments for k-th phase quantity
892 ! (derived from -ipolzn-; used only in 'lpcoef')
893 ! cbiga(n) bessel function ratio capital-a-sub-n (ref. 2, eq. 2)
894 ! ( complex version )
895 ! cior complex index of refraction with negative
896 ! imaginary part (van de hulst convention)
897 ! cioriv 1 / cior
898 ! coeff ( 2n + 1 ) / ( n ( n + 1 ) )
899 ! fn floating point version of index in loop performing
900 ! mie series summation
901 ! lita,litb(n) mie coefficients -an-, -bn-, saved in arrays for
902 ! use in calculating legendre moments *pmom*
903 ! maxtrm max. possible no. of terms in mie series
904 ! mm + 1 and - 1, alternately.
905 ! mim magnitude of imaginary refractive index
906 ! mre real part of refractive index
907 ! maxang max. possible value of input variable -numang-
908 ! nangd2 (numang+1)/2 ( no. of angles in 0-90 deg; anyang=f )
909 ! noabs true, sphere non-absorbing (determined by -mimcut-)
910 ! np1dn ( n + 1 ) / n
911 ! npquan highest-numbered phase quantity for which moments are
912 ! to be calculated (the largest digit in -ipolzn-
913 ! if ipolzn .ne. 0)
914 ! ntrm no. of terms in mie series
915 ! pass1 true on first entry, false thereafter; for self-test
916 ! pin(j) angular function little-pi-sub-n ( ref. 2, eq. 3 )
917 ! at j-th angle
918 ! pinm1(j) little-pi-sub-(n-1) ( see -pin- ) at j-th angle
919 ! psinm1 ricatti-bessel function psi-sub-(n-1), argument -xx-
920 ! psin ricatti-bessel function psi-sub-n of argument -xx-
921 ! ( ref. 1, p. 11 ff. )
922 ! rbiga(n) bessel function ratio capital-a-sub-n (ref. 2, eq. 2)
923 ! ( real version, for when imag refrac index = 0 )
924 ! rioriv 1 / mre
925 ! rn 1 / n
926 ! rtmp (real) temporary variable
927 ! sp(j) s+ for j-th angle ( ref. 2, p. 1507 )
928 ! sm(j) s- for j-th angle ( ref. 2, p. 1507 )
929 ! sps(j) s+ for (numang+1-j)-th angle ( anyang=false )
930 ! sms(j) s- for (numang+1-j)-th angle ( anyang=false )
931 ! taun angular function little-tau-sub-n ( ref. 2, eq. 4 )
932 ! at j-th angle
933 ! tcoef n ( n+1 ) ( 2n+1 ) (for summing tforw,tback series)
934 ! twonp1 2n + 1
935 ! yesang true if scattering amplitudes are to be calculated
936 ! zetnm1 ricatti-bessel function zeta-sub-(n-1) of argument
937 ! -xx- ( ref. 2, eq. 17 )
938 ! zetn ricatti-bessel function zeta-sub-n of argument -xx-
939 !
940 ! ----------------------------------------------------------------------
941 ! -------- i / o specifications for subroutine miev0 -----------------
942 ! ----------------------------------------------------------------------
943 implicit none
944 logical anyang, perfct, prnt(*)
945 integer ipolzn, momdim, numang, nmom
946 real*8 gqsc, mimcut, pmom( 0:momdim, * ), qext, qsca, &
947 xmu(*), xx
948 complex*16 crefin, sforw, sback, s1(*), s2(*), tforw(*), &
949 tback(*)
950 integer maxang,mxang2,maxtrm
951 real*8 onethr
952 ! ----------------------------------------------------------------------
953 !
954 parameter ( maxang = 501, mxang2 = maxang/2 + 1 )
955 !
956 ! ** note -- maxtrm = 10100 is neces-
957 ! ** sary to do some of the test probs,
958 ! ** but 1100 is sufficient for most
959 ! ** conceivable applications
960 parameter ( maxtrm = 1100 )
961 parameter ( onethr = 1./3. )
962 !
963 logical anysav, calcmo(4), noabs, ok, pass1, persav, yesang
964 integer npquan
965 integer i,j,n,nmosav,iposav,numsav,ntrm,nangd2
966 real*8 mim, mimsav, mre, mm, np1dn
967 real*8 rioriv,xmusav,xxsav,sq,fn,rn,twonp1,tcoef, coeff
968 real*8 xinv,psinm1,chinm1,psin,chin,rtmp,taun
969 real*8 rbiga( maxtrm ), pin( maxang ), pinm1( maxang )
970 complex*16 an, bn, anm1, bnm1, anp, bnp, anpm, bnpm, cresav, &
971 cior, cioriv, ctmp, zet, zetnm1, zetn
972 complex*16 cbiga( maxtrm ), lita( maxtrm ), litb( maxtrm ), &
973 sp( maxang ), sm( maxang ), sps( mxang2 ), sms( mxang2 )
974 equivalence ( cbiga, rbiga )
975 save pass1
976 sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
977 data pass1 / .true. /
978 !
979 !
980 if ( pass1 ) then
981 ! ** save certain user input values
982 xxsav = xx
983 cresav = crefin
984 mimsav = mimcut
985 persav = perfct
986 anysav = anyang
987 nmosav = nmom
988 iposav = ipolzn
989 numsav = numang
990 xmusav = xmu( 1 )
991 ! ** reset input values for test case
992 xx = 10.0
993 crefin = ( 1.5, - 0.1 )
994 perfct = .false.
995 mimcut = 0.0
996 anyang = .true.
997 numang = 1
998 xmu( 1 )= - 0.7660444
999 nmom = 1
1000 ipolzn = - 1
1001 !
1002 end if
1003 ! ** check input and calculate
1004 ! ** certain variables from input
1005 !
1006 10 call ckinmi( numang, maxang, xx, perfct, crefin, momdim, &
1007 nmom, ipolzn, anyang, xmu, calcmo, npquan )
1008 !
1009 if ( perfct .and. xx .le. 0.1 ) then
1010 ! ** use totally-reflecting
1011 ! ** small-particle limit
1012 !
1013 call small1 ( xx, numang, xmu, qext, qsca, gqsc, sforw, &
1014 sback, s1, s2, tforw, tback, lita, litb )
1015 ntrm = 2
1016 go to 200
1017 end if
1018 !
1019 if ( .not.perfct ) then
1020 !
1021 cior = crefin
1022 if ( dimag( cior ) .gt. 0.0 ) cior = dconjg( cior )
1023 mre = dble( cior )
1024 mim = - dimag( cior )
1025 noabs = mim .le. mimcut
1026 cioriv = 1.0 / cior
1027 rioriv = 1.0 / mre
1028 !
1029 if ( xx * dmax1( 1.d0, cdabs(cior) ) .le. 0.d1 ) then
1030 !
1031 ! ** use general-refractive-index
1032 ! ** small-particle limit
1033 ! ** ( ref. 2, p. 1508 )
1034 !
1035 call small2 ( xx, cior, .not.noabs, numang, xmu, qext, &
1036 qsca, gqsc, sforw, sback, s1, s2, tforw, &
1037 tback, lita, litb )
1038 ntrm = 2
1039 go to 200
1040 end if
1041 !
1042 end if
1043 !
1044 nangd2 = ( numang + 1 ) / 2
1045 yesang = numang .gt. 0
1046 ! ** estimate number of terms in mie series
1047 ! ** ( ref. 2, p. 1508 )
1048 if ( xx.le.8.0 ) then
1049 ntrm = xx + 4. * xx**onethr + 1.
1050 else if ( xx.lt.4200. ) then
1051 ntrm = xx + 4.05 * xx**onethr + 2.
1052 else
1053 ntrm = xx + 4. * xx**onethr + 2.
1054 end if
1055 if ( ntrm+1 .gt. maxtrm ) &
1056 call errmsg( 'miev0--parameter maxtrm too small', .true. )
1057 !
1058 ! ** calculate logarithmic derivatives of
1059 ! ** j-bessel-fcn., big-a-sub-(1 to ntrm)
1060 if ( .not.perfct ) &
1061 call biga( cior, xx, ntrm, noabs, yesang, rbiga, cbiga )
1062 !
1063 ! ** initialize ricatti-bessel functions
1064 ! ** (psi,chi,zeta)-sub-(0,1) for upward
1065 ! ** recurrence ( ref. 1, eq. 19 )
1066 xinv = 1.0 / xx
1067 psinm1 = dsin( xx )
1068 chinm1 = dcos( xx )
1069 psin = psinm1 * xinv - chinm1
1070 chin = chinm1 * xinv + psinm1
1071 zetnm1 = dcmplx( psinm1, chinm1 )
1072 zetn = dcmplx( psin, chin )
1073 ! ** initialize previous coeffi-
1074 ! ** cients for -gqsc- series
1075 anm1 = ( 0.0, 0.0 )
1076 bnm1 = ( 0.0, 0.0 )
1077 ! ** initialize angular function little-pi
1078 ! ** and sums for s+, s- ( ref. 2, p. 1507 )
1079 if ( anyang ) then
1080 do 60 j = 1, numang
1081 pinm1( j ) = 0.0
1082 pin( j ) = 1.0
1083 sp ( j ) = ( 0.0, 0.0 )
1084 sm ( j ) = ( 0.0, 0.0 )
1085 60 continue
1086 else
1087 do 70 j = 1, nangd2
1088 pinm1( j ) = 0.0
1089 pin( j ) = 1.0
1090 sp ( j ) = ( 0.0, 0.0 )
1091 sm ( j ) = ( 0.0, 0.0 )
1092 sps( j ) = ( 0.0, 0.0 )
1093 sms( j ) = ( 0.0, 0.0 )
1094 70 continue
1095 end if
1096 ! ** initialize mie sums for efficiencies, etc.
1097 qsca = 0.0
1098 gqsc = 0.0
1099 sforw = ( 0., 0. )
1100 sback = ( 0., 0. )
1101 tforw( 1 ) = ( 0., 0. )
1102 tback( 1 ) = ( 0., 0. )
1103 !
1104 !
1105 ! --------- loop to sum mie series -----------------------------------
1106 !
1107 mm = + 1.0
1108 do 100 n = 1, ntrm
1109 ! ** compute various numerical coefficients
1110 fn = n
1111 rn = 1.0 / fn
1112 np1dn = 1.0 + rn
1113 twonp1 = 2 * n + 1
1114 coeff = twonp1 / ( fn * ( n + 1 ) )
1115 tcoef = twonp1 * ( fn * ( n + 1 ) )
1116 !
1117 ! ** calculate mie series coefficients
1118 if ( perfct ) then
1119 ! ** totally-reflecting case
1120 !
1121 an = ( ( fn*xinv ) * psin - psinm1 ) / &
1122 ( ( fn*xinv ) * zetn - zetnm1 )
1123 bn = psin / zetn
1124 !
1125 else if ( noabs ) then
1126 ! ** no-absorption case
1127 !
1128 an = ( ( rioriv*rbiga(n) + ( fn*xinv ) ) * psin - psinm1 ) &
1129 / ( ( rioriv*rbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
1130 bn = ( ( mre * rbiga(n) + ( fn*xinv ) ) * psin - psinm1 ) &
1131 / ( ( mre * rbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
1132 else
1133 ! ** absorptive case
1134 !
1135 an = ( ( cioriv * cbiga(n) + ( fn*xinv ) ) * psin - psinm1 ) &
1136 /( ( cioriv * cbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
1137 bn = ( ( cior * cbiga(n) + ( fn*xinv ) ) * psin - psinm1 ) &
1138 /( ( cior * cbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
1139 qsca = qsca + twonp1 * ( sq( an ) + sq( bn ) )
1140 !
1141 end if
1142 ! ** save mie coefficients for *pmom* calculation
1143 lita( n ) = an
1144 litb( n ) = bn
1145 ! ** increment mie sums for non-angle-
1146 ! ** dependent quantities
1147 !
1148 sforw = sforw + twonp1 * ( an + bn )
1149 tforw( 1 ) = tforw( 1 ) + tcoef * ( an - bn )
1150 sback = sback + ( mm * twonp1 ) * ( an - bn )
1151 tback( 1 ) = tback( 1 ) + ( mm * tcoef ) * ( an + bn )
1152 gqsc = gqsc + ( fn - rn ) * dble( anm1 * dconjg( an ) &
1153 + bnm1 * dconjg( bn ) ) &
1154 + coeff * dble( an * dconjg( bn ) )
1155 !
1156 if ( yesang ) then
1157 ! ** put mie coefficients in form
1158 ! ** needed for computing s+, s-
1159 ! ** ( ref. 2, p. 1507 )
1160 anp = coeff * ( an + bn )
1161 bnp = coeff * ( an - bn )
1162 ! ** increment mie sums for s+, s-
1163 ! ** while upward recursing
1164 ! ** angular functions little pi
1165 ! ** and little tau
1166 if ( anyang ) then
1167 ! ** arbitrary angles
1168 !
1169 ! ** vectorizable loop
1170 do 80 j = 1, numang
1171 rtmp = ( xmu( j ) * pin( j ) ) - pinm1( j )
1172 taun = fn * rtmp - pinm1( j )
1173 sp( j ) = sp( j ) + anp * ( pin( j ) + taun )
1174 sm( j ) = sm( j ) + bnp * ( pin( j ) - taun )
1175 pinm1( j ) = pin( j )
1176 pin( j ) = ( xmu( j ) * pin( j ) ) + np1dn * rtmp
1177 80 continue
1178 !
1179 else
1180 ! ** angles symmetric about 90 degrees
1181 anpm = mm * anp
1182 bnpm = mm * bnp
1183 ! ** vectorizable loop
1184 do 90 j = 1, nangd2
1185 rtmp = ( xmu( j ) * pin( j ) ) - pinm1( j )
1186 taun = fn * rtmp - pinm1( j )
1187 sp ( j ) = sp ( j ) + anp * ( pin( j ) + taun )
1188 sms( j ) = sms( j ) + bnpm * ( pin( j ) + taun )
1189 sm ( j ) = sm ( j ) + bnp * ( pin( j ) - taun )
1190 sps( j ) = sps( j ) + anpm * ( pin( j ) - taun )
1191 pinm1( j ) = pin( j )
1192 pin( j ) = ( xmu( j ) * pin( j ) ) + np1dn * rtmp
1193 90 continue
1194 !
1195 end if
1196 end if
1197 ! ** update relevant quantities for next
1198 ! ** pass through loop
1199 mm = - mm
1200 anm1 = an
1201 bnm1 = bn
1202 ! ** upward recurrence for ricatti-bessel
1203 ! ** functions ( ref. 1, eq. 17 )
1204 !
1205 zet = ( twonp1 * xinv ) * zetn - zetnm1
1206 zetnm1 = zetn
1207 zetn = zet
1208 psinm1 = psin
1209 psin = dble( zetn )
1210 100 continue
1211 !
1212 ! ---------- end loop to sum mie series --------------------------------
1213 !
1214 !
1215 qext = 2. / xx**2 * dble( sforw )
1216 if ( perfct .or. noabs ) then
1217 qsca = qext
1218 else
1219 qsca = 2. / xx**2 * qsca
1220 end if
1221 !
1222 gqsc = 4. / xx**2 * gqsc
1223 sforw = 0.5 * sforw
1224 sback = 0.5 * sback
1225 tforw( 2 ) = 0.5 * ( sforw + 0.25 * tforw( 1 ) )
1226 tforw( 1 ) = 0.5 * ( sforw - 0.25 * tforw( 1 ) )
1227 tback( 2 ) = 0.5 * ( sback + 0.25 * tback( 1 ) )
1228 tback( 1 ) = 0.5 * ( - sback + 0.25 * tback( 1 ) )
1229 !
1230 if ( yesang ) then
1231 ! ** recover scattering amplitudes
1232 ! ** from s+, s- ( ref. 1, eq. 11 )
1233 if ( anyang ) then
1234 ! ** vectorizable loop
1235 do 110 j = 1, numang
1236 s1( j ) = 0.5 * ( sp( j ) + sm( j ) )
1237 s2( j ) = 0.5 * ( sp( j ) - sm( j ) )
1238 110 continue
1239 !
1240 else
1241 ! ** vectorizable loop
1242 do 120 j = 1, nangd2
1243 s1( j ) = 0.5 * ( sp( j ) + sm( j ) )
1244 s2( j ) = 0.5 * ( sp( j ) - sm( j ) )
1245 120 continue
1246 ! ** vectorizable loop
1247 do 130 j = 1, nangd2
1248 s1( numang+1 - j ) = 0.5 * ( sps( j ) + sms( j ) )
1249 s2( numang+1 - j ) = 0.5 * ( sps( j ) - sms( j ) )
1250 130 continue
1251 end if
1252 !
1253 end if
1254 ! ** calculate legendre moments
1255 200 if ( nmom.gt.0 ) &
1256 call lpcoef ( ntrm, nmom, ipolzn, momdim, calcmo, npquan, &
1257 lita, litb, pmom )
1258 !
1259 if ( dimag(crefin) .gt. 0.0 ) then
1260 ! ** take complex conjugates
1261 ! ** of scattering amplitudes
1262 sforw = dconjg( sforw )
1263 sback = dconjg( sback )
1264 do 210 i = 1, 2
1265 tforw( i ) = dconjg( tforw(i) )
1266 tback( i ) = dconjg( tback(i) )
1267 210 continue
1268 !
1269 do 220 j = 1, numang
1270 s1( j ) = dconjg( s1(j) )
1271 s2( j ) = dconjg( s2(j) )
1272 220 continue
1273 !
1274 end if
1275 !
1276 if ( pass1 ) then
1277 ! ** compare test case results with
1278 ! ** correct answers and abort if bad
1279 !
1280 call testmi ( qext, qsca, gqsc, sforw, sback, s1, s2, &
1281 tforw, tback, pmom, momdim, ok )
1282 if ( .not. ok ) then
1283 prnt(1) = .false.
1284 prnt(2) = .false.
1285 call miprnt( prnt, xx, perfct, crefin, numang, xmu, qext, &
1286 qsca, gqsc, nmom, ipolzn, momdim, calcmo, &
1287 pmom, sforw, sback, tforw, tback, s1, s2 )
1288 call errmsg( 'miev0 -- self-test failed', .true. )
1289 end if
1290 ! ** restore user input values
1291 xx = xxsav
1292 crefin = cresav
1293 mimcut = mimsav
1294 perfct = persav
1295 anyang = anysav
1296 nmom = nmosav
1297 ipolzn = iposav
1298 numang = numsav
1299 xmu(1) = xmusav
1300 pass1 = .false.
1301 go to 10
1302 !
1303 end if
1304 !
1305 if ( prnt(1) .or. prnt(2) ) &
1306 call miprnt( prnt, xx, perfct, crefin, numang, xmu, qext, &
1307 qsca, gqsc, nmom, ipolzn, momdim, calcmo, &
1308 pmom, sforw, sback, tforw, tback, s1, s2 )
1309 !
1310 return
1311 !
1312 end subroutine miev0
1313 !****************************************************************************
1314 subroutine ckinmi( numang, maxang, xx, perfct, crefin, momdim, &
1315 nmom, ipolzn, anyang, xmu, calcmo, npquan )
1316 !
1317 ! check for bad input to 'miev0' and calculate -calcmo,npquan-
1318 !
1319 implicit none
1320 logical perfct, anyang, calcmo(*)
1321 integer numang, maxang, momdim, nmom, ipolzn, npquan
1322 real*8 xx, xmu(*)
1323 integer i,l,j,ip
1324 complex*16 crefin
1325 !
1326 character*4 string
1327 logical inperr
1328 !
1329 inperr = .false.
1330 !
1331 if ( numang.gt.maxang ) then
1332 call errmsg( 'miev0--parameter maxang too small', .true. )
1333 inperr = .true.
1334 end if
1335 if ( numang.lt.0 ) call wrtbad( 'numang', inperr )
1336 if ( xx.lt.0. ) call wrtbad( 'xx', inperr )
1337 if ( .not.perfct .and. dble(crefin).le.0. ) &
1338 call wrtbad( 'crefin', inperr )
1339 if ( momdim.lt.1 ) call wrtbad( 'momdim', inperr )
1340 !
1341 if ( nmom.ne.0 ) then
1342 if ( nmom.lt.0 .or. nmom.gt.momdim ) call wrtbad('nmom',inperr)
1343 if ( iabs(ipolzn).gt.4444 ) call wrtbad( 'ipolzn', inperr )
1344 npquan = 0
1345 do 5 l = 1, 4
1346 calcmo( l ) = .false.
1347 5 continue
1348 if ( ipolzn.ne.0 ) then
1349 ! ** parse out -ipolzn- into its digits
1350 ! ** to find which phase quantities are
1351 ! ** to have their moments calculated
1352 !
1353 write( string, '(i4)' ) iabs(ipolzn)
1354 do 10 j = 1, 4
1355 ip = ichar( string(j:j) ) - ichar( '0' )
1356 if ( ip.ge.1 .and. ip.le.4 ) calcmo( ip ) = .true.
1357 if ( ip.eq.0 .or. (ip.ge.5 .and. ip.le.9) ) &
1358 call wrtbad( 'ipolzn', inperr )
1359 npquan = max0( npquan, ip )
1360 10 continue
1361 end if
1362 end if
1363 !
1364 if ( anyang ) then
1365 ! ** allow for slight imperfections in
1366 ! ** computation of cosine
1367 do 20 i = 1, numang
1368 if ( xmu(i).lt.-1.00001 .or. xmu(i).gt.1.00001 ) &
1369 call wrtbad( 'xmu', inperr )
1370 20 continue
1371 else
1372 do 22 i = 1, ( numang + 1 ) / 2
1373 if ( xmu(i).lt.-0.00001 .or. xmu(i).gt.1.00001 ) &
1374 call wrtbad( 'xmu', inperr )
1375 22 continue
1376 end if
1377 !
1378 if ( inperr ) &
1379 call errmsg( 'miev0--input error(s). aborting...', .true. )
1380 !
1381 if ( xx.gt.20000.0 .or. dble(crefin).gt.10.0 .or. &
1382 dabs( dimag(crefin) ).gt.10.0 ) call errmsg( &
1383 'miev0--xx or crefin outside tested range', .false. )
1384 !
1385 return
1386 end subroutine ckinmi
1387 !***********************************************************************
1388 subroutine lpcoef ( ntrm, nmom, ipolzn, momdim, calcmo, npquan, &
1389 a, b, pmom )
1390 !
1391 ! calculate legendre polynomial expansion coefficients (also
1392 ! called moments) for phase quantities ( ref. 5 formulation )
1393 !
1394 ! input: ntrm number terms in mie series
1395 ! nmom, ipolzn, momdim 'miev0' arguments
1396 ! calcmo flags calculated from -ipolzn-
1397 ! npquan defined in 'miev0'
1398 ! a, b mie series coefficients
1399 !
1400 ! output: pmom legendre moments ('miev0' argument)
1401 !
1402 ! *** notes ***
1403 !
1404 ! (1) eqs. 2-5 are in error in dave, appl. opt. 9,
1405 ! 1888 (1970). eq. 2 refers to m1, not m2; eq. 3 refers to
1406 ! m2, not m1. in eqs. 4 and 5, the subscripts on the second
1407 ! term in square brackets should be interchanged.
1408 !
1409 ! (2) the general-case logic in this subroutine works correctly
1410 ! in the two-term mie series case, but subroutine 'lpco2t'
1411 ! is called instead, for speed.
1412 !
1413 ! (3) subroutine 'lpco1t', to do the one-term case, is never
1414 ! called within the context of 'miev0', but is included for
1415 ! complete generality.
1416 !
1417 ! (4) some improvement in speed is obtainable by combining the
1418 ! 310- and 410-loops, if moments for both the third and fourth
1419 ! phase quantities are desired, because the third phase quantity
1420 ! is the real part of a complex series, while the fourth phase
1421 ! quantity is the imaginary part of that very same series. but
1422 ! most users are not interested in the fourth phase quantity,
1423 ! which is related to circular polarization, so the present
1424 ! scheme is usually more efficient.
1425 !
1426 implicit none
1427 logical calcmo(*)
1428 integer ipolzn, momdim, nmom, ntrm, npquan
1429 real*8 pmom( 0:momdim, * )
1430 complex*16 a(*), b(*)
1431 !
1432 ! ** specification of local variables
1433 !
1434 ! am(m) numerical coefficients a-sub-m-super-l
1435 ! in dave, eqs. 1-15, as simplified in ref. 5.
1436 !
1437 ! bi(i) numerical coefficients b-sub-i-super-l
1438 ! in dave, eqs. 1-15, as simplified in ref. 5.
1439 !
1440 ! bidel(i) 1/2 bi(i) times factor capital-del in dave
1441 !
1442 ! cm,dm() arrays c and d in dave, eqs. 16-17 (mueller form),
1443 ! calculated using recurrence derived in ref. 5
1444 !
1445 ! cs,ds() arrays c and d in ref. 4, eqs. a5-a6 (sekera form),
1446 ! calculated using recurrence derived in ref. 5
1447 !
1448 ! c,d() either -cm,dm- or -cs,ds-, depending on -ipolzn-
1449 !
1450 ! evenl true for even-numbered moments; false otherwise
1451 !
1452 ! idel 1 + little-del in dave
1453 !
1454 ! maxtrm max. no. of terms in mie series
1455 !
1456 ! maxmom max. no. of non-zero moments
1457 !
1458 ! nummom number of non-zero moments
1459 !
1460 ! recip(k) 1 / k
1461 !
1462 integer maxtrm,maxmom,mxmom2,maxrcp
1463 parameter ( maxtrm = 1102, maxmom = 2*maxtrm, mxmom2 = maxmom/2, &
1464 maxrcp = 4*maxtrm + 2 )
1465 real*8 am( 0:maxtrm ), bi( 0:mxmom2 ), bidel( 0:mxmom2 ), &
1466 recip( maxrcp )
1467 complex*16 cm( maxtrm ), dm( maxtrm ), cs( maxtrm ), ds( maxtrm ), &
1468 c( maxtrm ), d( maxtrm )
1469 integer k,j,l,nummom,ld2,idel,m,i,mmax,imax
1470 real*8 thesum
1471 equivalence ( c, cm ), ( d, dm )
1472 logical pass1, evenl
1473 save pass1, recip
1474 data pass1 / .true. /
1475 !
1476 !
1477 if ( pass1 ) then
1478 !
1479 do 1 k = 1, maxrcp
1480 recip( k ) = 1.0 / k
1481 1 continue
1482 pass1 = .false.
1483 !
1484 end if
1485 !
1486 do 5 j = 1, max0( 1, npquan )
1487 do 5 l = 0, nmom
1488 pmom( l, j ) = 0.0
1489 5 continue
1490 !
1491 if ( ntrm.eq.1 ) then
1492 call lpco1t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
1493 return
1494 else if ( ntrm.eq.2 ) then
1495 call lpco2t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
1496 return
1497 end if
1498 !
1499 if ( ntrm+2 .gt. maxtrm ) &
1500 call errmsg( 'lpcoef--parameter maxtrm too small', .true. )
1501 !
1502 ! ** calculate mueller c, d arrays
1503 cm( ntrm+2 ) = ( 0., 0. )
1504 dm( ntrm+2 ) = ( 0., 0. )
1505 cm( ntrm+1 ) = ( 1. - recip( ntrm+1 ) ) * b( ntrm )
1506 dm( ntrm+1 ) = ( 1. - recip( ntrm+1 ) ) * a( ntrm )
1507 cm( ntrm ) = ( recip(ntrm) + recip(ntrm+1) ) * a( ntrm ) &
1508 + ( 1. - recip(ntrm) ) * b( ntrm-1 )
1509 dm( ntrm ) = ( recip(ntrm) + recip(ntrm+1) ) * b( ntrm ) &
1510 + ( 1. - recip(ntrm) ) * a( ntrm-1 )
1511 !
1512 do 10 k = ntrm-1, 2, -1
1513 cm( k ) = cm( k+2 ) - ( 1. + recip(k+1) ) * b( k+1 ) &
1514 + ( recip(k) + recip(k+1) ) * a( k ) &
1515 + ( 1. - recip(k) ) * b( k-1 )
1516 dm( k ) = dm( k+2 ) - ( 1. + recip(k+1) ) * a( k+1 ) &
1517 + ( recip(k) + recip(k+1) ) * b( k ) &
1518 + ( 1. - recip(k) ) * a( k-1 )
1519 10 continue
1520 cm( 1 ) = cm( 3 ) + 1.5 * ( a( 1 ) - b( 2 ) )
1521 dm( 1 ) = dm( 3 ) + 1.5 * ( b( 1 ) - a( 2 ) )
1522 !
1523 if ( ipolzn.ge.0 ) then
1524 !
1525 do 20 k = 1, ntrm + 2
1526 c( k ) = ( 2*k - 1 ) * cm( k )
1527 d( k ) = ( 2*k - 1 ) * dm( k )
1528 20 continue
1529 !
1530 else
1531 ! ** compute sekera c and d arrays
1532 cs( ntrm+2 ) = ( 0., 0. )
1533 ds( ntrm+2 ) = ( 0., 0. )
1534 cs( ntrm+1 ) = ( 0., 0. )
1535 ds( ntrm+1 ) = ( 0., 0. )
1536 !
1537 do 30 k = ntrm, 1, -1
1538 cs( k ) = cs( k+2 ) + ( 2*k + 1 ) * ( cm( k+1 ) - b( k ) )
1539 ds( k ) = ds( k+2 ) + ( 2*k + 1 ) * ( dm( k+1 ) - a( k ) )
1540 30 continue
1541 !
1542 do 40 k = 1, ntrm + 2
1543 c( k ) = ( 2*k - 1 ) * cs( k )
1544 d( k ) = ( 2*k - 1 ) * ds( k )
1545 40 continue
1546 !
1547 end if
1548 !
1549 !
1550 if( ipolzn.lt.0 ) nummom = min0( nmom, 2*ntrm - 2 )
1551 if( ipolzn.ge.0 ) nummom = min0( nmom, 2*ntrm )
1552 if ( nummom .gt. maxmom ) &
1553 call errmsg( 'lpcoef--parameter maxtrm too small', .true. )
1554 !
1555 ! ** loop over moments
1556 do 500 l = 0, nummom
1557 ld2 = l / 2
1558 evenl = mod( l,2 ) .eq. 0
1559 ! ** calculate numerical coefficients
1560 ! ** a-sub-m and b-sub-i in dave
1561 ! ** double-sums for moments
1562 if( l.eq.0 ) then
1563 !
1564 idel = 1
1565 do 60 m = 0, ntrm
1566 am( m ) = 2.0 * recip( 2*m + 1 )
1567 60 continue
1568 bi( 0 ) = 1.0
1569 !
1570 else if( evenl ) then
1571 !
1572 idel = 1
1573 do 70 m = ld2, ntrm
1574 am( m ) = ( 1. + recip( 2*m-l+1 ) ) * am( m )
1575 70 continue
1576 do 75 i = 0, ld2-1
1577 bi( i ) = ( 1. - recip( l-2*i ) ) * bi( i )
1578 75 continue
1579 bi( ld2 ) = ( 2. - recip( l ) ) * bi( ld2-1 )
1580 !
1581 else
1582 !
1583 idel = 2
1584 do 80 m = ld2, ntrm
1585 am( m ) = ( 1. - recip( 2*m+l+2 ) ) * am( m )
1586 80 continue
1587 do 85 i = 0, ld2
1588 bi( i ) = ( 1. - recip( l+2*i+1 ) ) * bi( i )
1589 85 continue
1590 !
1591 end if
1592 ! ** establish upper limits for sums
1593 ! ** and incorporate factor capital-
1594 ! ** del into b-sub-i
1595 mmax = ntrm - idel
1596 if( ipolzn.ge.0 ) mmax = mmax + 1
1597 imax = min0( ld2, mmax - ld2 )
1598 if( imax.lt.0 ) go to 600
1599 do 90 i = 0, imax
1600 bidel( i ) = bi( i )
1601 90 continue
1602 if( evenl ) bidel( 0 ) = 0.5 * bidel( 0 )
1603 !
1604 ! ** perform double sums just for
1605 ! ** phase quantities desired by user
1606 if( ipolzn.eq.0 ) then
1607 !
1608 do 110 i = 0, imax
1609 ! ** vectorizable loop (cray)
1610 thesum = 0.0
1611 do 100 m = ld2, mmax - i
1612 thesum = thesum + am( m ) * &
1613 ( dble( c(m-i+1) * dconjg( c(m+i+idel) ) ) &
1614 + dble( d(m-i+1) * dconjg( d(m+i+idel) ) ) )
1615 100 continue
1616 pmom( l,1 ) = pmom( l,1 ) + bidel( i ) * thesum
1617 110 continue
1618 pmom( l,1 ) = 0.5 * pmom( l,1 )
1619 go to 500
1620 !
1621 end if
1622 !
1623 if ( calcmo(1) ) then
1624 do 160 i = 0, imax
1625 ! ** vectorizable loop (cray)
1626 thesum = 0.0
1627 do 150 m = ld2, mmax - i
1628 thesum = thesum + am( m ) * &
1629 dble( c(m-i+1) * dconjg( c(m+i+idel) ) )
1630 150 continue
1631 pmom( l,1 ) = pmom( l,1 ) + bidel( i ) * thesum
1632 160 continue
1633 end if
1634 !
1635 !
1636 if ( calcmo(2) ) then
1637 do 210 i = 0, imax
1638 ! ** vectorizable loop (cray)
1639 thesum = 0.0
1640 do 200 m = ld2, mmax - i
1641 thesum = thesum + am( m ) * &
1642 dble( d(m-i+1) * dconjg( d(m+i+idel) ) )
1643 200 continue
1644 pmom( l,2 ) = pmom( l,2 ) + bidel( i ) * thesum
1645 210 continue
1646 end if
1647 !
1648 !
1649 if ( calcmo(3) ) then
1650 do 310 i = 0, imax
1651 ! ** vectorizable loop (cray)
1652 thesum = 0.0
1653 do 300 m = ld2, mmax - i
1654 thesum = thesum + am( m ) * &
1655 ( dble( c(m-i+1) * dconjg( d(m+i+idel) ) ) &
1656 + dble( c(m+i+idel) * dconjg( d(m-i+1) ) ) )
1657 300 continue
1658 pmom( l,3 ) = pmom( l,3 ) + bidel( i ) * thesum
1659 310 continue
1660 pmom( l,3 ) = 0.5 * pmom( l,3 )
1661 end if
1662 !
1663 !
1664 if ( calcmo(4) ) then
1665 do 410 i = 0, imax
1666 ! ** vectorizable loop (cray)
1667 thesum = 0.0
1668 do 400 m = ld2, mmax - i
1669 thesum = thesum + am( m ) * &
1670 ( dimag( c(m-i+1) * dconjg( d(m+i+idel) ) ) &
1671 + dimag( c(m+i+idel) * dconjg( d(m-i+1) ) ))
1672 400 continue
1673 pmom( l,4 ) = pmom( l,4 ) + bidel( i ) * thesum
1674 410 continue
1675 pmom( l,4 ) = - 0.5 * pmom( l,4 )
1676 end if
1677 !
1678 500 continue
1679 !
1680 !
1681 600 return
1682 end subroutine lpcoef
1683 !*********************************************************************
1684 subroutine lpco1t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
1685 !
1686 ! calculate legendre polynomial expansion coefficients (also
1687 ! called moments) for phase quantities in special case where
1688 ! no. terms in mie series = 1
1689 !
1690 ! input: nmom, ipolzn, momdim 'miev0' arguments
1691 ! calcmo flags calculated from -ipolzn-
1692 ! a(1), b(1) mie series coefficients
1693 !
1694 ! output: pmom legendre moments
1695 !
1696 implicit none
1697 logical calcmo(*)
1698 integer ipolzn, momdim, nmom,nummom,l
1699 real*8 pmom( 0:momdim, * ),sq,a1sq,b1sq
1700 complex*16 a(*), b(*), ctmp, a1b1c
1701 sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
1702 !
1703 !
1704 a1sq = sq( a(1) )
1705 b1sq = sq( b(1) )
1706 a1b1c = a(1) * dconjg( b(1) )
1707 !
1708 if( ipolzn.lt.0 ) then
1709 !
1710 if( calcmo(1) ) pmom( 0,1 ) = 2.25 * b1sq
1711 if( calcmo(2) ) pmom( 0,2 ) = 2.25 * a1sq
1712 if( calcmo(3) ) pmom( 0,3 ) = 2.25 * dble( a1b1c )
1713 if( calcmo(4) ) pmom( 0,4 ) = 2.25 *dimag( a1b1c )
1714 !
1715 else
1716 !
1717 nummom = min0( nmom, 2 )
1718 ! ** loop over moments
1719 do 100 l = 0, nummom
1720 !
1721 if( ipolzn.eq.0 ) then
1722 if( l.eq.0 ) pmom( l,1 ) = 1.5 * ( a1sq + b1sq )
1723 if( l.eq.1 ) pmom( l,1 ) = 1.5 * dble( a1b1c )
1724 if( l.eq.2 ) pmom( l,1 ) = 0.15 * ( a1sq + b1sq )
1725 go to 100
1726 end if
1727 !
1728 if( calcmo(1) ) then
1729 if( l.eq.0 ) pmom( l,1 ) = 2.25 * ( a1sq + b1sq / 3. )
1730 if( l.eq.1 ) pmom( l,1 ) = 1.5 * dble( a1b1c )
1731 if( l.eq.2 ) pmom( l,1 ) = 0.3 * b1sq
1732 end if
1733 !
1734 if( calcmo(2) ) then
1735 if( l.eq.0 ) pmom( l,2 ) = 2.25 * ( b1sq + a1sq / 3. )
1736 if( l.eq.1 ) pmom( l,2 ) = 1.5 * dble( a1b1c )
1737 if( l.eq.2 ) pmom( l,2 ) = 0.3 * a1sq
1738 end if
1739 !
1740 if( calcmo(3) ) then
1741 if( l.eq.0 ) pmom( l,3 ) = 3.0 * dble( a1b1c )
1742 if( l.eq.1 ) pmom( l,3 ) = 0.75 * ( a1sq + b1sq )
1743 if( l.eq.2 ) pmom( l,3 ) = 0.3 * dble( a1b1c )
1744 end if
1745 !
1746 if( calcmo(4) ) then
1747 if( l.eq.0 ) pmom( l,4 ) = - 1.5 * dimag( a1b1c )
1748 if( l.eq.1 ) pmom( l,4 ) = 0.0
1749 if( l.eq.2 ) pmom( l,4 ) = 0.3 * dimag( a1b1c )
1750 end if
1751 !
1752 100 continue
1753 !
1754 end if
1755 !
1756 return
1757 end subroutine lpco1t
1758 !********************************************************************
1759 subroutine lpco2t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
1760 !
1761 ! calculate legendre polynomial expansion coefficients (also
1762 ! called moments) for phase quantities in special case where
1763 ! no. terms in mie series = 2
1764 !
1765 ! input: nmom, ipolzn, momdim 'miev0' arguments
1766 ! calcmo flags calculated from -ipolzn-
1767 ! a(1-2), b(1-2) mie series coefficients
1768 !
1769 ! output: pmom legendre moments
1770 !
1771 implicit none
1772 logical calcmo(*)
1773 integer ipolzn, momdim, nmom,l,nummom
1774 real*8 pmom( 0:momdim, * ),sq,pm1,pm2,a2sq,b2sq
1775 complex*16 a(*), b(*)
1776 complex*16 a2c, b2c, ctmp, ca, cac, cat, cb, cbc, cbt, cg, ch
1777 sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
1778 !
1779 !
1780 ca = 3. * a(1) - 5. * b(2)
1781 cat= 3. * b(1) - 5. * a(2)
1782 cac = dconjg( ca )
1783 a2sq = sq( a(2) )
1784 b2sq = sq( b(2) )
1785 a2c = dconjg( a(2) )
1786 b2c = dconjg( b(2) )
1787 !
1788 if( ipolzn.lt.0 ) then
1789 ! ** loop over sekera moments
1790 nummom = min0( nmom, 2 )
1791 do 50 l = 0, nummom
1792 !
1793 if( calcmo(1) ) then
1794 if( l.eq.0 ) pmom( l,1 ) = 0.25 * ( sq(cat) + &
1795 (100./3.) * b2sq )
1796 if( l.eq.1 ) pmom( l,1 ) = (5./3.) * dble( cat * b2c )
1797 if( l.eq.2 ) pmom( l,1 ) = (10./3.) * b2sq
1798 end if
1799 !
1800 if( calcmo(2) ) then
1801 if( l.eq.0 ) pmom( l,2 ) = 0.25 * ( sq(ca) + &
1802 (100./3.) * a2sq )
1803 if( l.eq.1 ) pmom( l,2 ) = (5./3.) * dble( ca * a2c )
1804 if( l.eq.2 ) pmom( l,2 ) = (10./3.) * a2sq
1805 end if
1806 !
1807 if( calcmo(3) ) then
1808 if( l.eq.0 ) pmom( l,3 ) = 0.25 * dble( cat*cac + &
1809 (100./3.)*b(2)*a2c )
1810 if( l.eq.1 ) pmom( l,3 ) = 5./6. * dble( b(2)*cac + &
1811 cat*a2c )
1812 if( l.eq.2 ) pmom( l,3 ) = 10./3. * dble( b(2) * a2c )
1813 end if
1814 !
1815 if( calcmo(4) ) then
1816 if( l.eq.0 ) pmom( l,4 ) = -0.25 * dimag( cat*cac + &
1817 (100./3.)*b(2)*a2c )
1818 if( l.eq.1 ) pmom( l,4 ) = -5./6. * dimag( b(2)*cac + &
1819 cat*a2c )
1820 if( l.eq.2 ) pmom( l,4 ) = -10./3. * dimag( b(2) * a2c )
1821 end if
1822 !
1823 50 continue
1824 !
1825 else
1826 !
1827 cb = 3. * b(1) + 5. * a(2)
1828 cbt= 3. * a(1) + 5. * b(2)
1829 cbc = dconjg( cb )
1830 cg = ( cbc*cbt + 10.*( cac*a(2) + b2c*cat) ) / 3.
1831 ch = 2.*( cbc*a(2) + b2c*cbt )
1832 !
1833 ! ** loop over mueller moments
1834 nummom = min0( nmom, 4 )
1835 do 100 l = 0, nummom
1836 !
1837 if( ipolzn.eq.0 .or. calcmo(1) ) then
1838 if( l.eq.0 ) pm1 = 0.25 * sq(ca) + sq(cb) / 12. &
1839 + (5./3.) * dble(ca*b2c) + 5.*b2sq
1840 if( l.eq.1 ) pm1 = dble( cb * ( cac/6. + b2c ) )
1841 if( l.eq.2 ) pm1 = sq(cb)/30. + (20./7.) * b2sq &
1842 + (2./3.) * dble( ca * b2c )
1843 if( l.eq.3 ) pm1 = (2./7.) * dble( cb * b2c )
1844 if( l.eq.4 ) pm1 = (40./63.) * b2sq
1845 if ( calcmo(1) ) pmom( l,1 ) = pm1
1846 end if
1847 !
1848 if( ipolzn.eq.0 .or. calcmo(2) ) then
1849 if( l.eq.0 ) pm2 = 0.25*sq(cat) + sq(cbt) / 12. &
1850 + (5./3.) * dble(cat*a2c) + 5.*a2sq
1851 if( l.eq.1 ) pm2 = dble( cbt * ( dconjg(cat)/6. + a2c) )
1852 if( l.eq.2 ) pm2 = sq(cbt)/30. + (20./7.) * a2sq &
1853 + (2./3.) * dble( cat * a2c )
1854 if( l.eq.3 ) pm2 = (2./7.) * dble( cbt * a2c )
1855 if( l.eq.4 ) pm2 = (40./63.) * a2sq
1856 if ( calcmo(2) ) pmom( l,2 ) = pm2
1857 end if
1858 !
1859 if( ipolzn.eq.0 ) then
1860 pmom( l,1 ) = 0.5 * ( pm1 + pm2 )
1861 go to 100
1862 end if
1863 !
1864 if( calcmo(3) ) then
1865 if( l.eq.0 ) pmom( l,3 ) = 0.25 * dble( cac*cat + cg + &
1866 20.*b2c*a(2) )
1867 if( l.eq.1 ) pmom( l,3 ) = dble( cac*cbt + cbc*cat + &
1868 3.*ch ) / 12.
1869 if( l.eq.2 ) pmom( l,3 ) = 0.1 * dble( cg + (200./7.) * &
1870 b2c * a(2) )
1871 if( l.eq.3 ) pmom( l,3 ) = dble( ch ) / 14.
1872 if( l.eq.4 ) pmom( l,3 ) = 40./63. * dble( b2c * a(2) )
1873 end if
1874 !
1875 if( calcmo(4) ) then
1876 if( l.eq.0 ) pmom( l,4 ) = 0.25 * dimag( cac*cat + cg + &
1877 20.*b2c*a(2) )
1878 if( l.eq.1 ) pmom( l,4 ) = dimag( cac*cbt + cbc*cat + &
1879 3.*ch ) / 12.
1880 if( l.eq.2 ) pmom( l,4 ) = 0.1 * dimag( cg + (200./7.) * &
1881 b2c * a(2) )
1882 if( l.eq.3 ) pmom( l,4 ) = dimag( ch ) / 14.
1883 if( l.eq.4 ) pmom( l,4 ) = 40./63. * dimag( b2c * a(2) )
1884 end if
1885 !
1886 100 continue
1887 !
1888 end if
1889 !
1890 return
1891 end subroutine lpco2t
1892 !*********************************************************************
1893 subroutine biga( cior, xx, ntrm, noabs, yesang, rbiga, cbiga )
1894 !
1895 ! calculate logarithmic derivatives of j-bessel-function
1896 !
1897 ! input : cior, xx, ntrm, noabs, yesang (defined in 'miev0')
1898 !
1899 ! output : rbiga or cbiga (defined in 'miev0')
1900 !
1901 ! internal variables :
1902 !
1903 ! confra value of lentz continued fraction for -cbiga(ntrm)-,
1904 ! used to initialize downward recurrence.
1905 ! down = true, use down-recurrence. false, do not.
1906 ! f1,f2,f3 arithmetic statement functions used in determining
1907 ! whether to use up- or down-recurrence
1908 ! ( ref. 2, eqs. 6-8 )
1909 ! mre real refractive index
1910 ! mim imaginary refractive index
1911 ! rezinv 1 / ( mre * xx ); temporary variable for recurrence
1912 ! zinv 1 / ( cior * xx ); temporary variable for recurrence
1913 !
1914 implicit none
1915 logical down, noabs, yesang
1916 integer ntrm,n
1917 real*8 mre, mim, rbiga(*), xx, rezinv, rtmp, f1,f2,f3
1918 ! complex*16 cior, ctmp, confra, cbiga(*), zinv
1919 complex*16 cior, ctmp, cbiga(*), zinv
1920 f1( mre ) = - 8.0 + mre**2 * ( 26.22 + mre * ( - 0.4474 &
1921 + mre**3 * ( 0.00204 - 0.000175 * mre ) ) )
1922 f2( mre ) = 3.9 + mre * ( - 10.8 + 13.78 * mre )
1923 f3( mre ) = - 15.04 + mre * ( 8.42 + 16.35 * mre )
1924 !
1925 ! ** decide whether 'biga' can be
1926 ! ** calculated by up-recurrence
1927 mre = dble( cior )
1928 mim = dabs( dimag( cior ) )
1929 if ( mre.lt.1.0 .or. mre.gt.10.0 .or. mim.gt.10.0 ) then
1930 down = .true.
1931 else if ( yesang ) then
1932 down = .true.
1933 if ( mim*xx .lt. f2( mre ) ) down = .false.
1934 else
1935 down = .true.
1936 if ( mim*xx .lt. f1( mre ) ) down = .false.
1937 end if
1938 !
1939 zinv = 1.0 / ( cior * xx )
1940 rezinv = 1.0 / ( mre * xx )
1941 !
1942 if ( down ) then
1943 ! ** compute initial high-order 'biga' using
1944 ! ** lentz method ( ref. 1, pp. 17-20 )
1945 !
1946 ctmp = confra( ntrm, zinv, xx )
1947 !
1948 ! *** downward recurrence for 'biga'
1949 ! *** ( ref. 1, eq. 22 )
1950 if ( noabs ) then
1951 ! ** no-absorption case
1952 rbiga( ntrm ) = dble( ctmp )
1953 do 25 n = ntrm, 2, - 1
1954 rbiga( n-1 ) = (n*rezinv) &
1955 - 1.0 / ( (n*rezinv) + rbiga( n ) )
1956 25 continue
1957 !
1958 else
1959 ! ** absorptive case
1960 cbiga( ntrm ) = ctmp
1961 do 30 n = ntrm, 2, - 1
1962 cbiga( n-1 ) = (n*zinv) - 1.0 / ( (n*zinv) + cbiga( n ) )
1963 30 continue
1964 !
1965 end if
1966 !
1967 else
1968 ! *** upward recurrence for 'biga'
1969 ! *** ( ref. 1, eqs. 20-21 )
1970 if ( noabs ) then
1971 ! ** no-absorption case
1972 rtmp = dsin( mre*xx )
1973 rbiga( 1 ) = - rezinv &
1974 + rtmp / ( rtmp*rezinv - dcos( mre*xx ) )
1975 do 40 n = 2, ntrm
1976 rbiga( n ) = - ( n*rezinv ) &
1977 + 1.0 / ( ( n*rezinv ) - rbiga( n-1 ) )
1978 40 continue
1979 !
1980 else
1981 ! ** absorptive case
1982 ctmp = cdexp( - dcmplx(0.d0,2.d0) * cior * xx )
1983 cbiga( 1 ) = - zinv + (1.-ctmp) / ( zinv * (1.-ctmp) - &
1984 dcmplx(0.d0,1.d0)*(1.+ctmp) )
1985 do 50 n = 2, ntrm
1986 cbiga( n ) = - (n*zinv) + 1.0 / ((n*zinv) - cbiga( n-1 ))
1987 50 continue
1988 end if
1989 !
1990 end if
1991 !
1992 return
1993 end subroutine biga
1994 !**********************************************************************
1995 complex*16 function confra( n, zinv, xx )
1996 !
1997 ! compute bessel function ratio capital-a-sub-n from its
1998 ! continued fraction using lentz method ( ref. 1, pp. 17-20 )
1999 !
2000 ! zinv = reciprocal of argument of capital-a
2001 !
2002 ! i n t e r n a l v a r i a b l e s
2003 ! ------------------------------------
2004 !
2005 ! cak term in continued fraction expansion of capital-a
2006 ! ( ref. 1, eq. 25 )
2007 ! capt factor used in lentz iteration for capital-a
2008 ! ( ref. 1, eq. 27 )
2009 ! cdenom denominator in -capt- ( ref. 1, eq. 28b )
2010 ! cnumer numerator in -capt- ( ref. 1, eq. 28a )
2011 ! cdtd product of two successive denominators of -capt-
2012 ! factors ( ref. 1, eq. 34c )
2013 ! cntn product of two successive numerators of -capt-
2014 ! factors ( ref. 1, eq. 34b )
2015 ! eps1 ill-conditioning criterion
2016 ! eps2 convergence criterion
2017 ! kk subscript k of -cak- ( ref. 1, eq. 25b )
2018 ! kount iteration counter ( used only to prevent runaway )
2019 ! maxit max. allowed no. of iterations
2020 ! mm + 1 and - 1, alternately
2021 !
2022 implicit none
2023 integer n,maxit,mm,kk,kount
2024 real*8 xx,eps1,eps2
2025 complex*16 zinv
2026 complex*16 cak, capt, cdenom, cdtd, cnumer, cntn
2027 ! data eps1 / 1.e - 2 /, eps2 / 1.e - 8 /
2028 data eps1 / 1.d-2 /, eps2 / 1.d-8 /
2029 data maxit / 10000 /
2030 !
2031 ! *** ref. 1, eqs. 25a, 27
2032 confra = ( n + 1 ) * zinv
2033 mm = - 1
2034 kk = 2 * n + 3
2035 cak = ( mm * kk ) * zinv
2036 cdenom = cak
2037 cnumer = cdenom + 1.0 / confra
2038 kount = 1
2039 !
2040 20 kount = kount + 1
2041 if ( kount.gt.maxit ) &
2042 call errmsg( 'confra--iteration failed to converge$', .true.)
2043 !
2044 ! *** ref. 2, eq. 25b
2045 mm = - mm
2046 kk = kk + 2
2047 cak = ( mm * kk ) * zinv
2048 ! *** ref. 2, eq. 32
2049 if ( cdabs( cnumer/cak ).le.eps1 &
2050 .or. cdabs( cdenom/cak ).le.eps1 ) then
2051 !
2052 ! ** ill-conditioned case -- stride
2053 ! ** two terms instead of one
2054 !
2055 ! *** ref. 2, eqs. 34
2056 cntn = cak * cnumer + 1.0
2057 cdtd = cak * cdenom + 1.0
2058 confra = ( cntn / cdtd ) * confra
2059 ! *** ref. 2, eq. 25b
2060 mm = - mm
2061 kk = kk + 2
2062 cak = ( mm * kk ) * zinv
2063 ! *** ref. 2, eqs. 35
2064 cnumer = cak + cnumer / cntn
2065 cdenom = cak + cdenom / cdtd
2066 kount = kount + 1
2067 go to 20
2068 !
2069 else
2070 ! ** well-conditioned case
2071 !
2072 ! *** ref. 2, eqs. 26, 27
2073 capt = cnumer / cdenom
2074 confra = capt * confra
2075 ! ** check for convergence
2076 ! ** ( ref. 2, eq. 31 )
2077 !
2078 if ( dabs( dble(capt) - 1.0 ).ge.eps2 &
2079 .or. dabs( dimag(capt) ) .ge.eps2 ) then
2080 !
2081 ! *** ref. 2, eqs. 30a-b
2082 cnumer = cak + 1.0 / cnumer
2083 cdenom = cak + 1.0 / cdenom
2084 go to 20
2085 end if
2086 end if
2087 !
2088 return
2089 !
2090 end function confra
2091 !********************************************************************
2092 subroutine miprnt( prnt, xx, perfct, crefin, numang, xmu, &
2093 qext, qsca, gqsc, nmom, ipolzn, momdim, &
2094 calcmo, pmom, sforw, sback, tforw, tback, &
2095 s1, s2 )
2096 !
2097 ! print scattering quantities of a single particle
2098 !
2099 implicit none
2100 logical perfct, prnt(*), calcmo(*)
2101 integer ipolzn, momdim, nmom, numang,i,m,j
2102 real*8 gqsc, pmom( 0:momdim, * ), qext, qsca, xx, xmu(*)
2103 real*8 fi1,fi2,fnorm
2104 complex*16 crefin, sforw, sback, tforw(*), tback(*), s1(*), s2(*)
2105 character*22 fmt
2106 !
2107 !
2108 if ( perfct ) write ( *, '(''1'',10x,a,1p,e11.4)' ) &
2109 'perfectly conducting case, size parameter =', xx
2110 if ( .not.perfct ) write ( *, '(''1'',10x,3(a,1p,e11.4))' ) &
2111 'refractive index: real ', dble(crefin), &
2112 ' imag ', dimag(crefin), ', size parameter =', xx
2113 !
2114 if ( prnt(1) .and. numang.gt.0 ) then
2115 !
2116 write ( *, '(/,a)' ) &
2117 ' cos(angle) ------- s1 --------- ------- s2 ---------'// &
2118 ' --- s1*conjg(s2) --- i1=s1**2 i2=s2**2 (i1+i2)/2'// &
2119 ' deg polzn'
2120 do 10 i = 1, numang
2121 fi1 = dble( s1(i) ) **2 + dimag( s1(i) )**2
2122 fi2 = dble( s2(i) ) **2 + dimag( s2(i) )**2
2123 write( *, '( i4, f10.6, 1p,10e11.3 )' ) &
2124 i, xmu(i), s1(i), s2(i), s1(i)*dconjg(s2(i)), &
2125 fi1, fi2, 0.5*(fi1+fi2), (fi2-fi1)/(fi2+fi1)
2126 10 continue
2127 !
2128 end if
2129 !
2130 !
2131 if ( prnt(2) ) then
2132 !
2133 write ( *, '(/,a,9x,a,17x,a,17x,a,/,(0p,f7.2, 1p,6e12.3) )' ) &
2134 ' angle', 's-sub-1', 't-sub-1', 't-sub-2', &
2135 0.0, sforw, tforw(1), tforw(2), &
2136 180., sback, tback(1), tback(2)
2137 write ( *, '(/,4(a,1p,e11.4))' ) &
2138 ' efficiency factors, extinction:', qext, &
2139 ' scattering:', qsca, &
2140 ' absorption:', qext-qsca, &
2141 ' rad. pressure:', qext-gqsc
2142 !
2143 if ( nmom.gt.0 ) then
2144 !
2145 write( *, '(/,a)' ) ' normalized moments of :'
2146 if ( ipolzn.eq.0 ) write ( *, '(''+'',27x,a)' ) 'phase fcn'
2147 if ( ipolzn.gt.0 ) write ( *, '(''+'',33x,a)' ) &
2148 'm1 m2 s21 d21'
2149 if ( ipolzn.lt.0 ) write ( *, '(''+'',33x,a)' ) &
2150 'r1 r2 r3 r4'
2151 !
2152 fnorm = 4. / ( xx**2 * qsca )
2153 do 20 m = 0, nmom
2154 write ( *, '(a,i4)' ) ' moment no.', m
2155 do 20 j = 1, 4
2156 if( calcmo(j) ) then
2157 write( fmt, 98 ) 24 + (j-1)*13
2158 write ( *,fmt ) fnorm * pmom(m,j)
2159 end if
2160 20 continue
2161 end if
2162 !
2163 end if
2164 !
2165 return
2166 !
2167 98 format( '( ''+'', t', i2, ', 1p,e13.4 )' )
2168 end subroutine miprnt
2169 !**************************************************************************
2170 subroutine small1 ( xx, numang, xmu, qext, qsca, gqsc, sforw, &
2171 sback, s1, s2, tforw, tback, a, b )
2172 !
2173 ! small-particle limit of mie quantities in totally reflecting
2174 ! limit ( mie series truncated after 2 terms )
2175 !
2176 ! a,b first two mie coefficients, with numerator and
2177 ! denominator expanded in powers of -xx- ( a factor
2178 ! of xx**3 is missing but is restored before return
2179 ! to calling program ) ( ref. 2, p. 1508 )
2180 !
2181 implicit none
2182 integer numang,j
2183 real*8 gqsc, qext, qsca, xx, xmu(*)
2184 real*8 twothr,fivthr,fivnin,sq,rtmp
2185 complex*16 a( 2 ), b( 2 ), sforw, sback, s1(*), s2(*), &
2186 tforw(*), tback(*)
2187 !
2188 parameter ( twothr = 2./3., fivthr = 5./3., fivnin = 5./9. )
2189 complex*16 ctmp
2190 sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
2191 !
2192 !
2193 a( 1 ) = dcmplx ( 0.d0, twothr * ( 1. - 0.2 * xx**2 ) ) &
2194 / dcmplx ( 1.d0 - 0.5 * xx**2, twothr * xx**3 )
2195 !
2196 b( 1 ) = dcmplx ( 0.d0, - ( 1. - 0.1 * xx**2 ) / 3. ) &
2197 / dcmplx ( 1.d0 + 0.5 * xx**2, - xx**3 / 3. )
2198 !
2199 a( 2 ) = dcmplx ( 0.d0, xx**2 / 30. )
2200 b( 2 ) = dcmplx ( 0.d0, - xx**2 / 45. )
2201 !
2202 qsca = 6. * xx**4 * ( sq( a(1) ) + sq( b(1) ) &
2203 + fivthr * ( sq( a(2) ) + sq( b(2) ) ) )
2204 qext = qsca
2205 gqsc = 6. * xx**4 * dble( a(1) * dconjg( a(2) + b(1) ) &
2206 + ( b(1) + fivnin * a(2) ) * dconjg( b(2) ) )
2207 !
2208 rtmp = 1.5 * xx**3
2209 sforw = rtmp * ( a(1) + b(1) + fivthr * ( a(2) + b(2) ) )
2210 sback = rtmp * ( a(1) - b(1) - fivthr * ( a(2) - b(2) ) )
2211 tforw( 1 ) = rtmp * ( b(1) + fivthr * ( 2.*b(2) - a(2) ) )
2212 tforw( 2 ) = rtmp * ( a(1) + fivthr * ( 2.*a(2) - b(2) ) )
2213 tback( 1 ) = rtmp * ( b(1) - fivthr * ( 2.*b(2) + a(2) ) )
2214 tback( 2 ) = rtmp * ( a(1) - fivthr * ( 2.*a(2) + b(2) ) )
2215 !
2216 do 10 j = 1, numang
2217 s1( j ) = rtmp * ( a(1) + b(1) * xmu(j) + fivthr * &
2218 ( a(2) * xmu(j) + b(2) * ( 2.*xmu(j)**2 - 1. )) )
2219 s2( j ) = rtmp * ( b(1) + a(1) * xmu(j) + fivthr * &
2220 ( b(2) * xmu(j) + a(2) * ( 2.*xmu(j)**2 - 1. )) )
2221 10 continue
2222 ! ** recover actual mie coefficients
2223 a( 1 ) = xx**3 * a( 1 )
2224 a( 2 ) = xx**3 * a( 2 )
2225 b( 1 ) = xx**3 * b( 1 )
2226 b( 2 ) = xx**3 * b( 2 )
2227 !
2228 return
2229 end subroutine small1
2230 !*************************************************************************
2231 subroutine small2 ( xx, cior, calcqe, numang, xmu, qext, qsca, &
2232 gqsc, sforw, sback, s1, s2, tforw, tback, &
2233 a, b )
2234 !
2235 ! small-particle limit of mie quantities for general refractive
2236 ! index ( mie series truncated after 2 terms )
2237 !
2238 ! a,b first two mie coefficients, with numerator and
2239 ! denominator expanded in powers of -xx- ( a factor
2240 ! of xx**3 is missing but is restored before return
2241 ! to calling program ) ( ref. 2, p. 1508 )
2242 !
2243 ! ciorsq square of refractive index
2244 !
2245 implicit none
2246 logical calcqe
2247 integer numang,j
2248 real*8 gqsc, qext, qsca, xx, xmu(*)
2249 real*8 twothr,fivthr,sq,rtmp
2250 complex*16 a( 2 ), b( 2 ), cior, sforw, sback, s1(*), s2(*), &
2251 tforw(*), tback(*)
2252 !
2253 parameter ( twothr = 2./3., fivthr = 5./3. )
2254 complex*16 ctmp, ciorsq
2255 sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
2256 !
2257 !
2258 ciorsq = cior**2
2259 ctmp = dcmplx( 0.d0, twothr ) * ( ciorsq - 1.0 )
2260 a(1) = ctmp * ( 1.0 - 0.1 * xx**2 + (ciorsq/350. + 1./280.)*xx**4) &
2261 / ( ciorsq + 2.0 + ( 1.0 - 0.7 * ciorsq ) * xx**2 &
2262 - ( ciorsq**2/175. - 0.275 * ciorsq + 0.25 ) * xx**4 &
2263 + xx**3 * ctmp * ( 1.0 - 0.1 * xx**2 ) )
2264 !
2265 b(1) = (xx**2/30.) * ctmp * ( 1.0 + (ciorsq/35. - 1./14.) *xx**2 ) &
2266 / ( 1.0 - ( ciorsq/15. - 1./6. ) * xx**2 )
2267 !
2268 a(2) = ( 0.1 * xx**2 ) * ctmp * ( 1.0 - xx**2 / 14. ) &
2269 / ( 2. * ciorsq + 3. - ( ciorsq/7. - 0.5 ) * xx**2 )
2270 !
2271 qsca = 6. * xx**4 * ( sq(a(1)) + sq(b(1)) + fivthr * sq(a(2)) )
2272 gqsc = 6. * xx**4 * dble( a(1) * dconjg( a(2) + b(1) ) )
2273 qext = qsca
2274 if ( calcqe ) qext = 6. * xx * dble( a(1) + b(1) + fivthr * a(2) )
2275 !
2276 rtmp = 1.5 * xx**3
2277 sforw = rtmp * ( a(1) + b(1) + fivthr * a(2) )
2278 sback = rtmp * ( a(1) - b(1) - fivthr * a(2) )
2279 tforw( 1 ) = rtmp * ( b(1) - fivthr * a(2) )
2280 tforw( 2 ) = rtmp * ( a(1) + 2. * fivthr * a(2) )
2281 tback( 1 ) = tforw( 1 )
2282 tback( 2 ) = rtmp * ( a(1) - 2. * fivthr * a(2) )
2283 !
2284 do 10 j = 1, numang
2285 s1( j ) = rtmp * ( a(1) + ( b(1) + fivthr * a(2) ) * xmu(j) )
2286 s2( j ) = rtmp * ( b(1) + a(1) * xmu(j) + fivthr * a(2) &
2287 * ( 2. * xmu(j)**2 - 1. ) )
2288 10 continue
2289 ! ** recover actual mie coefficients
2290 a( 1 ) = xx**3 * a( 1 )
2291 a( 2 ) = xx**3 * a( 2 )
2292 b( 1 ) = xx**3 * b( 1 )
2293 b( 2 ) = ( 0., 0. )
2294 !
2295 return
2296 end subroutine small2
2297 !***********************************************************************
2298 subroutine testmi ( qext, qsca, gqsc, sforw, sback, s1, s2, &
2299 tforw, tback, pmom, momdim, ok )
2300 !
2301 ! compare mie code test case results with correct answers
2302 ! and return ok=false if even one result is inaccurate.
2303 !
2304 ! the test case is : mie size parameter = 10
2305 ! refractive index = 1.5 - 0.1 i
2306 ! scattering angle = 140 degrees
2307 ! 1 sekera moment
2308 !
2309 ! results for this case may be found among the test cases
2310 ! at the end of reference (1).
2311 !
2312 ! *** note *** when running on some computers, esp. in single
2313 ! precision, the 'accur' criterion below may have to be relaxed.
2314 ! however, if 'accur' must be set larger than 10**-3 for some
2315 ! size parameters, your computer is probably not accurate
2316 ! enough to do mie computations for those size parameters.
2317 !
2318 implicit none
2319 integer momdim,m,n
2320 real*8 qext, qsca, gqsc, pmom( 0:momdim, * )
2321 complex*16 sforw, sback, s1(*), s2(*), tforw(*), tback(*)
2322 logical ok, wrong
2323 !
2324 real*8 accur, testqe, testqs, testgq, testpm( 0:1 )
2325 complex*16 testsf, testsb,tests1,tests2,testtf(2), testtb(2)
2326 data testqe / 2.459791 /, testqs / 1.235144 /, &
2327 testgq / 1.139235 /, testsf / ( 61.49476, -3.177994 ) /, &
2328 testsb / ( 1.493434, 0.2963657 ) /, &
2329 tests1 / ( -0.1548380, -1.128972) /, &
2330 tests2 / ( 0.05669755, 0.5425681) /, &
2331 testtf / ( 12.95238, -136.6436 ), ( 48.54238, 133.4656 ) /, &
2332 testtb / ( 41.88414, -15.57833 ), ( 43.37758, -15.28196 )/, &
2333 testpm / 227.1975, 183.6898 /
2334 real*8 calc,exact
2335 ! data accur / 1.e-5 /
2336 data accur / 1.e-4 /
2337 wrong( calc, exact ) = dabs( (calc - exact) / exact ) .gt. accur
2338 !
2339 !
2340 ok = .true.
2341 if ( wrong( qext,testqe ) ) &
2342 call tstbad( 'qext', abs((qext - testqe) / testqe), ok )
2343 if ( wrong( qsca,testqs ) ) &
2344 call tstbad( 'qsca', abs((qsca - testqs) / testqs), ok )
2345 if ( wrong( gqsc,testgq ) ) &
2346 call tstbad( 'gqsc', abs((gqsc - testgq) / testgq), ok )
2347 !
2348 if ( wrong( dble(sforw), dble(testsf) ) .or. &
2349 wrong( dimag(sforw), dimag(testsf) ) ) &
2350 call tstbad( 'sforw', cdabs((sforw - testsf) / testsf), ok )
2351 !
2352 if ( wrong( dble(sback), dble(testsb) ) .or. &
2353 wrong( dimag(sback), dimag(testsb) ) ) &
2354 call tstbad( 'sback', cdabs((sback - testsb) / testsb), ok )
2355 !
2356 if ( wrong( dble(s1(1)), dble(tests1) ) .or. &
2357 wrong( dimag(s1(1)), dimag(tests1) ) ) &
2358 call tstbad( 's1', cdabs((s1(1) - tests1) / tests1), ok )
2359 !
2360 if ( wrong( dble(s2(1)), dble(tests2) ) .or. &
2361 wrong( dimag(s2(1)), dimag(tests2) ) ) &
2362 call tstbad( 's2', cdabs((s2(1) - tests2) / tests2), ok )
2363 !
2364 do 20 n = 1, 2
2365 if ( wrong( dble(tforw(n)), dble(testtf(n)) ) .or. &
2366 wrong( dimag(tforw(n)), dimag(testtf(n)) ) ) &
2367 call tstbad( 'tforw', cdabs( (tforw(n) - testtf(n)) / &
2368 testtf(n) ), ok )
2369 if ( wrong( dble(tback(n)), dble(testtb(n)) ) .or. &
2370 wrong( dimag(tback(n)), dimag(testtb(n)) ) ) &
2371 call tstbad( 'tback', cdabs( (tback(n) - testtb(n)) / &
2372 testtb(n) ), ok )
2373 20 continue
2374 !
2375 do 30 m = 0, 1
2376 if ( wrong( pmom(m,1), testpm(m) ) ) &
2377 call tstbad( 'pmom', dabs( (pmom(m,1)-testpm(m)) / &
2378 testpm(m) ), ok )
2379 30 continue
2380 !
2381 return
2382 !
2383 end subroutine testmi
2384 !**************************************************************************
2385 subroutine errmsg( messag, fatal )
2386 !
2387 ! print out a warning or error message; abort if error
2388 !
2389 USE module_peg_util, only: peg_message, peg_error_fatal
2390
2391 implicit none
2392 logical fatal, once
2393 character*80 msg
2394 character*(*) messag
2395 integer maxmsg, nummsg
2396 save maxmsg, nummsg, once
2397 data nummsg / 0 /, maxmsg / 100 /, once / .false. /
2398 !
2399 !
2400 if ( fatal ) then
2401 ! write ( *, '(2a)' ) ' ******* error >>>>>> ', messag
2402 ! stop
2403 write( msg, '(a)' ) &
2404 'FASTJ mie fatal error ' // &
2405 messag
2406 call peg_message( lunerr, msg )
2407 call peg_error_fatal( lunerr, msg )
2408 end if
2409 !
2410 nummsg = nummsg + 1
2411 if ( nummsg.gt.maxmsg ) then
2412 ! if ( .not.once ) write ( *,99 )
2413 if ( .not.once )then
2414 write( msg, '(a)' ) &
2415 'FASTJ mie: too many warning messages -- no longer printing '
2416 call peg_message( lunerr, msg )
2417 end if
2418 once = .true.
2419 else
2420 msg = 'FASTJ mie warning ' // messag
2421 call peg_message( lunerr, msg )
2422 ! write ( *, '(2a)' ) ' ******* warning >>>>>> ', messag
2423 endif
2424 !
2425 return
2426 !
2427 ! 99 format( ///,' >>>>>> too many warning messages -- ', &
2428 ! 'they will no longer be printed <<<<<<<', /// )
2429 end subroutine errmsg
2430 !********************************************************************
2431 subroutine wrtbad ( varnam, erflag )
2432 !
2433 ! write names of erroneous variables
2434 !
2435 ! input : varnam = name of erroneous variable to be written
2436 ! ( character, any length )
2437 !
2438 ! output : erflag = logical flag, set true by this routine
2439 ! ----------------------------------------------------------------------
2440 USE module_peg_util, only: peg_message
2441
2442 implicit none
2443 character*(*) varnam
2444 logical erflag
2445 integer maxmsg, nummsg
2446 character*80 msg
2447 save nummsg, maxmsg
2448 data nummsg / 0 /, maxmsg / 50 /
2449 !
2450 !
2451 nummsg = nummsg + 1
2452 ! write ( *, '(3a)' ) ' **** input variable ', varnam, &
2453 ! ' in error ****'
2454 msg = 'FASTJ mie input variable in error ' // varnam
2455 call peg_message( lunerr, msg )
2456 erflag = .true.
2457 if ( nummsg.eq.maxmsg ) &
2458 call errmsg ( 'too many input variable errors. aborting...$', .true. )
2459 return
2460 !
2461 end subroutine wrtbad
2462 !******************************************************************
2463 subroutine tstbad( varnam, relerr, ok )
2464 !
2465 ! write name (-varnam-) of variable failing self-test and its
2466 ! percent error from the correct value. return ok = false.
2467 !
2468 implicit none
2469 character*(*) varnam
2470 logical ok
2471 real*8 relerr
2472 !
2473 !
2474 ok = .false.
2475 write( *, '(/,3a,1p,e11.2,a)' ) &
2476 ' output variable ', varnam,' differed by', 100.*relerr, &
2477 ' per cent from correct value. self-test failed.'
2478 return
2479 !
2480 end subroutine tstbad
2481 !-----------------------------------------------------------------------
2482 end module module_fastj_mie