anchor.diff:           !
            maxitems = max(nitems,maxitems)
            !
  -         !cycle
  +         if (vibtemperature_do) then 
  +            call readf(dtemp)
  +            call readf(ji)
  +            !
  +            if (info>0.and.nint(ji)==0) then
  +              write(out,"('It is illegal to use not J-sorted states with vib-temperatures')")
  +              stop 'States file must be ordered by J to use with Tvib'
  +            endif
  +            !
  +            if (info>0) cycle
  +            !
  +            if (nint(ji)==0) then 
  +              Nvib_states = Nvib_states + 1
  +            else 
  +              info = 1
  +            endif
  +            !
  +         endif
            !
            !10  exit
         end do
  
positive.code:    integer(ik),allocatable :: indexi_RAM(:),indexf_RAM(:)
     integer(ik),allocatable :: ileveli_RAM(:),ilevelf_RAM(:)
     integer(ik) :: nswap_,nswap,iswap,iswap_,iomp
     !
     integer(ik),allocatable :: nchars_quanta(:)  ! max number of characters used for each quantum number
     !
  
negative.code:                !
               case ('ABSORPTION')
                 !
                 !abscoef=cmcoef*acoef*gtot(ilevelf)*exp(-beta*energyi)*(1.0_rk-exp(-beta*tranfreq))/(tranfreq**2*partfunc)
                 !
                 if (.not.vibtemperature_do) then
                    abscoef=cmcoef*acoef*gtot(ilevelf)*exp(-beta*energyi)*(1.0_rk-exp(-beta*tranfreq))/(tranfreq**2*partfunc)
                 else 
                    ! split into a product of vib and rot parts 
                    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 
                 !
               case ('emission','EMISSION','EMISS','emiss')
                 !
                 ! emission coefficient [Ergs/mol/Sr]
                 !
                 if (.not.vibtemperature_do) then
                    abscoef=emcoef*acoef*gtot(ilevelf)/real(2*ji+1,rk)*real(2*jf+1,rk)*exp(-beta*energyf)*tranfreq/(partfunc)
                 else 
                    ! split into a product of vib and rot parts 
                    abscoef=emcoef*acoef*gtot(ilevelf)/real(2*ji+1,rk)*real(2*jf+1,rk)*exp(-c2/temp*ener_rot)*exp(-c2/temp_vib*ener_vib)*tranfreq/(partfunc)
                    !
                 endif 
                 !
               case ('LIFETIME')
                 !
  