anchor.diff:       !
      endif
      !
  -   allocate(quantum_numbers(maxitems,iline),stat=info)
  -   call ArrayStart('quantum_numbers',info,size(quantum_numbers),kind(quantum_numbers))
  -   !
  -   if (verbose>=3) call MemoryReport
  -   !
  -   quantum_numbers = ""
  -   !
  -   indices = 0
  -   !
  -   energies = 0
  -   Jrot = 0
  -   gtot = 0
  -   !
  -   ! 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 (specttype(1:4)=='PART') then
  -     !
  -     dtemp=maxtemp/real(npoints,rk)
  -     !
  -     if (ipartf<0.or.ipartf>3) stop "illegal partfunc component, can be only 0,1,2,3"
  -     !
  -     if (verbose>=1) print("('0,1,2,3 stand for PF, 1st and 2d moments and Cp:')")
  -     !
  -     if (verbose>=1) print('("!",5x,a4,( 1x,'//npoints_fmt//'( 5x,"T=",f8.2,5x) ) )'),'  J ',(i*dtemp,i=1,npoints)
  -     !
  -     allocate(pf(0:3,npoints),stat=info)
  -     call ArrayStart('pf',info,size(pf),kind(pf))
  -     !
  -     pf = 0 
  -     !
  -   endif
  -   j0 = 0
  -   i = 0
  -   !
  -   ! start reading the energies 
  -   !
  -   do
  -      !
  -      call read_line(eof,enunit) ; if (eof) exit
  -      i = i + 1
  -      !
  -      call readi(itemp)
  -      call readf(energy)
  -      !
  -      if (energy>enermax) then
  -        i = i - 1
  -        cycle
  -      endif
  -      !
  -      energies(i) = energy
  -      !
  -      call readi(gtot(i))
  -      call readf(jrot(i))
  -      !
  -      do kitem = 5,nitems
  -         !
  -         call reada(quantum_numbers(kitem-4,i))
  -         !
  -      enddo
  -      !
  -      if (filter.and.i==1) then
  -        !
  -        do ifilter = 1,Nfilters
  -          !
  -          if (lower(ifilter)%i<1.or.lower(ifilter)%i>nitems-4.or.upper(ifilter)%i<1.or.upper(ifilter)%i>nitems-4) then
  -            !
  -            write(out,"('wrong filter indices, upper or lower: ',2i)") upper(ifilter)%i+4,lower(ifilter)%i+4
  -            print*,quantum_numbers(:,i)
  -            stop 'wrong filter indices, upper or lower'
  -            !
  -          endif
  -          !
  -        enddo
  -        !
  -      endif 
  -      !
  -      indices(itemp) = i
  -      !
  -      j = 2*int(jrot(i))+1
  -      energy = energies(i)
  -      !
  -      if (partfunc_do) then
  -        if (j/=j0) then 
  -           !
  -           if (specttype(1:4)=='PART') then 
  -             if (verbose>=1.and.npoints<1000) print('("!",f8.1,3(1x,'//npoints_fmt//'es20.8))'),jrot(i),pf(ipartf,1:npoints)
  -           else 
  -             if (verbose>=4) print('("|",f8.1,1x,es16.8)'),jrot(i),partfunc
  -           endif
  -           !
  -        endif
  -        !
  -        partfunc=partfunc+gtot(i)*exp(-beta*energy)
  -        !
  -        j0 = j
  -        !
  +   ! open and read broenning files
  +   !
  +   do i =1,Nspecies
  +     !
  +     if ( trim(species(i)%filename)/="" ) then
  +       !
  +       write(ioname, '(a)') 'brodenning file'
  +       !
  +       call IOstart(trim(ioname),bunit)
  +       open(unit=bunit,file=trim(species(i)%filename),action='read',status='old')
  +       !
  +       Jmax = 0 
  +       !
  +       do
  +         !
  +         ! Scan and find Jmax
  +         read(bunit,*,end=14) gamma_,n_,J_
  +         !
  +         Jmax = max(Jmax,J_)
  +         !
  +         cycle 
  +         14  exit
  +         !
  +       enddo
  +       !
  +       if (Jmax>0) then
  +         !
  +         allocate(species(i)%gammaQN(0:Jmax,-1:1),stat=info)
  +         call ArrayStart('gammaQN',info,size(species(i)%gammaQN),kind(species(i)%gammaQN))
  +         species(i)%gammaQN(:,:) = species(i)%gamma
  +         !
  +         allocate(species(i)%nQN(0:Jmax,-1:1),stat=info)
  +         call ArrayStart('nQN',info,size(species(i)%nQN),kind(species(i)%nQN))
  +         species(i)%nQN(:,:) = species(i)%N
  +         !
  +        else
  +         !
  +         write(out,"('Jmax = 0 in ',a)") trim(species(i)%filename)
  +         stop 'Jmax = 0 in a broadening file'
  +         !
          endif
          !
  -       if (specttype(1:4)=='PART') then
  -         !
  -         do itemp = 1,npoints
  -           !
  -           if (energy>enermax) cycle
  -           !
  -           temp0 = real(itemp,rk)*dtemp
  -           !
  -           beta0 = planck*vellgt/(boltz*temp0)
  -           !
  -           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
  -           !
  -           pf(3,itemp) = ( pf(2,itemp)/pf(0,itemp)-(pf(1,itemp)/pf(0,itemp))**2 )
  -           !
  -         enddo
  -         !
  -       endif
  -       !
  -   end do
  -   close(enunit)
  -   !
  -   ! print out the computed part. function objects and finish
  -   !
  -   if (specttype(1:4)=='PART') then
  -     !
  -     j0rk = real(j0-1)*0.5_rk
  -     !
  -     if (verbose>=0) print('("!",f8.1,3(1x,'//npoints_fmt//'es20.8/))'),j0rk,pf(ipartf,1:npoints)
  -     !
  -     if (verbose>=0) print('("!0",8x,'//npoints_fmt//'es20.8)'),pf(0,1:npoints)
  -     if (verbose>=0) print('("!1",8x,'//npoints_fmt//'es20.8)'),pf(1,1:npoints)
  -     if (verbose>=0) print('("!2",8x,'//npoints_fmt//'es20.8)'),pf(2,1:npoints)
  -     if (verbose>=0) print('("!3",8x,'//npoints_fmt//'es20.8)'),pf(3,1:npoints)
  -     !
  -     if (verbose>=3) print('(/a,1x,es16.8/)'),'! partition function value is',partfunc
  -     !
  -     if (verbose>=0) then 
  -       print('(a)'),"Partition function"
  -       do itemp = 1,npoints
  -         temp0 = real(itemp,rk)*dtemp
  -         print('(f12.2,1x,es20.8)'),temp0,pf(0,itemp)
  +       rewind(bunit)
  +       !
  +       ! read in the QN-dependent broadening parameters
  +       !
  +       do
  +         !
  +         ! model specific read
  +         !
  +         select case ( trim(species(i)%model) )
  +           !
  +         case ('J','A0')
  +           !
  +           read(bunit,*,end=15) gamma_,n_,J
  +           species(i)%gammaQN(J,:) = gamma_
  +           species(i)%nQN(J,:) = n_
  +           !
  +         case ('JJ','A1')
  +           !
  +           read(bunit,*,end=15) gamma_,n_,Jpp,Jp
  +           !
  +           if (abs(Jp-Jpp)>1) stop 'Jp-Jpp in broadening file'
  +           !
  +           species(i)%gammaQN(Jpp,Jp-Jpp) = gamma_
  +           species(i)%nQN(Jpp,Jp-Jpp) = n_
  +           !
  +         case ('JJKK-X')
  +           !
  +           read(bunit,*,end=15) gamma_,n_,Jpp,Jp,Kp,Kpp
  +           !
  +           if (abs(Jp-Jpp)>1) stop 'Jp-Jpp in broadening file'
  +           if (Kp>Jp.or.Kpp>Jpp) stop 'K > J in broadening file'
  +           !
  +           !species(i)%gammaQN(Jpp,Jp-Jpp,Kp,Kpp) = gamma_
  +           !species(i)%nQN(Jpp,Jp-Jpp,Kp,Kpp) = n_
  +           !
  +         case default 
  +           !
  +           write(out,"('The broadening model',a,' is not implemented yet')") trim(species(i)%model)
  +           stop 'Illegal broadening model'
  +           !
  +         end select 
  +         !
  +         cycle 
  +         15  exit
  +         !
          enddo
  -     endif 
  -     !
  -     stop
  -     !
  -   else
  -     !
  -     if (verbose>=4) print('("|",f9.1,1x,es16.8)'),real(j0-1)*0.5,partfunc
  -     !
  -   endif
  -   !
  -   if (verbose>=2.or.(partfunc_do.and.verbose>0)) print('(1x,a,1x,es16.8/)'),'! partition function value is',partfunc
  -   !
  +       !
  +     endif
  +     !
  +     close(bunit)
  +     ! 
  +   enddo
      !
      ! Intensities
      !
  
positive.code:   integer(ik)    :: Nspecies = 0, Nfilters = 0
    type(speciesT) :: species(nspecies_max)
    !
    logical :: partfunc_do = .true., filter = .false., histogram = .false.
    !
    contains
    !
  
negative.code:             !
              alpha = halfwidth0
              gamma = halfwidth
              !
              !call voi_quad(npoints,ib,ie,freq,abscoef,intens,dfreq,tranfreq,alpha,gamma,nquad,abciss,weight)
              !
              sigma = alpha/ln22
              !
  