anchor.diff:            call readf(enermax)
             !
          case('GAUSSIAN','GAUSS','DOPPL','DOPPLER','RECT','BOX','BIN','STICKS','STICK','GAUS0','DOPP0',&
  -            'LOREN','LORENTZIAN','LORENTZ','MAX','VOIGT','HITRAN','PSEUDO','PSE-ROCCO','PSE-LIU')
  +            'LOREN','LORENTZIAN','LORENTZ','MAX','VOIGT','HITRAN','PSEUDO','PSE-ROCCO','PSE-LIU','VOI-QUAD')
             !
             proftype = trim(w)
             !
             if (trim(w(1:5))=="LOREN") ioffset = 500
  -          if (trim(w(1:5))=="VOIGT") ioffset = 500
  +          if (trim(w(1:))=="VOI") ioffset = 500
             if (trim(w(1:3))=="PSE") ioffset = 500
             offset = 25.0_rk
             !
  
positive.code:                !
                 tranfreq_i = freq(ipoint)
                 !
                 intens(ipoint)=intens(ipoint)+voigt_humlicek(tranfreq_i,tranfreq,halfwidth0,halfwidth)*abscoef
                 !
              enddo
  
negative.code:                intens(ipoint)=intens(ipoint)+abscoef*(intens1*eta1+intens2*eta2)
                 !
              enddo
              !$omp end parallel do             
              !
          case ('VOI-Q') ! VOIGT-QUADRATURES
              !
              alpha = halfwidth0
              gamma = halfwidth
              !
              sigma = alpha/ln22
              !
              y = gamma/sigma
              !
              dx2 = dfreq*0.5_rk/sigma
              !
              do iquad=1,nquad
                 !
                 xi = abciss(iquad)
                 !
                 x1 = (freq(ib)-tranfreq)/sigma-xi
                 x2 = (freq(ie)-tranfreq)/sigma-xi
                 !
                 L1 = atan( x1/y )
                 L2 = atan( x2/y )
                 !
                 bnormq(iquad) = 1.0_rk/(L2-L1)*pi
                 !
              enddo
              !
              Lorentz = 0
              !
              !$omp parallel do private(ipoint,xp,iquad,xi,x1,x2,L1,L2,Lorentz,voigt_) shared(intens) schedule(dynamic)
              do ipoint=ib,ie
                 !
                 xp = (freq(ipoint)-tranfreq)/sigma
                 !
                 do iquad=1,nquad
                    !
                    xi = abciss(iquad)
                    !
                    x1 = xp-dx2-xi
                    x2 = xp+dx2-xi
                    !
                    L1 = atan( x1/y )
                    L2 = atan( x2/y )
                    !
                    Lorentz(iquad) = (L2-L1)/pi
                    !
                 enddo
                 !
                 voigt_ = sum(Lorentz(1:nquad)*weight(1:nquad)*bnormq(1:nquad))/dfreq
                 !
                 intens(ipoint)=intens(ipoint)+abscoef*voigt_
                 !
              enddo
              !$omp end parallel do
              !
          case ('VOIGT')
              !
  