258 lines
8.7 KiB
Fortran
258 lines
8.7 KiB
Fortran
SUBROUTINE OUTPRI
|
|
C =================
|
|
C
|
|
C Output on unit 6 (printer)
|
|
C
|
|
INCLUDE 'IMPLIC.FOR'
|
|
INCLUDE 'BASICS.FOR'
|
|
INCLUDE 'ATOMIC.FOR'
|
|
INCLUDE 'MODELQ.FOR'
|
|
INCLUDE 'ARRAY1.FOR'
|
|
common/grdpra/GRD(MDEPTH),pra(mdepth),pgs0(mdepth),ANTP(MDEPTH)
|
|
dimension aes(mlevel,mlevel),bes(mlevel),poplte(mlevel),
|
|
* bfab(mlevel,mdepth)
|
|
C
|
|
C ************ Print emergent radiation field on unit 13, namely
|
|
C
|
|
C FREQ(IJ) - value of frequency
|
|
C FLUX(IJ) - emergent flux, precisely the second moment H(freq)
|
|
C at the surface, in ergs/cm**2/s/sterad/Hz
|
|
C FH(IJ) - Eddington factor f(H), ie the ratio H/J, J is the
|
|
C mean intensity of radiation
|
|
C
|
|
WRITE(6,600) ITER-1
|
|
FLTT=SIG4P*TEFF**4
|
|
TOTF=0.
|
|
DO IJ=1,NFREQ
|
|
IJP=IJ
|
|
IF(ispodf.eq.0) IJP=JIK(IJ)
|
|
IF(IJX(IJP).NE.-1) THEN
|
|
WRITE(13,602) FREQ(IJP),FLUX(IJP),FH(IJP)
|
|
TOTF=TOTF+FLUX(IJP)*W(IJP)
|
|
FLAM=FLUX(IJP)*FREQ(IJP)*FREQ(IJP)/2.997925E18
|
|
write(14,614) 2.997925e18/freq(ijp),flam
|
|
END IF
|
|
END DO
|
|
WRITE(6,603) TOTF
|
|
C
|
|
C ************ For partial opacity table, print electron
|
|
C densities - actual, and that from opacity table
|
|
C
|
|
if(ioptab.ne.0) call eldenc
|
|
C
|
|
C ************ Print basic model parameters, namely
|
|
C
|
|
C ID - depth index
|
|
C DM(ID) - mass-depth variable (in g/cm**2)
|
|
C TEMP(ID) - temperature (in K)
|
|
C ELEC(ID) - electron density (cm**-3)
|
|
C AN - total particle number density (cm**-3)
|
|
C DENS(ID) - mass density (g/cm**3)
|
|
C P - total gas pressure (cgs)
|
|
C GR - radiative acceleration (cgs)
|
|
C FLTOT(ID)- total (integrated over frequencies) radiative flux
|
|
C
|
|
if(idisk.eq.0) then
|
|
WRITE(6,611)
|
|
else
|
|
WRITE(6,613)
|
|
end if
|
|
DO IJ=1,NFREQE
|
|
IJT=IJFR(IJ)
|
|
CALL OPACF1(IJT)
|
|
DO ID=1,ND
|
|
ABSOEX(IJ,ID)=ABSO1(ID)
|
|
END DO
|
|
END DO
|
|
DO ID=1,ND
|
|
C
|
|
C contributions from explicit (linearized) frequencies to the
|
|
C flux and radiation pressure
|
|
C
|
|
GRP=0.
|
|
FLEX=0.
|
|
IF(NFREQE.GT.0) THEN
|
|
DO IJ=1,NFREQE
|
|
RAD0(IJ)=RADEX(IJ,ID)
|
|
FK0(IJ)=FAKEX(IJ,ID)
|
|
ABSO0(IJ)=ABSOEX(IJ,ID)
|
|
IJT=IJFR(IJ)
|
|
WD0C=W(IJT)
|
|
IF(ID.EQ.1) THEN
|
|
FLUXW=FH(IJT)*RAD0(IJ)-HEXTRD(IJT)
|
|
IF(.NOT.LSKIP(ID,IJT)) GRP=GRP+W(ijt)*FLUXW*ABSO0(IJ)
|
|
FLEX=FLEX+WD0C*FLUXW
|
|
ELSE
|
|
RADM(IJ)=RADEX(IJ,ID-1)
|
|
FKM(IJ)=FAKEX(IJ,ID-1)
|
|
ABSOM(IJ)=ABSOEX(IJ,ID-1)
|
|
FRD=FK0(IJ)*RAD0(IJ)-FKM(IJ)*RADM(IJ)
|
|
IF(.NOT.LSKIP(ID,IJFR(IJ))) GRP=GRP+W(ijt)*FRD
|
|
DTAUM=(ABSO0(IJ)*DENS1(ID)+ABSOM(IJ)*DENS1(ID-1))*
|
|
* DELDM(ID-1)
|
|
FLEX=FLEX+WD0C*FRD/DTAUM
|
|
END IF
|
|
END DO
|
|
END IF
|
|
GRAD(ID)=GRP+FPRD(ID)
|
|
if(ifryb.gt.0) GRAD(ID)=GRD(ID)
|
|
IF(ID.EQ.1) THEN
|
|
GRAD(ID)=GRAD(ID)/DENS(ID)
|
|
ELSE
|
|
GRAD(ID)=GRAD(ID)/(DM(ID)-DM(ID-1))
|
|
END IF
|
|
FLTOT(ID)=FLEX
|
|
C
|
|
C other quantities
|
|
C
|
|
AN=DENS(ID)/WMM(ID)+ELEC(ID)
|
|
P=AN*TEMP(ID)*BOLK
|
|
IF(ID.LT.ND.AND.GRAD(ID).GT.0.)
|
|
* GR=LOG10(GRAD(ID)*4.1916825D-10)
|
|
FLTO=FLTOT(ID)+FLFIX(ID)+flxc(id)
|
|
flto=flrd(id)+flxc(id)
|
|
ptotal(id)=pgs(id)+pradt(id)
|
|
if(idisk.eq.0) then
|
|
WRITE(6,612) ID,DM(ID),TROSS(ID),TEMP(ID),ELEC(ID),DENS(ID),
|
|
* P,GR,
|
|
* flrd(id)/fltt,flxc(id)/fltt,flto/fltt
|
|
else
|
|
IF(ID.EQ.1) THEN
|
|
GRV=QGRAV*ZD(ID)
|
|
pgint=pgs(1)/dens(1)*dm(1)
|
|
ptint=ptotal(1)/dens(1)*dm(1)
|
|
ELSE
|
|
GRV=QGRAV*(ZD(ID)+ZD(ID-1))*HALF
|
|
pgint=pgint+(dm(id)-dm(id-1))*(pgs(id)/dens(id)+
|
|
* pgs(id-1)/dens(id-1))*half
|
|
ptint=ptint+(dm(id)-dm(id-1))*(ptotal(id)/dens(id)+
|
|
* ptotal(id-1)/dens(id-1))*half
|
|
END IF
|
|
GRVL=0.
|
|
IF(GRV.GT.0.) GRVL=LOG10(GRV)
|
|
HMECH=SIG4P*TEFF**4*(UN-THETAV(ID))
|
|
WRITE(6,622) ID,DM(ID),TROSS(ID),TEMP(ID),ELEC(ID),
|
|
* dens(id),pgs(id),
|
|
* flxc(id)/FLTO,FLTO,HMECH,hmech/flto,ZD(ID),GRVL,GR
|
|
622 format(i4,1p2e10.2,0pf10.1,1p10e10.2)
|
|
wbar=wbarm/dm(nd)
|
|
if(p.gt.0.) alpg=omeg32*wbar*dens(id)*viscd(id)/p
|
|
pto=ptotal(id)
|
|
if(pto.gt.0.) alpt=omeg32*wbar*dens(id)*viscd(id)/pto
|
|
disip=viscd(id)*dens(id)*edisc
|
|
write(98,698) id,dm(id),zd(id),abrosd(id),temp(id),pgint,
|
|
* ptotal(id),p,dens(id),disip,alpg,alpt
|
|
|
|
if(id.eq.nd) write(98,698) id,edisc,viscd(id),dens(id),p,
|
|
* omeg32,wbar
|
|
698 format(i3,1p11d11.4)
|
|
end if
|
|
END DO
|
|
C
|
|
if(idisk.eq.1) then
|
|
if(pgint.gt.0.) alpgav=omeg32*wbar/pgint*dm(nd)
|
|
if(ptint.gt.0.) alptav=omeg32*wbar/ptint*dm(nd)
|
|
write(6,606) omeg32,wbar,alpgav,alptav
|
|
606 format(//
|
|
* ' omega*3/2 ',1pe10.2/
|
|
* ' wbar ',1pe10.2/
|
|
* ' equivalent alpha for Pg ',1pe10.2/
|
|
* ' equivalent alpha for Ptot',1pe10.2)
|
|
end if
|
|
C
|
|
IF(.NOT.LTE) THEN
|
|
C
|
|
C ************ Print b-factors on unit 12 and
|
|
C "absolute" b-factors on unit 22
|
|
c the traditional b-factors are already computed (BFAC),
|
|
c here we compute the absolute ones
|
|
c
|
|
LTE=.TRUE.
|
|
DO ID=1,ND
|
|
CALL WNSTOR(ID)
|
|
CALL SABOLF(ID)
|
|
CALL RATMAL(ID,AES,BES)
|
|
CALL LEVSOL(AES,BES,POPLTE,IIFOR,NLEVEL,0)
|
|
DO I=1,NLEVEL
|
|
BFAB(I,ID)=1.
|
|
IF(POPLTE(I).GT.0.) BFAB(I,ID)=POPUL(I,ID)/POPLTE(I)
|
|
END DO
|
|
END DO
|
|
LTE=.FALSE.
|
|
idlte=idlt0
|
|
C
|
|
NUMP=NLEVEL+3
|
|
IF(IFMOL.GT.0) NUMP=NLEVEL+4
|
|
IF(IDISK.EQ.0) THEN
|
|
NUMPAR=NUMP
|
|
IF(IFMOL.GT.0) NUMPAR=-NUMPAR
|
|
WRITE(12,701) ND,NUMPAR
|
|
WRITE(12,702) (DM(ID),ID=1,ND)
|
|
WRITE(22,701) ND,NUMPAR
|
|
WRITE(22,704) (DM(ID),ID=1,ND)
|
|
DO ID=1,ND
|
|
IF(IFMOL.EQ.0) THEN
|
|
WRITE(12,703) TEMP(ID),ELEC(ID),DENS(ID),
|
|
* (BFAC(J,ID),J=1,NLEVEL)
|
|
WRITE(22,703) TEMP(ID),ELEC(ID),DENS(ID),
|
|
* (BFAB(J,ID),J=1,NLEVEL)
|
|
ELSE
|
|
WRITE(12,703) TEMP(ID),ELEC(ID),DENS(ID),TOTN(ID),
|
|
* (BFAC(J,ID),J=1,NLEVEL)
|
|
WRITE(22,703) TEMP(ID),ELEC(ID),DENS(ID),TOTN(ID),
|
|
* (BFAB(J,ID),J=1,NLEVEL)
|
|
END IF
|
|
END DO
|
|
ELSE
|
|
NUMPAR=NUMP+1
|
|
IF(IFMOL.GT.0) NUMPAR=-NUMPAR
|
|
WRITE(12,701) ND,NUMPAR
|
|
WRITE(12,702) (DM(ID),ID=1,ND)
|
|
WRITE(22,701) ND,NUMPAR
|
|
WRITE(22,704) (DM(ID),ID=1,ND)
|
|
DO ID=1,ND
|
|
IF(IFMOL.EQ.0) THEN
|
|
WRITE(12,703) TEMP(ID),ELEC(ID),DENS(ID),ZD(ID),
|
|
* (BFAC(J,ID),J=1,NLEVEL)
|
|
WRITE(22,703) TEMP(ID),ELEC(ID),DENS(ID),ZD(ID),
|
|
* (BFAB(J,ID),J=1,NLEVEL)
|
|
ELSE
|
|
WRITE(12,703) TEMP(ID),ELEC(ID),DENS(ID),TOTN(ID),
|
|
* ZD(ID),(BFAC(J,ID),J=1,NLEVEL)
|
|
WRITE(22,703) TEMP(ID),ELEC(ID),DENS(ID),TOTN(ID),
|
|
* ZD(ID),(BFAB(J,ID),J=1,NLEVEL)
|
|
END IF
|
|
END DO
|
|
END IF
|
|
END IF
|
|
C
|
|
600 FORMAT(/' ************************************'/
|
|
* ' FINAL RESULTS:'/' '/
|
|
* ' MODEL QUANTITIES IN',I3,'. ITERATION'/
|
|
* ' ************************************'/)
|
|
c 601 FORMAT(' IJ',8X,'FREQ',11X,'LAMBDA',8X,'FLUX',9X,
|
|
c * 'FH',4X,'LOG(FNU)',3X,'LOG(FLAM)'/)
|
|
602 FORMAT(1PE15.8,1PE12.4,0PF7.3)
|
|
603 FORMAT(' TOTAL SURFACE FLUX',1PD15.8)
|
|
611 FORMAT(/' ----------------------'/
|
|
* ' FINAL MODEL ATMOSPHERE'/
|
|
* ' ----------------------'/
|
|
* ' ID MASS',6X,'TAUROSS',5X,'TEMP',7X,'NE',9X,'DENS',
|
|
* 6X,'P_gas',4X,'LOG(G_rad)',3x,'RAD/TOT',3x,'CON/TOT',
|
|
* 2x,'(RAD+CON)/TOT'/)
|
|
612 FORMAT(1H ,I3,1P2E11.3,0PF10.1,1P6E11.3,3E13.5)
|
|
613 FORMAT(/' ---------------------'/
|
|
* ' FINAL DISK RING MODEL'/
|
|
* ' ---------------------'/
|
|
*' ID MASS',4X,'TAUROSS',5X,'TEMP',7X,'NE',7X,'RHO',
|
|
* 7X,'PGAS'5X,'CON/TOT RAD.FLX DISSIP',2X,
|
|
* 'FLX/DISSIP',4X,'Z',7X,
|
|
* 'LOG G',2X,'LOG G(RAD)'/)
|
|
614 FORMAT(F15.3,1pe15.3)
|
|
701 FORMAT(2I5)
|
|
702 FORMAT(1P8E10.3)
|
|
703 FORMAT(1P5E15.6)
|
|
704 FORMAT(1P6E13.6)
|
|
RETURN
|
|
END
|