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