211 lines
5.8 KiB
Fortran
211 lines
5.8 KiB
Fortran
SUBROUTINE OPADD(MODE,ID,FR,ABAD,EMAD,SCAD)
|
|
C ===========================================
|
|
C
|
|
C Additional opacities
|
|
C This is basically user-supplied procedure; here are some more
|
|
C important non-standard opacity sources, namely
|
|
C Rayleigh scattering, H- opacity, H2+ opacity, and additional
|
|
C opacity of He I and He II.
|
|
C Inclusion of these opacities is contolled by switches transmitted
|
|
C by COMMON/OPCPAR - see description in START.
|
|
C
|
|
C Input parameters:
|
|
C MODE - controls the nature and the amount of calculations
|
|
C = -1 - (OPADD called from START) evaluation of relevant
|
|
C depth-dependent quantities (usually photoionization
|
|
C cross-sections, but also possibly other), which are
|
|
C stored in array CROS
|
|
C = 0 - evaluation of an additional opacity, emissivity, and
|
|
C scattering - for procedure OPAC0
|
|
C ID - depth index
|
|
C FR - frequency
|
|
C
|
|
C Output:
|
|
C
|
|
C ABAD - absorption coefficient (at frequency FR and depth ID)
|
|
C EMAD - emission coefficient (at frequency FR and depth ID)
|
|
C SCAD - scattering coefficient (at frequency FR and depth ID)
|
|
C
|
|
C
|
|
INCLUDE 'PARAMS.FOR'
|
|
INCLUDE 'MODELP.FOR'
|
|
PARAMETER (FRAYH = 2.463E15,
|
|
* FRAYHe = 5.150E15,
|
|
* FRAYH2 = 2.922E15,
|
|
* CLS = 2.997925e18)
|
|
C
|
|
AB0=0.
|
|
AB1=0.
|
|
ABAD=0.
|
|
EMAD=0.
|
|
SCAD=0.
|
|
C
|
|
if(iath.gt.0) then
|
|
N0HN=NFIRST(IELH)
|
|
NKH=NKA(IATH)
|
|
C
|
|
IF(MODE.GE.0) THEN
|
|
T=TEMP(ID)
|
|
ANE=ELEC(ID)
|
|
HKT=HK/T
|
|
T32=1./T/SQRT(T)
|
|
END IF
|
|
anh=dens(id)/(wmm(id)*ytot(id))
|
|
anhe=rrr(id,1,2)
|
|
C
|
|
IT=NLEVEL
|
|
C
|
|
C -----------------------
|
|
C HI Rayleigh scattering
|
|
C -----------------------
|
|
C
|
|
IF(IRSCT.NE.0.AND.IOPHLI.NE.1.AND.IOPHLI.NE.2) THEN
|
|
X=1.D0/(CLS/MIN(FR,FRAYH))**2
|
|
SG=(5.799E-13+(1.422E-6+2.784*X)*X)*X*X
|
|
c ABAD=POPUL(N0HN,ID)*SG
|
|
SCAD=POPUL(N0HN,ID)*SG
|
|
scad=anh*sg
|
|
END IF
|
|
IF(IOPHMI.NE.0) THEN
|
|
C
|
|
C ----------------------------
|
|
C H- bound-free and free-free
|
|
C ----------------------------
|
|
C Note: IOPHMI must not by taken non-zero if H- is considered
|
|
C explicitly, because H- opacity would be taken twice
|
|
C
|
|
SG=SBFHMI(FR)
|
|
XHM=8762.9/T
|
|
SB=1.0353E-16*T32*EXP(XHM)*POPUL(N0HN,ID)*ANE*SG
|
|
SF=SFFHMI(POPUL(N0HN,ID),FR,T)*ANE
|
|
AB0=SB+SF
|
|
END IF
|
|
C
|
|
C -----------------------
|
|
C He I Rayleigh scattering
|
|
C -----------------------
|
|
C
|
|
IF(IRSCHE.NE.0.AND.MODE.GE.0) THEN
|
|
X=(CLS/MIN(FR,FRAYHe))**2
|
|
CS=5.484E-14/X/X*(1.+(2.44E5+5.94E10/(X-2.90E5))/X)**2
|
|
sg=anhe*cs
|
|
c abad=abad+sg
|
|
scad=scad+sg
|
|
END IF
|
|
C
|
|
C -----------------------
|
|
C H2 Rayleigh scattering
|
|
C -----------------------
|
|
C
|
|
IF(IRSCH2.NE.0.AND.MODE.GE.0.AND.IFMOL.GT.0) THEN
|
|
X=(CLS/MIN(FR,FRAYH2))**2
|
|
X2=1./X/X
|
|
CS=(8.14E-13+1.28E-6/X+1.61*X2)*X2
|
|
sg=cs*anh2(id)
|
|
c abad=abad+sg
|
|
scad=scad+sg
|
|
END IF
|
|
C
|
|
IF(IOPH2P.GT.0.AND.IFMOL.GT.0.and.
|
|
* t.lt.tmolim.and.fr.lt.3.28e15) THEN
|
|
C
|
|
C -----------------------------
|
|
C H2+ bound-free and free-free
|
|
C -----------------------------
|
|
C
|
|
X=FR*1.E-15
|
|
SG1=(-7.342E-3+(-2.409+(1.028+(-4.23E-1+
|
|
* (1.224E-1-1.351E-2*X)*X)*X)*X)*X)*1.602E-12/BOLK
|
|
IT=IT+1
|
|
X=LOG(FR)
|
|
SG2=-3.0233E3+(3.7797E2+(-1.82496E1+(3.9207E-1-
|
|
* 3.1672E-3*X)*X)*X)*X
|
|
X2=-SG1/T+SG2
|
|
SB=0.
|
|
IF(X2.GT.-150.) SB=POPUL(N0HN,ID)*POPUL(NKH,ID)*EXP(X2)
|
|
AB0=AB0+SB
|
|
END IF
|
|
end if
|
|
C
|
|
C -----------------------------
|
|
C He- free-free
|
|
C -----------------------------
|
|
C
|
|
if(mode.ge.0.and.iophem.gt.0) then
|
|
A=3.397D-46+(-5.216D-31+7.039D-15/FR)/FR
|
|
B=-4.116D-42+(1.067D-26+8.135D-11/FR)/FR
|
|
C=5.081D-37+(-8.724D-23-5.659D-8/FR)/FR
|
|
cs=a*t+b+c/t
|
|
sg=anhe*ane*cs
|
|
ab0=ab0+sg
|
|
end if
|
|
C
|
|
C -----------------------------
|
|
C H2- free-free
|
|
C -----------------------------
|
|
C
|
|
IF(IOPH2M.NE.0.AND.MODE.GE.0.AND.IFMOL.GT.0.AND.T.LT.TMOLIM) THEN
|
|
call h2minus(t,anh2(id),ane,fr,oph2)
|
|
ab1=ab1+oph2
|
|
END IF
|
|
C
|
|
C -----------------------------
|
|
C CH and OH continuuum opacity
|
|
C -----------------------------
|
|
C
|
|
if(mode.ge.0.and.ifmol.gt.0.and.t.lt.tmolim) then
|
|
if(iopch.gt.0) ab0=ab0+sbfch(fr,t)*anch(id)
|
|
if(iopoh.gt.0) ab0=ab0+sbfoh(fr,t)*anoh(id)
|
|
C
|
|
C ---------------------------
|
|
C CIA H2-H2 opacity
|
|
C ---------------------------
|
|
C
|
|
if(ioh2h2.gt.0) then
|
|
call cia_h2h2(t,anh2(id),fr,oph2)
|
|
ab1=ab1+oph2
|
|
end if
|
|
C
|
|
C ---------------------------
|
|
C CIA H2-He opacity
|
|
C ---------------------------
|
|
C
|
|
if(ioh2he.gt.0) then
|
|
call cia_h2he(t,anh2(id),anhe,fr,oph2)
|
|
ab1=ab1+oph2
|
|
end if
|
|
C
|
|
C ---------------------------
|
|
C CIA H2-H opacity
|
|
C ---------------------------
|
|
C
|
|
if(ioh2h1.gt.0) then
|
|
call cia_h2h(t,anh2(id),anh,fr,oph2)
|
|
ab1=ab1+oph2
|
|
end if
|
|
C
|
|
C ---------------------------
|
|
C CIA H-He opacity
|
|
C ---------------------------
|
|
C
|
|
if(iohhe.gt.0) then
|
|
call cia_hhe(t,anh,anhe,fr,oph2)
|
|
ab1=ab1+oph2
|
|
end if
|
|
end if
|
|
C
|
|
C ----------------------------------------------
|
|
C The user may supply more opacity sources here:
|
|
C ----------------------------------------------
|
|
C
|
|
C Finally, actual absorption and emission coefficients
|
|
|
|
IF(MODE.LT.0) RETURN
|
|
X=EXP(-HKT*FR)
|
|
X1=1.-X
|
|
BNX=BN*(FR*1.E-15)**3*X
|
|
ABAD=ABAD+X1*AB0+AB1
|
|
EMAD=EMAD+BNX*(AB0+AB1/X1)
|
|
RETURN
|
|
END
|