anchor.diff:                 !
                  if (.not.vibtemperature_do) then
                     abscoef=cmcoef*acoef*gtot(ilevelf)*exp(-beta*energyi)*(1.0_rk-exp(-beta*tranfreq))/(tranfreq**2*partfunc)
  -               else 
  +               elseif(vibpopulation_do) then
  +                  ! apply vibrational population
  +                  read(quantum_numbers(QN%dens_col-4,ileveli),*) ndensity
  +                  abscoef=cmcoef*acoef*gtot(ilevelf)*exp(-c2/temp*ener_rot)*ndensity*&
  +                          (1.0_rk-exp(-c2/temp_vib*tranfreq))/(tranfreq**2*partfunc)
                     ! split into a product of vib and rot parts 
  +               else
                     abscoef=cmcoef*acoef*gtot(ilevelf)*exp(-c2/temp*ener_rot)*exp(-c2/temp_vib*ener_vib)*&
                             (1.0_rk-exp(-c2/temp_vib*tranfreq))/(tranfreq**2*partfunc)
                  endif 
  
positive.code:                 QN%Rotcol(2) = QN%Rotcol(1)
                  !
                endif
                !
              case ("VIB","VIBRATIONAL")
                !
  
negative.code:        !
         !
         if (abscoef>1.0D-99) then 
            write(h_fmt,"(a)") "(i2,i1,f12.6,es10.3,e10.3,f5.4,f5.3,f10.4,f4.2,f8.6)"
         else
            write(h_fmt,"(a)") "(i2,i1,f12.6,es10.2E3,e10.3,f5.4,f5.3,f10.4,f4.2,f8.6)"
         endif
         !
         if (stick_hitran) then
  