anchor.diff:          !
           Nintens(ilog) = Nintens(ilog)+1
           !
  -        if (abscoef<thresh) cycle
  +        !
  +        ! HITAN cutoff model as an intensity threshold
  +        if (trim(cutoff_model)=="HITRAN") then
  +           if (tranfreq<=nu_crit) then
  +             cutoff = S_crit*tranfreq/nu_crit*tanh(0.5_rk*c2/temp*tranfreq)
  +           else
  +             cutoff = S_crit
  +           endif
  +        else
  +           cutoff = thresh
  +        endif
  +        !
  +        if (abscoef<cutoff) cycle
           !
           ! if only stick spectrum needed
           if (proftype(1:5)=='STICK') then
  
positive.code:    integer(ik) :: indexf_,indexi_
     real(rk)    :: beta,ln2,ln22,dtemp,dfreq,temp0,beta0,intband,dpwcoef,x0,tranfreq,tranfreq_i,abscoef,dfreq_,xp,xm,de,lor,b,lor2,dfreq_2,halfwidth0,dnu_half,dxp
     real(rk)    :: cmcoef,emcoef,energy,energyf,energyi,elow,jf,ji,acoef,j0rk,bnorm,f,eta1,eta2,intens1,intens2,Va0,gammaV,a,wg,d
     real(rk)    :: sigma,alpha,gamma,y,x1,x2,voigt_,dx2,xi,L1,L2,acoef_
     integer(ik) :: Jmax,Jp,Jpp,Kp,Kpp,J_,ispecies,alloc_p
     real(rk)    :: gamma_,n_
     character(len=cl) :: ioname
  
negative.code:   integer(ik)   :: GNS=1,npoints=1001,nchar=1,nfiles=1,ipartf=0,verbose=3,ioffset = 10,iso=1
    real(rk)      :: temp=298.0,partfunc=-1.0,freql=-small_,freqr= 20000.0,thresh=1.0d-90,halfwidth=1e-2,meanmass=1.0,maxtemp=10000.0
    real(rk)      :: voigt_gamma = 0.05, voigt_n = 0.44, offset = -25.0, pressure = 1.0_rk 
    real(rk)      :: enermax = 1e6, abscoef_thresh = 1.0d-50, abundance = 1.0d0
    real(rk)      :: S_crit = 1e-29      ! cm/molecule, HITRAN cut-off paramater 
    real(rk)      :: nu_crit = 2000.0d0  ! cm-1, HITRAN cut-off paramater 
    character(len=cl) :: cutoff_model = "NONE"
    integer(ik)   :: nquad = 20      ! Number of quadrature points 
    !
    character(len=cl) :: specttype="absorption",enrfilename="none",intfilename(nfiles_max),proftype="DOPPL",output="output"
  