! **********************************************************************
  program main
 
! *****
! Will use f90 netCDF
  use netcdf

! Has max dimensions of the output indices and extinction spectra vectors
  include 'param1.inc'

  integer :: nopr,iout
  character(len=5) :: fout

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

! See rdwork.f90
  integer :: noprwork,noprtest,noprlist,noprind,noprwave,noprsize,noprext
  integer :: itest,iset,iwrind,iwrsize,iwrext
  integer :: iwave
  real :: w1,w2,dw
  real :: r1,r2
  real :: den1,rad1,sig1
  real :: den2,rad2,sig2
  integer :: icompare,iset1,iset2
  integer :: imix,imix1,imix2

! See rdlistrn.f90
  integer :: nsets
  integer,parameter :: nfmax=100
  character(len=100) :: listrn(nfmax)

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

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

! See calcsized.f90
  integer,parameter :: ndist2=200
  integer :: ndist
  real :: radr(ndist2),sized(ndist2)

! See calcext.f90
  real :: bext(nwavemax),babs(nwavemax)
  real :: bsca(nwavemax),asym(nwavemax)
  real :: omega(nwavemax),back(nwavemax)

! *****
! Output ascii file f.out
   iout=6
  fout='f.out'
  open(iout,form='formatted',file='f.out',status='unknown')
 
  write(iout,fmt=10)
  10 format(/,2x,"main: Hitran_Refractive_Indices_program")

! *************************
! Read in the specifications of the input refractive index and output
! subdirectory to which output text and ps graphs are written to
! Information is supplied by the user in the    directory.dat   file
   nopr=1
  call rdinoutdir(nopr,iout,indir,outdir,fout)
  if (nopr .eq. 1) then
!  stop
  end if

! *************************
! Information is supplied by the user in the    work.dat   file
   noprwork=1
  call rdwork(noprwork,&
  noprtest,noprlist,noprind,&
  noprwave,noprsize,noprext,&
  iout,indir,outdir,fout,&
  itest,iset,&
  iwrind,iwrsize,iwrext,&
  iwave,w1,w2,dw,r1,r2,&
  den1,rad1,sig1,&
  den2,rad2,sig2,&
  icompare,iset1,iset2,&
  imix,imix1,imix2)

! *************************
! Read in the listing of refractive index sets
! from the   indices.dat  ascii listing
  call rdlistrn(noprlist,iout,indir,nsets,listrn)

! *************************
! Will test the routine to make sure you get the desired result
! if (itest eq 1) then begin
!  test1,noprtest,$
!  iout,indir,outdir,fout
! endif

! *************************
! Compare two sets of indices
! if (icompare eq 1) then begin
!  compareindx,nopr,iout,$
!  icompare,iset1,iset2,$
!  nsets,listrn,lstrn,$
!  indir,outdir,iwave,w1,w2
! endif

! *************************
! Compare mixing rules
! if (imix eq 1) then begin
!  comparemix,nopr,iout,$
!  imix,imix1,imix2,$
!  nsets,listrn,lstrn,$
!  indir,outdir,iwave,w1,w2
! endif

! *************************
  write(0,fmt=100)
  100 format(2x," Will call rdindices.f90")

! Read in the indices of refraction for set iset
! Will interpolate to the output wave spectra grid with nwave points
! For igraphi will graph the real and imaginary indices of refraction
  call rdindices(noprind,iout,&
  indir,iset,nsets,listrn,lstrn,&
  nlines,wcm,wavelength,rndat,ridat,&
  outdir,iwave,iwrind,lset,llset,ibothpos,ndifferent)

  write(0,fmt=110)
  110 format(2x," Called rdindices.f90")

! *************************
  write(0,fmt=200)
  200 format(2x," Will call calcwave.f90")

! Obtain the wavenumber (or wavelength) scale
! Will use interpolation to obtain rn and ri indices
  call calcwave(noprwave,iout,&
  nlines,wcm,wavelength,rndat,ridat,&
  iwave,w1,w2,dw,&
  nwave,wave,rn,ri,&
  outdir,iwrind,lset,lstrn,&
  ibothpos)

  write(0,fmt=210)
  210 format(2x," Called calcwave.f90")

! *************************
! Calculate the size distribution with ndist points
  call calcsized(noprsize,iout,&
  r1,r2,&
  den1,rad1,sig1,&
  den2,rad2,sig2,&
  ndist,radr,sized,&
  outdir,iwrsize,lset)

! *************************
! Calculate the spectra (extinction, absorption, scattering, single scattering
! albedo, and asymmetry parameter)
  call calcext(nopr,iout,&
  iwave,nwave,wave,rn,ri,&
  ndist,radr,sized,&
  bext,babs,bsca,asym,back,omega,&
  outdir,iwrext,lset)

! *************************
  stop
  end
