module params

  real,parameter:: Sunl=3.826e33, Sunl30=3.826e3, Bolsun=4.77
  real,parameter:: sigma=5.669e-5,pigreco=3.141593

  integer, parameter:: maxlib=50
  integer, save:: nages=150,ntab,ntab1,nlam
  integer, save, allocatable:: nmax(:) 
  real, save, allocatable::fluxes(:,:,:),ztab(:),ages(:,:),wave(:)
  
end module params

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

Program READFLUX
  
  use params  
  implicit none
  
  real,allocatable::  HNU(:),norma(:)
  integer,allocatable:: h_int(:),w_int(:)
  Character(len=80):: TITLE,titleb,fileasc(maxlib),fileasca,nomfil
  Character(len=200):: filein
  real:: teff,agelog,fflum,ssp_mass,ssp_zeta,freq,ssp_mass_orig
  integer:: jz,jjz,jok,il,ia,err_alloc,nlamold,err
  real, parameter::safe=1.0

!!!!!!!!!!!!!!!!!!!!
!!Example of file sspfile.inp : comments in the first 2 raws, then
!!in each raw index (1 or 0) and name:

!Names of SSPs to be used in increasing Z order.List ended by a comment line
!First column put 1 if ssp is to be used, 0 otherwise. 
!1    'Z0004.sal'
!1    'Z004.sal'
!1    'Z008.sal'
!1    'Z02.sal'
!1    'Z05.sal'
!# 
!!!!!!!!!!!!!!!!!!!!!!
  
!open file sspfile.inp from where the ssp files names are read
  nomfil='sspfile.inp'
  open (15,file=nomfil)
!skip first 2 raws
  read(15,*)
  read(15,*)
    
  jz=0
  do jjz=1,maxlib
  !read index and file name in each raw. Index=0 means ssp not used, otherwise index=1
     read(15,end=901,err=901,fmt=*)jok,fileasca  
     if(jok.eq.1) then
        jz=jz+1
        fileasc(jz)=fileasca
     endif
  enddo
901 close(15)
  ntab=jz  !number of metallicities

  allocate(norma(ntab),ztab(ntab), STAT=ERR_ALLOC)
  IF (ERR_ALLOC /= 0) THEN
       print *,"ALLOCATION ERROR in readflux"
  ENDIF
  norma=0;ztab=0
  norma=1. !not used 
  

! read SSP header
  do jz=1,ntab

      filein=fileasc(jz)
      open(15,file=filein,status='old',iostat=err)     
      if (err/=0) then
            print *,'ERROR in opening SSP file ',filein
            stop
      endif     
      read(15,'(a)') titleb ! information on IMF: mass limits and slope
      read(15,*) ssp_mass_orig ! ssp_mass written in ssp (IMF integral)

      read(15,*) ztab(jz)  ! Metallicity of ssp
      read(15,*) nlam      ! number of wavelenghts for SSP SED
     
         if (jz==1) nlamold=nlam
         if (nlamold/=nlam) then
             print *,'SSP files have different wavelenghts'
         endif
         nlamold=nlam
         close(15)
 
  enddo  !metallicity cicle
 
  allocate(wave(nlam),hnu(nlam),fluxes(nlam,nages,ntab),&
  w_int(nlam),h_int(nlam),ages(nages,ntab),nmax(ntab), STAT=ERR_ALLOC)
  IF (ERR_ALLOC /= 0) THEN
     print *, "ALLOCATION ERROR in readflux"
  ENDIF
  wave=0;hnu=0;fluxes=0;w_int=0;h_int=0;ages=0;nmax=0

! now read SSP SED
! ia=age index (1=oldest age, nmax youngest)

  do jz=1,ntab

      filein=fileasc(jz)
      open(15,file=filein,status='old',iostat=err)
      !now the header can be jumped
      read(15,'(a)') titleb ! information on IMF
      read(15,*) ; read(15,*) ; read(15,*) 
     
      read(15,'(20I6)') W_INT
      wave(:)=EXP(w_int(:)/10000.0) 
      ! transform from nanometer to angstrom 
      wave(:)=wave(:)*10.

      do ia=1,1000

          if(ia .gt. nages) then
               write(*,*) ' recompile with greater NAGES'
               stop
          endif
          read(15,'(a)',END=90)TITLE
          read(TITLE,*)TEFF,agelog,fflum
        
          ages(ia,jz)=agelog
        
          ! read the flux in 10**30 ergs/s/hz/Msun
          READ(15,'(20i5)') H_int
          hnu(:)=h_int(:)/1000.0-100
          where (hnu>(-range(hnu(1))+safe)*log(10.0)) 
               hnu=exp(hnu)
          elsewhere
               hnu=0
          end where
         ! SSP spectrum in 10^30  ergs/sec/Angstrom/Msun   
         fluxes(:,ia,jz)=4*norma(jz)*hnu(:)*2.99792458e18/wave(:)/wave(:)

         enddo

    90   continue
         nmax(jz)=ia-1
         close(15)

      enddo
  
      ages(:,:)=(10**ages(:,:))/1.e9  !SSP age in Gyr

      deallocate(w_int,h_int,hnu,norma)


!!!!!!!!Write all the input SSP files in new files, with wavelengths and
!luminosities in columns. 
!Change here if you want different formats 

do jz=1,ntab

      nomfil=trim(fileasc(jz))//'.col'
      open(1,file=nomfil)
      write(1,'("#1st column: wave[A], then Lum[10^30 erg/s/A/Msun]")')
      !2nd raw: write SSP ages in Gyr
      write(1,'("#Age[Gyr]=",1x,100(1pe10.2))') ages(1:nmax(jz),jz) 
      write(1,'("#")')     
      do il=1,nlam
        !wave in Ang, fluxes in 10^30 erg/s/A/Msun
        write(1,'(1x,100(1pe10.2))') wave(il), fluxes(il,1:nmax(jz),jz) 
      enddo
      close(1)

enddo 



end Program READFLUX