module ave_stelspec_mod
      
      implicit none
      
      contains
      
      subroutine ave_stelspec(STELLAR)

!==================================================================
!     
!     Purpose
!     -------
!     Average the chosen high resolution stellar spectrum over the
!     visible bands in the model.
!     
!     Authors
!     ------- 
!     Robin Wordsworth (2010).
!     Generalized to very late spectral types (and Brown dwarfs) Jeremy Leconte (2012)
!     Modified to account for any stellar spectrum file (Lucas Teinturier and Martin Turbet, 2023-2025)
!
!     Called by
!     ---------
!     setspv.F
!     
!     Calls
!     -----
!     none
!     
!==================================================================

      use radinc_h, only: L_NSPECTV
      use radcommon_h, only: BWNV, DWNV, tstellar
      use datafile_mod, only: datadir
      use callkeys_mod, only: stelbbody,stelTbb
      use ioipsl_getin_p_mod, only: getin_p

      implicit none

      real*8 STELLAR(L_NSPECTV)

      integer Nfine
      integer,parameter :: Nfineband=200
      integer ifine,band

      real,allocatable,save :: lam(:),stel_f(:) ! read by master thread
                                                ! but used by all threads
      real lamm,lamp
      real dl

      character(len=100) :: file_id,file_id_lam
      character(len=200) :: file_path,file_path_lam
      character(len=150) :: stelspec_file

      real lam_temp
      double precision stel_temp
      
      integer :: ios ! file opening/reading status
      logical :: file_exists

      STELLAR(:)=0.0

      print*,'enter ave_stellspec'
      if(stelbbody)then
         tstellar=stelTbb
	 Nfine=L_NSPECTV*Nfineband
	 do band=1,L_NSPECTV
	    lamm=10000.0/BWNV(band+1)
	    lamp=10000.0/BWNV(band)
	    dl=(lamp-lamm)/(Nfineband)
	    do ifine=1,Nfineband
	       lam_temp=lamm+(lamp-lamm)*(ifine-1.)/(Nfineband)
               call blackl(dble(lam_temp*1e-6),dble(tstellar),stel_temp)
	       STELLAR(band)=STELLAR(band)+stel_temp*dl
	    enddo	    
	 end do
         STELLAR(1:L_NSPECTV)=STELLAR(1:L_NSPECTV)/sum(STELLAR(1:L_NSPECTV))
      else !stelbbody
         ! look for a " tstellar= ..." option in def files
         tstellar = -1. ! default
         call getin_p("tstellar",tstellar) ! default path
         if (tstellar.eq.-1.) then
	   write(*,*)'Beware that startype is now deprecated, you should use '
	   write(*,*)'stelspec_file and tstellar to define the input stellar spectrum.'
	   write(*,*)'     '
           write(*,*)'Error: tstellar (effective stellar temperature) needs to be specified'
           write(*,*)'in callphys.def: tstellar=...'
           call abort_physic("ave_stelspec", "tstellar needs to be specified",1)
         end if
	 
         write(*,*) "Input stellar temperature is:"
         write(*,*) "tstellar = ",tstellar

         ! load high resolution stellar data
         ! look for a " stelspec_file= ..." option in def files
         stelspec_file = "None" ! default
         call getin_p("stelspec_file",stelspec_file) ! default path
	 
         write(*,*) "Input stellar spectrum file is:"
         write(*,*) "stelspec_file = ",trim(stelspec_file)
         write(*,*) 'Please use ',1,' and only ',1,' header line in ',trim(stelspec_file)

         ! Check the target file is there
         file_path = trim(datadir)//'/stellar_spectra/'//stelspec_file
         print*, 'stellar flux : ', file_path
         inquire(FILE=file_path,EXIST=file_exists)         
   
         if (.not.file_exists) THEN
	   write(*,*)'Beware that startype is now deprecated, you should use '
	   write(*,*)'stelspec_file and tstellar to define the input stellar spectrum.'
	   write(*,*)'     '
           write(*,*)'Error: cannot open stelspec_file file ', trim(stelspec_file)
           write(*,*)'It should be in :',trim(datadir),'/stellar_spectra/'
           write(*,*)'1) You can change the data directory in callphys.def'
           write(*,*)'   with:'
           write(*,*)'   datadir=/path/to/the/directory'
           write(*,*)'2) You can change the input stelspec_file file name in'
           write(*,*)'   callphys.def with:'
           write(*,*)'   stelspec_file=filename'
	   write(*,*)'You can check the online repository to search for '
	   write(*,*)'available stellar spectra here : '
	   write(*,*)'https://web.lmd.jussieu.fr/~lmdz/planets/generic/datagcm/stellar_spectra/'
           call abort_physic("ave_stelspec", "Unable to read stellar flux file", 1)
         end if

!$OMP MASTER
         ! Open the file
         OPEN(UNIT=110,FILE=file_path,STATUS='old',iostat=ios)
         ! Get number of line in the file
         READ(110,*) ! skip first line header just in case
         Nfine = 0
         do
           read(110,*,iostat=ios)
           if (ios<0) exit
           Nfine = Nfine + 1
         end do
         rewind(110) ! Rewind file after counting lines
         READ(110,*) ! skip first line header just in case

	 allocate(lam(Nfine),stel_f(Nfine))

         do ifine=1,Nfine
           read(110,*) lam(ifine), stel_f(ifine) ! lam [um] stel_f [per unit of wavelength] (integrated and normalized by Fat1AU)
         enddo

!$OMP END MASTER
!$OMP BARRIER
	 
         ! sum data by band
         band=1
	 Do while(lam(1).lt. real(10000.0/BWNV(band+1)))
	    if (band.gt.L_NSPECTV-1) exit
            band=band+1
	 enddo
	 dl=lam(2)-lam(1)
         STELLAR(band)=STELLAR(band)+stel_f(1)*dl
         do ifine = 2,Nfine
            if(lam(ifine) .gt. real(10000.0/BWNV(band)))then
               band=band-1
            endif
            if(band .lt. 1) exit
	    dl=lam(ifine)-lam(ifine-1)
            STELLAR(band)=STELLAR(band)+stel_f(ifine)*dl
         end do
	       
	 
         STELLAR(1:L_NSPECTV)=STELLAR(1:L_NSPECTV)/sum(STELLAR(1:L_NSPECTV))
!$OMP BARRIER
!$OMP MASTER
	 deallocate(lam)
	 deallocate(stel_f)
!$OMP END MASTER
!$OMP BARRIER         
      endif !stelbbody

      end subroutine ave_stelspec
      
      end module ave_stelspec_mod