522 lines
15 KiB
Fortran
522 lines
15 KiB
Fortran
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
|