C***********************************************************************C C***********************************************************************C C C C AA GGGGG BBBBBBB SSSSS PPPPPPP EEEEEEE CCCCC C C AAA GG GG BB BB SS SS PP PP EE E CC CC C C AAAA GG BB BB SS PP PP EE CC C C AA AA GG GGG BBBBBB SSSSS PPPPPP EEEE CC C C AAAAAA GG G BB BB SS PP EE CC C C AA AA GG GG BB BB SS SS PP EE E CC CC C C AAAA AAAA GGGGG BBBBBBB SSSSS PPPP EEEEEEE CCCCC C C C C***********************************************************************C C***********************************************************************C C Subroutine package for calculating the emitted spectrum for a C C simple wind model generated by the program module AGBDUST C C Author: H.-P. Gail, Institute for Theoretical Astrophysics C C University of Heidelberg C C e-mail: gail@ita.uni-heidelberg.de C C Version: July 2010 C C***********************************************************************C C***********************************************************************C C Radiative transfer in a circumstellar dust shell C C***********************************************************************C SUBROUTINE CALCSPECT(DLS,TE) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(KDIM=1000,NDIM=80,IFDIM=11) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT COMMON/ABUND/EPSHE,EPSC,EPSN,EPSO,EPSNE,EPSNA,EPSMG,EPSAL,EPSSI, * EPSS,EPSAR,EPSCA,EPSFE,EPSNI COMMON/CSTAR/DLSTERN,RSTERN,TEFF COMMON/CHEMTYPE/RCO C -- Stellar parameter -- DLSTERN=DLS*SUNL TEFF=TE RSTERN=DSQRT(DLSTERN/(4D0*PI*SIG*TEFF**4)) RCO=EPSC/EPSO C check for the presence of dust species CALL CHECKDU C set initial temperature CALL STARTEM C calculate model ITT=0 CALL RADMODIT(ITT,IERR) C colors CALL COLOR RETURN END C**********************************************************************C C Preparing dust data and check for the important dust species C C**********************************************************************C SUBROUTINE CHECKDU PARAMETER(KDIM=1000,NSPDIM=11,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) LOGICAL FLDU COMMON/RMOD/RD(KDIM),TG(KDIM),RHO(KDIM),VV(KDIM),DKAP(KDIM), * GAM(KDIM),FV(KDIM,IFDIM),VDR(KDIM,IFDIM),AV(KDIM,IFDIM),KL,KK COMMON/CMOD/TAU(KDIM),TAL(KDIM),TAH(KDIM),CHIHT(KDIM),KLU COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) COMMON/CFDU/FDU(KDIM,NSPDIM) COMMON/CSTAR/DLSTERN,RSTERN,TEFF COMMON/CTE/CHM,RCOS,RCOC COMMON/CHEMTYPE/RCO CHARACTER*3 ADU C Store dust condensation fractions as they are used for radiative transfer DO K=1,KK DO N=1,NSPDIM FDU(K,N)=1D-50 END DO IF(RCO.LE.RCOS) THEN C -- M-stars -- FDU(K,1)=FDU(K,1)+FV(K,1) ! olivine FDU(K,2)=FDU(K,2)+FV(K,3) ! pyroxene FDU(K,3)=FDU(K,3)+FV(K,5) ! quartz FDU(K,4)=FDU(K,4)+FV(K,9) ! corundum FDU(K,5)=FDU(K,5)+FV(K,10) ! hibonite FDU(K,6)=FDU(K,6)+FV(K,11) ! spinel FDU(K,7)=FDU(K,7)+FV(K,6) ! iron FDU(K,8)=FDU(K,8)+FV(K,7) ! wuestite ELSE IF(RCO.LE.RCOC) THEN C -- S-stars -- FDU(K,3)=FDU(K,3)+FV(K,2) ! quartz FDU(K,7)=FDU(K,7)+FV(K,1) ! iron ELSE C -- C-stars -- FDU(K,7)=FDU(K,7)+FV(K,3) ! iron FDU(K,9)=FDU(K,9)+FV(K,1) ! carbon FDU(K,10)=FDU(K,10)+FV(K,2) ! silicon carbide FDU(K,11)=FDU(K,11)+FV(K,4) ! magnesium sulfide END IF END IF END DO C -- position of inner boundary of dust shell -- DO K=1,KK S=0D0 DO L=1,NSPDIM S=S+FDU(K,L) END DO IF(S.LE.1D-7) THEN KLU=K END IF END DO IF(KLU.EQ.KK) THEN PRINT*,'no dust is formed' STOP END IF KLU=KLU+1 C -- Which dust species are important ? -- NDON=0 DO I=1,NSPDIM DO K=KLU,KK IF(FDU(K,I).GT.1D-10) THEN NDON=NDON+1 NDONP(NDON)=I GOTO 100 END IF END DO 100 CONTINUE END DO C -- Which species condenses first ? -- FM=0D0 DO L1=1,NDON L=NDONP(L1) IF(FDU(KLU,L).GT.FM) THEN FM=FDU(KLU,L) IFIRSTC=L END IF END DO RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C S u b r o u t i n e s f o r c a l c u l a t i n g a C C C C r a d i a t i v e t r a n s f e r m o d e l C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C**********************************************************************C C Model iteration C C**********************************************************************C SUBROUTINE RADMODIT(ITT,IERR) IMPLICIT REAL*8(A-H,O-Z) REAL*8 VJNA,VKNA LOGICAL FLDU PARAMETER(KDIM=1000,NDIM=80,NSPDIM=11,IFDIM=11) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT COMMON/RMOD/RD(KDIM),TG(KDIM),RHO(KDIM),VV(KDIM),DKAP(KDIM), * GAM(KDIM),FV(KDIM,IFDIM),VDR(KDIM,IFDIM),AV(KDIM,IFDIM),KL,KK COMMON/CSTAR/DLSTERN,RSTERN,TEFF COMMON/CMOD/TAU(KDIM),TAL(KDIM),TAH(KDIM),CHIHT(KDIM),KLU COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) COMMON/CBDJ/HTA(KDIM),BN(KDIM,NSPDIM),SNUE(KDIM),VJN(KDIM), * VHN(KDIM),VKN(KDIM),VJNA(KDIM,NDIM),VKNA(KDIM,NDIM) COMMON/CTEMP/TS(KDIM,NSPDIM),TGH(KDIM) COMMON/FREQ/DLA(NDIM),DNU(NDIM),GNU(NDIM),DLAMREF,NN,NREF COMMON/CFLUX/FLUX(NDIM,4),HFLUX,HS,SFLUX COMMON/CRADM/VJ(KDIM),VH(KDIM),VK(KDIM),BV(KDIM,NSPDIM), * SPHF(KDIM),CHIBSP(KDIM,NSPDIM),CHIJSP(KDIM,NSPDIM), * CHIHSP(KDIM,NSPDIM) COMMON/CHEMTYPE/RCO COMMON/CTIT/ERRMH,ERRMT,ERRMJ,ITTMAX,IWH COMMON/CFDU/FDU(KDIM,NSPDIM) DIMENSION TW(NSPDIM) CHARACTER*3 ADU C -- Hnue of the star -- CALL STARSP(TEFF) HKORR=0D0 C -- initial values for moments J and K C (.NE.0 to avoid division by 0 in the first calculation of fEdd) -- DO N=1,NN DO K=1,KK VJNA(K,N)=0D0 VKNA(K,N)=0D0 END DO END DO IF(ITT.EQ.-1) RETURN C----------------------------------------------------------------------C C Model iteration C C----------------------------------------------------------------------C IERR=0 1000 ITT=ITT+1 IF(ITT.GT.ITTMAX) THEN PRINT *,'No convergence of the temperature correction procedure' IERR=1 GOTO 100 END IF C -- radiative transfer -- CALL RADMOD(HKORR,ERRJ,EH,KH) C -- temperature correction -- CALL KORRTEMP(ERRT,KT,KS) IF(IWH.NE.0) THEN 10 FORMAT('It:',I5,17X,'flux error:',1PE11.3,3X,'k:',I4,3X) WRITE(6,10) ITT,EH,KH 20 FORMAT(11X,'maximum temperature correction :',1PE11.3,3X, * 'k:',I4,' ks:',I4) WRITE(6,20) ERRT,KT,KS END IF C -- check of accuracy -- IF(DABS(ERRT).GT.ERRMT) GOTO 1000 C----------------------------------------------------------------------C C program end C C----------------------------------------------------------------------C 100 CONTINUE C -- optical depths averages -- TAH(KLU)=0D0 TAL(KLU)=0D0 DO K=KLU+1,KK H=.5D0*(RD(K)-RD(K-1)) TAH(K)=TAH(K-1)+H*(RHO(K-1)*CHIHT(K-1)+RHO(K)*CHIHT(K)) TAL(K)=TAL(K-1)+H*RSTERN*RSTERN*(RHO(K)*CHIHT(K)/(RD(K)**2)+ * RHO(K-1)*CHIHT(K-1)/(RD(K-1)**2)) END DO CALL MKTAU(0) C -- show optical depth variables -- 30 FORMAT('tau(',F3.1,'):',1PE10.3,3X,' tauH:',E10.3,3X * ,' tauL:',E10.3) WRITE(6,30) DLA(NREF),TAU(KK),TAH(KK),TAL(KK) C -- show temperature at inner boundary -- c 80 FORMAT('first condensate: ',A3,3X,'at Tc:',F12.6,3X) c WRITE(6,80) ADU(IFIRSTC),TS(KLU,IFIRSTC) C -- model output -- OPEN(10,FILE='huel.sp') 2 FORMAT(1P20E13.5) 3 FORMAT('#',I5,1P3E13.5) WRITE(10,3) NN,TEFF,DLSTERN/SUNL,RCO DO N=1,NN WRITE(10,2) DLA(N),DNU(N),FLUX(N,1),FLUX(N,2),FLUX(N,3) END DO CLOSE(10) 40 FORMAT(2X,'r',13X,9('T_',A3,8X)) OPEN(10,FILE='huel.mo') WRITE(10,40) (ADU(NDONP(L1)),L1=1,NDON) DO K=KLU,KK DO L1=1,NDON L=NDONP(L1) IF(FDU(K,L).GT.0D0) THEN TW(L1)=TS(K,L) ELSE TW(L1)=0D0 END IF END DO X1=RD(K)*RD(K)*VH(K)/(RSTERN*RSTERN*HS) WRITE(10,2) RD(K)/RSTERN,(TW(L1),L1=1,NDON),TGH(K),X1, * VK(K)/VJ(K) END DO CLOSE(10) END C**********************************************************************C C Radiative transfer in the dust shell C C**********************************************************************C SUBROUTINE RADMOD(HKORR,ERRJ,EH,KH) IMPLICIT REAL*8(A-H,O-Z) REAL*8 VJNA,VKNA LOGICAL FLDU PARAMETER(KDIM=1000,NDIM=80,NSPDIM=11,IFDIM=11) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT COMMON/RMOD/RD(KDIM),TG(KDIM),RHO(KDIM),VV(KDIM),DKAP(KDIM), * GAM(KDIM),FV(KDIM,IFDIM),VDR(KDIM,IFDIM),AV(KDIM,IFDIM),KL,KK COMMON/CFDU/FDU(KDIM,NSPDIM) COMMON/CSTAR/DLSTERN,RSTERN,TEFF COMMON/CMOD/TAU(KDIM),TAL(KDIM),TAH(KDIM),CHIHT(KDIM),KLU COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) COMMON/CBDJ/HTA(KDIM),BN(KDIM,NSPDIM),SNUE(KDIM),VJN(KDIM), * VHN(KDIM),VKN(KDIM),VJNA(KDIM,NDIM),VKNA(KDIM,NDIM) COMMON/CTEMP/TS(KDIM,NSPDIM),TGH(KDIM) COMMON/FREQ/DLA(NDIM),DNU(NDIM),GNU(NDIM),DLAMREF,NN,NREF COMMON/CFLUX/FLUX(NDIM,4),HFLUX,HS,SFLUX COMMON/DUEXT/ABSDU(NDIM,NSPDIM),STRDU(NDIM,NSPDIM) COMMON/DUEXT0/ABSDU0(NDIM,NSPDIM),STRDU0(NDIM,NSPDIM) COMMON/CDUGR/XIC0,XISIC0,XIIRO0,XIC,XISIC,XIIRO COMMON/CRADM/VJ(KDIM),VH(KDIM),VK(KDIM),BV(KDIM,NSPDIM), * SPHF(KDIM),CHIBSP(KDIM,NSPDIM),CHIJSP(KDIM,NSPDIM), * CHIHSP(KDIM,NSPDIM) COMMON/CTIT/ERRMH,ERRMT,ERRMJ,ITTMAX,IWH DIMENSION DKA(NSPDIM),DKS(NSPDIM) CHARACTER*3 ADU HSTERN=DLSTERN/(16D0*PI**2*RSTERN**2) C -- initial values -- PRINT*,'radiative transfer' DO K=KLU,KK VJ(K)=0D0 VH(K)=0D0 VK(K)=0D0 DO L=1,NSPDIM BV(K,L)=0D0 CHIBSP(K,L)=1D-30 CHIJSP(K,L)=1D-30 CHIHSP(K,L)=1D-30 END DO CHIHT(K)=0D0 END DO C Temperaturmodifikation zur Ber"ucksichtung der Einstrahlung TBH=(4D0*PI/SIG*(HSTERN+HKORR))**.25D0 DO N=1,NN FLUX(N,2)=.25D0*BNU(DNU(N),TBH) END DO HKORR=0D0 KJ=0 NJ=0 ERRJ=0D0 C---------------------------- begin of frequency loop ------------------ DO 100 N=1,NN ANU=DNU(N) C -- absorption & scattering -- DO L1=1,NDON L=NDONP(L1) DKA(L)=ABSDU(N,L) DKS(L)=DMIN1(STRDU(N,L),1D3*ABSDU(N,L)) END DO C -- source function -- DO K=KLU,KK DO L1=1,NDON L=NDONP(L1) BN(K,L)=BNU(ANU,TS(K,L)) END DO END DO C -- solution of radiative transfer equation -- CALL SOLVERAD(DKA,DKS,N) C -- frequency integrations -- DO K=KLU,KK DUSTE=0D0 DO L1=1,NDON L=NDONP(L1) BV(K,L)=BV(K,L)+GNU(N)*BN(K,L) CHIBSP(K,L)=CHIBSP(K,L)+GNU(N)*BN(K,L)*DKA(L)*FDU(K,L) CHIJSP(K,L)=CHIJSP(K,L)+GNU(N)*VJN(K)*DKA(L)*FDU(K,L) CHIHT(K)=CHIHT(K)+GNU(N)*VHN(K)*(DKA(L)+DKS(L))*FDU(K,L) C Q/A for drift velocity CHIHSP(K,L)=CHIHSP(K,L)+ * GNU(N)*VHN(K)*(ABSDU0(N,L)+STRDU0(N,L)) END DO VJ(K)=VJ(K)+GNU(N)*VJN(K) VH(K)=VH(K)+GNU(N)*VHN(K) VK(K)=VK(K)+GNU(N)*VKN(K) C -- store J_nue, K_nue -- IF(RD(K).LT.100D0*RSTERN.AND.DLA(N).LT.50D0. * AND.VJNA(K,N).GT.1D-15) THEN X=VJN(K)/VJNA(K,N)-1D0 IF(DABS(X).GT.DABS(ERRJ)) THEN ERRJ=X KJ=K NJ=N END IF END IF FFJ=1.70D0 VJNA(K,N)=FFJ*VJN(K)+(1D0-FFJ)*VJNA(K,N) VKNA(K,N)=FFJ*VKN(K)+(1D0-FFJ)*VKNA(K,N) END DO C -- spectrum of the dust shell -- FLUX(N,1)=HFLUX*(RD(KK)/RSTERN)**2 HKORR=HKORR+SFLUX*GNU(N)*(RD(KK)/RSTERN)**2 C -- optical depths of dust shell -- CALL MKTAU(N) FLUX(N,3)=TAU(KK) C---------------------------- end of frequency loop -------------------- 100 CONTINUE WRITE(6,'(3HeJ:2X,F8.3,1PE11.3,2X,I5)') DLA(NJ),ERRJ,KJ C -- frequency averages -- DO K=KLU,KK CHIHT(K)=CHIHT(K)/VH(K) DO L1=1,NDON L=NDONP(L1) CHIBSP(K,L)=RHO(K)*CHIBSP(K,L)/BV(K,L) CHIJSP(K,L)=RHO(K)*CHIJSP(K,L)/VJ(K) CHIHSP(K,L)=CHIHSP(K,L)/VH(K) END DO END DO C -- flux error -- EH=0D0 FL=1D0/(RSTERN*RSTERN*HS) TAH(KLU-1)=0D0 DO K=KLU,KK TAH(K)=TAH(K-1)+.5D0*(RD(K)-RD(K-1))*(CHIHT(K)+CHIHT(K-1)) X=FL*RD(K)*RD(K)*VH(K)-1D0 IF(DABS(X).GT.DABS(EH)) THEN EH=X KH=K END IF END DO IF(IWH.NE.0) THEN 10 FORMAT(1X,'flux error:',1PE11.3,3X,'k:',I4,3X) WRITE(6,10) EH,KH END IF C -- sphericity factor -- SPHF(KK)=0D0 FSPHA=3D0-VJ(KK)/VK(KK) DO K=KK-1,KLU,-1 FSPHN=3D0-VJ(K)/VK(K) SPHF(K)=SPHF(K+1)+.5D0*(FSPHN+FSPHA)*DLOG(RD(K+1)/RD(K)) FSPHA=FSPHN END DO DO K=KLU,KK SPHF(K)=DEXP(SPHF(K)) END DO RETURN END C**********************************************************************C C Temperature correction C C**********************************************************************C SUBROUTINE KORRTEMP(ERRT,KT,KS) IMPLICIT REAL*8(A-H,O-Z) REAL*8 VJNA,VKNA LOGICAL FLDU PARAMETER(KDIM=1000,NDIM=80,NSPDIM=11,IFDIM=11) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT COMMON/RMOD/RD(KDIM),TG(KDIM),RHO(KDIM),VV(KDIM),DKAP(KDIM), * GAM(KDIM),FV(KDIM,IFDIM),VDR(KDIM,IFDIM),AV(KDIM,IFDIM),KL,KK COMMON/CFDU/FDU(KDIM,NSPDIM) COMMON/CSTAR/DLSTERN,RSTERN,TEFF COMMON/CMOD/TAU(KDIM),TAL(KDIM),TAH(KDIM),CHIHT(KDIM),KLU COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) COMMON/CTEMP/TS(KDIM,NSPDIM),TGH(KDIM) COMMON/CRADM/VJ(KDIM),VH(KDIM),VK(KDIM),BV(KDIM,NSPDIM), * SPHF(KDIM),CHIBSP(KDIM,NSPDIM),CHIJSP(KDIM,NSPDIM), * CHIHSP(KDIM,NSPDIM) COMMON/CBDJ/HTA(KDIM),BN(KDIM,NSPDIM),SNUE(KDIM),VJN(KDIM), * VHN(KDIM),VKN(KDIM),VJNA(KDIM,NDIM),VKNA(KDIM,NDIM) COMMON/FREQ/DLA(NDIM),DNU(NDIM),GNU(NDIM),DLAMREF,NN,NREF COMMON/CFLUX/FLUX(NDIM,4),HFLUX,HS,SFLUX COMMON/CTIT/ERRMH,ERRMT,ERRMJ,ITTMAX,IWH DIMENSION TW(NSPDIM) CHARACTER*3 ADU ERRT=0D0 FL=HS*RSTERN**2 HA=FL/(RD(KK)*RD(KK)) AA=-(VH(KK)-HA) DO 100 K=KK,KLU,-1 DO L1=1,NDON L=NDONP(L1) X=CHIJSP(K,L)/CHIBSP(K,L) DB1=X*AA DB2=X*VJ(K)-BV(K,L) TA=TS(K,L) TN4=TA**4+PI*(DB1+0.6D0*DB2)/SIG IF(TN4.GT.0D0) THEN TS(K,L)=DSQRT(DSQRT(TN4)) ELSE TS(K,L)=.9D0*TA END IF X=(TS(K,L)-TA)/TS(K,L) IF(DABS(X).GT.DABS(ERRT).AND.FDU(K,L).GT.1D-3) THEN IF(TS(K,L).GT.1D2.AND.L.LT.10) THEN ERRT=X KT=K KS=L END IF END IF END DO C --- gas temperature --- TN4=PI*(VJ(K)+AA)/SIG TGH(K)=DSQRT(DSQRT(TN4)) X1=RD(K)*RD(K)*VH(K)/FL C --- integration --- IF(K.GT.KLU) THEN HN=FL/(RD(K-1)*RD(K-1)) HL=RD(K)-RD(K-1) RFEDD=VJ(K)/(VK(K)*SPHF(K)) AA=AA-.5D0*RFEDD*(SPHF(K)*RHO(K)*CHIHT(K)*(VH(K)-HA)+ * SPHF(K-1)*RHO(K-1)*CHIHT(K-1)*(VH(K-1)-HN))*HL HA=HN END IF 100 CONTINUE IF(KLU.EQ.1) GOTO 110 C extrapolate dust temperatures to smaller r DO K=KLU-1,1,-1 RF=DSQRT(RD(KLU)/RD(K)) DO L1=1,NDON L=NDONP(L1) TS(K,L)=TS(KLU,L)*RF END DO TGH(K)=TGH(KLU)*RF END DO 110 DO K=1,KK TG(K)=TGH(K) END DO IF(IWH.NE.0) THEN 20 FORMAT(11X,'maximum temperature correction :',1PE11.3,3X, * 'k:',I4,' ks:',I4) WRITE(6,20) ERRT,KT,KS END IF C set dummy value for T for species not present DO K=KLU,KK DO N=1,NSPDIM TW(N)=TGH(K) END DO DO L1=1,NDON L=NDONP(L1) TW(L)=TS(K,L) END DO DO N=1,NSPDIM TS(K,N)=TW(N) END DO END DO DO K=KLU,KK DO L1=1,NDON L=NDONP(L1) IF(FDU(K,L).GT.1D-6) THEN TW(L1)=TS(K,L) ELSE TW(L1)=0D0 END IF END DO HH=FL/RD(K)**2 END DO RETURN END C**********************************************************************C C Temperature initialization C C**********************************************************************C SUBROUTINE STARTEM IMPLICIT REAL*8(A-H,O-Z) PARAMETER(KDIM=1000,NSPDIM=11,IFDIM=11) COMMON/RMOD/RD(KDIM),TG(KDIM),RHO(KDIM),VV(KDIM),DKAP(KDIM), * GAM(KDIM),FV(KDIM,IFDIM),VDR(KDIM,IFDIM),AV(KDIM,IFDIM),KL,KK COMMON/CMOD/TAU(KDIM),TAL(KDIM),TAH(KDIM),CHIHT(KDIM),KLU COMMON/CTEMP/TS(KDIM,NSPDIM),TGH(KDIM) DO K=KK,KLU,-1 IF(TG(K).LE.0D0) THEN PRINT*,'T<=0 not possible K:',K,' T:',TG(K) STOP END IF DO L=1,NSPDIM TS(K,L)=TG(K) END DO TGH(K)=TG(K) END DO RETURN END C**********************************************************************C C tau-scale C C**********************************************************************C SUBROUTINE MKTAU(N) IMPLICIT REAL*8(A-H,O-Z) LOGICAL FLDU PARAMETER(KDIM=1000,NDIM=80,NSPDIM=11,IFDIM=11) COMMON/RMOD/RD(KDIM),TG(KDIM),RHO(KDIM),VV(KDIM),DKAP(KDIM), * GAM(KDIM),FV(KDIM,IFDIM),VDR(KDIM,IFDIM),AV(KDIM,IFDIM),KL,KK COMMON/CFDU/FDU(KDIM,NSPDIM) COMMON/CMOD/TAU(KDIM),TAL(KDIM),TAH(KDIM),CHIHT(KDIM),KLU COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) COMMON/DUEXT/ABSDU(NDIM,NSPDIM),STRDU(NDIM,NSPDIM) COMMON/FREQ/DLA(NDIM),DNU(NDIM),GNU(NDIM),DLAMREF,NN,NREF CHARACTER*3 ADU IF(N.EQ.0) THEN NF=NREF ELSE NF=N END IF TAU(KLU)=0D0 DO K=KLU+1,KK DUSTA=0D0 DUSTS=0D0 DO L1=1,NDON L=NDONP(L1) DUSTA=DUSTA+ABSDU(NF,L)*FDU(K,L) DUSTS=DUSTS+STRDU(NF,L)*FDU(K,L) END DO DUSTE=DUSTA+DUSTS TAU(K)=TAU(K-1)+.5D0*(RD(K)-RD(K-1))*(RHO(K)+RHO(K-1))*DUSTE END DO RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C S u b r o u t i n e s f o r r a d i a t i v e t r a n s f e r C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C**********************************************************************C C Kirchhoff-Planck-function C C**********************************************************************C DOUBLE PRECISION FUNCTION BNU(ANU,TT) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT X=4.7995D-11*ANU/TT RRL=ANU/CL Y=3.973D-16*RRL*RRL*RRL IF(X.LT.1D-3) THEN BNU=Y/(X+.5D0*X*X) RETURN END IF IF(X.LT.15D0) THEN BNU=Y/(DEXP(X)-1D0) RETURN END IF IF(X.LT.200D0) THEN BNU=Y*DEXP(-X) RETURN END IF BNU=0D0 RETURN END C**********************************************************************C C Initializing the frequency scale C C**********************************************************************C SUBROUTINE FREQSCALE(DUSTFILE) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(NDIM=80) COMMON/FREQ/DLA(NDIM),DNU(NDIM),GNU(NDIM),DLAMREF,NN,NREF DIMENSION DL(33) DATA CL/2.99793D+10/ DATA DL/500D0,230D0,100D0,48D0,35D0,30D0,27D0,24D0, * 20D0,18D0,15D0,13D0,12D0,11.7D0,-11.53D0,-11.52D0,-11.51D0, * 11.3D0,11D0,10.5D0,9.5D0,8.5D0,7.5D0,5.0D0,3.4D0,3.1D0,2.8D0, * 1.5D0,1.0D0,.70D0,.50D0,.40D0,0D0/ CHARACTER*1 CH CHARACTER*20 FILENAM,DUSTFILE 2 FORMAT(A1,A20) C Read filename FILENAM=DUSTFILE OPEN(10,ERR=9010,FILE=FILENAM,STATUS='OLD') IF(INE.NE.0) GOTO 9010 10 READ(10,2) CH,FILENAM IF(CH.EQ.'!'.OR.CH.EQ.'c'.OR.CH.EQ.'C') GOTO 10 IF(CH.EQ.' ') GOTO 12 11 PRINT*,'--> Error in file DUSTLIST in line',IZ WRITE(6,2) CH,FILENAM STOP 12 CLOSE(10) OPEN(10,ERR=9010,FILE=FILENAM,STATUS='OLD') NN=1 100 READ(10,*) DLA(NN) IF(DLA(NN).NE.0D0) THEN NN=NN+1 GOTO 100 END IF CLOSE(10) NN=NN-1 N=1 L=1 GNU(1)=0D0 200 IF(N.EQ.NN) GOTO 400 IF(L.GT.NDIM-2) GOTO 9000 IF(DLA(N).LT.0D0) GOTO 300 C -- Integration with Simpsons rule -- DNU(L)=CL/(DLA(N)*.0001D0) DNU(L+2)=CL/(DABS(DLA(N+1))*.0001D0) DNU(L+1)=.5D0*(DNU(L)+DNU(L+2)) H=DNU(L+1)-DNU(L) GNU(L)=GNU(L)+H/3D0 GNU(L+1)=4D0*H/3D0 GNU(L+2)=H/3D0 L=L+2 N=N+1 GOTO 200 C -- Integration with trapezoidal rule -- 300 DNU(L)=-CL/(DLA(N)*1D-4) DNU(L+1)=CL/(DABS(DLA(N+1))*1D-4) H=DNU(L+1)-DNU(L) GNU(L)=GNU(L)+H/2D0 GNU(L+1)=H/2D0 L=L+1 N=N+1 GOTO 200 400 NN=L PRINT*,'Number of frequency grid-points:',NN DO N=1,NN DLA(N)=1D4*CL/DNU(N) END DO NREF=0 DO N=1,NN IF(DABS((DLA(N)-DLAMREF)/DLA(N)).LT.1D-4) NREF=N END DO IF(NREF.EQ.0) THEN PRINT*,'Reference wavelength not element of lambda-grid' STOP END IF RETURN 9000 PRINT*,'Too small dimension NDIM' STOP 9010 PRINT*,'In FREQSCALE: file >',FILENAM,'< does not exist' STOP END C**********************************************************************C C Hnue at the stellar surface: black body C C**********************************************************************C SUBROUTINE STARSP(TEFF) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(NDIM=80) COMMON/FREQ/DLA(NDIM),DNU(NDIM),GNU(NDIM),DLAMREF,NN,NREF COMMON/CFLUX/FLUX(NDIM,4),HFLUX,HS,SFLUX HS=0D0 DO N=1,NN FLUX(N,2)=.25D0*BNU(DNU(N),TEFF) HS=HS+FLUX(N,2)*GNU(N) END DO RETURN END C**********************************************************************C C External Radiation field C C**********************************************************************C SUBROUTINE ISRF(MISRF) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(NDIM=80) COMMON/FREQ/DLA(NDIM),DNU(NDIM),GNU(NDIM),DLAMREF,NN,NREF COMMON/CFLUX/FLUX(NDIM,4),HFLUX,HS,SFLUX GOTO(100,200,300),MISRF+1 C -- no Radiation from outside -- 100 WIS=0D0 GOTO 210 C -- diluted black body radiation -- 200 WIS=1D-14 210 TISRAD=1D+4 DO N=1,NN FLUX(N,4)=WIS*BNU(DNU(N),TISRAD) END DO RETURN C -- Modell of the interstellar radiation field -- 300 DO N=1,NN FLUX(N,4)=RADISM(DLA(N)) END DO RETURN END C**********************************************************************C C Interstellar radiation field Inue in the solar neighbourhood C C Mathis, Mezger, Panagia A & A 128, 212 C C (other models e.g. in Draine & Bertoldi ApJ 468, 269) C C**********************************************************************C DOUBLE PRECISION FUNCTION RADISM(DLAM) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT DIMENSION DL(44),DJ(44) DATA DL/.091,.100,.110,.130,.143,.180,.200,.210,.216,.230,.250, * .346,.435,.550,.700,.900,1.20,1.80,2.20,2.40,3.40,4.00, * 5.00,6.00,8.00,10.0,12.0,20.0,25.0,30.0,40.0,50.0,60.0, * 70.0,80.0,90.0,100.,150.,200.,300.,400.,600.,800.,1000./ DATA DJ/1.07D-2,1.47D-2,2.04D-2,2.05D-2,1.82D-2,1.24D-2,1.04D-2, * 9.61D-3,9.17D-3,8.25D-3,7.27D-3,1.30D-2,1.50D-2,1.57D-2, * 1.53D-2,1.32D-2,9.26D-3,4.06D-3,2.41D-3,1.89D-3,6.49D-4, * 3.79D-4,1.76D-4,9.21D-5,4.33D-5,3.26D-5,2.21D-5,7.13D-6, * 1.14D-5,1.70D-5,3.06D-5,4.08D-5,4.76D-5,4.51D-5,4.20D-5, * 3.62D-5,3.16D-5,1.36D-5,5.61D-6,1.39D-6,5.61D-7,1.32D-7, * 4.20D-8,2.14D-8/ DATA II/44/ ANU=1D4*CL/DLAM UM=CL/(4D0*PI*ANU*ANU) IF(DLAM.GT.1000D0) THEN C Extrapolation for lam>1000 mu RADISM=UM*DJ(II)*(1D3/DLAM)**1.5D0 ELSE CALL FIPOL(DL,DLAM,II,IO,IL) RADISM=UM*(DJ(IL)+(DJ(IO)-DJ(IL))*(DLAM-DL(IL))/(DL(IO)-DL(IL))) END IF RETURN END C***********************************************************************C C Solution of the radiative transfer equation C C***********************************************************************C SUBROUTINE SOLVERAD(DKA,DKS,N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(KDIM=1000,NDIM=80,NSPDIM=11,IFDIM=11) REAL*8 VJNA,VKNA COMMON/CBDJ/HTA(KDIM),BN(KDIM,NSPDIM),SNUE(KDIM),VJN(KDIM), * VHN(KDIM),VKN(KDIM),VJNA(KDIM,NDIM),VKNA(KDIM,NDIM) COMMON/CFEAU/SD(KDIM),DMU(KDIM),DMA(KDIM),U(KDIM),V(KDIM), * RR(KDIM),PR(KDIM),RRA(KDIM),PRA(KDIM) COMMON/RMOD/RD(KDIM),TG(KDIM),RHO(KDIM),VV(KDIM),DKAP(KDIM), * GAM(KDIM),FV(KDIM,IFDIM),VDR(KDIM,IFDIM),AV(KDIM,IFDIM),KL,KK COMMON/CFLUX/FLUX(NDIM,4),HFLUX,HS,SFLUX COMMON/CSTAR/DLSTERN,RSTERN,TEFF COMMON/CMOD/TAU(KDIM),TAL(KDIM),TAH(KDIM),CHIHT(KDIM),KLU COMMON/CSRP/KPS,KPC DIMENSION DKA(NSPDIM),DKS(NSPDIM) DO K=KLU,KK VJN(K)=0D0 VHN(K)=0D0 VKN(K)=0D0 DMA(K)=1D0 END DO RED=1D0 RS=RED*RSTERN AIP=4D0*FLUX(N,2)/RED**2 IPL=KLU IFL=1 C -- rays hitting the star -- IRD=-1 DO M=0,KPS P=M*0.999D0*RS/KPS CALL FEAUTRIER(AIP,P,IPL,IFL,IRD,DKA,DKS,N) END DO C -- rays between star and dust shell -- IRD=0 DO M=0,KPC P=RS+M*(RD(KLU-1)-RS)/KPC CALL FEAUTRIER(AIP,P,IPL,IFL,IRD,DKA,DKS,N) END DO C -- rays through the dust shell --- DO IPL=KLU,KK P=RD(IPL) CALL FEAUTRIER(AIP,P,IPL,IFL,IRD,DKA,DKS,N) END DO RETURN END C**********************************************************************C C Feautrier-Method C C**********************************************************************C SUBROUTINE FEAUTRIER(AIP,P,IPL,IFL,IRD,DKA,DKS,N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(KDIM=1000,NDIM=80,NSPDIM=11,IFDIM=11) REAL*8 VJNA,VKNA LOGICAL FLDU COMMON/RMOD/RD(KDIM),TG(KDIM),RHO(KDIM),VV(KDIM),DKAP(KDIM), * GAM(KDIM),FV(KDIM,IFDIM),VDR(KDIM,IFDIM),AV(KDIM,IFDIM),KL,KK COMMON/CFDU/FDU(KDIM,NSPDIM) COMMON/CMOD/TAU(KDIM),TAL(KDIM),TAH(KDIM),CHIHT(KDIM),KLU COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) COMMON/CBDJ/HTA(KDIM),BN(KDIM,NSPDIM),SNUE(KDIM),VJN(KDIM), * VHN(KDIM),VKN(KDIM),VJNA(KDIM,NDIM),VKNA(KDIM,NDIM) COMMON/CFEAU/SD(KDIM),DMU(KDIM),DMA(KDIM),U(KDIM),V(KDIM), * RR(KDIM),PR(KDIM),RRA(KDIM),PRA(KDIM) CHARACTER*3 ADU COMMON/CFLUX/FLUX(NDIM,4),HFLUX,HS,SFLUX COMMON/FREQ/DLA(NDIM),DNU(NDIM),GNU(NDIM),DLAMREF,NN,NREF DIMENSION DKA(NSPDIM),DKS(NSPDIM) DIMENSION DKEXT(KDIM),DERKE(KDIM) IF(IPL.EQ.KK) GOTO 100 C -- s and angle mu -- DO K=IPL,KK SD(K)=DSQRT(RD(K)*RD(K)-P*P) DMU(K)=SD(K)/RD(K) END DO C -- absorption, scattering, source function -- DO K=IPL,KK EMISS=0D0 DUSTA=0D0 DUSTS=0D0 DO L1=1,NDON L=NDONP(L1) DUSTA=DUSTA+DKA(L)*FDU(K,L) DUSTS=DUSTS+DKS(L)*FDU(K,L) EMISS=EMISS+DKA(L)*FDU(K,L)*BN(K,L) END DO DUSTA=DUSTA*RHO(K) DUSTS=0.5D0*DUSTS*RHO(K) C thermal emission of all components EMISS=EMISS*RHO(K) C scattering with dipole characterics FEDD=VKNA(K,N)/VJNA(K,N) SCATT=DUSTS*.375D0*((3D0-DMU(K)*DMU(K))*VJNA(K,N)+ * (3D0*DMU(K)*DMU(K)-1D0)*VKNA(K,N)) C isotropic scattering c SCATT=DUSTS*VJNA(K,N) DKEXT(K)=DUSTA+DUSTS SNUE(K)=(EMISS+SCATT)/DKEXT(K) END DO C -- optical depths intervals -- DO K=IPL+1,KK-1 DERKE(K)=(DKEXT(K+1)-DKEXT(K-1))/(RD(K+1)-RD(K-1)) END DO DO K=IPL,KK-1 IF(K.EQ.IPL) THEN HTA(K)=(SD(K+1)-SD(K))*(2D0*DKEXT(K)+DKEXT(K+1))/3D0 ELSE HTA(K)=.5D0*(SD(K+1)-SD(K))*(DKEXT(K)+DKEXT(K+1) * +(SD(K+1)-SD(K))*(DERKE(K+1)-DERKE(K))/6D0) END IF END DO IF(IPL.EQ.KK-1) GOTO 200 C -- inner boundary -- HH2=.5D0*HTA(IPL)**2 IF(IRD.NE.0) THEN U(IPL)=1D0/(1D0+HTA(IPL)+HH2) V(IPL)=(AIP*HTA(IPL)+HH2*SNUE(IPL))*U(IPL) ELSE U(IPL)=1D0/(1D0+HH2) V(IPL)=HH2*SNUE(IPL)*U(IPL) END IF C -- grid points -- DO K=IPL+1,KK-1 HL=HTA(K-1) HR=HTA(K) HM=HL+HR c XX=2D0/(HL*HM) c Q=-2D0/(HL*HR)-1D0+XX*U(K-1) c U(K)=-2D0/(HR*HM*Q) c V(K)=-(SNUE(K)+XX*V(K-1))/Q C Genauigkeit O(h^4) D3=(HR-HL)/(3.D0*HM) C Genauigkeit O(h^2) C D3=0.D0 XX=2.D0/(HL*HM)+D3 Q=-2.D0/(HL*HR)-1.D0+XX*U(K-1) U(K)=-(2.D0/(HR*HM)-D3)/Q V(K)=-(SNUE(K)+D3*(SNUE(K+1)-SNUE(K-1))+XX*V(K-1))/Q END DO C -- outer boundary -- AIM=FLUX(N,4) HH2=.5D0*HTA(KK-1)*HTA(KK-1) Q=1D0+HTA(KK-1)+HH2-U(KK-1) PR(KK)=(HTA(KK-1)*AIM+HH2*SNUE(KK)+V(KK-1))/Q C -- calculation downwards -- DO K=KK-1,IPL,-1 PR(K)=PR(K+1)*U(K)+V(K) END DO C -- calculating R -- IF(IRD.EQ.0) THEN RR(IPL)=0D0 ELSE RR(IPL)=AIP-PR(IPL) END IF DO K=IPL+1,KK-1 HL=HTA(K-1) HR=HTA(K) HM=HL+HR RR(K)=HR*PR(K-1)/(HL*HM)-(HR-HL)*PR(K)/(HL*HR)-HL*PR(K+1)/ * (HR*HM) END DO RR(KK)=PR(KK)-AIM GOTO 300 C -- the two outermost rays -- 100 PR(KK)=0D0 RR(KK)=0D0 GOTO 300 200 AIM=0D0 XX=.5D0*HTA(KK-1)*(SNUE(KK)+SNUE(KK-1)) PR(KK-1)=XX+AIM RR(KK-1)=0D0 PR(KK)=PR(KK-1) RR(KK)=XX C -- angular integrations -- 300 IF(IFL.EQ.1) THEN IFL=2 HFLUX=0D0 SFLUX=0D0 GOTO 400 END IF DO K=IPL,KK GMU=(DMA(K)-DMU(K))*.5D0 VJN(K)=VJN(K)+(PR(K)+PRA(K))*GMU VHN(K)=VHN(K)+(RR(K)*DMU(K)+RRA(K)*DMA(K))*GMU VKN(K)=VKN(K)+(PR(K)*DMU(K)*DMU(K)+PRA(K)*DMA(K)*DMA(K))*GMU END DO HFLUX=HFLUX+0.5D0*((PR(KK)+RR(KK))*DMU(KK)+ * (PRA(KK)+RRA(KK))*DMA(KK))*GMU IF(IRD.EQ.-1) SFLUX=SFLUX+0.5D0*((PR(IPL)-RR(IPL))*DMU(IPL)+ * (PRA(IPL)-RRA(IPL))*DMA(IPL))*GMU 400 DO K=IPL,KK PRA(K)=PR(K) RRA(K)=RR(K) DMA(K)=DMU(K) END DO RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C S u b r o u t i n e s f o r c a l c u l a t i n g c o l o r s C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C**********************************************************************C C Colors C C**********************************************************************C SUBROUTINE COLOR IMPLICIT REAL*8(A-H,O-Z) PARAMETER(NDIM=80) COMMON/FREQ/DLA(NDIM),DNU(NDIM),GNU(NDIM),DLAMREF,NN,NREF COMMON/CFLUX/FLUX(NDIM,4),HFLUX,HS,SFLUX COMMON/CCOL/COLU,COLB,COLV,COLR,COLI,COLJ,COLK,COLL,COLM,COLN, * COLJBB,COLHBB,COLKBB,COLLBB,COLLSBB,COLMBB DIMENSION DLU(13),FU(13),DLB(11),FB(11),DLV(15),FV(15), * DLR(23),FR(23),DLI(26),FI(26),DLJ(31),FJ(31),DLK(17),FK(17), * DLL(15),FL(15),DLM(23),FM(23),DLN(34),FN(34), * DLJBB(21),FJBB(21),DLHBB(20),FHBB(20),DLKBB(25),FKBB(25), * DLLBB(21),FLBB(21),DLSBB(20),FLSBB(20),DLMBB(17),FMBB(17) DIMENSION FA(NDIM),FIP(34),FRE(34) C---------------------------------------------------------------------C C UBVRIJKLMN Colors according to Johnson, ApJ 141, 923 C C---------------------------------------------------------------------C DATA DLU/.42,.41,.40,.39,.38,.37,.36,.35,.34,.33,.32,.31,.30/ DATA FU/.00,.01,.05,.36,.73,.97,1.0,.97,.93,.84 ,.61,.10,.00/ DATA DLB/.56,.54,.52,.50,.48,.46,.44,.42,.40,.38,.36/ DATA FB/.00,.04,.15,.36,.58,.79,.94,1.0,.92,.11,.00/ DATA DLV/.74,.72,.70,.68,.66,.64,.62,.60,.58,.56,.54,.52,.50, * .48,.46/ DATA FV/.00,.01,.01,.01,.02,.08,.20,.40,.62,.72,.98,.91,.38, * .02,.00/ DATA DLR/.96,.94,.92,.90,.88,.86,.84,.82,.80,.78,.76,.74,.72,.70, * .68,.66,.64,.62,.60,.58,.56,.54,.52/ DATA FR/.00,.01,.02,.04,.06,.11,.17,.31,.42,.57,.73,.85,.94,1.0, * .98,.94,.88,.79,.69,.50,.28,.06,.00/ DATA DLI/1.20,1.16,1.14,1.12,1.10,1.08,1.06,1.04,1.01,1.00,.98, * .96,.94,.92,.90,.88,.86,.84,.82,.80,.78,.76,.74,.72,.70,.68/ DATA FI/.00,.02,.03,.05,.08,.10,.15,.20,.28,.36,.47,.58,.71,.84, * .93,.98,1.0,.99,.98,.96,.76,.56,.36,.17,.01,.00/ DATA DLJ/1.56,1.54,1.52,1.50,1.48,1.46,1.44,1.42,1.40,1.38,1.36, * 1.34,1.32,1.30,1.28,1.26,1.24,1.22,1.20,1.18,1.16,1.14,1.12, * 1.10,1.08,1.06,1.04,1.02,1.00,.98,.96/ DATA FJ/.00,.02,.06,.09,.14,.27,.46,.60,.66,.70,.70,.68,.66,.63, * .63,.64,.75,.93,.85,.80,.78,.78,.85,.93,.62,.35,.16,.06,.03, * .02,.00/ DATA DLK/2.60,2.55,2.50,2.45,2.40,2.35,2.30,2.25,2.20,2.15, * 2.10,2.05,2.00,1.95,1.90,1.85,1.80/ DATA FK/.00,.08,.46,.84,.95,.95,.94,.96,.97,.95,.96,.98,1.0,.95, * .48,.10,.00/ DATA DLL/4.3,4.2,4.1,4.0,3.9,3.8,3.7,3.6,3.5,3.4,3.3,3.2,3.1, * 3.0,2.9/ DATA FL/.00,.04,.40,.65,.59,.65,.81,.89,1.0,1.0,1.0,.95,.68,.14, * .00/ DATA DLM/6.3,6.2,6.1,6.0,5.9,5.8,5.7,5.6,5.5,5.4,5.3,5.2,5.1, * 5.0,4.9,4.8,4.7,4.6,4.5,4.4,4.3,4.2,4.1/ DATA FM/.00,.02,.09,.19,.21,.26,.32,.49,.60,.68,.74,.79,.76,.78, * .84,.89,.92,.93,1.0,.48,.25,.04,.00/ DATA DLN/13.6,13.4,13.2,13.0,12.8,12.6,12.4,12.2,12.0,11.8, * 11.6,11.4,11.2,11.0,10.8,10.6,10.4,10.2,10.0,9.8,9.6,9.4,9.2, * 9.0,8.8,8.6,8.4,8.2,8.0,7.8,7.6,7.4,7.2,7.0/ DATA FN/.43,.46,.50,.54,.63,.71,.76,.76,.76,.71,.68,.65,.64,.61, * .60,.60,.65,.72,.73,.69,.63,.60,.68,.79,.92,.98,.89,.84,.76, * .70,.59,.34,.10,.00/ C---------------------------------------------------------------------C C JHKLM colors according to Bessell & Brett, PASP 106, 1134 C C---------------------------------------------------------------------C DATA DLJBB/1.44,1.42,1.40,1.38,1.36,1.34,1.32,1.30,1.28,1.26, * 1.24,1.22,1.20,1.18,1.16,1.14,1.12,1.10,1.08,1.06,1.04/ DATA FJBB/.00,.03,.07,0.11,.04,.68,.79,.91,.94,.88,.83,.81,.77, * .73,.63,.47,.32,.42,.11,.02,.00/ DATA DLHBB/1.84,1.82,1.80,1.78,1.76,1.74,1.72,1.70,1.68,1.66, * 1.64,1.62,1.60,1.58,1.56,1.54,1.52,1.50,1.48,1.46/ DATA FHBB/.00,.02,.52,.71,.84,.87,.95,.99,.99,.99,.99,.99,.99, * .95,.98,.94,.86,.44,.15,.00/ DATA DLKBB/2.48,2.44,2.40,2.38,2.34,2.32,2.30,2.28,2.26,2.24, * 2.22,2.20,2.18,2.16,2.14,2.12,2.10,2.08,2.06,2.04,2.02,2.00, * 1.98,1.96,1.94/ DATA FKBB/.00,.10,.64,.75,.84,.88,.91,.96,.97,.98,.96,.94,.95, * .94,.94,.90,.85,.77,.55,.74,.55,.30,.20,.12,.00/ DATA DLLBB/3.84,3.80,3.76,3.72,3.68,3.64,3.60,3.56,3.52,3.48, * 3.44,3.40,3.36,3.32,3.28,3.24,3.20,3.16,3.12,3.08,3.04/ DATA FLBB/.00,.04,.19,.65,.86,.84,.84,.84,.88,.85,.70,.61,.50, * .41,.61,.50,.30,.38,.09,.02,.00/ DATA DLSBB/4.20,4.16,4.12,4.08,4.04,4.00,3.96,3.92,3.88,3.84, * 3.80,3.76,3.72,3.68,3.64,3.60,3.56,3.52,3.48,3.44/ DATA FLSBB/.00,.24,.54,.63,.76,.97,.97,.97,.91,.86,.88,.86,.82, * .85,.91,.90,.80,.19,.02,.00/ DATA DLMBB/5.08,5.04,5.00,4.96,4.92,4.88,4.84,4.80,4.76,4.72, * 4.68,4.64,4.60,4.56,4.52,4.48,4.44/ DATA FMBB/.00,.03,.07,.23,.37,.44,.37,.33,.35,.16,.44,.50,.39, * .30,.34,.13,.00/ C -- shell spectrum -- DO N=1,NN FA(N)=FLUX(N,1) END DO C-------------------------- Johnson colors -----------------------------C C -- U -- CALL ICOL(COLU,DNU,FA,DLU,FIP,FU,FRE,NN,13) C -- B -- CALL ICOL(COLB,DNU,FA,DLB,FIP,FB,FRE,NN,11) C -- V -- CALL ICOL(COLV,DNU,FA,DLV,FIP,FV,FRE,NN,15) C -- R -- CALL ICOL(COLR,DNU,FA,DLR,FIP,FR,FRE,NN,23) C -- I -- CALL ICOL(COLI,DNU,FA,DLI,FIP,FI,FRE,NN,26) C -- J -- CALL ICOL(COLJ,DNU,FA,DLJ,FIP,FJ,FRE,NN,31) C -- K -- CALL ICOL(COLK,DNU,FA,DLK,FIP,FK,FRE,NN,17) C -- L -- CALL ICOL(COLL,DNU,FA,DLL,FIP,FL,FRE,NN,15) C -- M -- CALL ICOL(COLM,DNU,FA,DLM,FIP,FM,FRE,NN,23) C -- N -- CALL ICOL(COLN,DNU,FA,DLN,FIP,FN,FRE,NN,34) C------------------------ Bessell & Brett colors -----------------------C 1234 continue C -- J -- CALL ICOL(COLJBB,DNU,FA,DLJBB,FIP,FJBB,FRE,NN,21) C -- H -- CALL ICOL(COLHBB,DNU,FA,DLHBB,FIP,FHBB,FRE,NN,20) C -- K -- CALL ICOL(COLKBB,DNU,FA,DLKBB,FIP,FKBB,FRE,NN,25) C -- L -- CALL ICOL(COLLBB,DNU,FA,DLLBB,FIP,FLBB,FRE,NN,21) C -- L' -- CALL ICOL(COLLSBB,DNU,FA,DLSBB,FIP,FLSBB,FRE,NN,20) C -- M -- CALL ICOL(COLMBB,DNU,FA,DLMBB,FIP,FMBB,FRE,NN,17) COLJBB=-COLJBB-48.58D0-0.90D0 COLHBB=-COLHBB-48.58D0-1.37D0 COLKBB=-COLKBB-48.58D0-1.88D0 COLLBB=-COLLBB-48.58D0-2.77D0 COLLSBB=-COLLSBB-48.58D0-2.97D0 COLMBB=-COLMBB-48.58D0-3.42D0 RETURN END C----------------------------------------------------------------------C SUBROUTINE ICOL(CO,DNU,FA,DL,FN,FC,FREQ,NN,NF) C DNU,FA,NN : spectrum C DL,FC,NF : weight function C FN,FREQ : used for interpolation IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT DIMENSION DNU(NN),FA(NN),DL(NF),FN(NF),FC(NF),FREQ(NF) C Frequencies of the filter function DO N=1,NF FREQ(N)=CL/(DL(N)*1D-4) END DO CALL IPVEC(1,DNU,FA,FREQ,FN,NN,NF) CO=0D0 SF=0D0 DO N=1,NF-1 CO=CO+.5D0*(FREQ(N+1)-FREQ(N))*FC(N)*(FN(N+1)+FN(N)) SF=SF+.5D0*(FREQ(N+1)-FREQ(N))*FC(N) END DO CO=2.5D0*DLOG10(CO/SF) RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C S u b r o u t i n e s f o r d u s t e x t i n c t i o n C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C***********************************************************************C C Read table of Q_a, Q_s C C***********************************************************************C SUBROUTINE DUSTABS(DUSTFILE) IMPLICIT REAL*8(A-H,O-Z) LOGICAL FLDU PARAMETER(NDIM=80,NSPDIM=11,NRDIM=7000) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT COMMON/FREQ/DLA(NDIM),DNU(NDIM),GNU(NDIM),DLAMREF,NN,NREF COMMON/DUEXT0/ABSDU0(NDIM,NSPDIM),STRDU0(NDIM,NSPDIM) COMMON/COLI/V0OLI,A0OLI,EPSOLI,ALFOLI,ALFOLEX COMMON/CPYR/V0PYR,A0PYR,EPSPYR,ALFPYR,ALFPYEX COMMON/CQUA/V0QUA,A0QUA,EPSQUA,ALFQUA COMMON/CCOR/V0COR,A0COR,EPSCOR,ALFCOR COMMON/CHIB/V0HIB,A0HIB,EPSHIB,ALFHIB COMMON/CSPI/V0SPI,A0SPI,EPSSPI,ALFSPI COMMON/CIRO/V0IRO,A0IRO,EPSIRO,ALFIRO COMMON/CCAR/V0CAR,A0CAR,EPSCAR,ALFCAR COMMON/CSIC/V0SIC,A0SIC,EPSSIC,ALFSIC COMMON/CMGS/V0MGS,A0MGS,EPSMGS,ALFMGS COMMON/CWUE/V0WUE,A0WUE,EPSWUE,ALFWUMG,ALFWUFE COMMON/CAW/AWSIO,AWMG,AWFE,AWH2O,AWC2H2,AWSI,AWSIS,AWH2S,AWAL, * AWFEO,AWMGO,AWCA CHARACTER*3 ADU,SKU COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) DIMENSION IDU(NSPDIM),VALFA(NSPDIM),VV0(NSPDIM) CHARACTER*20 FILENAM,DUSTFILE CHARACTER*1 CH,CB 1 FORMAT(1P30E10.3) 2 FORMAT(A1,A3,A1,A20) 3 FORMAT(F8.5) 4 FORMAT('alpha=',F8.5) 5 FORMAT(A1,A12) DO N=1,NSPDIM IDU(N)=0 FLDU(N)=.TRUE. END dO FILENAM=DUSTFILE OPEN(11,ERR=9010,FILE=FILENAM,STATUS='OLD') 10 READ(11,5) CH,FILENAM IF(CH.EQ.'!'.OR.CH.EQ.'c'.OR.CH.EQ.'C') GOTO 10 IZ=0 11 IZ=IZ+1 C Read filenames READ(11,2) CH,SKU,CB,FILENAM IF(CH.EQ.'!'.OR.CH.EQ.'c'.OR.CH.EQ.'C') GOTO 11 IF(CH.EQ.'*') GOTO 13 IF(CH.EQ.' ') THEN C growth coefficient READ(11,3) ALF M=0 DO N=1,NSPDIM IF(SKU.EQ.ADU(N)) THEN M=N VALFA(N)=ALF END IF END DO IF(M.NE.0) GOTO 12 PRINT*,'--> Unknown species: ',SKU GOTO 11 12 IDU(M)=1 C Input optical data for specíes CALL READNK(M,FILENAM,RHODU,AWDU) VV0(M)=AWDU*PM/RHODU WRITE(6,4) VALFA(M) IF(VALFA(M).EQ.0D0) FLDU(M)=.FALSE. GOTO 11 END IF PRINT*,'--> Error in file >',DUSTFILE,'< in line',IZ WRITE(6,2) CH,SKU,CB,FILENAM STOP 13 CLOSE(11) C Check for missing data FL=0 DO N=1,NSPDIM IF(IDU(N).EQ.0) THEN PRINT*,'No data for species ',ADU(N),' supplied !!!' FL=1 END IF END DO IF(FL.NE.0) STOP DO N=1,NSPDIM IF(.NOT.FLDU(N)) PRINT*,'---> Species switched off: ',ADU(N) END DO C -- Reset some dust data -- DO N=1,NSPDIM IF(ADU(N).EQ.'Oli') THEN IF(VALFA(N).LT.1D0) ALFOLI=VALFA(N) V0OLI=VV0(N) END IF IF(ADU(N).EQ.'Pyr') THEN IF(VALFA(N).LT.1D0) ALFPYR=VALFA(N) V0PYR=VV0(N) END IF IF(ADU(N).EQ.'Car') THEN IF(VALFA(N).LT.1D0) ALFCAR=VALFA(N) V0CAR=VV0(N) END IF IF(ADU(N).EQ.'SiC') THEN IF(VALFA(N).LT.1D0) ALFSIC=VALFA(N) V0SIC=VV0(N) END IF IF(ADU(N).EQ.'Irn') THEN IF(VALFA(N).LT.1D0) ALFIRO=VALFA(N) V0IRO=VV0(N) END IF IF(ADU(N).EQ.'MgS') THEN IF(VALFA(N).LT.1D0) ALFMGS=VALFA(N) V0MGS=VV0(N) END IF IF(ADU(N).EQ.'Cor') THEN IF(VALFA(N).LT.1D0) ALFCOR=VALFA(N) V0COR=VV0(N) END IF IF(ADU(N).EQ.'Hib') THEN IF(VALFA(N).LT.1D0) ALFHIB=VALFA(N) V0HIB=VV0(N) END IF IF(ADU(N).EQ.'Spi') THEN IF(VALFA(N).LT.1D0) ALFSPI=VALFA(N) V0SPI=VV0(N) END IF IF(ADU(N).EQ.'Qtz') THEN IF(VALFA(N).LT.1D0) ALFQUA=VALFA(N) V0QUA=VV0(N) END IF IF(ADU(N).EQ.'Wue') THEN IF(VALFA(N).LT.1D0) ALFWUMG=VALFA(N) IF(VALFA(N).LT.1D0) ALFWUFE=ALFWUMG*DSQRT(AWFE/AWMG) V0WUE=VV0(N) END IF END DO C -- write extinction data -- OPEN(10,FILE='abs.da') WRITE(10,'(2X,3Hlam,7X,11(3HQa(,A3,1H),3X,3HQs(,A3,1H),3X))') * (ADU(J),ADU(J),J=1,NSPDIM) DO N=1,NN WRITE(10,1) DLA(N),(ABSDU0(N,I),STRDU0(N,I),I=1,NSPDIM) END DO CLOSE(10) RETURN 9010 PRINT*,'In DUSTABS: file >',FILENAM,'< does not exist' STOP END C***********************************************************************C C Read tables for species C C***********************************************************************C SUBROUTINE READNK(NSPDU,FILENAM,RHODU,AWDU) IMPLICIT REAL*8(A-H,O-Z) LOGICAL FLDU PARAMETER(NDIM=80,NSPDIM=11,NRDIM=7000) COMMON/DUEXT0/ABSDU0(NDIM,NSPDIM),STRDU0(NDIM,NSPDIM) COMMON/FREQ/DLA(NDIM),DNU(NDIM),GNU(NDIM),DLAMREF,NN,NREF CHARACTER*3 ADU COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) CHARACTER*1 TEXT CHARACTER*20 FILENAM DIMENSION TEXT(80),VL(NRDIM),QAB(NRDIM),QSC(NRDIM),VNU(NRDIM) DATA CL/2.99793D+10/ 1 FORMAT(1P20E10.3) 2 FORMAT(80A1) 3 FORMAT(24X,I5) 4 FORMAT(46X,F5.2) 5 FORMAT('Range:',F9.2,2X,'to',F7.3,2X,'mu') 6 FORMAT(7X,F6.6,18X,F9.3) WRITE(6,'(10H----------,/,8Hspezies:,1X,A3)') ADU(NSPDU) OPEN(10,ERR=9010,FILE=FILENAM,STATUS='OLD') READ(10,3) NZ DO I=1,NZ READ(10,2) (TEXT(J),J=1,80) IF(I.EQ.1) WRITE(6,2) (TEXT(J),J=1,80) END DO C Density, Molecular weight of dust material READ(10,6) RHODU,AWDU C particle radius READ(10,4) RDU C Number of wavelength grid-points READ(10,3) NLDU IF(NLDU.GE.NRDIM) THEN PRINT*,'Too many input data:',NLDU STOP END IF READ(10,2) (TEXT(J),J=1,80) C wavelength (in mu) and complex index of refraction DO I=1,NLDU READ(10,*) VL(I),DN,DK,QAB(I),QSC(I) END DO CLOSE(10) C Sorting: Wavelength decreasing with increasing index IF(VL(1).GT.VL(2)) GOTO 100 DO J=1,NLDU/2 H1=VL(J) H2=QAB(J) H3=QSC(J) VL(J)=VL(NLDU-J+1) QAB(J)=QAB(NLDU-J+1) QSC(J)=QSC(NLDU-J+1) VL(NLDU-J+1)=H1 QAB(NLDU-J+1)=H2 QSC(NLDU-J+1)=H3 END DO 100 WRITE(6,5) VL(1),VL(NLDU) DO I=1,NLDU VNU(I)=CL/(1D-4*VL(I)) END DO C Interpolation R0=1D-5 VDUGES=R0**3 DO N=1,NN CALL FIPOL(VNU,DNU(N),NLDU,IO,IL) G1=(VNU(IO)-DNU(N))/(VNU(IO)-VNU(IL)) G2=1D0-G1 QA=G1*QAB(IL)+G2*QAB(IO) QS=G1*QSC(IL)+G2*QSC(IO) C Conversion of units ABSDU0(N,NSPDU)=QA*1D4 STRDU0(N,NSPDU)=QS*1D16*VDUGES END DO RETURN 9010 PRINT*,'In READNK: file >',FILENAM,'< does not exist' STOP END C**********************************************************************C C Reset dust extinction for new abundances C C**********************************************************************C SUBROUTINE SETDUSTABS IMPLICIT REAL*8(A-H,O-Z) LOGICAL FLDU PARAMETER(KDIM=1000,NDIM=80,NSPDIM=11) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT COMMON/ABUND/EPSHE,EPSC,EPSN,EPSO,EPSNE,EPSNA,EPSMG,EPSAL,EPSSI, * EPSS,EPSAR,EPSCA,EPSFE,EPSNI COMMON/COLI/V0OLI,A0OLI,EPSOLI,ALFOLI,ALFOLEX COMMON/CPYR/V0PYR,A0PYR,EPSPYR,ALFPYR,ALFPYEX COMMON/CQUA/V0QUA,A0QUA,EPSQUA,ALFQUA COMMON/CIRO/V0IRO,A0IRO,EPSIRO,ALFIRO COMMON/CCAR/V0CAR,A0CAR,EPSCAR,ALFCAR COMMON/CSIC/V0SIC,A0SIC,EPSSIC,ALFSIC COMMON/CMGS/V0MGS,A0MGS,EPSMGS,ALFMGS COMMON/CWUE/V0WUE,A0WUE,EPSWUE,ALFWUMG,ALFWUFE COMMON/CCOR/V0COR,A0COR,EPSCOR,ALFCOR COMMON/CHIB/V0HIB,A0HIB,EPSHIB,ALFHIB COMMON/CSPI/V0SPI,A0SPI,EPSSPI,ALFSPI COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) COMMON/FREQ/DLA(NDIM),DNU(NDIM),GNU(NDIM),DLAMREF,NN,NREF COMMON/DUEXT/ABSDU(NDIM,NSPDIM),STRDU(NDIM,NSPDIM) COMMON/DUEXT0/ABSDU0(NDIM,NSPDIM),STRDU0(NDIM,NSPDIM) CHARACTER*3 ADU C -- Olivine -- DKOLIDUF=.75D0*EPSSI*V0OLI/(1.4D0*PM) C -- Pyroxene -- DKPYRDUF=.75D0*EPSSI*V0PYR/(1.4D0*PM) C -- Quartz -- DKQUADUF=.75D0*EPSSI*V0QUA/(1.4D0*PM) C -- Carbon -- DKCARDUF=.75D0*EPSC*V0CAR/(1.4D0*PM) C -- SiC -- DKSICDUF=.75D0*EPSSI*V0SIC/(1.4D0*PM) C -- iron -- DKIRODUF=.75D0*EPSFE*V0IRO/(1.4D0*PM) C -- MgS --- DKMGSDUF=.75D0*EPSS*V0MGS/(1.4D0*PM) C -- Wuestite --- DKWUEDUF=.75D0*EPSMG*V0WUE/(1.4D0*PM) C -- Corundum --- DKCORDUF=.75D0*EPSAL*V0COR/(1.4D0*PM) C -- Hibonite --- DKHIBDUF=.75D0*EPSCA*V0HIB/(1.4D0*PM) C -- Spinel --- DKSPIDUF=.75D0*EPSAL*V0SPI/(1.4D0*PM) C Multiplication of extinction coefficients with element abundance C dependent factors DO N=1,NN ABSDU(N,1)=DKOLIDUF*ABSDU0(N,1) STRDU(N,1)=DKOLIDUF*STRDU0(N,1) ABSDU(N,2)=DKPYRDUF*ABSDU0(N,2) STRDU(N,2)=DKPYRDUF*STRDU0(N,2) ABSDU(N,3)=DKQUADUF*ABSDU0(N,3) STRDU(N,3)=DKQUADUF*STRDU0(N,3) ABSDU(N,4)=DKCORDUF*ABSDU0(N,4) STRDU(N,4)=DKCORDUF*STRDU0(N,4) ABSDU(N,5)=DKHIBDUF*ABSDU0(N,5) STRDU(N,5)=DKHIBDUF*STRDU0(N,5) ABSDU(N,6)=DKSPIDUF*ABSDU0(N,6) STRDU(N,6)=DKSPIDUF*STRDU0(N,6) ABSDU(N,7)=DKIRODUF*ABSDU0(N,7) STRDU(N,7)=DKIRODUF*STRDU0(N,7) ABSDU(N,8)=DKWUEDUF*ABSDU0(N,8) STRDU(N,8)=DKWUEDUF*STRDU0(N,8) ABSDU(N,9)=DKCARDUF*ABSDU0(N,9) STRDU(N,9)=DKCARDUF*STRDU0(N,9) ABSDU(N,10)=DKSICDUF*ABSDU0(N,10) STRDU(N,10)=DKSICDUF*STRDU0(N,10) ABSDU(N,11)=DKMGSDUF*ABSDU0(N,11) STRDU(N,11)=DKMGSDUF*STRDU0(N,11) END DO RETURN END C**********************************************************************C C Intializing the progam module for calculation of spectra C C**********************************************************************C SUBROUTINE INIRAD(DUSTFILE) IMPLICIT REAL*8 (A-H,O-Z) PARAMETER(KDIM=1000,NDIM=80,IPDIM=13,IFDIM=11,NSPDIM=11) COMMON/FREQ/DLA(NDIM),DNU(NDIM),GNU(NDIM),DLAMREF,NN,NREF COMMON/DUEXT/ABSDU(NDIM,NSPDIM),STRDU(NDIM,NSPDIM) COMMON/CSRP/KPS,KPC COMMON/CTIT/ERRMH,ERRMT,ERRMJ,ITTMAX,IWH CHARACTER*20 DUSTFILE C Reference wavelength (in mu) DLAMREF=1.0D0 C make frequency scale CALL FREQSCALE(DUSTFILE) C input dust extinction data and make extinktion coefficients CALL DUSTABS(DUSTFILE) C make extinction coefficients for given abundances c CALL SETDUSTABS C -- outer boundary condition of the radiation field -- C MISRF=0: No external radiation field C =1: Diluted black-body radiation C =2: Model of the interstellar radiation field MISRF=2 CALL ISRF(MISRF) C --Parameter for the solution of the radiative transfer problem -- C Number of rays hitting the star KPS=8 C Number of rays between star and dust shell KPC=15 C -- Shows information on iteration if IWH<>0 -- IWH=1 C -- Parameter for temperature iteration -- C Maximum error for flux ERRMH=.010D0 C Maximum error for temperature ERRMT=3D-5 C Maximum error for scattering ERRMJ=3D-3 C Maximum number of iteration steps for flux iteration ITTMAX=2000 RETURN END C**********************************************************************C C Interpolation in ascending order C C**********************************************************************C SUBROUTINE IPVEC(KLU,XA,YA,XN,YN,KA,KN) IMPLICIT REAL*8(A-H,O-Z) DIMENSION XA(1),YA(1),XN(1),YN(1) C upper limit: old: KA new: KN K=KLU KL=KLU 100 IF(XN(K).LE.XA(KLU)) THEN YN(K)=YA(KLU) K=K+1 GOTO 100 END IF 200 IF(XA(KL).LT.XN(K)) THEN IF(KL.LT.KA) THEN KL=KL+1 GOTO 200 ELSE GOTO 300 END IF END IF YN(K)=YA(KL-1)+(YA(KL)-YA(KL-1))*(XN(K)-XA(KL-1))/ * (XA(KL)-XA(KL-1)) IF(K.LT.KN) THEN K=K+1 GOTO 200 END IF 300 IF(XN(K).GE.XA(KA)) THEN YN(K)=YA(KA) IF(K.LT.KN) THEN K=K+1 GOTO 300 END IF END IF RETURN END