SpectraRust/tlusty/extracted/princ.f
2026-03-19 14:05:33 +08:00

157 lines
5.1 KiB
Fortran

SUBROUTINE PRINC
C ================
C
C Auxiliary output routine, which enables printing
C more detailed information about chosen individual transitions
C
C Input: NCT - number of transitions for which information is
C printed
C ICTR(I) - index of the I-the considered transition
C INFR(I) - index of the characteristic frequency point for
C the I-th transition (i.e., standardly, the point just
C shortward of edge for continua; or the line center
C for lines)
C - if =0, then the program choses the characteristic
C frequency itself, in the standard manner
C
C The printed table contains for each transition the following functions
C of depth:
C
C depth index;
C optical depth at the characteristic frequency;
C partial radiation force (actually, acceleration) - g(rad) -
C produced by the transition
C b-factor of the lower level;
C b-factor of the upper level;
C upward radiative rate;
C downward radiative rate;
C mean intensity at the characteristic frequency;
C Planck function at the characteristic frequency;
C total source function at the characteristic frequency;
C net line (or continuum) source function at the characteristic frequency;
C net heating rate (i.e heating minus cooling rate);
C
INCLUDE 'IMPLIC.FOR'
INCLUDE 'BASICS.FOR'
INCLUDE 'ATOMIC.FOR'
INCLUDE 'MODELQ.FOR'
INCLUDE 'ALIPAR.FOR'
PARAMETER (NPTR=30,CCOR=0.09,SIXTH=UN/6.)
DIMENSION ST(NPTR,MDEPTH),TAU(NPTR,MDEPTH),
* ABST(NPTR,MDEPTH),EMIT(NPTR,MDEPTH),SCTR(NPTR,MDEPTH),
* INFR(NPTR),ICTR(NPTR),ABM(NPTR),PRF(MFREQ),
* DWF(MFREQ)
C DIMENSION RADM(MFREQ),ABSM(MFREQ),FAK0(MFREQ),RAD0(MFREQ)
c
READ(44,*,END=200,ERR=200) NCT
READ(44,*,END=200,ERR=200) (ICTR(I),I=1,NCT)
READ(44,*,END=200,ERR=200) (INFR(I),I=1,NCT)
c
DO IC=1,NCT
ITR=ICTR(IC)
IFR=INFR(IC)
IF(LINE(ITR)) THEN
IF(IFR.EQ.0)
* INFR(IC)=(IFR0(ITR)+IFR1(ITR))/2
ELSE
IF(IFR.EQ.0) INFR(IC)=IFR1(ITR)
END IF
IFR=INFR(IC)
CALL OPACF1(IFR)
DO ID=1,ND
ABST(IC,ID)=ABSO1(ID)
EMIT(IC,ID)=EMIS1(ID)
SCTR(IC,ID)=SCAT1(ID)
END DO
END DO
C
DO ID=1,ND
T=TEMP(ID)
ANE=ELEC(ID)
SQT=SQRT(T)
ANES=EXP(SIXTH*LOG(ANE))
AACOR=CCOR*ANES/SQT
CALL SABOLF(ID)
DO I=1,NLEVEL
POP(I)=POPUL(I,ID)
END DO
DO IC=1,NCT
ITR=ICTR(IC)
IFR=INFR(IC)
I=ILOW(ITR)
J=IUP(ITR)
IF(LINE(ITR)) THEN
CALL LINPRO(ITR,ID,PRF)
GG=G(I)/G(J)
SG=PRF(IFR)
ELSE
GG=ANE*SBF(I)*EXP(-HK*FREQ(IFR)/T)
QZ=IZ(IEL(I))
MW=MCDW(ITR)
CALL DWNFR(MW,NFREQ,FR0(ITR),AACOR,ANE,QZ,FREQ,DWF)
SG=CROSS(ITRA(J,I),IFR)*DWF(IFR)
END IF
IJE=IJEX(IFR)
ESCT=SCTR(IC,ID)*RAD(IFR,ID)
ST(IC,ID)=(EMIT(IC,ID)+ESCT)/ABST(IC,ID)
IF(ID.EQ.1) THEN
TAU(IC,ID)=HALF *ABST(IC,ID)/DENS(ID)*DM(ID)
ELSE
TAU(IC,ID)=TAU(IC,ID-1)+(ABST(IC,ID)/DENS(ID)+ABM(IC))*
* (DM(ID)-DM(ID-1))*HALF
END IF
ABM(IC)=ABST(IC,ID)/DENS(ID)
END DO
END DO
C
DO IC=1,NCT
ITR=ICTR(IC)
I=ILOW(ITR)
J=IUP(ITR)
K=NNEXT(IEL(I))
IFR=INFR(IC)
IJE=IJEX(IFR)
FR=FREQ(IFR)
FR15=FR*1.D-15
BNU=BN*FR15*FR15*FR15
WRITE(16,600) ITR,IFR,FR,2.997925d18/fr
DO ID=1,ND
T=TEMP(ID)
TK=BOLK*T
SB=2.0706D-16/T/SQRT(T)*G(I)/G(K)*EXP(ENION(I)/TK)
IF(J.LT.K)
* SJ=2.0706D-16/T/SQRT(T)*G(J)/G(K)*EXP(ENION(J)/TK)
PI=POPUL(I,ID)
X=EXP(-HK*FR/T)
PLTE=SB*ELEC(ID)*POPUL(K,ID)
BI=PI/PLTE
BJ=1.
IF(J.LT.K) BJ=POPUL(J,ID)/SJ/ELEC(ID)/POPUL(K,ID)
PLANCK=BNU/(UN/X-UN)
IF(LINE(ITR)) THEN
RD=RRD(ITR,ID)*G(I)/G(J)
GG=G(I)/G(J)*POPUL(J,ID)
ELSE
GG=PLTE*X
RD=RRD(ITR,ID)
END IF
SL=BNU*GG/(PI-GG)
ggrad=0.
heat=0.
WRITE(16,601) ID,TAU(IC,ID),GGRAD,
* BI,BJ,RRU(ITR,ID),RD,
* RAD(IFR,ID),PLANCK,ST(IC,ID),SL,
* HEAT
END DO
END DO
C
600 FORMAT(/
* ' PARAMETERS FOR TRANSITION',I5,' IFR =' ,I5,
* ' FREQ =',1PD15.5,' Wavelength =',0pf11.3/
* 1H ,8X,'TAU',6X,'GR ',6X,'B-I',6X,'B-J',6X,' RU',
* 6X,' RD',6X,'RAD',5X,'PLANCK',5X,'STOT',6X,'SL',
* 5X,'HEAT'/)
601 FORMAT(1H ,I3,1P13D9.2)
200 RETURN
END