SUBROUTINE OPCTAB(FR,IJ,ID,T,RHO,AB,SC,SCT,IGRAM)
C     =================================================
C
C     opacity for a given temperature and density computed 
C     by an interpolation of the precalculated opacity table
C
C     This is a simplified routine with all interpolations linear
C
C     Input:   FR  - frequency   (Hz)
C              T   - temperature (K)
C              RHO - density     (g cm^-3)
C      Outout: AB  - absorptive opacity  (per gram)
C              SC  - scattering opacity  (per gram) 
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      parameter (frray0 = 5.0872638d14)
C
      jf=ij
      frij=fr
c
      if(numtemp.eq.nd) then
         opac=absopac(id,1,jf)
         go to 10
      end if
c
      TL=LOG(T)
      DELTAT=(TL-TTAB1)/(TTAB2-TTAB1)*FLOAT(numtemp-1)
      JT = 1 + IDINT(DELTAT)
      JU = JT + 1
      IF(JT.LT.1) JT = 1
      IF(JT.GT.numtemp-1) JT = numtemp-1
      t1i=tempvec(jt)
      t2i=tempvec(jt+1)
      dti=(TL-T1i)/(T2i-T1i)
      if(deltat.lt.0.) dti = 0.d0
C
      if(numrh(1).ne.1) then
c
c        lower temperature
c
         numrho=numrh(jt)
         rtab1=rhomat(jt,1)
         rtab2=rhomat(jt,numrho)
         RL = LOG(RHO) 
         DELTAR=(RL-RTAB1)/(RTAB2-RTAB1)*FLOAT(numrho-1)
         JR = 1 + IDINT(DELTAR)
         IF(JR.LT.1) JR = 1
         IF(JR.GT.(numrho-1)) JR = numrho-1
         r1i=rhomat(jt,jr)
         r2i=rhomat(jt,jr+1)
         dri=(RL-R1i)/(R2i-R1i)
         if(deltar.lt.0.) dri = 0.d0
         opr1=absopac(jt,jr,jf)+
     *        dri*(absopac(jt,jr+1,jf)-absopac(jt,jr,jf))
c
c        higher temperature
c        
         ju=jt+1
         numrho=numrh(ju)
         rtab1=rhomat(ju,1)
         rtab2=rhomat(ju,numrho)
         RL = LOG(RHO)
         DELTAR=(RL-RTAB1)/(RTAB2-RTAB1)*FLOAT(numrho-1)
         JR = 1 + IDINT(DELTAR)
         IF(JR.LT.1) JR = 1
         IF(JR.GT.(numrho-1)) JR = numrho-1
         r1i=rhomat(ju,jr)
         r2i=rhomat(ju,jr+1)
         dri=(RL-R1i)/(R2i-R1i)
         if(deltar.lt.0.) dri = 0.d0
         opr2=absopac(ju,jr,jf)+
     *        dri*(absopac(ju,jr+1,jf)-absopac(ju,jr,jf))
         opac=opr1+(opr2-opr1)*dti
       else
         jr=1
         opac=absopac(jt,jr,jf)+(absopac(ju,jr,jf)-
     *        absopac(jt,jr,jf))*dti
      end if
   10 continue
      opac=exp(opac)
C
      AB=opac
C
C     ************************************************************
C     scattering
C     ************************************************************
C     1. Rayleigh scattering
C
      sct=0.
      if(ifrayl.lt.0) then
         sct=raysc(id)*(freq(jf)/frray0)**4
       else if(ifrayl.gt.0) then
         call rayleigh(1,ij,id,scr)
         sct=scr/dens(id)
      end if
      if(ioptab.lt.0) sct=sct+sige*elec(id)/rho     
c     sct=sc+sige*elec(id)/dens(id)
c     sct=sc
C
c     2. cloud scattering (not yet implemented)
c
      if(iter.le.0) return
c
c     ab=ab+abscld(id,jf)
c     sct=sct+scacld(id,jf,1)
c
      if(igram.eq.0) then
         ab = ab*rho
         sc = sc*rho
         sct=sct*rho
      end if
c
      RETURN
      END