anchor.diff:     !
      logical :: eof
      character(len=cl) :: w
  +   !
  +   integer(ik) :: imol
  +   real(rk)    :: gf,gi
  +   character(55) ch_q
      !
      ln2=log(2.0_rk)
      ln22 = ln2*2.0_rk
  
positive.code:             !
              alpha = halfwidth0
              gamma = halfwidth
              !
              sigma = alpha/ln22
              !
  
negative.code:    !
     forall(ipoint=1:npoints) freq(ipoint)=freql+real(ipoint-1,rk)*dfreq
     !
     ! open and count number of lines (levels) in the Energy files
     !
     if (trim(enrfilename)/="none".and..not.hitran) then
        !
        write(ioname, '(a)') 'Energy file'
        call IOstart(trim(ioname),enunit)
        open(unit=enunit,file=trim(enrfilename),action='read',status='old')
        i = 0
        iline = 0
        maxitems = 5
        !
        call input_options(echo_lines=.false.,error_flag=1)
        !
        do
           !
           call read_line(eof,enunit) ; if (eof) exit
           !
           if (nitems<4) then 
             call report ("In .states less than  than 4 columns, mixed up with .trans?"//trim(w),.true.)
           endif 
           !
           call readi(itemp)
           call readf(energy)
           !
           if (energy>enermax) cycle
           !
           iline = iline + 1
           !
           i = max(i,itemp)
           !
           maxitems = max(nitems,maxitems)
           !
           !cycle
           !
           !10  exit
        end do
        !
        nlines = iline
        !
        if (verbose>=3) write(out,"('Max number of energy records = ',i9)") maxitems
        !
        open(unit=enunit,file=trim(enrfilename),action='read',status='old')
        !
        rewind(enunit)
        !
        nlevels = i
        !
        maxitems = maxitems-4
        !
        if (verbose>=3) print*,"nlevels: (total)",nlevels," (selected)",nlines
        !
        allocate(energies(nlines),Jrot(nlines),gtot(nlines),indices(nlevels),stat=info)
        call ArrayStart('energies',info,size(energies),kind(energies))
        call ArrayStart('Jrot',info,size(Jrot),kind(Jrot))
        call ArrayStart('gtot',info,size(gtot),kind(gtot))
        call ArrayStart('indices',info,size(indices),kind(indices))
        !
        if (trim(specttype)=='LIFETIME') THEN 
          !
          ! allocate the matrix for the sum of the A-coeffs 
          !
          allocate(Asum(nlines),stat=info)
          call ArrayStart('Asum',info,size(Asum),kind(Asum))
          !
          Asum = -1.0_rk
          !
          write(my_fmt,'(a)')  '(1x,i11,1x,f12.4,1x,f5.1,1x,es16.8,5x)'
          !
        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
          !
          if (trim(enrfilename)=="none") then 
             !
             write(out,"('PARTITION FUNCTION ERROR: no energy file is specified')") 
             stop "no energy file is specified"
             !
          endif
          !
          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
             !
            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)
            enddo
          endif 
          !
          stop
          !
        else
          !
          if (verbose>=4) print('("|",f9.1,1x,es16.8)'),real(j0-1)*0.5,partfunc
          !
        endif
        !
     endif  ! end of processing the states files 
     !   
     if (verbose>=2.or.(partfunc_do.and.verbose>0)) print('(1x,a,1x,es16.8/)'),'! partition function value is',partfunc
     !
     ! prepare the quadratures (Gauss-Hermite)
     !
  