359 lines
10 KiB
Fortran
359 lines
10 KiB
Fortran
SUBROUTINE OPACF0(ID,NFRQ)
|
|
C ==========================
|
|
C
|
|
C Absorption, emission, and scattering coefficients
|
|
C at depth ID
|
|
C
|
|
C Input: ID opacity and emissivity is calculated for the
|
|
C depth point ID
|
|
C Output: ABSO - array of absorption coefficient
|
|
C EMIS - array of emission coefficient
|
|
C SCAT - array of scattering coefficient
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ODFPAR.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
PARAMETER (FRH=3.28805E15, PH2=2.815D29*2., EHB=157802.77355)
|
|
PARAMETER (CFF1=1.3727D-25,CFF2=4.3748D-10,CFF3=2.5993D-7)
|
|
PARAMETER (C14=2.99793D14)
|
|
PARAMETER (SGFF0 = 3.694D8)
|
|
common/hmolab/anh2(mdepth),anhm(mdepth)
|
|
DIMENSION FREDG(NLMX),S(NLMX),SUM(NLMX),PRF(MFREQL)
|
|
C
|
|
C initialize
|
|
c
|
|
C this part is analogous to TDPINI - for one depth only
|
|
C
|
|
T=TEMP(ID)
|
|
T1=UN/T
|
|
HKT1(ID)=HK*T1
|
|
HKT21(ID)=HKT1(ID)*T1
|
|
TK1(ID)=HKT1(ID)/H
|
|
SQT1(ID)=SQRT(T)
|
|
TEMP1(ID)=T1
|
|
CALL GFREE0(ID)
|
|
EMEL1(ID)=UN
|
|
LASER=ITER.GT.ITLAS
|
|
c
|
|
C this part is analogous to OPAINI - for one depth only
|
|
C
|
|
ANE=ELEC(ID)
|
|
ELEC1(ID)=UN/ANE
|
|
DENS1(ID)=UN/DENS(ID)
|
|
DENSI(ID)=DENS1(ID)
|
|
DENSIM(ID)=DENSI(ID)*WMM(ID)
|
|
ELSCAT(ID)=ANE*SIGE
|
|
if(izscal.eq.1) then
|
|
densi(id)=un
|
|
densim(id)=0.
|
|
end if
|
|
CALL DWNFR0(ID)
|
|
CALL WNSTOR(ID)
|
|
CALL SABOLF(ID)
|
|
c
|
|
c quantities for the bound-free opacity
|
|
c
|
|
DO IBFT=1,NTRANC
|
|
ITR=ITRBF(IBFT)
|
|
IF(INDEXP(ITR).NE.0) THEN
|
|
II=ILOW(ITR)
|
|
JJ=IUP(ITR)
|
|
IT=ITRA(JJ,II)
|
|
IE=IEL(II)
|
|
NKE=NNEXT(IE)
|
|
CORR=UN
|
|
IF(NKE.NE.JJ) CORR=G(NKE)/G(JJ)*
|
|
* EXP((ENION(NKE)-ENION(JJ))*TK1(ID))
|
|
ABTRA(ITR,ID)=POPUL(II,ID)
|
|
EMTRA(ITR,ID)=POPUL(JJ,ID)*ANE*SBF(II)*WOP(II,ID)*CORR
|
|
END IF
|
|
END DO
|
|
c
|
|
c quantities for the free-free opacity
|
|
c
|
|
IF(IELHM.GT.0) THEN
|
|
CFFN(ID)=POPUL(NFIRST(IELH),ID)*ANE
|
|
CFFT(ID)=CFF2-CFF3/T
|
|
END IF
|
|
SGFF=SGFF0/SQT1(ID)*ANE
|
|
DO ION=1,NION
|
|
SFF2(ION,ID)=EXP(FF(ION)*HKT1(ID))
|
|
SFF3(ION,ID)=POPUL(NNEXT(ION),ID)*CHARG2(ION)*SGFF
|
|
END DO
|
|
c
|
|
C this part is analogous to SGMER0 - for one depth only
|
|
C
|
|
IMER=0
|
|
DO II=1,NLEVEL
|
|
IF(IFWOP(II).LT.0) THEN
|
|
IMER=IMER+1
|
|
IMRG(II)=IMER
|
|
IIMER(IMER)=II
|
|
IE=IEL(II)
|
|
CH=IZ(IE)*IZ(IE)
|
|
FRCH(IMER)=FRH*CH
|
|
SGM0(IMER)=PH2*CH*CH
|
|
II0=NQUANT(II-1)+1
|
|
EX=EHB*CH*TEMP1(ID)
|
|
DO I=II0,NLMX
|
|
FREDG(I)=FRCH(IMER)*XI2(I)
|
|
EXI=EXP(EX*XI2(I))
|
|
S(I)=EXI*WNHINT(I,ID)*XI3(I)
|
|
SUM(I)=0.
|
|
END DO
|
|
SUM(NLMX)=S(NLMX)
|
|
DO I=NLMX-1,II0,-1
|
|
SUM(I)=SUM(I+1)+S(I)
|
|
END DO
|
|
DO I=1,II0-1
|
|
SUM(I)=SUM(II0)
|
|
END DO
|
|
SGEM=SGM0(IMER)/GMER(IMER,ID)
|
|
DO I=1,NLMX
|
|
SGMSUM(I,IMER,ID)=SUM(I)*SGEM
|
|
END DO
|
|
END IF
|
|
END DO
|
|
C
|
|
C initialization of the line opacity
|
|
C
|
|
IF(NFRQ.GT.NFREQC) THEN
|
|
DO 10 ITR=1,NTRANS
|
|
IF(.NOT.LINE(ITR)) GO TO 10
|
|
IF(INTMOD(ITR).EQ.0) GO TO 10
|
|
INDXA=IABS(INDEXP(ITR))
|
|
IJL0=IFR0(ITR)
|
|
IJL1=IFR1(ITR)
|
|
IF(ISPODF.GE.1) THEN
|
|
IJL0=KFR0(ITR)
|
|
IJL1=KFR1(ITR)
|
|
END IF
|
|
II=ILOW(ITR)
|
|
JJ=IUP(ITR)
|
|
IF(INDXA.LT.2.OR.INDXA.GT.4) THEN
|
|
CALL LINPRO(ITR,ID,PRF)
|
|
DO IJ=IJL0,IJL1
|
|
PRFLIN(ID,IJ)=real(PRF(IJ-IJL0+1))
|
|
END DO
|
|
END IF
|
|
10 CONTINUE
|
|
GG=G(II)/G(JJ)
|
|
IF(IFWOP(JJ).GE.0) THEN
|
|
PI=POPUL(II,ID)*WOP(JJ,ID)
|
|
PJ=POPUL(JJ,ID)*WOP(II,ID)*GG
|
|
ELSE
|
|
PI=POPUL(II,ID)
|
|
PJ=POPUL(JJ,ID)*WOP(II,ID)*G(II)/GMER(IMRG(JJ),ID)
|
|
END IF
|
|
ABTRA(ITR,ID)=PI
|
|
EMTRA(ITR,ID)=PJ*EXP(FR0(ITR)*HKT1(ID))
|
|
IF(LASER) THEN
|
|
qtt=0.
|
|
if(pi.ne.pj) QTT=PJ/(PI-PJ)*(EXP(FR0(ITR)*HKT1(ID))-UN)
|
|
lfr=fr0(itr).lt.frtabm.and.iadop(iatm(ii)).gt.0
|
|
IF(QTT.LT.0. .OR. QTT.GT.QTLAS .or. lfr) THEN
|
|
ABTRA(ITR,ID)=0.
|
|
EMTRA(ITR,ID)=0.
|
|
END IF
|
|
END IF
|
|
END IF
|
|
C
|
|
C ---------------------------------------------------------
|
|
C
|
|
C loop over frequency points
|
|
C
|
|
ICALL=1
|
|
DO IJ=1,NFRQ
|
|
if(icompt.gt.0) ELSCAT(ID)=ELEC(ID)*SIGEC(IJ)
|
|
ABSO(IJ)=ELSCAT(ID)
|
|
EMIS(IJ)=0.
|
|
SCAT(IJ)=ELSCAT(ID)
|
|
C
|
|
C basic frequency- and depth-dependent quantities
|
|
C
|
|
FR=FREQ(IJ)
|
|
FRINV=UN/FR
|
|
FR3INV=FRINV*FRINV*FRINV
|
|
XKF(ID)=EXP(-HKT1(ID)*FR)
|
|
XKF1(ID)=UN-XKF(ID)
|
|
XKFB(ID)=XKF(ID)*BNUE(IJ)
|
|
C
|
|
C ******** 1. bound-free contribution
|
|
C
|
|
DO 30 IBFT=1,NTRANC
|
|
ITR=ITRBF(IBFT)
|
|
II=ILOW(ITR)
|
|
JJ=IUP(ITR)
|
|
if(iadop(iatm(ii)).gt.0.and.fr.le.frtabm) go to 30
|
|
if(ifdiel.eq.0) then
|
|
SG=CROSS(IBFT,IJ)
|
|
else
|
|
SG=CROSSD(IBFT,IJ,ID)
|
|
endif
|
|
IF(IFWOP(II).LT.0) THEN
|
|
IMER=IMRG(II)
|
|
CALL SGMER1(FRINV,FR3INV,IMER,ID,SGME1)
|
|
SGMG(IMER,ID)=SGME1
|
|
SG=SGME1
|
|
END IF
|
|
IF(SG.LE.0.) GO TO 30
|
|
IF(MCDW(ITR).GT.0) THEN
|
|
IZZ=IZ(IEL(II))
|
|
CALL DWNFR1(FR,FR0(ITR),ID,IZZ,DW1)
|
|
DWF1(MCDW(ITR),ID)=DW1
|
|
SG=SG*DW1
|
|
END IF
|
|
EMISBF=SG*EMTRA(ITR,ID)
|
|
ABSO(IJ)=ABSO(IJ)+SG*ABTRA(ITR,ID)
|
|
EMIS(IJ)=EMIS(IJ)+EMISBF
|
|
30 CONTINUE
|
|
C
|
|
C ******** 2. free-free contribution
|
|
C
|
|
DO 40 ION=1,NION
|
|
IT=ITRA(NNEXT(ION),NNEXT(ION))
|
|
if(iadop(iatm(nnext(ioon))).gt.0.and.fr.le.frtabm) go to 40
|
|
C
|
|
C hydrogenic with Gaunt factor = 1
|
|
C
|
|
IF(IT.EQ.1) THEN
|
|
SF1=SFF3(ION,ID)*FR3INV
|
|
SF2=SFF2(ION,ID)
|
|
IF(FR.LT.FF(ION)) SF2=UN/XKF(ID)
|
|
ABSOFF=SF1*SF2
|
|
ABSO(IJ)=ABSO(IJ)+ABSOFF
|
|
EMIS(IJ)=EMIS(IJ)+ABSOFF
|
|
C
|
|
C hydrogenic with exact Gaunt factor
|
|
C
|
|
ELSE IF(IT.EQ.2) THEN
|
|
SF1=SFF3(ION,ID)*FR3INV
|
|
SF2=SFF2(ION,ID)
|
|
IF(FR.LT.FF(ION)) SF2=UN/XKF(ID)
|
|
X=C14*CHARG2(ION)/FR
|
|
SF2=SF2-UN+GFREE1(ID,X)
|
|
ABSOFF=SF1*SF2
|
|
ABSO(IJ)=ABSO(IJ)+ABSOFF
|
|
EMIS(IJ)=EMIS(IJ)+ABSOFF
|
|
C
|
|
C H minus free-free opacity
|
|
C
|
|
ELSE IF(IT.EQ.3) THEN
|
|
c ABSOFF=(CFF1+CFFT(ID)*FRINV)*CFFN(ID)*FRINV
|
|
ABSOFF=SFFHMI(POPUL(NFIRST(IELH),ID),FR,TEMP(ID))*
|
|
* ELEC(ID)
|
|
ABSO(IJ)=ABSO(IJ)+ABSOFF
|
|
EMIS(IJ)=EMIS(IJ)+ABSOFF
|
|
C
|
|
C special evaluation of the cross-section
|
|
C
|
|
ELSE IF(IT.LT.0) THEN
|
|
ABSOFF=FFCROS(ION,IT,TEMP(ID),FR)*
|
|
* POPUL(NNEXT(ION),ID)*ELEC(ID)
|
|
ABSO(IJ)=ABSO(IJ)+ABSOFF
|
|
EMIS(IJ)=EMIS(IJ)+ABSOFF
|
|
END IF
|
|
40 CONTINUE
|
|
C
|
|
C ******** 3. - additional continuum opacity (OPADD)
|
|
C
|
|
IF(IOPADD.NE.0) THEN
|
|
CALL OPADD(0,ICALL,IJ,ID)
|
|
ABSO(IJ)=ABSO(IJ)+ABAD
|
|
EMIS(IJ)=EMIS(IJ)+EMAD
|
|
SCAT(IJ)=SCAT(IJ)+SCAD
|
|
END IF
|
|
C
|
|
C ******** 4. - opacity and emissivity in lines
|
|
C
|
|
IF(ISPODF.EQ.0) THEN
|
|
IF(IJLIN(IJ).GT.0) THEN
|
|
C
|
|
C the "primary" line at the given frequency
|
|
C
|
|
ITR=IJLIN(IJ)
|
|
iad=iadop(iatm(ilow(itr)))
|
|
if(iad.eq.0.or.(lfre.and.iad.gt.0)) then
|
|
SG=PRFLIN(ID,IJ)
|
|
ABSO(IJ)=ABSO(IJ)+SG*ABTRA(ITR,ID)
|
|
EMIS(IJ)=EMIS(IJ)+SG*EMTRA(ITR,ID)
|
|
end if
|
|
END IF
|
|
IF(NLINES(IJ).LE.0) GO TO 110
|
|
C
|
|
C the "overlapping" lines at the given frequency
|
|
C
|
|
DO 100 ILINT=1,NLINES(IJ)
|
|
ITR=ITRLIN(ILINT,IJ)
|
|
iad=iadop(iatm(ilow(itr)))
|
|
if(iad.gt.0.and..not.lfre) go to 100
|
|
if(linexp(itr)) go to 100
|
|
IJ0=IFR0(ITR)
|
|
DO IJT=IJ0,IFR1(ITR)
|
|
IF(FREQ(IJT).LE.FR) THEN
|
|
IJ0=IJT
|
|
GO TO 70
|
|
END IF
|
|
END DO
|
|
70 IJ1=IJ0-1
|
|
X=UN/(FREQ(IJ1)-FREQ(IJ0))
|
|
A1=(FR-FREQ(IJ0))*X
|
|
X=UN/(FREQ(IJ1)-FREQ(IJ0))
|
|
A2=(FREQ(IJ1)-FR)*X
|
|
SG=A1*PRFLIN(ID,IJ1)+A2*PRFLIN(ID,IJ0)
|
|
ABSO(IJ)=ABSO(IJ)+SG*ABTRA(ITR,ID)
|
|
EMIS(IJ)=EMIS(IJ)+SG*EMTRA(ITR,ID)
|
|
100 CONTINUE
|
|
110 CONTINUE
|
|
C
|
|
C Opacity sampling option
|
|
C
|
|
ELSE
|
|
IF(NLINES(IJ).LE.0) GO TO 400
|
|
DO 300 ILINT=1,NLINES(IJ)
|
|
ITR=ITRLIN(ILINT,IJ)
|
|
iad=iadop(iatm(ilow(itr)))
|
|
if(iad.gt.0.and..not.lfre) go to 300
|
|
KJ=IJ-IFR0(ITR)+KFR0(ITR)
|
|
INDXPA=IABS(INDEXP(ITR))
|
|
IF(INDXPA.NE.3 .AND. INDXPA.NE.4) THEN
|
|
DO ID=1,ND
|
|
SG=PRFLIN(ID,KJ)
|
|
ABSO(IJ)=ABSO(IJ)+SG*ABTRA(ITR,ID)
|
|
EMIS(IJ)=EMIS(IJ)+SG*EMTRA(ITR,ID)
|
|
END DO
|
|
ELSE
|
|
DO ID=1,ND
|
|
KJD=JIDI(ID)
|
|
SG=EXP(XJID(ID)*SIGFE(KJD,KJ)+(UN-XJID(ID))*
|
|
* SIGFE(KJD+1,KJ))
|
|
ABSO(IJ)=ABSO(IJ)+SG*ABTRA(ITR,ID)
|
|
EMIS(IJ)=EMIS(IJ)+SG*EMTRA(ITR,ID)
|
|
END DO
|
|
END IF
|
|
300 CONTINUE
|
|
400 CONTINUE
|
|
END IF
|
|
C
|
|
C ----------------------------
|
|
C total opacity and emissivity
|
|
C ----------------------------
|
|
C
|
|
ABSO(IJ)=ABSO(IJ)-EMIS(IJ)*XKF(ID)
|
|
EMIS(IJ)=EMIS(IJ)*XKFB(ID)
|
|
c
|
|
c --------------------------------------------------------
|
|
c contribution from precalculated background opacity table
|
|
c --------------------------------------------------------`
|
|
c
|
|
if(ioptab.gt.0) then
|
|
call opact1(ij)
|
|
end if
|
|
c
|
|
END DO
|
|
RETURN
|
|
END
|