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     nu:     array of frequencies
c     table:  absorptive opacities in cm^2/gm
c     (NOTE:  Quantities in absorption.tab are in log_e)
C
      INCLUDE 'IMPLIC.FOR'
      INCLUDE 'BASICS.FOR'
      INCLUDE 'MODELQ.FOR'
      common/intcfg/yint(mfreq),jgint(mfreq)
c
      dimension absort(mfhtab),frlt(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
         frlt(k)=log10(frgtab(k))
         do i = 1, nugtemp
            read(53,*) (hydcrs(i,j,k),j=1,nugele)
         end do
      end do
  501 format(40x,f17.14)
      close(53)
c
      FRGTB1 = log10(frgtab(1))
      FRGTB2 = log10(frgtab(nugfreq))
c
c     interpolate HYDCRS from tabular to actual frequencies
c
c     set up interpolation coefficients for frequency interpolation
c     by bisection
c
      frg1=frgtab(1)
      frg2=frgtab(nugfreq)
      do ij=1,nfreq
         jgint(ij)=0
         xint=freq(ij)
         if(xint.ge.frg1.and.xint.le.frg2) then
            jl=0
            ju=nugfreq+1
   40       continue           
            if(ju-jl.gt.1) then
               jm=(ju+jl)/2
               if((frg2.gt.frg1).eqv.(xint.gt.frgtab(jm))) then
                  jl=jm
                else
                  ju=jm
               end if
               go to 40
            end if
            j=jl
            if(j.eq.nugfreq) j=j-1
            if(j.eq.0) j=j+1
            jgint(ij)=j
            yint(ij)=un/log10(frgtab(j+1)/frgtab(j))
         end if
      end do
c      
      do it=1,nugtemp
         do ir=1,nugele
            do k=1,nugfreq
               absort(k)=log(hydcrs(it,ir,k))
            end do
            do ij=1,nfreq
               j=jgint(ij)
               hydcrs(it,ir,ij)=0.
               if(j.gt.0) then
                  rc=(absort(j+1)-absort(j))*yint(ij)
                  hcs=rc*log10(freq(ij)/frgtab(j))+absort(j)
                  hydcrs(it,ir,ij)=hcs
               end if        
            end do
         end do
      end do
c
      RETURN
      END