  subroutine rdmyhre_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=4733,mtempmax=5,mcases=7
  integer :: i,mn,nc
  integer :: ival,jval
  integer :: ncases,ntext,n70
  integer :: start(1),count(1),stride(1)
  integer :: start2(2),count2(2),stride2(2)
  integer :: start3(3),count3(3),stride3(3)
  real :: a1,a2
  real :: wghth2so4
  real :: vec1(mtempmax),vec2(mtempmax)
  real :: diff1,diff2
! Note that the ndat,mtemp is opposite to the ncdump dump
  real :: numlinesr(mcases),numlinesi(mcases)
  real :: temps(mcases,mtempmax)
  real :: ntemps(mcases)
  real :: wh2so4(mcases,mtempmax)
  real :: wcmr(ndat,mcases),wavelengthr(ndat,mcases)
  real :: wcmi(ndat,mcases),wavelengthi(ndat,mcases)
  real :: rn(ndat,mcases,mtempmax)
  real :: ri(ndat,mcases,mtempmax)
  character(len=30) :: fil
  character (len=16) :: llset

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

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

! ***
! The compound string
!       '1234567890123456'
  llset=' Myhre 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,"rdmyhre_h2so4: id1,filnm ",i8,2x,a100)
  end if

! ***
! Get dimension ids (work with variables part of ncdump)
  istat=nf90_inq_varid(id1,'ncases',idncases)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'ntempmax',idntempmax)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'nlinesr',idnlinesr)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'nlinesi',idnlinesi)
   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)

  if (nopr .eq. 1) then 
   write(iout,fmt=180) idncases,idntempmax,idnlinesr
   180 format(2x,"idncases,idntempmax,idnlinesr ",/,2x,3(1x,i4))
   write(iout,fmt=190) idnlinesi,idntext,idn70
   190 format(2x,"idnlinesi,idntext,idn70 ",/,2x,3(1x,i4))
  end if

! *****
! Get the dimensions
  istat=nf90_get_var(id1,idncases,ncases)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idntempmax,ntempmax)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idnlinesr,nlinesr)
   call handle(iout,istat)
  istat=nf90_get_var(id1,idnlinesi,nlinesi)
   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=200) ncases,ntempmax,nlinesr,nlinesi,ntext,n70
   200 format(/,2x,"ncases,ntempmax,nlinesr,nlinesi,ntext,n70 ",/,&
   2x,6(1x,i4))
  end if

! *****
! Get variable ids
  istat=nf90_inq_varid(id1,'numlinesr',idnumlinesr)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'numlinesi',idnumlinesi)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'temps',idtemps)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'ntemps',idntemps)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'wh2so4',idwh2so4)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'wcmr',idwcmr)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'wavelengthr',idwavelengthr)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'rn',idrn)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'wcmi',idwcmi)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'wavelengthi',idwavelengthi)
   call handle(iout,istat)
  istat=nf90_inq_varid(id1,'ri',idri)
   call handle(iout,istat)

  if (nopr .eq. 1) then
   write(iout,fmt=210) idnumlinesr,idnumlinesi,idtemps,idntemps,idwh2so4
   210 format(/,2x,"idnumlinesr,idnumlinesi,idtemps,idntemps,idwh2so4, ",/,&
   2x,5(1x,i4))
   write(iout,fmt=220) idwcmr,idwavelengthr,idrn
   220 format(/,2x,"idwcmr,idwavelengthr,idrn ",/,2x,3(1x,i4))
   write(iout,fmt=230) idwcmi,idwavelengthi,idri
   230 format(/,2x,"idwcmi,idwavelengthi,idri ",/,2x,3(1x,i4))
  end if

! ****
! Get the data

! The number of lines for real indices
   start(1)=1
   count(1)=ncases
   stride(1)=1
  istat=nf90_get_var(id1,idnumlinesr,numlinesr,start,count,stride)

! The number of lines for imaginary indices
   start(1)=1
   count(1)=ncases
   stride(1)=1
  istat=nf90_get_var(id1,idnumlinesi,numlinesi,start,count,stride)

! ****
! The temperatures
   start2(1)=1
   start2(2)=1
   count2(1)=mcases
   count2(2)=mtempmax
   stride2(1)=1
   stride2(2)=1
  istat=nf90_get_var(id1,idtemps,temps,start2,count2,stride2)

! The number of temperatures for each case
   start(1)=1
   count(1)=mcases
   stride(1)=1
  istat=nf90_get_var(id1,idntemps,ntemps,start,count,stride)

! The wh2so4 values
   start2(1)=1
   start2(2)=1
   count2(1)=mcases
   count2(2)=mtempmax
   stride2(1)=1
   stride2(1)=1
  istat=nf90_get_var(id1,idwh2so4,wh2so4,start2,count2,stride2)

! ****
! The wcmr values
   start2(1)=1
   start2(2)=1
   count2(1)=ndat
   count2(2)=mcases
   stride2(1)=1
   stride2(2)=1
  istat=nf90_get_var(id1,idwcmr,wcmr,start2,count2,stride2)

! The wavelengthr values
   start2(1)=1
   start2(2)=1
   count2(1)=ndat
   count2(2)=mcases
   stride2(1)=1
   stride2(2)=1
  istat=nf90_get_var(id1,idwavelengthr,wavelengthr,start2,count2,stride2)

! Read in the real indices
! Note that count3(1) and count3(2) etc are opposite to the ncdump dump
   start3(1)=1
   start3(2)=1
   start3(3)=1
   count3(1)=ndat
   count3(2)=mcases
   count3(3)=mtempmax
   stride3(1)=1
   stride3(2)=1
   stride3(3)=1
  istat=nf90_get_var(id1,idrn,rn,start3,count3,stride3)
   call handle(iout,istat)

! ****
! The wcmi values
   start2(1)=1
   start2(2)=1
   count2(1)=ndat
   count2(2)=mcases
   stride2(1)=1
   stride2(2)=1
  istat=nf90_get_var(id1,idwcmi,wcmi,start2,count2,stride2)
   call handle(iout,istat)

! The wavelengthi values
   start2(1)=1
   start2(2)=1
   count2(1)=ndat
   count2(2)=mcases
   stride2(1)=1
   stride2(2)=1
  istat=nf90_get_var(id1,idwavelengthi,wavelengthi,start2,count2,stride2)
   call handle(iout,istat)

! Read in the imaginary indices
! Note that count3(1) and count3(2) etc are opposite to the ncdump dump
   start3(1)=1
   start3(2)=1
   start3(3)=1
   count3(1)=ndat
   count3(2)=mcases
   count3(3)=mtempmax
   stride3(1)=1
   stride3(2)=1
   stride3(3)=1
  istat=nf90_get_var(id1,idri,ri,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
  
! Will specify the wghth2so4 and temp pair
! First, will write out the possibilities

! The wh2so4 and temps values
  if (nopr .eq. 1) then
   write(iout,fmt=320) 
   320 format(/,2x,"rdmyhre_h2so4: i,j,wh2so4(i,j),temps(i,j)")
   do i=1,mcases
   do j=1,mtempmax
    write(iout,fmt=330) i,j,wh2so4(i,j),temps(i,j)
    330 format(2x,2(1x,i4),2x,2(1x,f10.4))
   end do
   end do
  end if

! **
! Write out to the screen
   write(0,fmt=360) 
   360 format(/,2x,"rdmyhre_h2so4: i,j,wh2so4(i,j),temps(i,j)")
   do i=1,mcases
   do j=1,mtempmax
    vec1(j)=wh2so4(i,j)
    vec2(j)=temps(i,j)
   end do
    write(0,fmt=365) i,(vec1(j),j=1,5)
    365 format(2x,1(1x,i3),5(1x,f10.4))
    write(0,fmt=367) (vec2(j),j=1,5)
    367 format(2x,4x,5(1x,f10.4))
   end do

! Specify the temp and wghth2so4 values
  write(0,fmt=160)
  160 format(2x,"Specify wghth2so4 value (e.g. 72.0) ")
  read(5,*) wghth2so4

  write(0,fmt=165)
  165 format(2x,"Specify temp value (e.g. 233.0) ")
  read(5,*) temp

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

    diffmin=1.0e3
    ival=-99
    jval=-99
   do i=1,mcases
   do j=1,mtempmax
    diff1=wh2so4(i,j)-wghth2so4
    diff2=temps(i,j)-temp
    a1=sqrt((diff1*diff1)+(diff2*diff2))
    if (a1 .lt. diffmin) then 
     ival=i
     jval=j
     diffmin=a1
    end if
   end do
   end do

  if ((ival .lt. 0) .or. (jval .lt. 0)) then
   write(iout,fmt=350) ival,jval
   350 format(/,2x,"rdmyhre_h2so4: ival,jval ",2(1x,i4)," will stop")
   stop
  end if

  if (nopr .eq. 1) then
   write(iout,fmt=410) ival,jval,wh2so4(ival,jval),temps(ival,jval)
   410 format(/,2x,"rdmyhre_h2so4:ival,jval,wh2so4,temps, ",/,&
   2x,2(1x,i4),2(1x,f10.4))
  end if

! *****
! Put values into the output arrays
! rnval(nlinesr,ncases,ntempmax)
! rival(nlinesr,ncases,ntempmax)
   nlines=nlinesr
  do i=1,nlines
   wcm(i)=wcmr(i,ival)
   wavelength(i)=wavelengthr(i,ival)
   a1=rn(i,ival,jval)
   a2=ri(i,ival,jval)
   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
    800 format(/,2x,"rdmyhre_h2so4: temp,wghth2so4 ",/,2x,2(1x,f10.4),/,&
    2x,"rdmyhre_h2so4: 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
