
 subroutine get_solar_spectrum( &
    nsmax,                 ns_elas,     nwave,      &
    dnu,                   nwave_aux,   wavenr_aux, &
    solar_irradiance_aux)

!-----------------------------------------------------------
! Subroutine:
!   get_solar_spectrum
!
! Purpose:
!   Reads MODTRAN solar spectrum given on 1cm-1 grid and
!   returns this reference solar spectrum on the desired
!   wavenr grid in [photons cm-2 (cm-1)-1 s-1]
!
! Output:
!    solar_irradiance_aux 
!    on wavenr_aux grid
!-------------------------------------------------------------------

  implicit none

  integer,parameter::dp=selected_real_kind(8) !double precision

  integer,intent(in):: &
    nsmax, &
    ns_elas, &
    nwave, &  !wavenr grid
    nwave_aux !extended wavenr grid

  integer:: & 
    istat,          &
    l,              &
    n,              &
    nstart,         & !1cm-1 index corresponding to wavenr_aux(1) 
    nend,           & !1cm-1 " to wavenr_aux(nwave_aux)
    nwave_1cm,      &
    sun_data_switch
  
  real(kind=dp):: &
    nu_min,     &
    nu_max,     &
    wavenr_1cm_min, &
    wavenr_1cm_max, &
    a, &
    b, &
    filter_radius, &
    scale_factor

  real(kind=dp),dimension(:),allocatable:: &
    wavenr_1cm, &
    solar_irradiance_1cm

  integer,parameter:: &
    ndata = 49951

  real(kind=dp),intent(in):: &
    dnu,        & !wavenr binwidth [cm-1]
    wavenr_aux(nwave_aux)

  real(kind=dp),intent(out):: &
    solar_irradiance_aux(nwave_aux)

  real(kind=dp)::nu0,dnu0

!---------------------------------------------------------
! Read extraterrestrial solar irradiance [W cm-2 cm]
! on 1cm-1 wavenr grid 
!
! sun_data_switch = 1 chkur     Chance & Spurr
!                 = 2 newkur    Kurucz
!                 = 3 cebchkur  
!                 = 4 thkur
!                 = 5 vandeelen2005, Eq.(20), for testing
!                 = 6 constant, for testing
!---------------------------------------------------------

  sun_data_switch = 1 !chkur
! sun_data_switch = 6 

! wavenr range of input 
  nu_min =    50.d0 !cm-1 = wavelength 200000 nm
  nu_max = 50000.d0 !cm-1 = wavelength 200 nm 

! filter radius of resample function
  filter_radius = 2.d0

! wavenr range of 1 cm-1 spectrum
  wavenr_1cm_min = wavenr_aux(1)         - filter_radius*dnu
  wavenr_1cm_max = wavenr_aux(nwave_aux) + filter_radius*dnu

  if(wavenr_1cm_min.lt.nu_min) & 
    write(*,*) 'ietsun.f90: Uhoh... outside 1cm-1 wavenr range'
  if(wavenr_1cm_max.gt.nu_max) &
    write(*,*) 'getsun.f90: Uhoh... outside 1cm-1 wavenr range'

! read solar irradiance spectrum from file 
! on 1cm-1 wavenr grid

  nwave_1cm = int(wavenr_1cm_max-wavenr_1cm_min)+2
  
  allocate( &
    wavenr_1cm(nwave_1cm),          &
    solar_irradiance_1cm(nwave_1cm), &
    stat=istat)

!-------------------------------------------------------------
! Read one of MODTRAN solar spectra
!-------------------------------------------------------------

  if(sun_data_switch.eq.1) &
     open(1,file='../data/chkur.dat',status='old')
  if(sun_data_switch.eq.2) &
     open(1,file='../data/newkur.dat',status='old')
  if(sun_data_switch.eq.3) &
     open(1,file='../data/cebchkur.dat',status='old')
  if(sun_data_switch.eq.4) &
     open(1,file='../data/thkur.dat',status='old')
  if((sun_data_switch.eq.5).or.(sun_data_switch.eq.6))then
     open(1,file='../data/thkur.dat',status='old') !dummy
  endif

  read(1,*)
  read(1,*)

  n=1
  do l=1,ndata
     read(1,*) a,b
     if(a.ge.wavenr_1cm_min)then
         wavenr_1cm(n)           = a
         solar_irradiance_1cm(n) = b
         if(abs(wavenr_1cm(n)-wavenr_aux(1)).lt.1.d-4) nstart=n
         if(abs(wavenr_1cm(n)-wavenr_aux(nwave_aux)).lt.1.d-4) nend=n
         n=n+1
         if(n.gt.nwave_1cm) goto 11
     endif
   enddo
11 continue

  close(1)
 
!--------------------------------------------------------------
! Artificial solar spectrum for testing purposes 
! with one solar Fraunhofer line
!-------------------------------------------------------------

  if(sun_data_switch.eq.5)then  

     nu0  = 25414.d0 !cm-1
     dnu0 = 70.d0 !cm-1
 
     do n=1,nwave_1cm
        if(abs(wavenr_1cm(n)-nu0).lt.dnu0)then
            solar_irradiance_1cm(n) = 1000.d0*(1.d0 - &
               0.82d0*(1.d0-abs(wavenr_1cm(n)-nu0)/dnu0))
        else 
            solar_irradiance_1cm(n) = 1000.d0 * 1.d0
        endif
     enddo

!    open(10,file='../output/spectrum_1cm.dat',status='unknown')
!    write(10,'(a21)') '# spectrum 1cm-1 grid'
!    write(10,'(a8,a12)') '# w [nm]','F [ ]'
!    do n=nwave_1cm,1,-1
!       write(10,'(f8.4,f12.4)') &
!        1.d7/wavenr_1cm(n),solar_irradiance_1cm(n)
!    enddo
!    close(10)

  endif

!--------------------------------------------------------------
! Artificial solar continuum for testing purposes 
!-------------------------------------------------------------

  if(sun_data_switch.eq.6)then  
 
     do n=1,nwave_1cm
        solar_irradiance_1cm(n) = 1000.d0
     enddo

  endif

!----------------------------------------------------------
! Resample solar spectrum on 1cm-1 grid to solar spectrum 
! on wavenr_aux grid  
!----------------------------------------------------------

! 1 cm-1 binwidth / dnu binwidth [cm-1] 
  scale_factor = 1.d0/dnu

  call resample( &
       nsmax,          ns_elas,      dnu,                  &
       filter_radius,  nstart,       nend,                 &
       nwave_1cm,      wavenr_1cm,   solar_irradiance_1cm, &
       nwave_aux,      wavenr_aux,   solar_irradiance_aux)

!----------------------------------------------------------
! Test
!----------------------------------------------------------
! do n=1,nwave_aux
!    solar_irradiance_aux(n) = 1.d0
!    if(n.le.nsmax)         solar_irradiance_aux(n) = 0.d0
!    if(n.ge.nsmax+nwave+1) solar_irradiance_aux(n) = 0.d0
! enddo
!----------------------------------------------------------

  deallocate( &
    wavenr_1cm,            &
    solar_irradiance_1cm)

  return
  end subroutine get_solar_spectrum
