  subroutine rdshettle(nopr,iout,indir,outdir,lstrn,iset,&
  lset,llset,nlines,wcm,wavelength,rndat,ridat,iwave)

! ******
! Read in the data from the netCDF file

! ****
! Will use netCDF f90 routines
  use netcdf

! param1.inc gives the maximum size of indices and extinction vectors
  include 'param1.inc'

! ****
! See rdinoutdir.f90
  character(len=60) :: indir,outdir
  integer :: nopr,iout,iset

! See main.f90
  integer :: nlines,iwave
  real :: wcm(nwavemax),wavelength(nwavemax)
  real :: rndat(nwavemax),ridat(nwavemax)
  character(len=100) :: lstrn,lset,filnm

! Used here
  integer :: icheck,id1,istat,noprf,ilook,mset
  integer :: start(1),count(1),stride(1)
  integer,parameter :: ndat=5000
  real :: wcmdat(ndat),wavedat(ndat),rnval(ndat),rival(ndat)
  character(len=30) :: fil
  character (len=70) :: c70
  character (len=16) :: files(18),llset

! ***
! Read in the data from the netCDF file
   lset=lstrn

! ***
! The input netCDF file
!     '123456789012345678901234567890'
  fil='afcrl1987_shettle.nc          '

! Obtain the full pathname of the input netCDF file
   noprf=0
  call getfilnm(noprf,iout,fil,indir,filnm)
  if (noprf .eq. 1) then
   stop
  end if

! ***
! Open the netCDF file
  istat=nf90_open(filnm,nf90_nowrite,id1)
   call handle(iout,istat)

  if (nopr .eq. 1) then
   write(iout,fmt=100) id1,filnm
   100 format(/,2x,"rdshettle: id1,filnm ",i8,2x,a100)
  end if

! ***
! Get dimension ids (work with variables part of ncdump)
  istat=nf90_inq_varid(id1,'ntext',idntext)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'n80',idn80)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'n16',idn16)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'ncomp',idncomp)
   call handle(iout,istat)

  istat=nf90_inq_varid(id1,'nlines1',idnlines1)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'nlines2',idnlines2)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'nlines3',idnlines3)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'nlines4',idnlines4)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'nlines5',idnlines5)
   call handle(iout,istat)

! *****
! Get the dimensions
  istat=nf90_get_var(id1,idnlines1,nlines1)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idnlines2,nlines2)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idnlines3,nlines3)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idnlines4,nlines4)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idnlines5,nlines5)
   call handle(iout,istat)

  istat=nf90_get_var(id1,idntext,ntext)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idn80,n80)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idn16,n16)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idncomp,ncomp)
   call handle(iout,istat)
  
  if (nopr .eq. 1) then
   write(iout,fmt=110) nlines1,nlines2,nlines3,nlines4,nlines5
   110 format(/,2x,"rdshettle:  nlines1,nlines2,nlines3,nlines4,nlines5 ",/,&
   2x,5(1x,i4))
   write(iout,fmt=120) ntext,n80,n16,ncomp
   120 format(/,2x,"rdshettle:  ntext,n80,n16,ncomp",/,&
   2x,4(1x,i4))
  end if

! ***
! Obtain the listing of the compounds
! Will ask the user to choose which compound to work with

! Specify by hand (have code read in values)
!          '1234567890123456'
  files(1)='Water           '
  files(2)='Ice             '
  files(3)='NaCl            ' 
  files(4)='Seasalt         '
  files(5)='Water_soluble   '
  files(6)='Ammonium        '
  files(7)='Carbonaceous    '
  files(8)='Volcanic_dust   '
  files(9)='H2SO4 215K      '
  files(10)='H2SO4 300K     '
  files(11)='Meteor dust    '
  files(12)='Quartz oray    '
  files(13)='Quartz eray    '
  files(14)='Hematite oray  '
  files(15)='Hematite eray  '
  files(16)='Sand oray      '
  files(17)='Sand eray      '
  files(18)='Dust like      '

! *
! Write out the compounds included in the netCDF file
  write(iout,fmt=130) ncomp
  130 format(/,2x,"rdshettle: ncomp ",i3)

  write(0,fmt=140)
  140 format(2x," Will select mset value for a specific compound")
  write(0,fmt=145)
  145 format(2x," i, compound")

  do i=1,ncomp
   j=i-1
   write(0,fmt=150) j,files(i)
   150 format(2x,i3,2x,a16)
  end do

  write(0,fmt=160)
  160 format(2x,"Specify mset (0-17) e.g. mset=7 for volcanic dust ")

  read(5,*) mset
  mset=int(mset)

! The compound string
  llset=files(mset+1)

  write(0,fmt=170) mset,llset
  170 format(2x," mset,llset ",i4,2x,a16)

! ***
! Original text from Shettle's ascii file
! ncdf_varget,nc_id,idntext,ntext
! ncdf_varget,nc_id,idtext,text
! text=string(text)

! ***
! Obtain the specific indices

  if ((mset .ge. 0) .and. (mset .le. 3)) then 
     nlines=nlines1
    istat=nf90_inq_varid(id1,'wcm1',idwcm)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'wavelength1',idwave)
     call handle(iout,istat)
   if (mset .eq. 0) then 
    istat=nf90_inq_varid(id1,'watern',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'wateri',idri)
     call handle(iout,istat)
   end if
   if (mset .eq. 1) then
    istat=nf90_inq_varid(id1,'icen',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'icei',idri)
     call handle(iout,istat)
   end if
   if (mset .eq. 2) then
    istat=nf90_inq_varid(id1,'ncln',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'ncli',idri)
     call handle(iout,istat)
   end if
   if (mset .eq. 3) then
    istat=nf90_inq_varid(id1,'seasaltn',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'seasalti',idri)
     call handle(iout,istat)
   end if
  end if

  if ((mset .ge. 4) .and. (mset .le. 6)) then
     nlines=nlines2
    istat=nf90_inq_varid(id1,'wcm2',idwcm)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'wavelength2',idwave)
     call handle(iout,istat)
   if (mset .eq. 4) then 
    istat=nf90_inq_varid(id1,'solublen',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'solublei',idri)
     call handle(iout,istat)
   end if
   if (mset .eq. 5) then 
    istat=nf90_inq_varid(id1,'ammoniumn',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'ammoniumi',idri)
     call handle(iout,istat)
   end if
   if (mset .eq. 6) then 
    istat=nf90_inq_varid(id1,'carbonaceousn',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'carbonaceousi',idri)
     call handle(iout,istat)
   end if
  end if

  if ((mset .ge. 7) .and. (mset .le. 10)) then 
     nlines=nlines3
    istat=nf90_inq_varid(id1,'wcm3',idwcm)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'wavelength3',idwave)
     call handle(iout,istat)
   if (mset .eq. 7) then 
    istat=nf90_inq_varid(id1,'volcanic_dustn',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'volcanic_dusti',idri)
     call handle(iout,istat)
   end if
   if (mset .eq. 8) then 
    istat=nf90_inq_varid(id1,'h2so4215n',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'h2so4215i',idri)
     call handle(iout,istat)
   end if
   if (mset .eq. 9) then 
    istat=nf90_inq_varid(id1,'h2so4300n',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'h2so4300i',idri)
     call handle(iout,istat)
   end if
   if (mset .eq. 10) then 
    istat=nf90_inq_varid(id1,'meteor_dustn',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'meteor_dusti',idri)
     call handle(iout,istat)
   end if
  end if

  if ((mset .ge. 11) .and. (mset .le. 14)) then
     nlines=nlines4
    istat=nf90_inq_varid(id1,'wcm4',idwcm)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'wavelength4',idwave)
     call handle(iout,istat)
   if (mset .eq. 11) then
    istat=nf90_inq_varid(id1,'quartz_orayn',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'quartz_orayi',idri)
     call handle(iout,istat)
   end if
   if (mset .eq. 12) then 
    istat=nf90_inq_varid(id1,'quartz_erayn',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'quartz_erayi',idri)
     call handle(iout,istat)
   end if
   if (mset .eq. 13) then 
    istat=nf90_inq_varid(id1,'hematite_orayn',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'hematite_orayi',idri)
     call handle(iout,istat)
   end if
   if (mset .eq. 14) then
    istat=nf90_inq_varid(id1,'hematite_erayn',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'hematite_erayi',idri)
     call handle(iout,istat)
   end if
  end if

  if ((mset .ge. 15) .and. (mset .le. 17)) then 
     nlines=nlines5
    istat=nf90_inq_varid(id1,'wcm5',idwcm)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'wavelength5',idwave)
     call handle(iout,istat)
   if (mset .eq. 15) then
    istat=nf90_inq_varid(id1,'sand_orayn',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'sand_orayi',idri)
     call handle(iout,istat)
   end if
   if (mset .eq. 16) then 
    istat=nf90_inq_varid(id1,'sand_erayn',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'sand_erayi',idri)
     call handle(iout,istat)
   end if
   if (mset .eq. 17) then 
    istat=nf90_inq_varid(id1,'dust_liken',idrn)
     call handle(iout,istat)
    istat=nf90_inq_varid(id1,'dust_likei',idri)
     call handle(iout,istat)
   end if
  end if

! ****
!  Read in the real and imaginary indices
   start(1)=1
   count(1)=nlines
   stride(1)=1
  istat=nf90_get_var(id1,idwcm,wcmdat,start,count,stride)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idwave,wavedat,start,count,stride)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idrn,rnval,start,count,stride)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idri,rival,start,count,stride)
   call handle(iout,istat)

! ****
! Exit define mode
  istat=nf90_enddef(id1)

! Close the data set
  istat=nf90_close(id1)

! *****
! Put values into the output arrays
  do i=1,nlines
   wcm(i)=wcmdat(i)
   wavelength(i)=wavedat(i)
   rndat(i)=rnval(i)
   ridat(i)=rival(i)
  end do

! *****
  if (nopr .eq. 1) then

    write(iout,fmt=800) mset,files(mset)
    800 format(/,2x,"rdshettle: work with mset, files(mset)",/,&
    2x,i3,2x,a16,/,&
    2x,"rdshettle: i,files(i)")
   do i=1,ncomp 
    write(iout,810) i,files(i)
    810 format(2x,i3,2x,a70)
   end do

    iwrtext=0
   if (iwrtext .eq. 1) then 
    write(iout,fmt=820)
    820 format(/,2x,"rdshettle: text(i)")
    do i=1,ntext
     write(iout,fmt=830) text(i)
     830 format(2x,a70)
    end do
   end if

    write(iout,fmt=840)
    840 format(/,2x,"rdshettle: i,wcm,wavelength,rndat,ridat")
   do i=1,nlines 
    write(iout,fmt=850) i,wcm(i),wavelength(i),rndat(i),ridat(i)
    850 format(2x,i4,2(1x,f10.4),1p,2(1x,e10.3))
   end do

  end if

! ******
  return
  end


