anchor.diff:  
                  elseif(vibpopulation_do) then
                     ! apply vibrational population
  -                  read(quantum_numbers(QN%dens_col-4,ilevelf),*) ndensity
  +                  !read(quantum_numbers(QN%dens_col-4,ilevelf),*) ndensity
                     !
                     abscoef=emcoef*acoef*gtot(ilevelf)*exp(-c2/temp*ener_rot)*ndensity*&
                             tranfreq/(partfunc)
  
positive.code:         !
          ivib = 0 
          i = -imin
          !
          do
             !
  
negative.code:               !           
                if(vibpopulation_do) then
                   ! apply vibrational population keeping the rotational part as Boltzmann
                   !read(quantum_numbers(QN%dens_col-4,i),*) ndensity
                   !
                   ndensity = dens_vib(ivib)
                   !
                   partfunc=partfunc+gtot(i)*exp(-c2/temp*ener_rot)*ndensity
                else 
                   ! split PF into a product of vib and rot parts 
  