  subroutine calcwave(nopr,iout,&
  nlines,wcm,wavelength,rndat,ridat,&
  iwave,w1,w2,dw,&
  nwave,wave,rn,ri,&
  outdir,iwrind,lset,lstrn,&
  ibothpos)

! *************************
  integer :: nopr,iout

! See rdinoutdir.f90
  character(len=60) :: indir,outdir

! See rdwork.f90
  integer :: noprwork,noprtest,noprlist,noprind,noprwave,noprsize,noprext
  integer :: iwrind
  integer :: iwave
  real :: w1,w2,dw

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

! See calcwave.f90
! wave is wavelength or wcm
  integer :: nwave
  real :: wave(nwavemax)
  real :: rn(nwavemax),ri(nwavemax)

! used here
  integer :: mn,iorder
  integer :: j1,j2,jskip,ioops,iok,i1,i2
  real :: dwd,deriv1,deriv2,a1,a2

! *************************
! Check on input
   iwrsp=1
  if (iwrsp .eq. 1) then
   write(iout,fmt=10) nlines
   10 format(/,2x,"calcwave: nlines ",i4,/,&
   2x,"i,wcm(i),wavelength(i),rndat(i),ridat(i)")
   do i=1,nlines
    write(iout,fmt=20) i,wcm(i),wavelength(i),rndat(i),ridat(i)
    20 format(2x,i4,2(1x,f10.4),1p,2(1x,e10.3))
   end do
!  stop
  end if

! *************************
! For case with only 1 index of refraction
  if (nlines .eq. 1) then
   write(iout,fmt=100)
   100 format(/,2x,"calcwave: there is only one wavelength")
     nwave=nlines
    rn(1)=rndat(0)
    ri(1)=ridat(0)
   if (iwave .eq. 2) then
    wave(1)=wavelength(1)
   end if
   if (iwave .eq. 1) then
    wave(1)=wcm(1)
   end if
   go to 1000
  end if

! *************************
! For wavenumber scale
  if (iwave .eq. 1) then

! Check that wavenumbers are OK
    if ((w1 .lt. 100.0) .or. (w2 .lt. 100.0)) then
     write(iout,fmt=200)
     200 format(/,2x,"calcwave: iwave=1, expect w1 and w2 > 100.0 cm-1",&
     /,2x,"calcwave: will stop")
     stop
    end if

! There will be nwave spectral points in the extinction spectrum
     a1=(w2-w1)/dw
     nwave=int(a1)

! Wave is in wavenumbers
    do i=1,nwave 
     wave(i)=w1+((i-1)*dw)
    end do

! wavenumber increases from wcm(0) to wcm(1)
    if (wcm(2) .ge. wcm(1)) then
     iorder=1
     j1=1
     j2=nlines
     jskip=1
      ioops=0
     if (wcm(1) .gt. wave(1)) then 
      ioops=1
     end if
     if (wave(nwave) .gt. wcm(nlines)) then
      ioops=1
     end if
    end if

! wavenumber decreases from wcm(0) to wcm(1)
    if (wcm(1) .ge. wcm(2)) then
     iorder=2
     j2=1
     j1=nlines
     jskip=-1
      ioops=0
     if (wcm(nlines) .gt. wave(1)) then 
      ioops=1
     end if
     if (wave(nwave) .gt. wcm(1)) then 
      ioops=1
     end if
    end if

    if (ioops .eq. 1) then 
     write(iout,fmt=210) iwave,wave(1),wave(nwave),wcm(1),wcm(nlines),&
     wavelength(1),wavelength(nlines)
     210 format(/,2x,&
     2x,"calcwave: refractive index set does not have data ",/,&
     2x,"calcwave: for the full wave grid ",/,&
     2x,"calcwave: iwave ",i4,/,&
     2x,"calcwave: wave(1),wave(nwave) ",/,2x,2(1x,f10.4),/,&
     2x,"calcwave: wcm(1),wcm(nlines) ",/,2x,2(1x,f10.4),/,&
     2x,"calcwave: wavelength(1),wavelength(nlines) ",/,2x,2(1x,f10.4))
     stop
    end if

! **
! Loop over the output grid wave with nwave points
   do i=1,nwave 

! **
! wavenumber increases as j increases 
    if (iorder .eq. 1) then

      iok=0
     do j=j1,j2,jskip 
      if (wcm(j) .gt. wave(i)) then 
       iok=1
       i1=j-1
       i2=j
       j1=i1
       dw2=wave(i)-wcm(i1)
       go to 220
      end if
     end do
     220 mn=0

     if (iok .eq. 1) then
      dwd=wcm(i2)-wcm(i1)
      deriv1=(rndat(i2)-rndat(i1))/dwd
      deriv2=(ridat(i2)-ridat(i1))/dwd
       a1=rndat(i1)+(deriv1*dw2)
      if (a1 .lt. 0.0) then 
       a1=0.0
      end if
       a2=ridat(i1)+(deriv2*dw2)
      if (a2 .lt. 0.0) then 
       a2=0.0
      end if
      rn(i)=a1
      ri(i)=a2
     end if
     
    end if
! iorder=1

! **
! wavenumber decreases
    if (iorder .eq. 2) then

      iok=0
     do j=j1,j2,jskip 
      if (wcm(j) .gt. wave(i)) then
       iok=1
       i1=j
       i2=j+1
       j1=i2
       dw2=wave(i)-wcm(i2)
       go to 230
      end if
     end do
     230 mn=0

     if (iok .eq. 1) then 
      dwd=wcm(i2)-wcm(i1)
      deriv1=(rndat(i2)-rndat(i1))/dwd
      deriv2=(ridat(i2)-ridat(i1))/dwd
       a1=rndat(i2)+(deriv1*dw2)
      if (a1 .lt. 0.0) then 
       a1=0.0
      end if
       a2=ridat(i2)+(deriv2*dw2)
      if (a2 .lt. 0.0) then 
       a2=0.0
      end if
      rn(i)=a1
      ri(i)=a2
     end if
     
    end if
! iorder=1

   end do
! loop over nwave wave

  end if
! iwave=1

! *************************
! For wavelength scale
  if (iwave .eq. 2) then

! Check that wavelengths are OK
    if ((w1 .ge. 100.0) .or. (w2 .ge. 100.0)) then
     write(iout,fmt=300)
     300 format(/,2x,"calcwave: iwave=2, expect w1 and w2 < 100 microns",/,&
     2x,"Will stop in calcwave")
     stop
    end if

! There will be nwave spectral points in the extinction spectrum
     a1=(w2-w1)/dw
     nwave=int(a1)

! Wave is in microns
    do i=1,nwave
     wave(i)=w1+((i-1)*dw)
    end do

! Check end points

! wavelength increases from wavelength(1) to wavelength(2)
    if (wavelength(2) .ge. wavelength(1)) then
     iorder=1
     j1=1
     j2=nlines
     jskip=1
      ioops=0
     if (wave(1) .lt. wavelength(1)) then
      ioops=1
     end if
     if (wave(nwave) .gt. wavelength(nlines)) then
      ioops=1
     end if
    end if

! wavelength decreases from wavelength(1) to wavelength(2)
    if (wavelength(1) .ge. wavelength(2)) then
     iorder=2
     j2=1
     j1=nlines
     jskip=-1
      ioops=0
     if (wave(1) .lt. wavelength(nlines)) then
      ioops=1
     end if
     if (wave(nwave) .gt. wavelength(1)) then
      ioops=1
     end if
    end if

    if (ioops .eq. 1) then 
     write(iout,fmt=310) iwave,wave(1),wave(nwave),wcm(1),wcm(nlines),& 
     wavelength(1),wavelength(nlines)
     310 format(/,2x,&
     2x,"calcwave: refractive index set does not have data ",/,&
     2x,"calcwave: for the specified wave grid, will stop ",/,&
     2x,"calcwave: iwave ",i4,/,&
     2x,"calcwave: wave(1),wave(nwave) ",/,2x,2(1x,f10.4),/,&
     2x,"calcwave: wcm(1),wcm(nlines) ",/,2x,2(1x,f10.4),/,&
     2x,"calcwave: wavelength(1),wavelength(nlines) ",/,2x,2(1x,f10.4))
     stop
    end if

! **
! Loop over the output grid wave with nwave points
   do i=1,nwave 

! **
! wavelength increases as j increases
    if (iorder .eq. 1) then 

      iok=0
     do j=j1,j2,jskip 
      if (wavelength(j) .gt. wave(i)) then
       iok=1
       i1=j-1
       i2=j
       if (i1 .lt. 0) then
        stop
       end if
       j1=i1
       dw2=wave(i)-wavelength(i1)
       go to 320
      end if
     end do
     320 mn=0

     if (iok .eq. 1) then 
      dwd=wavelength(i2)-wavelength(i1)
      deriv1=(rndat(i2)-rndat(i1))/dwd
      deriv2=(ridat(i2)-ridat(i1))/dwd
       a1=rndat(i1)+(deriv1*dw2)
      if (a1 .lt. 0.0) then 
       a1=0.0
      end if
       a2=ridat(i1)+(deriv2*dw2)
      if (a2 .lt. 0.0) then 
       a2=0.0
      end if
      rn(i)=a1
      ri(i)=a2
     end if
     
    end if
! iorder=1

! **
! wavelength decreases 
    if (iorder .eq. 2) then

      iok=0
     do j=j1,j2,jskip 
      if (wavelength(j) .gt. wave(i)) then 
       iok=1
       i1=j
       i2=j+1
       j1=i2
       dw2=wave(i)-wavelength(i2)
       go to 330
      end if
     end do
     330 mn=0

     if (iok .eq. 1) then 
      dwd=wavelength(i2)-wavelength(i1)
      deriv1=(rndat(i2)-rndat(i1))/dwd
      deriv2=(ridat(i2)-ridat(i1))/dwd
       a1=rndat(i2)+(deriv1*dw2)
      if (a1 .lt. 0.0) then 
       a1=0.0
      end if
       a2=ridat(i2)+(deriv2*dw2)
      if (a2 .lt. 0.0) then
       a2=0.0
      end if
      rn(i)=a1
      ri(i)=a2
     end if
     
    end if
! iorder=2
! **

   end do
  end if

! *************************
! For case with one wavelength
  1000 mn=0

! *************************
   if (nopr .eq. 1) then 
    write(iout,fmt=800)
    800 format(/,2x,"calcwave: i,wave(i),rn(i),ri(i)")
    do i=1,nwave 
     write(iout,fmt=810) i,wave(i),rn(i),ri(i)
     810 format(2x,i4,2x,1p,4(1x,e10.3))
    end do
   end if

! *************************
! Write out the interpolated indices to ascii and netCDF files
   if (iwrind .eq. 1) then 
    call wrindices_interp(nopr,iout,outdir,&
    nwave,wave,rn,ri,iwave,&
    lset)
   end if

! *************************
  if (ibothpos .eq. 0) then
   write(iout,fmt=900)
   900 format(/,2x," calcwave: ibothpos=0, real or imaginary indices are zero",&
   2x,"calcwave: will stop")
   stop
  end if

! *************************
  return
  end
