subroutine calcext(nopr,iout,&
  iwave,nwave,wave,rn,ri,&
  ndist,radr,sized,&
  bext,babs,bsca,asym,back,omega,&
  outdir,iwrext,lset)

! *************************
  integer :: nopr,iout,iwrext
  character(len=100) :: lset

! See calcwave.f90
  integer :: nwave
  integer :: iwave
  real :: wave(nwavemax)
  real :: rn(nwavemax),ri(nwavemax)

! See calcsized.f90
  integer,parameter :: ndist2=200
  integer :: ndist
  real :: radr(ndist2),sized(ndist2)

! See calcext.f90
  real :: bext(nwavemax),babs(nwavemax)
  real :: bsca(nwavemax),asym(nwavemax)
  real :: omega(nwavemax),back(nwavemax)

! See rdinoutdir.f90
  character(len=60) :: outdir

! Used here
  integer,parameter :: nang=20
  complex :: refrel,s1(200),s2(200)
  real :: dr(ndist2)

! *************************
   pi=3.14159265
   const=sqrt(2.00*pi)
   constx=2.00*pi*1.0e-4
   refmed=1.00

! ****
! For the size distribution
   volume=0.0
   area=0.0
   rave=0.0
   totden=0.0

! ****
! The radii increments
   do i=1,ndist 
    i1=i
    i2=i+1
    dr(i)=radr(i2)-radr(i1)
   end do
    i3=ndist
    i2=ndist-1
   dr(i3)=dr(i2)

! ****
! Loop over spectra grid
  do i=1,nwave 

! Initialize for each spectra grid point
   bext(i)=0.0
   babs(i)=0.0
   bsca(i)=0.0
   asym(i)=0.0
   back(i)=0.0
   omega(i)=0.0

! For wavenumbers
   if (iwave .eq. 1) then
    wavcm=wave(i)
   end if
! For wavelength in microns, calculate wavcm
   if (iwave .eq. 2) then
    wavcm=1.0e4/wave(i)
   end if

! The complex index of refraction
   refre=rn(i)
   refim=ri(i)
!  refrel=(cmpl(refre,refim))/refmed
   refrel=cmplx(refre,refim)

! Loop over the size distribution
   do j=1,ndist 

! The Mie x parameter (2 pi Radius / Wavelength)
   x=constx*radr(j)*wavcm

! Use the Bohren and Huffman BHMIE routine
   call bhmie(x,refrel,nang,s1,s2,qext,qsca,qback,gfac)
 
! The Q absorption efficieny factor 
   qabs=qext-qsca

! Calculate the extinction, scattering, absorption coefficient
! 1.0e-3 converts cm-1 to km-1
    rd2=radr(j)*radr(j)
    weight=pi*rd2*sized(j)*dr(j)*1.0e-3

! Add to the output arrays
   bext(i)=bext(i)+(weight*qext)
   babs(i)=babs(i)+(weight*qabs)
   bsca(i)=bsca(i)+(weight*qsca)
   back(i)=back(i)+(weight*qback)
   asym(i)=asym(i)+(weight*qsca*gfac)

  end do
! Loop over the size distribution is done

! ***
! The asymmetry factor
   asym(i)=asym(i)/bsca(i)

! The single scattering albedo
   omega(i)=bsca(i)/bext(i)

  end do
! Loop over spectra grid points is done

! *************************
  if (nopr .eq. 1) then 
   write(iout,fmt=500)
   500 format(/,2x,"calcext: bext,babs,bsca, are in 1/km units ",/,&
   2x,"calcext: asym,omega are unitless ",/,&
   2x,"calcext: i,wave(i),bext(i),babs(i),bsca(i),asym(i),omega(i) ")
   do i=1,nwave 
    write(iout,fmt=510) i,wave(i),bext(i),babs(i),bsca(i),asym(i),omega(i)
    510 format(1(1x,i3),1(1x,f10.4),1p,5(1x,e10.3))
   end do
  end if

! *************************
! Write the spectra out to ascii and netCDF files
   if (iwrext .eq. 1) then 
    call wrext(nopr,iout,outdir,iwave,&
    nwave,wave,bext,babs,bsca,asym,omega)
   end if

! *************************
  return
  end