!*********************************************************************************
! This computer software was prepared by Battelle Memorial Institute, hereinafter
! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of 
! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
!
! Module to Compute Aerosol Optical Properties
! * Author: Jerome D. Fast
! * Originators of parts of code:
!   Rahul A. Zaveri, Jim Barnard, Richard C. Easter, William I. Gustafson Jr.
! Last update: April 2007
!
! Contact:
! Jerome D. Fast, PhD
! Staff Scientist
! Pacific Northwest National Laboratory
! P.O. Box 999, MSIN K9-30
! Richland, WA, 99352
! Phone: (509) 372-6116
! Email: Jerome.Fast@pnl.gov
!
! Please report any bugs or problems to Jerome Fast, the WRF-chem implmentation
! team leader for PNNL
!
! Terms of Use:
!  1) This module may not be included in  commerical package or used for any
!     commercial applications without the consent of the PNNL contact. 
!  2) This module is provided to the WRF modeling community; however, no portion
!     of it can be used separately or in another code without the consent of the
!     PNNL contact.
!  3) This module may be used for research, educational, and non-profit purposes
!     only.  Any other usage must be first approved by the PNNL contact.
!  4) Publications resulting from the usage of this module must use one the
!     reference below for proper acknowledgment.
!
! References: 
! * Fast, J.D., W.I. Gustafson Jr., R.C. Easter, R.A. Zaveri, J.C. Barnard, E.G.
!   Chapman, G.A. Grell, and S.E. Peckham (2005), Evolution of ozone, particulates,
!   and aerosol direct radiative forcing in the vicinity of Houston using a fully- 
!   coupled meteorology-chemistry-aerosol model. JGR, 111, doi:10.1029/2005JD006721.
!
! Contact Jerome Fast for updates on the status of manuscripts under review.  
!
! Additional information:
! *  www.pnl.gov/atmos_sciences/Jdf/wrfchem.html
!
! Support:
! Funding for this code development was provided by the U.S. Department of Energy
! under the auspices of Atmospheric Science Program of the Office of Biological
! and Environmental Research the PNNL Laboratory Research and Directed Research
! and Development program.
!**********************************************************************************
	module module_optical_averaging
	
        implicit none
 	integer, parameter :: lunerr = -1
	
	contains

!----------------------------------------------------------------------------------
! Aerosol optical properties computed using three methods (option_method):
! 1) volume averaging mixing rule: method that assumes internal-mixing of aerosol
!    composition that averages the refractive indices for each size bin
! 2) Maxwell-Garnett mixing rule: method that randomly distributes black carbon
!    within a particle
! 3) shell-core: method that assumes a "core" composed of black carbon surrounded 
!    by a "shell" composed of all other compositions
!
! There are two Mie routines included (option_mie):
! 1) subroutine mieaer: Employs a Chebyshev economization (Fast et al. 2006, Ghan
!    et al. (2001) so that full Mie computations are called only once and then
!    expansion coeffiecients are used for subsequent times to save CPU.  This
!    method is somewhat less accurate than full Mie calculation.
! 2) subroutine mieaer_sc: Full Mie calculation at each time step that also 
!    permits computation of shell-core method.
! 
! Sectional and modal size distributions are treated similary, but there is 
! separate code currrently to handle differences between MOSAIC and MADE/SORGAM.
!
! Methodology for sectional:
! * 3-D arrays for refractive index, wet radius, and aerosol number produced by
!   optical_prep_sectional are then passed into mieaer_sectional
! * subroutine mieaer produces vertical profiles of aerosol optical properties for
!   4 wavelengths that are put into 3-D arrays and passed back up to chem_driver.F
! * tauaer*, waer*, gaer* passed to module_ra_gsfcsw.F
! * tauaer*, waer*, gaer*, l2-l7 passed to module_phot_fastj.F
! Methodology for modal:
! * similar to sectional, except divide modal mass into discrete size bins first
! * currently assume same 8 size bins as MOSAIC, but other bins are possible
!   
! THIS CODE IS STILL BEING TESTED.  USERS ARE ENCOURAGED TO USE ONLY
! AER_OP_OPT=1
!
      subroutine optical_averaging(id,ktau,dtstep,config_flags,         &
                 nbin_o,haveaer,option_method,option_mie,chem,dz8w,alt, &
                 h2oai,h2oaj,                                           &
                 tauaer1,tauaer2,tauaer3,tauaer4,                       &
                 gaer1,gaer2,gaer3,gaer4,                               &
                 waer1,waer2,waer3,waer4,                               &
                 bscoef1,bscoef2,bscoef3,bscoef4,                       &
                 l2aer,l3aer,l4aer,l5aer,l6aer,l7aer,                   &
                 ids,ide, jds,jde, kds,kde,                             &
                 ims,ime, jms,jme, kms,kme,                             &
                 its,ite, jts,jte, kts,kte                              )
!----------------------------------------------------------------------------------
   USE module_configure
   USE module_state_description
   USE module_model_constants
!  USE module_data_mosaic_therm, only: nbin_a, nbin_a_maxd      
!
   INTEGER,      INTENT(IN   ) :: id,                                  &
                                  ids,ide, jds,jde, kds,kde,           &
                                  ims,ime, jms,jme, kms,kme,           &
                                  its,ite, jts,jte, kts,kte
   INTEGER,      INTENT(IN   ) :: ktau, nbin_o
      REAL,      INTENT(IN   ) :: dtstep
!
! array that holds all advected chemical species
!
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),             &
         INTENT(IN ) ::  chem
!
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                       &
         INTENT(IN ) ::  dz8w, alt, h2oai, h2oaj
!
! arrays that hold the aerosol optical properties
!
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                       &
         INTENT(INOUT ) ::                                             &
           tauaer1, tauaer2, tauaer3, tauaer4,                         &
           gaer1, gaer2, gaer3, gaer4,                                 &
           waer1, waer2, waer3, waer4,                                 &
           bscoef1, bscoef2, bscoef3, bscoef4
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:4 ),                  &
         INTENT(INOUT ) ::                                             &
           l2aer, l3aer, l4aer, l5aer, l6aer, l7aer
!
   TYPE(grid_config_rec_type),  INTENT(IN   )    :: config_flags
   LOGICAL, INTENT(IN) :: haveaer
!
! local variables
!
   real, dimension( its:ite, kts:kte, jts:jte, 1:nbin_o ) ::           &
           radius_wet, number_bin, radius_core
   real, dimension( 1:nbin_o, kts:kte) ::                              &
           radius_wet_col, number_bin_col, radius_core_col
   complex, dimension( its:ite, kts:kte, jts:jte, 1:nbin_o ) ::        &
           refindx, refindx_core, refindx_shell
   complex, dimension( 1:nbin_o, kts:kte) ::                           &
           refindx_col, refindx_core_col, refindx_shell_col
   real, dimension( kts:kte ) :: dz
!
   integer isecfrm0, nspint
   integer iclm, jclm, k, isize
   integer option_method, option_mie
   parameter ( nspint = 4 ) ! number of spectral interval bands 
   real, dimension( nspint, kts:kte ) ::                               &
           sizeaer,extaer,waer,gaer,tauaer,bscoef
   real, dimension( nspint, kts:kte ) ::                               &
           l2, l3, l4, l5, l6, l7
   real refr
   real fv
   complex aa, bb
!  save :: sizeaer,extaer,waer,gaer,tauaer,bscoef
!  save :: l2,l3,l4,l5,l6,l7
!----------------------------------------------------------------------------------
!
   chem_select: SELECT CASE(config_flags%chem_opt)
!
   CASE (RADM2SORG,RADM2SORG_KPP)
     call optical_prep_modal(nbin_o, chem, alt,                        &
          h2oai, h2oaj, refindx, radius_wet, number_bin,               &
          radius_core, refindx_core, refindx_shell,                    &
          ids,ide, jds,jde, kds,kde,                                   &
          ims,ime, jms,jme, kms,kme,                                   &
          its,ite, jts,jte, kts,kte                                    )
!    call mieaer_modal()
!
   CASE (CBMZ_MOSAIC_4BIN, CBMZ_MOSAIC_8BIN,                           &
         CBMZ_MOSAIC_4BIN_AQ, CBMZ_MOSAIC_8BIN_AQ)
     call optical_prep_sectional(nbin_o, chem, alt,                    &
          refindx, radius_wet, number_bin,                             &
          radius_core, refindx_core, refindx_shell,                    &
          ids,ide, jds,jde, kds,kde,                                   &
          ims,ime, jms,jme, kms,kme,                                   &
          its,ite, jts,jte, kts,kte                                    )
   END SELECT chem_select
     do jclm = jts, jte
     do iclm = its, ite
       do k = kts, kte
          dz(k) = dz8w(iclm, k, jclm)   ! cell depth (m)
       end do
       do k = kts, kte
       do isize = 1, nbin_o
          number_bin_col(isize,k) = number_bin(iclm,k,jclm,isize)
          radius_wet_col(isize,k) = radius_wet(iclm,k,jclm,isize)
          refindx_col(isize,k)    = refindx(iclm,k,jclm,isize)
          refr=real(refindx_col(isize,k))
          radius_core_col(isize,k)   = radius_core(iclm,k,jclm,isize)
          refindx_core_col(isize,k)  = refindx_core(iclm,k,jclm,isize)
          refindx_shell_col(isize,k) = refindx_shell(iclm,k,jclm,isize)

! JCB, Feb. 20, 2008:  in the case of shell/core and the use of the Mie 
! routine, set the refractive index of the shell used in the printout
! equal to the actual refractive index of the shell
          if(option_method.eq.3.and.option_mie.eq.2) &
               refindx_col(isize,k) = refindx_shell(iclm,k,jclm,isize) ! JCB

! JCB, Feb. 20, 2008:  set core radius = 0 for very small cores; this
! prevents problems with full-blown Mie calculations that do not deal
! well with very small cores.  For very small cores, the amount of
! absorption is negligible, and therefore setting the core radius to zero
! has virtually no effect on calculated optical properties
          if(radius_wet_col(isize,k) < 1e-20) then
               radius_core_col(isize,k)=0.0             
          else if(radius_core_col(isize,k)/radius_wet_col(isize,k)**3.le.0.0001) then
               radius_core_col(isize,k)=0.0  ! JCB
          end if
       enddo
       enddo
       isecfrm0 = ifix ( (ktau-1) * dtstep )

       if (option_method .eq. 2) then
          do k = kts, kte
          do isize = 1, nbin_o
           fv = (radius_core_col(isize,k)/radius_wet_col(isize,k))**3 ! volume fraction
           aa=(refindx_core_col(isize,k)**2+2.0*refindx_shell(iclm,k,jclm,isize)**2)
           bb=fv*(refindx_core_col(isize,k)**2-refindx_shell(iclm,k,jclm,isize)**2)
           refindx_col(isize,k)= refindx_shell(iclm,k,jclm,isize)*sqrt((aa+2.0*bb)/(aa-bb))
           refr=real(refindx_col(isize,k))
          enddo
          enddo
       endif

       if (option_method .le. 2) then
          do k = kts, kte
          do isize = 1, nbin_o
             radius_core_col(isize,k) = 0.0
             refindx_core_col(isize,k) = cmplx(0.0,0.0)
           enddo
          enddo
      endif
!
!!$       if(id.eq.1.and.iclm.eq.84.and.jclm.eq.52) then
!!$         print*,'jdf printout 1'
!!$         do isize = 1, nbin_o
!!$            write(*,888) isize,number_bin_col(isize,1),                &
!!$                   radius_wet_col(isize,1),radius_core_col(isize,1),   &
!!$                   real(refindx_col(isize,1)),                         &
!!$                   imag(refindx_col(isize,1)),                         &
!!$                   real(refindx_core_col(isize,1)),                    &
!!$                   imag(refindx_core_col(isize,1)),dz(1)
!!$         enddo
!!$       endif
!!$       if(id.eq.2.and.iclm.eq.59.and.jclm.eq.63) then
!!$         print*,'jdf printout 2'
!!$         do isize = 1, nbin_o
!!$            write(*,888) isize,number_bin_col(isize,1),                &
!!$                   radius_wet_col(isize,1),radius_core_col(isize,1),   &
!!$                   real(refindx_col(isize,1)),                         &
!!$                   imag(refindx_col(isize,1)),                         &
!!$                   real(refindx_core_col(isize,1)),                    &
!!$                   imag(refindx_core_col(isize,1)),dz(1)
!!$         enddo
!!$       endif
!!$  888  format(i3,9e12.5)
!
       if (option_mie .eq. 1) then 
          call mieaer(id, iclm, jclm, nbin_o,                          &
             number_bin_col, radius_wet_col, refindx_col,              &
             dz, isecfrm0, kte,                                        &
             sizeaer, extaer, waer, gaer, tauaer,                      &
             l2, l3, l4, l5, l6, l7, bscoef                            )
       endif
       if (option_mie .ge. 2 .and. option_method .le. 2) then 
          call mieaer_sc(id, iclm, jclm, nbin_o,                       &
             number_bin_col, radius_wet_col, refindx_col,              &
             radius_core_col, refindx_core_col,                        &
             dz, isecfrm0, kte,                                        &
             sizeaer, extaer, waer, gaer, tauaer,                      &
             l2, l3, l4, l5, l6, l7, bscoef                            )
       endif
       if (option_mie .ge. 2 .and. option_method .eq. 3) then 
          call mieaer_sc(id, iclm, jclm, nbin_o,                       &
             number_bin_col, radius_wet_col, refindx_shell_col,        &
             radius_core_col, refindx_core_col,                        &
             dz, isecfrm0, kte,                                        &
             sizeaer, extaer, waer, gaer, tauaer,                      &
             l2, l3, l4, l5, l6, l7, bscoef                            )
       endif
!
       do k=kts,kte
         tauaer1(iclm,k,jclm) = tauaer(1,k)
         tauaer2(iclm,k,jclm) = tauaer(2,k)
         tauaer3(iclm,k,jclm) = tauaer(3,k)
         tauaer4(iclm,k,jclm) = tauaer(4,k)
         gaer1(iclm,k,jclm)   = gaer(1,k)
         gaer2(iclm,k,jclm)   = gaer(2,k)
         gaer3(iclm,k,jclm)   = gaer(3,k)
         gaer4(iclm,k,jclm)   = gaer(4,k)
         waer1(iclm,k,jclm)   = waer(1,k)
         waer2(iclm,k,jclm)   = waer(2,k)
         waer3(iclm,k,jclm)   = waer(3,k)
         waer4(iclm,k,jclm)   = waer(4,k)
         bscoef1(iclm,k,jclm) = bscoef(1,k)
         bscoef2(iclm,k,jclm) = bscoef(2,k)
         bscoef3(iclm,k,jclm) = bscoef(3,k)
         bscoef4(iclm,k,jclm) = bscoef(4,k)
         l2aer(iclm,k,jclm,1) = l2(1,k)
         l2aer(iclm,k,jclm,2) = l2(2,k)
         l2aer(iclm,k,jclm,3) = l2(3,k)
         l2aer(iclm,k,jclm,4) = l2(4,k)
         l3aer(iclm,k,jclm,1) = l3(1,k)
         l3aer(iclm,k,jclm,2) = l3(2,k)
         l3aer(iclm,k,jclm,3) = l3(3,k)
         l3aer(iclm,k,jclm,4) = l3(4,k)
         l4aer(iclm,k,jclm,1) = l4(1,k)
         l4aer(iclm,k,jclm,2) = l4(2,k)
         l4aer(iclm,k,jclm,3) = l4(3,k)
         l4aer(iclm,k,jclm,4) = l4(4,k)
         l5aer(iclm,k,jclm,1) = l5(1,k)
         l5aer(iclm,k,jclm,2) = l5(2,k)
         l5aer(iclm,k,jclm,3) = l5(3,k)
         l5aer(iclm,k,jclm,4) = l5(4,k)
         l6aer(iclm,k,jclm,1) = l6(1,k)
         l6aer(iclm,k,jclm,2) = l6(2,k)
         l6aer(iclm,k,jclm,3) = l6(3,k)
         l6aer(iclm,k,jclm,4) = l6(4,k)
         l7aer(iclm,k,jclm,1) = l7(1,k)
         l7aer(iclm,k,jclm,2) = l7(2,k)
         l7aer(iclm,k,jclm,3) = l7(3,k)
         l7aer(iclm,k,jclm,4) = l7(4,k)
       enddo
!!$       if(id.eq.1.and.iclm.eq.84.and.jclm.eq.52) then
!!$          write(*,889) sizeaer(1,1),sizeaer(2,1),sizeaer(3,1),sizeaer(4,1)
!!$          write(*,889) extaer(1,1),extaer(2,1),extaer(3,1),extaer(4,1)
!!$          write(*,889) waer(1,1),waer(2,1),waer(3,1),waer(4,1)
!!$          write(*,889) gaer(1,1),gaer(2,1),gaer(3,1),gaer(4,1)
!!$          write(*,889) tauaer(1,1),tauaer(2,1),tauaer(3,1),tauaer(4,1)
!!$          write(*,889) bscoef(1,1),bscoef(2,1),bscoef(3,1),bscoef(4,1)
!!$          write(*,889) l2(1,1),l2(2,1),l2(3,1),l2(4,1)
!!$          write(*,889) l3(1,1),l3(2,1),l3(3,1),l3(4,1)
!!$          write(*,889) l4(1,1),l4(2,1),l4(3,1),l4(4,1)
!!$          write(*,889) l5(1,1),l5(2,1),l5(3,1),l5(4,1)
!!$          write(*,889) l6(1,1),l6(2,1),l6(3,1),l6(4,1)
!!$          write(*,889) l7(1,1),l7(2,1),l7(3,1),l7(4,1)
!!$       endif
!!$       if(id.eq.2.and.iclm.eq.59.and.jclm.eq.63) then
!!$          write(*,889) sizeaer(1,1),sizeaer(2,1),sizeaer(3,1),sizeaer(4,1)
!!$          write(*,889) extaer(1,1),extaer(2,1),extaer(3,1),extaer(4,1)
!!$          write(*,889) waer(1,1),waer(2,1),waer(3,1),waer(4,1)
!!$          write(*,889) gaer(1,1),gaer(2,1),gaer(3,1),gaer(4,1)
!!$          write(*,889) tauaer(1,1),tauaer(2,1),tauaer(3,1),tauaer(4,1)
!!$          write(*,889) bscoef(1,1),bscoef(2,1),bscoef(3,1),bscoef(4,1)
!!$          write(*,889) l2(1,1),l2(2,1),l2(3,1),l2(4,1)
!!$          write(*,889) l3(1,1),l3(2,1),l3(3,1),l3(4,1)
!!$          write(*,889) l4(1,1),l4(2,1),l4(3,1),l4(4,1)
!!$          write(*,889) l5(1,1),l5(2,1),l5(3,1),l5(4,1)
!!$          write(*,889) l6(1,1),l6(2,1),l6(3,1),l6(4,1)
!!$          write(*,889) l7(1,1),l7(2,1),l7(3,1),l7(4,1)
!!$       endif
!!$  889  format(4e12.5)
     enddo
     enddo
!
      return
!
      end subroutine optical_averaging
!----------------------------------------------------------------------------------
! This subroutine computes volume-averaged refractive index and wet radius needed 
! by the mie calculations. Aerosol number is also passed into the mie calculations
! in terms of other units.
!
      subroutine optical_prep_sectional(nbin_o, chem, alt,       &
        refindx, radius_wet, number_bin,                         &
        radius_core, refindx_core,  refindx_shell,               &
        ids,ide, jds,jde, kds,kde,                               &
        ims,ime, jms,jme, kms,kme,                               &
        its,ite, jts,jte, kts,kte                                )
!
   USE module_configure
!  USE module_state_description
   USE module_model_constants
   USE module_data_mosaic_asect
   USE module_data_mosaic_other
   USE module_state_description, only:  param_first_scalar
!
   INTEGER, INTENT(IN   ) :: its,ite, jts,jte, kts,kte, nbin_o
   INTEGER, INTENT(IN   ) :: ims,ime, jms,jme, kms,kme
   INTEGER, INTENT(IN   ) :: ids,ide, jds,jde, kds,kde

   REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),             &
         INTENT(IN ) ::  chem
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                       &
         INTENT(IN ) ::  alt
   REAL, DIMENSION( its:ite, kts:kte, jts:jte, 1:nbin_o),              &
         INTENT(OUT ) ::                                               &
           radius_wet, number_bin, radius_core
   COMPLEX, DIMENSION( its:ite, kts:kte, jts:jte, 1:nbin_o),           &
         INTENT(OUT ) ::                                               &
           refindx, refindx_core, refindx_shell
!
! local variables
!
   integer i, j, k, l, isize, itype, iphase
   integer p1st
   complex  ref_index_nh4so4 , ref_index_lvcite  , ref_index_nh4hso4, &
            ref_index_nh4msa , ref_index_nh4no3  , ref_index_nh4cl , &
            ref_index_nacl   , ref_index_nano3   , ref_index_na2so4, &
            ref_index_na3hso4, ref_index_nahso4  , ref_index_namsa,  &
            ref_index_caso4  , ref_index_camsa2  , ref_index_cano3,  &
            ref_index_cacl2  , ref_index_caco3   , ref_index_h2so4,  &
            ref_index_hhso4  , ref_index_hno3    , ref_index_hcl,    &
            ref_index_msa    , ref_index_oc      , ref_index_bc,     &
            ref_index_oin    , ref_index_aro1    , ref_index_aro2,   &
            ref_index_alk1   , ref_index_ole1    , ref_index_api1,   &
            ref_index_api2   , ref_index_im1     , ref_index_im2,    &
            ref_index_h2o    , ri_dum            , ri_ave_a
   real  dens_so4  , dens_no3  , dens_cl   , dens_msa  , dens_co3 ,  &
         dens_nh4  , dens_na   , dens_ca   , dens_oin  , dens_oc  ,  &
         dens_bc   , dens_aro1 , dens_aro2 , dens_alk1 , dens_ole1,  &
         dens_api1 , dens_api2 , dens_lim1 , dens_lim2 , dens_h2o
   real  mass_so4  , mass_no3  , mass_cl   , mass_msa  , mass_co3 ,  &
         mass_nh4  , mass_na   , mass_ca   , mass_oin  , mass_oc  ,  &
         mass_bc   , mass_aro1 , mass_aro2 , mass_alk1 , mass_ole1,  &
         mass_api1 , mass_api2 , mass_lim1 , mass_lim2 , mass_h2o
   real  vol_so4   , vol_no3   , vol_cl    , vol_msa   , vol_co3  ,  &
         vol_nh4   , vol_na    , vol_ca    , vol_oin   , vol_oc   ,  &
         vol_bc    , vol_aro1  , vol_aro2  , vol_alk1  , vol_ole1 ,  &
         vol_api1  , vol_api2  , vol_lim1  , vol_lim2  , vol_h2o
   real  conv1a, conv1b
   real  mass_dry_a, mass_wet_a, vol_dry_a , vol_wet_a , vol_shell,  &
         dp_dry_a  , dp_wet_a  , num_a     , dp_bc_a
   real  num_a_lo  , num_a_hi
   real  refr
!
! Define refractive indicies
! * assume na and cl are the same as nacl
! * assume so4, no3, and nh4 are the same as nh4no3
! * assume ca and co3 are the same as caco3
! * assume msa is just msa
! Further work:
! * to be more precise, need to compute electrolytes to apportion
!   so4, no3, nh4, na, cl, msa, ca, co3 among various componds
!   as was done previously in module_mosaic_therm.F
!
      ref_index_nh4so4 = cmplx(1.52,0.)
      ref_index_lvcite = cmplx(1.50,0.)
      ref_index_nh4hso4= cmplx(1.47,0.)
      ref_index_nh4msa = cmplx(1.50,0.)     ! assumed
      ref_index_nh4no3 = cmplx(1.50,0.)
      ref_index_nh4cl  = cmplx(1.50,0.)
      ref_index_nacl   = cmplx(1.45,0.)
      ref_index_nano3  = cmplx(1.50,0.)
      ref_index_na2so4 = cmplx(1.50,0.)
      ref_index_na3hso4= cmplx(1.50,0.)
      ref_index_nahso4 = cmplx(1.50,0.)
      ref_index_namsa  = cmplx(1.50,0.)     ! assumed
      ref_index_caso4  = cmplx(1.56,0.006)
      ref_index_camsa2 = cmplx(1.56,0.006)  ! assumed
      ref_index_cano3  = cmplx(1.56,0.006)
      ref_index_cacl2  = cmplx(1.52,0.006)
      ref_index_caco3  = cmplx(1.68,0.006)
      ref_index_h2so4  = cmplx(1.43,0.)
      ref_index_hhso4  = cmplx(1.43,0.)
      ref_index_hno3   = cmplx(1.50,0.)
      ref_index_hcl    = cmplx(1.50,0.)
      ref_index_msa    = cmplx(1.43,0.)     ! assumed
      ref_index_oc     = cmplx(1.45,0.)  ! JCB, Feb. 20, 2008: no complex part?
! JCB, Feb. 20, 2008:  set the refractive index of BC equal to the midpoint
! of ranges given in Bond and Bergstrom, Light absorption by carboneceous
! particles: an investigative review 2006, Aerosol Sci. and Tech., 40:27-67.
!      ref_index_bc     = cmplx(1.82,0.74), old value
      ref_index_bc     = cmplx(1.85,0.71)
      ref_index_oin    = cmplx(1.55,0.006)  ! JCB, Feb. 20, 2008: if  "other inorganics" includes dust, then this refractive index should be wavelength depedent
      ref_index_aro1   = cmplx(1.45,0.)
      ref_index_aro2   = cmplx(1.45,0.)
      ref_index_alk1   = cmplx(1.45,0.)
      ref_index_ole1   = cmplx(1.45,0.)
      ref_index_api1   = cmplx(1.45,0.)
      ref_index_api2   = cmplx(1.45,0.)
      ref_index_im1   = cmplx(1.45,0.)
      ref_index_im2   = cmplx(1.45,0.)
      ref_index_h2o    = cmplx(1.33,0.)
!
! densities in g/cc
!
      dens_so4   = 1.8        ! used
      dens_no3   = 1.8        ! used
      dens_cl    = 2.2        ! used
      dens_msa   = 1.8        ! used
      dens_co3   = 2.6        ! used
      dens_nh4   = 1.8        ! used
      dens_na    = 2.2        ! used
      dens_ca    = 2.6        ! used
      dens_oin   = 2.6        ! used
      dens_oc    = 1.0        ! used
! JCB, Feb. 20, 2008:  the density of BC is updated to reflect values
! published by Bond and Bergstrom, Light absorption by carboneceous
! particles: an investigative review 2006, Aerosol Sci. and Tech., 40:27-67.
!      dens_bc    = 1.7        ! used, old value
      dens_bc    =  1.8        !  midpoint of Bond and Bergstrom value
      dens_aro1  = 1.0
      dens_aro2  = 1.0
      dens_alk1  = 1.0
      dens_ole1  = 1.0
      dens_api1  = 1.0
      dens_api2  = 1.0
      dens_lim1  = 1.0
      dens_lim2  = 1.0
      dens_h2o   = 1.0
!
      p1st = param_first_scalar
!
      do isize = 1, nbin_o
      do j = jts, jte
      do k = kts, kte
      do i = its, ite
        refindx(i,k,j,isize)=0.0
        radius_wet(i,k,j,isize)=0.0
        number_bin(i,k,j,isize)=0.0
        radius_core(i,k,j,isize)=0.0
        refindx_core(i,k,j,isize)=0.0
        refindx_shell(i,k,j,isize)=0.0
      enddo
      enddo
      enddo
      enddo
!
! units:
! * mass     - g/cc(air)
! * number   - #/cc(air)
! * volume   - cc(air)/cc(air)
! * diameter - cm
!
      itype=1
      iphase=1
      do j = jts, jte
      do k = kts, kte
      do i = its, ite
        do isize = 1, nbin_o
          mass_so4 = 0.0
          mass_no3 = 0.0
          mass_nh4 = 0.0
          mass_oin = 0.0
          mass_oc = 0.0
          mass_bc = 0.0
          mass_na = 0.0
          mass_cl = 0.0
          mass_msa = 0.0
          mass_co3 = 0.0
          mass_ca = 0.0
          mass_h2o = 0.0
! convert ug / kg dry air to g / cc  air
          conv1a = (1.0/alt(i,k,j)) * 1.0e-12
! convert # / kg dry air to # / cc  air
          conv1b = (1.0/alt(i,k,j)) * 1.0e-6
          l=lptr_so4_aer(isize,itype,iphase)
          if (l .ge. p1st)  mass_so4= chem(i,k,j,l)*conv1a
          l=lptr_no3_aer(isize,itype,iphase)
          if (l .ge. p1st)  mass_no3= chem(i,k,j,l)*conv1a
          l=lptr_nh4_aer(isize,itype,iphase)
          if (l .ge. p1st)  mass_nh4= chem(i,k,j,l)*conv1a
          l=lptr_oin_aer(isize,itype,iphase)
          if (l .ge. p1st)  mass_oin= chem(i,k,j,l)*conv1a
          l=lptr_oc_aer(isize,itype,iphase)
          if (l .ge. p1st)  mass_oc= chem(i,k,j,l)*conv1a
          l=lptr_bc_aer(isize,itype,iphase)
          if (l .ge. p1st)  mass_bc= chem(i,k,j,l)*conv1a
          l=lptr_na_aer(isize,itype,iphase)
          if (l .ge. p1st)  mass_na= chem(i,k,j,l)*conv1a
          l=lptr_cl_aer(isize,itype,iphase)
          if (l .ge. p1st)  mass_cl= chem(i,k,j,l)*conv1a
          l=lptr_msa_aer(isize,itype,iphase)
          if (l .ge. p1st)  mass_msa= chem(i,k,j,l)*conv1a
          l=lptr_co3_aer(isize,itype,iphase)
          if (l .ge. p1st)  mass_co3= chem(i,k,j,l)*conv1a
          l=lptr_ca_aer(isize,itype,iphase)
          if (l .ge. p1st)  mass_ca= chem(i,k,j,l)*conv1a
          l=waterptr_aer(isize,itype)
          if (l .ge. p1st)  mass_h2o= chem(i,k,j,l)*conv1a
          l=numptr_aer(isize,itype,iphase)
          if (l .ge. p1st)  num_a= chem(i,k,j,l)*conv1b
!
          vol_so4 = mass_so4 / dens_so4
          vol_no3 = mass_no3 / dens_no3
          vol_nh4 = mass_nh4 / dens_nh4
          vol_oin = mass_oin / dens_oin
          vol_oc  = mass_oc  / dens_oc
          vol_bc  = mass_bc  / dens_bc
          vol_na  = mass_na  / dens_na
          vol_cl  = mass_cl  / dens_cl
          vol_msa = mass_msa / dens_msa
          vol_co3 = mass_co3 / dens_co3
          vol_ca  = mass_ca  / dens_ca
!         mass_h2o = 0.0 ! testing purposes only
          vol_h2o = mass_h2o / dens_h2o
          mass_dry_a = mass_so4 + mass_no3 + mass_nh4 + mass_oin + &
                       mass_oc  + mass_bc  + mass_na  + mass_cl  + &
                       mass_msa + mass_co3 + mass_ca
          mass_wet_a = mass_dry_a + mass_h2o 
          vol_dry_a  = vol_so4  + vol_no3  + vol_nh4  + vol_oin  + &
                       vol_oc   + vol_bc   + vol_na   + vol_cl   + &
                       vol_msa  + vol_co3  + vol_ca
          vol_wet_a  = vol_dry_a + vol_h2o
          vol_shell  = vol_wet_a - vol_bc
! jdf: Adjustment of aerosol number if it falls outside of reasonable bounds.
! This is necessary since advection scheme will cause the mass and number
! prognostic equations to give inconsistent values, usually at sharp gradients
! of these quantities.
          num_a_lo=1.90985*vol_dry_a/(dlo_sect(isize,itype)**3)
          num_a_hi=1.90985*vol_dry_a/(dhi_sect(isize,itype)**3)
          if(num_a.gt.num_a_lo) then
!           if(i.eq.84.and.j.eq.52) then
!             print*,'jdf adjust lower',k,isize,vol_dry_a,num_a,num_a_lo,num_a_hi
!           endif
            num_a=num_a_lo
          elseif(num_a.lt.num_a_hi) then
!           if(i.eq.84.and.j.eq.52) then
!             print*,'jdf adjust higher',k,isize,vol_dry_a,num_a,num_a_lo,num_a_hi
!           endif
            num_a=num_a_hi
          endif
!
          dp_dry_a   = (1.90985*vol_dry_a/num_a)**0.3333333
          dp_wet_a   = (1.90985*vol_wet_a/num_a)**0.3333333
          dp_bc_a    = (1.90985*vol_bc/num_a)**0.3333333
          ri_dum     = (0.0,0.0)
          ri_dum     = (ref_index_nh4so4 * mass_so4 / dens_so4) +  &
                       (ref_index_nh4no3 * mass_no3 / dens_no3) +  &
                       (ref_index_nh4no3 * mass_nh4 / dens_nh4) +  &
                       (ref_index_oin    * mass_oin / dens_oin) +  &
                       (ref_index_oc     * mass_oc  / dens_oc) +   &
                       (ref_index_bc     * mass_bc  / dens_bc) +   &
                       (ref_index_nacl   * mass_na  / dens_na) +   &
                       (ref_index_nacl   * mass_cl  / dens_cl) +   &
                       (ref_index_msa    * mass_msa / dens_msa) +  &
                       (ref_index_caco3  * mass_ca  / dens_ca) +   &
                       (ref_index_caco3  * mass_co3 / dens_co3) +  &
                       (ref_index_h2o    * mass_h2o / dens_h2o) 
          ri_ave_a   = ri_dum/vol_wet_a
          ri_dum     = (ref_index_nh4so4 * mass_so4 / dens_so4) +  &
                       (ref_index_nh4no3 * mass_no3 / dens_no3) +  &
                       (ref_index_nh4no3 * mass_nh4 / dens_nh4) +  &
                       (ref_index_oin    * mass_oin / dens_oin) +  &
                       (ref_index_oc     * mass_oc  / dens_oc) +   &
                       (ref_index_nacl   * mass_na  / dens_na) +   &
                       (ref_index_nacl   * mass_cl  / dens_cl) +   &
                       (ref_index_msa    * mass_msa / dens_msa) +  &
                       (ref_index_caco3  * mass_ca  / dens_ca) +   &
                       (ref_index_caco3  * mass_co3 / dens_co3) +  &
                       (ref_index_h2o    * mass_h2o / dens_h2o) 
          if(dp_wet_a/2.0 .lt. dlo_sect(isize,itype)/2.0) then
            refindx(i,k,j,isize)    = (1.5,0.0)
            radius_wet(i,k,j,isize) =dlo_sect(isize,itype)/2.0
            number_bin(i,k,j,isize) =num_a
            radius_core(i,k,j,isize) =0.0
            refindx_core(i,k,j,isize) = ref_index_bc
            refindx_shell(i,k,j,isize) = ref_index_oin
          elseif(vol_shell .lt. 1.0e-20) then
            refindx(i,k,j,isize)    = (1.5,0.0)
            radius_wet(i,k,j,isize) =dlo_sect(isize,itype)/2.0
            number_bin(i,k,j,isize) =num_a
            radius_core(i,k,j,isize) =0.0
            refindx_core(i,k,j,isize) = ref_index_bc
            refindx_shell(i,k,j,isize) = ref_index_oin
          else
            refindx(i,k,j,isize)    =ri_ave_a
            radius_wet(i,k,j,isize) =dp_wet_a/2.0
            number_bin(i,k,j,isize) =num_a
            radius_core(i,k,j,isize) =dp_bc_a/2.0
            refindx_core(i,k,j,isize) =ref_index_bc
            refindx_shell(i,k,j,isize) =ri_dum/vol_shell
          endif
!
          refr=real(refindx(i,k,j,isize))
        enddo
        enddo
        enddo
      enddo

      return

      end subroutine optical_prep_sectional
!
!----------------------------------------------------------------------------------
! This subroutine computes volume-averaged refractive index and wet radius needed 
! by the mie calculations. Aerosol number is also passed into the mie calculations
! in terms of other units.
!
      subroutine optical_prep_modal(nbin_o, chem, alt,           &
        h2oai, h2oaj, refindx, radius_wet, number_bin,           &
        radius_core, refindx_core, refindx_shell,                &
        ids,ide, jds,jde, kds,kde,                               &
        ims,ime, jms,jme, kms,kme,                               &
        its,ite, jts,jte, kts,kte                                )
!
   USE module_configure
!  USE module_state_description
   USE module_model_constants
   USE module_state_description, only:  param_first_scalar
   USE module_data_sorgam
!
   INTEGER, INTENT(IN   ) :: its,ite, jts,jte, kts,kte, nbin_o
   INTEGER, INTENT(IN   ) :: ims,ime, jms,jme, kms,kme
   INTEGER, INTENT(IN   ) :: ids,ide, jds,jde, kds,kde

   REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ),             &
         INTENT(IN ) ::  chem
   REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),                       &
         INTENT(IN ) ::  alt, h2oai, h2oaj
   REAL, DIMENSION( its:ite, kts:kte, jts:jte, 1:nbin_o),              &
         INTENT(OUT ) ::                                               &
           radius_wet, number_bin, radius_core
   COMPLEX, DIMENSION( its:ite, kts:kte, jts:jte, 1:nbin_o),           &
         INTENT(OUT ) ::                                               &
           refindx, refindx_core, refindx_shell
!
! local variables
!
   integer i, j, k, l, isize, itype, iphase
   integer p1st
   complex  ref_index_nh4so4 , ref_index_lvcite  , ref_index_nh4hso4, &
            ref_index_nh4msa , ref_index_nh4no3  , ref_index_nh4cl , &
            ref_index_nacl   , ref_index_nano3   , ref_index_na2so4, &
            ref_index_na3hso4, ref_index_nahso4  , ref_index_namsa,  &
            ref_index_caso4  , ref_index_camsa2  , ref_index_cano3,  &
            ref_index_cacl2  , ref_index_caco3   , ref_index_h2so4,  &
            ref_index_hhso4  , ref_index_hno3    , ref_index_hcl,    &
            ref_index_msa    , ref_index_oc      , ref_index_bc,     &
            ref_index_oin    , ref_index_aro1    , ref_index_aro2,   &
            ref_index_alk1   , ref_index_ole1    , ref_index_api1,   &
            ref_index_api2   , ref_index_im1     , ref_index_im2,    &
            ref_index_h2o    , ri_dum            , ri_ave_a
   real  dens_so4  , dens_no3  , dens_cl   , dens_msa  , dens_co3 ,  &
         dens_nh4  , dens_na   , dens_ca   , dens_oin  , dens_oc  ,  &
         dens_bc   , dens_aro1 , dens_aro2 , dens_alk1 , dens_ole1,  &
         dens_api1 , dens_api2 , dens_lim1 , dens_lim2 , dens_h2o
   real  mass_so4  , mass_no3  , mass_cl   , mass_msa  , mass_co3 ,  &
         mass_nh4  , mass_na   , mass_ca   , mass_oin  , mass_oc  ,  &
         mass_bc   , mass_aro1 , mass_aro2 , mass_alk1 , mass_ole1,  &
         mass_api1 , mass_api2 , mass_lim1 , mass_lim2 , mass_h2o
   real  mass_so4i , mass_no3i , mass_cli  , mass_msai , mass_co3i,  &
         mass_nh4i , mass_nai  , mass_cai  , mass_oini , mass_oci ,  &
         mass_bci  , mass_aro1i, mass_aro2i, mass_alk1i, mass_ole1i, &
         mass_ba1i , mass_ba2i,  mass_ba3i , mass_ba4i , mass_pai,   &
         mass_h2oi
   real  mass_so4j , mass_no3j , mass_clj  , mass_msaj , mass_co3j,  &
         mass_nh4j , mass_naj  , mass_caj  , mass_oinj , mass_ocj ,  &
         mass_bcj  , mass_aro1j, mass_aro2j, mass_alk1j, mass_ole1j, &
         mass_ba1j , mass_ba2j,  mass_ba3j , mass_ba4j , mass_paj,   &
         mass_h2oj
   real  mass_antha, mass_seas, mass_soil
   real  num_ai, num_aj, num_ac,  vol_ai, vol_aj, vol_ac
   real  vol_so4   , vol_no3   , vol_cl    , vol_msa   , vol_co3  ,  &
         vol_nh4   , vol_na    , vol_ca    , vol_oin   , vol_oc   ,  &
         vol_bc    , vol_aro1  , vol_aro2  , vol_alk1  , vol_ole1 ,  &
         vol_api1  , vol_api2  , vol_lim1  , vol_lim2  , vol_h2o
   real  conv1a, conv1b
   real  mass_dry_a, mass_wet_a, vol_dry_a , vol_wet_a , vol_shell,  &
         dp_dry_a  , dp_wet_a  , num_a     , dp_bc_a
   real  ifac, jfac, cfac
   real  refr
   real  dgnum_um,drydens,duma,dlo_um,dhi_um,dgmin,sixpi,ss1,ss2,ss3,dtemp
   integer  iflag
   real, dimension(1:nbin_o) :: xnum_secti,xnum_sectj,xnum_sectc
   real, dimension(1:nbin_o) :: xmas_secti,xmas_sectj,xmas_sectc
   real, dimension(1:nbin_o) :: xdia_um, xdia_cm
!
!  real  sginin,sginia,sginic from module_data_sorgam.F
! 
! Mass from modal distribution is divided into individual sections before
! being passed back into the Mie routine.
! * currently use the same size bins as 8 default MOSAIC size bins
! * dlo_um and dhi_um define the lower and upper bounds of individual sections
!   used to compute optical properties
! * sigmas for 3 modes taken from module_sorgan_data.F
! * these parameters are needed by sect02 that is called later
! * sginin=1.7, sginia=2.0, sginic=2.5
!
      sixpi=6.0/3.14159265359
      dlo_um=0.0390625
      dhi_um=10.0
      drydens=1.8
      iflag=2
      duma=1.0
      dgmin=1.0e-07 ! in (cm)
      dtemp=dlo_um
      do isize=1,nbin_o
        xdia_um(isize)=(dtemp+dtemp*2.0)/2.0
        dtemp=dtemp*2.0
      enddo
!
! Define refractive indicies
! * assume na and cl are the same as nacl
! * assume so4, no3, and nh4 are the same as nh4no3
! * assume ca and co3 are the same as caco3
! * assume msa is just msa
! Further work:
! * to be more precise, need to compute electrolytes to apportion
!   so4, no3, nh4, na, cl, msa, ca, co3 among various componds
!   as was done previously in module_mosaic_therm.F
!
      ref_index_nh4so4 = cmplx(1.52,0.)
      ref_index_lvcite = cmplx(1.50,0.)
      ref_index_nh4hso4= cmplx(1.47,0.)
      ref_index_nh4msa = cmplx(1.50,0.)     ! assumed
      ref_index_nh4no3 = cmplx(1.50,0.)
      ref_index_nh4cl  = cmplx(1.50,0.)
      ref_index_nacl   = cmplx(1.45,0.)
      ref_index_nano3  = cmplx(1.50,0.)
      ref_index_na2so4 = cmplx(1.50,0.)
      ref_index_na3hso4= cmplx(1.50,0.)
      ref_index_nahso4 = cmplx(1.50,0.)
      ref_index_namsa  = cmplx(1.50,0.)     ! assumed
      ref_index_caso4  = cmplx(1.56,0.006)
      ref_index_camsa2 = cmplx(1.56,0.006)  ! assumed
      ref_index_cano3  = cmplx(1.56,0.006)
      ref_index_cacl2  = cmplx(1.52,0.006)
      ref_index_caco3  = cmplx(1.68,0.006)
      ref_index_h2so4  = cmplx(1.43,0.)
      ref_index_hhso4  = cmplx(1.43,0.)
      ref_index_hno3   = cmplx(1.50,0.)
      ref_index_hcl    = cmplx(1.50,0.)
      ref_index_msa    = cmplx(1.43,0.)     ! assumed
      ref_index_oc     = cmplx(1.45,0.)  ! JCB, Feb. 20, 2008: no complex part?
! JCB, Feb. 20, 2008:  set the refractive index of BC equal to the
! midpoint of ranges given in Bond and Bergstrom, Light absorption by
! carboneceous particles: an investigative review 2006, Aerosol Sci.
! and Tech., 40:27-67.
!     ref_index_bc     = cmplx(1.82,0.74) old value
      ref_index_bc     = cmplx(1.85,0.71)
      ref_index_oin    = cmplx(1.55,0.006)  ! JCB, Feb. 20, 2008: if  "other inorganics" includes dust, then this refractive index should be wavelength depedent
      ref_index_aro1   = cmplx(1.45,0.)
      ref_index_aro2   = cmplx(1.45,0.)
      ref_index_alk1   = cmplx(1.45,0.)
      ref_index_ole1   = cmplx(1.45,0.)
      ref_index_api1   = cmplx(1.45,0.)
      ref_index_api2   = cmplx(1.45,0.)
      ref_index_im1   = cmplx(1.45,0.)
      ref_index_im2   = cmplx(1.45,0.)
      ref_index_h2o    = cmplx(1.33,0.)
!
! densities in g/cc
!
      dens_so4   = 1.8        ! used
      dens_no3   = 1.8        ! used
      dens_cl    = 2.2        ! used
      dens_msa   = 1.8        ! used
      dens_co3   = 2.6        ! used
      dens_nh4   = 1.8        ! used
      dens_na    = 2.2        ! used
      dens_ca    = 2.6        ! used
      dens_oin   = 2.6        ! used
      dens_oc    = 1.0        ! used
! JCB, Feb. 20, 2008:  the density of BC is updated to reflect values
! published by Bond and Bergstrom, Light absorption by carboneceous
! particles: an investigative review 2006, Aerosol Sci. and Tech., 40:27-67.
!      dens_bc    = 1.7        ! used, old value
      dens_bc    =  1.8        !  midpoint of Bond and Bergstrom value
      dens_aro1  = 1.0
      dens_aro2  = 1.0
      dens_alk1  = 1.0
      dens_ole1  = 1.0
      dens_api1  = 1.0
      dens_api2  = 1.0
      dens_lim1  = 1.0
      dens_lim2  = 1.0
      dens_h2o   = 1.0
!
      p1st = param_first_scalar
!
      do isize = 1, nbin_o
      do j = jts, jte
      do k = kts, kte
      do i = its, ite
        refindx(i,k,j,isize)=0.0
        radius_wet(i,k,j,isize)=0.0
        number_bin(i,k,j,isize)=0.0
        radius_core(i,k,j,isize)=0.0
        refindx_core(i,k,j,isize)=0.0
        refindx_shell(i,k,j,isize)=0.0
      enddo
      enddo
      enddo
      enddo
!
! units:
! * mass     - g/cc(air)
! * number   - #/cc(air)
! * volume   - cc(air)/cc(air)
! * diameter - cm
!
      itype=1
      iphase=1
      do j = jts, jte
      do k = kts, kte
      do i = its, ite
        mass_so4i = 0.0
        mass_so4j = 0.0
        mass_no3i = 0.0
        mass_no3j = 0.0
        mass_nh4i = 0.0
        mass_nh4j = 0.0
        mass_oini = 0.0
        mass_oinj = 0.0
        mass_aro1i = 0.0
        mass_aro1j = 0.0
        mass_aro2i = 0.0
        mass_aro2j = 0.0
        mass_alk1i = 0.0
        mass_alk1j = 0.0
        mass_ole1i = 0.0
        mass_ole1j = 0.0
        mass_ba1i = 0.0
        mass_ba1j = 0.0
        mass_ba2i = 0.0
        mass_ba2j = 0.0
        mass_ba3i = 0.0
        mass_ba3j = 0.0
        mass_ba4i = 0.0
        mass_ba4j = 0.0
        mass_pai = 0.0
        mass_paj = 0.0
        mass_oci = 0.0
        mass_ocj = 0.0
        mass_bci = 0.0
        mass_bcj = 0.0
        mass_cai = 0.0
        mass_caj = 0.0
        mass_co3i = 0.0
        mass_co3j = 0.0
        mass_msai = 0.0
        mass_msaj = 0.0
        mass_h2oi = 0.0
        mass_h2oj = 0.0
        mass_antha = 0.0
        mass_seas = 0.0
        mass_soil = 0.0
        vol_aj = 0.0
        vol_ai = 0.0
        vol_ac = 0.0
        num_aj = 0.0
        num_ai = 0.0
        num_ac = 0.0

! convert ug / kg dry air to g / cc  air
        conv1a = (1.0/alt(i,k,j)) * 1.0e-12
! convert # / kg dry air to # / cc  air
        conv1b = (1.0/alt(i,k,j)) * 1.0e-6

! Accumulation mode...
        l=lptr_so4_aer(1,1,iphase)
        if (l .ge. p1st)  mass_so4j= chem(i,k,j,l)*conv1a
        l=lptr_no3_aer(1,1,iphase)
        if (l .ge. p1st)  mass_no3j= chem(i,k,j,l)*conv1a
        l=lptr_nh4_aer(1,1,iphase)
        if (l .ge. p1st)  mass_nh4j= chem(i,k,j,l)*conv1a
        l=lptr_p25_aer(1,1,iphase)
        if (l .ge. p1st)  mass_oinj= chem(i,k,j,l)*conv1a
        l=lptr_orgaro1_aer(1,1,iphase)
        if (l .ge. p1st)  mass_aro1j= chem(i,k,j,l)*conv1a
        l=lptr_orgaro2_aer(1,1,iphase)
        if (l .ge. p1st)  mass_aro2j= chem(i,k,j,l)*conv1a
        l=lptr_orgalk_aer(1,1,iphase)
        if (l .ge. p1st)  mass_alk1j= chem(i,k,j,l)*conv1a
        l=lptr_orgole_aer(1,1,iphase)
        if (l .ge. p1st)  mass_ole1j= chem(i,k,j,l)*conv1a
        l=lptr_orgba1_aer(1,1,iphase)
        if (l .ge. p1st)  mass_ba1j= chem(i,k,j,l)*conv1a
        l=lptr_orgba2_aer(1,1,iphase)
        if (l .ge. p1st)  mass_ba2j= chem(i,k,j,l)*conv1a
        l=lptr_orgba3_aer(1,1,iphase)
        if (l .ge. p1st)  mass_ba3j= chem(i,k,j,l)*conv1a
        l=lptr_orgba4_aer(1,1,iphase)
        if (l .ge. p1st)  mass_ba4j= chem(i,k,j,l)*conv1a
        l=lptr_orgpa_aer(1,1,iphase)
        if (l .ge. p1st)  mass_paj= chem(i,k,j,l)*conv1a
        l=lptr_ec_aer(1,1,iphase)
        if (l .ge. p1st)  mass_bcj= chem(i,k,j,l)*conv1a
        l=numptr_aer(1,1,iphase)
        if (l .ge. p1st)  num_aj= chem(i,k,j,l)*conv1b
        mass_h2oj= h2oaj(i,k,j) * 1.0e-12
        mass_ocj=mass_aro1j+mass_aro2j+mass_alk1j+mass_ole1j+ &
                 mass_ba1j+mass_ba2j+mass_ba3j+mass_ba4j+mass_paj

! Aitken mode...
        l=lptr_so4_aer(1,2,iphase)
        if (l .ge. p1st)  mass_so4i= chem(i,k,j,l)*conv1a
        l=lptr_no3_aer(1,2,iphase)
        if (l .ge. p1st)  mass_no3i= chem(i,k,j,l)*conv1a
        l=lptr_nh4_aer(1,2,iphase)
        if (l .ge. p1st)  mass_nh4i= chem(i,k,j,l)*conv1a
        l=lptr_p25_aer(1,2,iphase)
        if (l .ge. p1st)  mass_oini= chem(i,k,j,l)*conv1a
        l=lptr_orgaro1_aer(1,2,iphase)
        if (l .ge. p1st)  mass_aro1i= chem(i,k,j,l)*conv1a
        l=lptr_orgaro2_aer(1,2,iphase)
        if (l .ge. p1st)  mass_aro2i= chem(i,k,j,l)*conv1a
        l=lptr_orgalk_aer(1,2,iphase)
        if (l .ge. p1st)  mass_alk1i= chem(i,k,j,l)*conv1a
        l=lptr_orgole_aer(1,2,iphase)
        if (l .ge. p1st)  mass_ole1i= chem(i,k,j,l)*conv1a
        l=lptr_orgba1_aer(1,2,iphase)
        if (l .ge. p1st)  mass_ba1i= chem(i,k,j,l)*conv1a
        l=lptr_orgba2_aer(1,2,iphase)
        if (l .ge. p1st)  mass_ba2i= chem(i,k,j,l)*conv1a
        l=lptr_orgba3_aer(1,2,iphase)
        if (l .ge. p1st)  mass_ba3i= chem(i,k,j,l)*conv1a
        l=lptr_orgba4_aer(1,2,iphase)
        if (l .ge. p1st)  mass_ba4i= chem(i,k,j,l)*conv1a
        l=lptr_orgpa_aer(1,2,iphase)
        if (l .ge. p1st)  mass_pai= chem(i,k,j,l)*conv1a
        l=lptr_ec_aer(1,2,iphase)
        if (l .ge. p1st)  mass_bci= chem(i,k,j,l)*conv1a
        l=numptr_aer(1,2,iphase)
        if (l .ge. p1st)  num_ai= chem(i,k,j,l)*conv1b
        mass_h2oi= h2oai(i,k,j) * 1.0e-12
        mass_oci=mass_aro1i+mass_aro2i+mass_alk1i+mass_ole1i+ &
                 mass_ba1i+mass_ba2i+mass_ba3i+mass_ba4i+mass_pai

! Coarse mode...
        l=lptr_anth_aer(1,3,iphase)
        if (l .ge. p1st)  mass_antha= chem(i,k,j,l)*conv1a
        l=lptr_seas_aer(1,3,iphase)
        if (l .ge. p1st)  mass_seas= chem(i,k,j,l)*conv1a
        l=lptr_soil_aer(1,3,iphase)
        if (l .ge. p1st)  mass_soil= chem(i,k,j,l)*conv1a
        l=numptr_aer(1,3,iphase)
        if (l .ge. p1st)  num_ac= chem(i,k,j,l)*conv1b

        vol_ai = (mass_so4i/dens_so4)+(mass_no3i/dens_no3)+ &
                 (mass_nh4i/dens_nh4)+(mass_oini/dens_oin)+ &
                 (mass_aro1i/dens_oc)+(mass_alk1i/dens_oc)+ &
                 (mass_ole1i/dens_oc)+(mass_ba1i/dens_oc)+  &
                 (mass_ba2i/dens_oc)+(mass_ba3i/dens_oc)+   &
                 (mass_ba4i/dens_oc)+(mass_pai/dens_oc)+    &
                 (mass_aro2i/dens_oc)+(mass_bci/dens_bc)
        vol_aj = (mass_so4j/dens_so4)+(mass_no3j/dens_no3)+ &
                 (mass_nh4j/dens_nh4)+(mass_oinj/dens_oin)+ &
                 (mass_aro1j/dens_oc)+(mass_alk1j/dens_oc)+ &
                 (mass_ole1j/dens_oc)+(mass_ba1j/dens_oc)+  &
                 (mass_ba2j/dens_oc)+(mass_ba3j/dens_oc)+   &
                 (mass_ba4j/dens_oc)+(mass_paj/dens_oc)+    &
                 (mass_aro2j/dens_oc)+(mass_bcj/dens_bc)
        vol_ac = (mass_antha/dens_oin)+ &
                 (mass_seas*(22.9897/58.4428)/dens_na)+ &
                 (mass_seas*(35.4270/58.4428)/dens_cl)+ &
                 (mass_soil/dens_oin)
!
! Now divide mass into sections which is done by sect02:
! * xmas_secti is for aiken mode
! * xmas_sectj is for accumulation mode
! * xmas_sectc is for coarse mode
! * sect02 expects input in um
! * pass in generic mass of 1.0 just to get a percentage distribution of mass among bins
!
        ss1=alog(sginin)
        ss2=exp(ss1*ss1*36.0/8.0) 
        ss3=(sixpi*vol_ai/(num_ai*ss2))**0.3333333
        dgnum_um=max(dgmin,ss3)*1.0e+04
        call sect02(dgnum_um,sginin,drydens,iflag,duma,nbin_o,dlo_um,dhi_um, &
          xnum_secti,xmas_secti)
        ss1=alog(sginia)
        ss2=exp(ss1*ss1*36.0/8.0) 
        ss3=(sixpi*vol_aj/(num_aj*ss2))**0.3333333
        dgnum_um=max(dgmin,ss3)*1.0e+04
        call sect02(dgnum_um,sginia,drydens,iflag,duma,nbin_o,dlo_um,dhi_um, &
          xnum_sectj,xmas_sectj)
        ss1=alog(sginic)
        ss2=exp(ss1*ss1*36.0/8.0) 
        ss3=(sixpi*vol_ac/(num_ac*ss2))**0.3333333
        dgnum_um=max(dgmin,ss3)*1.0e+04
        call sect02(dgnum_um,sginic,drydens,iflag,duma,nbin_o,dlo_um,dhi_um, &
          xnum_sectc,xmas_sectc)

        do isize = 1, nbin_o
          xdia_cm(isize)=xdia_um(isize)*1.0e-04
          mass_so4 = mass_so4i*xmas_secti(isize) + mass_so4j*xmas_sectj(isize)
          mass_no3 = mass_no3i*xmas_secti(isize) + mass_no3j*xmas_sectj(isize)
          mass_nh4 = mass_nh4i*xmas_secti(isize) + mass_nh4j*xmas_sectj(isize)
          mass_oin = mass_oini*xmas_secti(isize) + mass_oinj*xmas_sectj(isize) + &
                     mass_soil*xmas_sectc(isize) + mass_antha*xmas_sectc(isize)
          mass_oc = (mass_pai+mass_aro1i+mass_aro2i+mass_alk1i+mass_ole1i+ &
                     mass_ba1i+mass_ba2i+mass_ba3i+mass_ba4i)*xmas_secti(isize) + &
                    (mass_paj+mass_aro1j+mass_aro2j+mass_alk1j+mass_ole1j+ &
                     mass_ba1j+mass_ba2j+mass_ba3j+mass_ba4j)*xmas_sectj(isize)
          mass_bc  = mass_bci*xmas_secti(isize) + mass_bcj*xmas_sectj(isize)
          mass_na = mass_seas*xmas_sectc(isize)*(22.9897/58.4428)
          mass_cl = mass_seas*xmas_sectc(isize)*(35.4270/58.4428)
          mass_h2o = mass_h2oi*xmas_secti(isize) + mass_h2oj*xmas_sectj(isize)
!         mass_h2o = 0.0 ! testing purposes only
          vol_so4 = mass_so4 / dens_so4
          vol_no3 = mass_no3 / dens_no3
          vol_nh4 = mass_nh4 / dens_nh4
          vol_oin = mass_oin / dens_oin
          vol_oc  = mass_oc  / dens_oc
          vol_bc  = mass_bc  / dens_bc
          vol_na  = mass_na  / dens_na
          vol_cl  = mass_cl  / dens_cl
          vol_h2o = mass_h2o / dens_h2o
          mass_dry_a = mass_so4 + mass_no3 + mass_nh4 + mass_oin + &
                       mass_oc  + mass_bc  + mass_na  + mass_cl
          mass_wet_a = mass_dry_a + mass_h2o 
          vol_dry_a  = vol_so4  + vol_no3  + vol_nh4  + vol_oin  + &
                       vol_oc   + vol_bc   + vol_na   + vol_cl
          vol_wet_a  = vol_dry_a + vol_h2o
          vol_shell  = vol_wet_a - vol_bc
          num_a      = vol_wet_a / (0.52359877*xdia_cm(isize)*xdia_cm(isize)*xdia_cm(isize))
          ri_dum     = (0.0,0.0)
          ri_dum     = (ref_index_nh4so4 * mass_so4 / dens_so4) +  &
                       (ref_index_nh4no3 * mass_no3 / dens_no3) +  &
                       (ref_index_nh4no3 * mass_nh4 / dens_nh4) +  &
                       (ref_index_oin    * mass_oin / dens_oin) +  &
                       (ref_index_oc     * mass_oc  / dens_oc) +   &
                       (ref_index_bc     * mass_bc  / dens_bc) +   &
                       (ref_index_nacl   * mass_na  / dens_na) +   &
                       (ref_index_nacl   * mass_cl  / dens_cl) +   &
                       (ref_index_h2o    * mass_h2o / dens_h2o) 
!
! for some reason MADE/SORGAM occasionally produces zero aerosols so
! need to add a check here to avoid divide by zero
!
          IF(num_a .lt. 1.0e-20 .or. vol_wet_a .lt. 1.0e-20) then
            dp_dry_a   = 0.0
            dp_wet_a   = 0.0
            dp_bc_a    = 0.0
            ri_ave_a   = 0.0
            ri_dum     = 0.0
          else
            dp_dry_a   = (1.90985*vol_dry_a/num_a)**0.3333333
            dp_wet_a   = (1.90985*vol_wet_a/num_a)**0.3333333
            dp_bc_a    = (1.90985*vol_bc/num_a)**0.3333333
            ri_ave_a   = ri_dum/vol_wet_a
            ri_dum   = (ref_index_nh4so4 * mass_so4 / dens_so4) +  &
                       (ref_index_nh4no3 * mass_no3 / dens_no3) +  &
                       (ref_index_nh4no3 * mass_nh4 / dens_nh4) +  &
                       (ref_index_oin    * mass_oin / dens_oin) +  &
                       (ref_index_oc     * mass_oc  / dens_oc) +   &
                       (ref_index_nacl   * mass_na  / dens_na) +   &
                       (ref_index_nacl   * mass_cl  / dens_cl) +   &
                       (ref_index_h2o    * mass_h2o / dens_h2o) 
          endif
          if(dp_wet_a/2.0 .lt. dlo_um*1.0e-4/2.0) then
            refindx(i,k,j,isize)    = (1.5,0.0)
            radius_wet(i,k,j,isize) =dlo_um*1.0e-4/2.0
            number_bin(i,k,j,isize) =num_a
            radius_core(i,k,j,isize) =0.0
            refindx_core(i,k,j,isize) = ref_index_bc
            refindx_shell(i,k,j,isize) = ref_index_oin
          elseif(vol_shell .lt. 1.0e-20) then
            refindx(i,k,j,isize)    = (1.5,0.0)
            radius_wet(i,k,j,isize) =dlo_um*1.0e-4/2.0
            number_bin(i,k,j,isize) =num_a
            radius_core(i,k,j,isize) =0.0
            refindx_core(i,k,j,isize) = ref_index_bc
            refindx_shell(i,k,j,isize) = ref_index_oin
          else
            refindx(i,k,j,isize)    =ri_ave_a
            radius_wet(i,k,j,isize) =dp_wet_a/2.0
            number_bin(i,k,j,isize) =num_a
            radius_core(i,k,j,isize) =dp_bc_a/2.0
            refindx_core(i,k,j,isize) =ref_index_bc
            refindx_shell(i,k,j,isize) =ri_dum/vol_shell
          endif
!
          refr=real(refindx(i,k,j,isize))
        enddo
      enddo
      enddo
      enddo

      return

      end subroutine optical_prep_modal
!
!***********************************************************************
! <1.> subr mieaer
! Purpose:  calculate aerosol optical depth, single scattering albedo,
!   asymmetry factor, extinction, Legendre coefficients, and average aerosol
!   size. parameterizes aerosol coefficients using chebychev polynomials
!   requires double precision on 32-bit machines
!   uses Wiscombe's (1979) mie scattering code
! INPUT
!       id -- grid id number
!   	iclm, jclm -- i,j of grid column being processed
!       nbin_a -- number of bins
!	number_bin(nbin_a,kmaxd) --   number density in layer, #/cm^3
!	radius_wet(nbin_a,kmaxd) -- wet radius, cm
!    	refindx(nbin_a,kmaxd) --volume averaged complex index of refraction
!	dz -- depth of individual cells in column, m
!	isecfrm0 - time from start of run, sec
!	lpar -- number of grid cells in vertical (via module_fastj_cmnh)
!   kmaxd -- predefined maximum allowed levels from module_data_mosaic_other
!            passed here via module_fastj_cmnh
! OUTPUT: saved in module_fastj_cmnmie
!   real tauaer  ! aerosol optical depth
!        waer    ! aerosol single scattering albedo
!        gaer    ! aerosol asymmetery factor
!        extaer  ! aerosol extinction
!        l2,l3,l4,l5,l6,l7 ! Legendre coefficients, numbered 0,1,2,......
!        sizeaer ! average wet radius
!	 bscoef ! aerosol backscatter coefficient with units km-1 * steradian -1  JCB 2007/02/01
!***********************************************************************
        subroutine mieaer( &
	          id, iclm, jclm, nbin_a,   &
              number_bin_col, radius_wet_col, refindx_col,   &
              dz, isecfrm0, lpar, &
              sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7,bscoef)  ! added bscoef JCB 2007/02/01

!       USE module_data_mosaic_other, only : kmaxd
!       USE module_data_mosaic_therm, ONLY: nbin_a_maxd
        USE module_peg_util, only:  peg_message
        
        IMPLICIT NONE
!   subr arguments        
!jdf
        integer,parameter :: nspint = 4 ! Num of spectral intervals across
                                        ! solar spectrum for FAST-J
        integer, intent(in) :: lpar
!jdf    real, dimension (nspint, kmaxd+1),intent(out) :: sizeaer,extaer,waer,gaer,tauaer
!jdf    real, dimension (nspint, kmaxd+1),intent(out) :: l2,l3,l4,l5,l6,l7
!jdf    real, dimension (nspint, kmaxd+1),intent(out) :: bscoef  !JCB 2007/02/01
        real, dimension (nspint, lpar+1),intent(out) :: sizeaer,extaer,waer,gaer,tauaer
        real, dimension (nspint, lpar+1),intent(out) :: l2,l3,l4,l5,l6,l7
        real, dimension (nspint, lpar+1),intent(out) :: bscoef  !JCB 2007/02/01
        real, dimension (nspint),save :: wavmid !cm
        data wavmid     &
            / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
!jdf
    integer, intent(in) :: id, iclm, jclm, nbin_a, isecfrm0
!jdf    real, intent(in), dimension(nbin_a, kmaxd) :: number_bin_col
!jdf    real, intent(inout), dimension(nbin_a, kmaxd) :: radius_wet_col
!jdf    complex, intent(in) :: refindx_col(nbin_a, kmaxd)
    real, intent(in), dimension(nbin_a, lpar+1) :: number_bin_col
    real, intent(inout), dimension(nbin_a, lpar+1) :: radius_wet_col
    complex, intent(in) :: refindx_col(nbin_a, lpar+1)
    real, intent(in)    :: dz(lpar)		

!local variables
	real weighte, weights
! various bookeeping variables
	integer ltype ! total number of indicies of refraction
	parameter (ltype = 1)  ! bracket refractive indices based on information from Rahul, 2002/11/07
	real x
	real thesum ! for normalizing things
	real sizem ! size in microns
	integer kcallmieaer
!
	integer m, j, nc, klevel
        real pext	    ! parameterized specific extinction (cm2/g)
        real pasm       ! parameterized asymmetry factor
        real ppmom2     ! 2 Lengendre expansion coefficient (numbered 0,1,2,...)
        real ppmom3     ! 3     ...
        real ppmom4     ! 4     ...
        real ppmom5     ! 5     ...
        real ppmom6     ! 6     ...
        real ppmom7     ! 7     ...
        real sback2     ! JCB 2007/02/01 sback*conjg(sback)

        integer ns 	    ! Spectral loop index
        integer i  	    ! Longitude loop index
        integer k       ! Level loop index
        real pscat      !scattering cross section

        integer prefr,prefi,nrefr,nrefi,nr,ni
        save nrefr,nrefi
        parameter (prefr=7,prefi=7)

        complex*16 sforw,sback,tforw(2),tback(2)
        integer numang,nmom,ipolzn,momdim
        real*8 pmom(0:7,1)
        logical perfct,anyang,prnt(2)
        logical, save :: first
        data first/.true./
        integer, parameter ::  nsiz=200,nlog=30,ncoef=50
!
        real p2(nsiz),p3(nsiz),p4(nsiz),p5(nsiz)
        real p6(nsiz),p7(nsiz)
!
        real*8 xmu(1)
        data xmu/1./,anyang/.false./
        data numang/0/
        complex*16 s1(1),s2(1)
        real*8 mimcut
        data perfct/.false./,mimcut/0.0/
        data nmom/7/,ipolzn/0/,momdim/7/
        data prnt/.false.,.false./
!     coefficients for parameterizing aerosol radiative properties
!     in terms of refractive index and wet radius
      real, save :: extp(ncoef,prefr,prefi,nspint)	! specific extinction
      real, save :: albp(ncoef,prefr,prefi,nspint)	! single scat alb
      real, save :: asmp(ncoef,prefr,prefi,nspint)	! asymmetry factor
      real, save :: ascat(ncoef,prefr,prefi,nspint)	! scattering efficiency, JCB 2004/02/09
      real, save :: pmom2(ncoef,prefr,prefi,nspint)	! phase function expansion, #2
      real, save :: pmom3(ncoef,prefr,prefi,nspint)	! phase function expansion, #3
      real, save :: pmom4(ncoef,prefr,prefi,nspint)	! phase function expansion, #4
      real, save :: pmom5(ncoef,prefr,prefi,nspint)	! phase function expansion, #5
      real, save :: pmom6(ncoef,prefr,prefi,nspint)	! phase function expansion, #6
      real, save :: pmom7(ncoef,prefr,prefi,nspint) ! phase function expansion, #7
      real, save :: sback2p(ncoef,prefr,prefi,nspint) ! JCB 2007/02/01 - backscatter
!--------------
      real cext(ncoef),casm(ncoef),cpmom2(ncoef)
      real cscat(ncoef)  ! JCB 2004/02/09
      real cpmom3(ncoef),cpmom4(ncoef),cpmom5(ncoef)
      real cpmom6(ncoef),cpmom7(ncoef)
      real cpsback2p(ncoef) ! JCB 2007/02/09  - backscatter
      integer itab,jtab
      real ttab,utab

!     nsiz = number of wet particle sizes
!     crefin = complex refractive index
      integer n
      real*8 thesize 	! 2 pi radpart / waveleng = size parameter
      real*8 qext(nsiz) ! array of extinction efficiencies
      real*8 qsca(nsiz) ! array of scattering efficiencies
      real*8 gqsc(nsiz) ! array of asymmetry factor * scattering efficiency
      real asymm(nsiz)  ! array of asymmetry factor
      real scat(nsiz)   ! JCB 2004/02/09
      real sb2(nsiz)     ! JCB 2007/02/01 - 4*abs(sback)^2/(size parameter)^2 backscattering efficiency
!     specabs = absorption coeff / unit dry mass
!     specscat = scattering coeff / unit dry mass
      complex*16 crefin,crefd,crefw
      save crefw
      real, save :: rmin,rmax   ! min, max aerosol size bin
      data rmin /0.005e-4/ 	! rmin in cm. 5e-3 microns min allowable size
      data rmax /50.e-4 /	! rmax in cm. 50 microns, big particle, max allowable size

      real bma,bpa

      real, save :: xrmin,xrmax,xr
      real rs(nsiz) ! surface mode radius (cm)
      real xrad ! normalized aerosol radius
      real ch(ncoef) ! chebychev polynomial

      real, save :: rhoh2o ! density of liquid water (g/cm3)
      data rhoh2o/1./

      real refr     ! real part of refractive index
      real refi     ! imaginary part of refractive index
      real qextr4(nsiz)          !  extinction, real*4

      real refrmin ! minimum of real part of refractive index
      real refrmax ! maximum of real part of refractive index
      real refimin ! minimum of imag part of refractive index
      real refimax ! maximum of imag part of refractive index
      real drefr ! increment in real part of refractive index
      real drefi ! increment in imag part of refractive index
      real, save :: refrtab(prefr) ! table of real refractive indices for aerosols
      real, save :: refitab(prefi) ! table of imag refractive indices for aerosols
	complex specrefndx(ltype) ! refractivr indices
      real pie,third

      integer irams, jrams
! diagnostic declarations
      integer kcallmieaer2
      integer ibin
      character*150 msg

#if (defined(CHEM_DBG_I) && defined(CHEM_DBG_J) && defined(CHEM_DBG_K))
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ec  diagnostics
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ec run_out.25 has aerosol physical parameter info for bins 1-8
!ec and vertical cells 1 to kmaxd.
!        ilaporte = 33
!        jlaporte = 34
	if (iclm .eq. CHEM_DBG_I) then
	  if (jclm .eq. CHEM_DBG_J) then
!   initial entry
	   if (kcallmieaer2 .eq. 0) then
	      write(*,9099)iclm, jclm
 9099	format('for cell i = ', i3, 2x, 'j = ', i3)	
              write(*,9100)
 9100     format(   &
               'isecfrm0', 3x, 'i', 3x, 'j', 3x,'k', 3x,   &
      	       'ibin', 3x,   &
               'refindx_col(ibin,k)', 3x,   &
               'radius_wet_col(ibin,k)', 3x,   &
               'number_bin_col(ibin,k)'   &
               )
	   end if
!ec output for run_out.25
            do k = 1, lpar
            do ibin = 1, nbin_a
              write(*, 9120)   &
                 isecfrm0,iclm, jclm, k, ibin,   &
                 refindx_col(ibin,k),   &
                 radius_wet_col(ibin,k),   &
                 number_bin_col(ibin,k)
9120    format( i7,3(2x,i4),2x,i4, 4x, 4(e14.6,2x))
            end do
            end do
	kcallmieaer2 = kcallmieaer2 + 1
        end if
        end if
!ec end print of aerosol physical parameter diagnostics	
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
#endif
!
! assign fast-J wavelength, these values are in cm
!      wavmid(1)=0.30e-4
!      wavmid(2)=0.40e-4
!      wavmid(3)=0.60e-4
!      wavmid(4)=0.999e-04
!
      pie=4.*atan(1.)
      third=1./3.
      if(first)then
        first=.false.

!       parameterize aerosol radiative properties in terms of
!       relative humidity, surface mode wet radius, aerosol species,
!       and wavelength

!       first find min,max of real and imaginary parts of refractive index

!       real and imaginary parts of water refractive index

        crefw=cmplx(1.33,0.0)
        refrmin=real(crefw)
        refrmax=real(crefw)
! change Rahul's imaginary part of the refractive index from positive to negative
        refimin=-imag(crefw)
        refimax=-imag(crefw)

!       aerosol mode loop
        specrefndx(1)=cmplx(1.85,-0.71) ! max values from Bond and Bergstrom
!
	
	do i=1,ltype ! loop over all possible refractive indices

!             real and imaginary parts of aerosol refractive index

              refrmin=amin1(refrmin,real(specrefndx(ltype)))
              refrmax=amax1(refrmax,real(specrefndx(ltype)))
              refimin=amin1(refimin,aimag(specrefndx(ltype)))
              refimax=amax1(refimax,aimag(specrefndx(ltype)))

           enddo

         drefr=(refrmax-refrmin)
         if(drefr.gt.1.e-4)then
            nrefr=prefr
            drefr=drefr/(nrefr-1)
         else
            nrefr=1
         endif

         drefi=(refimax-refimin)
         if(drefi.gt.1.e-4)then
            nrefi=prefi
            drefi=drefi/(nrefi-1)
         else
            nrefi=1
         endif

!
               bma=0.5*alog(rmax/rmin) ! JCB
               bpa=0.5*alog(rmax*rmin) ! JCB
	xrmin=alog(rmin)
	xrmax=alog(rmax)

!        wavelength loop

         do 200 ns=1,nspint

!           calibrate parameterization with range of refractive indices

            do 120 nr=1,nrefr
            do 120 ni=1,nrefi

               refrtab(nr)=refrmin+(nr-1)*drefr
! JCB, Feb. 20, 2008:   here we fix a conceptual problem associated with
! the old Ghan code. See Steve Ghan's updated method;  Ghan and Zaveri,
! Parameterization of optical properties for hydrated internally mixed
! aerosol, J. Geophys. Res., 112, D10201, doi:10.1029/2006JD007927, 2007.
! JCB      refitab(ni)=refimin+(ni-1)*drefi, old code
! JCB, note that the complex part of the refractive index is < 0
		refitab(ni)=refimin/0.2*(0.2**real(ni))  ! JCB change, slightly different parameters than Ghan and Zaveri
		if(ni.eq.nrefi)refitab(nrefi)=0.0  ! JCB change
               crefd=cmplx(refrtab(nr),refitab(ni))

!              mie calculations of optical efficiencies

               do n=1,nsiz
                  xr=cos(pie*(float(n)-0.5)/float(nsiz))
                  rs(n)=exp(xr*bma+bpa)

!                 size parameter and weighted refractive index

                  thesize=2.*pie*rs(n)/wavmid(ns)
                  thesize=min(thesize,10000.d0)

                  call miev0(thesize,crefd,perfct,mimcut,anyang,   &
                     numang,xmu,nmom,ipolzn,momdim,prnt,   &
                     qext(n),qsca(n),gqsc(n),pmom,sforw,sback,s1,   &
                     s2,tforw,tback )
		  qextr4(n)=qext(n)
!                 qabs(n)=qext(n)-qsca(n) ! not necessary anymore JCB 2004/02/09
                  scat(n)=qsca(n) ! JCB 2004/02/09
                  asymm(n)=gqsc(n)/qsca(n) ! assume always greater than zero
! coefficients of phase function expansion; note modification by JCB of miev0 coefficients
		  p2(n)=pmom(2,1)/pmom(0,1)*5.0
		  p3(n)=pmom(3,1)/pmom(0,1)*7.0
		  p4(n)=pmom(4,1)/pmom(0,1)*9.0
		  p5(n)=pmom(5,1)/pmom(0,1)*11.0
		  p6(n)=pmom(6,1)/pmom(0,1)*13.0
		  p7(n)=pmom(7,1)/pmom(0,1)*15.0
! backscattering efficiency, Bohren and Huffman, page 122
! as stated by Bohren and Huffman, this is 4*pie times what is should be
! may need to be smoothed - a very rough function - for the time being we won't apply smoothing
! and let the integration over the size distribution be the smoothing
		  sb2(n)=4.0*sback*dconjg(sback)/(thesize*thesize) ! JCB 2007/02/01  
		enddo
  100          continue
!
               call fitcurv(rs,qextr4,extp(1,nr,ni,ns),ncoef,nsiz)
               call fitcurv(rs,scat,ascat(1,nr,ni,ns),ncoef,nsiz) ! JCB 2004/02/07 - scattering efficiency
               call fitcurv(rs,asymm,asmp(1,nr,ni,ns),ncoef,nsiz)
               call fitcurv_nolog(rs,p2,pmom2(1,nr,ni,ns),ncoef,nsiz)
	       call fitcurv_nolog(rs,p3,pmom3(1,nr,ni,ns),ncoef,nsiz)
	       call fitcurv_nolog(rs,p4,pmom4(1,nr,ni,ns),ncoef,nsiz)
               call fitcurv_nolog(rs,p5,pmom5(1,nr,ni,ns),ncoef,nsiz)
	       call fitcurv_nolog(rs,p6,pmom6(1,nr,ni,ns),ncoef,nsiz)
               call fitcurv_nolog(rs,p7,pmom7(1,nr,ni,ns),ncoef,nsiz)
               call fitcurv(rs,sb2,sback2p(1,nr,ni,ns),ncoef,nsiz) ! JCB 2007/02/01 - backscattering efficiency             

  120       continue
  200    continue

      endif
! begin level loop
	do 2000 klevel=1,lpar
! sum densities for normalization
	thesum=0.0
	do m=1,nbin_a
	thesum=thesum+number_bin_col(m,klevel)
	enddo
! Begin spectral loop
      do 1000 ns=1,nspint

!        aerosol optical properties

               tauaer(ns,klevel)=0.
               waer(ns,klevel)=0.
               gaer(ns,klevel)=0.
	       sizeaer(ns,klevel)=0.0
		extaer(ns,klevel)=0.0
		l2(ns,klevel)=0.0
		l3(ns,klevel)=0.0
		l4(ns,klevel)=0.0
		l5(ns,klevel)=0.0
		l6(ns,klevel)=0.0
		l7(ns,klevel)=0.0
		bscoef(ns,klevel)=0.0  ! JCB 2007/02/01 - backscattering coefficient
		if(thesum.le.1e-21)goto 1000 ! set everything = 0 if no aerosol !wig changed 0.0 to 1e-21, 31-Oct-2005

! loop over the bins
               do m=1,nbin_a ! nbin_a is number of bins
! check to see if there's any aerosol
!jdf   		if(number_bin_col(m,klevel).le.1e-21)goto 70  ! no aerosol !wig changed 0.0 to 1e-21, 31-Oct-2005
! here's the size
		sizem=radius_wet_col(m,klevel) ! radius in cm
! check limits of particle size
! rce 2004-dec-07 - use klevel in write statements
		if(radius_wet_col(m,klevel).le.rmin)then
		  radius_wet_col(m,klevel)=rmin
		  write( msg, '(a, 5i4,1x, e11.4)' )	&
		  'mieaer: radius_wet set to rmin,'  //	&
		  'id,i,j,k,m,rm(m,k)', id, iclm, jclm, klevel, m, radius_wet_col(m,klevel)
		  call peg_message( lunerr, msg )
!		   write(6,'('' particle size too small '')')
		endif
!
		if(radius_wet_col(m,klevel).gt.rmax)then
		   write( msg, '(a, 5i4,1x, e11.4)' )	&
                'mieaer: radius_wet set to rmax,'  //	&
                'id,i,j,k,m,rm(m,k)', &
                id, iclm, jclm, klevel, m, radius_wet_col(m,klevel)
           call peg_message( lunerr, msg )
           radius_wet_col(m,klevel)=rmax
!           write(6,'('' particle size too large '')')
		endif
!
		x=alog(radius_wet_col(m,klevel)) ! radius in cm
!
	          crefin=refindx_col(m,klevel)
                  refr=real(crefin)
! change Rahul's imaginary part of the index of refraction from positive to negative
                  refi=-imag(crefin)
!
		xrad=x

		thesize=2.0*pie*exp(x)/wavmid(ns)
! normalize size parameter
                   xrad=(2*xrad-xrmax-xrmin)/(xrmax-xrmin)
! retain this diagnostic code
		  if(abs(refr).gt.10.0.or.abs(refr).le.0.001)then
                     write ( msg, '(a,1x, e14.5)' )  &
           'mieaer /refr/ outside range 1e-3 - 10 ' //  &
           'refr= ', refr
         	     call peg_message( lunerr, msg )
!                      print *,'refr=',refr
                      call exit(1)
                  endif
                  if(abs(refi).gt.10.)then
                     write ( msg, '(a,1x, e14.5)' )  &
           'mieaer /refi/ >10 '  //  &
            'refi', refi
         	     call peg_message( lunerr, msg )                  
!                      print *,'refi=',refi
                      call exit(1)
                  endif

! interpolate coefficients linear in refractive index
! first call calcs itab,jtab,ttab,utab
                  itab=0
                  call binterp(extp(1,1,1,ns),ncoef,nrefr,nrefi,   &
                               refr,refi,refrtab,refitab,itab,jtab,   &
                               ttab,utab,cext)

! JCB 2004/02/09  -- new code for scattering cross section
                  call binterp(ascat(1,1,1,ns),ncoef,nrefr,nrefi,   &
                               refr,refi,refrtab,refitab,itab,jtab,   &
                               ttab,utab,cscat)
                  call binterp(asmp(1,1,1,ns),ncoef,nrefr,nrefi,   &
                               refr,refi,refrtab,refitab,itab,jtab,   &
                               ttab,utab,casm)
                  call binterp(pmom2(1,1,1,ns),ncoef,nrefr,nrefi,   &
                               refr,refi,refrtab,refitab,itab,jtab,   &
                               ttab,utab,cpmom2)
                  call binterp(pmom3(1,1,1,ns),ncoef,nrefr,nrefi,   &
                               refr,refi,refrtab,refitab,itab,jtab,   &
                               ttab,utab,cpmom3)
                  call binterp(pmom4(1,1,1,ns),ncoef,nrefr,nrefi,   &
                               refr,refi,refrtab,refitab,itab,jtab,   &
                               ttab,utab,cpmom4)
                  call binterp(pmom5(1,1,1,ns),ncoef,nrefr,nrefi,   &
                               refr,refi,refrtab,refitab,itab,jtab,   &
                               ttab,utab,cpmom5)
                  call binterp(pmom6(1,1,1,ns),ncoef,nrefr,nrefi,   &
                               refr,refi,refrtab,refitab,itab,jtab,   &
                               ttab,utab,cpmom6)
                  call binterp(pmom7(1,1,1,ns),ncoef,nrefr,nrefi,   &
                               refr,refi,refrtab,refitab,itab,jtab,   &
                               ttab,utab,cpmom7)
                  call binterp(sback2p(1,1,1,ns),ncoef,nrefr,nrefi,   &
                               refr,refi,refrtab,refitab,itab,jtab,   &
                               ttab,utab,cpsback2p)

!                 chebyshev polynomials

                  ch(1)=1.
                  ch(2)=xrad
                  do nc=3,ncoef
                     ch(nc)=2.*xrad*ch(nc-1)-ch(nc-2)
                  enddo
!                 parameterized optical properties

                     pext=0.5*cext(1)
                     do nc=2,ncoef
                        pext=pext+ch(nc)*cext(nc)
                     enddo
                    pext=exp(pext)
	
! JCB 2004/02/09 -- for scattering efficiency
                  pscat=0.5*cscat(1)
                  do nc=2,ncoef
                     pscat=pscat+ch(nc)*cscat(nc)
                  enddo
                  pscat=exp(pscat)
!
                  pasm=0.5*casm(1)
                  do nc=2,ncoef
                     pasm=pasm+ch(nc)*casm(nc)
                  enddo
                  pasm=exp(pasm)
!
                  ppmom2=0.5*cpmom2(1)
                  do nc=2,ncoef
                     ppmom2=ppmom2+ch(nc)*cpmom2(nc)
                  enddo
		  if(ppmom2.le.0.0)ppmom2=0.0
!
                  ppmom3=0.5*cpmom3(1)
                  do nc=2,ncoef
                     ppmom3=ppmom3+ch(nc)*cpmom3(nc)
                  enddo
!		  ppmom3=exp(ppmom3)  ! no exponentiation required
		  if(ppmom3.le.0.0)ppmom3=0.0
!
                  ppmom4=0.5*cpmom4(1)
                  do nc=2,ncoef
                     ppmom4=ppmom4+ch(nc)*cpmom4(nc)
                  enddo
		  if(ppmom4.le.0.0.or.sizem.le.0.03e-04)ppmom4=0.0
!
                  ppmom5=0.5*cpmom5(1)
                  do nc=2,ncoef
                     ppmom5=ppmom5+ch(nc)*cpmom5(nc)
                  enddo
		  if(ppmom5.le.0.0.or.sizem.le.0.03e-04)ppmom5=0.0
!
                  ppmom6=0.5*cpmom6(1)
                  do nc=2,ncoef
                     ppmom6=ppmom6+ch(nc)*cpmom6(nc)
                  enddo
		  if(ppmom6.le.0.0.or.sizem.le.0.03e-04)ppmom6=0.0
!
                  ppmom7=0.5*cpmom7(1)
                  do nc=2,ncoef
                     ppmom7=ppmom7+ch(nc)*cpmom7(nc)
                  enddo
		  if(ppmom7.le.0.0.or.sizem.le.0.03e-04)ppmom7=0.0
!
                  sback2=0.5*cpsback2p(1) ! JCB 2007/02/01 - backscattering efficiency
                  do nc=2,ncoef
                     sback2=sback2+ch(nc)*cpsback2p(nc)
                  enddo
                     sback2=exp(sback2)
		  if(sback2.le.0.0)sback2=0.0
!
!
! weights:
	weighte=pext*pie*exp(x)**2 ! JCB, extinction cross section
	weights=pscat*pie*exp(x)**2 ! JCB, scattering cross section
       	tauaer(ns,klevel)=tauaer(ns,klevel)+weighte*number_bin_col(m,klevel)  ! must be multiplied by deltaZ
!	extaer(ns,klevel)=extaer(ns,klevel)+pext*number_bin_col(m,klevel) ! not the true extinction, JCB 2007/02/01
	sizeaer(ns,klevel)=sizeaer(ns,klevel)+exp(x)*10000.0*   &
      	number_bin_col(m,klevel)
        waer(ns,klevel)=waer(ns,klevel)+weights*number_bin_col(m,klevel) !JCB
        gaer(ns,klevel)=gaer(ns,klevel)+pasm*weights*   &
        number_bin_col(m,klevel) !JCB
! need weighting by scattering cross section ?  JCB 2004/02/09
	l2(ns,klevel)=l2(ns,klevel)+weights*ppmom2*number_bin_col(m,klevel)
	l3(ns,klevel)=l3(ns,klevel)+weights*ppmom3*number_bin_col(m,klevel)
	l4(ns,klevel)=l4(ns,klevel)+weights*ppmom4*number_bin_col(m,klevel)
	l5(ns,klevel)=l5(ns,klevel)+weights*ppmom5*number_bin_col(m,klevel)
	l6(ns,klevel)=l6(ns,klevel)+weights*ppmom6*number_bin_col(m,klevel)
	l7(ns,klevel)=l7(ns,klevel)+weights*ppmom7*number_bin_col(m,klevel)	
! convert backscattering efficiency to backscattering coefficient, units (cm)^-1
	bscoef(ns,klevel)=bscoef(ns,klevel)+pie*exp(x)**2*sback2*number_bin_col(m,klevel) ! JCB 2007/02/01 - backscatter coefficient,

        end do ! end of nbin_a loop
! take averages - weighted by cross section - new code JCB 2004/02/09
	sizeaer(ns,klevel)=sizeaer(ns,klevel)/thesum
	gaer(ns,klevel)=gaer(ns,klevel)/waer(ns,klevel) ! JCB removed *3 factor 2/9/2004
! because factor is applied in subroutine opmie, file zz01fastj_mod.f
	l2(ns,klevel)=l2(ns,klevel)/waer(ns,klevel)
	l3(ns,klevel)=l3(ns,klevel)/waer(ns,klevel)
	l4(ns,klevel)=l4(ns,klevel)/waer(ns,klevel)
	l5(ns,klevel)=l5(ns,klevel)/waer(ns,klevel)
	l6(ns,klevel)=l6(ns,klevel)/waer(ns,klevel)
	l7(ns,klevel)=l7(ns,klevel)/waer(ns,klevel)
! backscatter coef, divide by 4*Pie to get units of (km*ster)^-1 JCB 2007/02/01
	bscoef(ns,klevel)=bscoef(ns,klevel)*1.0e5  ! JCB 2007/02/01 - units are now (km)^-1
	extaer(ns,klevel)=tauaer(ns,klevel)*1.0e5  ! JCB 2007/02/01 - now true extincion, units (km)^-1
! this must be last!! JCB 2007/02/01
	waer(ns,klevel)=waer(ns,klevel)/tauaer(ns,klevel) ! JCB

 70   continue ! bail out if no aerosol;go on to next wavelength bin

 1000 continue  ! end of wavelength loop

2000   continue  ! end of klevel loop
!
! before returning, multiply tauaer by depth of individual cells.
! tauaer is in cm-1, dz in m; multiply dz by 100 to convert from m to cm.
  	do ns = 1, nspint
  	do klevel = 1, lpar
  	   tauaer(ns,klevel) = tauaer(ns,klevel) * dz(klevel)* 100.   
        end do
        end do 	

#if (defined(CHEM_DBG_I) && defined(CHEM_DBG_J) && defined(CHEM_DBG_K))
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ec  fastj diagnostics
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ec run_out.30 has aerosol optical info for cells 1 to kmaxd.
!        ilaporte = 33
!         jlaporte = 34
 	if (iclm .eq. CHEM_DBG_I) then
 	  if (jclm .eq. CHEM_DBG_J) then
!   initial entry
 	   if (kcallmieaer .eq. 0) then
 	       write(*,909) CHEM_DBG_I, CHEM_DBG_J
 909	format(	' for cell i = ', i3, ' j = ', i3)          	
               write(*,910)
 910     format(   &
               'isecfrm0', 3x, 'i', 3x, 'j', 3x,'k', 3x,   &
      	        'dzmfastj', 8x,   &
               'tauaer(1,k)',1x, 'tauaer(2,k)',1x,'tauaer(3,k)',3x,   &
               'tauaer(4,k)',5x,   &
               'waer(1,k)', 7x, 'waer(2,k)', 7x,'waer(3,k)', 7x,   &
               'waer(4,k)', 7x,   &
               'gaer(1,k)', 7x, 'gaer(2,k)', 7x,'gaer(3,k)', 7x,   &
               'gaer(4,k)', 7x,   &
               'extaer(1,k)',5x, 'extaer(2,k)',5x,'extaer(3,k)',5x,   &
               'extaer(4,k)',5x,   &
               'sizeaer(1,k)',4x, 'sizeaer(2,k)',4x,'sizeaer(3,k)',4x,   &
               'sizeaer(4,k)'  )
 	   end if
!ec output for run_out.30
         do k = 1, lpar
         write(*, 912)   &
           isecfrm0,iclm, jclm, k,   &
           dz(k) ,   &
           (tauaer(n,k),   n=1,4),   &
           (waer(n,k),     n=1,4),   &
           (gaer(n,k),     n=1,4),   &
           (extaer(n,k),   n=1,4),   &
           (sizeaer(n,k),  n=1,4)
 912    format( i7,3(2x,i4),2x,21(e14.6,2x))
         end do
 	kcallmieaer = kcallmieaer + 1
        end if
         end if
!ec end print of fastj diagnostics	
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
#endif

      return
      end subroutine mieaer
!****************************************************************
      subroutine fitcurv(rs,yin,coef,ncoef,maxm)

!     fit y(x) using Chebychev polynomials
!     wig 7-Sep-2004: Removed dependency on pre-determined maximum
!                     array size and replaced with f90 array info.

      USE module_peg_util, only:  peg_message

      IMPLICIT NONE
!      integer nmodes, nrows, maxm, ncoef
!      parameter (nmodes=500,nrows=8)
      integer, intent(in) :: maxm, ncoef

!      real rs(nmodes),yin(nmodes),coef(ncoef)
!      real x(nmodes),y(nmodes)
      real, dimension(ncoef) :: coef
      real, dimension(:) :: rs, yin
      real x(size(rs)),y(size(yin))

      integer m
      real xmin, xmax
      character*80 msg

!!$      if(maxm.gt.nmodes)then
!!$        write ( msg, '(a, 1x,i6)' )  &
!!$           'FASTJ mie nmodes too small in fitcurv, '  //  &
!!$           'maxm ', maxm
!!$          call peg_message( lunerr, msg )
!!$!        write(*,*)'nmodes too small in fitcurv',maxm
!!$        call exit(1)
!!$      endif

      do 100 m=1,maxm
      x(m)=alog(rs(m))
      y(m)=alog(yin(m))
  100 continue

      xmin=x(1)
      xmax=x(maxm)
      do 110 m=1,maxm
      x(m)=(2*x(m)-xmax-xmin)/(xmax-xmin)
  110 continue

      call chebft(coef,ncoef,maxm,y)

      return
      end subroutine fitcurv                        
!**************************************************************
      subroutine fitcurv_nolog(rs,yin,coef,ncoef,maxm)

!     fit y(x) using Chebychev polynomials
!     wig 7-Sep-2004: Removed dependency on pre-determined maximum
!                     array size and replaced with f90 array info.

      USE module_peg_util, only:  peg_message
      IMPLICIT NONE

!      integer nmodes, nrows, maxm, ncoef
!      parameter (nmodes=500,nrows=8)
      integer, intent(in) :: maxm, ncoef

!      real rs(nmodes),yin(nmodes),coef(ncoef)
      real, dimension(:) :: rs, yin
      real, dimension(ncoef) :: coef(ncoef)
      real x(size(rs)),y(size(yin))

      integer m
      real xmin, xmax
      character*80 msg
           
!!$      if(maxm.gt.nmodes)then
!!$        write ( msg, '(a,1x, i6)' )  &
!!$           'FASTJ mie nmodes too small in fitcurv '  //  &
!!$           'maxm ', maxm
!!$          call peg_message( lunerr, msg )
!!$!        write(*,*)'nmodes too small in fitcurv',maxm
!!$        call exit(1)
!!$      endif

      do 100 m=1,maxm
      x(m)=alog(rs(m))
      y(m)=yin(m) ! note, no "alog" here
  100 continue

      xmin=x(1)
      xmax=x(maxm)
      do 110 m=1,maxm
      x(m)=(2*x(m)-xmax-xmin)/(xmax-xmin)
  110 continue

      call chebft(coef,ncoef,maxm,y)

      return
      end subroutine fitcurv_nolog                        
!************************************************************************
      subroutine chebft(c,ncoef,n,f)
!     given a function f with values at zeroes x_k of Chebychef polynomial
!     T_n(x), calculate coefficients c_j such that
!     f(x)=sum(k=1,n) c_k t_(k-1)(y) - 0.5*c_1
!     where y=(x-0.5*(xmax+xmin))/(0.5*(xmax-xmin))
!     See Numerical Recipes, pp. 148-150.

      IMPLICIT NONE
      real pi
      integer ncoef, n
      parameter (pi=3.14159265)
      real c(ncoef),f(n)

! local variables      
      real fac, thesum
      integer j, k
      
      fac=2./n
      do j=1,ncoef
         thesum=0
         do k=1,n
            thesum=thesum+f(k)*cos((pi*(j-1))*((k-0.5)/n))
         enddo
         c(j)=fac*thesum
      enddo
      return
      end subroutine chebft             
!*************************************************************************
      subroutine binterp(table,km,im,jm,x,y,xtab,ytab,ix,jy,t,u,out)

!     bilinear interpolation of table
!
      implicit none
      integer im,jm,km
      real table(km,im,jm),xtab(im),ytab(jm),out(km)
      integer i,ix,ip1,j,jy,jp1,k
      real x,dx,t,y,dy,u,tu,  tuc,tcu,tcuc

      if(ix.gt.0)go to 30
      if(im.gt.1)then
        do i=1,im
          if(x.lt.xtab(i))go to 10
        enddo
   10   ix=max0(i-1,1)
        ip1=min0(ix+1,im)
        dx=(xtab(ip1)-xtab(ix))
        if(abs(dx).gt.1.e-20)then
           t=(x-xtab(ix))/(xtab(ix+1)-xtab(ix))
        else
           t=0
        endif
      else
        ix=1
        ip1=1
        t=0
      endif
      if(jm.gt.1)then
        do j=1,jm
          if(y.lt.ytab(j))go to 20
        enddo
   20   jy=max0(j-1,1)
        jp1=min0(jy+1,jm)
        dy=(ytab(jp1)-ytab(jy))
        if(abs(dy).gt.1.e-20)then
           u=(y-ytab(jy))/dy
        else
           u=0
        endif
      else
        jy=1
        jp1=1
        u=0
      endif
   30 continue
      jp1=min(jy+1,jm)
      ip1=min(ix+1,im)
      tu=t*u
      tuc=t-tu
      tcuc=1-tuc-u
      tcu=u-tu
      do k=1,km
         out(k)=tcuc*table(k,ix,jy)+tuc*table(k,ip1,jy)   &
               +tu*table(k,ip1,jp1)+tcu*table(k,ix,jp1)
      enddo
      return
      end subroutine binterp                                            
!***************************************************************
      subroutine  miev0 ( xx, crefin, perfct, mimcut, anyang,   &
                          numang, xmu, nmom, ipolzn, momdim, prnt,   &
                          qext, qsca, gqsc, pmom, sforw, sback, s1,   &
                          s2, tforw, tback )
!
!    computes mie scattering and extinction efficiencies; asymmetry
!    factor;  forward- and backscatter amplitude;  scattering
!    amplitudes for incident polarization parallel and perpendicular
!    to the plane of scattering, as functions of scattering angle;
!    coefficients in the legendre polynomial expansions of either the
!    unpolarized phase function or the polarized phase matrix;
!    and some quantities needed in polarized radiative transfer.
!
!      calls :  biga, ckinmi, small1, small2, testmi, miprnt,
!               lpcoef, errmsg
!
!      i n t e r n a l   v a r i a b l e s
!      -----------------------------------
!
!  an,bn           mie coefficients  little-a-sub-n, little-b-sub-n
!                     ( ref. 1, eq. 16 )
!  anm1,bnm1       mie coefficients  little-a-sub-(n-1),
!                     little-b-sub-(n-1);  used in -gqsc- sum
!  anp             coeffs. in s+ expansion ( ref. 2, p. 1507 )
!  bnp             coeffs. in s- expansion ( ref. 2, p. 1507 )
!  anpm            coeffs. in s+ expansion ( ref. 2, p. 1507 )
!                     when  mu  is replaced by  - mu
!  bnpm            coeffs. in s- expansion ( ref. 2, p. 1507 )
!                     when  mu  is replaced by  - mu
!  calcmo(k)       true, calculate moments for k-th phase quantity
!                     (derived from -ipolzn-; used only in 'lpcoef')
!  cbiga(n)        bessel function ratio capital-a-sub-n (ref. 2, eq. 2)
!                     ( complex version )
!  cior            complex index of refraction with negative
!                     imaginary part (van de hulst convention)
!  cioriv          1 / cior
!  coeff           ( 2n + 1 ) / ( n ( n + 1 ) )
!  fn              floating point version of index in loop performing
!                     mie series summation
!  lita,litb(n)    mie coefficients -an-, -bn-, saved in arrays for
!                     use in calculating legendre moments *pmom*
!  maxtrm          max. possible no. of terms in mie series
!  mm              + 1 and  - 1,  alternately.
!  mim             magnitude of imaginary refractive index
!  mre             real part of refractive index
!  maxang          max. possible value of input variable -numang-
!  nangd2          (numang+1)/2 ( no. of angles in 0-90 deg; anyang=f )
!  noabs           true, sphere non-absorbing (determined by -mimcut-)
!  np1dn           ( n + 1 ) / n
!  npquan          highest-numbered phase quantity for which moments are
!                     to be calculated (the largest digit in -ipolzn-
!                     if  ipolzn .ne. 0)
!  ntrm            no. of terms in mie series
!  pass1           true on first entry, false thereafter; for self-test
!  pin(j)          angular function little-pi-sub-n ( ref. 2, eq. 3 )
!                     at j-th angle
!  pinm1(j)        little-pi-sub-(n-1) ( see -pin- ) at j-th angle
!  psinm1          ricatti-bessel function psi-sub-(n-1), argument -xx-
!  psin            ricatti-bessel function psi-sub-n of argument -xx-
!                     ( ref. 1, p. 11 ff. )
!  rbiga(n)        bessel function ratio capital-a-sub-n (ref. 2, eq. 2)
!                     ( real version, for when imag refrac index = 0 )
!  rioriv          1 / mre
!  rn              1 / n
!  rtmp            (real) temporary variable
!  sp(j)           s+  for j-th angle  ( ref. 2, p. 1507 )
!  sm(j)           s-  for j-th angle  ( ref. 2, p. 1507 )
!  sps(j)          s+  for (numang+1-j)-th angle ( anyang=false )
!  sms(j)          s-  for (numang+1-j)-th angle ( anyang=false )
!  taun            angular function little-tau-sub-n ( ref. 2, eq. 4 )
!                     at j-th angle
!  tcoef           n ( n+1 ) ( 2n+1 ) (for summing tforw,tback series)
!  twonp1          2n + 1
!  yesang          true if scattering amplitudes are to be calculated
!  zetnm1          ricatti-bessel function  zeta-sub-(n-1) of argument
!                     -xx-  ( ref. 2, eq. 17 )
!  zetn            ricatti-bessel function  zeta-sub-n of argument -xx-
!
! ----------------------------------------------------------------------
! --------  i / o specifications for subroutine miev0  -----------------
! ----------------------------------------------------------------------
      implicit none
      logical  anyang, perfct, prnt(*)
      integer  ipolzn, momdim, numang, nmom
      real*8     gqsc, mimcut, pmom( 0:momdim, * ), qext, qsca,   &
               xmu(*), xx
      complex*16  crefin, sforw, sback, s1(*), s2(*), tforw(*),   &
               tback(*)
      integer maxang,mxang2,maxtrm
      real*8 onethr
! ----------------------------------------------------------------------
!
      parameter ( maxang = 501, mxang2 = maxang/2 + 1 )
!
!                                  ** note --  maxtrm = 10100  is neces-
!                                  ** sary to do some of the test probs,
!                                  ** but 1100 is sufficient for most
!                                  ** conceivable applications
      parameter ( maxtrm = 1100 )
      parameter ( onethr = 1./3. )
!
      logical   anysav, calcmo(4), noabs, ok, persav, yesang
      integer   npquan
      integer i,j,n,nmosav,iposav,numsav,ntrm,nangd2
      real*8      mim, mimsav, mre, mm, np1dn
      real*8 rioriv,xmusav,xxsav,sq,fn,rn,twonp1,tcoef, coeff
      real*8 xinv,psinm1,chinm1,psin,chin,rtmp,taun
      real*8      rbiga( maxtrm ), pin( maxang ), pinm1( maxang )
      complex*16   an, bn, anm1, bnm1, anp, bnp, anpm, bnpm, cresav,   &
                cior, cioriv, ctmp, zet, zetnm1, zetn
      complex*16   cbiga( maxtrm ), lita( maxtrm ), litb( maxtrm ),   &
                sp( maxang ), sm( maxang ), sps( mxang2 ), sms( mxang2 )
      equivalence  ( cbiga, rbiga )
      logical, save :: pass1
      data  pass1 / .true. /
      sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
!
!
      if ( pass1 )  then
!                                   ** save certain user input values
         xxsav  = xx
         cresav = crefin
         mimsav = mimcut
         persav = perfct
         anysav = anyang
         nmosav = nmom
         iposav = ipolzn
         numsav = numang
         xmusav = xmu( 1 )
!                              ** reset input values for test case
         xx      = 10.0
         crefin  = ( 1.5, - 0.1 )
         perfct  = .false.
         mimcut  = 0.0
         anyang  = .true.
         numang  = 1
         xmu( 1 )= - 0.7660444
         nmom    = 1
         ipolzn  = - 1
!
      end if
!                                        ** check input and calculate
!                                        ** certain variables from input
!
   10 call  ckinmi( numang, maxang, xx, perfct, crefin, momdim,   &
                    nmom, ipolzn, anyang, xmu, calcmo, npquan )
!
      if ( perfct .and. xx .le. 0.1 )  then
!                                            ** use totally-reflecting
!                                            ** small-particle limit
!
         call  small1 ( xx, numang, xmu, qext, qsca, gqsc, sforw,   &
                        sback, s1, s2, tforw, tback, lita, litb )
         ntrm = 2
         go to 200
      end if
!
      if ( .not.perfct )  then
!
         cior = crefin
         if ( dimag( cior ) .gt. 0.0 )  cior = dconjg( cior )
         mre =     dble( cior )
         mim =  - dimag( cior )
         noabs = mim .le. mimcut
         cioriv = 1.0 / cior
         rioriv = 1.0 / mre
!
         if ( xx * dmax1( 1.d0, cdabs(cior) ) .le. 0.d1 ) then
!
!                                    ** use general-refractive-index
!                                    ** small-particle limit
!                                    ** ( ref. 2, p. 1508 )
!
            call  small2 ( xx, cior, .not.noabs, numang, xmu, qext,   &
                           qsca, gqsc, sforw, sback, s1, s2, tforw,   &
                           tback, lita, litb )
            ntrm = 2
            go to 200
         end if
!
      end if
!
      nangd2 = ( numang + 1 ) / 2
      yesang = numang .gt. 0
!                              ** estimate number of terms in mie series
!                              ** ( ref. 2, p. 1508 )
      if ( xx.le.8.0 )  then
         ntrm = xx + 4. * xx**onethr + 1.
      else if ( xx.lt.4200. )  then
         ntrm = xx + 4.05 * xx**onethr + 2.
      else
         ntrm = xx + 4. * xx**onethr + 2.
      end if
      if ( ntrm+1 .gt. maxtrm )   &
           call errmsg( 'miev0--parameter maxtrm too small', .true. )
!
!                            ** calculate logarithmic derivatives of
!                            ** j-bessel-fcn., big-a-sub-(1 to ntrm)
      if ( .not.perfct )   &
           call  biga( cior, xx, ntrm, noabs, yesang, rbiga, cbiga )
!
!                            ** initialize ricatti-bessel functions
!                            ** (psi,chi,zeta)-sub-(0,1) for upward
!                            ** recurrence ( ref. 1, eq. 19 )
      xinv = 1.0 / xx
      psinm1   = dsin( xx )
      chinm1   = dcos( xx )
      psin = psinm1 * xinv - chinm1
      chin = chinm1 * xinv + psinm1
      zetnm1 = dcmplx( psinm1, chinm1 )
      zetn   = dcmplx( psin, chin )
!                                     ** initialize previous coeffi-
!                                     ** cients for -gqsc- series
      anm1 = ( 0.0, 0.0 )
      bnm1 = ( 0.0, 0.0 )
!                             ** initialize angular function little-pi
!                             ** and sums for s+, s- ( ref. 2, p. 1507 )
      if ( anyang )  then
         do  60  j = 1, numang
            pinm1( j ) = 0.0
            pin( j )   = 1.0
            sp ( j ) = ( 0.0, 0.0 )
            sm ( j ) = ( 0.0, 0.0 )
   60    continue
      else
         do  70  j = 1, nangd2
            pinm1( j ) = 0.0
            pin( j )   = 1.0
            sp ( j ) = ( 0.0, 0.0 )
            sm ( j ) = ( 0.0, 0.0 )
            sps( j ) = ( 0.0, 0.0 )
            sms( j ) = ( 0.0, 0.0 )
   70    continue
      end if
!                         ** initialize mie sums for efficiencies, etc.
      qsca = 0.0
      gqsc = 0.0
      sforw      = ( 0., 0. )
      sback      = ( 0., 0. )
      tforw( 1 ) = ( 0., 0. )
      tback( 1 ) = ( 0., 0. )
!
!
! ---------  loop to sum mie series  -----------------------------------
!
      mm = + 1.0
      do  100  n = 1, ntrm
!                           ** compute various numerical coefficients
         fn     = n
         rn     = 1.0 / fn
         np1dn  = 1.0 + rn
         twonp1 = 2 * n + 1
         coeff  = twonp1 / ( fn * ( n + 1 ) )
         tcoef  = twonp1 * ( fn * ( n + 1 ) )
!
!                              ** calculate mie series coefficients
         if ( perfct )  then
!                                   ** totally-reflecting case
!
            an = ( ( fn*xinv ) * psin - psinm1 ) /   &
                 ( ( fn*xinv ) * zetn - zetnm1 )
            bn = psin / zetn
!
         else if ( noabs )  then
!                                      ** no-absorption case
!
            an =  ( ( rioriv*rbiga(n) + ( fn*xinv ) ) * psin - psinm1 )   &
                / ( ( rioriv*rbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
            bn =  ( (  mre * rbiga(n) + ( fn*xinv ) ) * psin - psinm1 )   &
                / ( (  mre * rbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
         else
!                                       ** absorptive case
!
            an = ( ( cioriv * cbiga(n) + ( fn*xinv ) ) * psin - psinm1 )   &
                /( ( cioriv * cbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
            bn = ( (   cior * cbiga(n) + ( fn*xinv ) ) * psin - psinm1 )   &
                /( (   cior * cbiga(n) + ( fn*xinv ) ) * zetn - zetnm1 )
            qsca = qsca + twonp1 * ( sq( an ) + sq( bn ) )
!
         end if
!                       ** save mie coefficients for *pmom* calculation
         lita( n ) = an
         litb( n ) = bn
!                            ** increment mie sums for non-angle-
!                            ** dependent quantities
!
         sforw      = sforw      + twonp1 * ( an + bn )
         tforw( 1 ) = tforw( 1 ) + tcoef  * ( an - bn )
         sback      = sback      + ( mm * twonp1 ) * ( an - bn )
         tback( 1 ) = tback( 1 ) + ( mm * tcoef )  * ( an + bn )
         gqsc = gqsc + ( fn - rn ) * dble( anm1 * dconjg( an )   &
                                         + bnm1 * dconjg( bn ) )   &
                + coeff * dble( an * dconjg( bn ) )
!
         if ( yesang )  then
!                                      ** put mie coefficients in form
!                                      ** needed for computing s+, s-
!                                      ** ( ref. 2, p. 1507 )
            anp = coeff * ( an + bn )
            bnp = coeff * ( an - bn )
!                                      ** increment mie sums for s+, s-
!                                      ** while upward recursing
!                                      ** angular functions little pi
!                                      ** and little tau
            if ( anyang )  then
!                                         ** arbitrary angles
!
!                                              ** vectorizable loop
               do  80  j = 1, numang
                  rtmp = ( xmu( j ) * pin( j ) ) - pinm1( j )
                  taun =  fn * rtmp - pinm1( j )
                  sp( j )  = sp( j ) + anp * ( pin( j ) + taun )
                  sm( j )  = sm( j ) + bnp * ( pin( j ) - taun )
                  pinm1( j ) = pin( j )
                  pin( j ) = ( xmu( j ) * pin( j ) ) + np1dn * rtmp
   80          continue
!
            else
!                                  ** angles symmetric about 90 degrees
               anpm = mm * anp
               bnpm = mm * bnp
!                                          ** vectorizable loop
               do  90  j = 1, nangd2
                  rtmp = ( xmu( j ) * pin( j ) ) - pinm1( j )
                  taun =  fn * rtmp - pinm1( j )
                  sp ( j ) = sp ( j ) +  anp * ( pin( j ) + taun )
                  sms( j ) = sms( j ) + bnpm * ( pin( j ) + taun )
                  sm ( j ) = sm ( j ) +  bnp * ( pin( j ) - taun )
                  sps( j ) = sps( j ) + anpm * ( pin( j ) - taun )
                  pinm1( j ) = pin( j )
                  pin( j ) = ( xmu( j ) * pin( j ) ) + np1dn * rtmp
   90          continue
!
            end if
         end if
!                          ** update relevant quantities for next
!                          ** pass through loop
         mm   =  - mm
         anm1 = an
         bnm1 = bn
!                           ** upward recurrence for ricatti-bessel
!                           ** functions ( ref. 1, eq. 17 )
!
         zet    = ( twonp1 * xinv ) * zetn - zetnm1
         zetnm1 = zetn
         zetn   = zet
         psinm1 = psin
         psin   = dble( zetn )
  100 continue
!
! ---------- end loop to sum mie series --------------------------------
!
!
      qext = 2. / xx**2 * dble( sforw )
      if ( perfct .or. noabs )  then
         qsca = qext
      else
         qsca = 2. / xx**2 * qsca
      end if
!
      gqsc = 4. / xx**2 * gqsc
      sforw = 0.5 * sforw
      sback = 0.5 * sback
      tforw( 2 ) = 0.5 * (   sforw + 0.25 * tforw( 1 ) )
      tforw( 1 ) = 0.5 * (   sforw - 0.25 * tforw( 1 ) )
      tback( 2 ) = 0.5 * (   sback + 0.25 * tback( 1 ) )
      tback( 1 ) = 0.5 * ( - sback + 0.25 * tback( 1 ) )
!
      if ( yesang )  then
!                                ** recover scattering amplitudes
!                                ** from s+, s- ( ref. 1, eq. 11 )
         if ( anyang )  then
!                                         ** vectorizable loop
            do  110  j = 1, numang
               s1( j ) = 0.5 * ( sp( j ) + sm( j ) )
               s2( j ) = 0.5 * ( sp( j ) - sm( j ) )
  110       continue
!
         else
!                                         ** vectorizable loop
            do  120  j = 1, nangd2
               s1( j ) = 0.5 * ( sp( j ) + sm( j ) )
               s2( j ) = 0.5 * ( sp( j ) - sm( j ) )
  120       continue
!                                         ** vectorizable loop
            do  130  j = 1, nangd2
               s1( numang+1 - j ) = 0.5 * ( sps( j ) + sms( j ) )
               s2( numang+1 - j ) = 0.5 * ( sps( j ) - sms( j ) )
  130       continue
         end if
!
      end if
!                                         ** calculate legendre moments
  200 if ( nmom.gt.0 )   &
           call lpcoef ( ntrm, nmom, ipolzn, momdim, calcmo, npquan,   &
                         lita, litb, pmom )
!
      if ( dimag(crefin) .gt. 0.0 )  then
!                                         ** take complex conjugates
!                                         ** of scattering amplitudes
         sforw = dconjg( sforw )
         sback = dconjg( sback )
         do  210  i = 1, 2
            tforw( i ) = dconjg( tforw(i) )
            tback( i ) = dconjg( tback(i) )
  210    continue
!
         do  220  j = 1, numang
            s1( j ) = dconjg( s1(j) )
            s2( j ) = dconjg( s2(j) )
  220    continue
!
      end if
!
      if ( pass1 )  then
!                             ** compare test case results with
!                             ** correct answers and abort if bad
!
         call  testmi ( qext, qsca, gqsc, sforw, sback, s1, s2,   &
                        tforw, tback, pmom, momdim, ok )
         if ( .not. ok )  then
            prnt(1) = .false.
            prnt(2) = .false.
            call  miprnt( prnt, xx, perfct, crefin, numang, xmu, qext,   &
                          qsca, gqsc, nmom, ipolzn, momdim, calcmo,   &
                          pmom, sforw, sback, tforw, tback, s1, s2 )
            call errmsg( 'miev0 -- self-test failed', .true. )
         end if
!                                       ** restore user input values
         xx     = xxsav
         crefin = cresav
         mimcut = mimsav
         perfct = persav
         anyang = anysav
         nmom   = nmosav
         ipolzn = iposav
         numang = numsav
         xmu(1) = xmusav
         pass1 = .false.
         go to 10
!
      end if
!
      if ( prnt(1) .or. prnt(2) )   &
         call  miprnt( prnt, xx, perfct, crefin, numang, xmu, qext,   &
                       qsca, gqsc, nmom, ipolzn, momdim, calcmo,   &
                       pmom, sforw, sback, tforw, tback, s1, s2 )
!
      return
!
      end subroutine  miev0 
!****************************************************************************                                           
      subroutine  ckinmi( numang, maxang, xx, perfct, crefin, momdim,   &
                          nmom, ipolzn, anyang, xmu, calcmo, npquan )
!
!        check for bad input to 'miev0' and calculate -calcmo,npquan-
!
      implicit none
      logical  perfct, anyang, calcmo(*)
      integer  numang, maxang, momdim, nmom, ipolzn, npquan
      real*8    xx, xmu(*)
      integer i,l,j,ip
      complex*16  crefin
!
      character*4  string
      logical  inperr
!
      inperr = .false.
!
      if ( numang.gt.maxang )  then
         call errmsg( 'miev0--parameter maxang too small', .true. )
         inperr = .true.
      end if
      if ( numang.lt.0 )  call  wrtbad( 'numang', inperr )
      if ( xx.lt.0. )     call  wrtbad( 'xx', inperr )
      if ( .not.perfct .and. dble(crefin).le.0. )   &
           call wrtbad( 'crefin', inperr )
      if ( momdim.lt.1 )  call wrtbad( 'momdim', inperr )
!
      if ( nmom.ne.0 )  then
         if ( nmom.lt.0 .or. nmom.gt.momdim ) call wrtbad('nmom',inperr)
         if ( iabs(ipolzn).gt.4444 )  call  wrtbad( 'ipolzn', inperr )
         npquan = 0
         do 5  l = 1, 4
            calcmo( l ) = .false.
    5    continue
         if ( ipolzn.ne.0 )  then
!                                 ** parse out -ipolzn- into its digits
!                                 ** to find which phase quantities are
!                                 ** to have their moments calculated
!
            write( string, '(i4)' )  iabs(ipolzn)
            do 10  j = 1, 4
               ip = ichar( string(j:j) ) - ichar( '0' )
               if ( ip.ge.1 .and. ip.le.4 )  calcmo( ip ) = .true.
               if ( ip.eq.0 .or. (ip.ge.5 .and. ip.le.9) )   &
                    call  wrtbad( 'ipolzn', inperr )
               npquan = max0( npquan, ip )
   10       continue
         end if
      end if
!
      if ( anyang )  then
!                                ** allow for slight imperfections in
!                                ** computation of cosine
          do  20  i = 1, numang
             if ( xmu(i).lt.-1.00001 .or. xmu(i).gt.1.00001 )   &
                  call wrtbad( 'xmu', inperr )
   20     continue
      else
          do  22  i = 1, ( numang + 1 ) / 2
             if ( xmu(i).lt.-0.00001 .or. xmu(i).gt.1.00001 )   &
                  call wrtbad( 'xmu', inperr )
   22     continue
      end if
!
      if ( inperr )   &
           call errmsg( 'miev0--input error(s).  aborting...', .true. )
!
      if ( xx.gt.20000.0 .or. dble(crefin).gt.10.0 .or.   &
           dabs( dimag(crefin) ).gt.10.0 )  call  errmsg(   &
           'miev0--xx or crefin outside tested range', .false. )
!
      return
      end subroutine  ckinmi
!***********************************************************************                                                   
      subroutine  lpcoef ( ntrm, nmom, ipolzn, momdim, calcmo, npquan,   &
                           a, b, pmom )
!
!         calculate legendre polynomial expansion coefficients (also
!         called moments) for phase quantities ( ref. 5 formulation )
!
!     input:  ntrm                    number terms in mie series
!             nmom, ipolzn, momdim    'miev0' arguments
!             calcmo                  flags calculated from -ipolzn-
!             npquan                  defined in 'miev0'
!             a, b                    mie series coefficients
!
!     output: pmom                   legendre moments ('miev0' argument)
!
!     *** notes ***
!
!         (1)  eqs. 2-5 are in error in dave, appl. opt. 9,
!         1888 (1970).  eq. 2 refers to m1, not m2;  eq. 3 refers to
!         m2, not m1.  in eqs. 4 and 5, the subscripts on the second
!         term in square brackets should be interchanged.
!
!         (2)  the general-case logic in this subroutine works correctly
!         in the two-term mie series case, but subroutine  'lpco2t'
!         is called instead, for speed.
!
!         (3)  subroutine  'lpco1t', to do the one-term case, is never
!         called within the context of 'miev0', but is included for
!         complete generality.
!
!         (4)  some improvement in speed is obtainable by combining the
!         310- and 410-loops, if moments for both the third and fourth
!         phase quantities are desired, because the third phase quantity
!         is the real part of a complex series, while the fourth phase
!         quantity is the imaginary part of that very same series.  but
!         most users are not interested in the fourth phase quantity,
!         which is related to circular polarization, so the present
!         scheme is usually more efficient.
!
      implicit none
      logical  calcmo(*)
      integer  ipolzn, momdim, nmom, ntrm, npquan
      real*8    pmom( 0:momdim, * )
      complex*16  a(*), b(*)
!
!           ** specification of local variables
!
!      am(m)       numerical coefficients  a-sub-m-super-l
!                     in dave, eqs. 1-15, as simplified in ref. 5.
!
!      bi(i)       numerical coefficients  b-sub-i-super-l
!                     in dave, eqs. 1-15, as simplified in ref. 5.
!
!      bidel(i)    1/2 bi(i) times factor capital-del in dave
!
!      cm,dm()     arrays c and d in dave, eqs. 16-17 (mueller form),
!                     calculated using recurrence derived in ref. 5
!
!      cs,ds()     arrays c and d in ref. 4, eqs. a5-a6 (sekera form),
!                     calculated using recurrence derived in ref. 5
!
!      c,d()       either -cm,dm- or -cs,ds-, depending on -ipolzn-
!
!      evenl       true for even-numbered moments;  false otherwise
!
!      idel        1 + little-del  in dave
!
!      maxtrm      max. no. of terms in mie series
!
!      maxmom      max. no. of non-zero moments
!
!      nummom      number of non-zero moments
!
!      recip(k)    1 / k
!
      integer maxtrm,maxmom,mxmom2,maxrcp
      parameter  ( maxtrm = 1102, maxmom = 2*maxtrm, mxmom2 = maxmom/2,   &
                   maxrcp = 4*maxtrm + 2 )
      real*8      am( 0:maxtrm ), bi( 0:mxmom2 ), bidel( 0:mxmom2 )
      real*8, save :: recip( maxrcp )
      complex*16 cm( maxtrm ), dm( maxtrm ), cs( maxtrm ), ds( maxtrm ),   &
                 c( maxtrm ), d( maxtrm )
      integer k,j,l,nummom,ld2,idel,m,i,mmax,imax
      real*8 thesum
      equivalence  ( c, cm ),  ( d, dm )
      logical evenl
      logical, save :: pass1
      data  pass1 / .true. /
!
!
      if ( pass1 )  then
!
         do  1  k = 1, maxrcp
            recip( k ) = 1.0 / k
    1    continue
         pass1 = .false.
!
      end if
!
      do  5  j = 1, max0( 1, npquan )
         do  5  l = 0, nmom
            pmom( l, j ) = 0.0
    5 continue
!
      if ( ntrm.eq.1 )  then
         call  lpco1t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
         return
      else if ( ntrm.eq.2 )  then
         call  lpco2t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
         return
      end if
!
      if ( ntrm+2 .gt. maxtrm )   &
           call errmsg( 'lpcoef--parameter maxtrm too small', .true. )
!
!                                     ** calculate mueller c, d arrays
      cm( ntrm+2 ) = ( 0., 0. )
      dm( ntrm+2 ) = ( 0., 0. )
      cm( ntrm+1 ) = ( 1. - recip( ntrm+1 ) ) * b( ntrm )
      dm( ntrm+1 ) = ( 1. - recip( ntrm+1 ) ) * a( ntrm )
      cm( ntrm ) = ( recip(ntrm) + recip(ntrm+1) ) * a( ntrm )   &
                   + ( 1. - recip(ntrm) ) * b( ntrm-1 )
      dm( ntrm ) = ( recip(ntrm) + recip(ntrm+1) ) * b( ntrm )   &
                   + ( 1. - recip(ntrm) ) * a( ntrm-1 )
!
      do  10  k = ntrm-1, 2, -1
         cm( k ) = cm( k+2 ) - ( 1. + recip(k+1) ) * b( k+1 )   &
                             + ( recip(k) + recip(k+1) ) * a( k )   &
                             + ( 1. - recip(k) ) * b( k-1 )
         dm( k ) = dm( k+2 ) - ( 1. + recip(k+1) ) * a( k+1 )   &
                             + ( recip(k) + recip(k+1) ) * b( k )   &
                             + ( 1. - recip(k) ) * a( k-1 )
   10 continue
      cm( 1 ) = cm( 3 ) + 1.5 * ( a( 1 ) - b( 2 ) )
      dm( 1 ) = dm( 3 ) + 1.5 * ( b( 1 ) - a( 2 ) )
!
      if ( ipolzn.ge.0 )  then
!
         do  20  k = 1, ntrm + 2
            c( k ) = ( 2*k - 1 ) * cm( k )
            d( k ) = ( 2*k - 1 ) * dm( k )
   20    continue
!
      else
!                                    ** compute sekera c and d arrays
         cs( ntrm+2 ) = ( 0., 0. )
         ds( ntrm+2 ) = ( 0., 0. )
         cs( ntrm+1 ) = ( 0., 0. )
         ds( ntrm+1 ) = ( 0., 0. )
!
         do  30  k = ntrm, 1, -1
            cs( k ) = cs( k+2 ) + ( 2*k + 1 ) * ( cm( k+1 ) - b( k ) )
            ds( k ) = ds( k+2 ) + ( 2*k + 1 ) * ( dm( k+1 ) - a( k ) )
   30    continue
!
         do  40  k = 1, ntrm + 2
            c( k ) = ( 2*k - 1 ) * cs( k )
            d( k ) = ( 2*k - 1 ) * ds( k )
   40    continue
!
      end if
!
!
      if( ipolzn.lt.0 )  nummom = min0( nmom, 2*ntrm - 2 )
      if( ipolzn.ge.0 )  nummom = min0( nmom, 2*ntrm )
      if ( nummom .gt. maxmom )   &
           call errmsg( 'lpcoef--parameter maxtrm too small', .true. )
!
!                               ** loop over moments
      do  500  l = 0, nummom
         ld2 = l / 2
         evenl = mod( l,2 ) .eq. 0
!                                    ** calculate numerical coefficients
!                                    ** a-sub-m and b-sub-i in dave
!                                    ** double-sums for moments
         if( l.eq.0 )  then
!
            idel = 1
            do  60  m = 0, ntrm
               am( m ) = 2.0 * recip( 2*m + 1 )
   60       continue
            bi( 0 ) = 1.0
!
         else if( evenl )  then
!
            idel = 1
            do  70  m = ld2, ntrm
               am( m ) = ( 1. + recip( 2*m-l+1 ) ) * am( m )
   70       continue
            do  75  i = 0, ld2-1
               bi( i ) = ( 1. - recip( l-2*i ) ) * bi( i )
   75       continue
            bi( ld2 ) = ( 2. - recip( l ) ) * bi( ld2-1 )
!
         else
!
            idel = 2
            do  80  m = ld2, ntrm
               am( m ) = ( 1. - recip( 2*m+l+2 ) ) * am( m )
   80       continue
            do  85  i = 0, ld2
               bi( i ) = ( 1. - recip( l+2*i+1 ) ) * bi( i )
   85       continue
!
         end if
!                                     ** establish upper limits for sums
!                                     ** and incorporate factor capital-
!                                     ** del into b-sub-i
         mmax = ntrm - idel
         if( ipolzn.ge.0 )  mmax = mmax + 1
         imax = min0( ld2, mmax - ld2 )
         if( imax.lt.0 )  go to 600
         do  90  i = 0, imax
            bidel( i ) = bi( i )
   90    continue
         if( evenl )  bidel( 0 ) = 0.5 * bidel( 0 )
!
!                                    ** perform double sums just for
!                                    ** phase quantities desired by user
         if( ipolzn.eq.0 )  then
!
            do  110  i = 0, imax
!                                           ** vectorizable loop (cray)
               thesum = 0.0
               do  100  m = ld2, mmax - i
                  thesum = thesum + am( m ) *   &
                            ( dble( c(m-i+1) * dconjg( c(m+i+idel) ) )   &
                            + dble( d(m-i+1) * dconjg( d(m+i+idel) ) ) )
  100          continue
               pmom( l,1 ) = pmom( l,1 ) + bidel( i ) * thesum
  110       continue
            pmom( l,1 ) = 0.5 * pmom( l,1 )
            go to 500
!
         end if
!
         if ( calcmo(1) )  then
            do  160  i = 0, imax
!                                           ** vectorizable loop (cray)
               thesum = 0.0
               do  150  m = ld2, mmax - i
                  thesum = thesum + am( m ) *   &
                              dble( c(m-i+1) * dconjg( c(m+i+idel) ) )
  150          continue
               pmom( l,1 ) = pmom( l,1 ) + bidel( i ) * thesum
  160       continue
         end if
!
!
         if ( calcmo(2) )  then
            do  210  i = 0, imax
!                                           ** vectorizable loop (cray)
               thesum = 0.0
               do  200  m = ld2, mmax - i
                  thesum = thesum + am( m ) *   &
                              dble( d(m-i+1) * dconjg( d(m+i+idel) ) )
  200          continue
               pmom( l,2 ) = pmom( l,2 ) + bidel( i ) * thesum
  210       continue
         end if
!
!
         if ( calcmo(3) )  then
            do  310  i = 0, imax
!                                           ** vectorizable loop (cray)
               thesum = 0.0
               do  300  m = ld2, mmax - i
                  thesum = thesum + am( m ) *   &
                            ( dble( c(m-i+1) * dconjg( d(m+i+idel) ) )   &
                            + dble( c(m+i+idel) * dconjg( d(m-i+1) ) ) )
  300          continue
               pmom( l,3 ) = pmom( l,3 ) + bidel( i ) * thesum
  310       continue
            pmom( l,3 ) = 0.5 * pmom( l,3 )
         end if
!
!
         if ( calcmo(4) )  then
            do  410  i = 0, imax
!                                           ** vectorizable loop (cray)
               thesum = 0.0
               do  400  m = ld2, mmax - i
                  thesum = thesum + am( m ) *   &
                            ( dimag( c(m-i+1) * dconjg( d(m+i+idel) ) )   &
                            + dimag( c(m+i+idel) * dconjg( d(m-i+1) ) ))
  400          continue
               pmom( l,4 ) = pmom( l,4 ) + bidel( i ) * thesum
  410       continue
            pmom( l,4 ) = - 0.5 * pmom( l,4 )
         end if
!
  500 continue
!
!
  600 return
      end subroutine  lpcoef
!*********************************************************************                                                    
      subroutine  lpco1t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
!
!         calculate legendre polynomial expansion coefficients (also
!         called moments) for phase quantities in special case where
!         no. terms in mie series = 1
!
!        input:  nmom, ipolzn, momdim     'miev0' arguments
!                calcmo                   flags calculated from -ipolzn-
!                a(1), b(1)               mie series coefficients
!
!        output: pmom                     legendre moments
!
      implicit none
      logical  calcmo(*)
      integer  ipolzn, momdim, nmom,nummom,l
      real*8    pmom( 0:momdim, * ),sq,a1sq,b1sq
      complex*16  a(*), b(*), ctmp, a1b1c
      sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
!
!
      a1sq = sq( a(1) )
      b1sq = sq( b(1) )
      a1b1c = a(1) * dconjg( b(1) )
!
      if( ipolzn.lt.0 )  then
!
         if( calcmo(1) )  pmom( 0,1 ) = 2.25 * b1sq
         if( calcmo(2) )  pmom( 0,2 ) = 2.25 * a1sq
         if( calcmo(3) )  pmom( 0,3 ) = 2.25 * dble( a1b1c )
         if( calcmo(4) )  pmom( 0,4 ) = 2.25 *dimag( a1b1c )
!
      else
!
         nummom = min0( nmom, 2 )
!                                   ** loop over moments
         do  100  l = 0, nummom
!
            if( ipolzn.eq.0 )  then
               if( l.eq.0 )  pmom( l,1 ) = 1.5 * ( a1sq + b1sq )
               if( l.eq.1 )  pmom( l,1 ) = 1.5 * dble( a1b1c )
               if( l.eq.2 )  pmom( l,1 ) = 0.15 * ( a1sq + b1sq )
               go to 100
            end if
!
            if( calcmo(1) )  then
               if( l.eq.0 )  pmom( l,1 ) = 2.25 * ( a1sq + b1sq / 3. )
               if( l.eq.1 )  pmom( l,1 ) = 1.5 * dble( a1b1c )
               if( l.eq.2 )  pmom( l,1 ) = 0.3 * b1sq
            end if
!
            if( calcmo(2) )  then
               if( l.eq.0 )  pmom( l,2 ) = 2.25 * ( b1sq + a1sq / 3. )
               if( l.eq.1 )  pmom( l,2 ) = 1.5 * dble( a1b1c )
               if( l.eq.2 )  pmom( l,2 ) = 0.3 * a1sq
            end if
!
            if( calcmo(3) )  then
               if( l.eq.0 )  pmom( l,3 ) = 3.0 * dble( a1b1c )
               if( l.eq.1 )  pmom( l,3 ) = 0.75 * ( a1sq + b1sq )
               if( l.eq.2 )  pmom( l,3 ) = 0.3 * dble( a1b1c )
            end if
!
            if( calcmo(4) )  then
               if( l.eq.0 )  pmom( l,4 ) = - 1.5 * dimag( a1b1c )
               if( l.eq.1 )  pmom( l,4 ) = 0.0
               if( l.eq.2 )  pmom( l,4 ) = 0.3 * dimag( a1b1c )
            end if
!
  100    continue
!
      end if
!
      return
      end subroutine  lpco1t 
!********************************************************************                                                  
      subroutine  lpco2t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
!
!         calculate legendre polynomial expansion coefficients (also
!         called moments) for phase quantities in special case where
!         no. terms in mie series = 2
!
!        input:  nmom, ipolzn, momdim     'miev0' arguments
!                calcmo                   flags calculated from -ipolzn-
!                a(1-2), b(1-2)           mie series coefficients
!
!        output: pmom                     legendre moments
!
      implicit none
      logical  calcmo(*)
      integer  ipolzn, momdim, nmom,l,nummom
      real*8    pmom( 0:momdim, * ),sq,pm1,pm2,a2sq,b2sq
      complex*16  a(*), b(*)
      complex*16  a2c, b2c, ctmp, ca, cac, cat, cb, cbc, cbt, cg, ch
      sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
!
!
      ca = 3. * a(1) - 5. * b(2)
      cat= 3. * b(1) - 5. * a(2)
      cac = dconjg( ca )
      a2sq = sq( a(2) )
      b2sq = sq( b(2) )
      a2c = dconjg( a(2) )
      b2c = dconjg( b(2) )
!
      if( ipolzn.lt.0 )  then
!                                   ** loop over sekera moments
         nummom = min0( nmom, 2 )
         do  50  l = 0, nummom
!
            if( calcmo(1) )  then
               if( l.eq.0 ) pmom( l,1 ) = 0.25 * ( sq(cat) +   &
                                                   (100./3.) * b2sq )
               if( l.eq.1 ) pmom( l,1 ) = (5./3.) * dble( cat * b2c )
               if( l.eq.2 ) pmom( l,1 ) = (10./3.) * b2sq
            end if
!
            if( calcmo(2) )  then
               if( l.eq.0 ) pmom( l,2 ) = 0.25 * ( sq(ca) +   &
                                                   (100./3.) * a2sq )
               if( l.eq.1 ) pmom( l,2 ) = (5./3.) * dble( ca * a2c )
               if( l.eq.2 ) pmom( l,2 ) = (10./3.) * a2sq
            end if
!
            if( calcmo(3) )  then
               if( l.eq.0 ) pmom( l,3 ) = 0.25 * dble( cat*cac +   &
                                                 (100./3.)*b(2)*a2c )
               if( l.eq.1 ) pmom( l,3 ) = 5./6. * dble( b(2)*cac +   &
                                                        cat*a2c )
               if( l.eq.2 ) pmom( l,3 ) = 10./3. * dble( b(2) * a2c )
            end if
!
            if( calcmo(4) )  then
               if( l.eq.0 ) pmom( l,4 ) = -0.25 * dimag( cat*cac +   &
                                                 (100./3.)*b(2)*a2c )
               if( l.eq.1 ) pmom( l,4 ) = -5./6. * dimag( b(2)*cac +   &
                                                        cat*a2c )
               if( l.eq.2 ) pmom( l,4 ) = -10./3. * dimag( b(2) * a2c )
            end if
!
   50    continue
!
      else
!
         cb = 3. * b(1) + 5. * a(2)
         cbt= 3. * a(1) + 5. * b(2)
         cbc = dconjg( cb )
         cg = ( cbc*cbt + 10.*( cac*a(2) + b2c*cat) ) / 3.
         ch = 2.*( cbc*a(2) + b2c*cbt )
!
!                                   ** loop over mueller moments
         nummom = min0( nmom, 4 )
         do  100  l = 0, nummom
!
            if( ipolzn.eq.0 .or. calcmo(1) )  then
               if( l.eq.0 ) pm1 = 0.25 * sq(ca) + sq(cb) / 12.   &
                                  + (5./3.) * dble(ca*b2c) + 5.*b2sq
               if( l.eq.1 ) pm1 = dble( cb * ( cac/6. + b2c ) )
               if( l.eq.2 ) pm1 = sq(cb)/30. + (20./7.) * b2sq   &
                                  + (2./3.) * dble( ca * b2c )
               if( l.eq.3 ) pm1 = (2./7.) * dble( cb * b2c )
               if( l.eq.4 ) pm1 = (40./63.) * b2sq
               if ( calcmo(1) )  pmom( l,1 ) = pm1
            end if
!
            if( ipolzn.eq.0 .or. calcmo(2) )  then
               if( l.eq.0 ) pm2 = 0.25*sq(cat) + sq(cbt) / 12.   &
                                  + (5./3.) * dble(cat*a2c) + 5.*a2sq
               if( l.eq.1 ) pm2 = dble( cbt * ( dconjg(cat)/6. + a2c) )
               if( l.eq.2 ) pm2 = sq(cbt)/30. + (20./7.) * a2sq   &
                                  + (2./3.) * dble( cat * a2c )
               if( l.eq.3 ) pm2 = (2./7.) * dble( cbt * a2c )
               if( l.eq.4 ) pm2 = (40./63.) * a2sq
               if ( calcmo(2) )  pmom( l,2 ) = pm2
            end if
!
            if( ipolzn.eq.0 )  then
               pmom( l,1 ) = 0.5 * ( pm1 + pm2 )
               go to 100
            end if
!
            if( calcmo(3) )  then
               if( l.eq.0 ) pmom( l,3 ) = 0.25 * dble( cac*cat + cg +   &
                                                       20.*b2c*a(2) )
               if( l.eq.1 ) pmom( l,3 ) = dble( cac*cbt + cbc*cat +   &
                                                3.*ch ) / 12.
               if( l.eq.2 ) pmom( l,3 ) = 0.1 * dble( cg + (200./7.) *   &
                                                      b2c * a(2) )
               if( l.eq.3 ) pmom( l,3 ) = dble( ch ) / 14.
               if( l.eq.4 ) pmom( l,3 ) = 40./63. * dble( b2c * a(2) )
            end if
!
            if( calcmo(4) )  then
               if( l.eq.0 ) pmom( l,4 ) = 0.25 * dimag( cac*cat + cg +   &
                                                        20.*b2c*a(2) )
               if( l.eq.1 ) pmom( l,4 ) = dimag( cac*cbt + cbc*cat +   &
                                                 3.*ch ) / 12.
               if( l.eq.2 ) pmom( l,4 ) = 0.1 * dimag( cg + (200./7.) *   &
                                                       b2c * a(2) )
               if( l.eq.3 ) pmom( l,4 ) = dimag( ch ) / 14.
               if( l.eq.4 ) pmom( l,4 ) = 40./63. * dimag( b2c * a(2) )
            end if
!
  100    continue
!
      end if
!
      return
      end subroutine  lpco2t 
!*********************************************************************                                                  
      subroutine  biga( cior, xx, ntrm, noabs, yesang, rbiga, cbiga )
!
!        calculate logarithmic derivatives of j-bessel-function
!
!     input :  cior, xx, ntrm, noabs, yesang  (defined in 'miev0')
!
!    output :  rbiga or cbiga  (defined in 'miev0')
!
!    internal variables :
!
!       confra     value of lentz continued fraction for -cbiga(ntrm)-,
!                     used to initialize downward recurrence.
!       down       = true, use down-recurrence.  false, do not.
!       f1,f2,f3   arithmetic statement functions used in determining
!                     whether to use up-  or down-recurrence
!                     ( ref. 2, eqs. 6-8 )
!       mre        real refractive index
!       mim        imaginary refractive index
!       rezinv     1 / ( mre * xx ); temporary variable for recurrence
!       zinv       1 / ( cior * xx ); temporary variable for recurrence
!
      implicit none
      logical  down, noabs, yesang
      integer  ntrm,n
      real*8    mre, mim, rbiga(*), xx, rezinv, rtmp, f1,f2,f3
!      complex*16  cior, ctmp, confra, cbiga(*), zinv
      complex*16  cior, ctmp,  cbiga(*), zinv
      f1( mre ) =  - 8.0 + mre**2 * ( 26.22 + mre * ( - 0.4474   &
                   + mre**3 * ( 0.00204 - 0.000175 * mre ) ) )
      f2( mre ) = 3.9 + mre * ( - 10.8 + 13.78 * mre )
      f3( mre ) =  - 15.04 + mre * ( 8.42 + 16.35 * mre )
!
!                                  ** decide whether 'biga' can be
!                                  ** calculated by up-recurrence
      mre =  dble( cior )
      mim =  dabs( dimag( cior ) )
      if ( mre.lt.1.0 .or. mre.gt.10.0 .or. mim.gt.10.0 )  then
         down = .true.
      else if ( yesang )  then
         down = .true.
         if ( mim*xx .lt. f2( mre ) )  down = .false.
      else
         down = .true.
         if ( mim*xx .lt. f1( mre ) )  down = .false.
      end if
!
      zinv  = 1.0 / ( cior * xx )
      rezinv = 1.0 / ( mre * xx )
!
      if ( down )  then
!                          ** compute initial high-order 'biga' using
!                          ** lentz method ( ref. 1, pp. 17-20 )
!
         ctmp = confra( ntrm, zinv, xx )
!
!                                   *** downward recurrence for 'biga'
!                                   *** ( ref. 1, eq. 22 )
         if ( noabs )  then
!                                            ** no-absorption case
            rbiga( ntrm ) = dble( ctmp )
            do  25  n = ntrm, 2, - 1
               rbiga( n-1 ) = (n*rezinv)   &
                               - 1.0 / ( (n*rezinv) + rbiga( n ) )
   25       continue
!
         else
!                                            ** absorptive case
            cbiga( ntrm ) = ctmp
            do  30  n = ntrm, 2, - 1
               cbiga( n-1 ) = (n*zinv) - 1.0 / ( (n*zinv) + cbiga( n ) )
   30       continue
!
         end if
!
      else
!                              *** upward recurrence for 'biga'
!                              *** ( ref. 1, eqs. 20-21 )
         if ( noabs )  then
!                                            ** no-absorption case
            rtmp = dsin( mre*xx )
            rbiga( 1 ) =  - rezinv   &
                           + rtmp / ( rtmp*rezinv - dcos( mre*xx ) )
            do  40  n = 2, ntrm
               rbiga( n ) = - ( n*rezinv )   &
                             + 1.0 / ( ( n*rezinv ) - rbiga( n-1 ) )
   40       continue
!
         else
!                                                ** absorptive case
            ctmp = cdexp( - dcmplx(0.d0,2.d0) * cior * xx )
            cbiga( 1 ) = - zinv + (1.-ctmp) / ( zinv * (1.-ctmp) -   &
                           dcmplx(0.d0,1.d0)*(1.+ctmp) )
            do  50  n = 2, ntrm
               cbiga( n ) = - (n*zinv) + 1.0 / ((n*zinv) - cbiga( n-1 ))
   50       continue
         end if
!
      end if
!
      return
      end subroutine  biga  
!**********************************************************************                                                   
      complex*16 function  confra( n, zinv, xx )
!
!         compute bessel function ratio capital-a-sub-n from its
!         continued fraction using lentz method ( ref. 1, pp. 17-20 )
!
!         zinv = reciprocal of argument of capital-a
!
!    i n t e r n a l    v a r i a b l e s
!    ------------------------------------
!
!    cak      term in continued fraction expansion of capital-a
!                ( ref. 1, eq. 25 )
!    capt     factor used in lentz iteration for capital-a
!                ( ref. 1, eq. 27 )
!    cdenom   denominator in -capt-  ( ref. 1, eq. 28b )
!    cnumer   numerator   in -capt-  ( ref. 1, eq. 28a )
!    cdtd     product of two successive denominators of -capt-
!                factors  ( ref. 1, eq. 34c )
!    cntn     product of two successive numerators of -capt-
!                factors  ( ref. 1, eq. 34b )
!    eps1     ill-conditioning criterion
!    eps2     convergence criterion
!    kk       subscript k of -cak-  ( ref. 1, eq. 25b )
!    kount    iteration counter ( used only to prevent runaway )
!    maxit    max. allowed no. of iterations
!    mm       + 1  and - 1, alternately
!
      implicit none
      integer   n,mm,kk,kount
      integer, save :: maxit
      data  maxit / 10000 /
      real*8     xx
      real*8, save :: eps1,eps2
      data  eps1 / 1.d-2 /, eps2 / 1.d-8 /
      complex*16   zinv
      complex*16   cak, capt, cdenom, cdtd, cnumer, cntn
!
!                                      *** ref. 1, eqs. 25a, 27
      confra = ( n + 1 ) * zinv
      mm     =  - 1
      kk     = 2 * n + 3
      cak    = ( mm * kk ) * zinv
      cdenom = cak
      cnumer = cdenom + 1.0 / confra
      kount  = 1
!
   20 kount = kount + 1
      if ( kount.gt.maxit )   &
           call errmsg( 'confra--iteration failed to converge$', .true.)
!
!                                         *** ref. 2, eq. 25b
      mm  =  - mm
      kk  = kk + 2
      cak = ( mm * kk ) * zinv
!                                         *** ref. 2, eq. 32
      if (      cdabs( cnumer/cak ).le.eps1   &
           .or. cdabs( cdenom/cak ).le.eps1 )  then
!
!                                  ** ill-conditioned case -- stride
!                                  ** two terms instead of one
!
!                                         *** ref. 2, eqs. 34
         cntn   = cak * cnumer + 1.0
         cdtd   = cak * cdenom + 1.0
         confra = ( cntn / cdtd ) * confra
!                                             *** ref. 2, eq. 25b
         mm  =  - mm
         kk  = kk + 2
         cak = ( mm * kk ) * zinv
!                                         *** ref. 2, eqs. 35
         cnumer = cak + cnumer / cntn
         cdenom = cak + cdenom / cdtd
         kount  = kount + 1
         go to 20
!
      else
!                                ** well-conditioned case
!
!                                        *** ref. 2, eqs. 26, 27
         capt   = cnumer / cdenom
         confra = capt * confra
!                                    ** check for convergence
!                                    ** ( ref. 2, eq. 31 )
!
         if (      dabs( dble(capt) - 1.0 ).ge.eps2   &
              .or. dabs( dimag(capt) )      .ge.eps2 )  then
!
!                                        *** ref. 2, eqs. 30a-b
            cnumer = cak + 1.0 / cnumer
            cdenom = cak + 1.0 / cdenom
            go to 20
         end if
      end if
!
      return
!
      end function confra
!********************************************************************      
      subroutine  miprnt( prnt, xx, perfct, crefin, numang, xmu,   &
                          qext, qsca, gqsc, nmom, ipolzn, momdim,   &
                          calcmo, pmom, sforw, sback, tforw, tback,   &
                          s1, s2 )
!
!         print scattering quantities of a single particle
!
      implicit none
      logical  perfct, prnt(*), calcmo(*)
      integer  ipolzn, momdim, nmom, numang,i,m,j
      real*8    gqsc, pmom( 0:momdim, * ), qext, qsca, xx, xmu(*)
      real*8 fi1,fi2,fnorm
      complex*16  crefin, sforw, sback, tforw(*), tback(*), s1(*), s2(*)
      character*22  fmt
!
!
      if ( perfct )  write ( *, '(''1'',10x,a,1p,e11.4)' )   &
                      'perfectly conducting case, size parameter =', xx
      if ( .not.perfct )  write ( *, '(''1'',10x,3(a,1p,e11.4))' )   &
                        'refractive index:  real ', dble(crefin),   &
                   '  imag ', dimag(crefin), ',   size parameter =', xx
!
      if ( prnt(1) .and. numang.gt.0 )  then
!
         write ( *, '(/,a)' )   &
          '    cos(angle)  ------- s1 ---------  ------- s2 ---------'//   &
          '  --- s1*conjg(s2) ---   i1=s1**2   i2=s2**2  (i1+i2)/2'//   &
          '  deg polzn'
         do  10  i = 1, numang
            fi1 = dble( s1(i) ) **2 + dimag( s1(i) )**2
            fi2 = dble( s2(i) ) **2 + dimag( s2(i) )**2
            write( *, '( i4, f10.6, 1p,10e11.3 )'   )   &
                    i, xmu(i), s1(i), s2(i), s1(i)*dconjg(s2(i)),   &
                    fi1, fi2, 0.5*(fi1+fi2), (fi2-fi1)/(fi2+fi1)
   10    continue
!
      end if
!
!
      if ( prnt(2) )  then
!
         write ( *, '(/,a,9x,a,17x,a,17x,a,/,(0p,f7.2, 1p,6e12.3) )' )   &
                 '  angle', 's-sub-1', 't-sub-1', 't-sub-2',   &
                     0.0,     sforw,    tforw(1),  tforw(2),   &
                    180.,     sback,    tback(1),  tback(2)
         write ( *, '(/,4(a,1p,e11.4))' )   &
                 ' efficiency factors,  extinction:', qext,   &
                                    '   scattering:', qsca,   &
                                    '   absorption:', qext-qsca,   &
                                 '   rad. pressure:', qext-gqsc
!
         if ( nmom.gt.0 )  then
!
            write( *, '(/,a)' )  ' normalized moments of :'
            if ( ipolzn.eq.0 ) write ( *, '(''+'',27x,a)' ) 'phase fcn'
            if ( ipolzn.gt.0 )  write ( *, '(''+'',33x,a)' )   &
               'm1           m2          s21          d21'
            if ( ipolzn.lt.0 )  write ( *, '(''+'',33x,a)' )   &
               'r1           r2           r3           r4'
!
            fnorm = 4. / ( xx**2 * qsca )
            do  20  m = 0, nmom
               write ( *, '(a,i4)' )  '      moment no.', m
               do 20  j = 1, 4
                  if( calcmo(j) )  then
                     write( fmt, 98 )  24 + (j-1)*13
                     write ( *,fmt )  fnorm * pmom(m,j)
                  end if
   20       continue
         end if
!
      end if
!
      return
!
   98 format( '( ''+'', t', i2, ', 1p,e13.4 )' )
      end subroutine  miprnt  
!**************************************************************************                                            
      subroutine  small1 ( xx, numang, xmu, qext, qsca, gqsc, sforw,   &
                           sback, s1, s2, tforw, tback, a, b )
!
!       small-particle limit of mie quantities in totally reflecting
!       limit ( mie series truncated after 2 terms )
!
!        a,b       first two mie coefficients, with numerator and
!                  denominator expanded in powers of -xx- ( a factor
!                  of xx**3 is missing but is restored before return
!                  to calling program )  ( ref. 2, p. 1508 )
!
      implicit none
      integer  numang,j
      real*8    gqsc, qext, qsca, xx, xmu(*)
      real*8 twothr,fivthr,fivnin,sq,rtmp
      complex*16  a( 2 ), b( 2 ), sforw, sback, s1(*), s2(*),   &
               tforw(*), tback(*)
!
      parameter  ( twothr = 2./3., fivthr = 5./3., fivnin = 5./9. )
      complex*16    ctmp
      sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
!
!
      a( 1 ) = dcmplx ( 0.d0, twothr * ( 1. - 0.2 * xx**2 ) )   &
             / dcmplx ( 1.d0 - 0.5 * xx**2, twothr * xx**3 )
!
      b( 1 ) = dcmplx ( 0.d0, - ( 1. - 0.1 * xx**2 ) / 3. )   &
             / dcmplx ( 1.d0 + 0.5 * xx**2, - xx**3 / 3. )
!
      a( 2 ) = dcmplx ( 0.d0,   xx**2 / 30. )
      b( 2 ) = dcmplx ( 0.d0, - xx**2 / 45. )
!
      qsca = 6. * xx**4 * ( sq( a(1) ) + sq( b(1) )   &
                            + fivthr * ( sq( a(2) ) + sq( b(2) ) ) )
      qext = qsca
      gqsc = 6. * xx**4 * dble( a(1) * dconjg( a(2) + b(1) )   &
                          + ( b(1) + fivnin * a(2) ) * dconjg( b(2) ) )
!
      rtmp = 1.5 * xx**3
      sforw      = rtmp * ( a(1) + b(1) + fivthr * ( a(2) + b(2) ) )
      sback      = rtmp * ( a(1) - b(1) - fivthr * ( a(2) - b(2) ) )
      tforw( 1 ) = rtmp * ( b(1) + fivthr * ( 2.*b(2) - a(2) ) )
      tforw( 2 ) = rtmp * ( a(1) + fivthr * ( 2.*a(2) - b(2) ) )
      tback( 1 ) = rtmp * ( b(1) - fivthr * ( 2.*b(2) + a(2) ) )
      tback( 2 ) = rtmp * ( a(1) - fivthr * ( 2.*a(2) + b(2) ) )
!
      do  10  j = 1, numang
         s1( j ) = rtmp * ( a(1) + b(1) * xmu(j) + fivthr *   &
                    ( a(2) * xmu(j) + b(2) * ( 2.*xmu(j)**2 - 1. )) )
         s2( j ) = rtmp * ( b(1) + a(1) * xmu(j) + fivthr *   &
                    ( b(2) * xmu(j) + a(2) * ( 2.*xmu(j)**2 - 1. )) )
   10 continue
!                                     ** recover actual mie coefficients
      a( 1 ) = xx**3 * a( 1 )
      a( 2 ) = xx**3 * a( 2 )
      b( 1 ) = xx**3 * b( 1 )
      b( 2 ) = xx**3 * b( 2 )
!
      return
      end subroutine  small1 
!*************************************************************************                                                 
      subroutine  small2 ( xx, cior, calcqe, numang, xmu, qext, qsca,   &
                           gqsc, sforw, sback, s1, s2, tforw, tback,   &
                           a, b )
!
!       small-particle limit of mie quantities for general refractive
!       index ( mie series truncated after 2 terms )
!
!        a,b       first two mie coefficients, with numerator and
!                  denominator expanded in powers of -xx- ( a factor
!                  of xx**3 is missing but is restored before return
!                  to calling program )  ( ref. 2, p. 1508 )
!
!        ciorsq    square of refractive index
!
      implicit none
      logical  calcqe
      integer  numang,j
      real*8    gqsc, qext, qsca, xx, xmu(*)
      real*8 twothr,fivthr,sq,rtmp
      complex*16  a( 2 ), b( 2 ), cior, sforw, sback, s1(*), s2(*),   &
               tforw(*), tback(*)
!
      parameter  ( twothr = 2./3., fivthr = 5./3. )
      complex*16  ctmp, ciorsq
      sq( ctmp ) = dble( ctmp )**2 + dimag( ctmp )**2
!
!
      ciorsq = cior**2
      ctmp = dcmplx( 0.d0, twothr ) * ( ciorsq - 1.0 )
      a(1) = ctmp * ( 1.0 - 0.1 * xx**2 + (ciorsq/350. + 1./280.)*xx**4)   &
             / ( ciorsq + 2.0 + ( 1.0 - 0.7 * ciorsq ) * xx**2   &
                 - ( ciorsq**2/175. - 0.275 * ciorsq + 0.25 ) * xx**4   &
                 + xx**3 * ctmp * ( 1.0 - 0.1 * xx**2 ) )
!
      b(1) = (xx**2/30.) * ctmp * ( 1.0 + (ciorsq/35. - 1./14.) *xx**2 )   &
             / ( 1.0 - ( ciorsq/15. - 1./6. ) * xx**2 )
!
      a(2) = ( 0.1 * xx**2 ) * ctmp * ( 1.0 - xx**2 / 14. )   &
             / ( 2. * ciorsq + 3. - ( ciorsq/7. - 0.5 ) * xx**2 )
!
      qsca = 6. * xx**4 * ( sq(a(1)) + sq(b(1)) + fivthr * sq(a(2)) )
      gqsc = 6. * xx**4 * dble( a(1) * dconjg( a(2) + b(1) ) )
      qext = qsca
      if ( calcqe ) qext = 6. * xx * dble( a(1) + b(1) + fivthr * a(2) )
!
      rtmp = 1.5 * xx**3
      sforw      = rtmp * ( a(1) + b(1) + fivthr * a(2) )
      sback      = rtmp * ( a(1) - b(1) - fivthr * a(2) )
      tforw( 1 ) = rtmp * ( b(1) - fivthr * a(2) )
      tforw( 2 ) = rtmp * ( a(1) + 2. * fivthr * a(2) )
      tback( 1 ) = tforw( 1 )
      tback( 2 ) = rtmp * ( a(1) - 2. * fivthr * a(2) )
!
      do  10  j = 1, numang
         s1( j ) = rtmp * ( a(1) + ( b(1) + fivthr * a(2) ) * xmu(j) )
         s2( j ) = rtmp * ( b(1) + a(1) * xmu(j) + fivthr * a(2)   &
                            * ( 2. * xmu(j)**2 - 1. ) )
   10 continue
!                                     ** recover actual mie coefficients
      a( 1 ) = xx**3 * a( 1 )
      a( 2 ) = xx**3 * a( 2 )
      b( 1 ) = xx**3 * b( 1 )
      b( 2 ) = ( 0., 0. )
!
      return
      end subroutine  small2  
!***********************************************************************                                                 
      subroutine  testmi ( qext, qsca, gqsc, sforw, sback, s1, s2,   &
                           tforw, tback, pmom, momdim, ok )
!
!         compare mie code test case results with correct answers
!         and return  ok=false  if even one result is inaccurate.
!
!         the test case is :  mie size parameter = 10
!                             refractive index   = 1.5 - 0.1 i
!                             scattering angle = 140 degrees
!                             1 sekera moment
!
!         results for this case may be found among the test cases
!         at the end of reference (1).
!
!         *** note *** when running on some computers, esp. in single
!         precision, the 'accur' criterion below may have to be relaxed.
!         however, if 'accur' must be set larger than 10**-3 for some
!         size parameters, your computer is probably not accurate
!         enough to do mie computations for those size parameters.
!
      implicit none
      integer momdim,m,n
      real*8    qext, qsca, gqsc, pmom( 0:momdim, * )
      complex*16  sforw, sback, s1(*), s2(*), tforw(*), tback(*)
      logical  ok, wrong
!
      real*8    accur, testqe, testqs, testgq, testpm( 0:1 )
      complex*16 testsf, testsb,tests1,tests2,testtf(2), testtb(2)
      data   testqe / 2.459791 /,  testqs / 1.235144 /,   &
             testgq / 1.139235 /,  testsf / ( 61.49476, -3.177994 ) /,   &
             testsb / ( 1.493434, 0.2963657 ) /,   &
             tests1 / ( -0.1548380, -1.128972) /,   &
             tests2 / ( 0.05669755, 0.5425681) /,   &
             testtf / ( 12.95238, -136.6436 ), ( 48.54238, 133.4656 ) /,   &
             testtb / ( 41.88414, -15.57833 ), ( 43.37758, -15.28196 )/,   &
             testpm / 227.1975, 183.6898 /
      real*8 calc,exact
!      data   accur / 1.e-5 /
      data   accur / 1.e-4 /
      wrong( calc, exact ) = dabs( (calc - exact) / exact ) .gt. accur
!
!
      ok = .true.
      if ( wrong( qext,testqe ) )   &
           call  tstbad( 'qext', abs((qext - testqe) / testqe), ok )
      if ( wrong( qsca,testqs ) )   &
           call  tstbad( 'qsca', abs((qsca - testqs) / testqs), ok )
      if ( wrong( gqsc,testgq ) )   &
           call  tstbad( 'gqsc', abs((gqsc - testgq) / testgq), ok )
!
      if ( wrong(  dble(sforw),  dble(testsf) ) .or.   &
           wrong( dimag(sforw), dimag(testsf) ) )   &
           call  tstbad( 'sforw', cdabs((sforw - testsf) / testsf), ok )
!
      if ( wrong(  dble(sback),  dble(testsb) ) .or.   &
           wrong( dimag(sback), dimag(testsb) ) )   &
           call  tstbad( 'sback', cdabs((sback - testsb) / testsb), ok )
!
      if ( wrong(  dble(s1(1)),  dble(tests1) ) .or.   &
           wrong( dimag(s1(1)), dimag(tests1) ) )   &
           call  tstbad( 's1', cdabs((s1(1) - tests1) / tests1), ok )
!
      if ( wrong(  dble(s2(1)),  dble(tests2) ) .or.   &
           wrong( dimag(s2(1)), dimag(tests2) ) )   &
           call  tstbad( 's2', cdabs((s2(1) - tests2) / tests2), ok )
!
      do  20  n = 1, 2
         if ( wrong(  dble(tforw(n)),  dble(testtf(n)) ) .or.   &
              wrong( dimag(tforw(n)), dimag(testtf(n)) ) )   &
              call  tstbad( 'tforw', cdabs( (tforw(n) - testtf(n)) /   &
                                           testtf(n) ), ok )
         if ( wrong(  dble(tback(n)),  dble(testtb(n)) ) .or.   &
              wrong( dimag(tback(n)), dimag(testtb(n)) ) )   &
              call  tstbad( 'tback', cdabs( (tback(n) - testtb(n)) /   &
                                            testtb(n) ), ok )
   20 continue
!
      do  30  m = 0, 1
         if ( wrong( pmom(m,1), testpm(m) ) )   &
              call  tstbad( 'pmom', dabs( (pmom(m,1)-testpm(m)) /   &
                                         testpm(m) ), ok )
   30 continue
!
      return
!
      end subroutine  testmi  
!**************************************************************************                                              
      subroutine  errmsg( messag, fatal )
!
!        print out a warning or error message;  abort if error
!
      USE module_peg_util, only:  peg_message, peg_error_fatal

      implicit none
      logical       fatal
      logical, save :: once
      data once / .false. /
      character*80 msg
      character*(*) messag
      integer, save :: maxmsg, nummsg
      data nummsg / 0 /,  maxmsg / 100 /
!
!
      if ( fatal )  then
!         write ( *, '(2a)' )  ' ******* error >>>>>>  ', messag
!         stop
   	  write( msg, '(a)' )   &
                  'optical averaging mie fatal error ' //   &
                  messag                  
                  call peg_message( lunerr, msg )             
                  call peg_error_fatal( lunerr, msg )
      end if
!
      nummsg = nummsg + 1
      if ( nummsg.gt.maxmsg )  then
!         if ( .not.once )  write ( *,99 )
	 if ( .not.once )then
	    write( msg, '(a)' )   &
             'optical averaging mie: too many warning messages -- no longer printing '
            call peg_message( lunerr, msg )
         end if    
         once = .true.
      else
         msg =   'optical averaging mie warning '  // messag
         call peg_message( lunerr, msg )  
!         write ( *, '(2a)' )  ' ******* warning >>>>>>  ', messag
      endif
!
      return
!
!   99 format( ///,' >>>>>>  too many warning messages --  ',   &
!         'they will no longer be printed  <<<<<<<', /// )
      end subroutine  errmsg  
!********************************************************************                     
      subroutine  wrtbad ( varnam, erflag )
!
!          write names of erroneous variables
!
!      input :   varnam = name of erroneous variable to be written
!                         ( character, any length )
!
!      output :  erflag = logical flag, set true by this routine
! ----------------------------------------------------------------------
      USE module_peg_util, only:  peg_message
      
      implicit none
      character*(*)  varnam
      logical        erflag
      character*80 msg
      integer, save :: maxmsg, nummsg
      data  nummsg / 0 /,  maxmsg / 50 /     
!
!
      nummsg = nummsg + 1
!      write ( *, '(3a)' )  ' ****  input variable  ', varnam,   &
!                           '  in error  ****'
        msg = 'optical averaging mie input variable in error ' // varnam                   
      call peg_message( lunerr, msg )
      erflag = .true.
      if ( nummsg.eq.maxmsg )   &     
         call  errmsg ( 'too many input variable errors.  aborting...$', .true. )
      return
!
      end subroutine  wrtbad 
!******************************************************************                        
      subroutine  tstbad( varnam, relerr, ok )
!
!       write name (-varnam-) of variable failing self-test and its
!       percent error from the correct value.  return  ok = false.
!
      implicit none
      character*(*)  varnam
      logical        ok
      real*8          relerr
!
!
      ok = .false.
      write( *, '(/,3a,1p,e11.2,a)' )   &
             ' output variable  ', varnam,'  differed by', 100.*relerr,   &
             '  per cent from correct value.  self-test failed.'
      return
!
      end subroutine  tstbad                      
!******************************************************************
!
      subroutine sect02(dgnum_um,sigmag,drydens,iflag,duma,nbin,dlo_um,dhi_um, &
        xnum_sect,xmas_sect)
!
!   user specifies a single log-normal mode and a set of section boundaries
!   prog calculates mass and number for each section
!
        implicit none
        REAL, DIMENSION(nbin), INTENT(OUT) :: xnum_sect, xmas_sect
        integer iflag, n, nbin 
        real &
          dgnum, dgnum_um, dhi, dhi_um, dlo, dlo_um,   &
          drydens, dstar, duma, dumfrac, dx,   &
          sigmag, sumnum, summas,   &
          sx, sxroot2, thi, tlo, vtot,   &
          x0, x3, xhi, xlo, xmtot, xntot, xvtot
        real dlo_sect(nbin), dhi_sect(nbin)
!       real erfc_num_recipes
        real pi
        parameter (pi = 3.1415926536)
!
        if (iflag .le. 1) then
            xntot = duma
        else
            xmtot = duma
        end if
!   compute total volume and number for mode
        dgnum = dgnum_um*1.0e-4
        sx = alog( sigmag )
        x0 = alog( dgnum )
        x3 = x0 + 3.*sx*sx
        dstar = dgnum * exp(1.5*sx*sx)
        if (iflag .le. 1) then
            xvtot = xntot*(pi/6.0)*dstar*dstar*dstar
            xmtot = xvtot*drydens*1.0e12
        else
            xvtot = xmtot/(drydens*1.0e12)
            xntot = xvtot/((pi/6.0)*dstar*dstar*dstar)
        end if
!   compute section boundaries
        dlo = dlo_um*1.0e-4
        dhi = dhi_um*1.0e-4
        xlo = log( dlo )
        xhi = log( dhi )
        dx = (xhi - xlo)/nbin
        do n = 1, nbin
            dlo_sect(n) = exp( xlo + dx*(n-1) )
            dhi_sect(n) = exp( xlo + dx*n )
        end do
!   compute modal "working" parameters including total num/vol/mass
        dgnum = dgnum_um*1.0e-4
        sx = alog( sigmag )
        x0 = alog( dgnum )
        x3 = x0 + 3.*sx*sx
        dstar = dgnum * exp(1.5*sx*sx)
        if (iflag .le. 1) then
            xvtot = xntot*(pi/6.0)*dstar*dstar*dstar
            xmtot = xvtot*drydens*1.0e12
        else
            xvtot = xmtot/(drydens*1.0e12)
            xntot = xvtot/((pi/6.0)*dstar*dstar*dstar)
        end if
!   compute number and mass for each section
        sxroot2 = sx * sqrt( 2.0 )
        sumnum = 0.
        summas = 0.
!       write(22,*)
!       write(22,*) 'dgnum_um, sigmag = ', dgnum_um, sigmag
!       write(22,*) 'drydens =', drydens
!       write(22,*) 'ntot (#/cm3), mtot (ug/m3) = ', xntot, xmtot
!        write(22,9220)
!9220    format( /   &
!        '  n   dlo(um)   dhi(um)       number         mass' / )
!9225    format(   i3, 2f10.6, 2(1pe13.4) )
!9230    format( / 'sum over all sections  ', 2(1pe13.4) )
!9231    format(   'modal totals           ', 2(1pe13.4) )
        do n = 1, nbin
            xlo = alog( dlo_sect(n) )
            xhi = alog( dhi_sect(n) )
            tlo = (xlo - x0)/sxroot2
            thi = (xhi - x0)/sxroot2
            if (tlo .le. 0.) then
                dumfrac = 0.5*( erfc_num_recipes(-thi) - erfc_num_recipes(-tlo) )
            else
                dumfrac = 0.5*( erfc_num_recipes(tlo) - erfc_num_recipes(thi) )
            end if
            xnum_sect(n) = xntot*dumfrac
            tlo = (xlo - x3)/sxroot2
            thi = (xhi - x3)/sxroot2
            if (tlo .le. 0.) then
                dumfrac = 0.5*( erfc_num_recipes(-thi) - erfc_num_recipes(-tlo) )
            else
                dumfrac = 0.5*( erfc_num_recipes(tlo) - erfc_num_recipes(thi) )
            end if
            xmas_sect(n) = xmtot*dumfrac
            sumnum = sumnum + xnum_sect(n)
            summas = summas + xmas_sect(n)
!           write(22,9225) n, 1.e4*dlo_sect(n), 1.e4*dhi_sect(n),   &
!               xnum_sect(n), xmas_sect(n)
        end do
!       write(22,9230) sumnum, summas
!       write(22,9231) xntot, xmtot

      end subroutine  sect02
!-----------------------------------------------------------------------
        real function erfc_num_recipes( x )
!
!   from press et al, numerical recipes, 1990, page 164
!
        implicit none
        real x
        double precision erfc_dbl, dum, t, z
        z = abs(x)
        t = 1.0/(1.0 + 0.5*z)
!       erfc_num_recipes =
!     &   t*exp( -z*z - 1.26551223 + t*(1.00002368 + t*(0.37409196 +
!     &   t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 +
!     &                                    t*(-1.13520398 +
!     &   t*(1.48851587 + t*(-0.82215223 + t*0.17087277 )))))))))
        dum =  ( -z*z - 1.26551223 + t*(1.00002368 + t*(0.37409196 +   &
          t*(0.09678418 + t*(-0.18628806 + t*(0.27886807 +   &
                                           t*(-1.13520398 +   &
          t*(1.48851587 + t*(-0.82215223 + t*0.17087277 )))))))))
        erfc_dbl = t * exp(dum)
        if (x .lt. 0.0) erfc_dbl = 2.0d0 - erfc_dbl
        erfc_num_recipes = erfc_dbl
        return

        end function erfc_num_recipes
!-----------------------------------------------------------------------

!****************************************************************************
! <1.> subr mieaer_sc
! Purpose:  calculate aerosol optical depth, single scattering albedo,
!   asymmetry factor, extinction, Legendre coefficients, and average aerosol
!   size. parameterizes aerosol coefficients using a full-blown Mie code
! Calculates these properties for either (1) an aerosol internally mixed in a
! shell/core configuration or (2) internally mixed aerosol represented by
! volume averaging of refractive indices
!  Uses the ACKMIE code developed eons ago by Tom Ackerman  (Ackerman and Toon, 1981:
! absorption of visible radiation in the atmosphere containing mixtures of absorbing and 
! non-absorbing particles, Appl. Opt., 20, 3661-3668.
!
! INPUT
!       id -- grid id number
!   	iclm, jclm -- i,j of grid column being processed
!       nbin_a -- number of bins
!	number_bin_col(nbin_a,kmaxd) --   number density in layer, #/cm^3
!	radius_wet_col(nbin_a,kmaxd) -- wet radius, shell, cm
!       radius_core_col(nbin_a,kmaxd) -- core radius, cm; NOTE:
!		if this is set to zero, the code will assumed a volume averaging
!		of refractive indices
!    	refindx_col(nbin_a,kmaxd) -- volume complex index of refraction for shell, or
!       	volume averaged complex index of refraction for the whole aerosol
!		in volume averaged mode
!	refindx_core_col(nbin_a,kmaxd) -- complex index of refraction for core
!	dz -- depth of individual cells in column, m
!	isecfrm0 - time from start of run, sec
!	lpar -- number of grid cells in vertical (via module_fastj_cmnh)
!   kmaxd -- predefined maximum allowed levels from module_data_mosaic_other
!            passed here via module_fastj_cmnh
! OUTPUT: saved in module_fastj_cmnmie
!   real tauaer  ! aerosol optical depth
!        waer    ! aerosol single scattering albedo
!        gaer    ! aerosol asymmetery factor
!        extaer  ! aerosol extinction
!        l2,l3,l4,l5,l6,l7 ! Legendre coefficients, numbered 0,1,2,......
!        sizeaer ! average wet radius
!	 bscoef ! aerosol backscatter coefficient with units km-1 * steradian -1  JCB 2007/02/01
!***********************************************************************

        subroutine mieaer_sc( &
	          id, iclm, jclm, nbin_a,   &
              number_bin_col, radius_wet_col, refindx_col,   &
              radius_core_col, refindx_core_col, &  ! jcb, 2007/07/25; for shell/core implementation, set radius_cor_col=0 for volume-average configuration
              dz, isecfrm0, lpar, &
              sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7,bscoef)  ! added bscoef JCB 2007/02/01
!
	USE module_data_mosaic_other, only : kmaxd
	USE module_data_mosaic_therm, only : nbin_a_maxd
	USE module_peg_util, only : peg_message

        IMPLICIT NONE
!   subr arguments
!jdf
        integer,parameter :: nspint = 4 ! Num of spectral intervals across
                                        ! solar spectrum for FAST-J
        integer, intent(in) :: lpar
!jdf    real, dimension (nspint, kmaxd+1),intent(out) :: sizeaer,extaer,waer,gaer,tauaer
!jdf    real, dimension (nspint, kmaxd+1),intent(out) :: l2,l3,l4,l5,l6,l7
!jdf    real, dimension (nspint, kmaxd+1),intent(out) :: bscoef  !JCB 2007/02/01
        real, dimension (nspint, lpar+1),intent(out) :: sizeaer,extaer,waer,gaer,tauaer
        real, dimension (nspint, lpar+1),intent(out) :: l2,l3,l4,l5,l6,l7
        real, dimension (nspint, lpar+1),intent(out) :: bscoef  !JCB 2007/02/01
        real, dimension (nspint),save :: wavmid !cm
        data wavmid     &
            / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
!jdf
    integer, intent(in) :: id, iclm, jclm, nbin_a, isecfrm0
!jdf    real, intent(in), dimension(nbin_a, kmaxd) :: number_bin_col
!jdf    real, intent(inout), dimension(nbin_a, kmaxd) :: radius_wet_col
!jdf    complex, intent(in) :: refindx_col(nbin_a, kmaxd)
    real, intent(in), dimension(nbin_a, lpar+1) :: number_bin_col
    real, intent(inout), dimension(nbin_a, lpar+1) :: radius_wet_col, radius_core_col  ! jcb  2007/07/25
    complex, intent(in) :: refindx_col(nbin_a, lpar+1), refindx_core_col(nbin_a,lpar+1)  ! jcb  2007/07/25, 
    real, intent(in)    :: dz(lpar)
	real thesum, sum ! for normalizing things and testing
!
      integer m,l,j,nl,ll,nc,klevel
      integer       ns, &       ! Spectral loop index
                    i,  &      ! Longitude loop index
                    k         ! Level loop index

      real*8 dp_wet_a,dp_core_a
      complex*16 ri_shell_a,ri_core_a
      real*8 qextc,qscatc,qbackc,extc,scatc,backc,gscac
      real*8 vlambc
      integer n,kkk,jjj
      real*8 pmom(0:7,1)
      real weighte, weights, pscat
      real pie,sizem
      real ratio
!
	real,save ::rmin,rmax  ! min, max aerosol size bin
!       data rmin /0.005e-04/   ! rmin in cm, 5e-3 microns min allowable size
!       data rmax /50.0e-04/    ! rmax in cm. 50 microns, big particle, max allowable size
	data rmin /0.010e-04/   ! rmin in cm, 5e-3 microns min allowable size
	data rmax /7.0e-04/    ! rmax in cm. 50 microns, big particle, max allowable size
! diagnostic declarations
      integer kcallmieaer2
      integer ibin
      character*150 msg

#if (defined(CHEM_DBG_I) && defined(CHEM_DBG_J) && defined(CHEM_DBG_K))
!cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ec  diagnostics
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ec run_out.25 has aerosol physical parameter info for bins 1-8
!ec and vertical cells 1 to kmaxd.
!        ilaporte = 33
!        jlaporte = 34
	if (iclm .eq. CHEM_DBG_I) then
	  if (jclm .eq. CHEM_DBG_J) then
!   initial entry
	   if (kcallmieaer2 .eq. 0) then
	      write(*,9099)iclm, jclm
 9099	format('for cell i = ', i3, 2x, 'j = ', i3)	
              write(*,9100)
 9100     format(   &
               'isecfrm0', 3x, 'i', 3x, 'j', 3x,'k', 3x,   &
      	       'ibin', 3x,   &
               'refindx_col(ibin,k)', 3x,   &
               'radius_wet_col(ibin,k)', 3x,   &
               'number_bin_col(ibin,k)'   &
               )
	   end if
!ec output for run_out.25
            do k = 1, lpar
            do ibin = 1, nbin_a
              write(*, 9120)   &
                 isecfrm0,iclm, jclm, k, ibin,   &
                 refindx_col(ibin,k),   &
                 radius_wet_col(ibin,k),   &
                 number_bin_col(ibin,k)
9120    format( i7,3(2x,i4),2x,i4, 4x, 4(e14.6,2x))
            end do
            end do
	kcallmieaer2 = kcallmieaer2 + 1
        end if
        end if
!ec end print of aerosol physical parameter diagnostics	
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
#endif
!
! loop over levels
	do 2000 klevel=1,lpar
	thesum=0.0
	do m=1,nbin_a
	thesum=thesum+number_bin_col(m,klevel)
	enddo
      pie=4.*atan(1.)
! Begin spectral loop
      do 1000 ns=1,nspint
!        aerosol optical properties
               tauaer(ns,klevel)=0.
               waer(ns,klevel)=0.
               gaer(ns,klevel)=0.
	       sizeaer(ns,klevel)=0.0
		extaer(ns,klevel)=0.0
		l2(ns,klevel)=0.0
		l3(ns,klevel)=0.0
		l4(ns,klevel)=0.0
		l5(ns,klevel)=0.0
		l6(ns,klevel)=0.0
		l7(ns,klevel)=0.0
		bscoef(ns,klevel)=0.0
		if(thesum.le.1.e-21)goto 1000  ! set everything = 0 if no aerosol ! wig changed 0.0 to 1e-21

! loop over the bins, nbin_a is the number of bins
               do m=1,nbin_a
! check to see if there's any aerosol
!jdf   		if(number_bin_col(m,klevel).le.1e-21)goto 70  ! no aerosol wig changed 0.0 to 1e-21, 31-Oct-2005
! here's the size
		sizem=radius_wet_col(m,klevel) ! radius in cm
                ratio=radius_core_col(m,klevel)/radius_wet_col(m,klevel)
! check limits of particle size
! rce 2004-dec-07 - use klevel in write statements
		if(radius_wet_col(m,klevel).le.rmin)then
		  radius_wet_col(m,klevel)=rmin
                  radius_core_col(m,klevel)=rmin*ratio
		  write( msg, '(a, 5i4,1x, e11.4)' )	&
		  'mieaer_sc: radius_wet set to rmin,'  //	&
		  'id,i,j,k,m,rm(m,k)', id, iclm, jclm, klevel, m, radius_wet_col(m,klevel)
		  call peg_message( lunerr, msg )
!		   write(6,'('' particle size too small '')')
		endif
!
		if(radius_wet_col(m,klevel).gt.rmax)then
		   write( msg, '(a, 5i4,1x, e11.4)' )	&
                'mieaer_sc: radius_wet set to rmax,'  //	&
                'id,i,j,k,m,rm(m,k)', &
                id, iclm, jclm, klevel, m, radius_wet_col(m,klevel)
           call peg_message( lunerr, msg )
           radius_wet_col(m,klevel)=rmax
           radius_core_col(m,klevel)=rmax*ratio
!           write(6,'('' particle size too large '')')
		endif
!
                  ri_shell_a=dcmplx(real(refindx_col(m,klevel)),abs(aimag(refindx_col(m,klevel)))) ! need positive complex part of refractive index here
                  ri_core_a=dcmplx(real(refindx_core_col(m,klevel)),abs(aimag(refindx_core_col(m,klevel))))  ! need positive complex part of refractive index here
!
		dp_wet_a= 2.0*radius_wet_col(m,klevel)*1.0e04  ! radius_wet is in cm,dp_wet_a should be in microns
		dp_core_a=2.0*radius_core_col(m,klevel)*1.0e04
		vlambc=wavmid(ns)*1.0e04
!
!	write(6,*)dp_wet_a
!	write(6,*)ri_shell_a
!	write(6,*)vlambc
		call miedriver(dp_wet_a,dp_core_a,ri_shell_a,ri_core_a, vlambc, &
         	qextc,qscatc,gscac,extc,scatc,qbackc,backc,pmom)
! check, note that pmom(1,1)/pmom(0,1) is indeed the asymmetry parameter as calculated by Tom's code, jcb, July 7, 2007
! correct in the Rayleigh limit, July 3, 2007: jcb
!     	write(6,*)
!     	do ii=0,7
!	write(6,*)pmom(ii,1),pmom(ii,1)/pmom(0,1)
!	enddo
!   	write(6,*)qextc,qscatc,gscac,extc,scatc
!	write(6,*)
!
	weighte=extc*1.0e-08 ! extinction cross section, converted to cm^2
	weights=scatc*1.0e-08 ! scattering cross section, converted to cm^2
       	tauaer(ns,klevel)=tauaer(ns,klevel)+weighte* &
        number_bin_col(m,klevel)  ! must be multiplied by deltaZ
	sizeaer(ns,klevel)=sizeaer(ns,klevel)+radius_wet_col(m,klevel)*10000.0*  &
        number_bin_col(m,klevel)
        waer(ns,klevel)=waer(ns,klevel)+weights*number_bin_col(m,klevel)
        gaer(ns,klevel)=gaer(ns,klevel)+gscac*weights*  &
        number_bin_col(m,klevel)
	l2(ns,klevel)=l2(ns,klevel)+weights*pmom(2,1)/pmom(0,1)*5.0*number_bin_col(m,klevel)
	l3(ns,klevel)=l3(ns,klevel)+weights*pmom(3,1)/pmom(0,1)*7.0*number_bin_col(m,klevel)
	l4(ns,klevel)=l4(ns,klevel)+weights*pmom(4,1)/pmom(0,1)*9.0*number_bin_col(m,klevel)
	l5(ns,klevel)=l5(ns,klevel)+weights*pmom(5,1)/pmom(0,1)*11.0*number_bin_col(m,klevel)
	l6(ns,klevel)=l6(ns,klevel)+weights*pmom(6,1)/pmom(0,1)*13.0*number_bin_col(m,klevel)
	l7(ns,klevel)=l7(ns,klevel)+weights*pmom(7,1)/pmom(0,1)*15.0*number_bin_col(m,klevel)
! the 4*pi gives the correct value in the Rayleigh limit compared with the old core, which we assume is correct
	bscoef(ns,klevel)=bscoef(ns,klevel)+backc*1.0e-08*number_bin_col(m,klevel)*4.0*pie ! converting cross-section from microns ^2 to cm^2, 4*pie needed
2001	continue
        end do ! end of nbin loop
! take averages
	sizeaer(ns,klevel)=sizeaer(ns,klevel)/thesum
	gaer(ns,klevel)=gaer(ns,klevel)/waer(ns,klevel)
	l2(ns,klevel)=l2(ns,klevel)/waer(ns,klevel)
	l3(ns,klevel)=l3(ns,klevel)/waer(ns,klevel)
	l4(ns,klevel)=l4(ns,klevel)/waer(ns,klevel)
	l5(ns,klevel)=l5(ns,klevel)/waer(ns,klevel)
	l6(ns,klevel)=l6(ns,klevel)/waer(ns,klevel)
	l7(ns,klevel)=l7(ns,klevel)/waer(ns,klevel)
!	write(6,*)ns,klevel,l4(ns,klevel)
! this is beta, not beta/(4*pie)
	bscoef(ns,klevel)=bscoef(ns,klevel)*1.0e5 ! unit (km)^-1
! SSA checked by comparson with Travis and Hansen, get exact result
	waer(ns,klevel)=waer(ns,klevel)/tauaer(ns,klevel)  ! must be last
	extaer(ns,klevel)=tauaer(ns,klevel)*1.0e5 ! unit (km)^-1
 70   continue ! end of nbin_a loop
 1000 continue  ! end of wavelength loop
2000   continue  ! end of klevel loop
! before returning, multiply tauaer by depth of individual cells.
! tauaer is in cm-1, dz in m; multiply dz by 100 to convert from m to cm.
  	do ns = 1, nspint
  	do klevel = 1, lpar
  	   tauaer(ns,klevel) = tauaer(ns,klevel) * dz(klevel)* 100.   
        end do
        end do 	
!
#if (defined(CHEM_DBG_I) && defined(CHEM_DBG_J) && defined(CHEM_DBG_K))
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ec  fastj diagnostics
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
!ec run_out.30 has aerosol optical info for cells 1 to kmaxd.
!        ilaporte = 33
!         jlaporte = 34
 	if (iclm .eq. CHEM_DBG_I) then
 	  if (jclm .eq. CHEM_DBG_J) then
!   initial entry
 	   if (kcallmieaer .eq. 0) then
 	       write(*,909) CHEM_DBG_I, CHEM_DBG_J
 909	format(	' for cell i = ', i3, ' j = ', i3)          	
               write(*,910)
 910     format(   &
               'isecfrm0', 3x, 'i', 3x, 'j', 3x,'k', 3x,   &
      	        'dzmfastj', 8x,   &
               'tauaer(1,k)',1x, 'tauaer(2,k)',1x,'tauaer(3,k)',3x,   &
               'tauaer(4,k)',5x,   &
               'waer(1,k)', 7x, 'waer(2,k)', 7x,'waer(3,k)', 7x,   &
               'waer(4,k)', 7x,   &
               'gaer(1,k)', 7x, 'gaer(2,k)', 7x,'gaer(3,k)', 7x,   &
               'gaer(4,k)', 7x,   &
               'extaer(1,k)',5x, 'extaer(2,k)',5x,'extaer(3,k)',5x,   &
               'extaer(4,k)',5x,   &
               'sizeaer(1,k)',4x, 'sizeaer(2,k)',4x,'sizeaer(3,k)',4x,   &
               'sizeaer(4,k)'  )
 	   end if
!ec output for run_out.30
         do k = 1, lpar
         write(*, 912)   &
           isecfrm0,iclm, jclm, k,   &
           dz(k) ,   &
           (tauaer(n,k),   n=1,4),   &
           (waer(n,k),     n=1,4),   &
           (gaer(n,k),     n=1,4),   &
           (extaer(n,k),   n=1,4),   &
           (sizeaer(n,k),  n=1,4)
 912    format( i7,3(2x,i4),2x,21(e14.6,2x))
         end do
 	kcallmieaer = kcallmieaer + 1
        end if
         end if
!ec end print of fastj diagnostics	
!ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
#endif
!
      return
      end subroutine mieaer_sc
!
	subroutine miedriver(dp_wet_a,dp_core_a,ri_shell_a,ri_core_a, vlambc, &
        qextc,qscatc,gscac,extc,scatc,qbackc,backc,pmom)
! MOSAIC INPUTS
!	dp_wet_a = diameter (cm) of aerosol
!	dp_core_a = diameter (cm) of the aerosol's core
!	ri_shell_a = refractive index (complex) of shell
!	ri_core_a = refractve index (complex ) of core (usually assumed to be LAC)
!	vlambc = wavelength of calculation (um, convert to cm)
! MOSAIC outputs
! 	qextc = scattering efficiency
!	qscac = scattering efficiency
!	gscac = asymmetry parameter
!	extc = extinction cross section  (cm^2)
!	scac = scattering cross section (cm^2)
! drives concentric sphere program
! /*---------------------------------------------------------------*/
! /* INPUTS:                                                       */
! /*---------------------------------------------------------------*/

!   VLAMBc: Wavelength of the radiation
!   NRGFLAGc: Flag to indicate a number density of volume radius
!   RGc: Number (RGN = Rm) or volume (RGV) weighted mean radius of
!        the particle size distribution
!   SIGMAGc: Geometric standard deviation of the distribution
!   SHELRc: Real part of the index of refraction for the shell
!   SHELIc: Imaginary part of the index of refraction for the shell
!   RINc: Inner core radius as a fraction of outer shell radius
!   CORERc: Real part of the index of refraction for the core
!   COREIc: Imaginary part of the index of refraction for the core
!   NANG: Number of scattering angles between 0 and 90 degrees,
!         inclusive

! /*---------------------------------------------------------------*/
! /* OUTPUTS:                                                      */
! /*---------------------------------------------------------------*/

!   QEXTc: Extinction efficiency of the particle
!   QSCAc: Scattering efficiency of the particle
!   QBACKc: Backscatter efficiency of the particle
!   EXTc: Extinction cross section of the particle
!   SCAc: Scattering cross section of the particle
!   BACKc: Backscatter cross section of the particle
!   GSCA: Asymmetry parameter of the particles phase function
!   ANGLES(NAN): Scattering angles in degrees
!   S1R(NAN): Real part of the amplitude scattering matrix
!   S1C(NAN): Complex part of the amplitude scattering matrix
!   S2R(NAN): Real part of the amplitude scattering matrix
!   S2C(NAN): Complex part of the amplitude scattering matrix
!   S11N: Normalization coefficient of the scattering matrix
!   S11(NAN): S11 scattering coefficients
!   S12(NAN): S12 scattering coefficients
!   S33(NAN): S33 scattering coefficients
!   S34(NAN): S34 scattering coefficients
!   SPOL(NAN): Degree of polarization of unpolarized, incident light
!   SP(NAN): Phase function
!
! NOTE: NAN=2*NANG-1 is the number of scattering angles between
!       0 and 180 degrees, inclusive.
! /*---------------------------------------------------------------*/
      REAL*8 VLAMBc,RGcmin,RGcmax,RGc,SIGMAGc,SHELRc,SHELIc
      REAL*8 RINc,CORERc,COREIc
      INTEGER*4 NRGFLAGc,NANG
      REAL*8 QEXTc,QSCATc,QBACKc,EXTc,SCATc,BACKc,GSCAc
      REAL*8 ANGLESc(200),S1R(200),S1C(200),S2R(200),S2C(200)
      REAL*8 S11N,S11(200),S12(200),S33(200),S34(200),SPOL(200),SP(200)
      real*8 pmom(0:7,1)
      real*8 dp_wet_a,dp_core_a
      complex*16 ri_shell_a,ri_core_a
!
	nang=2 ! only one angle
	nrgflagc=0 ! size distribution
!
	rgc=dp_wet_a/2.0 ! radius of particle
	rinc=dp_core_a/dp_wet_a ! fraction of radius that is the core
	rgcmin=0.001
	rgcmax=5.0
	sigmagc=1.0  ! no particle size dispersion
	shelrc=real(ri_shell_a)
	shelic=aimag(ri_shell_a)
	corerc=real(ri_core_a)
	coreic=aimag(ri_core_a)
       CALL ACKMIEPARTICLE( VLAMBc,NRGFLAGc,RGcmin,RGcmax, &
                  RGc,SIGMAGc,SHELRc, &
                  SHELIc, RINc,CORERc,COREIc,NANG,QEXTc,QSCATc, &
                  QBACKc, EXTc,SCATc,BACKc, GSCAc, &
               ANGLESc,S1R,S1C,S2R,S2C,S11N,S11,S12,S33,S34,SPOL,SP,pmom)  ! jcb
!	write(6,1010)rgc,qextc,qscatc,qscatc/qextc,gscac
1010	format(5f20.12)
1020	format(2f12.6)
	end subroutine miedriver
!
!     /*--------------------------------------------------------*/
!     /* The Toon-Ackerman SUBROUTINE DMIESS for calculating the*/
!     /* scatter off of a coated sphere of some sort.           */
!     /* Toon and Ackerman, Applied Optics, Vol. 20, Pg. 3657   */
!     /*--------------------------------------------------------*/

!**********************************************

      SUBROUTINE DMIESS(  RO,      RFR,     RFI,     THETD,     JX,  &
                         QEXT,    QSCAT,   CTBRQS,  ELTRMX,    PIE, &
                         TAU,     CSTHT,   SI2THT,  ACAP, QBS, IT,  &
                         LL,      R,       RE2,     TMAG2,     WVNO, an,bn, ntrm  )
!
! **********************************************************************
!    THIS SUBROUTINE COMPUTES MIE SCATTERING BY A STRATIFIED SPHERE,
!    I.E. A PARTICLE CONSISTING OF A SPHERICAL CORE SURROUNDED BY A
!    SPHERICAL SHELL.  THE BASIC CODE USED WAS THAT DESCRIBED IN THE
!    REPORT: SUBROUTINES FOR COMPUTING THE PARAMETERS OF THE
!    ELECTROMAGNETIC RADIATION SCATTERED BY A SPHERE J.V. DAVE,
!    I B M SCIENTIFIC CENTER, PALO ALTO , CALIFORNIA.
!    REPORT NO. 320 - 3236 .. MAY 1968 .
!
!    THE MODIFICATIONS FOR STRATIFIED SPHERES ARE DESCRIBED IN
!        TOON AND ACKERMAN, APPL. OPTICS, IN PRESS, 1981
!
!    THE PARAMETERS IN THE CALLING STATEMENT ARE DEFINED AS FOLLOWS :
!      RO IS THE OUTER (SHELL) RADIUS;
!      R  IS THE CORE RADIUS;
!      RFR, RFI  ARE THE REAL AND IMAGINARY PARTS OF THE SHELL INDEX
!          OF REFRACTION IN THE FORM (RFR - I* RFI);
!      RE2, TMAG2  ARE THE INDEX PARTS FOR THE CORE;
!          ( WE ASSUME SPACE HAS UNIT INDEX. )
!      THETD(J): ANGLE IN DEGREES BETWEEN THE DIRECTIONS OF THE INCIDENT
!          AND THE SCATTERED RADIATION.  THETD(J) IS< OR= 90.0
!          IF THETD(J) SHOULD HAPPEN TO BE GREATER THAN 90.0, ENTER WITH
!          SUPPLEMENTARY VALUE, SEE COMMENTS BELOW ON ELTRMX;
!      JX: TOTAL NUMBER OF THETD FOR WHICH THE COMPUTATIONS ARE
!          REQUIRED.  JX SHOULD NOT EXCEED IT UNLESS THE DIMENSIONS
!          STATEMENTS ARE APPROPRIATEDLY MODIFIED;
!
!      THE DEFINITIONS FOR THE FOLLOWING SYMBOLS CAN BE FOUND IN LIGHT
!          SCATTERING BY SMALL PARTICLES, H.C.VAN DE HULST, JOHN WILEY
!          SONS, INC., NEW YORK, 1957.
!      QEXT: EFFIECIENCY FACTOR FOR EXTINCTION,VAN DE HULST,P.14 127.
!      QSCAT: EFFIECINCY FACTOR FOR SCATTERING,V.D. HULST,P.14 127.
!      CTBRQS: AVERAGE(COSINE THETA) * QSCAT,VAN DE HULST,P.128
!      ELTRMX(I,J,K): ELEMENTS OF THE TRANSFORMATION MATRIX F,V.D.HULST
!          ,P.34,45 125. I=1: ELEMENT M SUB 2..I=2: ELEMENT M SUB 1..
!          I = 3: ELEMENT S SUB 21.. I = 4: ELEMENT D SUB 21..
!      ELTRMX(I,J,1) REPRESENTS THE ITH ELEMENT OF THE MATRIX FOR
!          THE ANGLE THETD(J).. ELTRMX(I,J,2) REPRESENTS THE ITH ELEMENT
!          OF THE MATRIX FOR THE ANGLE 180.0 - THETD(J) ..
!      QBS IS THE BACK SCATTER CROSS SECTION.
!
!      IT: IS THE DIMENSION OF THETD, ELTRMX, CSTHT, PIE, TAU, SI2THT,
!          IT MUST CORRESPOND EXACTLY TO THE SECOND DIMENSION OF ELTRMX.
!      LL: IS THE DIMENSION OF ACAP
!          IN THE ORIGINAL PROGRAM THE DIMENSION OF ACAP WAS 7000.
!          FOR CONSERVING SPACE THIS SHOULD BE NOT MUCH HIGHER THAN
!          THE VALUE, N=1.1*(NREAL**2 + NIMAG**2)**.5 * X + 1
!      WVNO: 2*PIE / WAVELENGTH
!
!    ALSO THE SUBROUTINE COMPUTES THE CAPITAL A FUNCTION BY MAKING USE O
!    DOWNWARD RECURRENCE RELATIONSHIP.
!
!      TA(1): REAL PART OF WFN(1).  TA(2): IMAGINARY PART OF WFN(1).
!      TA(3): REAL PART OF WFN(2).  TA(4): IMAGINARY PART OF WFN(2).
!      TB(1): REAL PART OF FNA.     TB(2): IMAGINARY PART OF FNA.
!      TC(1): REAL PART OF FNB.     TC(2): IMAGINARY PART OF FNB.
!      TD(1): REAL PART OF FNAP.    TD(2): IMAGINARY PART OF FNAP.
!      TE(1): REAL PART OF FNBP.    TE(2): IMAGINARY PART OF FNBP.
!      FNAP, FNBP  ARE THE PRECEDING VALUES OF FNA, FNB RESPECTIVELY.
! **********************************************************************

!     /*--------------------------------------------------------------*/
!     /* Initially, make all types undefined.                         */
!     /*--------------------------------------------------------------*/

!      IMPLICIT UNDEFINED(A-Z)

!     /*--------------------------------------------------------*/
!     /* Dimension statements.                                  */
!     /*--------------------------------------------------------*/

      INTEGER*4  JX, IT, LL

      REAL*8     RO, RFR, RFI, THETD(IT), QEXT, QSCAT, CTBRQS, &
                ELTRMX(4,IT,2), PIE(3,IT), TAU(3,IT), CSTHT(IT), &
                SI2THT(IT), QBS, R,  RE2, TMAG2, WVNO

      COMPLEX*16 ACAP(LL)

!     /*--------------------------------------------------------*/
!     /* Variables used in the calculations below.              */
!     /*--------------------------------------------------------*/

      INTEGER*4  IFLAG, J, K, M, N, NN, NMX1, NMX2

      REAL*8     T(5), TA(4), TB(2), TC(2), TD(2), TE(2), X, &
                RX, X1, Y1, X4, Y4, SINX1, SINX4, COSX1, COSX4, &
                EY1, E2Y1, EY4, EY1MY4, EY1PY4, AA, BB, &
                CC, DD, DENOM, REALP, AMAGP, QBSR, QBSI, RMM, &
                PIG, RXP4
!
      COMPLEX*16 FNAP,      FNBP,      W,  &
                FNA,       FNB,       RF,           RRF, &
                RRFX,      WM1,       FN1,          FN2,   &
                TC1,       TC2,       WFN(2),       Z(4), &
                K1,        K2,        K3,    &
                RC,        U(8),      DH1, &
                DH2,       DH4,       P24H24,       P24H21, &
                PSTORE,    HSTORE,    DUMMY,        DUMSQ
! jcb
      complex*16 an(500),bn(500) ! a,b Mie coefficients, jcb  Hansen and Travis, eqn 2.44
      integer*4 ntrm
!
!     /*--------------------------------------------------------*/
!     /* Define the common block.                               */
!     /*--------------------------------------------------------*/

      COMMON / WARRAY / W(3,9000)

!
!      EQUIVALENCE   (FNA,TB(1)),(FNB,TC(1)),(FNAP,TD(1)),(FNBP,TE(1))
!
!   IF THE CORE IS SMALL SCATTERING IS COMPUTED FOR THE SHELL ONLY
!

!     /*--------------------------------------------------------*/
!     /* Begin the Mie calculations.                            */
!     /*--------------------------------------------------------*/
      IFLAG = 1
      ntrm=0 ! jcb
      IF ( R/RO .LT. 1.0D-06 )   IFLAG = 2                              
      IF ( JX .LE. IT )   GO TO 20
         WRITE( 6,7 )                                                   
         WRITE( 6,6 )
         STOP 30                                                        
   20 RF =  CMPLX( RFR,  -RFI )
      RC =  CMPLX( RE2,-TMAG2 )
      X  =  RO * WVNO
      K1 =  RC * WVNO
      K2 =  RF * WVNO
      K3 =  CMPLX( WVNO, 0.0D0 )                                          
      Z(1) =  K2 * RO
      Z(2) =  K3 * RO                                                   
      Z(3) =  K1 * R
      Z(4) =  K2 * R                                                    
      X1   =  DREAL( Z(1) )
      Y1   =  DIMAG( Z(1) )                                              
      X4   =  DREAL( Z(4) )
      Y4   =  DIMAG( Z(4) )                                              
      RRF  =  1.0D0 / RF
      RX   =  1.0D0 / X                                                   
      RRFX =  RRF * RX
      T(1) =  ( X**2 ) * ( RFR**2 + RFI**2 )                            
      T(1) =  DSQRT( T(1) )
      NMX1 =  1.30D0* T(1)
!
      IF ( NMX1 .LE. LL-1 )   GO TO 21
         WRITE(6,8)
         STOP 32                                                        
   21 NMX2 = T(1) * 1.2
	nmx1=min(nmx1+5,150)  ! jcb
	nmx2=min(nmx2+5,135)  ! jcb
!   	write(6,*)x,nmx1,nmx2,ll  ! jcb
!	stop
      IF ( NMX1 .GT.  150 )   GO TO 22
!        NMX1 = 150
!        NMX2 = 135
!
   22 ACAP( NMX1+1 )  =  ( 0.0D0,0.0D0 )
      IF ( IFLAG .EQ. 2 )   GO TO 26
         DO 29   N = 1,3                                                
   29    W( N,NMX1+1 )  =  ( 0.0D0,0.0D0 )
   26 CONTINUE                                                          
      DO 23   N = 1,NMX1
         NN = NMX1 - N + 1
         ACAP(NN) = (NN+1)*RRFX - 1.0D0 / ((NN+1)*RRFX + ACAP(NN+1))
         IF ( IFLAG .EQ. 2 )   GO TO 23                                 
            DO 31   M = 1,3
   31       W( M,NN ) = (NN+1) / Z(M+1)  -   &
                        1.0D0 / ((NN+1) / Z(M+1) + W( M,NN+1 ))
   23 CONTINUE                                                          
!
      DO 30   J = 1,JX                                                  
      IF ( THETD(J) .LT. 0.0D0 )  THETD(J) =  DABS( THETD(J) )
      IF ( THETD(J) .GT. 0.0D0 )  GO TO 24                                
      CSTHT(J)  = 1.0D0
      SI2THT(J) = 0.0D0                                                   
      GO TO 30
   24 IF ( THETD(J) .GE. 90.0D0 )  GO TO 25                               
      T(1)      =  ( 3.14159265359 * THETD(J) ) / 180.0D0
      CSTHT(J)  =  DCOS( T(1) )                                          
      SI2THT(J) =  1.0D0 - CSTHT(J)**2
      GO TO 30                                                          
   25 IF ( THETD(J) .GT. 90.0 )  GO TO 28
      CSTHT(J)  =  0.0D0                                               
      SI2THT(J) =  1.0D0
      GO TO 30                                                          
   28 WRITE( 6,5 )  THETD(J)
      WRITE( 6,6 )                                                      
      STOP 34
   30 CONTINUE                                                          
!
      DO 35  J = 1,JX                                                   
      PIE(1,J) =  0.0D0
      PIE(2,J) =  1.0D0
      TAU(1,J) =  0.0D0
      TAU(2,J) =  CSTHT(J)                                              
   35 CONTINUE
!
! INITIALIZATION OF HOMOGENEOUS SPHERE
!
      T(1)   =  DCOS(X)
      T(2)   =  DSIN(X)                                                  
      WM1    =  CMPLX( T(1),-T(2) )
      WFN(1) =  CMPLX( T(2), T(1) )                                     
      TA(1)  =  T(2)
      TA(2)  =  T(1)                                                    
      WFN(2) =  RX * WFN(1) - WM1
      TA(3)  =  DREAL(WFN(2))                                           
      TA(4)  =  DIMAG(WFN(2))
!
	n=1 ! jcb, bug???
      IF ( IFLAG .EQ. 2 )   GO TO 560
      N = 1                                                             
!
! INITIALIZATION PROCEDURE FOR STRATIFIED SPHERE BEGINS HERE
!
      SINX1   =  DSIN( X1 )                                              
      SINX4   =  DSIN( X4 )
      COSX1   =  DCOS( X1 )                                              
      COSX4   =  DCOS( X4 )
      EY1     =  DEXP( Y1 )                                              
      E2Y1    =  EY1 * EY1
      EY4     =  DEXP( Y4 )                                              
      EY1MY4  =  DEXP( Y1 - Y4 )
      EY1PY4  =  EY1 * EY4                                              
      EY1MY4  =  DEXP( Y1 - Y4 )
      AA  =  SINX4 * ( EY1PY4 + EY1MY4 )                                
      BB  =  COSX4 * ( EY1PY4 - EY1MY4 )
      CC  =  SINX1 * ( E2Y1 + 1.0D0 )
      DD  =  COSX1 * ( E2Y1 - 1.0D0 )
      DENOM   =  1.0D0  +  E2Y1 * (4.0D0*SINX1*SINX1 - 2.0D0 + E2Y1)    
      REALP   =  ( AA * CC  +  BB * DD ) / DENOM
      AMAGP   =  ( BB * CC  -  AA * DD ) / DENOM
      DUMMY   =  CMPLX( REALP, AMAGP )
      AA  =  SINX4 * SINX4 - 0.5D0                                        
      BB  =  COSX4 * SINX4
      P24H24  =  0.5D0 + CMPLX( AA,BB ) * EY4 * EY4                       
      AA  =  SINX1 * SINX4  -  COSX1 * COSX4
      BB  =  SINX1 * COSX4  +  COSX1 * SINX4                            
      CC  =  SINX1 * SINX4  +  COSX1 * COSX4
      DD  = -SINX1 * COSX4  +  COSX1 * SINX4                            
      P24H21  =  0.5D0 * CMPLX( AA,BB ) * EY1 * EY4  + &
                0.5D0 * CMPLX( CC,DD ) * EY1MY4
      DH4  =  Z(4) / (1.0D0 + (0.0D0,1.0D0) * Z(4))  -  1.0D0 / Z(4)
      DH1  =  Z(1) / (1.0D0 + (0.0D0,1.0D0) * Z(1))  -  1.0D0 / Z(1)        
      DH2  =  Z(2) / (1.0D0 + (0.0D0,1.0D0) * Z(2))  -  1.0D0 / Z(2)
      PSTORE  =  ( DH4 + N / Z(4) )  *  ( W(3,N) + N / Z(4) )           
      P24H24  =  P24H24 / PSTORE
      HSTORE  =  ( DH1 + N / Z(1) )  *  ( W(3,N) + N / Z(4) )           
      P24H21  =  P24H21 / HSTORE
      PSTORE  =  ( ACAP(N) + N / Z(1) )  /  ( W(3,N) + N / Z(4) )       
      DUMMY   =  DUMMY * PSTORE
      DUMSQ   =  DUMMY * DUMMY                                          
!
! NOTE:  THE DEFINITIONS OF U(I) IN THIS PROGRAM ARE NOT THE SAME AS
!        THE USUBI DEFINED IN THE ARTICLE BY TOON AND ACKERMAN.  THE
!        CORRESPONDING TERMS ARE:
!          USUB1 = U(1)                       USUB2 = U(5)
!          USUB3 = U(7)                       USUB4 = DUMSQ
!          USUB5 = U(2)                       USUB6 = U(3)
!          USUB7 = U(6)                       USUB8 = U(4)
!          RATIO OF SPHERICAL BESSEL FTN TO SPHERICAL HENKAL FTN = U(8)
!
      U(1) =  K3 * ACAP(N)  -  K2 * W(1,N)
      U(2) =  K3 * ACAP(N)  -  K2 * DH2                                 
      U(3) =  K2 * ACAP(N)  -  K3 * W(1,N)
      U(4) =  K2 * ACAP(N)  -  K3 * DH2                                 
      U(5) =  K1 *  W(3,N)  -  K2 * W(2,N)
      U(6) =  K2 *  W(3,N)  -  K1 * W(2,N)                              
      U(7) =  ( 0.0D0,-1.0D0 )  *  ( DUMMY * P24H21 - P24H24 )
      U(8) =  TA(3) / WFN(2)                                            
!
      FNA  =  U(8) * ( U(1)*U(5)*U(7)  +  K1*U(1)  -  DUMSQ*K3*U(5) ) / &
                    ( U(2)*U(5)*U(7)  +  K1*U(2)  -  DUMSQ*K3*U(5) )
      FNB  =  U(8) * ( U(3)*U(6)*U(7)  +  K2*U(3)  -  DUMSQ*K2*U(6) ) / &
                    ( U(4)*U(6)*U(7)  +  K2*U(4)  -  DUMSQ*K2*U(6) )
!
!  Explicit equivalences added by J. Francis

      TB(1) = DREAL(FNA)
      TB(2) = DIMAG(FNA)
      TC(1) = DREAL(FNB)
      TC(2) = DIMAG(FNB)
      GO TO 561                                                         
  560 TC1  =  ACAP(1) * RRF  +  RX
      TC2  =  ACAP(1) * RF   +  RX                                      
      FNA  =  ( TC1 * TA(3)  -  TA(1) ) / ( TC1 * WFN(2)  -  WFN(1) )
      FNB  =  ( TC2 * TA(3)  -  TA(1) ) / ( TC2 * WFN(2)  -  WFN(1) )   
      TB(1) = DREAL(FNA)
      TB(2) = DIMAG(FNA)
      TC(1) = DREAL(FNB)
      TC(2) = DIMAG(FNB)
!
  561 CONTINUE
! jcb
	ntrm=ntrm+1
	an(n)=fna
	bn(n)=fnb
!	write(6,1010)ntrm,n,an(n),bn(n)
1010	format(2i5,4e15.6)
! jcb
      FNAP = FNA
      FNBP = FNB                                                        
      TD(1) = DREAL(FNAP)
      TD(2) = DIMAG(FNAP)
      TE(1) = DREAL(FNBP)
      TE(2) = DIMAG(FNBP)
      T(1) = 1.50D0
!
!    FROM HERE TO THE STATMENT NUMBER 90, ELTRMX(I,J,K) HAS
!    FOLLOWING MEANING:
!    ELTRMX(1,J,K): REAL PART OF THE FIRST COMPLEX AMPLITUDE.
!    ELTRMX(2,J,K): IMAGINARY PART OF THE FIRST COMPLEX AMPLITUDE.
!    ELTRMX(3,J,K): REAL PART OF THE SECOND COMPLEX AMPLITUDE.
!    ELTRMX(4,J,K): IMAGINARY PART OF THE SECOND COMPLEX AMPLITUDE.
!    K = 1 : FOR THETD(J) AND K = 2 : FOR 180.0 - THETD(J)
!    DEFINITION OF THE COMPLEX AMPLITUDE: VAN DE HULST,P.125.
!
      TB(1) = T(1) * TB(1)                                              
      TB(2) = T(1) * TB(2)
      TC(1) = T(1) * TC(1)                                              
      TC(2) = T(1) * TC(2)
!      DO 60 J = 1,JX
!          ELTRMX(1,J,1) = TB(1) * PIE(2,J) + TC(1) * TAU(2,J)
!          ELTRMX(2,J,1) = TB(2) * PIE(2,J) + TC(2) * TAU(2,J)
!          ELTRMX(3,J,1) = TC(1) * PIE(2,J) + TB(1) * TAU(2,J)
!          ELTRMX(4,J,1) = TC(2) * PIE(2,J) + TB(2) * TAU(2,J)
!          ELTRMX(1,J,2) = TB(1) * PIE(2,J) - TC(1) * TAU(2,J)
!          ELTRMX(2,J,2) = TB(2) * PIE(2,J) - TC(2) * TAU(2,J)
!          ELTRMX(3,J,2) = TC(1) * PIE(2,J) - TB(1) * TAU(2,J)
!          ELTRMX(4,J,2) = TC(2) * PIE(2,J) - TB(2) * TAU(2,J)
   60 CONTINUE
!
      QEXT   = 2.0D0 * ( TB(1) + TC(1))
      QSCAT  = ( TB(1)**2 + TB(2)**2 + TC(1)**2 + TC(2)**2 ) / 0.75D0     
      CTBRQS = 0.0D0
      QBSR   = -2.0D0*(TC(1) - TB(1))                                     
      QBSI   = -2.0D0*(TC(2) - TB(2))
      RMM    = -1.0D0                                                     
      N = 2
   65 T(1) = 2*N - 1   ! start of loop, JCB
      T(2) =   N - 1
      T(3) = 2*N + 1
      DO 70  J = 1,JX
          PIE(3,J) = ( T(1)*PIE(2,J)*CSTHT(J) - N*PIE(1,J) ) / T(2) 
          TAU(3,J) = CSTHT(J) * ( PIE(3,J) - PIE(1,J) ) - &
                               T(1)*SI2THT(J)*PIE(2,J) + TAU(1,J)
   70 CONTINUE
!
! HERE SET UP HOMOGENEOUS SPHERE
!
      WM1    =  WFN(1)
      WFN(1) =  WFN(2)                                                  
      TA(1)  =  DREAL(WFN(1))
      TA(2)  =  DIMAG(WFN(1))                                            
      WFN(2) =  T(1) * RX * WFN(1)  -  WM1
      TA(3)  =  DREAL(WFN(2))                                            
      TA(4)  =  DIMAG(WFN(2))
!
      IF ( IFLAG .EQ. 2 )   GO TO 1000
!
! HERE SET UP STRATIFIED SPHERE
!
      DH2  =  - N / Z(2)  +  1.0D0 / ( N / Z(2) - DH2 )
      DH4  =  - N / Z(4)  +  1.0D0 / ( N / Z(4) - DH4 )                   
      DH1  =  - N / Z(1)  +  1.0D0 / ( N / Z(1) - DH1 )
      PSTORE  =  ( DH4 + N / Z(4) )  *  ( W(3,N) + N / Z(4) )           
      P24H24  =  P24H24 / PSTORE
      HSTORE  =  ( DH1 + N / Z(1) )  *  ( W(3,N) + N / Z(4) )           
      P24H21  =  P24H21 / HSTORE
      PSTORE  =  ( ACAP(N) + N / Z(1) )  /  ( W(3,N) + N / Z(4) )       
      DUMMY   =  DUMMY * PSTORE
      DUMSQ   =  DUMMY * DUMMY                                          
!
      U(1) =  K3 * ACAP(N)  -  K2 * W(1,N)                              
      U(2) =  K3 * ACAP(N)  -  K2 * DH2
      U(3) =  K2 * ACAP(N)  -  K3 * W(1,N)
      U(4) =  K2 * ACAP(N)  -  K3 * DH2
      U(5) =  K1 *  W(3,N)  -  K2 * W(2,N)                              
      U(6) =  K2 *  W(3,N)  -  K1 * W(2,N)
      U(7) =  ( 0.0D0,-1.0D0 )  *  ( DUMMY * P24H21 - P24H24 )              
      U(8) =  TA(3) / WFN(2)
!
      FNA  =  U(8) * ( U(1)*U(5)*U(7)  +  K1*U(1)  -  DUMSQ*K3*U(5) ) / &
                    ( U(2)*U(5)*U(7)  +  K1*U(2)  -  DUMSQ*K3*U(5) )
      FNB  =  U(8) * ( U(3)*U(6)*U(7)  +  K2*U(3)  -  DUMSQ*K2*U(6) ) / &
                    ( U(4)*U(6)*U(7)  +  K2*U(4)  -  DUMSQ*K2*U(6) )
      TB(1) = DREAL(FNA)
      TB(2) = DIMAG(FNA)
      TC(1) = DREAL(FNB)
      TC(2) = DIMAG(FNB)
!
 1000 CONTINUE                                                          
      TC1  =  ACAP(N) * RRF  +  N * RX
      TC2  =  ACAP(N) * RF   +  N * RX                                  
      FN1  =  ( TC1 * TA(3)  -  TA(1) ) /  ( TC1 * WFN(2) - WFN(1) )
      FN2  =  ( TC2 * TA(3)  -  TA(1) ) /  ( TC2 * WFN(2) - WFN(1) )    
      M    =  WVNO * R
      IF ( N .LT. M )   GO TO 1002                                      
      IF ( IFLAG .EQ. 2 )   GO TO 1001
      IF ( ABS(  ( FN1-FNA ) / FN1  )  .LT. 1.0D-09   .AND.    &
          ABS(  ( FN2-FNB ) / FN2  )  .LT. 1.0D-09  )     IFLAG = 2
      IF ( IFLAG .EQ. 1 )   GO TO 1002                                  
 1001 FNA  =  FN1
      FNB  =  FN2                                                       
      TB(1) = DREAL(FNA)
      TB(2) = DIMAG(FNA)
      TC(1) = DREAL(FNB)
      TC(2) = DIMAG(FNB)
!
 1002 CONTINUE
! jcb
 	ntrm=ntrm+1
 	an(n)=fna
	bn(n)=fnb
!	write(6,1010)ntrm,n,an(n),bn(n)
! jcb
      T(5)  =  N
      T(4)  =  T(1) / ( T(5) * T(2) )                                   
      T(2)  =  (  T(2) * ( T(5) + 1.0D0 )  ) / T(5)
!
      CTBRQS  =  CTBRQS  +  T(2) * ( TD(1) * TB(1)  +  TD(2) * TB(2)  &
                        +           TE(1) * TC(1)  +  TE(2) * TC(2) )  &
                        +  T(4) * ( TD(1) * TE(1)  +  TD(2) * TE(2) )
      QEXT    =   QEXT  +  T(3) * ( TB(1) + TC(1) )                     
      T(4)    =  TB(1)**2 + TB(2)**2 + TC(1)**2 + TC(2)**2
      QSCAT   =  QSCAT  +  T(3) * T(4)                                  
      RMM     =  -RMM
      QBSR    =  QBSR + T(3)*RMM*(TC(1) - TB(1))                        
      QBSI    =  QBSI + T(3)*RMM*(TC(2) - TB(2))
!
      T(2)    =  N * (N+1)
      T(1)    =  T(3) / T(2)                                            
      K = (N/2)*2
!      DO 80 J = 1,JX
!       ELTRMX(1,J,1)=ELTRMX(1,J,1)+T(1)*(TB(1)*PIE(3,J)+TC(1)*TAU(3,J))
!       ELTRMX(2,J,1)=ELTRMX(2,J,1)+T(1)*(TB(2)*PIE(3,J)+TC(2)*TAU(3,J))
!       ELTRMX(3,J,1)=ELTRMX(3,J,1)+T(1)*(TC(1)*PIE(3,J)+TB(1)*TAU(3,J))
!       ELTRMX(4,J,1)=ELTRMX(4,J,1)+T(1)*(TC(2)*PIE(3,J)+TB(2)*TAU(3,J))
!      IF ( K .EQ. N )  THEN
!       ELTRMX(1,J,2)=ELTRMX(1,J,2)+T(1)*(-TB(1)*PIE(3,J)+TC(1)*TAU(3,J))
!       ELTRMX(2,J,2)=ELTRMX(2,J,2)+T(1)*(-TB(2)*PIE(3,J)+TC(2)*TAU(3,J))
!       ELTRMX(3,J,2)=ELTRMX(3,J,2)+T(1)*(-TC(1)*PIE(3,J)+TB(1)*TAU(3,J))
!       ELTRMX(4,J,2)=ELTRMX(4,J,2)+T(1)*(-TC(2)*PIE(3,J)+TB(2)*TAU(3,J))
!      ELSE
!       ELTRMX(1,J,2)=ELTRMX(1,J,2)+T(1)*(TB(1)*PIE(3,J)-TC(1)*TAU(3,J))
!       ELTRMX(2,J,2)=ELTRMX(2,J,2)+T(1)*(TB(2)*PIE(3,J)-TC(2)*TAU(3,J))
!       ELTRMX(3,J,2)=ELTRMX(3,J,2)+T(1)*(TC(1)*PIE(3,J)-TB(1)*TAU(3,J))
!       ELTRMX(4,J,2)=ELTRMX(4,J,2)+T(1)*(TC(2)*PIE(3,J)-TB(2)*TAU(3,J))
!      END IF
!   80 CONTINUE
!
!      IF ( T(4) .LT. 1.0D-14 )   GO TO 100         ! bail out of loop
      IF ( T(4) .LT. 1.0D-10 .OR. N .GE. NMX2)   GO TO 100         ! bail out of loop, JCB
      N = N + 1
!      DO 90 J = 1,JX
!         PIE(1,J)  =   PIE(2,J)
!         PIE(2,J)  =   PIE(3,J)
!         TAU(1,J)  =  TAU(2,J)
!         TAU(2,J)  =  TAU(3,J)
   90 CONTINUE                                                          
      FNAP  =  FNA
      FNBP  =  FNB                                                      
      TD(1) = DREAL(FNAP)
      TD(2) = DIMAG(FNAP)
      TE(1) = DREAL(FNBP)
      TE(2) = DIMAG(FNBP)
      IF ( N .LE. NMX2 )   GO TO 65
         WRITE( 6,8 )                                                   
         STOP 36
  100 CONTINUE
!      DO 120 J = 1,JX
!      DO 120 K = 1,2
!         DO  115  I= 1,4
!         T(I)  =  ELTRMX(I,J,K)
!  115    CONTINUE
!         ELTRMX(2,J,K)  =      T(1)**2  +  T(2)**2
!         ELTRMX(1,J,K)  =      T(3)**2  +  T(4)**2
!         ELTRMX(3,J,K)  =  T(1) * T(3)  +  T(2) * T(4)
!         ELTRMX(4,J,K)  =  T(2) * T(3)  -  T(4) * T(1)
!  120 CONTINUE
      T(1)    =  2.0D0 * RX**2
      QEXT    =   QEXT * T(1)                                           
      QSCAT   =  QSCAT * T(1)
      CTBRQS  =  2.0D0 * CTBRQS * T(1)                                    
!
! QBS IS THE BACK SCATTER CROSS SECTION
!
      PIG   = DACOS(-1.0D0)                                                
      RXP4  = RX*RX/(4.0D0*PIG)
      QBS   = RXP4*(QBSR**2 + QBSI**2)                                  
!
    5 FORMAT( 10X,' THE VALUE OF THE SCATTERING ANGLE IS GREATER THAN 90.0 DEGREES. IT IS ', E15.4 )
    6 FORMAT( // 10X, 'PLEASE READ COMMENTS.' // )
    7 FORMAT( // 10X, 'THE VALUE OF THE ARGUMENT JX IS GREATER THAN IT')
    8 FORMAT( // 10X, 'THE UPPER LIMIT FOR ACAP IS NOT ENOUGH. SUGGEST GET DETAILED OUTPUT AND MODIFY SUBROUTINE' // )
!
      RETURN
      END SUBROUTINE DMIESS
!
! /*****************************************************************/
! /* SUBROUTINE ACKMIEPARICLE                                      */
! /*****************************************************************/

! THIS PROGRAM COMPUTES SCATTERING PROPERTIES FOR DISTRIBUTIONS OF
! PARTICLES COMPOSED OF A CORE OF ONE MATERIAL AND A SHELL OF ANOTHER.

! /*---------------------------------------------------------------*/
! /* INPUTS:                                                       */
! /*---------------------------------------------------------------*/

!   VLAMBc: Wavelength of the radiation
!   NRGFLAGc: Flag to indicate a number density of volume radius
!   RGc: Geometric mean radius of the particle distribution
!   SIGMAGc: Geometric standard deviation of the distribution
!   SHELRc: Real part of the index of refraction for the shell
!   SHELIc: Imaginary part of the index of refraction for the shell
!   RINc: Inner core radius as a fraction of outer shell radius
!   CORERc: Real part of the index of refraction for the core
!   COREIc: Imaginary part of the index of refraction for the core
!   NANG: Number of scattering angles between 0 and 90 degrees,
!         inclusive

! /*---------------------------------------------------------------*/
! /* OUTPUTS:                                                      */
! /*---------------------------------------------------------------*/

!   QEXTc: Extinction efficiency of the particle
!   QSCAc: Scattering efficiency of the particle
!   QBACKc: Backscatter efficiency of the particle
!   EXTc: Extinction cross section of the particle
!   SCAc: Scattering cross section of the particle
!   BACKc: Backscatter cross section of the particle
!   GSCA: Asymmetry parameter of the particle's phase function
!   ANGLES(NAN): Scattering angles in degrees
!   S1R(NAN): Real part of the amplitude scattering matrix
!   S1C(NAN): Complex part of the amplitude scattering matrix
!   S2R(NAN): Real part of the amplitude scattering matrix
!   S2C(NAN): Complex part of the amplitude scattering matrix
!   S11N: Normalization coefficient of the scattering matrix
!   S11(NAN): S11 scattering coefficients
!   S12(NAN): S12 scattering coefficients
!   S33(NAN): S33 scattering coefficients
!   S34(NAN): S34 scattering coefficients
!   SPOL(NAN): Degree of polarization of unpolarized, incident light
!   SP(NAN): Phase function
!
! NOTE: NAN=2*NANG-1 is the number of scattering angles between
!       0 and 180 degrees, inclusive.
! /*---------------------------------------------------------------*/

      SUBROUTINE ACKMIEPARTICLE( VLAMBc,NRGFLAGc,RGcmin,RGcmax, &
                  RGc,SIGMAGc,SHELRc, &
                  SHELIc, RINc,CORERc,COREIc,NANG,QEXTc,QSCATc, &
                  QBACKc, EXTc,SCATc,BACKc, GSCAc, &
               ANGLESc,S1R,S1C,S2R,S2C,S11N,S11,S12,S33,S34,SPOL,SP,pmom)  ! jcb

!     /*--------------------------------------------------------*/
!     /* Set reals to 8 bytes, i.e., double precision.          */
!     /*--------------------------------------------------------*/

      IMPLICIT REAL*8 (A-H, O-Z)

!     /*--------------------------------------------------------*/
!     /* Parameter statements.                                  */
!     /*--------------------------------------------------------*/

      integer*4 mxnang

      PARAMETER(MXNANG=501)

!     /*--------------------------------------------------------*/
!     /* Dimension statements.                                  */
!     /*--------------------------------------------------------*/

      REAL*8 VLAMBc,RGcmin,RGcmax,RGc,SIGMAGc,SHELRc,SHELIc
      REAL*8 RINc,CORERc,COREIc
      INTEGER*4 NANG,NRGFLAGc,NSCATH
      REAL*8 QEXTc,QSCATc,QBACKc,EXTc,SCATc,BACKc,GSCAc
      REAL*8 ANGLESc(*),S1R(*),S1C(*),S2R(*),S2C(*)
      REAL*8 S11N,S11(*),S12(*),S33(*),S34(*),SPOL(*),SP(*)

!     /*--------------------------------------------------------*/
!     /* Define the types of the common block.                  */
!     /*--------------------------------------------------------*/

      INTEGER*4 IPHASEmie

      REAL*8 ALAMB, RGmin, RGmax, RGV, SIGMAG, RGCFRAC, RFRS,RFIS, RFRC, RFIC

! for calculating the Legendre coefficient, jcb
        complex*16 an(500),bn(500) ! a,b Mie coefficients, jcb  Hansen and Travis, eqn 2.44
	integer*4 ntrmj,ntrm,nmom,ipolzn,momdim
	real*8 pmom(0:7,1)

!     /*--------------------------------------------------------*/
!     /* Set reals to 8 bytes, i.e., double precision.          */
!     /*--------------------------------------------------------*/

!      IMPLICIT REAL*8 (A-H, O-Z)

!     /*--------------------------------------------------------*/
!     /* Input common block for scattering calculations.        */
!     /*--------------------------------------------------------*/

!jdf  COMMON / PHASE  / IPHASEmie

!jdf  COMMON / INPUTS / ALAMB, RGmin, RGmax, RGV, SIGMAG, &
!jdf                    RGCFRAC, RFRS,  RFIS,  RFRC,  RFIC

!     /*--------------------------------------------------------*/
!     /* Output common block for scattering calculations.       */
!     /*--------------------------------------------------------*/

!jdf  COMMON / OUTPUTS / QEXT, QSCAT, QBS, EXT, SCAT, BSCAT, ASY

!     /*--------------------------------------------------------*/
!     /* Arrays to hold the results of the scattering and       */
!     /* moment calculations.                                   */
!     /*--------------------------------------------------------*/

      REAL*8 COSPHI(2*MXNANG-1), SCTPHS(2*MXNANG-1)

!     /*--------------------------------------------------------*/
!     /* Copy the input parameters into the common block INPUTS */
!     /*--------------------------------------------------------*/

         NSCATH  = NANG

         IPHASEmie = 0
         ALAMB   = VLAMBc
         RGmin   = RGcmin
         RGmax   = RGcmax
         RGV     = RGc
         SIGMAG  = SIGMAGc
         RGCFRAC = RINc
         RFRS    = SHELRc
         RFIS    = SHELIc
         RFRC    = CORERc
         RFIC    = COREIc
 
!     /*--------------------------------------------------------*/
!     /* Calculate the particle scattering properties for the   */
!     /* given wavelength, particle distribution and indices of */
!     /* refraction of inner and outer material.                */
!     /*--------------------------------------------------------*/

         CALL PFCNPARTICLE(NSCATH, COSPHI, SCTPHS, &
            ANGLESc,S1R,S1C,S2R,S2C,S11N,S11,S12,S33,S34,SPOL,SP, an, bn, ntrm,  & ! jcb
            ALAMB,RGmin,RGmax,RGV,SIGMAG,RGCFRAC,RFRS,RFIS,RFRC,RFIC,           & ! jdf
            QEXT,QSCAT,QBS,EXT,SCAT,BSCAT,ASY,                                  & ! jdf
            IPHASEmie)                                                            ! jdf

!     /*--------------------------------------------------------*/
!     /* If IPHASE = 1, then the full phase function was        */
!     /* calculated; now, go calculate its moments.             */
!     /*--------------------------------------------------------*/

!        IF (IPHASE .gt. 0) CALL DISMOM (NSCATA,COSPHI,SCTPHS,RMOMS)

!     /*--------------------------------------------------------*/
!     /* Copy the variables in the common block OUTPUTS to the  */
!     /* variable addresses passed into this routine.           */
!     /*--------------------------------------------------------*/

         QEXTc  = QEXT
         QSCATc = QSCAT
         QBACKc = QBS

         EXTc  = EXT
         SCATc = SCAT
         BACKc = BSCAT

         GSCAc  = ASY
! jcb
!	ntrmj = number of terms in Mie series, jcb
	nmom= 7 ! largest Legendre coefficient to calculate 0:7  (8 total), jcb
	ipolzn=0 ! unpolarized light, jcb
	momdim=7 ! dimension of pmom, pmom(0:7), jcb
! a, b = Mie coefficients
! pmom = output of Legendre coefficients, pmom(0:7)
!	write(6,*)ntrm
!	do ii=1,ntrm
!	write(6,1030)ii,an(ii),bn(ii)
1030	format(i5,4e15.6)
!	enddo

	call lpcoefjcb(ntrm,nmom,ipolzn,momdim,an,bn,pmom)
!	do ii=0,7
!	write(6,1040)ii,pmom(ii,1),pmom(ii,1)/pmom(0,1)
!1040	format(i5,2e15.6)
!	enddo
!     /*--------------------------------------------------------*/
!     /* FORMAT statements.                                     */
!     /*--------------------------------------------------------*/

  107 FORMAT ( ///, 1X, I6, ' IS AN INVALID MEAN RADIUS FLAG')

!     /*--------------------------------------------------------*/
!     /* DONE with this subroutine so exit.                     */
!     /*--------------------------------------------------------*/

  END SUBROUTINE ACKMIEPARTICLE

! /*****************************************************************/
! /* SUBROUTINE PFCNPARTICLE                                       */
! /*****************************************************************/
!
! THIS SUBROUTINE COMPUTES THE PHASE FUNCTION SCTPHS(I) AT NSCATA
! ANGLES BETWEEN 0.0 AND 180.0 DEGREES SPECIFIED BY COSPHI(I) WHICH
! CONTAINS THE ANGLE COSINES. THESE VALUES ARE RETURNED TO FUNCTION
! PHASFN FOR AZIMUTHAL AVERAGING.
! INPUT DATA FOR THIS ROUTINE IS PASSED THROUGH COMMON /SIZDIS/
! AND CONSISTS OF THE FOLLOWING PARAMETERS
! NEWSD  = 1 IF SIZE DIS VALUES HAVE NOT PREVIOUSLY BEEN USED IN
!          THIS ROUTINE,  = 0 OTHERWISE.
! RGV    = GEOMETRIC MEAN RADIUS FOR THE VOLUME DISTRIBUTION OF THE
!          SPECIFIED PARTICLES
! SIGMAG = GEOMETRIC STANDARD DEVIATION
! RFR,I  = REAL AND IMAGINARY INDEX OF REFRACTION OF PARTICLES
! ALAMB  = WAVELENGTH AT WHICH CALCULATIONS ARE TO BE PERFORMED
!
! /*---------------------------------------------------------------*/

      SUBROUTINE PFCNPARTICLE( NSCATH, COSPHI, SCTPHS, &
          ANGLESc,S1R,S1C,S2R,S2C,S11N,S11,S12,S33,S34,SPOL,SP,an,bn,ntrm, & ! jcb
          ALAMB,RGmin,RGmax,RGV,SIGMAG,RGCFRAC,RFRS,RFIS,RFRC,RFIC,        & ! jdf
          QEXT,QSCAT,QBS,EXT,SCAT,BSCAT,ASY,                               & ! jdf
          IPHASEmie)                                                         ! jdf

!     /*--------------------------------------------------------*/
!     /* Set reals to 8 bytes, i.e., double precision.          */
!     /*--------------------------------------------------------*/

      IMPLICIT REAL*8 (A-H, O-Z)

!     /*--------------------------------------------------------*/
!     /* Parameter statements.                                  */
!     /*--------------------------------------------------------*/

	integer*4 MXNANG, MXNWORK, JX,LL,IT,IT2

      PARAMETER(MXNANG=501, MXNWORK=500000)

!     /*--------------------------------------------------------*/
!     /* Dimension statements.                                  */
!     /*--------------------------------------------------------*/

      REAL*8 ANGLESc(*),S1R(*),S1C(*),S2R(*),S2C(*)
      REAL*8 S11N,S11(*),S12(*),S33(*),S34(*),SPOL(*),SP(*)

!     /*--------------------------------------------------------*/
!     /* Define the types of the common block.                  */
!     /*--------------------------------------------------------*/

      INTEGER*4 IPHASEmie,NSCATH

      REAL*8 ALAMB, RGmin, RGmax, RGV, SIGMAG, RGCFRAC, RFRS, &
           RFIS, RFRC, RFIC

      complex*16 an(500),bn(500)
      integer*4 ntrm

!     /*--------------------------------------------------------*/
!     /* Set reals to 8 bytes, i.e., double precision.          */
!     /*--------------------------------------------------------*/

!      IMPLICIT REAL*8 (A-H, O-Z)

!     /*--------------------------------------------------------*/
!     /* Input common block for scattering calculations.        */
!     /*--------------------------------------------------------*/

!jdf  COMMON / PHASE  / IPHASEmie

!jdf  COMMON / INPUTS / ALAMB, RGmin, RGmax, RGV, SIGMAG, &
!jdf                   RGCFRAC, RFRS,  RFIS,  RFRC,  RFIC

!     /*--------------------------------------------------------*/
!     /* Output common block for scattering calculations.       */
!     /*--------------------------------------------------------*/

!jdf  COMMON / OUTPUTS / QEXT, QSCAT, QBS, EXT, SCAT, BSCAT, ASY

!     /*--------------------------------------------------------*/
!     /* Arrays to perform the scattering calculations and to   */
!     /* hold the subsequent results.                           */
!     /*--------------------------------------------------------*/

      REAL*8  THETA(MXNANG), ELTRMX(4,MXNANG,2), PII(3,MXNANG), &
             TAU(3,MXNANG), CSTHT(MXNANG),      SI2THT(MXNANG)

      REAL*8   ROUT, RFRO, RFIO, DQEXT, DQSCAT, CTBRQS, DQBS, &
              RIN,  RFRI, RFII, WNUM
 
      COMPLEX*16  ACAP(MXNWORK)
 
      REAL*8  COSPHI(2*MXNANG-1), SCTPHS(2*MXNANG-1)

      INTEGER J, JJ, NINDEX, NSCATA

!     /*--------------------------------------------------------*/
!     /* Obvious variable initializations.                      */
!     /*--------------------------------------------------------*/

         PIE    = DACOS( -1.0D0 )

!     /*--------------------------------------------------------*/
!     /* Maximum number of scattering angles between 0 and 90   */
!     /* degrees, inclusive.                                    */
!     /*--------------------------------------------------------*/

         IT = MXNANG

!     /*--------------------------------------------------------*/
!     /* Maximum number of scattering angles between 0 and 180  */
!     /* degrees, inclusive.                                    */
!     /*--------------------------------------------------------*/

         IT2 = 2 * IT - 1

!     /*--------------------------------------------------------*/
!     /* Dimension of the work array ACAP.                      */
!     /*--------------------------------------------------------*/

         LL = MXNWORK
 
!     /*--------------------------------------------------------*/
!     /* NSCATA is the actual user-requested number of          */
!     /* scattering angles between 0 and 90 degrees, inclusive. */
!     /*--------------------------------------------------------*/

         NSCATA = 2 * NSCATH - 1

!     /*--------------------------------------------------------*/
!     /* If the user did not request a phase function, then we  */
!     /* can set NSCATA and NSCATH to 0.                        */
!     /*--------------------------------------------------------*/

         IF ( IPHASEmie  .le.  0 )  then
              NSCATH  =  0
              NSCATA  =  0
         ENDIF

!     /*--------------------------------------------------------*/
!     /* Check to make sure that the user-requested number of   */
!     /* scattering angles does not excede the current maximum  */
!     /* limit.                                                 */
!     /*--------------------------------------------------------*/

         IF ( NSCATA .gt. IT2  .OR.  NSCATH .gt. IT)  then
              WRITE( 6,105 )  NSCATA, NSCATH, IT2, IT
              STOP 11
         ENDIF

!     /*--------------------------------------------------------*/
!     /* Subroutine SCATANGLES was added by EEC[0495] in order  */
!     /* to facilitate changing the scattering angle locations  */
!     /* output by the Ackerman and Toon Mie code.              */
!     /*--------------------------------------------------------*/

!         CALL SCATANGLES(NSCATH,THETA,COSPHI)

!     /*--------------------------------------------------------*/
!     /* COMPUTE SCATTERING PROPERTIES OF THE PARTICLE.         */
!     /*--------------------------------------------------------*/

!     /*--------------------------------------------------------*/
!     /* DMIESS expects a wavenumber.                           */
!     /*--------------------------------------------------------*/

      WNUM = (2.D0*PIE) / ALAMB

!     /*--------------------------------------------------------*/
!     /* DMIESS assignments of the indices of refraction of the */
!     /* core and shell materials.                              */
!     /*--------------------------------------------------------*/

      RFRO = RFRS
      RFIO = RFIS
      RFRI = RFRC
      RFII = RFIC

!     /*--------------------------------------------------------*/
!     /* DMIESS core and shell radii.                           */
!     /*--------------------------------------------------------*/

       ROUT = RGV
       RIN  = RGCFRAC * ROUT

!     /*--------------------------------------------------------*/
!     /* Scattering angles are symmetric about 90 degrees.      */
!     /*--------------------------------------------------------*/

      IF ( NSCATH  .eq.  0.0 )  THEN
           JX  =  1
      ELSE
           JX   = NSCATH
      ENDIF

!     /*--------------------------------------------------------*/
!     /* Compute the scattering properties for this particle.   */
!     /*--------------------------------------------------------*/

      CALL DMIESS(  ROUT,    RFRO,    RFIO,    THETA,   JX,  &
                   DQEXT,   DQSCAT,  CTBRQS,  ELTRMX,  PII, &
                   TAU,     CSTHT,   SI2THT,  ACAP,    DQBS,  IT, &
                   LL,      RIN,     RFRI,    RFII,    WNUM, an, bn, ntrm   )  ! jcb
 
!     /*--------------------------------------------------------*/
!     /* Compute total cross-sectional area of the particle.    */
!     /*--------------------------------------------------------*/

      X = PIE * RGV * RGV

!     /*--------------------------------------------------------*/
!     /* Assign the final extinction efficiency.                */
!     /*--------------------------------------------------------*/

      QEXT = DQEXT

!     /*--------------------------------------------------------*/
!     /* Compute total extinction cross-section due to particle.*/
!     /*--------------------------------------------------------*/

      EXT = DQEXT * X

!     /*--------------------------------------------------------*/
!     /* Assign the final scattering efficiency.                */
!     /*--------------------------------------------------------*/

      QSCAT = DQSCAT

!     /*--------------------------------------------------------*/
!     /* Compute total scattering cross-section due to particle.*/
!     /*--------------------------------------------------------*/

      SCAT = DQSCAT * X

!     /*--------------------------------------------------------*/
!     /* Assign the final backscatter efficiency.               */
!     /*--------------------------------------------------------*/

      QBS = DQBS

!     /*--------------------------------------------------------*/
!     /* Compute backscatter due to particle.                   */
!     /*--------------------------------------------------------*/

      BSCAT = DQBS * X

!     /*--------------------------------------------------------*/
!     /* Compute asymmetry parameter due to particle.           */
!     /*--------------------------------------------------------*/

      ASY = (CTBRQS * X) / SCAT

!     /*--------------------------------------------------------*/
!     /* If IPHASE is 1, compute the phase function.            */
!     /* S33 and S34 matrix elements are normalized by S11. S11 */
!     /* is normalized to 1.0 in the forward direction.  The    */
!     /* variable SPOL is the degree of polarization for        */
!     /* incident unpolarized light.                            */
!     /*--------------------------------------------------------*/

      IF ( IPHASEmie  .gt.  0 )  THEN

        DO 355 J=1,NSCATA

          IF (J .LE. JX)  THEN
             JJ = J
             NINDEX = 1
          ELSE
             JJ = NSCATA - J + 1
             NINDEX = 2
          ENDIF

          ANGLESc(J) = COSPHI(J)

          S1R(J) = ELTRMX(1,JJ,NINDEX)
          S1C(J) = ELTRMX(2,JJ,NINDEX)
          S2R(J) = ELTRMX(3,JJ,NINDEX)
          S2C(J) = ELTRMX(4,JJ,NINDEX)

          S11(J) = 0.5D0*(S1R(J)**2+S1C(J)**2+S2R(J)**2+S2C(J)**2)
          S12(J) = 0.5D0*(S2R(J)**2+S2C(J)**2-S1R(J)**2-S1C(J)**2)
          S33(J) = S2R(J)*S1R(J) + S2C(J)*S1C(J)
          S34(J) = S2R(J)*S1C(J) - S1R(J)*S2C(J)

          SPOL(J) = -S12(J) / S11(J)

          SP(J) = (4.D0*PIE)*(S11(J) / (SCAT*WNUM**2))

  355   CONTINUE

!        /*-----------------------------------------------------*/
!        /* DONE with the phase function so exit the IF.        */
!        /*-----------------------------------------------------*/

      ENDIF

!     /*--------------------------------------------------------*/
!     /* END of the computations so exit the routine.           */
!     /*--------------------------------------------------------*/

      RETURN

!     /*--------------------------------------------------------*/
!     /* FORMAT statements.                                     */
!     /*--------------------------------------------------------*/

  100 FORMAT ( 7X, I3 )
  105 FORMAT ( ///, 1X,'NUMBER OF ANGLES SPECIFIED =',2I6, / &
                  10X,'EXCEEDS ARRAY DIMENSIONS =',2I6 )

  120 FORMAT (/10X,'INTEGRATED VOLUME',           T40,'=',1PE14.5,/ &
              15X,'PERCENT VOLUME IN CORE',      T40,'=',0PF10.5,/ &
              15X,'PERCENT VOLUME IN SHELL',     T40,'=',0PF10.5,/ &
              10X,'INTEGRATED SURFACE AREA',     T40,'=',1PE14.5,/ &
              10X,'INTEGRATED NUMBER DENSITY',   T40,'=',1PE14.5 )
  125 FORMAT ( 10X,'CORE RADIUS COMPUTED FROM :', /, 20X, 9A8, /  )

  150 FORMAT ( ///,1X,'* * * WARNING * * *', / &
              10X,'PHASE FUNCTION CALCULATION MAY NOT HAVE CONVERGED'/ &
              10X,'VALUES OF S1 AT NSDI-1 AND NSDI ARE :', 2E14.6, /  &
              10X,'VALUE OF X AT NSDI =', E14.6 )
 
!     /*--------------------------------------------------------*/
!     /* DONE with this subroutine so exit.                     */
!     /*--------------------------------------------------------*/

      END SUBROUTINE PFCNPARTICLE

! /*****************************************************************/
! /*****************************************************************/
	subroutine  lpcoefjcb ( ntrm, nmom, ipolzn, momdim,a, b, pmom )
!
!         calculate legendre polynomial expansion coefficients (also
!         called moments) for phase quantities ( ref. 5 formulation )
!
!     input:  ntrm                    number terms in mie series
!             nmom, ipolzn, momdim    'miev0' arguments
!             a, b                    mie series coefficients
!
!     output: pmom                   legendre moments ('miev0' argument)
!
!     *** notes ***
!
!         (1)  eqs. 2-5 are in error in dave, appl. opt. 9,
!         1888 (1970).  eq. 2 refers to m1, not m2;  eq. 3 refers to
!         m2, not m1.  in eqs. 4 and 5, the subscripts on the second
!         term in square brackets should be interchanged.
!
!         (2)  the general-case logic in this subroutine works correctly
!         in the two-term mie series case, but subroutine  'lpco2t'
!         is called instead, for speed.
!
!         (3)  subroutine  'lpco1t', to do the one-term case, is never
!         called within the context of 'miev0', but is included for
!         complete generality.
!
!         (4)  some improvement in speed is obtainable by combining the
!         310- and 410-loops, if moments for both the third and fourth
!         phase quantities are desired, because the third phase quantity
!         is the real part of a complex series, while the fourth phase
!         quantity is the imaginary part of that very same series.  but
!         most users are not interested in the fourth phase quantity,
!         which is related to circular polarization, so the present
!         scheme is usually more efficient.
!
      implicit none
      integer  ipolzn, momdim, nmom, ntrm
      real*8    pmom( 0:momdim,1 )  ! the ",1" dimension is for historical reasons
      complex*16  a(500), b(500)
!
!           ** specification of local variables
!
!      am(m)       numerical coefficients  a-sub-m-super-l
!                     in dave, eqs. 1-15, as simplified in ref. 5.
!
!      bi(i)       numerical coefficients  b-sub-i-super-l
!                     in dave, eqs. 1-15, as simplified in ref. 5.
!
!      bidel(i)    1/2 bi(i) times factor capital-del in dave
!
!      cm,dm()     arrays c and d in dave, eqs. 16-17 (mueller form),
!                     calculated using recurrence derived in ref. 5
!
!      cs,ds()     arrays c and d in ref. 4, eqs. a5-a6 (sekera form),
!                     calculated using recurrence derived in ref. 5
!
!      c,d()       either -cm,dm- or -cs,ds-, depending on -ipolzn-
!
!      evenl       true for even-numbered moments;  false otherwise
!
!      idel        1 + little-del  in dave
!
!      maxtrm      max. no. of terms in mie series
!
!      maxmom      max. no. of non-zero moments
!
!      nummom      number of non-zero moments
!
!      recip(k)    1 / k
!
      integer maxtrm,maxmom,mxmom2,maxrcp
      parameter  ( maxtrm = 1102, maxmom = 2*maxtrm, mxmom2 = maxmom/2, maxrcp = 4*maxtrm + 2 )
      real*8      am( 0:maxtrm ), bi( 0:mxmom2 ), bidel( 0:mxmom2 ), recip( maxrcp )
      complex*16 cm( maxtrm ), dm( maxtrm ), cs( maxtrm ), ds( maxtrm )
      integer k,j,l,nummom,ld2,idel,m,i,mmax,imax
      real*8 sum
      logical    pass1, evenl
      save  pass1, recip
      data  pass1 / .true. /
!
!
      if ( pass1 )  then
!
         do  1  k = 1, maxrcp
            recip( k ) = 1.0 / k
    1    continue
         pass1 = .false.
!
      end if
!
         do  l = 0, nmom
            pmom( l, 1 ) = 0.0
	 enddo
! these will never be called
!      if ( ntrm.eq.1 )  then
!         call  lpco1t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
!         return
!      else if ( ntrm.eq.2 )  then
!         call  lpco2t ( nmom, ipolzn, momdim, calcmo, a, b, pmom )
!         return
!      end if
!
      if ( ntrm+2 .gt. maxtrm ) &
          write(6,1010)
1010    format( ' lpcoef--parameter maxtrm too small' )
!
!                                     ** calculate mueller c, d arrays
      cm( ntrm+2 ) = ( 0., 0. )
      dm( ntrm+2 ) = ( 0., 0. )
      cm( ntrm+1 ) = ( 1. - recip( ntrm+1 ) ) * b( ntrm )
      dm( ntrm+1 ) = ( 1. - recip( ntrm+1 ) ) * a( ntrm )
      cm( ntrm ) = ( recip(ntrm) + recip(ntrm+1) ) * a( ntrm ) &
                  + ( 1. - recip(ntrm) ) * b( ntrm-1 )
      dm( ntrm ) = ( recip(ntrm) + recip(ntrm+1) ) * b( ntrm ) &
                  + ( 1. - recip(ntrm) ) * a( ntrm-1 )
!
      do  10  k = ntrm-1, 2, -1
         cm( k ) = cm( k+2 ) - ( 1. + recip(k+1) ) * b( k+1 ) &
                            + ( recip(k) + recip(k+1) ) * a( k ) &
                            + ( 1. - recip(k) ) * b( k-1 )
         dm( k ) = dm( k+2 ) - ( 1. + recip(k+1) ) * a( k+1 ) &
                            + ( recip(k) + recip(k+1) ) * b( k ) &
                            + ( 1. - recip(k) ) * a( k-1 )
   10 continue
      cm( 1 ) = cm( 3 ) + 1.5 * ( a( 1 ) - b( 2 ) )
      dm( 1 ) = dm( 3 ) + 1.5 * ( b( 1 ) - a( 2 ) )
!
      if ( ipolzn.ge.0 )  then
!
         do  20  k = 1, ntrm + 2
            cm( k ) = ( 2*k - 1 ) * cm( k )
            dm( k ) = ( 2*k - 1 ) * dm( k )
   20    continue
!
      else
!                                    ** compute sekera c and d arrays
         cs( ntrm+2 ) = ( 0., 0. )
         ds( ntrm+2 ) = ( 0., 0. )
         cs( ntrm+1 ) = ( 0., 0. )
         ds( ntrm+1 ) = ( 0., 0. )
!
         do  30  k = ntrm, 1, -1
            cs( k ) = cs( k+2 ) + ( 2*k + 1 ) * ( cm( k+1 ) - b( k ) )
            ds( k ) = ds( k+2 ) + ( 2*k + 1 ) * ( dm( k+1 ) - a( k ) )
   30    continue
!
         do  40  k = 1, ntrm + 2
            cm( k ) = ( 2*k - 1 ) * cs( k )
            dm( k ) = ( 2*k - 1 ) * ds( k )
   40    continue
!
      end if
!
!
      if( ipolzn.lt.0 )  nummom = min0( nmom, 2*ntrm - 2 )
      if( ipolzn.ge.0 )  nummom = min0( nmom, 2*ntrm )
      if ( nummom .gt. maxmom ) &
          write(6,1020)
1020	format( ' lpcoef--parameter maxtrm too small')
!
!                               ** loop over moments
      do  500  l = 0, nummom
         ld2 = l / 2
         evenl = mod( l,2 ) .eq. 0
!                                    ** calculate numerical coefficients
!                                    ** a-sub-m and b-sub-i in dave
!                                    ** double-sums for moments
         if( l.eq.0 )  then
!
            idel = 1
            do  60  m = 0, ntrm
               am( m ) = 2.0 * recip( 2*m + 1 )
   60       continue
            bi( 0 ) = 1.0
!
         else if( evenl )  then
!
            idel = 1
            do  70  m = ld2, ntrm
               am( m ) = ( 1. + recip( 2*m-l+1 ) ) * am( m )
   70       continue
            do  75  i = 0, ld2-1
               bi( i ) = ( 1. - recip( l-2*i ) ) * bi( i )
   75       continue
            bi( ld2 ) = ( 2. - recip( l ) ) * bi( ld2-1 )
!
         else
!
            idel = 2
            do  80  m = ld2, ntrm
               am( m ) = ( 1. - recip( 2*m+l+2 ) ) * am( m )
   80       continue
            do  85  i = 0, ld2
               bi( i ) = ( 1. - recip( l+2*i+1 ) ) * bi( i )
   85       continue
!
         end if
!                                     ** establish upper limits for sums
!                                     ** and incorporate factor capital-
!                                     ** del into b-sub-i
         mmax = ntrm - idel
         if( ipolzn.ge.0 )  mmax = mmax + 1
         imax = min0( ld2, mmax - ld2 )
         if( imax.lt.0 )  go to 600
         do  90  i = 0, imax
            bidel( i ) = bi( i )
   90    continue
         if( evenl )  bidel( 0 ) = 0.5 * bidel( 0 )
!
!                                    ** perform double sums just for
!                                    ** phase quantities desired by user
         if( ipolzn.eq.0 )  then
!
            do  110  i = 0, imax
!                                           ** vectorizable loop (cray)
               sum = 0.0
               do  100  m = ld2, mmax - i
                  sum = sum + am( m ) * &
                           ( dble( cm(m-i+1) * dconjg( cm(m+i+idel) ) ) &
                           + dble( dm(m-i+1) * dconjg( dm(m+i+idel) ) ) )
  100          continue
               pmom( l,1 ) = pmom( l,1 ) + bidel( i ) * sum
  110       continue
            pmom( l,1 ) = 0.5 * pmom( l,1 )
            go to 500
!
         end if
!
  500 continue
!
!
  600 return
      end subroutine lpcoefjcb
!
	end module module_optical_averaging