  subroutine rdtisdale(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
! Notice that the exact number of elements in rnval,rival works
  integer,parameter :: ndat=1685,nwt=8
  integer :: i,j1,j2,mn
  integer :: start(1),count(1),stride(1)
  integer :: start2(2),count2(2),stride2(2)
  real :: diff,diffmin
! Note that the ndat,nwt is opposite to the ncdump dump
  real :: wcmdat(ndat),wavedat(ndat)
  real :: rnval(ndat,nwt),rival(ndat,nwt)
  real :: h2so4wt(nwt),wh2so4
  real :: dh,dh2so4,derivn,derivi
  character(len=30) :: fil
  character (len=16) :: llset

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

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

! Obtain the full pathname of the input netCDF file
   noprf=1
  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,"rdtisdale: id1,filnm ",i8,2x,a100)
  end if

! ***
! Get dimension ids
  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,'nlines',idnlines)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'nweights',idnweights)
   call handle(iout,istat)

! *****
! Get the dimensions
  istat=nf90_get_var(id1,idnlines,nlines)
   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,idnweights,nweights)
   call handle(iout,istat)

! *****
! Get variable ids
  istat=nf90_inq_varid(id1,'h2so4wt',idh2so4wt)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'wcm',idwcm)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'wavelength',idwave)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'rn8',idrn)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'ri8',idri)
   call handle(iout,istat)

  if (nopr .eq. 1) then
   write(iout,fmt=500) nlines,nweights,ntext,n80
   500 format(/,2x,"nlines,nweights,ntext,n80",/,2x,4(2x,i4))
   write(iout,fmt=510) idh2so4wt,idwcm,idwave,idrn,idri
   510 format(/,2x,"idh2so4wt,idwcm,idwave,idrn,idri",/,2x,5(2x,i4))
  end if

! ****
! The weight percent h2so4 values
   start(1)=1
   count(1)=nwt
   stride(1)=1
  istat=nf90_get_var(id1,idh2so4wt,h2so4wt,start,count,stride)

! ****
! Read in the wcmdat and wavedat values
   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)

!  Read in the real and imaginary indices
! Note that count2(1) and count2(2) are opposite to the ncdump dump
   start2(1)=1
   start2(2)=1
   count2(1)=nlines
   count2(2)=nwt
   stride2(1)=1
   stride2(2)=1
  istat=nf90_get_var(id1,idrn,rnval,start2,count2,stride2)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idri,rival,start2,count2,stride2)
   call handle(iout,istat)

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

! Close the data set
  istat=nf90_close(id1)

! ****
! Specify the weight percent h2so4
   write(0,fmt=130)
   130 format(/,2x,"The weight percent h2so4 values to choose from")
  do i=1,nwt
   write(0,fmt=135) h2so4wt(i)
   135 format(2x,f10.4)
  end do

  write(0,fmt=160)
  160 format(2x,"Specify wh2so4 ")

  read(5,*) wh2so4

! The compound string
!       '1234567890123456'
  llset=' Tisdale H2SO4  '

  write(0,fmt=170) wh2so4
  170 format(2x," wh2so4 ",f10.4)

! *****
! Find the indices for the interpolation
   j1=-99
   j2=-99
    dh=0.0
   if (wh2so4 .lt. h2so4wt(1)) then
    j1=1
    j2=2
    dh=wh2so4-h2so4wt(1)
    go to 600 
   end if
   if (wh2so4 .gt. h2so4wt(nwt)) then
    j1=nwt-1
    j2=nwt
    dh=wh2so4-h2so4wt(1)
    go to 600 
   end if
  do i=1,nwt-1
    i2=i+1
   if ((wh2so4 .ge. h2so4wt(i)) .and. (wh2so4 .lt. h2so4wt(i2))) then
    j1=i
    j2=i2
    dh=wh2so4-h2so4wt(i)
    go to 600 
   end if
  end do
  600 mn=0

! Check
  if ((j1 .le. 0) .or. (j2 .le. 0)) then 
   write(0,fmt=610)
   610 format(2x,"j1 or j2 le 0, will stop")
   stop
  end if

! Put values into the output arrays
   dh2so4=h2so4wt(j2)-h2so4wt(j1)
  do i=1,nlines
   wcm(i)=wcmdat(i)
   wavelength(i)=wavedat(i)
   derivn=(rnval(i,j2)-rnval(i,j1))/dh2so4
   derivi=(rival(i,j2)-rival(i,j1))/dh2so4
   a1=rnval(i,j1)+(dh*derivn)
   a2=rival(i,j1)+(dh*derivi)
   if (a1 .le. 0.0) then 
    a1=0.0
   endif
   if (a2 .le. 0.0) then 
    a2=0.0
   endif
   rndat(i)=a1
   ridat(i)=a2
  end do

! *****
  if (nopr .eq. 1) then
    write(iout,fmt=840) wh2so4,j1,j2,h2so4wt(j1),h2so4wt(j2),dh
    840 format(/,2x,"rdtisdale:  wh2so4,j1,j2,h2so4wt(j1),h2so4wt(j2),dh ",/,&
    2x,f10.4,2(1x,i3),3(1x,f10.4),/,&
    2x,"rdtisdale: 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,f12.4),1p,2(1x,e10.3))
   end do
  end if

! ******
  return
  end

