program spectrum_BYTe
 !

 integer isym,nsym,npoints,ipoint,j,jf,ji,ilevel,ilevelf,ileveli,igamma,igammaf,igammai,info,j0,i
 integer ib,ie,ifreq,nfiles,j_t,tau_,v(6)
 !
 integer  nlevels,nn,itemp

 real(8) temp,temp0,freql,freqr,gns(10),halfwidth,meanmass,partfunc,dfreq,dfreq_,acoef,abscoef,&
         energy,energyf,energyi,tranfreq,cmcoef,beta,ln2,pi,dpwcoef,absthresh,linestr,emcoef,beta0,exp_coef
 real(8) :: thresh = 1e-16
 real(8),allocatable :: freq(:),intens(:),energies(:),pf(:,:)
 integer,allocatable :: sym(:),n(:,:),jrot(:),krot(:),taurot(:),l(:,:),symv(:)

 character(4) gamma,gammaf,gammai,symlab(10)
 character(60) intfilename(200),enrfilename,swapfile
 character(30) proftype,specttype
 character(2)  :: ch_t2
 character(1)  :: ch_t1
 character(95) :: ch_t95
 character(16)  :: ch_npform
 !
 logical :: swap = .false.,partfunc_do=.false.


 real(8),parameter :: planck=6.6260693d-27,avogno=6.0221415d+23,vellgt=2.99792458d+10,boltz=1.380658d-16

!read number of intensity files
 read*,nfiles

if (nfiles>200) then 
  print('("Too many files (>200):",i8)'),nfiles
  stop 'Too many files'
endif

!read intensities filenames
do i = 1,nfiles
 read*,intfilename(i)
 !print('(1x,50a)'),intfilename(i)
enddo

!read energies filename
 read*,enrfilename
 !print('(1x,50a)'),enrfilename

!read swapfile filename
 read*,swapfile
 !print('(1x,50a)'),swapfile

!read number of irreps
 read*,nsym

!read nuclear spin statistics weights
 read*,(symlab(isym),gns(isym),isym=1,nsym)

!read temperature, partition function
 read*,temp,partfunc

!read frequency range and number of points
 read*,freql,freqr,npoints

!read type of profiling
 read*,proftype

!read type of spectra
 read*,specttype

!read gaussian half-width or molecular mean mass
 select case (proftype(1:5))
    case ('gauss'); read*,halfwidth
    case ('doppl'); read*,meanmass
    case ('stick'); read*,absthresh
 end select

!read the threshold
 read*,thresh

!compute constants
 ln2=log(2.0d0)
 pi=acos(-1.0d0)
 beta=planck*vellgt/(boltz*temp)
 cmcoef=avogno/(8.0d0*pi*vellgt)
 emcoef=avogno*planck*vellgt/(4.0d0*pi)
 dpwcoef=sqrt(2.0*ln2*boltz*avogno)/vellgt

!set up freq. grid
 allocate(freq(npoints),stat=info); if (info/=0) stop 'error: freq is out of memory'
 allocate(intens(npoints),stat=info); if (info/=0) stop 'error: intens is out of memory'
 dfreq=(freqr-freql)/real(npoints-1,8)
 forall(ipoint=1:npoints) freq(ipoint)=freql+real(ipoint-1,8)*dfreq

!count number of lines (levels) in the Energy files
 open(unit=1,file=trim(enrfilename))
 i = 0
 do
    !
    read(1,"(i8)",end=10) igammai
    !
    i = i + 1
    !
    cycle
    !
10  exit
 end do
 !
 rewind(1)
 !
 nlevels = i
 !
 allocate(sym(nlevels),symv(nlevels),energies(nlevels),Jrot(nlevels),n(nlevels,4),l(nlevels,3:4),stat=info)
 if (info/=0) stop 'error: sym,energies,Jrot,n,l are out of memory'
 allocate(krot(nlevels),taurot(nlevels),stat=info)
 if (info/=0) stop 'error: krot,taurot is out of memory'
 !
 ! In case of specttype = partfunc the partition function will be computed for a series of temperatures. 
 ! npoints stands for the number of temperature steps, and the frequency limits are used to set the maximal temperature
 !
 if (trim(specttype)=='partfunc') then
   !
   write(ch_npform,"(i4,'(f9.2,10x)')") npoints
   ch_npform = adjustl(ch_npform)
   !
   dfreq=temp/real(npoints,8)
   !
   print('("!","   J ",2x,'//ch_npform//')'),(i*dfreq,i=1,npoints)
   allocate(pf(0:2,npoints),stat=info)
   if (info/=0) stop 'error: pf'
   !
   write(ch_npform,"(i4,'(es18.8,1x)')") npoints
   ch_npform = adjustl(ch_npform)
   !
 endif
 !
 ! prepare computing the partition function
 if (partfunc<=0) then
    partfunc=0.0 ; j0=0
    partfunc_do = .true.
    if (trim(specttype)=='partfunc') pf = 0
 endif
 !
 ! start reading the energies 
 !
 do i=1,nlevels
    !
    read(1,*) nn,jrot(i),sym(i),ilevel,energies(i),n(i,1:4),l(i,3:4),tau_,j_t,krot(i),taurot(i),v(1:6),symv(i)
    !
    j = jrot(i)
    igamma = sym(i)
    energy = energies(i)
    !
    if (partfunc_do) then
      if (j/=j0) then 
         !
         if (trim(specttype)=='partfunc') then 
           print('("!",i4,2x,'//ch_npform//')'),j0,pf(0,1:npoints)
         else 
           print('("#",i4,1x,es16.8)'),j0,partfunc
         endif
         !
      endif
      !
      ! symmetry number
      ! partition function
      !
      partfunc=partfunc+real((2*j+1),8)*gns(igamma)*exp(-beta*energy)
      !
      j0 = j
     endif
     !
     if (trim(specttype)=='partfunc') then
       !
       do itemp = 1,npoints
         !
         if (energy>freqr) cycle
         !
         temp0 = real(itemp,8)*dfreq
         !
         beta0 = planck*vellgt/(boltz*temp0)
         !
         ! apart from the partition function, the first two derivatieves are also computed:
         !
         pf(0,itemp) = pf(0,itemp) + real((2*j+1),8)*                  gns(igamma)*exp(-beta0*energy)
         pf(1,itemp) = pf(1,itemp) + real((2*j+1),8)*(beta0*energy)   *gns(igamma)*exp(-beta0*energy)
         pf(2,itemp) = pf(2,itemp) + real((2*j+1),8)*(beta0*energy)**2*gns(igamma)*exp(-beta0*energy)
         !
       enddo
       !
     endif
     !
 end do
 close(1)
 !
 ! print out the computed part. function objects and finish
 !
 if (trim(specttype)=='partfunc') then
   !
   print('("!0",i4,1x,'//ch_npform//')'),j0,pf(0,1:npoints)
   print('("!1",i4,1x,'//ch_npform//')'),j0,pf(1,1:npoints)
   print('("!2",i4,1x,'//ch_npform//')'),j0,pf(2,1:npoints)
   !
   print('(1x,a,1x,es16.8)'),'! partition function value is',partfunc
   !
   stop
   !
 else
   !
   print('("#",i4,1x,es16.8)'),j0,partfunc
   !
 endif
 !
 print('(1x,a,1x,es16.8)'),'# partition function value is',partfunc
 !
 intens=0.0
 !
 ! the selected by the frequency limits transitions will be stored to the temporary file 'swapfile', 
 ! which can be later used in stead of the Transition files - this will make the simualtions faster
 !
 if (all(trim(swapfile)/=(/"null","NULL","none","NONE"/))) then
   open(unit=11,file=trim(swapfile),action='write',status='replace')
   swap = .true.
 endif
 !
 !start loop over all transitions
 !
 do i = 1,nfiles
   !
   open(unit=1,file=trim(intfilename(i)))
   do
      !   read new line from intensities file
      !
      read(1,*,end=20) ilevelf,ileveli,acoef
      !
      energyf = energies(ilevelf)
      energyi = energies(ileveli)
      !
      jf = jrot(ilevelf)
      ji = jrot(ileveli)
      !
      igammaf = sym(ilevelf)
      igammai = sym(ileveli)
      !
      tranfreq = energyf-energyi
      !
      if (tranfreq<freql.or.tranfreq>freqr) cycle 
      !
      if (swap) then 
        write(11,"(2i12,2x,es16.8)"),ilevelf,ileveli,acoef
      endif 
      !
      !   half width for Doppler profiling
      if (proftype(1:5)=='doppl') halfwidth=dpwcoef*sqrt(temp/meanmass)*tranfreq

      !   if transition frequency is out of selected range
      if (tranfreq>freqr+10.0*halfwidth.or.tranfreq<freql-10.0*halfwidth) cycle
      !
      if (tranfreq<epsilon(1.0d0)) cycle
      !
      select case (trim(specttype))
        !
      case default
        !
        print('(a,2x,a)'),'Illegal key:',trim(specttype)
        stop 'Illegal specttype-key'
        !
      case ('absorption','ABSORPTION','absorpt','ABSORPT')
        !
        ! compute absorption coefficient [cm/mol]
        !
        abscoef=cmcoef*acoef*gns(igammaf)*real(2*jf+1,8)*exp(-beta*energyi)*(1.0d0-exp(-beta*tranfreq))/(tranfreq**2*partfunc)
        !
      case ('emission','EMISSION','EMISS','emiss')
        !
        ! emission coefficient [Ergs/mol/Sr]
        !
        abscoef=emcoef*acoef*gns(igammaf)*real(2*jf+1,8)*exp(-beta*energyf)*tranfreq/(partfunc)
        !
      end select
      !
      ! if only stick spectrum needed
      if (proftype(1:5)=='stick') then
        !
        if (abscoef>absthresh) then 
          print('((1x,2es16.8,2x,a3,1x,i4,2x,4i4,2x,2i4,1x,a3,3x,2i4,1x,f12.4," <- ",&
                                 a3,1x,i4,2x,4i4,2x,2i4,1x,a3,3x,2i4,1x,f12.4))'), &
                   tranfreq,abscoef,symlab(sym(ilevelf)),jrot(ilevelf),n(ilevelf,:),l(ilevelf,:),&
                                    symlab(symv(ilevelf)),krot(ilevelf),taurot(ilevelf),energyf,&
                                    symlab(sym(ileveli)),jrot(ileveli),n(ileveli,:),l(ileveli,:),&
                                    symlab(symv(ileveli)),krot(ileveli),taurot(ileveli),energyi
        endif
        cycle
      end if
      !
      !   normalize absorption coefficient
      !
      abscoef=abscoef*sqrt(ln2/pi)/halfwidth
      !
      if (abscoef<thresh) cycle 
      !
      exp_coef = -ln2/halfwidth**2
      !
      ib =  max(nint( ( tranfreq-halfwidth*10.0-freql)/dfreq )+1,1)
      ie =  min(nint( ( tranfreq+halfwidth*10.0-freql)/dfreq )+1,npoints)
      !
      !$omp parallel do private(ipoint,dfreq_) shared(intens) schedule(dynamic)
      do ipoint=ib,ie
         dfreq_=freq(ipoint)-tranfreq
         if (abs(dfreq_)>halfwidth*10.0) cycle
         intens(ipoint)=intens(ipoint)+abscoef*exp(exp_coef*dfreq_**2)
      enddo
      !$omp end parallel do
      !
      !print('(2(1x,es16.8),i4,i4,f11.4,f11.4,f11.4,es16.8,i6,i6)'),freq(1),intens(1),ji,jf,energyi,energyf,tranfreq,linestr,ileveli,ilevelf
      !
      !ji_ = ji ; jf_ = jf
      !
      cycle
   20  exit
   enddo
   !
 enddo 

 close(1)

 if (proftype(1:5)=='doppl'.or.proftype(1:5)=='gauss') then 
   !
   !print out the generated intensities convoluted with a given profile
   !
   print('(2(1x,g16.8))'),(freq(ipoint),intens(ipoint),ipoint=1,npoints)
   !
 endif 

end program spectrum_BYTe