SUBROUTINE OPACFD(IJ) C ===================== C C Absorption and emission coefficients, and their derivatives C C This procedure is very similar to OPACF1, the only differences is C the evaluation of derivatives C C Input: C IJ - depth index C Output: C ABSO1 - array of absorption coefficient C EMIS1 - array of emission coefficient C SCAT1 - array of scattering coefficient C Dxxy - array of derivatives of xx (=AB for absorption, =EM for C emission) coefficient wrt y (=T for temperature, =N for C electron density) C INCLUDE 'IMPLIC.FOR' INCLUDE 'BASICS.FOR' INCLUDE 'ATOMIC.FOR' INCLUDE 'MODELQ.FOR' INCLUDE 'ODFPAR.FOR' INCLUDE 'ALIPAR.FOR' INCLUDE 'ARRAY1.FOR' INCLUDE 'ITERAT.FOR' PARAMETER (C14=2.99793D14, CFF1=1.3727D-25) PARAMETER (DELT=1.D-3, DELR=1.D-3) common/hmolab/anh2(mdepth),anhm(mdepth) DIMENSION DABP0(MLEVEL),DEMP0(MLEVEL) DIMENSION DABN1A(MDEPTH),DEMN1A(MDEPTH),DSCN1A(MDEPTH) common/dsctva/dsct1(mdepth),dscn1(mdepth) common/rhoder/drhodt(mdepth) c if(ioptab.lt.0) then call opactd(ij) return end if C C initialize c DO ID=1,ND ELSCAT(ID)=ELEC(ID)*SIGEC(IJ) END DO C DO ID=1,ND ABSO1(ID)=0. EMIS1(ID)=0. SCAT1(ID)=ELSCAT(ID) DABT1(ID)=0. DEMT1(ID)=0. DABN1(ID)=SIGEC(IJ) DEMN1(ID)=0. ABSFF(ID)=0. DABFT(ID)=0. DABFN(ID)=0. DO II=1,NLEVEL DABP1(II,ID)=0. DEMP1(II,ID)=0. END DO END DO C C basic frequency- and depth-dependent quantities C FR=FREQ(IJ) lfre=fr.gt.frtabm 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 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 II=ILOW(ITR) 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 if(iifix(iatm(ii)).le.0) then DEMT1(ID)=DEMT1(ID)+EMISBF*DEMLT(ITR,ID) DEMN1(ID)=DEMN1(ID)+EMISBF*ELEC1(ID) DEMP1(JJ,ID)=DEMP1(JJ,ID)+EMISBF*POPINV(JJ,ID) if(ipzero(ii,id).eq.0) DABP1(II,ID)=DABP1(II,ID)+SGD end if END DO END IF END DO C C ******** 1b. bound-free contribution - with dielectronic rec. C else DO IBFT=1,NTRANC ITR=ITRBF(IBFT) II=ILOW(ITR) SG=CROSS(IBFT,IJ) iad=iadop(iatm(ii)) lcomop=iad.eq.0.or.(lfre.and.iad.gt.0) 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 if(iifix(iatm(ii)).le.0) then DEMT1(ID)=DEMT1(ID)+EMISBF*DEMLT(ITR,ID) DEMN1(ID)=DEMN1(ID)+EMISBF*ELEC1(ID) DEMP1(JJ,ID)=DEMP1(JJ,ID)+EMISBF*POPINV(JJ,ID) if(ipzero(ii,id).eq.0) DABP1(II,ID)=DABP1(II,ID)+SGD end if end if END DO END IF END DO end if C C ******** 2. free-free contribution C DO 40 ION=1,NION II=NNEXT(ION) IT=ITRA(II,II) C C hydrogenic (with Gaunt factor = 1 for IT=1; exact for IT=2) C (derivative of Gaunt factor wrt T is neglected) C iad=iadop(iatm(nnext(ion))) if(iad.gt.0.and..not.lfre) go to 40 IF(IT.LE.2) THEN DO ID=1,ND SF1=SFF3(ION,ID)*FR3INV SF2=SFF2(ION,ID) DSF2=DSFF(ION,ID) IF(FR.LT.FF(ION)) THEN SF2=UN/XKF(ID) DSF2=(HKT1(ID)*FR+HALF)*TEMP1(ID) END IF IF(IT.EQ.2) THEN X=C14*CHARG2(ION)/FR SF2=SF2-UN+GFREE1(ID,X) ELSE IF(IT.EQ.3) THEN CALL GFREED(ID,FR,CHARG2(ION),GFR,DGFR) SF2=SF2-UN+GFR DSF2=DSF2-(DGFR-(GFR-UN)*TEMP1(ID)*HALF)/SF2 END IF ABSOFF=SF1*SF2 ABSFF(ID)=ABSFF(ID)+ABSOFF if(iifix(iatm(ii)).eq.0) then DABFT(ID)=DABFT(ID)-ABSOFF*DSF2 DABFN(ID)=DABFN(ID)+ABSOFF*ELEC1(ID) DABPP=ABSOFF*POPINV(II,ID) DABP1(II,ID)=DABP1(II,ID)+DABPP DEMP1(II,ID)=DEMP1(II,ID)+DABPP end if END DO C C H minus free-free opacity C (all derivatives are neglected) C ELSE IF(IT.EQ.3) THEN DO ID=1,ND ABSOFF=SFFHMI(POPUL(NFIRST(IELH),ID),FR,TEMP(ID))* * ELEC(ID) ABSFF(ID)=ABSFF(ID)+ABSOFF END DO C C special evaluation of the cross-section C (all derivatives are neglected) C ELSE IF(IT.LT.0) THEN DO ID=1,ND ABSOFF=FFCROS(ION,IT,TEMP(ID),FR)* * POPUL(NNEXT(ION),ID)*ELEC(ID) ABSFF(ID)=ABSFF(ID)+ABSOFF END DO END IF 40 CONTINUE C ******** 3. - additional 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 DABT1(ID)=DABT1(ID)+DAT DEMT1(ID)=DEMT1(ID)+DET DABN1(ID)=DABN1(ID)+DAN DEMN1(ID)=DEMN1(ID)+DEN END DO END IF C C ----------------------- C total continuum opacity C ----------------------- C DO ID=1,ND ABSO1(ID)=ABSO1(ID)+ABSFF(ID) DABT1(ID)=DABT1(ID)+DABFT(ID) DABN1(ID)=DABN1(ID)+DABFN(ID) EMIS1(ID)=EMIS1(ID)+ABSFF(ID) DEMT1(ID)=DEMT1(ID)+DABFT(ID) DEMN1(ID)=DEMN1(ID)+DABFN(ID) END DO C C ******** 4. - opacity and emissivity in lines C LASER=ITER.GT.ITLAS IF(ISPODF.EQ.0) THEN IF(IJLIN(IJ).GT.0) THEN ITR=IJLIN(IJ) iad=iadop(iatm(ilow(itr))) if(iad.eq.0.or.(lfre.and.iad.gt.0)) then C C the "primary" line at the given frequency C ITR=IJLIN(IJ) II=ILOW(ITR) JJ=IUP(ITR) DO 50 ID=1,ND SG=PRFLIN(ID,IJ) SGPI=SG*ABTRA(ITR,ID) IF(SGPI.LE.0.AND.LASER) GO TO 50 SGPJ=SG*EMTRA(ITR,ID) ABSO1(ID)=ABSO1(ID)+SGPI EMIS1(ID)=EMIS1(ID)+SGPJ if(iifix(iatm(ii)).gt.0) go to 50 DEMT1(ID)=DEMT1(ID)+SGPJ*DEMLT(ITR,ID) DABP1(II,ID)=DABP1(II,ID)+SGPI*POPINV(II,ID) DEMP1(JJ,ID)=DEMP1(JJ,ID)+SGPJ*POPINV(JJ,ID) 50 CONTINUE end if ENDIF IF(NLINES(IJ).LE.0) GO TO 200 C C the "overlapping" lines at the given frequency C DO 100 ILINT=1,NLINES(IJ) ITR=ITRLIN(ILINT,IJ) if(linexp(itr)) goto 100 II=ILOW(ITR) JJ=IUP(ITR) IJ0=IFR0(ITR) iad=iadop(iatm(ii)) if(iad.gt.0.and..not.lfre) go to 100 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 80 ID=1,ND SG=A1*PRFLIN(ID,IJ1)+A2*PRFLIN(ID,IJ0) SGPI=SG*ABTRA(ITR,ID) IF(SGPI.LE.0.AND.LASER) GO TO 80 SGPJ=SG*EMTRA(ITR,ID) ABSO1(ID)=ABSO1(ID)+SGPI EMIS1(ID)=EMIS1(ID)+SGPJ if(iifix(iatm(ii)).gt.0) go to 80 DEMT1(ID)=DEMT1(ID)+SGPJ*DEMLT(ITR,ID) DABP1(II,ID)=DABP1(II,ID)+SGPI*POPINV(II,ID) DEMP1(JJ,ID)=DEMP1(JJ,ID)+SGPJ*POPINV(JJ,ID) 80 CONTINUE 100 CONTINUE 200 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) II=ILOW(ITR) JJ=IUP(ITR) iad=iadop(iatm(ii)) 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 310 ID=1,ND SGPI=PRFLIN(ID,KJ)*ABTRA(ITR,ID) IF(SGPI.LE.0.AND.LASER) GO TO 310 SGPJ=PRFLIN(ID,KJ)*EMTRA(ITR,ID) ABSO1(ID)=ABSO1(ID)+SGPI EMIS1(ID)=EMIS1(ID)+SGPJ if(iifix(iatm(ii)).gt.0) go to 310 DEMT1(ID)=DEMT1(ID)+SGPJ*DEMLT(ITR,ID) DABP1(II,ID)=DABP1(II,ID)+SGPI*POPINV(II,ID) DEMP1(JJ,ID)=DEMP1(JJ,ID)+SGPJ*POPINV(JJ,ID) 310 CONTINUE ELSE DO 320 ID=1,ND KJD=JIDI(ID) SG=EXP(XJID(ID)*SIGFE(KJD,KJ)+(UN-XJID(ID))* * SIGFE(KJD+1,KJ)) SGPI=SG*ABTRA(ITR,ID) IF(SGPI.LE.0.AND.LASER) GO TO 320 SGPJ=SG*EMTRA(ITR,ID) ABSO1(ID)=ABSO1(ID)+SGPI EMIS1(ID)=EMIS1(ID)+SGPJ if(iifix(iatm(ii)).gt.0) go to 320 DEMT1(ID)=DEMT1(ID)+SGPJ*DEMLT(ITR,ID) DABP1(II,ID)=DABP1(II,ID)+SGPI*POPINV(II,ID) DEMP1(JJ,ID)=DEMP1(JJ,ID)+SGPJ*POPINV(JJ,ID) 320 CONTINUE END IF 300 CONTINUE 400 CONTINUE END IF C c Lyman alpha and beta quasimolecular opacity c call quasim(ij) C C ------------------------------------------ C total opacity, emissivity, and derivatives C ------------------------------------------ C DO ID=1,ND DEMT1(ID)=DEMT1(ID)+EMIS1(ID)*FR*HKT21(ID) ABSO1(ID)=ABSO1(ID)-EMIS1(ID)*XKF(ID)+SCAT1(ID) DABN1(ID)=DABN1(ID)-DEMN1(ID)*XKF(ID) DABT1(ID)=DABT1(ID)-DEMT1(ID)*XKF(ID) EMIS1(ID)=EMIS1(ID)*XKFB(ID) DEMN1(ID)=DEMN1(ID)*XKFB(ID) DEMT1(ID)=DEMT1(ID)*XKFB(ID) DO II=1,NLEVEL DABP1(II,ID)=DABP1(II,ID)-DEMP1(II,ID)*XKF(ID) DEMP1(II,ID)=DEMP1(II,ID)*XKFB(ID) END DO absot(id)=abso1(id) END DO c IF(IOPLYM.GT.0) CALL LYMLIN(IJ) if(ifprd.gt.0) call prd(ij) C C --------------------------------------------- C C derivatives in the linearized explicit levels C IF(NLVEXP.LT.NLEVEL) THEN DO ID=1,ND DO I=1,NLVEXP DABP0(I)=0. DEMP0(I)=0. END DO DO I=1,NLEVEL if(iifix(iatm(i)).eq.0) then II=IIEXP(I) IF(II.GT.0) THEN DABP0(II)=DABP0(II)+DABP1(I,ID) DEMP0(II)=DEMP0(II)+DEMP1(I,ID) ELSE IF(II.LT.0) THEN DABP0(-II)=DABP0(-II)+DABP1(I,ID)*PP(I,ID) DEMP0(-II)=DEMP0(-II)+DEMP1(I,ID)*PP(I,ID) ELSE JJ=IIEXP(ILTREF(I,ID)) if(jj.gt.0) then DABP0(JJ)=DABP0(JJ)+DABP1(I,ID)*PP(I,ID) DEMP0(JJ)=DEMP0(JJ)+DEMP1(I,ID)*PP(I,ID) endif IF(IABS(IMODL(I)).LE.5) THEN DABT1(ID)=DABT1(ID)+DABP1(I,ID)*PT(I,ID) DEMT1(ID)=DEMT1(ID)+DEMP1(I,ID)*PT(I,ID) DABN1(ID)=DABN1(ID)+DABP1(I,ID)*PN(I,ID) DEMN1(ID)=DEMN1(ID)+DEMP1(I,ID)*PN(I,ID) END IF END IF end if END DO DO II=1,NLVEXP DABP1(II,ID)=DABP0(II) DEMP1(II,ID)=DEMP0(II) END DO END DO END IF c c contribution from the background opacity table c if(ioptab.gt.0) then c imodf=0 FR=FREQ(IJ) if(fr.lt.frtabm) then DO ID=1,ND T=TEMP(ID) T1=T*(UN+DELT) RHO=DENS(ID) RHO1=RHO*(UN+DELR) PLAN=XKFB(ID)/XKF1(ID) DPLAN=PLAN/XKF1(ID)*HKT1(ID)*FR/T CALL OPCTAB(FR,IJ,ID,T,RHO,AB,SC,SCT,imodf) CALL OPCTAB(FR,IJ,ID,T1,RHO,AB1,SC1,SCT1,imodf) CALL OPCTAB(FR,IJ,ID,T,RHO1,AB2,SC2,SCT2,imodf) ABSO1(ID)=ABSO1(id)+AB EMIS1(ID)=EMIS1(ID)+AB*PLAN scat1(id)=scat1(id)+sct c c derivatives w.r.t. temperature c DABTAB=(AB1-AB)/T/DELT dabtab=0. DABT1(ID)=DABT1(ID)+DABTAB DEMT1(ID)=DEMT1(ID)+AB*DPLAN+DABTAB*PLAN DSCT1(ID)=DSCT1(ID)+(SCT1-SCT)/T/DELT dabt1(id)=dabt1(id)+dsct1(id) c c derivatives w.r.t. density c DABN1A(ID)=(AB2-AB)/RHO/DELR DEMN1A(ID)=DABN1(ID)*PLAN c DSCN1A(ID)=(SCT2-SCT)/RHO/DELR c dabn1A(id)=dabn1a(id)+dscn1a(id) DABN1A(ID)=0. DEMN1A(ID)=0. DSCN1A(ID)=0. c c modify derivatives in case density is not a state parameter c IF(INHE.LE.0) THEN DABT1(ID)=DABT1(ID)+DABN1A(ID)*DRHODT(ID) DEMT1(ID)=DEMT1(ID)+DEMN1A(ID)*DRHODT(ID) DSCT1(ID)=DSCT1(ID)+DSCN1A(ID)*DRHODT(ID) ELSE DABN1(ID)=DABN1(ID)+DABN1A(ID) DEMN1(ID)=DEMN1(ID)+DEMN1A(ID) DSCN1(ID)=DSCN1(ID)+DSCN1A(ID) END IF c if(ifryb.gt.5) then abso1(id)=abso1(id)/dens(id) emis1(id)=emis1(id)/dens(id) scat1(id)=scat1(id)/dens(id) dabt1(id)=dabt1(id)/dens(id) demt1(id)=demt1(id)/dens(id) dsct1(id)=dsct1(id)/dens(id) end if c END DO end if 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)/dens(id) end do id=1 c if(mod(ij,1000).le.3) c * write(*,*) '+++++++opacfd',ij,abso1(id),absot(id) end if C C store quantities for explicit (linearized) frequencies C IF(IJEX(IJ).LE.0) RETURN IJE=IJEX(IJ) DO ID=1,ND ABSOEX(IJE,ID)=ABSO1(ID) EMISEX(IJE,ID)=EMIS1(ID) SCATEX(IJE,ID)=SCAT1(ID) DABTEX(IJE,ID)=DABT1(ID) DEMTEX(IJE,ID)=DEMT1(ID) DABNEX(IJE,ID)=DABN1(ID) DEMNEX(IJE,ID)=DEMN1(ID) DABMEX(IJE,ID)=DABM1(ID) DEMMEX(IJE,ID)=DEMM1(ID) DO II=1,NLVEXP DRCHEX(II,IJE,ID)=DABP1(II,ID) DRETEX(II,IJE,ID)=DEMP1(II,ID) END DO END DO C RETURN END