105 lines
2.8 KiB
Fortran
105 lines
2.8 KiB
Fortran
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
|