!************************************************************************
! This computer software was prepared by Battelle Memorial Institute,
! hereinafter the Contractor, under Contract No. DE-AC05-76RL0 1830 with
! the Department of Energy (DOE). NEITHER THE GOVERNMENT NOR THE
! CONTRACTOR MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY
! LIABILITY FOR THE USE OF THIS SOFTWARE.
!
! Module to Compute Aerosol Optical Properties
! * Primary investigator: James C. Barnard
! * Co-investigators: Rahul A. Zaveri, Richard C. Easter, William I.
! Gustafson Jr.
! Last update: February 2009
!
! Contact:
! Jerome D. Fast, PhD
! Staff Scientist
! Pacific Northwest National Laboratory
! P.O. Box 999, MSIN K9-30
! Richland, WA, 99352
! Phone: (509) 372-6116
! Email: Jerome.Fast@pnl.gov
!
! Please report any bugs or problems to Jerome Fast, the WRF-chem
! implmentation team leader for PNNL.
!
! Terms of Use:
! 1) Users are requested to consult the primary author prior to
! modifying this module or incorporating it or its submodules in
! another code. This is meant to ensure that the any linkages and/or
! assumptions will not adversely affect the operation of this module.
! 2) The source code in this module is intended for research and
! educational purposes. Users are requested to contact the primary
! author regarding the use of the code for any commercial application.
! 3) Users preparing publications resulting from the usage of this code
! are requested to cite one or more of the references below
! (depending on the application) for proper acknowledgement.
!
! Note that the aerosol optical properites are currently tied to the use
! of Fast-J and MOSAIC. Future modifications will make the calculations
! more generic so the code is not tied to the photolysis scheme and the
! code can be used for both modal and sectional treatments.
!
! References:
! * Fast, J.D., W.I. Gustafson Jr., R.C. Easter, R.A. Zaveri, J.C.
! Barnard, E.G. Chapman, G.A. Grell, and S.E. Peckham (2005), Evolution
! of ozone, particulates, and aerosol direct radiative forcing in the
! vicinity of Houston using a fully-coupled meteorology-chemistry-
! aerosol model, J. Geophys. Res., 111, D21305,
! doi:10.1029/2005JD006721.
! * Barnard, J.C., J.D. Fast, G. Paredes-Miranda, W.P. Arnott, and
! A. Laskin (2010), Technical note: evaluation of the WRF-Chem
! "aerosol chemical to aerosol optical properties" module using data
! from the MILAGRO campaign, Atmos. Chem. Phys., 10, 7325-7340,
! doi:10.5194/acp-10-7325-2010.
!
! Contact Jerome Fast for updates on the status of manuscripts under
! review.
!
! Additional information:
! * www.pnl.gov/atmospheric/research/wrf-chem
!
! Support:
! Funding for adapting Fast-J was provided by the U.S. Department of
! Energy under the auspices of Atmospheric Science Program of the Office
! of Biological and Environmental Research the PNNL Laboratory Research
! and Directed Research and Development program.
!************************************************************************
module module_fastj_mie
! jcb 2007-feb-01 added capability to calculate aerosol backscatter
! units (km)^-1
! need to divide by 4*pie to get units of (km*ster)-1
! jcb 2007-feb-01 now calculate true extinctionm units km^-1
! rce 2005-apr-22 - want lunerr = -1 for wrf-chem 3d;
! also, define lunerr here, and it's available everywhere
integer, parameter :: lunerr = -1
contains
!***********************************************************************
! <1.> subr mieaer
! Purpose: calculate aerosol optical depth, single scattering albedo,
! asymmetry factor, extinction, Legendre coefficients, and average aerosol
! size. parameterizes aerosol coefficients using chebychev polynomials
! requires double precision on 32-bit machines
! uses Wiscombe's (1979) mie scattering code
! INPUT
! id -- grid id number
! iclm, jclm -- i,j of grid column being processed
! nbin_a -- number of bins
! number_bin(nbin_a,kmaxd) -- number density in layer, #/cm^3
! radius_wet(nbin_a,kmaxd) -- wet radius, cm
! refindx(nbin_a,kmaxd) --volume averaged complex index of refraction
! dz -- depth of individual cells in column, m
! isecfrm0 - time from start of run, sec
! lpar -- number of grid cells in vertical (via module_fastj_cmnh)
! kmaxd -- predefined maximum allowed levels from module_data_mosaic_other
! passed here via module_fastj_cmnh
! OUTPUT: saved in module_fastj_cmnmie
! real tauaer ! aerosol optical depth
! waer ! aerosol single scattering albedo
! gaer ! aerosol asymmetery factor
! extaer ! aerosol extinction
! l2,l3,l4,l5,l6,l7 ! Legendre coefficients, numbered 0,1,2,......
! sizeaer ! average wet radius
! bscoef ! aerosol backscatter coefficient with units km-1 * steradian -1 JCB 2007/02/01
!----------------------------------------------------------------------
subroutine mieaer( &
id, iclm, jclm, nbin_a, &
number_bin, radius_wet, refindx, &
dz, isecfrm0, lpar, &
sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7,bscoef) ! added bscoef JCB 2007/02/01
USE module_data_mosaic_other, only : kmaxd
USE module_data_mosaic_therm, ONLY: nbin_a_maxd
USE module_peg_util, only: peg_error_fatal, peg_message
IMPLICIT NONE
! subr arguments
!jdf
integer,parameter :: nspint = 4 ! Num of spectral intervals across
! solar spectrum for FAST-J
integer, intent(in) :: lpar
real, dimension (nspint, kmaxd+1),intent(out) :: sizeaer,extaer,waer,gaer,tauaer
real, dimension (nspint, kmaxd+1),intent(out) :: l2,l3,l4,l5,l6,l7
real, dimension (nspint, kmaxd+1),intent(out) :: bscoef !JCB 2007/02/01
real, dimension (nspint),save :: wavmid !cm
data wavmid &
/ 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
!jdf
integer, intent(in) :: id, iclm, jclm, nbin_a, isecfrm0
real, intent(in), dimension(nbin_a_maxd, kmaxd) :: number_bin
real, intent(inout), dimension(nbin_a_maxd, kmaxd) :: radius_wet
complex, intent(in) :: refindx(nbin_a_maxd, kmaxd)
real, intent(in) :: dz(lpar)
!local variables
real weighte, weights
! various bookeeping variables
integer ltype ! total number of indicies of refraction
parameter (ltype = 1) ! bracket refractive indices based on information from Rahul, 2002/11/07
real x
real thesum ! for normalizing things
real sizem ! size in microns
integer kcallmieaer
!
integer m, j, nc, klevel
real pext ! parameterized specific extinction (cm2/g)
real pasm ! parameterized asymmetry factor
real ppmom2 ! 2 Lengendre expansion coefficient (numbered 0,1,2,...)
real ppmom3 ! 3 ...
real ppmom4 ! 4 ...
real ppmom5 ! 5 ...
real ppmom6 ! 6 ...
real ppmom7 ! 7 ...
real sback2 ! JCB 2007/02/01 sback*conjg(sback)
integer ns ! Spectral loop index
integer i ! Longitude loop index
integer k ! Level loop index
real pscat !scattering cross section
integer prefr,prefi,nrefr,nrefi,nr,ni
save nrefr,nrefi
parameter (prefr=7,prefi=7)
complex*16 sforw,sback,tforw(2),tback(2)
integer numang,nmom,ipolzn,momdim
real*8 pmom(0:7,1)
logical perfct,anyang,prnt(2)
logical first
integer, parameter :: nsiz=200,nlog=30,ncoef=50
!
real p2(nsiz),p3(nsiz),p4(nsiz),p5(nsiz)
real p6(nsiz),p7(nsiz)
!
real*8 xmu(1)
data xmu/1./,anyang/.false./
data numang/0/
complex*16 s1(1),s2(1)
data first/.true./
save first
real*8 mimcut
data perfct/.false./,mimcut/0.0/
data nmom/7/,ipolzn/0/,momdim/7/
data prnt/.false.,.false./
! coefficients for parameterizing aerosol radiative properties
! in terms of refractive index and wet radius
real extp(ncoef,prefr,prefi,nspint) ! specific extinction
real albp(ncoef,prefr,prefi,nspint) ! single scat alb
real asmp(ncoef,prefr,prefi,nspint) ! asymmetry factor
real ascat(ncoef,prefr,prefi,nspint) ! scattering efficiency, JCB 2004/02/09
real pmom2(ncoef,prefr,prefi,nspint) ! phase function expansion, #2
real pmom3(ncoef,prefr,prefi,nspint) ! phase function expansion, #3
real pmom4(ncoef,prefr,prefi,nspint) ! phase function expansion, #4
real pmom5(ncoef,prefr,prefi,nspint) ! phase function expansion, #5
real pmom6(ncoef,prefr,prefi,nspint) ! phase function expansion, #6
real pmom7(ncoef,prefr,prefi,nspint) ! phase function expansion, #7
real sback2p(ncoef,prefr,prefi,nspint) ! JCB 2007/02/01 - backscatter
save :: extp,albp,asmp,ascat,pmom2,pmom3,pmom4,pmom5,pmom6,pmom7,sback2p ! JCB 2007/02/01 - added sback2p
!--------------
real cext(ncoef),casm(ncoef),cpmom2(ncoef)
real cscat(ncoef) ! JCB 2004/02/09
real cpmom3(ncoef),cpmom4(ncoef),cpmom5(ncoef)
real cpmom6(ncoef),cpmom7(ncoef)
real cpsback2p(ncoef) ! JCB 2007/02/09 - backscatter
integer itab,jtab
real ttab,utab
! nsiz = number of wet particle sizes
! crefin = complex refractive index
integer n
real*8 thesize ! 2 pi radpart / waveleng = size parameter
real*8 qext(nsiz) ! array of extinction efficiencies
real*8 qsca(nsiz) ! array of scattering efficiencies
real*8 gqsc(nsiz) ! array of asymmetry factor * scattering efficiency
real asymm(nsiz) ! array of asymmetry factor
real scat(nsiz) ! JCB 2004/02/09
real sb2(nsiz) ! JCB 2007/02/01 - 4*abs(sback)^2/(size parameter)^2 backscattering efficiency
! specabs = absorption coeff / unit dry mass
! specscat = scattering coeff / unit dry mass
complex*16 crefin,crefd,crefw
save crefw
real, save :: rmin,rmax ! min, max aerosol size bin
data rmin /0.005e-4/ ! rmin in cm. 5e-3 microns min allowable size
data rmax /50.e-4 / ! rmax in cm. 50 microns, big particle, max allowable size
real bma,bpa
real, save :: xrmin,xrmax,xr
real rs(nsiz) ! surface mode radius (cm)
real xrad ! normalized aerosol radius
real ch(ncoef) ! chebychev polynomial
real, save :: rhoh2o ! density of liquid water (g/cm3)
data rhoh2o/1./
real refr ! real part of refractive index
real refi ! imaginary part of refractive index
real qextr4(nsiz) ! extinction, real*4
real refrmin ! minimum of real part of refractive index
real refrmax ! maximum of real part of refractive index
real refimin ! minimum of imag part of refractive index
real refimax ! maximum of imag part of refractive index
real drefr ! increment in real part of refractive index
real drefi ! increment in imag part of refractive index
real, save :: refrtab(prefr) ! table of real refractive indices for aerosols
real, save :: refitab(prefi) ! table of imag refractive indices for aerosols
complex specrefndx(ltype) ! refractivr indices
real pie,third
integer irams, jrams
! diagnostic declarations
integer kcallmieaer2
integer ibin
character*150 msg
#if (defined(CHEM_DBG_I) && defined(CHEM_DBG_J) && defined(CHEM_DBG_K))
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ec diagnostics
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ec run_out.25 has aerosol physical parameter info for bins 1-8
!ec and vertical cells 1 to kmaxd.
! ilaporte = 33
! jlaporte = 34
if (iclm .eq. CHEM_DBG_I) then
if (jclm .eq. CHEM_DBG_J) then
! initial entry
if (kcallmieaer2 .eq. 0) then
write(*,9099)iclm, jclm
9099 format('for cell i = ', i3, 2x, 'j = ', i3)
write(*,9100)
9100 format( &
'isecfrm0', 3x, 'i', 3x, 'j', 3x,'k', 3x, &
'ibin', 3x, &
'refindx(ibin,k)', 3x, &
'radius_wet(ibin,k)', 3x, &
'number_bin(ibin,k)' &
)
end if
!ec output for run_out.25
do k = 1, lpar
do ibin = 1, nbin_a
write(*, 9120) &
isecfrm0,iclm, jclm, k, ibin, &
refindx(ibin,k), &
radius_wet(ibin,k), &
number_bin(ibin,k)
9120 format( i7,3(2x,i4),2x,i4, 4x, 4(e14.6,2x))
end do
end do
kcallmieaer2 = kcallmieaer2 + 1
end if
end if
!ec end print of aerosol physical parameter diagnostics
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
#endif
!
! assign fast-J wavelength, these values are in cm
! wavmid(1)=0.30e-4
! wavmid(2)=0.40e-4
! wavmid(3)=0.60e-4
! wavmid(4)=0.999e-04
!
pie=4.*atan(1.)
third=1./3.
if(first)then
first=.false.
! parameterize aerosol radiative properties in terms of
! relative humidity, surface mode wet radius, aerosol species,
! and wavelength
! first find min,max of real and imaginary parts of refractive index
! real and imaginary parts of water refractive index
crefw=cmplx(1.33,0.0)
refrmin=real(crefw)
refrmax=real(crefw)
! change Rahul's imaginary part of the refractive index from positive to negative
refimin=-imag(crefw)
refimax=-imag(crefw)
! aerosol mode loop
specrefndx(1)=cmplx(1.82,-0.74) ! max values from Rahul, 7 Nov 2002
!
do i=1,ltype ! loop over all possible refractive indices
! real and imaginary parts of aerosol refractive index
refrmin=amin1(refrmin,real(specrefndx(ltype)))
refrmax=amax1(refrmax,real(specrefndx(ltype)))
refimin=amin1(refimin,aimag(specrefndx(ltype)))
refimax=amax1(refimax,aimag(specrefndx(ltype)))
enddo
drefr=(refrmax-refrmin)
if(drefr.gt.1.e-4)then
nrefr=prefr
drefr=drefr/(nrefr-1)
else
nrefr=1
endif
drefi=(refimax-refimin)
if(drefi.gt.1.e-4)then
nrefi=prefi
drefi=drefi/(nrefi-1)
else
nrefi=1
endif
!
bma=0.5*alog(rmax/rmin) ! JCB
bpa=0.5*alog(rmax*rmin) ! JCB
xrmin=alog(rmin)
xrmax=alog(rmax)
! wavelength loop
do 200 ns=1,nspint
! calibrate parameterization with range of refractive indices
do 120 nr=1,nrefr
do 120 ni=1,nrefi
refrtab(nr)=refrmin+(nr-1)*drefr
refitab(ni)=refimin+(ni-1)*drefi
crefd=cmplx(refrtab(nr),refitab(ni))
! mie calculations of optical efficiencies
do n=1,nsiz
xr=cos(pie*(float(n)-0.5)/float(nsiz))
rs(n)=exp(xr*bma+bpa)
! size parameter and weighted refractive index
thesize=2.*pie*rs(n)/wavmid(ns)
thesize=min(thesize,10000.d0)
call miev0(thesize,crefd,perfct,mimcut,anyang, &
numang,xmu,nmom,ipolzn,momdim,prnt, &
qext(n),qsca(n),gqsc(n),pmom,sforw,sback,s1, &
s2,tforw,tback )
qextr4(n)=qext(n)
! qabs(n)=qext(n)-qsca(n) ! not necessary anymore JCB 2004/02/09
scat(n)=qsca(n) ! JCB 2004/02/09
asymm(n)=gqsc(n)/qsca(n) ! assume always greater than zero
! coefficients of phase function expansion; note modification by JCB of miev0 coefficients
p2(n)=pmom(2,1)/pmom(0,1)*5.0
p3(n)=pmom(3,1)/pmom(0,1)*7.0
p4(n)=pmom(4,1)/pmom(0,1)*9.0
p5(n)=pmom(5,1)/pmom(0,1)*11.0
p6(n)=pmom(6,1)/pmom(0,1)*13.0
p7(n)=pmom(7,1)/pmom(0,1)*15.0
! backscattering efficiency, Bohren and Huffman, page 122
! as stated by Bohren and Huffman, this is 4*pie times what is should be
! may need to be smoothed - a very rough function - for the time being we won't apply smoothing
! and let the integration over the size distribution be the smoothing
sb2(n)=4.0*sback*dconjg(sback)/(thesize*thesize) ! JCB 2007/02/01
enddo
100 continue
!
call fitcurv(rs,qextr4,extp(1,nr,ni,ns),ncoef,nsiz)
call fitcurv(rs,scat,ascat(1,nr,ni,ns),ncoef,nsiz) ! JCB 2004/02/07 - scattering efficiency
call fitcurv(rs,asymm,asmp(1,nr,ni,ns),ncoef,nsiz)
call fitcurv_nolog(rs,p2,pmom2(1,nr,ni,ns),ncoef,nsiz)
call fitcurv_nolog(rs,p3,pmom3(1,nr,ni,ns),ncoef,nsiz)
call fitcurv_nolog(rs,p4,pmom4(1,nr,ni,ns),ncoef,nsiz)
call fitcurv_nolog(rs,p5,pmom5(1,nr,ni,ns),ncoef,nsiz)
call fitcurv_nolog(rs,p6,pmom6(1,nr,ni,ns),ncoef,nsiz)
call fitcurv_nolog(rs,p7,pmom7(1,nr,ni,ns),ncoef,nsiz)
call fitcurv(rs,sb2,sback2p(1,nr,ni,ns),ncoef,nsiz) ! JCB 2007/02/01 - backscattering efficiency
120 continue
200 continue
endif
! begin level loop
do 2000 klevel=1,lpar
! sum densities for normalization
thesum=0.0
do m=1,nbin_a
thesum=thesum+number_bin(m,klevel)
enddo
! Begin spectral loop
do 1000 ns=1,nspint
! aerosol optical properties
tauaer(ns,klevel)=0.
waer(ns,klevel)=0.
gaer(ns,klevel)=0.
sizeaer(ns,klevel)=0.0
extaer(ns,klevel)=0.0
l2(ns,klevel)=0.0
l3(ns,klevel)=0.0
l4(ns,klevel)=0.0
l5(ns,klevel)=0.0
l6(ns,klevel)=0.0
l7(ns,klevel)=0.0
bscoef(ns,klevel)=0.0 ! JCB 2007/02/01 - backscattering coefficient
if(thesum.le.1e-21)goto 1000 ! set everything = 0 if no aerosol !wig changed 0.0 to 1e-21, 31-Oct-2005
! loop over the bins
do m=1,nbin_a ! nbin_a is number of bins
! check to see if there's any aerosol
if(number_bin(m,klevel).le.1e-21)goto 70 ! no aerosol !wig changed 0.0 to 1e-21, 31-Oct-2005
! here's the size
sizem=radius_wet(m,klevel) ! radius in cm
! check limits of particle size
! rce 2004-dec-07 - use klevel in write statements
if(radius_wet(m,klevel).le.rmin)then
radius_wet(m,klevel)=rmin
write( msg, '(a, 5i4,1x, e11.4)' ) &
'FASTJ mie: radius_wet set to rmin,' // &
'id,i,j,k,m,rm(m,k)', id, iclm, jclm, klevel, m, radius_wet(m,klevel)
call peg_message( lunerr, msg )
! write(6,'('' particle size too small '')')
endif
!
if(radius_wet(m,klevel).gt.rmax)then
write( msg, '(a, 5i4,1x, e11.4)' ) &
'FASTJ mie: radius_wet set to rmax,' // &
'id,i,j,k,m,rm(m,k)', &
id, iclm, jclm, klevel, m, radius_wet(m,klevel)
call peg_message( lunerr, msg )
radius_wet(m,klevel)=rmax
! write(6,'('' particle size too large '')')
endif
!
x=alog(radius_wet(m,klevel)) ! radius in cm
!
crefin=refindx(m,klevel)
refr=real(crefin)
! change Rahul's imaginary part of the index of refraction from positive to negative
refi=-imag(crefin)
!
xrad=x
thesize=2.0*pie*exp(x)/wavmid(ns)
! normalize size parameter
xrad=(2*xrad-xrmax-xrmin)/(xrmax-xrmin)
! retain this diagnostic code
if(abs(refr).gt.10.0.or.abs(refr).le.0.001)then
write ( msg, '(a,1x, e14.5)' ) &
'FASTJ mie /refr/ outside range 1e-3 - 10 ' // &
'refr= ', refr
! print *,'refr=',refr
call peg_error_fatal( lunerr, msg )
endif
if(abs(refi).gt.10.)then
write ( msg, '(a,1x, e14.5)' ) &
'FASTJ mie /refi/ >10 ' // &
'refi', refi
! print *,'refi=',refi
call peg_error_fatal( lunerr, msg )
endif
! interpolate coefficients linear in refractive index
! first call calcs itab,jtab,ttab,utab
itab=0
call binterp(extp(1,1,1,ns),ncoef,nrefr,nrefi, &
refr,refi,refrtab,refitab,itab,jtab, &
ttab,utab,cext)
! JCB 2004/02/09 -- new code for scattering cross section
call binterp(ascat(1,1,1,ns),ncoef,nrefr,nrefi, &
refr,refi,refrtab,refitab,itab,jtab, &
ttab,utab,cscat)
call binterp(asmp(1,1,1,ns),ncoef,nrefr,nrefi, &
refr,refi,refrtab,refitab,itab,jtab, &
ttab,utab,casm)
call binterp(pmom2(1,1,1,ns),ncoef,nrefr,nrefi, &
refr,refi,refrtab,refitab,itab,jtab, &
ttab,utab,cpmom2)
call binterp(pmom3(1,1,1,ns),ncoef,nrefr,nrefi, &
refr,refi,refrtab,refitab,itab,jtab, &
ttab,utab,cpmom3)
call binterp(pmom4(1,1,1,ns),ncoef,nrefr,nrefi, &
refr,refi,refrtab,refitab,itab,jtab, &
ttab,utab,cpmom4)
call binterp(pmom5(1,1,1,ns),ncoef,nrefr,nrefi, &
refr,refi,refrtab,refitab,itab,jtab, &
ttab,utab,cpmom5)
call binterp(pmom6(1,1,1,ns),ncoef,nrefr,nrefi, &
refr,refi,refrtab,refitab,itab,jtab, &
ttab,utab,cpmom6)
call binterp(pmom7(1,1,1,ns),ncoef,nrefr,nrefi, &
refr,refi,refrtab,refitab,itab,jtab, &
ttab,utab,cpmom7)
call binterp(sback2p(1,1,1,ns),ncoef,nrefr,nrefi, &
refr,refi,refrtab,refitab,itab,jtab, &
ttab,utab,cpsback2p)
! chebyshev polynomials
ch(1)=1.
ch(2)=xrad
do nc=3,ncoef
ch(nc)=2.*xrad*ch(nc-1)-ch(nc-2)
enddo
! parameterized optical properties
pext=0.5*cext(1)
do nc=2,ncoef
pext=pext+ch(nc)*cext(nc)
enddo
pext=exp(pext)
! JCB 2004/02/09 -- for scattering efficiency
pscat=0.5*cscat(1)
do nc=2,ncoef
pscat=pscat+ch(nc)*cscat(nc)
enddo
pscat=exp(pscat)
!
pasm=0.5*casm(1)
do nc=2,ncoef
pasm=pasm+ch(nc)*casm(nc)
enddo
pasm=exp(pasm)
!
ppmom2=0.5*cpmom2(1)
do nc=2,ncoef
ppmom2=ppmom2+ch(nc)*cpmom2(nc)
enddo
if(ppmom2.le.0.0)ppmom2=0.0
!
ppmom3=0.5*cpmom3(1)
do nc=2,ncoef
ppmom3=ppmom3+ch(nc)*cpmom3(nc)
enddo
! ppmom3=exp(ppmom3) ! no exponentiation required
if(ppmom3.le.0.0)ppmom3=0.0
!
ppmom4=0.5*cpmom4(1)
do nc=2,ncoef
ppmom4=ppmom4+ch(nc)*cpmom4(nc)
enddo
if(ppmom4.le.0.0.or.sizem.le.0.03e-04)ppmom4=0.0
!
ppmom5=0.5*cpmom5(1)
do nc=2,ncoef
ppmom5=ppmom5+ch(nc)*cpmom5(nc)
enddo
if(ppmom5.le.0.0.or.sizem.le.0.03e-04)ppmom5=0.0
!
ppmom6=0.5*cpmom6(1)
do nc=2,ncoef
ppmom6=ppmom6+ch(nc)*cpmom6(nc)
enddo
if(ppmom6.le.0.0.or.sizem.le.0.03e-04)ppmom6=0.0
!
ppmom7=0.5*cpmom7(1)
do nc=2,ncoef
ppmom7=ppmom7+ch(nc)*cpmom7(nc)
enddo
if(ppmom7.le.0.0.or.sizem.le.0.03e-04)ppmom7=0.0
!
sback2=0.5*cpsback2p(1) ! JCB 2007/02/01 - backscattering efficiency
do nc=2,ncoef
sback2=sback2+ch(nc)*cpsback2p(nc)
enddo
sback2=exp(sback2)
if(sback2.le.0.0)sback2=0.0
!
!
! weights:
weighte=pext*pie*exp(x)**2 ! JCB, extinction cross section
weights=pscat*pie*exp(x)**2 ! JCB, scattering cross section
tauaer(ns,klevel)=tauaer(ns,klevel)+weighte*number_bin(m,klevel) ! must be multiplied by deltaZ
! extaer(ns,klevel)=extaer(ns,klevel)+pext*number_bin(m,klevel) ! not the true extinction, JCB 2007/02/01
sizeaer(ns,klevel)=sizeaer(ns,klevel)+exp(x)*10000.0* &
number_bin(m,klevel)
waer(ns,klevel)=waer(ns,klevel)+weights*number_bin(m,klevel) !JCB
gaer(ns,klevel)=gaer(ns,klevel)+pasm*weights* &
number_bin(m,klevel) !JCB
! need weighting by scattering cross section ? JCB 2004/02/09
l2(ns,klevel)=l2(ns,klevel)+weights*ppmom2*number_bin(m,klevel)
l3(ns,klevel)=l3(ns,klevel)+weights*ppmom3*number_bin(m,klevel)
l4(ns,klevel)=l4(ns,klevel)+weights*ppmom4*number_bin(m,klevel)
l5(ns,klevel)=l5(ns,klevel)+weights*ppmom5*number_bin(m,klevel)
l6(ns,klevel)=l6(ns,klevel)+weights*ppmom6*number_bin(m,klevel)
l7(ns,klevel)=l7(ns,klevel)+weights*ppmom7*number_bin(m,klevel)
! convert backscattering efficiency to backscattering coefficient, units (cm)^-1
bscoef(ns,klevel)=bscoef(ns,klevel)+pie*exp(x)**2*sback2*number_bin(m,klevel) ! JCB 2007/02/01 - backscatter coefficient,
end do ! end of nbin_a loop
! take averages - weighted by cross section - new code JCB 2004/02/09
sizeaer(ns,klevel)=sizeaer(ns,klevel)/thesum
gaer(ns,klevel)=gaer(ns,klevel)/waer(ns,klevel) ! JCB removed *3 factor 2/9/2004
! because factor is applied in subroutine opmie, file zz01fastj_mod.f
l2(ns,klevel)=l2(ns,klevel)/waer(ns,klevel)
l3(ns,klevel)=l3(ns,klevel)/waer(ns,klevel)
l4(ns,klevel)=l4(ns,klevel)/waer(ns,klevel)
l5(ns,klevel)=l5(ns,klevel)/waer(ns,klevel)
l6(ns,klevel)=l6(ns,klevel)/waer(ns,klevel)
l7(ns,klevel)=l7(ns,klevel)/waer(ns,klevel)
! backscatter coef, divide by 4*Pie to get units of (km*ster)^-1 JCB 2007/02/01
bscoef(ns,klevel)=bscoef(ns,klevel)*1.0e5 ! JCB 2007/02/01 - units are now (km)^-1
extaer(ns,klevel)=tauaer(ns,klevel)*1.0e5 ! JCB 2007/02/01 - now true extincion, units (km)^-1
! this must be last!! JCB 2007/02/01
waer(ns,klevel)=waer(ns,klevel)/tauaer(ns,klevel) ! JCB
70 continue ! bail out if no aerosol;go on to next wavelength bin
1000 continue ! end of wavelength loop
2000 continue ! end of klevel loop
!
! before returning, multiply tauaer by depth of individual cells.
! tauaer is in cm-1, dz in m; multiply dz by 100 to convert from m to cm.
do ns = 1, nspint
do klevel = 1, lpar
tauaer(ns,klevel) = tauaer(ns,klevel) * dz(klevel)* 100.
end do
end do
#if (defined(CHEM_DBG_I) && defined(CHEM_DBG_J) && defined(CHEM_DBG_K))
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ec fastj diagnostics
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ec run_out.30 has aerosol optical info for cells 1 to kmaxd.
! ilaporte = 33
! jlaporte = 34
if (iclm .eq. CHEM_DBG_I) then
if (jclm .eq. CHEM_DBG_J) then
! initial entry
if (kcallmieaer .eq. 0) then
write(*,909) CHEM_DBG_I, CHEM_DBG_J
909 format( ' for cell i = ', i3, ' j = ', i3)
write(*,910)
910 format( &
'isecfrm0', 3x, 'i', 3x, 'j', 3x,'k', 3x, &
'dzmfastj', 8x, &
'tauaer(1,k)',1x, 'tauaer(2,k)',1x,'tauaer(3,k)',3x, &
'tauaer(4,k)',5x, &
'waer(1,k)', 7x, 'waer(2,k)', 7x,'waer(3,k)', 7x, &
'waer(4,k)', 7x, &
'gaer(1,k)', 7x, 'gaer(2,k)', 7x,'gaer(3,k)', 7x, &
'gaer(4,k)', 7x, &
'extaer(1,k)',5x, 'extaer(2,k)',5x,'extaer(3,k)',5x, &
'extaer(4,k)',5x, &
'sizeaer(1,k)',4x, 'sizeaer(2,k)',4x,'sizeaer(3,k)',4x, &
'sizeaer(4,k)' )
end if
!ec output for run_out.30
do k = 1, lpar
write(*, 912) &
isecfrm0,iclm, jclm, k, &
dz(k) , &
(tauaer(n,k), n=1,4), &
(waer(n,k), n=1,4), &
(gaer(n,k), n=1,4), &
(extaer(n,k), n=1,4), &
(sizeaer(n,k), n=1,4)
912 format( i7,3(2x,i4),2x,21(e14.6,2x))
end do
kcallmieaer = kcallmieaer + 1
end if
end if
!ec end print of fastj diagnostics
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
#endif
return
end subroutine mieaer
!****************************************************************
subroutine fitcurv(rs,yin,coef,ncoef,maxm)
! fit y(x) using Chebychev polynomials
! wig 7-Sep-2004: Removed dependency on pre-determined maximum
! array size and replaced with f90 array info.
USE module_peg_util, only: peg_message
IMPLICIT NONE
! integer nmodes, nrows, maxm, ncoef
! parameter (nmodes=500,nrows=8)
integer, intent(in) :: maxm, ncoef
! real rs(nmodes),yin(nmodes),coef(ncoef)
! real x(nmodes),y(nmodes)
real, dimension(ncoef) :: coef
real, dimension(:) :: rs, yin
real x(size(rs)),y(size(yin))
integer m
real xmin, xmax
character*80 msg
!!$ if(maxm.gt.nmodes)then
!!$ write ( msg, '(a, 1x,i6)' ) &
!!$ 'FASTJ mie nmodes too small in fitcurv, ' // &
!!$ 'maxm ', maxm
!!$! write(*,*)'nmodes too small in fitcurv',maxm
!!$ call peg_error_fatal( lunerr, msg )
!!$ endif
do 100 m=1,maxm
x(m)=alog(rs(m))
y(m)=alog(yin(m))
100 continue
xmin=x(1)
xmax=x(maxm)
do 110 m=1,maxm
x(m)=(2*x(m)-xmax-xmin)/(xmax-xmin)
110 continue
call chebft(coef,ncoef,maxm,y)
return
end subroutine fitcurv
!**************************************************************
subroutine fitcurv_nolog(rs,yin,coef,ncoef,maxm)
! fit y(x) using Chebychev polynomials
! wig 7-Sep-2004: Removed dependency on pre-determined maximum
! array size and replaced with f90 array info.
USE module_peg_util, only: peg_message
IMPLICIT NONE
! integer nmodes, nrows, maxm, ncoef
! parameter (nmodes=500,nrows=8)
integer, intent(in) :: maxm, ncoef
! real rs(nmodes),yin(nmodes),coef(ncoef)
real, dimension(:) :: rs, yin
real, dimension(ncoef) :: coef(ncoef)
real x(size(rs)),y(size(yin))
integer m
real xmin, xmax
character*80 msg
!!$ if(maxm.gt.nmodes)then
!!$ write ( msg, '(a,1x, i6)' ) &
!!$ 'FASTJ mie nmodes too small in fitcurv ' // &
!!$ 'maxm ', maxm
!!$! write(*,*)'nmodes too small in fitcurv',maxm
!!$ call peg_error_fatal( lunerr, msg )
!!$ endif
do 100 m=1,maxm
x(m)=alog(rs(m))
y(m)=yin(m) ! note, no "alog" here
100 continue
xmin=x(1)
xmax=x(maxm)
do 110 m=1,maxm
x(m)=(2*x(m)-xmax-xmin)/(xmax-xmin)
110 continue
call chebft(coef,ncoef,maxm,y)
return
end subroutine fitcurv_nolog
!************************************************************************
subroutine chebft(c,ncoef,n,f)
! given a function f with values at zeroes x_k of Chebychef polynomial
! T_n(x), calculate coefficients c_j such that
! f(x)=sum(k=1,n) c_k t_(k-1)(y) - 0.5*c_1
! where y=(x-0.5*(xmax+xmin))/(0.5*(xmax-xmin))
! See Numerical Recipes, pp. 148-150.
IMPLICIT NONE
real pi
integer ncoef, n
parameter (pi=3.14159265)
real c(ncoef),f(n)
! local variables
real fac, thesum
integer j, k
fac=2./n
do j=1,ncoef
thesum=0
do k=1,n
thesum=thesum+f(k)*cos((pi*(j-1))*((k-0.5)/n))
enddo
c(j)=fac*thesum
enddo
return
end subroutine chebft
!*************************************************************************
subroutine binterp(table,km,im,jm,x,y,xtab,ytab,ix,jy,t,u,out)
! bilinear interpolation of table
!
implicit none
integer im,jm,km
real table(km,im,jm),xtab(im),ytab(jm),out(km)
integer i,ix,ip1,j,jy,jp1,k
real x,dx,t,y,dy,u,tu, tuc,tcu,tcuc
if(ix.gt.0)go to 30
if(im.gt.1)then
do i=1,im
if(x.lt.xtab(i))go to 10
enddo
10 ix=max0(i-1,1)
ip1=min0(ix+1,im)
dx=(xtab(ip1)-xtab(ix))
if(abs(dx).gt.1.e-20)then
t=(x-xtab(ix))/(xtab(ix+1)-xtab(ix))
else
t=0
endif
else
ix=1
ip1=1
t=0
endif
if(jm.gt.1)then
do j=1,jm
if(y.lt.ytab(j))go to 20
enddo
20 jy=max0(j-1,1)
jp1=min0(jy+1,jm)
dy=(ytab(jp1)-ytab(jy))
if(abs(dy).gt.1.e-20)then
u=(y-ytab(jy))/dy
else
u=0
endif
else
jy=1
jp1=1
u=0
endif
30 continue
jp1=min(jy+1,jm)
ip1=min(ix+1,im)
tu=t*u
tuc=t-tu
tcuc=1-tuc-u
tcu=u-tu
do k=1,km
out(k)=tcuc*table(k,ix,jy)+tuc*table(k,ip1,jy) &
+tu*table(k,ip1,jp1)+tcu*table(k,ix,jp1)
enddo
return
end subroutine binterp
!***************************************************************
subroutine miev0 ( xx, crefin, perfct, mimcut, anyang, &
numang, xmu, nmom, ipolzn, momdim, prnt, &
qext, qsca, gqsc, pmom, sforw, sback, s1, &
s2, tforw, tback )
!
! computes mie scattering and extinction efficiencies; asymmetry
! factor; forward- and backscatter amplitude; scattering
! amplitudes for incident polarization parallel and perpendicular
! to the plane of scattering, as functions of scattering angle;
! coefficients in the legendre polynomial expansions of either the
! unpolarized phase function or the polarized phase matrix;
! and some quantities needed in polarized radiative transfer.
!
! calls : biga, ckinmi, small1, small2, testmi, miprnt,
! lpcoef, errmsg
!
! i n t e r n a l v a r i a b l e s
! -----------------------------------
!
! an,bn mie coefficients little-a-sub-n, little-b-sub-n
! ( ref. 1, eq. 16 )
! anm1,bnm1 mie coefficients little-a-sub-(n-1),
! little-b-sub-(n-1); used in -gqsc- sum
! anp coeffs. in s+ expansion ( ref. 2, p. 1507 )
! bnp coeffs. in s- expansion ( ref. 2, p. 1507 )
! anpm coeffs. in s+ expansion ( ref. 2, p. 1507 )
! when mu is replaced by - mu
! bnpm coeffs. in s- expansion ( ref. 2, p. 1507 )
! when mu is replaced by - mu
! calcmo(k) true, calculate moments for k-th phase quantity
! (derived from -ipolzn-; used only in 'lpcoef')
! cbiga(n) bessel function ratio capital-a-sub-n (ref. 2, eq. 2)
! ( complex version )
! cior complex index of refraction with negative
! imaginary part (van de hulst convention)
! cioriv 1 / cior
! coeff ( 2n + 1 ) / ( n ( n + 1 ) )
! fn floating point version of index in loop performing
! mie series summation
! lita,litb(n) mie coefficients -an-, -bn-, saved in arrays for
! use in calculating legendre moments *pmom*
! maxtrm max. possible no. of terms in mie series
! mm + 1 and - 1, alternately.
! mim magnitude of imaginary refractive index
! mre real part of refractive index
! maxang max. possible value of input variable -numang-
! nangd2 (numang+1)/2 ( no. of angles in 0-90 deg; anyang=f )
! noabs true, sphere non-absorbing (determined by -mimcut-)
! np1dn ( n + 1 ) / n
! npquan highest-numbered phase quantity for which moments are
! to be calculated (the largest digit in -ipolzn-
! if ipolzn .ne. 0)
! ntrm no. of terms in mie series
! pass1 true on first entry, false thereafter; for self-test
! pin(j) angular function little-pi-sub-n ( ref. 2, eq. 3 )
! at j-th angle
! pinm1(j) little-pi-sub-(n-1) ( see -pin- ) at j-th angle
! psinm1 ricatti-bessel function psi-sub-(n-1), argument -xx-
! psin ricatti-bessel function psi-sub-n of argument -xx-
! ( ref. 1, p. 11 ff. )
! rbiga(n) bessel function ratio capital-a-sub-n (ref. 2, eq. 2)
! ( real version, for when imag refrac index = 0 )
! rioriv 1 / mre
! rn 1 / n
! rtmp (real) temporary variable
! sp(j) s+ for j-th angle ( ref. 2, p. 1507 )
! sm(j) s- for j-th angle ( ref. 2, p. 1507 )
! sps(j) s+ for (numang+1-j)-th angle ( anyang=false )
! sms(j) s- for (numang+1-j)-th angle ( anyang=false )
! taun angular function little-tau-sub-n ( ref. 2, eq. 4 )
! at j-th angle
! tcoef n ( n+1 ) ( 2n+1 ) (for summing tforw,tback series)
! twonp1 2n + 1
! yesang true if scattering amplitudes are to be calculated
! zetnm1 ricatti-bessel function zeta-sub-(n-1) of argument
! -xx- ( ref. 2, eq. 17 )
! zetn ricatti-bessel function zeta-sub-n of argument -xx-
!
! ----------------------------------------------------------------------
! -------- i / o specifications for subroutine miev0 -----------------
! ----------------------------------------------------------------------
implicit none
logical anyang, perfct, prnt(*)
integer ipolzn, momdim, numang, nmom
real*8 gqsc, mimcut, pmom( 0:momdim, * ), qext, qsca, &
xmu(*), xx
complex*16 crefin, sforw, sback, s1(*), s2(*), tforw(*), &
tback(*)
integer maxang,mxang2,maxtrm
real*8 onethr
! ----------------------------------------------------------------------
!
parameter ( maxang = 501, mxang2 = maxang/2 + 1 )
!
! ** note -- maxtrm = 10100 is neces-
! ** sary to do some of the test probs,
! ** but 1100 is sufficient for most
! ** conceivable applications
parameter ( maxtrm = 1100 )
parameter ( onethr = 1./3. )
!
logical anysav, calcmo(4), noabs, ok, pass1, persav, yesang
integer npquan
integer i,j,n,nmosav,iposav,numsav,ntrm,nangd2
real*8 mim, mimsav, mre, mm, np1dn
real*8 rioriv,xmusav,xxsav,sq,fn,rn,twonp1,tcoef, coeff
real*8 xinv,psinm1,chinm1,psin,chin,rtmp,taun
real*8 rbiga( maxtrm ), pin( maxang ), pinm1( maxang )
complex*16 an, bn, anm1, bnm1, anp, bnp, anpm, bnpm, cresav, &
cior, cioriv, ctmp, zet, zetnm1, zetn
complex*16 cbiga( maxtrm ), lita( maxtrm ), litb( maxtrm ), &
sp( maxang ), sm( maxang ), sps( mxang2 ), sms( mxang2 )
equivalence ( cbiga, rbiga )
save pass1
sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
data pass1 / .true. /
!
!
if ( pass1 ) then
! ** save certain user input values
xxsav = xx
cresav = crefin
mimsav = mimcut
persav = perfct
anysav = anyang
nmosav = nmom
iposav = ipolzn
numsav = numang
xmusav = xmu( 1 )
! ** reset input values for test case
xx = 10.0
crefin = ( 1.5, - 0.1 )
perfct = .false.
mimcut = 0.0
anyang = .true.
numang = 1
xmu( 1 )= - 0.7660444
nmom = 1
ipolzn = - 1
!
end if
! ** check input and calculate
! ** certain variables from input
!
10 call ckinmi( numang, maxang, xx, perfct, crefin, momdim, &
nmom, ipolzn, anyang, xmu, calcmo, npquan )
!
if ( perfct .and. xx .le. 0.1 ) then
! ** use totally-reflecting
! ** small-particle limit
!
call small1 ( xx, numang, xmu, qext, qsca, gqsc, sforw, &
sback, s1, s2, tforw, tback, lita, litb )
ntrm = 2
go to 200
end if
!
if ( .not.perfct ) then
!
cior = crefin
if ( dimag( cior ) .gt. 0.0 ) cior = dconjg( cior )
mre = dble( cior )
mim = - dimag( cior )
noabs = mim .le. mimcut
cioriv = 1.0 / cior
rioriv = 1.0 / mre
!
if ( xx * dmax1( 1.d0, cdabs(cior) ) .le. 0.d1 ) then
!
! ** use general-refractive-index
! ** small-particle limit
! ** ( ref. 2, p. 1508 )
!
call small2 ( xx, cior, .not.noabs, numang, xmu, qext, &
qsca, gqsc, sforw, sback, s1, s2, tforw, &
tback, lita, litb )
ntrm = 2
go to 200
end if
!
end if
!
nangd2 = ( numang + 1 ) / 2
yesang = numang .gt. 0
! ** estimate number of terms in mie series
! ** ( ref. 2, p. 1508 )
if ( xx.le.8.0 ) then
ntrm = xx + 4. * xx**onethr + 1.
else if ( xx.lt.4200. ) then
ntrm = xx + 4.05 * xx**onethr + 2.
else
ntrm = xx + 4. * xx**onethr + 2.
end if
if ( ntrm+1 .gt. maxtrm ) &
call errmsg( 'miev0--parameter maxtrm too small', .true. )
!
! ** calculate logarithmic derivatives of
! ** j-bessel-fcn., big-a-sub-(1 to ntrm)
if ( .not.perfct ) &
call biga( cior, xx, ntrm, noabs, yesang, rbiga, cbiga )
!
! ** initialize ricatti-bessel functions
! ** (psi,chi,zeta)-sub-(0,1) for upward
! ** recurrence ( ref. 1, eq. 19 )
xinv = 1.0 / xx
psinm1 = dsin( xx )
chinm1 = dcos( xx )
psin = psinm1 * xinv - chinm1
chin = chinm1 * xinv + psinm1
zetnm1 = dcmplx( psinm1, chinm1 )
zetn = dcmplx( psin, chin )
! ** initialize previous coeffi-
! ** cients for -gqsc- series
anm1 = ( 0.0, 0.0 )
bnm1 = ( 0.0, 0.0 )
! ** initialize angular function little-pi
! ** and sums for s+, s- ( ref. 2, p. 1507 )
if ( anyang ) then
do 60 j = 1, numang
pinm1( j ) = 0.0
pin( j ) = 1.0
sp ( j ) = ( 0.0, 0.0 )
sm ( j ) = ( 0.0, 0.0 )
60 continue
else
do 70 j = 1, nangd2
pinm1( j ) = 0.0
pin( j ) = 1.0
sp ( j ) = ( 0.0, 0.0 )
sm ( j ) = ( 0.0, 0.0 )
sps( j ) = ( 0.0, 0.0 )
sms( j ) = ( 0.0, 0.0 )
70 continue
end if
! ** initialize mie sums for efficiencies, etc.
qsca = 0.0
gqsc = 0.0
sforw = ( 0., 0. )
sback = ( 0., 0. )
tforw( 1 ) = ( 0., 0. )
tback( 1 ) = ( 0., 0. )
!
!
! --------- loop to sum mie series -----------------------------------
!
mm = + 1.0
do 100 n = 1, ntrm
! ** compute various numerical coefficients
fn = n
rn = 1.0 / fn
np1dn = 1.0 + rn
twonp1 = 2 * n + 1
coeff = twonp1 / ( fn * ( n + 1 ) )
tcoef = twonp1 * ( fn * ( n + 1 ) )
!
! ** calculate mie series coefficients
if ( perfct ) then
! ** totally-reflecting case
!
an = ( ( fn*xinv ) * psin - psinm1 ) / &
( ( fn*xinv ) * zetn - zetnm1 )
bn = psin / zetn
!
else if ( noabs ) then
! ** no-absorption case
!
an = ( ( rioriv*rbiga(n) + ( fn*xinv ) ) * psin - psinm1 ) &
/ ( ( rioriv*rbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
bn = ( ( mre * rbiga(n) + ( fn*xinv ) ) * psin - psinm1 ) &
/ ( ( mre * rbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
else
! ** absorptive case
!
an = ( ( cioriv * cbiga(n) + ( fn*xinv ) ) * psin - psinm1 ) &
/( ( cioriv * cbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
bn = ( ( cior * cbiga(n) + ( fn*xinv ) ) * psin - psinm1 ) &
/( ( cior * cbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
qsca = qsca + twonp1 * ( sq( an ) + sq( bn ) )
!
end if
! ** save mie coefficients for *pmom* calculation
lita( n ) = an
litb( n ) = bn
! ** increment mie sums for non-angle-
! ** dependent quantities
!
sforw = sforw + twonp1 * ( an + bn )
tforw( 1 ) = tforw( 1 ) + tcoef * ( an - bn )
sback = sback + ( mm * twonp1 ) * ( an - bn )
tback( 1 ) = tback( 1 ) + ( mm * tcoef ) * ( an + bn )
gqsc = gqsc + ( fn - rn ) * dble( anm1 * dconjg( an ) &
+ bnm1 * dconjg( bn ) ) &
+ coeff * dble( an * dconjg( bn ) )
!
if ( yesang ) then
! ** put mie coefficients in form
! ** needed for computing s+, s-
! ** ( ref. 2, p. 1507 )
anp = coeff * ( an + bn )
bnp = coeff * ( an - bn )
! ** increment mie sums for s+, s-
! ** while upward recursing
! ** angular functions little pi
! ** and little tau
if ( anyang ) then
! ** arbitrary angles
!
! ** vectorizable loop
do 80 j = 1, numang
rtmp = ( xmu( j ) * pin( j ) ) - pinm1( j )
taun = fn * rtmp - pinm1( j )
sp( j ) = sp( j ) + anp * ( pin( j ) + taun )
sm( j ) = sm( j ) + bnp * ( pin( j ) - taun )
pinm1( j ) = pin( j )
pin( j ) = ( xmu( j ) * pin( j ) ) + np1dn * rtmp
80 continue
!
else
! ** angles symmetric about 90 degrees
anpm = mm * anp
bnpm = mm * bnp
! ** vectorizable loop
do 90 j = 1, nangd2
rtmp = ( xmu( j ) * pin( j ) ) - pinm1( j )
taun = fn * rtmp - pinm1( j )
sp ( j ) = sp ( j ) + anp * ( pin( j ) + taun )
sms( j ) = sms( j ) + bnpm * ( pin( j ) + taun )
sm ( j ) = sm ( j ) + bnp * ( pin( j ) - taun )
sps( j ) = sps( j ) + anpm * ( pin( j ) - taun )
pinm1( j ) = pin( j )
pin( j ) = ( xmu( j ) * pin( j ) ) + np1dn * rtmp
90 continue
!
end if
end if
! ** update relevant quantities for next
! ** pass through loop
mm = - mm
anm1 = an
bnm1 = bn
! ** upward recurrence for ricatti-bessel
! ** functions ( ref. 1, eq. 17 )
!
zet = ( twonp1 * xinv ) * zetn - zetnm1
zetnm1 = zetn
zetn = zet
psinm1 = psin
psin = dble( zetn )
100 continue
!
! ---------- end loop to sum mie series --------------------------------
!
!
qext = 2. / xx**2 * dble( sforw )
if ( perfct .or. noabs ) then
qsca = qext
else
qsca = 2. / xx**2 * qsca
end if
!
gqsc = 4. / xx**2 * gqsc
sforw = 0.5 * sforw
sback = 0.5 * sback
tforw( 2 ) = 0.5 * ( sforw + 0.25 * tforw( 1 ) )
tforw( 1 ) = 0.5 * ( sforw - 0.25 * tforw( 1 ) )
tback( 2 ) = 0.5 * ( sback + 0.25 * tback( 1 ) )
tback( 1 ) = 0.5 * ( - sback + 0.25 * tback( 1 ) )
!
if ( yesang ) then
! ** recover scattering amplitudes
! ** from s+, s- ( ref. 1, eq. 11 )
if ( anyang ) then
! ** vectorizable loop
do 110 j = 1, numang
s1( j ) = 0.5 * ( sp( j ) + sm( j ) )
s2( j ) = 0.5 * ( sp( j ) - sm( j ) )
110 continue
!
else
! ** vectorizable loop
do 120 j = 1, nangd2
s1( j ) = 0.5 * ( sp( j ) + sm( j ) )
s2( j ) = 0.5 * ( sp( j ) - sm( j ) )
120 continue
! ** vectorizable loop
do 130 j = 1, nangd2
s1( numang+1 - j ) = 0.5 * ( sps( j ) + sms( j ) )
s2( numang+1 - j ) = 0.5 * ( sps( j ) - sms( j ) )
130 continue
end if
!
end if
! ** calculate legendre moments
200 if ( nmom.gt.0 ) &
call lpcoef ( ntrm, nmom, ipolzn, momdim, calcmo, npquan, &
lita, litb, pmom )
!
if ( dimag(crefin) .gt. 0.0 ) then
! ** take complex conjugates
! ** of scattering amplitudes
sforw = dconjg( sforw )
sback = dconjg( sback )
do 210 i = 1, 2
tforw( i ) = dconjg( tforw(i) )
tback( i ) = dconjg( tback(i) )
210 continue
!
do 220 j = 1, numang
s1( j ) = dconjg( s1(j) )
s2( j ) = dconjg( s2(j) )
220 continue
!
end if
!
if ( pass1 ) then
! ** compare test case results with
! ** correct answers and abort if bad
!
call testmi ( qext, qsca, gqsc, sforw, sback, s1, s2, &
tforw, tback, pmom, momdim, ok )
if ( .not. ok ) then
prnt(1) = .false.
prnt(2) = .false.
call miprnt( prnt, xx, perfct, crefin, numang, xmu, qext, &
qsca, gqsc, nmom, ipolzn, momdim, calcmo, &
pmom, sforw, sback, tforw, tback, s1, s2 )
call errmsg( 'miev0 -- self-test failed', .true. )
end if
! ** restore user input values
xx = xxsav
crefin = cresav
mimcut = mimsav
perfct = persav
anyang = anysav
nmom = nmosav
ipolzn = iposav
numang = numsav
xmu(1) = xmusav
pass1 = .false.
go to 10
!
end if
!
if ( prnt(1) .or. prnt(2) ) &
call miprnt( prnt, xx, perfct, crefin, numang, xmu, qext, &
qsca, gqsc, nmom, ipolzn, momdim, calcmo, &
pmom, sforw, sback, tforw, tback, s1, s2 )
!
return
!
end subroutine miev0
!****************************************************************************
subroutine ckinmi( numang, maxang, xx, perfct, crefin, momdim, &
nmom, ipolzn, anyang, xmu, calcmo, npquan )
!
! check for bad input to 'miev0' and calculate -calcmo,npquan-
!
implicit none
logical perfct, anyang, calcmo(*)
integer numang, maxang, momdim, nmom, ipolzn, npquan
real*8 xx, xmu(*)
integer i,l,j,ip
complex*16 crefin
!
character*4 string
logical inperr
!
inperr = .false.
!
if ( numang.gt.maxang ) then
call errmsg( 'miev0--parameter maxang too small', .true. )
inperr = .true.
end if
if ( numang.lt.0 ) call wrtbad( 'numang', inperr )
if ( xx.lt.0. ) call wrtbad( 'xx', inperr )
if ( .not.perfct .and. dble(crefin).le.0. ) &
call wrtbad( 'crefin', inperr )
if ( momdim.lt.1 ) call wrtbad( 'momdim', inperr )
!
if ( nmom.ne.0 ) then
if ( nmom.lt.0 .or. nmom.gt.momdim ) call wrtbad('nmom',inperr)
if ( iabs(ipolzn).gt.4444 ) call wrtbad( 'ipolzn', inperr )
npquan = 0
do 5 l = 1, 4
calcmo( l ) = .false.
5 continue
if ( ipolzn.ne.0 ) then
! ** parse out -ipolzn- into its digits
! ** to find which phase quantities are
! ** to have their moments calculated
!
write( string, '(i4)' ) iabs(ipolzn)
do 10 j = 1, 4
ip = ichar( string(j:j) ) - ichar( '0' )
if ( ip.ge.1 .and. ip.le.4 ) calcmo( ip ) = .true.
if ( ip.eq.0 .or. (ip.ge.5 .and. ip.le.9) ) &
call wrtbad( 'ipolzn', inperr )
npquan = max0( npquan, ip )
10 continue
end if
end if
!
if ( anyang ) then
! ** allow for slight imperfections in
! ** computation of cosine
do 20 i = 1, numang
if ( xmu(i).lt.-1.00001 .or. xmu(i).gt.1.00001 ) &
call wrtbad( 'xmu', inperr )
20 continue
else
do 22 i = 1, ( numang + 1 ) / 2
if ( xmu(i).lt.-0.00001 .or. xmu(i).gt.1.00001 ) &
call wrtbad( 'xmu', inperr )
22 continue
end if
!
if ( inperr ) &
call errmsg( 'miev0--input error(s). aborting...', .true. )
!
if ( xx.gt.20000.0 .or. dble(crefin).gt.10.0 .or. &
dabs( dimag(crefin) ).gt.10.0 ) call errmsg( &
'miev0--xx or crefin outside tested range', .false. )
!
return
end subroutine ckinmi
!***********************************************************************
subroutine lpcoef ( ntrm, nmom, ipolzn, momdim, calcmo, npquan, &
a, b, pmom )
!
! calculate legendre polynomial expansion coefficients (also
! called moments) for phase quantities ( ref. 5 formulation )
!
! input: ntrm number terms in mie series
! nmom, ipolzn, momdim 'miev0' arguments
! calcmo flags calculated from -ipolzn-
! npquan defined in 'miev0'
! a, b mie series coefficients
!
! output: pmom legendre moments ('miev0' argument)
!
! *** notes ***
!
! (1) eqs. 2-5 are in error in dave, appl. opt. 9,
! 1888 (1970). eq. 2 refers to m1, not m2; eq. 3 refers to
! m2, not m1. in eqs. 4 and 5, the subscripts on the second
! term in square brackets should be interchanged.
!
! (2) the general-case logic in this subroutine works correctly
! in the two-term mie series case, but subroutine 'lpco2t'
! is called instead, for speed.
!
! (3) subroutine 'lpco1t', to do the one-term case, is never
! called within the context of 'miev0', but is included for
! complete generality.
!
! (4) some improvement in speed is obtainable by combining the
! 310- and 410-loops, if moments for both the third and fourth
! phase quantities are desired, because the third phase quantity
! is the real part of a complex series, while the fourth phase
! quantity is the imaginary part of that very same series. but
! most users are not interested in the fourth phase quantity,
! which is related to circular polarization, so the present
! scheme is usually more efficient.
!
implicit none
logical calcmo(*)
integer ipolzn, momdim, nmom, ntrm, npquan
real*8 pmom( 0:momdim, * )
complex*16 a(*), b(*)
!
! ** specification of local variables
!
! am(m) numerical coefficients a-sub-m-super-l
! in dave, eqs. 1-15, as simplified in ref. 5.
!
! bi(i) numerical coefficients b-sub-i-super-l
! in dave, eqs. 1-15, as simplified in ref. 5.
!
! bidel(i) 1/2 bi(i) times factor capital-del in dave
!
! cm,dm() arrays c and d in dave, eqs. 16-17 (mueller form),
! calculated using recurrence derived in ref. 5
!
! cs,ds() arrays c and d in ref. 4, eqs. a5-a6 (sekera form),
! calculated using recurrence derived in ref. 5
!
! c,d() either -cm,dm- or -cs,ds-, depending on -ipolzn-
!
! evenl true for even-numbered moments; false otherwise
!
! idel 1 + little-del in dave
!
! maxtrm max. no. of terms in mie series
!
! maxmom max. no. of non-zero moments
!
! nummom number of non-zero moments
!
! recip(k) 1 / k
!
integer maxtrm,maxmom,mxmom2,maxrcp
parameter ( maxtrm = 1102, maxmom = 2*maxtrm, mxmom2 = maxmom/2, &
maxrcp = 4*maxtrm + 2 )
real*8 am( 0:maxtrm ), bi( 0:mxmom2 ), bidel( 0:mxmom2 ), &
recip( maxrcp )
complex*16 cm( maxtrm ), dm( maxtrm ), cs( maxtrm ), ds( maxtrm ), &
c( maxtrm ), d( maxtrm )
integer k,j,l,nummom,ld2,idel,m,i,mmax,imax
real*8 thesum
equivalence ( c, cm ), ( d, dm )
logical pass1, evenl
save pass1, recip
data pass1 / .true. /
!
!
if ( pass1 ) then
!
do 1 k = 1, maxrcp
recip( k ) = 1.0 / k
1 continue
pass1 = .false.
!
end if
!
do 5 j = 1, max0( 1, npquan )
do 5 l = 0, nmom
pmom( l, j ) = 0.0
5 continue
!
if ( ntrm.eq.1 ) then
call lpco1t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
return
else if ( ntrm.eq.2 ) then
call lpco2t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
return
end if
!
if ( ntrm+2 .gt. maxtrm ) &
call errmsg( 'lpcoef--parameter maxtrm too small', .true. )
!
! ** calculate mueller c, d arrays
cm( ntrm+2 ) = ( 0., 0. )
dm( ntrm+2 ) = ( 0., 0. )
cm( ntrm+1 ) = ( 1. - recip( ntrm+1 ) ) * b( ntrm )
dm( ntrm+1 ) = ( 1. - recip( ntrm+1 ) ) * a( ntrm )
cm( ntrm ) = ( recip(ntrm) + recip(ntrm+1) ) * a( ntrm ) &
+ ( 1. - recip(ntrm) ) * b( ntrm-1 )
dm( ntrm ) = ( recip(ntrm) + recip(ntrm+1) ) * b( ntrm ) &
+ ( 1. - recip(ntrm) ) * a( ntrm-1 )
!
do 10 k = ntrm-1, 2, -1
cm( k ) = cm( k+2 ) - ( 1. + recip(k+1) ) * b( k+1 ) &
+ ( recip(k) + recip(k+1) ) * a( k ) &
+ ( 1. - recip(k) ) * b( k-1 )
dm( k ) = dm( k+2 ) - ( 1. + recip(k+1) ) * a( k+1 ) &
+ ( recip(k) + recip(k+1) ) * b( k ) &
+ ( 1. - recip(k) ) * a( k-1 )
10 continue
cm( 1 ) = cm( 3 ) + 1.5 * ( a( 1 ) - b( 2 ) )
dm( 1 ) = dm( 3 ) + 1.5 * ( b( 1 ) - a( 2 ) )
!
if ( ipolzn.ge.0 ) then
!
do 20 k = 1, ntrm + 2
c( k ) = ( 2*k - 1 ) * cm( k )
d( k ) = ( 2*k - 1 ) * dm( k )
20 continue
!
else
! ** compute sekera c and d arrays
cs( ntrm+2 ) = ( 0., 0. )
ds( ntrm+2 ) = ( 0., 0. )
cs( ntrm+1 ) = ( 0., 0. )
ds( ntrm+1 ) = ( 0., 0. )
!
do 30 k = ntrm, 1, -1
cs( k ) = cs( k+2 ) + ( 2*k + 1 ) * ( cm( k+1 ) - b( k ) )
ds( k ) = ds( k+2 ) + ( 2*k + 1 ) * ( dm( k+1 ) - a( k ) )
30 continue
!
do 40 k = 1, ntrm + 2
c( k ) = ( 2*k - 1 ) * cs( k )
d( k ) = ( 2*k - 1 ) * ds( k )
40 continue
!
end if
!
!
if( ipolzn.lt.0 ) nummom = min0( nmom, 2*ntrm - 2 )
if( ipolzn.ge.0 ) nummom = min0( nmom, 2*ntrm )
if ( nummom .gt. maxmom ) &
call errmsg( 'lpcoef--parameter maxtrm too small', .true. )
!
! ** loop over moments
do 500 l = 0, nummom
ld2 = l / 2
evenl = mod( l,2 ) .eq. 0
! ** calculate numerical coefficients
! ** a-sub-m and b-sub-i in dave
! ** double-sums for moments
if( l.eq.0 ) then
!
idel = 1
do 60 m = 0, ntrm
am( m ) = 2.0 * recip( 2*m + 1 )
60 continue
bi( 0 ) = 1.0
!
else if( evenl ) then
!
idel = 1
do 70 m = ld2, ntrm
am( m ) = ( 1. + recip( 2*m-l+1 ) ) * am( m )
70 continue
do 75 i = 0, ld2-1
bi( i ) = ( 1. - recip( l-2*i ) ) * bi( i )
75 continue
bi( ld2 ) = ( 2. - recip( l ) ) * bi( ld2-1 )
!
else
!
idel = 2
do 80 m = ld2, ntrm
am( m ) = ( 1. - recip( 2*m+l+2 ) ) * am( m )
80 continue
do 85 i = 0, ld2
bi( i ) = ( 1. - recip( l+2*i+1 ) ) * bi( i )
85 continue
!
end if
! ** establish upper limits for sums
! ** and incorporate factor capital-
! ** del into b-sub-i
mmax = ntrm - idel
if( ipolzn.ge.0 ) mmax = mmax + 1
imax = min0( ld2, mmax - ld2 )
if( imax.lt.0 ) go to 600
do 90 i = 0, imax
bidel( i ) = bi( i )
90 continue
if( evenl ) bidel( 0 ) = 0.5 * bidel( 0 )
!
! ** perform double sums just for
! ** phase quantities desired by user
if( ipolzn.eq.0 ) then
!
do 110 i = 0, imax
! ** vectorizable loop (cray)
thesum = 0.0
do 100 m = ld2, mmax - i
thesum = thesum + am( m ) * &
( dble( c(m-i+1) * dconjg( c(m+i+idel) ) ) &
+ dble( d(m-i+1) * dconjg( d(m+i+idel) ) ) )
100 continue
pmom( l,1 ) = pmom( l,1 ) + bidel( i ) * thesum
110 continue
pmom( l,1 ) = 0.5 * pmom( l,1 )
go to 500
!
end if
!
if ( calcmo(1) ) then
do 160 i = 0, imax
! ** vectorizable loop (cray)
thesum = 0.0
do 150 m = ld2, mmax - i
thesum = thesum + am( m ) * &
dble( c(m-i+1) * dconjg( c(m+i+idel) ) )
150 continue
pmom( l,1 ) = pmom( l,1 ) + bidel( i ) * thesum
160 continue
end if
!
!
if ( calcmo(2) ) then
do 210 i = 0, imax
! ** vectorizable loop (cray)
thesum = 0.0
do 200 m = ld2, mmax - i
thesum = thesum + am( m ) * &
dble( d(m-i+1) * dconjg( d(m+i+idel) ) )
200 continue
pmom( l,2 ) = pmom( l,2 ) + bidel( i ) * thesum
210 continue
end if
!
!
if ( calcmo(3) ) then
do 310 i = 0, imax
! ** vectorizable loop (cray)
thesum = 0.0
do 300 m = ld2, mmax - i
thesum = thesum + am( m ) * &
( dble( c(m-i+1) * dconjg( d(m+i+idel) ) ) &
+ dble( c(m+i+idel) * dconjg( d(m-i+1) ) ) )
300 continue
pmom( l,3 ) = pmom( l,3 ) + bidel( i ) * thesum
310 continue
pmom( l,3 ) = 0.5 * pmom( l,3 )
end if
!
!
if ( calcmo(4) ) then
do 410 i = 0, imax
! ** vectorizable loop (cray)
thesum = 0.0
do 400 m = ld2, mmax - i
thesum = thesum + am( m ) * &
( dimag( c(m-i+1) * dconjg( d(m+i+idel) ) ) &
+ dimag( c(m+i+idel) * dconjg( d(m-i+1) ) ))
400 continue
pmom( l,4 ) = pmom( l,4 ) + bidel( i ) * thesum
410 continue
pmom( l,4 ) = - 0.5 * pmom( l,4 )
end if
!
500 continue
!
!
600 return
end subroutine lpcoef
!*********************************************************************
subroutine lpco1t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
!
! calculate legendre polynomial expansion coefficients (also
! called moments) for phase quantities in special case where
! no. terms in mie series = 1
!
! input: nmom, ipolzn, momdim 'miev0' arguments
! calcmo flags calculated from -ipolzn-
! a(1), b(1) mie series coefficients
!
! output: pmom legendre moments
!
implicit none
logical calcmo(*)
integer ipolzn, momdim, nmom,nummom,l
real*8 pmom( 0:momdim, * ),sq,a1sq,b1sq
complex*16 a(*), b(*), ctmp, a1b1c
sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
!
!
a1sq = sq( a(1) )
b1sq = sq( b(1) )
a1b1c = a(1) * dconjg( b(1) )
!
if( ipolzn.lt.0 ) then
!
if( calcmo(1) ) pmom( 0,1 ) = 2.25 * b1sq
if( calcmo(2) ) pmom( 0,2 ) = 2.25 * a1sq
if( calcmo(3) ) pmom( 0,3 ) = 2.25 * dble( a1b1c )
if( calcmo(4) ) pmom( 0,4 ) = 2.25 *dimag( a1b1c )
!
else
!
nummom = min0( nmom, 2 )
! ** loop over moments
do 100 l = 0, nummom
!
if( ipolzn.eq.0 ) then
if( l.eq.0 ) pmom( l,1 ) = 1.5 * ( a1sq + b1sq )
if( l.eq.1 ) pmom( l,1 ) = 1.5 * dble( a1b1c )
if( l.eq.2 ) pmom( l,1 ) = 0.15 * ( a1sq + b1sq )
go to 100
end if
!
if( calcmo(1) ) then
if( l.eq.0 ) pmom( l,1 ) = 2.25 * ( a1sq + b1sq / 3. )
if( l.eq.1 ) pmom( l,1 ) = 1.5 * dble( a1b1c )
if( l.eq.2 ) pmom( l,1 ) = 0.3 * b1sq
end if
!
if( calcmo(2) ) then
if( l.eq.0 ) pmom( l,2 ) = 2.25 * ( b1sq + a1sq / 3. )
if( l.eq.1 ) pmom( l,2 ) = 1.5 * dble( a1b1c )
if( l.eq.2 ) pmom( l,2 ) = 0.3 * a1sq
end if
!
if( calcmo(3) ) then
if( l.eq.0 ) pmom( l,3 ) = 3.0 * dble( a1b1c )
if( l.eq.1 ) pmom( l,3 ) = 0.75 * ( a1sq + b1sq )
if( l.eq.2 ) pmom( l,3 ) = 0.3 * dble( a1b1c )
end if
!
if( calcmo(4) ) then
if( l.eq.0 ) pmom( l,4 ) = - 1.5 * dimag( a1b1c )
if( l.eq.1 ) pmom( l,4 ) = 0.0
if( l.eq.2 ) pmom( l,4 ) = 0.3 * dimag( a1b1c )
end if
!
100 continue
!
end if
!
return
end subroutine lpco1t
!********************************************************************
subroutine lpco2t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
!
! calculate legendre polynomial expansion coefficients (also
! called moments) for phase quantities in special case where
! no. terms in mie series = 2
!
! input: nmom, ipolzn, momdim 'miev0' arguments
! calcmo flags calculated from -ipolzn-
! a(1-2), b(1-2) mie series coefficients
!
! output: pmom legendre moments
!
implicit none
logical calcmo(*)
integer ipolzn, momdim, nmom,l,nummom
real*8 pmom( 0:momdim, * ),sq,pm1,pm2,a2sq,b2sq
complex*16 a(*), b(*)
complex*16 a2c, b2c, ctmp, ca, cac, cat, cb, cbc, cbt, cg, ch
sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
!
!
ca = 3. * a(1) - 5. * b(2)
cat= 3. * b(1) - 5. * a(2)
cac = dconjg( ca )
a2sq = sq( a(2) )
b2sq = sq( b(2) )
a2c = dconjg( a(2) )
b2c = dconjg( b(2) )
!
if( ipolzn.lt.0 ) then
! ** loop over sekera moments
nummom = min0( nmom, 2 )
do 50 l = 0, nummom
!
if( calcmo(1) ) then
if( l.eq.0 ) pmom( l,1 ) = 0.25 * ( sq(cat) + &
(100./3.) * b2sq )
if( l.eq.1 ) pmom( l,1 ) = (5./3.) * dble( cat * b2c )
if( l.eq.2 ) pmom( l,1 ) = (10./3.) * b2sq
end if
!
if( calcmo(2) ) then
if( l.eq.0 ) pmom( l,2 ) = 0.25 * ( sq(ca) + &
(100./3.) * a2sq )
if( l.eq.1 ) pmom( l,2 ) = (5./3.) * dble( ca * a2c )
if( l.eq.2 ) pmom( l,2 ) = (10./3.) * a2sq
end if
!
if( calcmo(3) ) then
if( l.eq.0 ) pmom( l,3 ) = 0.25 * dble( cat*cac + &
(100./3.)*b(2)*a2c )
if( l.eq.1 ) pmom( l,3 ) = 5./6. * dble( b(2)*cac + &
cat*a2c )
if( l.eq.2 ) pmom( l,3 ) = 10./3. * dble( b(2) * a2c )
end if
!
if( calcmo(4) ) then
if( l.eq.0 ) pmom( l,4 ) = -0.25 * dimag( cat*cac + &
(100./3.)*b(2)*a2c )
if( l.eq.1 ) pmom( l,4 ) = -5./6. * dimag( b(2)*cac + &
cat*a2c )
if( l.eq.2 ) pmom( l,4 ) = -10./3. * dimag( b(2) * a2c )
end if
!
50 continue
!
else
!
cb = 3. * b(1) + 5. * a(2)
cbt= 3. * a(1) + 5. * b(2)
cbc = dconjg( cb )
cg = ( cbc*cbt + 10.*( cac*a(2) + b2c*cat) ) / 3.
ch = 2.*( cbc*a(2) + b2c*cbt )
!
! ** loop over mueller moments
nummom = min0( nmom, 4 )
do 100 l = 0, nummom
!
if( ipolzn.eq.0 .or. calcmo(1) ) then
if( l.eq.0 ) pm1 = 0.25 * sq(ca) + sq(cb) / 12. &
+ (5./3.) * dble(ca*b2c) + 5.*b2sq
if( l.eq.1 ) pm1 = dble( cb * ( cac/6. + b2c ) )
if( l.eq.2 ) pm1 = sq(cb)/30. + (20./7.) * b2sq &
+ (2./3.) * dble( ca * b2c )
if( l.eq.3 ) pm1 = (2./7.) * dble( cb * b2c )
if( l.eq.4 ) pm1 = (40./63.) * b2sq
if ( calcmo(1) ) pmom( l,1 ) = pm1
end if
!
if( ipolzn.eq.0 .or. calcmo(2) ) then
if( l.eq.0 ) pm2 = 0.25*sq(cat) + sq(cbt) / 12. &
+ (5./3.) * dble(cat*a2c) + 5.*a2sq
if( l.eq.1 ) pm2 = dble( cbt * ( dconjg(cat)/6. + a2c) )
if( l.eq.2 ) pm2 = sq(cbt)/30. + (20./7.) * a2sq &
+ (2./3.) * dble( cat * a2c )
if( l.eq.3 ) pm2 = (2./7.) * dble( cbt * a2c )
if( l.eq.4 ) pm2 = (40./63.) * a2sq
if ( calcmo(2) ) pmom( l,2 ) = pm2
end if
!
if( ipolzn.eq.0 ) then
pmom( l,1 ) = 0.5 * ( pm1 + pm2 )
go to 100
end if
!
if( calcmo(3) ) then
if( l.eq.0 ) pmom( l,3 ) = 0.25 * dble( cac*cat + cg + &
20.*b2c*a(2) )
if( l.eq.1 ) pmom( l,3 ) = dble( cac*cbt + cbc*cat + &
3.*ch ) / 12.
if( l.eq.2 ) pmom( l,3 ) = 0.1 * dble( cg + (200./7.) * &
b2c * a(2) )
if( l.eq.3 ) pmom( l,3 ) = dble( ch ) / 14.
if( l.eq.4 ) pmom( l,3 ) = 40./63. * dble( b2c * a(2) )
end if
!
if( calcmo(4) ) then
if( l.eq.0 ) pmom( l,4 ) = 0.25 * dimag( cac*cat + cg + &
20.*b2c*a(2) )
if( l.eq.1 ) pmom( l,4 ) = dimag( cac*cbt + cbc*cat + &
3.*ch ) / 12.
if( l.eq.2 ) pmom( l,4 ) = 0.1 * dimag( cg + (200./7.) * &
b2c * a(2) )
if( l.eq.3 ) pmom( l,4 ) = dimag( ch ) / 14.
if( l.eq.4 ) pmom( l,4 ) = 40./63. * dimag( b2c * a(2) )
end if
!
100 continue
!
end if
!
return
end subroutine lpco2t
!*********************************************************************
subroutine biga( cior, xx, ntrm, noabs, yesang, rbiga, cbiga )
!
! calculate logarithmic derivatives of j-bessel-function
!
! input : cior, xx, ntrm, noabs, yesang (defined in 'miev0')
!
! output : rbiga or cbiga (defined in 'miev0')
!
! internal variables :
!
! confra value of lentz continued fraction for -cbiga(ntrm)-,
! used to initialize downward recurrence.
! down = true, use down-recurrence. false, do not.
! f1,f2,f3 arithmetic statement functions used in determining
! whether to use up- or down-recurrence
! ( ref. 2, eqs. 6-8 )
! mre real refractive index
! mim imaginary refractive index
! rezinv 1 / ( mre * xx ); temporary variable for recurrence
! zinv 1 / ( cior * xx ); temporary variable for recurrence
!
implicit none
logical down, noabs, yesang
integer ntrm,n
real*8 mre, mim, rbiga(*), xx, rezinv, rtmp, f1,f2,f3
! complex*16 cior, ctmp, confra, cbiga(*), zinv
complex*16 cior, ctmp, cbiga(*), zinv
f1( mre ) = - 8.0 + mre**2 * ( 26.22 + mre * ( - 0.4474 &
+ mre**3 * ( 0.00204 - 0.000175 * mre ) ) )
f2( mre ) = 3.9 + mre * ( - 10.8 + 13.78 * mre )
f3( mre ) = - 15.04 + mre * ( 8.42 + 16.35 * mre )
!
! ** decide whether 'biga' can be
! ** calculated by up-recurrence
mre = dble( cior )
mim = dabs( dimag( cior ) )
if ( mre.lt.1.0 .or. mre.gt.10.0 .or. mim.gt.10.0 ) then
down = .true.
else if ( yesang ) then
down = .true.
if ( mim*xx .lt. f2( mre ) ) down = .false.
else
down = .true.
if ( mim*xx .lt. f1( mre ) ) down = .false.
end if
!
zinv = 1.0 / ( cior * xx )
rezinv = 1.0 / ( mre * xx )
!
if ( down ) then
! ** compute initial high-order 'biga' using
! ** lentz method ( ref. 1, pp. 17-20 )
!
ctmp = confra( ntrm, zinv, xx )
!
! *** downward recurrence for 'biga'
! *** ( ref. 1, eq. 22 )
if ( noabs ) then
! ** no-absorption case
rbiga( ntrm ) = dble( ctmp )
do 25 n = ntrm, 2, - 1
rbiga( n-1 ) = (n*rezinv) &
- 1.0 / ( (n*rezinv) + rbiga( n ) )
25 continue
!
else
! ** absorptive case
cbiga( ntrm ) = ctmp
do 30 n = ntrm, 2, - 1
cbiga( n-1 ) = (n*zinv) - 1.0 / ( (n*zinv) + cbiga( n ) )
30 continue
!
end if
!
else
! *** upward recurrence for 'biga'
! *** ( ref. 1, eqs. 20-21 )
if ( noabs ) then
! ** no-absorption case
rtmp = dsin( mre*xx )
rbiga( 1 ) = - rezinv &
+ rtmp / ( rtmp*rezinv - dcos( mre*xx ) )
do 40 n = 2, ntrm
rbiga( n ) = - ( n*rezinv ) &
+ 1.0 / ( ( n*rezinv ) - rbiga( n-1 ) )
40 continue
!
else
! ** absorptive case
ctmp = cdexp( - dcmplx(0.d0,2.d0) * cior * xx )
cbiga( 1 ) = - zinv + (1.-ctmp) / ( zinv * (1.-ctmp) - &
dcmplx(0.d0,1.d0)*(1.+ctmp) )
do 50 n = 2, ntrm
cbiga( n ) = - (n*zinv) + 1.0 / ((n*zinv) - cbiga( n-1 ))
50 continue
end if
!
end if
!
return
end subroutine biga
!**********************************************************************
complex*16 function confra( n, zinv, xx )
!
! compute bessel function ratio capital-a-sub-n from its
! continued fraction using lentz method ( ref. 1, pp. 17-20 )
!
! zinv = reciprocal of argument of capital-a
!
! i n t e r n a l v a r i a b l e s
! ------------------------------------
!
! cak term in continued fraction expansion of capital-a
! ( ref. 1, eq. 25 )
! capt factor used in lentz iteration for capital-a
! ( ref. 1, eq. 27 )
! cdenom denominator in -capt- ( ref. 1, eq. 28b )
! cnumer numerator in -capt- ( ref. 1, eq. 28a )
! cdtd product of two successive denominators of -capt-
! factors ( ref. 1, eq. 34c )
! cntn product of two successive numerators of -capt-
! factors ( ref. 1, eq. 34b )
! eps1 ill-conditioning criterion
! eps2 convergence criterion
! kk subscript k of -cak- ( ref. 1, eq. 25b )
! kount iteration counter ( used only to prevent runaway )
! maxit max. allowed no. of iterations
! mm + 1 and - 1, alternately
!
implicit none
integer n,maxit,mm,kk,kount
real*8 xx,eps1,eps2
complex*16 zinv
complex*16 cak, capt, cdenom, cdtd, cnumer, cntn
! data eps1 / 1.e - 2 /, eps2 / 1.e - 8 /
data eps1 / 1.d-2 /, eps2 / 1.d-8 /
data maxit / 10000 /
!
! *** ref. 1, eqs. 25a, 27
confra = ( n + 1 ) * zinv
mm = - 1
kk = 2 * n + 3
cak = ( mm * kk ) * zinv
cdenom = cak
cnumer = cdenom + 1.0 / confra
kount = 1
!
20 kount = kount + 1
if ( kount.gt.maxit ) &
call errmsg( 'confra--iteration failed to converge$', .true.)
!
! *** ref. 2, eq. 25b
mm = - mm
kk = kk + 2
cak = ( mm * kk ) * zinv
! *** ref. 2, eq. 32
if ( cdabs( cnumer/cak ).le.eps1 &
.or. cdabs( cdenom/cak ).le.eps1 ) then
!
! ** ill-conditioned case -- stride
! ** two terms instead of one
!
! *** ref. 2, eqs. 34
cntn = cak * cnumer + 1.0
cdtd = cak * cdenom + 1.0
confra = ( cntn / cdtd ) * confra
! *** ref. 2, eq. 25b
mm = - mm
kk = kk + 2
cak = ( mm * kk ) * zinv
! *** ref. 2, eqs. 35
cnumer = cak + cnumer / cntn
cdenom = cak + cdenom / cdtd
kount = kount + 1
go to 20
!
else
! ** well-conditioned case
!
! *** ref. 2, eqs. 26, 27
capt = cnumer / cdenom
confra = capt * confra
! ** check for convergence
! ** ( ref. 2, eq. 31 )
!
if ( dabs( dble(capt) - 1.0 ).ge.eps2 &
.or. dabs( dimag(capt) ) .ge.eps2 ) then
!
! *** ref. 2, eqs. 30a-b
cnumer = cak + 1.0 / cnumer
cdenom = cak + 1.0 / cdenom
go to 20
end if
end if
!
return
!
end function confra
!********************************************************************
subroutine miprnt( prnt, xx, perfct, crefin, numang, xmu, &
qext, qsca, gqsc, nmom, ipolzn, momdim, &
calcmo, pmom, sforw, sback, tforw, tback, &
s1, s2 )
!
! print scattering quantities of a single particle
!
implicit none
logical perfct, prnt(*), calcmo(*)
integer ipolzn, momdim, nmom, numang,i,m,j
real*8 gqsc, pmom( 0:momdim, * ), qext, qsca, xx, xmu(*)
real*8 fi1,fi2,fnorm
complex*16 crefin, sforw, sback, tforw(*), tback(*), s1(*), s2(*)
character*22 fmt
!
!
if ( perfct ) write ( *, '(''1'',10x,a,1p,e11.4)' ) &
'perfectly conducting case, size parameter =', xx
if ( .not.perfct ) write ( *, '(''1'',10x,3(a,1p,e11.4))' ) &
'refractive index: real ', dble(crefin), &
' imag ', dimag(crefin), ', size parameter =', xx
!
if ( prnt(1) .and. numang.gt.0 ) then
!
write ( *, '(/,a)' ) &
' cos(angle) ------- s1 --------- ------- s2 ---------'// &
' --- s1*conjg(s2) --- i1=s1**2 i2=s2**2 (i1+i2)/2'// &
' deg polzn'
do 10 i = 1, numang
fi1 = dble( s1(i) ) **2 + dimag( s1(i) )**2
fi2 = dble( s2(i) ) **2 + dimag( s2(i) )**2
write( *, '( i4, f10.6, 1p,10e11.3 )' ) &
i, xmu(i), s1(i), s2(i), s1(i)*dconjg(s2(i)), &
fi1, fi2, 0.5*(fi1+fi2), (fi2-fi1)/(fi2+fi1)
10 continue
!
end if
!
!
if ( prnt(2) ) then
!
write ( *, '(/,a,9x,a,17x,a,17x,a,/,(0p,f7.2, 1p,6e12.3) )' ) &
' angle', 's-sub-1', 't-sub-1', 't-sub-2', &
0.0, sforw, tforw(1), tforw(2), &
180., sback, tback(1), tback(2)
write ( *, '(/,4(a,1p,e11.4))' ) &
' efficiency factors, extinction:', qext, &
' scattering:', qsca, &
' absorption:', qext-qsca, &
' rad. pressure:', qext-gqsc
!
if ( nmom.gt.0 ) then
!
write( *, '(/,a)' ) ' normalized moments of :'
if ( ipolzn.eq.0 ) write ( *, '(''+'',27x,a)' ) 'phase fcn'
if ( ipolzn.gt.0 ) write ( *, '(''+'',33x,a)' ) &
'm1 m2 s21 d21'
if ( ipolzn.lt.0 ) write ( *, '(''+'',33x,a)' ) &
'r1 r2 r3 r4'
!
fnorm = 4. / ( xx**2 * qsca )
do 20 m = 0, nmom
write ( *, '(a,i4)' ) ' moment no.', m
do 20 j = 1, 4
if( calcmo(j) ) then
write( fmt, 98 ) 24 + (j-1)*13
write ( *,fmt ) fnorm * pmom(m,j)
end if
20 continue
end if
!
end if
!
return
!
98 format( '( ''+'', t', i2, ', 1p,e13.4 )' )
end subroutine miprnt
!**************************************************************************
subroutine small1 ( xx, numang, xmu, qext, qsca, gqsc, sforw, &
sback, s1, s2, tforw, tback, a, b )
!
! small-particle limit of mie quantities in totally reflecting
! limit ( mie series truncated after 2 terms )
!
! a,b first two mie coefficients, with numerator and
! denominator expanded in powers of -xx- ( a factor
! of xx**3 is missing but is restored before return
! to calling program ) ( ref. 2, p. 1508 )
!
implicit none
integer numang,j
real*8 gqsc, qext, qsca, xx, xmu(*)
real*8 twothr,fivthr,fivnin,sq,rtmp
complex*16 a( 2 ), b( 2 ), sforw, sback, s1(*), s2(*), &
tforw(*), tback(*)
!
parameter ( twothr = 2./3., fivthr = 5./3., fivnin = 5./9. )
complex*16 ctmp
sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
!
!
a( 1 ) = dcmplx ( 0.d0, twothr * ( 1. - 0.2 * xx**2 ) ) &
/ dcmplx ( 1.d0 - 0.5 * xx**2, twothr * xx**3 )
!
b( 1 ) = dcmplx ( 0.d0, - ( 1. - 0.1 * xx**2 ) / 3. ) &
/ dcmplx ( 1.d0 + 0.5 * xx**2, - xx**3 / 3. )
!
a( 2 ) = dcmplx ( 0.d0, xx**2 / 30. )
b( 2 ) = dcmplx ( 0.d0, - xx**2 / 45. )
!
qsca = 6. * xx**4 * ( sq( a(1) ) + sq( b(1) ) &
+ fivthr * ( sq( a(2) ) + sq( b(2) ) ) )
qext = qsca
gqsc = 6. * xx**4 * dble( a(1) * dconjg( a(2) + b(1) ) &
+ ( b(1) + fivnin * a(2) ) * dconjg( b(2) ) )
!
rtmp = 1.5 * xx**3
sforw = rtmp * ( a(1) + b(1) + fivthr * ( a(2) + b(2) ) )
sback = rtmp * ( a(1) - b(1) - fivthr * ( a(2) - b(2) ) )
tforw( 1 ) = rtmp * ( b(1) + fivthr * ( 2.*b(2) - a(2) ) )
tforw( 2 ) = rtmp * ( a(1) + fivthr * ( 2.*a(2) - b(2) ) )
tback( 1 ) = rtmp * ( b(1) - fivthr * ( 2.*b(2) + a(2) ) )
tback( 2 ) = rtmp * ( a(1) - fivthr * ( 2.*a(2) + b(2) ) )
!
do 10 j = 1, numang
s1( j ) = rtmp * ( a(1) + b(1) * xmu(j) + fivthr * &
( a(2) * xmu(j) + b(2) * ( 2.*xmu(j)**2 - 1. )) )
s2( j ) = rtmp * ( b(1) + a(1) * xmu(j) + fivthr * &
( b(2) * xmu(j) + a(2) * ( 2.*xmu(j)**2 - 1. )) )
10 continue
! ** recover actual mie coefficients
a( 1 ) = xx**3 * a( 1 )
a( 2 ) = xx**3 * a( 2 )
b( 1 ) = xx**3 * b( 1 )
b( 2 ) = xx**3 * b( 2 )
!
return
end subroutine small1
!*************************************************************************
subroutine small2 ( xx, cior, calcqe, numang, xmu, qext, qsca, &
gqsc, sforw, sback, s1, s2, tforw, tback, &
a, b )
!
! small-particle limit of mie quantities for general refractive
! index ( mie series truncated after 2 terms )
!
! a,b first two mie coefficients, with numerator and
! denominator expanded in powers of -xx- ( a factor
! of xx**3 is missing but is restored before return
! to calling program ) ( ref. 2, p. 1508 )
!
! ciorsq square of refractive index
!
implicit none
logical calcqe
integer numang,j
real*8 gqsc, qext, qsca, xx, xmu(*)
real*8 twothr,fivthr,sq,rtmp
complex*16 a( 2 ), b( 2 ), cior, sforw, sback, s1(*), s2(*), &
tforw(*), tback(*)
!
parameter ( twothr = 2./3., fivthr = 5./3. )
complex*16 ctmp, ciorsq
sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
!
!
ciorsq = cior**2
ctmp = dcmplx( 0.d0, twothr ) * ( ciorsq - 1.0 )
a(1) = ctmp * ( 1.0 - 0.1 * xx**2 + (ciorsq/350. + 1./280.)*xx**4) &
/ ( ciorsq + 2.0 + ( 1.0 - 0.7 * ciorsq ) * xx**2 &
- ( ciorsq**2/175. - 0.275 * ciorsq + 0.25 ) * xx**4 &
+ xx**3 * ctmp * ( 1.0 - 0.1 * xx**2 ) )
!
b(1) = (xx**2/30.) * ctmp * ( 1.0 + (ciorsq/35. - 1./14.) *xx**2 ) &
/ ( 1.0 - ( ciorsq/15. - 1./6. ) * xx**2 )
!
a(2) = ( 0.1 * xx**2 ) * ctmp * ( 1.0 - xx**2 / 14. ) &
/ ( 2. * ciorsq + 3. - ( ciorsq/7. - 0.5 ) * xx**2 )
!
qsca = 6. * xx**4 * ( sq(a(1)) + sq(b(1)) + fivthr * sq(a(2)) )
gqsc = 6. * xx**4 * dble( a(1) * dconjg( a(2) + b(1) ) )
qext = qsca
if ( calcqe ) qext = 6. * xx * dble( a(1) + b(1) + fivthr * a(2) )
!
rtmp = 1.5 * xx**3
sforw = rtmp * ( a(1) + b(1) + fivthr * a(2) )
sback = rtmp * ( a(1) - b(1) - fivthr * a(2) )
tforw( 1 ) = rtmp * ( b(1) - fivthr * a(2) )
tforw( 2 ) = rtmp * ( a(1) + 2. * fivthr * a(2) )
tback( 1 ) = tforw( 1 )
tback( 2 ) = rtmp * ( a(1) - 2. * fivthr * a(2) )
!
do 10 j = 1, numang
s1( j ) = rtmp * ( a(1) + ( b(1) + fivthr * a(2) ) * xmu(j) )
s2( j ) = rtmp * ( b(1) + a(1) * xmu(j) + fivthr * a(2) &
* ( 2. * xmu(j)**2 - 1. ) )
10 continue
! ** recover actual mie coefficients
a( 1 ) = xx**3 * a( 1 )
a( 2 ) = xx**3 * a( 2 )
b( 1 ) = xx**3 * b( 1 )
b( 2 ) = ( 0., 0. )
!
return
end subroutine small2
!***********************************************************************
subroutine testmi ( qext, qsca, gqsc, sforw, sback, s1, s2, &
tforw, tback, pmom, momdim, ok )
!
! compare mie code test case results with correct answers
! and return ok=false if even one result is inaccurate.
!
! the test case is : mie size parameter = 10
! refractive index = 1.5 - 0.1 i
! scattering angle = 140 degrees
! 1 sekera moment
!
! results for this case may be found among the test cases
! at the end of reference (1).
!
! *** note *** when running on some computers, esp. in single
! precision, the 'accur' criterion below may have to be relaxed.
! however, if 'accur' must be set larger than 10**-3 for some
! size parameters, your computer is probably not accurate
! enough to do mie computations for those size parameters.
!
implicit none
integer momdim,m,n
real*8 qext, qsca, gqsc, pmom( 0:momdim, * )
complex*16 sforw, sback, s1(*), s2(*), tforw(*), tback(*)
logical ok, wrong
!
real*8 accur, testqe, testqs, testgq, testpm( 0:1 )
complex*16 testsf, testsb,tests1,tests2,testtf(2), testtb(2)
data testqe / 2.459791 /, testqs / 1.235144 /, &
testgq / 1.139235 /, testsf / ( 61.49476, -3.177994 ) /, &
testsb / ( 1.493434, 0.2963657 ) /, &
tests1 / ( -0.1548380, -1.128972) /, &
tests2 / ( 0.05669755, 0.5425681) /, &
testtf / ( 12.95238, -136.6436 ), ( 48.54238, 133.4656 ) /, &
testtb / ( 41.88414, -15.57833 ), ( 43.37758, -15.28196 )/, &
testpm / 227.1975, 183.6898 /
real*8 calc,exact
! data accur / 1.e-5 /
data accur / 1.e-4 /
wrong( calc, exact ) = dabs( (calc - exact) / exact ) .gt. accur
!
!
ok = .true.
if ( wrong( qext,testqe ) ) &
call tstbad( 'qext', abs((qext - testqe) / testqe), ok )
if ( wrong( qsca,testqs ) ) &
call tstbad( 'qsca', abs((qsca - testqs) / testqs), ok )
if ( wrong( gqsc,testgq ) ) &
call tstbad( 'gqsc', abs((gqsc - testgq) / testgq), ok )
!
if ( wrong( dble(sforw), dble(testsf) ) .or. &
wrong( dimag(sforw), dimag(testsf) ) ) &
call tstbad( 'sforw', cdabs((sforw - testsf) / testsf), ok )
!
if ( wrong( dble(sback), dble(testsb) ) .or. &
wrong( dimag(sback), dimag(testsb) ) ) &
call tstbad( 'sback', cdabs((sback - testsb) / testsb), ok )
!
if ( wrong( dble(s1(1)), dble(tests1) ) .or. &
wrong( dimag(s1(1)), dimag(tests1) ) ) &
call tstbad( 's1', cdabs((s1(1) - tests1) / tests1), ok )
!
if ( wrong( dble(s2(1)), dble(tests2) ) .or. &
wrong( dimag(s2(1)), dimag(tests2) ) ) &
call tstbad( 's2', cdabs((s2(1) - tests2) / tests2), ok )
!
do 20 n = 1, 2
if ( wrong( dble(tforw(n)), dble(testtf(n)) ) .or. &
wrong( dimag(tforw(n)), dimag(testtf(n)) ) ) &
call tstbad( 'tforw', cdabs( (tforw(n) - testtf(n)) / &
testtf(n) ), ok )
if ( wrong( dble(tback(n)), dble(testtb(n)) ) .or. &
wrong( dimag(tback(n)), dimag(testtb(n)) ) ) &
call tstbad( 'tback', cdabs( (tback(n) - testtb(n)) / &
testtb(n) ), ok )
20 continue
!
do 30 m = 0, 1
if ( wrong( pmom(m,1), testpm(m) ) ) &
call tstbad( 'pmom', dabs( (pmom(m,1)-testpm(m)) / &
testpm(m) ), ok )
30 continue
!
return
!
end subroutine testmi
!**************************************************************************
subroutine errmsg( messag, fatal )
!
! print out a warning or error message; abort if error
!
USE module_peg_util, only: peg_message, peg_error_fatal
implicit none
logical fatal, once
character*80 msg
character*(*) messag
integer maxmsg, nummsg
save maxmsg, nummsg, once
data nummsg / 0 /, maxmsg / 100 /, once / .false. /
!
!
if ( fatal ) then
! write ( *, '(2a)' ) ' ******* error >>>>>> ', messag
! stop
write( msg, '(a)' ) &
'FASTJ mie fatal error ' // &
messag
call peg_message( lunerr, msg )
call peg_error_fatal( lunerr, msg )
end if
!
nummsg = nummsg + 1
if ( nummsg.gt.maxmsg ) then
! if ( .not.once ) write ( *,99 )
if ( .not.once )then
write( msg, '(a)' ) &
'FASTJ mie: too many warning messages -- no longer printing '
call peg_message( lunerr, msg )
end if
once = .true.
else
msg = 'FASTJ mie warning ' // messag
call peg_message( lunerr, msg )
! write ( *, '(2a)' ) ' ******* warning >>>>>> ', messag
endif
!
return
!
! 99 format( ///,' >>>>>> too many warning messages -- ', &
! 'they will no longer be printed <<<<<<<', /// )
end subroutine errmsg
!********************************************************************
subroutine wrtbad ( varnam, erflag )
!
! write names of erroneous variables
!
! input : varnam = name of erroneous variable to be written
! ( character, any length )
!
! output : erflag = logical flag, set true by this routine
! ----------------------------------------------------------------------
USE module_peg_util, only: peg_message
implicit none
character*(*) varnam
logical erflag
integer maxmsg, nummsg
character*80 msg
save nummsg, maxmsg
data nummsg / 0 /, maxmsg / 50 /
!
!
nummsg = nummsg + 1
! write ( *, '(3a)' ) ' **** input variable ', varnam, &
! ' in error ****'
msg = 'FASTJ mie input variable in error ' // varnam
call peg_message( lunerr, msg )
erflag = .true.
if ( nummsg.eq.maxmsg ) &
call errmsg ( 'too many input variable errors. aborting...$', .true. )
return
!
end subroutine wrtbad
!******************************************************************
subroutine tstbad( varnam, relerr, ok )
!
! write name (-varnam-) of variable failing self-test and its
! percent error from the correct value. return ok = false.
!
implicit none
character*(*) varnam
logical ok
real*8 relerr
!
!
ok = .false.
write( *, '(/,3a,1p,e11.2,a)' ) &
' output variable ', varnam,' differed by', 100.*relerr, &
' per cent from correct value. self-test failed.'
return
!
end subroutine tstbad
!-----------------------------------------------------------------------
end module module_fastj_mie