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