SpectraRust/synspec/extracted/ghydop.f
2026-03-19 14:05:33 +08:00

51 lines
1.3 KiB
Fortran

subroutine ghydop(id,i0,i1,pj,absoh,emish)
c ==========================================
c
c hydrogen opacity -- lines + pseudocontinuum from Gomez tables
c
INCLUDE 'PARAMS.FOR'
INCLUDE 'MODELP.FOR'
INCLUDE 'SYNTHP.FOR'
COMMON/GOMOPA/frgtab(mfhtab),wlgtab(mfhtab),hydopg(mfhtab,mdepth),
* nugfreq
dimension absoh(mfreq),emish(mfreq),pj(40)
c
frg1=frgtab(1)
frg2=frgtab(nugfreq)
do 20 ij=i0,i1
fr=freq(ij)
if(fr.lt.frg1.or.fr.gt.frg2) go to 20
wla=2.997925e18/fr
frl=log10(fr)
c
if(ij.eq.i0) igf=nugfreq
10 continue
if(wla.gt.wlgtab(igf)) then
igf=igf-1
go to 10
end if
ig0=igf
if(ig0.le.2) ig0=2
ig1=igf-1
abl=(hydopg(ig1,id)-hydopg(ig0,id))*(wla-wlgtab(ig0))/
* (wlgtab(ig1)-wlgtab(ig0))+hydopg(ig0,id)
c
ii=1
if(freq(ij).gt.8.22013e14) then
pp=pj(1)*2.
else
pp=pj(2)*8.
end if
c
F15=FR*1.E-15
XKF=EXP(-4.79928e-11*FR/TEMP(ID))
XKFB=XKF*1.4743E-2*F15*F15*F15
oph=exp(abl)*pp
absoh(ij)=absoh(ij)+oph
emish(ij)=emish(ij)+oph*xkfb/(1.-xkf)
20 continue
c
return
end