181 lines
5.4 KiB
Fortran
181 lines
5.4 KiB
Fortran
SUBROUTINE OPACTR(IJ)
|
|
C =====================
|
|
C
|
|
C Absorption and emission coefficients and their derivatives
|
|
C for the use by Rybicki variant (RYBSOL)
|
|
C
|
|
C This procedure is very similar to OPACT1, the only difference 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, =SC for scattering) coefficient wrt T
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
PARAMETER (DELT=1.D-2)
|
|
common/dsctva/dsct1(mdepth),dscn1(mdepth)
|
|
common/hmolab/anh2(mdepth),anhm(mdepth)
|
|
common/grdpra/GRD(MDEPTH),pra(mdepth),pgs0(mdepth),ANTP(MDEPTH)
|
|
|
|
DIMENSION ABSOPP(MFREQ,MDEPTH),SCATPP(MFREQ,MDEPTH),
|
|
* EMISPP(MFREQ,MDEPTH),AES(MLEVEL,MLEVEL),BES(MLEVEL),
|
|
* ELEC0(MDEPTH),DENS0(MDEPTH),AN0(MDEPTH),
|
|
* BFABS(MLEVEL,MDEPTH),ELERAT(MDEPTH),
|
|
* POPUL0(MLEVEL,MDEPTH),POPLTE(MLEVEL)
|
|
c
|
|
C compute opacities at T+DELTA(T) - for derivatives
|
|
C
|
|
IF(IJ.EQ.1) THEN
|
|
C
|
|
C first, sab\ve absolute b-factors (for evaluating
|
|
C derivatives w.r.t. T)
|
|
C
|
|
IF(.NOT.LTE) THEN
|
|
LTE0=LTE
|
|
DO ID=1,ND
|
|
CALL WNSTOR(ID)
|
|
CALL SABOLF(ID)
|
|
CALL RATMAL(ID,AES,BES)
|
|
CALL LEVSOL(AES,BES,POPLTE,IIFOR,NLEVEL,0)
|
|
ELERAT(ID)=ELEC(ID)/(DENS(ID)/WMM(ID)+ELEC(ID))
|
|
DO I=1,NLEVEL
|
|
POPUL0(I,ID)=POPUL(I,ID)
|
|
BFABS(I,ID)=1.
|
|
IF(POPLTE(I).GT.0.) BFABS(I,ID)=POPUL(I,ID)/POPLTE(I)
|
|
END DO
|
|
END DO
|
|
LTE=LTE0
|
|
END IF
|
|
C
|
|
C opacities at T+DELTA(T) - for derivatives
|
|
C for that, one needs also to estimate ELEC and DENS for T+DELTA(T)
|
|
C
|
|
DO ID=1,ND
|
|
TEMP(ID)=TEMP(ID)*(UN+DELT)
|
|
ELEC0(ID)=ELEC(ID)
|
|
DENS0(ID)=DENS(ID)
|
|
an0(id)=dens(id)/wmm(id)+elec(id)
|
|
END DO
|
|
CALL TDPINI
|
|
DO ID=1,ND
|
|
T=TEMP(ID)
|
|
if(ifprad.gt.0) then
|
|
if(idisk.eq.0) then
|
|
if(id.eq.1) then
|
|
pgs(id)=dm(id)*(grav-grd(id)*(un+4.*delt))
|
|
else
|
|
pgs(id)=pgs(id-1)+grav*(dm(id)-dm(id-1))
|
|
* -pck*grd(id)*(un+4.*delt)
|
|
end if
|
|
end if
|
|
end if
|
|
end do
|
|
|
|
do id=1,nd
|
|
t=temp(id)
|
|
an=pgs(id)/bolk/t
|
|
if(idisk.gt.0) an=antp(id)
|
|
c IF(LTE) THEN
|
|
CALL ELDENS(ID,T,AN,ANE,ENRG,ENTT,WM,1)
|
|
c ELSE
|
|
c ANE=AN*ELERAT(ID)
|
|
c END IF
|
|
RHO=WMM(ID)*(AN-ANE)
|
|
DENS(ID)=RHO
|
|
ELEC(ID)=ANE
|
|
CALL WNSTOR(ID)
|
|
IF(LTE.OR.IFMOL.GT.0) THEN
|
|
CALL STEQEQ(ID,POP,1)
|
|
ELSE
|
|
CALL SABOLF(ID)
|
|
CALL RATMAL(ID,AES,BES)
|
|
CALL LEVSOL(AES,BES,POPLTE,IIFOR,NLEVEL,0)
|
|
DO I=1,NLEVEL
|
|
POPUL(I,ID)=POPLTE(I)*BFABS(I,ID)
|
|
END DO
|
|
END IF
|
|
END DO
|
|
CALL OPAINI(1)
|
|
IOPLY0=IOPLYM
|
|
IOPLYM=0
|
|
DO IJP=1,NFREQ
|
|
CALL OPACF1(IJP)
|
|
DO ID=1,ND
|
|
ABSOPP(IJP,ID)=ABSO1(ID)/DENS(ID)
|
|
EMISPP(IJP,ID)=EMIS1(ID)/DENS(ID)
|
|
SCATPP(IJP,ID)=SCAT1(ID)/DENS(ID)
|
|
END DO
|
|
END DO
|
|
IOPLYM=IOPLY0
|
|
C
|
|
C reset the original structural parameters
|
|
C
|
|
DO ID=1,ND
|
|
TEMP(ID)=TEMP(ID)/(UN+DELT)
|
|
END DO
|
|
CALL TDPINI
|
|
c
|
|
IF(IDISK.EQ.0) THEN
|
|
DO ID=1,ND
|
|
ELEC(ID)=ELEC0(ID)
|
|
DENS(ID)=DENS0(ID)
|
|
T=TEMP(ID)
|
|
AN=DENS0(ID)/WMM(ID)+ELEC0(ID)
|
|
PGS(ID)=AN*BOLK*T
|
|
END DO
|
|
ELSE
|
|
call pgset(1)
|
|
DO ID=1,ND
|
|
T=TEMP(ID)
|
|
AN=ANTP(ID)
|
|
CALL ELDENS(ID,T,AN,ANE,ENRG,ENTT,WM,1)
|
|
RHO=WMM(ID)*(AN-ANE)
|
|
DENS(ID)=RHO
|
|
ELEC(ID)=ANE
|
|
PGS(ID)=AN*BOLK*T
|
|
END DO
|
|
END IF
|
|
DO ID=1,ND
|
|
CALL WNSTOR(ID)
|
|
CALL STEQEQ(ID,POP,1)
|
|
c DO I=1,NLEVEL
|
|
c POPUL(I,ID)=POPUL0(I,ID)
|
|
c END DO
|
|
END DO
|
|
CALL OPAINI(1)
|
|
END IF
|
|
C
|
|
C opacity at the original T, and derivatives
|
|
C
|
|
CALL OPACF1(IJ)
|
|
DO ID=1,ND
|
|
ABSO1(ID)=ABSO1(ID)/DENS(ID)
|
|
SCAT1(ID)=SCAT1(ID)/DENS(ID)
|
|
EMIS1(ID)=EMIS1(ID)/DENS(ID)
|
|
DABT1(ID)=(ABSOPP(IJ,ID)-ABSO1(ID))/T/DELT
|
|
DSCT1(ID)=(SCATPP(IJ,ID)-SCAT1(ID))/T/DELT
|
|
XKF(ID)=EXP(-HKT1(ID)*FREQ(IJ))
|
|
XKF1(ID)=UN-XKF(ID)
|
|
XKFB(ID)=XKF(ID)*BNUE(IJ)
|
|
IF(LTE.OR.IFMOL.GT.0) THEN
|
|
PLAN=XKFB(ID)/XKF1(ID)
|
|
DPLAN=PLAN/XKF1(ID)*HKT1(ID)*FREQ(IJ)/TEMP(ID)
|
|
DEMT1(ID)=(DABT1(ID)-DSCT1(ID))*PLAN+
|
|
* (ABSO1(ID)-SCAT1(ID))*DPLAN
|
|
ELSE
|
|
DEMT1(ID)=(EMISPP(IJ,ID)-EMIS1(ID))/T/DELT
|
|
END IF
|
|
END DO
|
|
C
|
|
RETURN
|
|
END
|