program spectrum
 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) :: isym,nsym,npoints,ipoint,j,jf,ji,ilevel,ilevelf,ileveli,igamma,igammaf,igammai,info,j0,i,ilevelf_,ileveli_
 integer(ik) :: quantaf(0:10),quantai(0:10),imode,nmodes,localf(0:10),locali(0:10) !,normalf(0:10),normali(0:10)
 integer(ik) :: ib,ie,nfiles,j_t,v(9),ilog,ienerlow
 integer(ik) :: verbose=3,smax=200,nchar
 integer(hik):: Nintens(-60:60)
 !
 integer(ik) :: nlevels,nn,itemp,ipartf,NNN
 real(rk)    :: int_vs_enerlow(0:12000)

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

 character(4) symlab(10),a_fmt
 character(1) invchar
 character(60) intfilename(200),enrfilename,swapfile
 character(30) proftype,specttype
 character(100),allocatable :: state_string(:)
 character(300) :: string_tmp
 character(len=200) :: my_fmt  !text variable containing formats for reads/writes
 character(len=9) :: npoints_fmt  !text variable containing formats for reads/writes
 !
 logical :: swap = .false.,partfunc_do=.false.,conv_units = .false. 


 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):",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 characters
  read*,nchar

 write(a_fmt,'(I3)') nchar
 !print*,a_fmt
 write(a_fmt,'("a",a3)') adjustl(a_fmt)
 !
 write(my_fmt,'(a,a4,a,a4,a)')  '(1x,2es16.8,1x,f5.1,1x,f12.4," <- ",f5.1,1x,f12.4,5x,',a_fmt,',"<-  ",',a_fmt,'))'
 !print*,my_fmt
 !
!read nuclear spin statistics weights
! read*,(symlab(isym),gns(isym),isym=1,nsym)

!read initial state quantum numbers 
! read*,(quantai(imode),imode=0,nmodes)

!read final state quantum numbers 
! read*,(quantaf(imode),imode=0,nmodes)

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

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

 write(npoints_fmt,'(i9)') 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','gauss-um'); read*,halfwidth
    case ('doppl','doppl-um'); read*,meanmass
    case ('stick','HITRAN','hitran','trunck'); read*,absthresh
    case ('rect ','bin  '); 
       read*,halfwidth
       halfwidth = min((freqr-freql)/real(npoints,8),halfwidth)
    case default ; stop 'illegal profile-type'
 end select

!read the threshold
 read*,thresh


!read the threshold
! read*,elow 


!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(6:8)=='-um') conv_units = .true.
 !
 if (proftype(1:5)=='doppl') then
   if (verbose>=4) print('("alpha = ",f18.8)'),dpwcoef
 endif 
 !
 if (proftype(1:5)=='gauss') then
     !
     if (dfreq>halfwidth*100.0) then 
       !
       if (verbose>=4) print('(" dnu  >> alpha")')
       !
     elseif(halfwidth>dfreq*100.0) then 
       !
       if (verbose>=4) print('(" dnu  << alpha")')
       !
     else
       !
       if (verbose>=4) print('(" dnu  ~ alpha")')
       !
    endif
 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,"(i8)",end=10) nn
    !
    i = i + 1
    !
    cycle
    !
10  exit
 end do
 !
 rewind(1)
 !
 nlevels = i
 !
 allocate(sym(nlevels),energies(nlevels),Jrot(nlevels),state_string(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
   !
   !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"
   !
   if (verbose>=1) print('("!",i1,a4,( 1x,'//npoints_fmt//'( 5x,"T=",f8.2,5x) ) )'),ipartf,'  J ',(i*dfreq,i=1,npoints)
   !
   allocate(pf(0:3,npoints),stat=info)
   if (info/=0) stop 'error: pf'
 endif
 !
 ! prepare computing the partition function
 if (partfunc<=0) then
    partfunc=0 ; j0=0
    partfunc_do = .true.
    if (trim(specttype)=='partfunc') pf = 0
 endif
 !
 ! start reading the energies 
 !
 do i=1,nlevels
    !
    !read(1,*) nn,energies(i),gtot(i),jrot(i),invchar,sym(i),NNN,n(i,1:7),j_t,krot(i),taurot(i),v(1:6),symv(i)
    !
    !read(1,*) nn,energies(i),gtot(i),jrot(i),string_tmp,state_string(i)(1:nchar)
    !
    read(1,"(a250)") string_tmp
    !
    read(string_tmp,*) nn,energies(i),gtot(i),jrot(i)
    read(string_tmp(42:),"("//a_fmt//")") state_string(i)(1:nchar)
    !
    !print("("//a_fmt//")"),adjustl(string_tmp(42:min(nchar,250)))
    !
    j = 2*jrot(i)+1
    igamma = sym(i)
    energy = energies(i)
    !
    if (partfunc_do) then
      if (j/=j0) then 
         !
         if (trim(specttype)=='partfunc') then 
           if (verbose>=1.and.npoints<1000) print('("!",i4,3(1x,'//npoints_fmt//'es20.8))'),j0,pf(ipartf,1:npoints)
         else 
           if (verbose>=4) 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
         !
         if (energy>freqr) cycle
         !
         temp0 = real(itemp,rk)*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),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)
         !
         pf(0,itemp) = pf(0,itemp) + gtot(i)*exp(-beta0*energy)
         pf(1,itemp) = pf(1,itemp) + gtot(i)*exp(-beta0*energy)*(beta0*energy)
         pf(2,itemp) = pf(2,itemp) + gtot(i)*exp(-beta0*energy)*(beta0*energy)**2
         !
         !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
   !
   if (verbose>=1) print('("!",i4,3(1x,'//npoints_fmt//'es20.8/))'),j0,pf(ipartf,1:npoints)
   !
   if (verbose>=1) print('("!0",i4,1x,'//npoints_fmt//'es20.8)'),j0,pf(0,1:npoints)
   if (verbose>=1) print('("!1",i4,1x,'//npoints_fmt//'es20.8)'),j0,pf(1,1:npoints)
   if (verbose>=1) print('("!2",i4,1x,'//npoints_fmt//'es20.8)'),j0,pf(2,1:npoints)
   if (verbose>=1) print('("!3",i4,1x,'//npoints_fmt//'es20.8)'),j0,pf(3,1:npoints)
   !
   if (verbose>=3) print('(1x,a,1x,es16.8)'),'! partition function value is',partfunc
   !
   if (verbose>=1) then 
     print('(a)'),"Partition function"
     do itemp = 1,npoints
       temp0 = real(itemp,rk)*dfreq
       print('(f12.2,1x,es20.8)'),temp0,pf(0,itemp)
     enddo
   endif 
   !
   stop
   !
 else
   !
   if (verbose>=4) print('("#",i4,1x,es16.8)'),j0,partfunc
   !
 endif
 !
 if (verbose>=2.or.partfunc_do) 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
 !
 intens=0.0
 intband = 0.0
 !
 !ji_ = 0 ; jf_ = 0
 !
 ilevelf_ = 0
 ileveli_ = 0
 !
 Nintens = 0
 int_vs_enerlow = 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
      !
      !read(1,*,end=20) ileveli,ilevelf,acoef
      !
      energyf = energies(ilevelf)
      energyi = energies(ileveli)
      !
      if (energyf-energyi<0) then 
        write(6,"(2i12,2x,3es16.8)"),ilevelf,ileveli,acoef,energyf,energyi
        stop 'wrong order of indices'
      endif 
      !
      if (elow<energyi) cycle
      !
      jf = jrot(ilevelf)
      ji = jrot(ileveli)
      !
      !igammaf = sym(ilevelf)
      !igammai = sym(ileveli)
      !
      tranfreq = energyf-energyi
      !
      !sprint*,energyf,energyi,tranfreq,freql,freqr
      !
      tranfreq0 = tranfreq
      if (conv_units) tranfreq0 = 10000.0_rk/tranfreq
      !
      !   half width for Doppler profiling
      if (proftype(1:5)=='doppl') then
        halfwidth=dpwcoef*tranfreq
        if (conv_units) halfwidth = 10000.0_rk/tranfreq*dpwcoef
      endif 
      !
      if (tranfreq0>freqr+10.0*halfwidth.or.tranfreq0<freql-10.0*halfwidth) cycle

      !
      if (swap) then 
        write(11,"(2i12,2x,es16.8)"),ilevelf,ileveli,acoef
      endif 
      !
      x0 = sqrt(ln2)/halfwidth*dfreq*0.5_rk
      !
      if (tranfreq0<epsilon(1.0_rk)) cycle
      !
      !localf(0) = symv(ilevelf)
      !locali(0) = symv(ileveli)
      !
      !localf(1:9) = n(ilevelf,1:9)
      !locali(1:9) = n(ileveli,1:9)
      !
      !do imode = 0,nmodes
      !  !
      !  if ( quantaf(imode)==-1 ) localf(imode) = -1
      !  if ( quantai(imode)==-1 ) locali(imode) = -1
      !  !
      !enddo
      !
      !if ( any( localf(:nmodes)/=quantaf(:nmodes)).or.any(locali(:nmodes)/=quantai(:nmodes) ) ) 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,rk)*exp(-beta*energyi)*(1.0_rk-exp(-beta*tranfreq))/(tranfreq**2*partfunc)
        !
        abscoef=cmcoef*acoef*gtot(ilevelf)*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)
        !
        abscoef=emcoef*acoef*gtot(ilevelf)/real(2*ji+1,rk)*real(2*jf+1,rk)*exp(-beta*energyf)*tranfreq/(partfunc)
        !
      end select
      !
      intband = intband + abscoef
      !
      ilog=log10(abscoef)
      !
      ilog = max(min(ubound(Nintens,dim=1),ilog),lbound(Nintens,dim=1))
      !
      Nintens(ilog) = Nintens(ilog)+1
      !
      ! build intensity vs lower energy
      ienerlow = nint(energyi)
      !
      ienerlow = max(min(ubound(int_vs_enerlow,dim=1),ienerlow),lbound(int_vs_enerlow,dim=1))
      !
      !int_vs_enerlow(ienerlow) = max(abscoef,int_vs_enerlow(ienerlow))
      !
      int_vs_enerlow(ienerlow) = (abscoef+int_vs_enerlow(ienerlow))
      !
      !
      !
      ! if only stick spectrum needed
      if (proftype(1:5)=='stick') then
         !
         !
         if (abscoef>absthresh)  print(my_fmt), &
         tranfreq,abscoef,jrot(ilevelf),energyf, jrot(ileveli),energyi,state_string(ilevelf)(1:nchar),&
         state_string(ileveli)(1:nchar)


         cycle
      end if
      !
      if (abscoef<thresh) cycle 
      !
      if (proftype(1:6)=='trunck') then
         if (abscoef>absthresh) write(1000+int(tranfreq/100),"(2i12,2x,es16.8,2x,f16.6)"),ilevelf,ileveli,acoef,tranfreq
         cycle
      end if
      !
      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 
          !
      case ('rect');
        !
        !$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
        !
      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)=intens(ipoint)+abscoef/dfreq
           !
        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)=='rect'.or.proftype(1:3)=='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 
 !
 if (verbose>=3) print('(/"Intensity distribution: ")')
 !
 do i = lbound(Nintens,dim=1),ubound(Nintens,dim=1)
   !
   if (verbose>=3) print('(i4,2x,i10)'), i,Nintens(i)
   !
 enddo
 !
 !do i = lbound(int_vs_enerlow,dim=1),ubound(int_vs_enerlow,dim=1)
 !  !
 !  if (verbose>=3) print('(f8.1,2x,es16.8)'), real(i),int_vs_enerlow(i)
 !  !
 !enddo
 !
end program spectrum