chem/module_fastj_mie.F

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