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

! ****
! 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,parameter :: ndat=180,ncomp=10
  integer :: id1,istat,noprf,mset
  integer :: start(1),count(1),stride(1)
  real :: wcmdat(ndat),wavedat(ndat),rnval(ndat),rival(ndat)
  character(len=30) :: fil
  character (len=16) :: files(ncomp),llset
  character (len=16) :: llsets(ncomp)
  character (len=30) :: filei(ncomp)

! **
! The different compositions
!          '1234567890123456'
  files(1)='SiO2 300K ab    '
  files(2)='SiO2 300K c     '
  files(3)='SiO2 551K ab    '
  files(4)='SiO2 551K c     '
  files(5)='SiO2 738K ab    '
  files(6)='SiO2 738K c     '
  files(7)='SiO2 833K ab    '
  files(8)='SiO2 833K c     '
  files(9)='SiO2 928K ab    '
  files(10)='SiO2 928K c     '

! Input netCDF files to choose from
!          '123456789012345678901234567890'
  filei(1)='zeidler_sio2_300k_ab.nc       '
  filei(2)='zeidler_sio2_300k_c.nc        '
  filei(3)='zeidler_sio2_551k_ab.nc       '
  filei(4)='zeidler_sio2_551k_c.nc        '
  filei(5)='zeidler_sio2_738k_ab.nc       '
  filei(6)='zeidler_sio2_738k_c.nc        '
  filei(7)='zeidler_sio2_833k_ab.nc       '
  filei(8)='zeidler_sio2_833k_c.nc        '
  filei(9)='zeidler_sio2_928k_ab.nc       '
  filei(10)='zeidler_sio2_928k_c.nc        '

!            '1234567890123456'
   llsets(1)='sio2 300K ab    '
   llsets(2)='sio2 300K c     ' 
   llsets(3)='sio2 551K ab    '
   llsets(4)='sio2 551K c     '
   llsets(5)='sio2 738K ab    '
   llsets(6)='sio2 738K c     '
   llsets(7)='sio2 833K ab    '
   llsets(8)='sio2 833K c     '
   llsets(9)='sio2 928K ab    '
   llsets(10)='sio2 928K c     '

! **
! User will select which set to work with
  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
   write(0,fmt=150) i,files(i)
   150 format(2x,i3,2x,a16)
  end do

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

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

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

! ***
! The input netCDF file
   fil=filei(mset)
   llset=llsets(mset)

   noprf=1
  call getfilnm(noprf,iout,fil,indir,filnm)
  if (noprf .eq. 1) then  
!  stop
  end if

! *****
! Open the netCDF file
   ifile=1
  istat=nf90_open(filnm,nf90_nowrite,id1)
  call handle(iout,istat)
 
! ***
! Get dimension ids
  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,idnlines,nlines)
   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)

  if (nopr .eq. 1) then
   write(iout,fmt=500) nlines,ntext,n70
   500 format(/,2x,"nlines ",i4,/,2x,"ntext ",i4,/,2x,"n70 ",i4)
  end if

! ***
! Get variable ids 
  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,'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=510) idwcm,idwave,idrn,idri
   510 format(/,2x," idwcm,idwave,idrn,idri ",/,2x,4(1x,i4))
  end if

! ***
! Get data
   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=200)
   200 format(/)
   write(iout,fmt=210)
   210 format(2x,"rdzeidler_sio2_hghtemp: text(i)")
   do i=1,ntext
!   write(iout,fmt=*) text(i)
   end do

   write(iout,fmt=200)
   write(iout,fmt=230)
   230 format(2x,"rdzeidler_sio2_hghtemp: i,wcm,wavelength,rndat and ridat")
   do i=1,nlines
    write(iout,fmt=240) i,wcm(i),wavelength(i),rndat(i),ridat(i)
    240 format(2x,i4,1p,2(1x,e10.3),2x,2(1x,e10.3))
   end do

  end if

! ******
  return
  end
