117 lines
3.0 KiB
Fortran
117 lines
3.0 KiB
Fortran
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
|