  subroutine rdniedziela_h2so4(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=2010,mtemp=8,nwt=6
  integer :: i,i1,i2,i3,i4,j1,j2,k1,k2,mn,nc
  integer :: ntemps,nnlines,ntext,n70,nwh2so4
  integer :: start(1),count(1),stride(1)
  integer :: start2(2),count2(2),stride2(2)
  integer :: start3(3),count3(3),stride3(3)
  real :: a1,a2
! Note that the ndat,mtemp is opposite to the ncdump dump
  real :: temps(mtemp),temp
  real :: numlines(mtemp)
  real :: wh2so4(mtemp,nwt),wghth2so4,ncomps(mtemp)
  real :: wcmdat(ndat,mtemp),wavel(ndat,mtemp)
  real :: rnval(ndat,mtemp,nwt),rival(ndat,mtemp,nwt)
  character(len=30) :: fil
  character (len=16) :: llset

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

! ***
! The input netCDF file
!     '123456789012345678901234567890'
  fil='niedziela_h2so4.nc            '
   noprf=1
  call getfilnm(noprf,iout,fil,indir,filnm)
  if (noprf .eq. 1) then
!  stop
  end if

! ***
! Specify the weight percent h2so4
  write(0,fmt=140)
  140 format(2x,"Specify temperature (e.g. 200 to 300 K) ")
  read(5,*) temp

  write(0,fmt=165) temp
  165 format(2x," temp ",f10.4)

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

! ***
! 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,"rdniedziela_h2so4: id1,filnm ",i8,2x,a100)
  end if

! ***
! Get dimension ids
  istat=nf90_inq_varid(id1,'ntemps',idntemps)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'nwh2so4',idnwh2so4)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'nlines',idnlines)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'ntext',idntext)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'n70',idn70)
   call handle(iout,istat)

! *****
! Get the dimensions
  istat=nf90_get_var(id1,idntemps,ntemps)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idnwh2so4,nwh2so4)
   call handle(iout,istat)
! notice that you use nnlines here
  istat=nf90_get_var(id1,idnlines,nnlines)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idntext,ntext)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idn70,n70)
   call handle(iout,istat)

! *****
! Get variable ids
  istat=nf90_inq_varid(id1,'temps',idtemps)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'numlines',idnumlines)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'wh2so4',idwh2so4)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'ncomps',idncomps)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'wcm',idwcm)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'wavelength',idwavel)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'rn',idrn)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'ri',idri)
   call handle(iout,istat)

  if (nopr .eq. 1) then
   write(iout,fmt=500) ntemps,nwh2so4,nnlines,ntext,n70
   500 format(/,2x,"ntemps,nwh2so4,nnlines,ntext,n70 ",/,2x,5(2x,i4))
   write(iout,fmt=510) idntemps,idnumlines,idwh2so4,idwcm,idwavel
   510 format(/,2x,"idntemps,idnumlines,idwh2so4,idwcm,idwavel ",/,&
   2x,5(1x,i4))
   write(iout,fmt=520) idrn,idri
   520 format(/,2x,"idnrn,idri ",/,2x,2(1x,i4))
  end if

! ****
! The temperatures
   start(1)=1
   count(1)=mtemp
   stride(1)=1
  istat=nf90_get_var(id1,idtemps,temps,start,count,stride)

! ****
! The number of lines for each of the ntemps temperatures
   start(1)=1
   count(1)=mtemp
   stride(1)=1
  istat=nf90_get_var(id1,idnumlines,numlines,start,count,stride)

! ****
! The weight percent h2so4 values
   start2(1)=1
   start2(2)=1
   count2(1)=mtemp
   count2(2)=nwt
   stride2(1)=1
   stride2(2)=1
  istat=nf90_get_var(id1,idwh2so4,wh2so4,start2,count2,stride2)

! ****
! The number of h2so4 values for each temperature
   start(1)=1
   count(1)=mtemp
   stride(1)=1
  istat=nf90_get_var(id1,idncomps,ncomps,start,count,stride)

! ****
! Read in the wcmdat and wavel values
! notice that you use nnlines here
   start2(1)=1
   start2(2)=1
   count2(1)=nnlines
   count2(2)=mtemp
   stride2(1)=1
   stride2(2)=1
  istat=nf90_get_var(id1,idwcm,wcmdat,start2,count2,stride2)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idwavel,wavel,start2,count2,stride2)
   call handle(iout,istat)

! Read in the real and imaginary indices
! Note that count3(1) and count3(2) etc are opposite to the ncdump dump
! notice that you use nnlines here
   start3(1)=1
   start3(2)=1
   start3(3)=1
   count3(1)=nnlines
   count3(2)=mtemp
   count3(3)=nwt
   stride3(1)=1
   stride3(2)=1
   stride3(3)=1
  istat=nf90_get_var(id1,idrn,rnval,start3,count3,stride3)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idri,rival,start3,count3,stride3)
   call handle(iout,istat)

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

! Close the data set
  istat=nf90_close(id1)

! *****
! Find the indices for the interpolation

! ***
! The temperature value
   k1=-99
  if (temp .lt. temps(1)) then
   k1=1
   go to 590
  end if
  if (temp .ge. temps(mtemp)) then
   k1=mtemp
   go to 590
  end if
  do i=1,mtemp-1
   i3=i
   i4=i+1
   if ((temp .ge. temps(i3)) .and. (temp .lt. temps(i4))) then
    k1=i
    go to 590
   end if
  end do
  590 mn=0

  if (k1 .lt. 0) then
   write(0,fmt=595)
   595 format(2x,"k1 lt 0 for temperature, will stop")
   stop
  end if

! The number of lines to work with
   nlines=numlines(k1)

  write(0,fmt=596) temp,k1,temps(k1),nlines
  596 format(2x,"temp,k1,temps(k1),nlines ",2x,f10.4,2x,i4,2x,f10.4,2x,i4)

! ***
! The number of h2so4 values to look at
    nc=ncomps(k1)
   write(0,fmt=605) nc
   605 format(2x,"nc, number of h2so4 values to look at ",i3)

! ***
! The weight percent h2so4
   j1=-99

   write(0,fmt=899) 
   899 format(2x,"Here are the wh2so4 values to choose from")
  do i=1,nc
   write(0,fmt=900) wh2so4(k1,i)
   900 format(2x,f10.4)
  end do

  write(0,fmt=160)
  160 format(2x,"Specify th wghth2so4 value ")
  read(5,*) wghth2so4

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

  if (wghth2so4 .lt. wh2so4(k1,1)) then
   j1=1
   go to 600 
  end if
  if (wghth2so4 .ge. wh2so4(k1,nc)) then
   j1=nc
   go to 600 
  end if
  do i=1,nc-1
    i3=i
    i4=i+1
    a1=wh2so4(k1,i3)
    a2=wh2so4(k1,i4)
   if ((wghth2so4 .ge. a1) .and. (wghth2so4 .lt. a2)) then
    j1=i
    go to 600 
   end if
  end do
  600 mn=0

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

  write(0,fmt=611) wghth2so4,j1,wh2so4(k1,j1),nlines
  611 format(2x,"wghth2so4,j1,wh2so4(k1,j1),nlines ",/,&
  2x,f10.4,2x,i4,2x,f10.4,2x,i4)

! Put values into the output arrays
  do i=1,nlines
   wcm(i)=wcmdat(i,k1)
   wavelength(i)=wavel(i,k1)
   a1=rnval(i,k1,j1)
   a2=rival(i,k1,j1)
   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=800) temp,wghth2so4,k1,temps(k1),j1,wh2so4(k1,j1)
    800 format(/,2x,"rdniedziela: temp,wghth2so4 ",/,2x,2(1x,f10.4),/,&
    2x,"rdniedziela: k1,temps(k1),j1,wh2so4(k1,j1) ",/,&
    2x,2(1x,i4,2x,f10.4),/,&
    2x,"rdniedziela: 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
