SUBROUTINE OPACF1(IJ) C ===================== C C Absorption, emission, and scattering coefficients C at frequency IJ and for all depths C C Input: IJ opacity and emissivity is calculated for the C frequency points with index IJ C Output: ABSO1 - array of absorption coefficient C EMIS1 - array of emission coefficient C SCAT1 - array of scattering coefficient (all scattering C mechanisms except electron scattering) C INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' INCLUDE 'ATOMIC.FOR' INCLUDE 'MODELQ.FOR' INCLUDE 'ODFPAR.FOR' INCLUDE 'ALIPAR.FOR' common/hmolab/anh2(mdepth),anhm(mdepth) common/ipricr/iprcrs,nprcrs PARAMETER (C14=2.99793D14, c10=c14*1.d-4,CFF1=1.3727D-25) dimension pold(mlevel),abtrh(mtrans) c if(ioptab.lt.0) then do id=1,nd abso1(id)=0. scat1(id)=0. emis1(id)=0 absot(id)=0. end do call opact1(ij) return end if C C initialize c elscat(id)=elec(id)*sige IF(ICOMPT.GT.0) THEN DO ID=1,ND ELSCAT(ID)=ELEC(ID)*SIGEC(IJ) END DO END IF C DO ID=1,ND c ABSO1(ID)=ELSCAT(ID) ABSO1(ID)=0. EMIS1(ID)=0. SCAT1(ID)=ELSCAT(ID) END DO C C basic frequency- and depth-dependent quantities C FR=FREQ(IJ) FRINV=UN/FR FR3INV=FRINV*FRINV*FRINV DO ID=1,ND XKF(ID)=EXP(-HKT1(ID)*FR) XKF1(ID)=UN-XKF(ID) XKFB(ID)=XKF(ID)*BNUE(IJ) END DO if(ielh.gt.0) n0hn=nfirst(ielh) al=2.997925e18/fr lpri=al.gt.1579.0.and.al.lt.1579.5 c lfre=freq(ij).gt.frtabm if(iprcrs.gt.0) then abso1(iprcrs)=0. do ii=nfirst(ielh),nlast(ielh) if(ii.ne.nprcrs+nfirst(ielh)-1) then pold(ii)=popul(ii,iprcrs) popul(ii,iprcrs)=0. do jj=ii+1,nnext(ielh) itrh=itra(ii,jj) abtrh(itrh)=abtra(itrh,iprcrs) abtra(itrh,iprcrs)=0. end do end if end do end if C C ******** 1a. bound-free contribution - without dielectronic rec. C if(ifdiel.eq.0) then DO IBFT=1,NTRANC ITR=ITRBF(IBFT) II=ILOW(ITR) iad=iadop(iatm(ii)) lcomop=iad.eq.0.or.(lfre.and.iad.gt.0) SG=CROSS(IBFT,IJ) IF(SG.GT.0..and.lcomop) THEN JJ=IUP(ITR) IZZ=IZ(IEL(II)) IMER=IMRG(II) DO ID=1,ND SGD=SG IF(MCDW(ITR).GT.0) THEN CALL DWNFR1(FR,FR0(ITR),ID,IZZ,DW1) DWF1(MCDW(ITR),ID)=DW1 SGD=SG*DW1 END IF IF(IFWOP(II).LT.0) THEN CALL SGMER1(FRINV,FR3INV,IMER,ID,SGME1) SGMG(IMER,ID)=SGME1 SGD=SGME1 END IF EMISBF=SGD*EMTRA(ITR,ID) ABSO1(ID)=ABSO1(ID)+SGD*ABTRA(ITR,ID) EMIS1(ID)=EMIS1(ID)+EMISBF c if(lpri.and.id.eq.40) c * write(65,621) id,ij,al,itr,ii,jj,typion(iel(ii)), c * ii-nfirst(iel(ii))+1,jj-nfirst(iel(ii))+1, c * popul(ii,id),sgd*abtra(itr,id),abso1(id),emisbf,emis1(id) c end if c 621 format('bf',i4,i6,f10.3,3i5,2x,a4,2x,2i4,1p5e14.7) END DO END IF END DO else C C ******** 1b. bound-free contribution - with dielectronic rec. C DO IBFT=1,NTRANC ITR=ITRBF(IBFT) II=ILOW(ITR) iad=iadop(iatm(ii)) lcomop=iad.eq.0.or.(lfre.and.iadop(iatm(ii)).gt.0) SG=CROSS(IBFT,IJ) IF(SG.GT.0..and.lcomop) THEN JJ=IUP(ITR) IZZ=IZ(IEL(II)) IMER=IMRG(II) DO ID=1,ND SG=CROSSD(IBFT,IJ,ID) IF(SG.GT.0.) THEN SGD=SG IF(MCDW(ITR).GT.0) THEN CALL DWNFR1(FR,FR0(ITR),ID,IZZ,DW1) DWF1(MCDW(ITR),ID)=DW1 SGD=SG*DW1 END IF IF(IFWOP(II).LT.0) THEN CALL SGMER1(FRINV,FR3INV,IMER,ID,SGME1) SGMG(IMER,ID)=SGME1 SGD=SGME1 END IF EMISBF=SGD*EMTRA(ITR,ID) ABSO1(ID)=ABSO1(ID)+SGD*ABTRA(ITR,ID) EMIS1(ID)=EMIS1(ID)+EMISBF END IF END DO END IF END DO end if C C ******** 2. free-free contribution C DO 40 ION=1,NION IT=ITRA(NNEXT(ION),NNEXT(ION)) iad=iadop(iatm(nnext(ion))) if(iad.gt.0.and..not.lfre) go to 40 C C hydrogenic with Gaunt factor = 1 C IF(IT.EQ.1) THEN DO ID=1,ND SF1=SFF3(ION,ID)*FR3INV SF2=SFF2(ION,ID) IF(FR.LT.FF(ION)) SF2=UN/XKF(ID) ABSOFF=SF1*SF2 ABSO1(ID)=ABSO1(ID)+ABSOFF EMIS1(ID)=EMIS1(ID)+ABSOFF c if(lpri.and.mod(id,20).eq.1) then c write(6,622) id,ij,ion,typion(ion),sf1,sf2,absoff,abso1(id) c end if c 622 format('ff',i4,i6,i4,2x,a4,2x,1p4e14.6) END DO C C hydrogenic with exact Gaunt factor C ELSE IF(IT.EQ.2) THEN DO ID=1,ND 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 ABSO1(ID)=ABSO1(ID)+ABSOFF EMIS1(ID)=EMIS1(ID)+ABSOFF c if(lpri.and.mod(id,20).eq.1) then c sgf=gfree1(id,x) c write(6,624) id,ij,ion,typion(ion),ff(ion),sgf,sf2,abso1(id) c end if c 624 format('ffh',i4,i6,i4,2x,a4,2x,1p4e14.6) END DO C C H minus free-free opacity C ELSE IF(IT.EQ.3) THEN DO ID=1,ND T=TEMP(ID) ANE=ELEC(ID) c ABSOFF=(CFF1+CFFT(ID)*FRINV)*CFFN(ID)*FRINV ABSOFF=SFFHMI(POPUL(N0HN,ID),FR,T)*ANE ABSO1(ID)=ABSO1(ID)+ABSOFF EMIS1(ID)=EMIS1(ID)+ABSOFF END DO C C special evaluation of the cross-section C ELSE IF(IT.LT.0) THEN DO ID=1,ND ABSOFF=FFCROS(ION,IT,TEMP(ID),FR)* * POPUL(NNEXT(ION),ID)*ELEC(ID) ABSO1(ID)=ABSO1(ID)+ABSOFF EMIS1(ID)=EMIS1(ID)+ABSOFF END DO END IF 40 CONTINUE C C ******** 3. - additional continuum opacity (OPADD) C IF(IOPADD.NE.0) THEN ICALL=1 DO ID=1,ND CALL OPADD(0,ICALL,IJ,ID) ABSO1(ID)=ABSO1(ID)+ABAD EMIS1(ID)=EMIS1(ID)+EMAD SCAT1(ID)=SCAT1(ID)+SCAD c if(lpri.and.mod(id,20).eq.1) then c if(lpri.and.id.eq.50) then c write(6,623) id,ij,abad,abso1(id),emis1(id),scat1(id) c write(6,623) id,ij,abad,emad,scad c write(*,*) 'elec',elec(id),sigec(ij),elscat(id) c end if c 623 format('ad',i4,i6,1p4e14.6) END DO 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 DO ID=1,ND SG=PRFLIN(ID,IJ) ABSO1(ID)=ABSO1(ID)+SG*ABTRA(ITR,ID) EMIS1(ID)=EMIS1(ID)+SG*EMTRA(ITR,ID) END DO end if ENDIF IF(NLINES(IJ).GT.0) THEN C C the "overlapping" lines at the given frequency C DO 100 ILINT=1,NLINES(IJ) ITR=ITRLIN(ILINT,IJ) if(linexp(itr)) go to 100 iad=iadop(iatm(ilow(itr))) if(iad.gt.0.and..not.lfre) 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 A1=(FR-FREQ(IJ0))/(FREQ(IJ1)-FREQ(IJ0)) A2=UN-A1 DO ID=1,ND SG=A1*PRFLIN(ID,IJ1)+A2*PRFLIN(ID,IJ0) ABSO1(ID)=ABSO1(ID)+SG*ABTRA(ITR,ID) EMIS1(ID)=EMIS1(ID)+SG*EMTRA(ITR,ID) END DO c if(lpri.and.mod(id,20).eq.1) write(6,648) ij,id,itr,abso1(id), c * emis1(id),sg,abtra(itr,id) c 648 format('lin1',i8,i4,i7,1p4e14.6) 100 CONTINUE END IF 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) ABSO1(ID)=ABSO1(ID)+SG*ABTRA(ITR,ID) EMIS1(ID)=EMIS1(ID)+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)) ABSO1(ID)=ABSO1(ID)+SG*ABTRA(ITR,ID) EMIS1(ID)=EMIS1(ID)+SG*EMTRA(ITR,ID) END DO ENDIF c if(lpri.and.mod(id,20).eq.1) write(6,649) ij,id,itr,abso1(id), c * emis1(id),sg,abtra(itr,id) c 649 format('linodf',i8,i4,i7,1p4e14.6) 300 CONTINUE 400 CONTINUE ENDIF C c Lyman alpha and beta quasimolecular opacity c call quasim(ij) c C ---------------------------- C total opacity and emissivity C ---------------------------- C DO ID=1,ND ABSO1(ID)=ABSO1(ID)-EMIS1(ID)*XKF(ID)+SCAT1(ID) EMIS1(ID)=EMIS1(ID)*XKFB(ID) absot(id)=abso1(id) c if(lpri.and.mod(id,20).eq.1) write(6,641) ij,id,abso1(id), c * emis1(id),scat1(id) c 641 format('opac1',i8,i4,1p3e14.6) END DO c c --------------------------------- c hydrogen pacity from Gomez tables c --------------------------------- c call ghydop(ij) c c approximate opacity in Lyman lines C if(ioplym.gt.0) call lymlin(ij) C c -------------------------------------------------------- c contribution from precalculated background opacity table c --------------------------------------------------------` c if(ioptab.gt.0) then call opact1(ij) end if C C if needed, evaluate the opacity per gram C if(izscal.eq.0) then do id=1,nd absot(id)=abso1(id)*dens1(id) end do end if c if(ifprd.gt.0) call prd(ij) c if(iprcrs.gt.0) then ih=nfirst(ielh)+nprcrs-1 crs=abso1(iprcrs)/(popul(ih,iprcrs)*g(ih)* * 0.0265*4.1347e-15) do ii=nfirst(ielh),nlast(ielh) if(ii.ne.ih) then popul(ii,iprcrs)=pold(ii) do jj=ii+1,nnext(ielh) itrh=itra(ii,jj) abtra(itrh,iprcrs)=abtrh(itrh) end do end if end do end if RETURN END