132 lines
3.9 KiB
Fortran
132 lines
3.9 KiB
Fortran
SUBROUTINE ROSSTD(IJ)
|
|
C =====================
|
|
C
|
|
C Rosseland mean opacity
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ITERAT.FOR'
|
|
INCLUDE 'ALIPAR.FOR'
|
|
c DIMENSION TAUR(MDEPTH)
|
|
dimension pld(mdepth),abpld(mdepth)
|
|
|
|
C
|
|
C for IJ > 0 - contribution from the frequency IJ to the Rosseland
|
|
C opacity integrals
|
|
C
|
|
IF(IJ.GT.0) THEN
|
|
if(ij.eq.1) then
|
|
do id=1,nd
|
|
pld(id)=0.
|
|
abpld(id)=0.
|
|
end do
|
|
end if
|
|
IF(IJX(IJ).GE.0) THEN
|
|
DO ID=1,ND
|
|
PLAN=XKFB(ID)/XKF1(ID)*W(IJ)
|
|
DPLAN=PLAN/XKF1(ID)*FREQ(IJ)*HKT21(ID)
|
|
ABROSD(ID)=ABROSD(ID)+DPLAN/ABSO1(ID)
|
|
SUMDPL(ID)=SUMDPL(ID)+DPLAN
|
|
pld(id)=pld(id)+plan
|
|
abpld(id)=abpld(id)+(abso1(id)-scat1(id))*pld(id)
|
|
END DO
|
|
END IF
|
|
ELSE
|
|
C
|
|
C for IJ=0 - evaluation of the Rosseland optical depth and
|
|
C the radiative equilibrium division point
|
|
C
|
|
ID=1
|
|
IDR=0
|
|
TAURS(ID)=HALF*DEDM1*ABROSD(ID)*DENS(ID)
|
|
DO ID=2,ND
|
|
DTAUR=DELDM(ID-1)*(ABROSD(ID)+ABROSD(ID-1))
|
|
TAURS(ID)=TAURS(ID-1)+DTAUR
|
|
IF(TAURS(ID).LE.TAUDIV) IDR=ID
|
|
END DO
|
|
C
|
|
C in the last iteration, output of The Rosseland opacity and
|
|
C optical depth; skip the rest
|
|
C
|
|
IF(LFIN) THEN
|
|
DO ID=1,ND
|
|
TROSS(ID)=TAURS(ID)
|
|
abpl=abpld(id)/pld(id)/dens(id)
|
|
pll=sig4p*4.*temp(id)**4
|
|
WRITE(11,611) ID,DM(ID),TAURS(ID),ABROSD(ID),
|
|
* TEMP(ID),ELEC(ID),DENS(ID),
|
|
* pld(id),pll,abpl
|
|
END DO
|
|
611 FORMAT(I4,2X,1P6E12.4,2x,3e12.4)
|
|
RETURN
|
|
END IF
|
|
C
|
|
C determination of the radiative equilibrium parameters
|
|
C REINT and REDIF;
|
|
C REINT=1 for all ID .le. ND-idlst
|
|
C REDIF=1 for Rosseland optical depth > taudiv (taur.gt.taudiv)
|
|
C typically 0.1 - 0.5;
|
|
C - this value has been empirically shown to yield
|
|
C the best convergence rate
|
|
C
|
|
IF(ITER.GT.ITNDRE) RETURN
|
|
if(ndre.gt.1) then
|
|
c write(6,600)
|
|
do id=1,nd
|
|
if(id.lt.ndre) then
|
|
reint(id)=1.
|
|
redif(id)=0.
|
|
else
|
|
reint(id)=0.
|
|
redif(id)=1.
|
|
end if
|
|
end do
|
|
idr=ndre
|
|
else if(ndre.eq.-1) then
|
|
write(6,600)
|
|
do id=1,nd
|
|
reint(id)=1.
|
|
redif(id)=0.
|
|
end do
|
|
redif(1)=1.
|
|
else
|
|
c
|
|
DO ID=1,ND
|
|
REINT(ID)=UN
|
|
REDIF(ID)=UN
|
|
IF(ID.GT.ND-IDLST) REINT(ID)=0.
|
|
IF(ID.LT.IDR) THEN
|
|
REDIF(ID)=0.
|
|
IF(MOD(NDRE,10).EQ.-1) THEN
|
|
REDIF(ID)=TAURS(ID)
|
|
ELSE IF(MOD(NDRE,10).EQ.-2) THEN
|
|
REDIF(ID)=TAURS(ID)*TAURS(ID)
|
|
END IF
|
|
END IF
|
|
IF(NDRE.LE.-10) THEN
|
|
REDIF(ID)=REDIF(ID)/SIG4P/TEFF**4
|
|
REINT(ID)=REINT(ID)/ABPLAD(ID)*DENS(ID)
|
|
END IF
|
|
END DO
|
|
ID=1
|
|
IF(MOD(NDRE,10).EQ.-5) REDIF(ID)=UN
|
|
IF(NDRE.LE.-10) THEN
|
|
REDIF(ID)=REDIF(ID)/SIG4P/TEFF**4
|
|
END IF
|
|
NDRE=1
|
|
if(iter.eq.1) then
|
|
WRITE(6,601) IDR,IDR+1,ND-idlst
|
|
WRITE(10,601) IDR,IDR+1,ND-idlst
|
|
endif
|
|
endif
|
|
WRITE(6,601) IDR,IDR+1,ND-idlst
|
|
600 FORMAT(/' id redif reint'/)
|
|
601 FORMAT(/' SCHEME OF RADIATIVE EQUIL. DETERMINED IN RESOLV'/
|
|
* ' ONLY INTEGRAL EQUATION FOR ID <= ',I3/
|
|
* ' BOTH FOR ',I5,' <= ID <= ',I3/)
|
|
END IF
|
|
RETURN
|
|
END
|