program spectrum_10to10
 implicit none
 !
 integer, parameter :: ik          = selected_int_kind(8)       ! "Normal" integers. This must map on
                                                                ! C "int" type, or the DX interface won't
                                                                ! work correctly.
 integer, parameter :: hik         = selected_int_kind(16)      ! "Pointer" integers - sufficient to store
                                                                ! memory address
 integer, parameter :: drk         = selected_real_kind(12,25)  ! "Double" reals and complex (complexi? :-)
 integer, parameter :: rk          = selected_real_kind(12,25)  ! "Normal" reals and complex (complexi? :-)
 integer, parameter :: ark         = selected_real_kind(25,32)  ! "Accurate" reals and complex (complexi? :-)
 integer(ik), parameter :: nsym = 5
 integer(ik), parameter :: nmodes = 9
 real(rk),    parameter :: gns(5) = (/5.0,5.0,2.0,3.0,3.0/)
 !
 integer(ik) :: npoints,ipoint,j,jf,ji,ilevel,ilevelf,ileveli,igamma,igammaf,igammai,info,j0,i,ilevelf_,ileveli_,ipolyad
 integer(ik) :: ib,ie,nfiles,j_t,v(nmodes)
 integer(ik) :: verbose=4
 !
 integer(ik) :: nlevels,nn,itemp,ipartf

 real(rk) :: temp,temp0,freql,freqr,halfwidth,meanmass,partfunc,dfreq,dfreq_,acoef,abscoef,coeff,&
         energy,energyf,energyi,tranfreq,cmcoef,beta,ln2,pi,dpwcoef,absthresh,emcoef,beta0,intband,tranfreq0
 real(rk) ::  de,x0,xp,xm
 real(rk) ::  thresh = 1e-16
 real(rk),allocatable :: freq(:),intens(:),energies(:),pf(:,:)
 integer(ik),allocatable :: sym(:),n(:,:),jrot(:),krot(:),taurot(:),l(:,:),symv(:),symr(:),gtot(:)

 character(4) :: symlab(5) = (/"A1","A2","E ","F1","F2"/)
 character(60) intfilename(200),enrfilename
 character(30) proftype,specttype
 !
 logical :: swap = .false.,partfunc_do=.false.
 character(len=130)  :: my_fmt ! contains format specification for intput/output


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

!read number of intensity files
 read*,nfiles

if (nfiles>200) then 
  print('("Too many files (>200):",i9)'),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 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

halfwidth = 0 ; meanmass = 0.0 ; absthresh = 0.0 ; 
!
!read gaussian half-width or molecular mean mass
 select case (trim(proftype))
    case ('gauss'); read*,halfwidth
    case ('doppl'); read*,meanmass
    case ('stick','HITRAN','hitran'); read*,absthresh
    case ('bin '); 
       read*,halfwidth
       halfwidth = min((freqr-freql)/real(npoints,8),halfwidth)
    case default ; stop 'illegal profile-type'
 end select

!read the threshold
 read*,thresh
 if (verbose>=5) print('(1x,"thresh = ",g18.7)'),thresh



!compute constants
 ln2=log(2.0_rk)
 pi=acos(-1.0_rk)
 !beta=planck*vellgt/(boltz*temp)
 !
 beta=c2/temp
 !
 !cmcoef=avogno/(8.0d0*pi*vellgt)
 !emcoef=avogno*planck*vellgt/(4.0d0*pi)
 !
 cmcoef=1.0_rk/(8.0d0*pi*vellgt)
 emcoef=1.0_rk*planck*vellgt/(4.0d0*pi)
 dfreq=(freqr-freql)/real(npoints-1,rk)
 !
 !   half width for Doppler profiling
 dpwcoef=sqrt(2.0*ln2*boltz*avogno)/vellgt
 dpwcoef = dpwcoef*sqrt(temp/meanmass)

 !
 if (proftype(1:5)=='doppl') then
   if (verbose>=5) print('("alpha = ",f18.8)'),dpwcoef
 endif 
 !
!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'
 forall(ipoint=1:npoints) freq(ipoint)=freql+real(ipoint-1,rk)*dfreq

!count number of lines (levels) in the Energy files
 open(unit=1,file=trim(enrfilename))
 i = 0
 do
    !
    read(1,"(i9)",end=10) nn
    !
    i = i + 1
    !
    cycle
    !
10  exit
 end do
 !
 rewind(1)
 !
 nlevels = i
 !
 allocate(sym(nlevels),symr(nlevels),symv(nlevels),energies(nlevels),Jrot(nlevels),n(nlevels,9),l(nlevels,2:4),krot(nlevels),&
          taurot(nlevels),gtot(nlevels),stat=info)
 if (info/=0) stop 'error: sym,energies,Jrot 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 as temperature limits
 !
 if (trim(specttype)=='partfunc') then
   !
   if (verbose>=6) print*,halfwidth,thresh 
   !
   ipartf = int(halfwidth,4)
   !
   dfreq=temp/real(npoints,rk)
   !
   if (ipartf<0.or.ipartf>3) stop "illegal partfunc component, can be only 0,1,2,3"
   !
   write(my_fmt,'("(1x,i1,a4,1x,",i7,"( 5x,2x,f8.2,5x))")') npoints
   if (verbose>=1) print(my_fmt),ipartf,'  J ',(i*dfreq,i=1,npoints)
   allocate(pf(0:3,npoints),stat=info)
   pf = 0 ; j0=0 ; partfunc_do = .true.
   if (info/=0) stop 'error: pf'
 endif
 !
 ! prepare computing the partition function
 if (partfunc<=0) then
    partfunc=0 ; j0=0
    partfunc_do = .true.
 endif
 !
 ! start reading the energies 
 !
 write(my_fmt, '("(1x,i4,3(1x,",i8,"es20.8))")') npoints
 !
 do i=1,nlevels
    !
    read(1,*) nn,energies(i),gtot(i),jrot(i),sym(i),n(i,1:9),symv(i),j_t,krot(i),taurot(i),symr(i),ipolyad,coeff,v(1:9)
    !
    j = jrot(i)
    igamma = sym(i)
    energy = energies(i)
    !
    if (partfunc_do) then
      if (j/=j0) then 
         !
         if (trim(specttype)=='partfunc') then 
           if (verbose>=1) print(my_fmt),j0,pf(ipartf,1:npoints)
         else 
           if (verbose>=3) print('("#",i4,1x,es16.8)'),j0,partfunc
         endif
         !
      endif
      !
      ! symmetry number
      ! partition function
      !
      !partfunc=partfunc+real((2*j+1),rk)*gns(igamma)*exp(-beta*energy)
      !
      partfunc=partfunc+gtot(i)*exp(-beta*energy)
      !
      j0 = j
     endif
     !
     if (trim(specttype)=='partfunc') then
       !
       do itemp = 1,npoints
         !
         temp0 = real(itemp,rk)*dfreq
         !
         beta0 = planck*vellgt/(boltz*temp0)
         !
         ! apart from the partition function, the first two moments and the heat capacity are also computed:
         !
         pf(0,itemp) = pf(0,itemp) + real((2*j+1),rk)*                  gns(igamma)*exp(-beta0*energy)
         pf(1,itemp) = pf(1,itemp) + real((2*j+1),rk)*(beta0*energy)   *gns(igamma)*exp(-beta0*energy)
         pf(2,itemp) = pf(2,itemp) + real((2*j+1),rk)*(beta0*energy)**2*gns(igamma)*exp(-beta0*energy)
         !
         !qp = temp0*pf(1,itemp)
         !qpp= temp0**2*pf(2,itemp)+2.0d0*qp
         pf(3,itemp) = ( pf(2,itemp)/pf(0,itemp)-(pf(1,itemp)/pf(0,itemp))**2 )
         !
       enddo
       !
     endif
     !
 end do
 close(1)
 !
 ! print out the computed part. function objects and finish
 !
 if (trim(specttype)=='partfunc') then
   !
   write(my_fmt, '("1x,i4,3(1x,",i8,"es20.8)")') npoints
   if (verbose>=1) print(my_fmt),j0,pf(ipartf,1:npoints)
   !
   if (verbose>=1) print(my_fmt),j0,pf(0,1:npoints)
   if (verbose>=1) print(my_fmt),j0,pf(1,1:npoints)
   if (verbose>=1) print(my_fmt),j0,pf(2,1:npoints)
   if (verbose>=1) print(my_fmt),j0,pf(3,1:npoints)
   !
   if (verbose>=2) print('(1x,a,1x,es16.8)'),'! partition function value is',partfunc
   !
   stop
   !
 else
   !
   if (verbose>=2) print('("#",i4,1x,es16.8)'),j0,partfunc
   !
 endif
 !
 if (verbose>=2) print('(1x,a,1x,es16.8)'),'# partition function value is',partfunc
 !
 intens=0.0
 !
 !start loop over all transitions
 !
 intens=0.0
 intband = 0.0
 !
 !ji_ = 0 ; jf_ = 0
 !
 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
      !
      !sprint*,energyf,energyi,tranfreq,freql,freqr
      !
      tranfreq0 = tranfreq
      !
      !   half width for Doppler profiling
      if (proftype(1:5)=='doppl') then
        halfwidth=dpwcoef*tranfreq
      endif 
      !
      !if (tranfreq0<freql.or.tranfreq0>freqr) cycle 
      !
      !   if transition frequency is out of selected range
      if (tranfreq0>freqr+10.0*halfwidth.or.tranfreq0<freql-10.0*halfwidth) cycle
      !
      x0 = sqrt(ln2)/halfwidth*dfreq*0.5_rk
      !
      if (tranfreq0<epsilon(1.0_rk)) cycle
      !
      if (gns(igammai)/=gns(igammaf)) then
         print('("gns are not the same for states ",2(1x,i9)," : ",2i7)'), ileveli,ilevel,gns(igammai),gns(igammaf)
         stop 'gns are not the same!'
      endif
      !
      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,rk)*exp(-beta*energyi)*(1.0_rk-exp(-beta*tranfreq))/(tranfreq**2*partfunc)
        !
        abscoef=cmcoef*acoef*gtot(ileveli)*exp(-beta*energyi)*(1.0_rk-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,rk)*exp(-beta*energyf)*tranfreq/(partfunc)
        !
      end select
      !
      intband = intband + abscoef
      !
      ! if only stick spectrum needed
      if (proftype(1:5)=='stick') then
        !
        if (abscoef>absthresh) &
          print('((1x,2es16.8,2x,a3,1x,i3,2x,9i3,1x,a3,1x,2i4,1x,a3,1x,f12.4," <- ",'//&
            &'a3,1x,i3,2x,9i3,1x,a3,1x,2i4,1x,a3,1x,f12.4))'),&
            tranfreq,abscoef,symlab(sym(ilevelf)),jrot(ilevelf),n(ilevelf,1:9),symlab(symv(ilevelf)),krot(ilevelf),taurot(ilevelf),&
            symlab(symr(ilevelf)),energyf,symlab(sym(ileveli)),jrot(ileveli),n(ileveli,1:9),symlab(symv(ileveli)),&
            krot(ileveli),taurot(ileveli),symlab(symr(ileveli)),energyi
        cycle
      end if
      !
      if (abscoef<thresh) cycle 
      !
      ib =  max(nint( ( tranfreq0-halfwidth*10.0-freql)/dfreq )+1,1)
      ie =  min(nint( ( tranfreq0+halfwidth*10.0-freql)/dfreq )+1,npoints)
      !
      !   normalize absorption coefficient
      !
      !read gaussian half-width or molecular mean mass
      select case (trim(proftype(1:5)))
          !
      case ('gauss','doppl')
          !
          !abscoef=abscoef*sqrt(ln2/pi)/halfwidth
          !
          !$omp parallel do private(ipoint,dfreq_,xp,xm,de) shared(intens) schedule(dynamic)
          do ipoint=ib,ie
             !
             dfreq_=freq(ipoint)-tranfreq0
             !if (abs(dfreq_)>halfwidth*10.0) cycle
             !
             xp = sqrt(ln2)/halfwidth*(dfreq_)+x0
             xm = sqrt(ln2)/halfwidth*(dfreq_)-x0
             !
             de = erf(xp)-erf(xm)
             !
             intens(ipoint)=intens(ipoint)+abscoef*0.5_rk/dfreq*de
             !
          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
          !
      case ('bin ');
        !
        !$omp parallel do private(ipoint,dfreq_) shared(intens) schedule(dynamic)
        do ipoint=ib,ie
           dfreq_=freq(ipoint)-tranfreq0
           !if (abs(dfreq_)>halfwidth*10.0) cycle
           intens(ipoint)=max(intens(ipoint),abscoef)
        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
          !
      end select
      !
      cycle
   20  exit
   enddo
   !
 enddo 

 close(1)

 if (proftype(1:5)=='doppl'.or.proftype(1:5)=='gauss'.or.proftype(1:4)=='bin') then 
   !
   !print out the generated intensities convoluted with a given profile
   !
   print('(2(1x,es16.8))'),(freq(ipoint),intens(ipoint),ipoint=1,npoints)
   !
   print('("Total intensity  (sum):",es16.8," (int):",es16.8)'), intband,sum(intens)*dfreq
   !
 else
   !
   print('("Total intensity = ",es16.8)'), intband
   !
 endif 
 !
end program spectrum_10to10