SUBROUTINE GOMINI
C     =================
C
C     Initialization and reading of the opacity table for thermal processe
C     and Rayleigh scattering
c     raytab: scattering opacities in cm^2/gm at 5.0872638d14 Hz (sodium D)
c     (NOTE: Quantities in rayleigh.tab are in log_e)
C
c     tempvec: array of temperatures
c     rhovec: array of densities (gm/cm^3)
c     nu:     array of frequencies
c     table:  absorptive opacities in cm^2/gm
c     (NOTE:  Quantities in absorption.tab are in log_e)
C
      INCLUDE 'PARAMS.FOR'
      INCLUDE 'MODELP.FOR'
      COMMON/GOMOPA/frgtab(mfhtab),wlgtab(mfhtab),hydopg(mfhtab,mdepth),
     *              nugfreq
      common/gompar/hglim,ihgom
      dimension temvec(mtabth),elevec(mtabeh),
     *          hydcrs(mtabth,mtabeh,mfhtab)
c
      if(ihgom.eq.0) return
C
      open(53,file='gomhyd.dat',status='old')
c
      read(53,*) nugfreq,nugtemp,nugele
      read(53,*)
      read(53,*) (temvec(i),i=1,nugtemp)
      read(53,*)
      read(53,*) (elevec(j),j=1,nugele)
      do it=1,nugtemp
         temvec(it)=log(temvec(it)*1.161e4)
      end do
c     write(6,600) ihgom,nugfreq,nugtemp,nugele
c 600 format(' ihgom,nugfr,nugt,nuge ',4i4)
c
      EGTAB1 = elevec(1)
      EGTAB2 = elevec(nugele)
      TGTAB1 = temvec(1)
      TGTAB2 = temvec(nugtemp)
c
      do k = 1, nugfreq
         read(53,501) eneev
         frgtab(k)=3.28805e15/13.595*eneev
         wlgtab(k)=2.997925e18/frgtab(k)
         do i = 1, nugtemp
            read(53,*) (hydcrs(i,j,k),j=1,nugele)
         end do
      end do
      frg1=frgtab(1)
      frg2=frgtab(nugfreq)
c
  501 format(40x,f17.14)
      close(53)
C
c     Interpolate to the actual temperature and electron density
c     at the individual depth points
C
      do 10 id=1,nd
         if(elec(id).lt.HGLIM) go to 10
         rl=log(elec(id))
         tl=log(temp(id))
c
         DELTAR=(RL-EGTAB1)/(EGTAB2-EGTAB1)*FLOAT(nugele-1)
         JR = 1 + IDINT(DELTAR)
         IF(JR.LT.1) JR = 1
         IF(JR.GT.(nugele-1)) JR = nugele-1
         r1i=elevec(jr)
         r2i=elevec(jr+1)
         dri=(RL-R1i)/(R2i-R1i)
         if(JR .eq. 1) dri = 0.d0
C
         DELTAT=(TL-TGTAB1)/(TGTAB2-TGTAB1)*FLOAT(nugtemp-1)
         JP = 1 + IDINT(DELTAT)
         IF(JP.LT.1) JP = 1
         IF(JP.GT.nugtemp-1) JP = nugtemp-1
         t1i=temvec(jp)
         t2i=temvec(jp+1)
         dti=(TL-T1i)/(T2i-T1i)
         if(JP .eq. 1) dti = 0.d0
C
c        loop over tabular frequencies
c
         do jf=1,nugfreq
            opr1=hydcrs(jp,jr,jf)+dti*
     *           (hydcrs(jp+1,jr,jf)-hydcrs(jp,jr,jf))
            opr2=hydcrs(jp,jr+1,jf)+dti*
     *           (hydcrs(jp+1,jr+1,jf)-hydcrs(jp,jr+1,jf))
            opac=opr1+dri*(opr2-opr1)
            hydopg(jf,id)=opac+log(0.02654*4.1347e-15)
         end do
   10 continue
      return
      end