C***********************************************************************C C***********************************************************************C C C C AA GGGGG BBBBBB DDDDDD UU UU SSSSS TTTTTTTT C C AAA GG GG BB BB DD DD UU UU SS SS T TT T C C AAAA GG BB BB DD DD UU UU SS TT C C AA AA GG GGG BBBBBB DD DD UU UU SSSSS TT C C AAAAAA GG G BB BB DD DD UU UU SS TT C C AA AA GG GG BB BB DD DD UU UU SS SS TT C C AAAA AAAA GGGGG BBBBBB DDDDDD UUUU SSSSS TTTT C C C C***********************************************************************C C***********************************************************************C C Subroutine package for calculating dust condensation in a simple C C wind model for AGB stars of all C/O abundance ratios 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 Master routine for calculating a wind model with dust formation C C***********************************************************************C SUBROUTINE ITTEMPL(STARM,STARL,TEFF,DMP,ROUT,VEXP,TAULUC,FOL,XOL, * FPY,XPY,FQU,FIR,FWU,XWU,FCO,FHI,FSP,FCA,FSC,FMS,VWIND,TAUH, * TAULFIN) PARAMETER(KDIM=1000,IFDIM=11,IPDIM=21) IMPLICIT REAL*8(A-H,O-Z) EXTERNAL MODELLO,MODELLC,MODELLS CHARACTER*20 DUSTFILE 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/RMOD/RD(KDIM),TG(KDIM),RHO(KDIM),VV(KDIM),DKAP(KDIM), * GAM(KDIM),FV(KDIM,IFDIM),VDR(KDIM,IFDIM),AV(KDIM,IFDIM),KL,KK DIMENSION PARM(IPDIM) COMMON/CTE/CHM,RCOS,RCOC COMMON/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM IF((EPSC+EPSO+EPSSI+EPSMG+EPSFE).GT.1D-2.OR.EPSHE.GT.0.5D0) THEN PRINT*,'Only hydrogen rich mixtures admitted' STOP END IF C -- input parameters for wind model -- PARM(1)=STARM*SUNM PARM(2)=STARL*SUNL PARM(3)=TEFF PARM(4)=VEXP*1D5 PARM(5)=DMP*SUNM/YEAR PARM(6)=TAULUC PARM(7)=ROUT RCO=EPSC/EPSO C -- M-stars -- IF(RCO.LE.RCOS) THEN C start model CALL ITTEMPLA(PARM,NDUSPO,MODELLO) C model iteration CALL ITTEMPRT(PARM,MODELLO) C -- return values -- FOL=PARM(NDUAM+1) XOL=PARM(NDUAM+2) FPY=PARM(NDUAM+3) XPY=PARM(NDUAM+4) FQU=PARM(NDUAM+5) FIR=PARM(NDUAM+6) FWU=PARM(NDUAM+7) XWU=PARM(NDUAM+8) FCO=PARM(NDUAM+9) FHI=PARM(NDUAM+10) FSP=PARM(NDUAM+11) VWIND=PARM(NDUSPO+8) TAUH=PARM(NDUSPO+9) TAULFIN=PARM(NDUSPO+10) FCA=0D0 FSC=0D0 FMS=0D0 RETURN END IF C -- S-stars -- IF(RCO.LE.RCOC) THEN CALL ITTEMPLA(PARM,NDUSPS,MODELLS) C model iteration CALL ITTEMPRT(PARM,MODELLS) C -- return values -- FIR=PARM(NDUAM+1) FQU=PARM(NDUAM+2) VWIND=PARM(NDUSPS+8) TAUH=PARM(NDUSPS+9) TAULFIN=PARM(NDUSPS+10) FOL=0D0 XOL=0D0 FPY=0D0 XPY=0D0 FWU=0D0 XWU=0D0 FCO=0D0 FSP=0D0 FCA=0D0 FSC=0D0 FMS=0D0 RETURN END IF C -- C-stars -- C start model CALL ITTEMPLA(PARM,NDUSPC,MODELLC) C model iteration CALL ITTEMPRT(PARM,MODELLC) C -- return values -- FCA=PARM(NDUAM+1) FSC=PARM(NDUAM+2) FIR=PARM(NDUAM+3) FMS=PARM(NDUAM+4) VWIND=PARM(NDUSPC+8) TAUH=PARM(NDUSPC+9) TAULFIN=PARM(NDUSPC+10) FOL=0D0 XOL=0D0 FPY=0D0 XPY=0D0 FQU=0D0 FWU=0D0 XWU=0D0 FCO=0D0 FSP=0D0 RETURN END C***********************************************************************C C Iteration of temperature stratification (Lucy-approximation) C C***********************************************************************C SUBROUTINE ITTEMPLA(PARM,NDUSPOC,DUSTFORM) PARAMETER(IPDIM=21) IMPLICIT REAL*8(A-H,O-Z) EXTERNAL DUSTFORM COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT COMMON/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM DIMENSION PARM(IPDIM) LOGICAL NEWMO 10 FORMAT(1X,'it:',I3,3X,'tauLi:',1PE11.4,3X,'tauLf:',E11.4,3X, * 'L:',E11.4,3X,'Teff:',E11.4,3X,'M:',E11.4,3X,'Mp:',E11.4) C -- initial values -- NEWMO=.TRUE. TAULH=PARM(6) CALL INITEM(PARM) TU=TAULH TO=TAULH TM=TAULH C -- bisection method to determine TAUL -- IM=1 IT=0 100 IT=IT+1 IF(IT.GT.100) GOTO 200 PARM(6)=TM CALL DUSTFORM(PARM,NEWMO) TAUL=PARM(NDUSPOC+10) IF(IT.GT.90) WRITE(6,10) IT,TM,TAUL,PARM(2)/SUNL, * PARM(3),PARM(1)/SUNM,PARM(5)*YEAR/SUNM IF(IM.NE.1) GOTO 110 C -- right boundary too small -- IF(TAUL.LT.0D0) THEN TU=TM TO=1.5D0*TO TM=TO ELSE IM=2 TM=TU END IF GOTO 100 110 IF(IM.NE.2) GOTO 120 C -- left boundary too big -- IF(TAUL.GE.0D0) THEN TU=.1D0*TU TM=TU ELSE IM=3 TM=.1D0*(TU+TO) END IF GOTO 100 120 IF((TO-TU)/TU.GT.1D-4.OR.TAUL.GT.1D-8) THEN C -- bisection -- IF(TAUL.GT.0D0) THEN TO=TM ELSE TU=TM END IF TM=.5D0*(TU+TO) GOTO 100 END IF C -- method successful -- PARM(NDUSPOC+10)=TM PARM(6)=TAULH RETURN 200 IF(DABS(TAUL).LT.1D-10) RETURN PRINT *,'Iteration of TAUL in ITTEMPL did not converge' STOP END C***********************************************************************C C Iteration of temperature stratification, solving radiative transfer C C***********************************************************************C SUBROUTINE ITTEMPRT(PARM,DUSTFORM) PARAMETER(KDIM=1000,NSPDIM=11,IFDIM=11,IPDIM=21) IMPLICIT REAL*8(A-H,O-Z) EXTERNAL DUSTFORM LOGICAL FLDU 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/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/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM COMMON/CSTAR/DLSTERN,RSTERN,TEFF COMMON/CHEMTYPE/RCO COMMON/CTIT/ERRMH,ERRMT,ERRMJ,ITTMAX,IWH COMMON/CTEMP/TS(KDIM,NSPDIM),TGH(KDIM) COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) COMMON/CFDU/FDU(KDIM,NSPDIM) COMMON/CMOD/TAU(KDIM),TAL(KDIM),TAH(KDIM),CHIHT(KDIM),KLU DIMENSION TSA(KDIM,NSPDIM),TGHA(KDIM),FVA(KDIM,NSPDIM), * RHOA(KDIM),CHIHA(KDIM),VVA(KDIM) DIMENSION PARM(IPDIM) LOGICAL NEWMO CHARACTER*3 ADU 10 FORMAT(1X,'it:',I3,3X,'tauLi:',1PE11.4,3X,'tauLf:',E11.4,3X, * 'L:',E11.4,3X,'Teff:',E11.4,3X,'M:',E11.4,3X,'Mp:',E11.4) C -- Stellar parameter -- DLSTERN=PARM(2) TEFF=PARM(3) RSTERN=DSQRT(DLSTERN/(4D0*PI*SIG*TEFF**4)) RCO=EPSC/EPSO C -- initial values -- NEWMO=.FALSE. CALL SETDUSTABS CALL CHECKDU ITR=ITTMAX ITTMAX=1 ITT=0 IT=0 100 IT=IT+1 IF(IT.GT.ITR) GOTO 200 WRITE(6,'(3Hit:,I5)') IT DO K=1,KK DO L=1,NDUSPO FVA(K,L)=FV(K,L) END DO RHOA(K)=RHO(K) VVA(K)=VV(K) END DO IF(IT.LT.40) GOTO 110 IF(IT.LE.200.AND.MOD(IT,10).EQ.0) GOTO 110 IF(IT.GT.400.AND.MOD(IT,20).EQ.0) GOTO 110 GOTO 120 110 CALL DUSTFORM(PARM,NEWMO) 120 FW=0.5D0 FR=0.9D0 FF=0.1D0 IF(IT.GT.100) FF=0.05D0 IF(IT.GE.600) FF=0.02D0 IF(IT.GE.600) FW=0.1D0 DO K=1,KK DO L=1,NDUSPO IF(FV(K,L).NE.0) FV(K,L)=FF*FV(K,L)+(1D0-FF)*FVA(K,L) END DO RHO(K)=FW*RHO(K)+(1D0-FW)*RHOA(K) END DO CALL CHECKDU IF(IT.EQ.1) THEN CALL STARTEM CALL RADMODIT(-1,IERR) END IF DO K=KK,KLU,-1 DO L1=1,NDON L=NDONP(L1) TSA(K,L)=TS(K,L) END DO TGHA(K)=TGH(K) CHIHA(K)=CHIHT(K) END DO CALL RADMOD(HKORR,ERRJ,EH,KH) C -- temperature correction -- CALL KORRTEMP(ERRT,KT,KS) DO K=KK,KLU,-1 DO L1=1,NDON L=NDONP(L1) TS(K,L)=FR*TS(K,L)+(1D0-FR)*TSA(K,L) END DO TGH(K)=FR*TGH(K)+(1D0-FR)*TGHA(K) CHIHT(K)=FR*CHIHT(K)+(1D0-FR)*CHIHA(K) VV(K)=FR*VV(K)+(1D0-FR)*VVA(K) END DO CALL MKTAU(0) PARM(6)=0D0 WRITE(6,'(11X,5Hvend:,F7.3,5H km/s,3X,4Htau:,F8.3)') * VV(KK)*1D-5,TAU(KK) IF(DABS(ERRT).GT.ERRMT.OR.DABS(ERRJ).GT.ERRMJ) THEN IETC=0 ELSE IETC=IETC+1 END IF IF(IETC.LT.3) GOTO 100 RETURN 200 PRINT *,'Iteration of Model in ITTEMPRT did not converge' RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C S u b r o u t i n e s f o r M - s t a r s C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C***********************************************************************C C Integration of wind equations (driver) C C for oxygen rich element mixture C C***********************************************************************C SUBROUTINE MODELLO(PARM,NEWMO) PARAMETER(KDIM=1000,IPDIM=21,IYDIM=IPDIM-4,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) 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/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/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/CCOR/V0COR,A0COR,EPSCOR,ALFCOR COMMON/CHIB/V0HIB,A0HIB,EPSHIB,ALFHIB COMMON/CSPI/V0SPI,A0SPI,EPSSPI,ALFSPI COMMON/CWUE/V0WUE,A0WUE,EPSWUE,ALFWUMG,ALFWUFE COMMON/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM DIMENSION YN(IYDIM),YA(IYDIM),DYDTN(IYDIM-1),DYDTA(IYDIM-1), * PARM(IPDIM) LOGICAL NEWMO 1 FORMAT(I8,1P20E10.3) C -- outer boundary -- STARR=DSQRT(PARM(2)/(4D0*PI*SIG*PARM(3)**4)) RMAX=PARM(7)*STARR C -- initial values -- RMIN=2.5D0*STARR CALL FIPOL(RD,RMIN,KK,KO,KL) YA(1)=RMIN YA(2)=PARM(4) YA(3)=PARM(6) YA(4)=0D0 YA(5)=A0OLI YA(6)=0.8D0 YA(7)=A0PYR YA(8)=0.8D0 YA(9)=A0QUA YA(10)=A0IRO YA(11)=A0WUE YA(12)=1D-4 YA(13)=A0COR YA(14)=A0HIB YA(15)=A0SPI IF(NEWMO) THEN CALL STOREO(YN,YA,PARM,-KL) END IF CALL DGLSO(DT,YA,DYDTA,PARM) C -- loop for time integration -- DT=1D5 DTA=DT KNEXT=KO ISTEP=0 1000 ISTEP=ISTEP+1 1010 CALL NEXTTO(DT,DTA,DTN,YN,YA,DYDTN,DYDTA,PARM) IF(MOD(ISTEP,10000).EQ.0) WRITE(6,1) ISTEP,DT,YN IF(DTN.LT.0.5D0*DT) THEN DT=DTN c PRINT *,'repeat last step:',ISTEP GOTO 1010 END IF CALL STOREO(YN,YA,PARM,KNEXT) C -- new time step-width -- DTA=DT IF(DTN.GT.1.2D0*DT) DT=1.2D0*DT IF(DTN.LT.DT) DT=DTN C -- save last step -- DO N=1,NDUSPO+6 YA(N)=YN(N) END DO DO N=1,NDUSPO+4 DYDTA(N)=DYDTN(N) END DO IF(YN(1).LT.RMAX) GOTO 1000 C -- store final values -- PARM(NDUAM+1)=4D0*PI*(YN(5)**3-A0OLI**3)*EPSOLI/(3D0*V0OLI*EPSSI) PARM(NDUAM+2)=YN(6) PARM(NDUAM+3)=4D0*PI*(YN(7)**3-A0PYR**3)*EPSPYR/(3D0*V0PYR*EPSSI) PARM(NDUAM+4)=YN(8) PARM(NDUAM+5)=4D0*PI*(YN(9)**3-A0QUA**3)*EPSQUA/(3D0*V0QUA*EPSSI) PARM(NDUAM+6)=4D0*PI*(YN(10)**3-A0IRO**3)*EPSIRO/(3D0*V0IRO*EPSFE) PARM(NDUAM+7)= (YN(11)**3-A0WUE**3)*EPSWUE/(V0WUE*EPSMG) PARM(NDUAM+8)=YN(12) PARM(NDUAM+9)=4D0*PI*(YN(13)**3-A0COR**3)*EPSCOR/(3D0*V0COR*EPSAL) PARM(NDUAM+10)=4D0*PI*(YN(14)**3-A0HIB**3)*EPSHIB/ * (3D0*V0HIB*EPSAL) PARM(NDUAM+11)=4D0*PI*(YN(15)**3-A0SPI**3)*EPSSPI/ * (3D0*V0SPI*EPSAL) PARM(NDUSPO+8)=YN(2)*1D-5 PARM(NDUSPO+9)=YN(4) PARM(NDUSPO+10)=YN(3) RETURN END C***********************************************************************C C Integration of wind equations (Stepper) C C Method: Adams - Bashforth C C for oxygen rich element mixtures C C***********************************************************************C SUBROUTINE NEXTTO(DT,DTA,DTN,YN,YA,DYDTN,DYDTA,PARM) PARAMETER(IPDIM=21,IYDIM=IPDIM-4) IMPLICIT REAL*8(A-H,O-Z) DIMENSION YN(IYDIM),YA(IYDIM),DYDTN(IYDIM-1),DYDTA(IYDIM-1), * PARM(IPDIM) COMMON/CTE/CHM,RCOS,RCOC COMMON/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM CALL DGLSO(DT,YA,DYDTN,PARM) YN(NDUSPO+5)=YA(NDUSPO+5) YN(NDUSPO+6)=YA(NDUSPO+6) DO N=1,NDUSPO+4 YN(N)=YA(N)+DT*(DYDTN(N)+.5D0*(DYDTN(N)-DYDTA(N))*DT/DTA) END DO DTN=1D30 DO N=1,NDUSPO+4 X=CHM*DABS(YN(N)/DYDTN(N)) IF((N.EQ.3).OR.(N.EQ.4)) X=1D30 IF((N.EQ.6).OR.(N.EQ.8).OR.(N.EQ.12)) THEN IF(DYDTN(N).GT.1D-30) THEN X=CHM*DABS((1D0-YN(N))/DYDTN(N)) ELSE X=1D30 END IF END IF IF(X.LT.DTN) DTN=X END DO RETURN END C***********************************************************************C C Right hand sides of the DEQ-System C C for oxygen rich element mixture C C***********************************************************************C C ---- Independent variables C Y(1) : r C Y(2) : v C Y(3) : tauL C Y(4) : tau* C Y(5) : aOli C Y(6) : xOli C Y(7) : aPyr C Y(8) : xPyr C Y(9) : aQua C Y(10): aFe C Y(11): aWu C Y(12): xWu C Y(13): aCo C Y(14): aHi C Y(15): aSp C Y(16): Tg C Y(17): Gamma C ---- Parameters C PARM(1) =STARM C PARM(2) =STARL C PARM(3) =TEFF C PARM(4) =VEXP C PARM(5) =DMP C PARM(6) =tauLini C PARM(7) =Rout C PARM(8) =FOL C PARM(9) =XOL C PARM(10)=FPY C PARM(11)=XPY C PARM(12)=FQU C PARM(13)=FIR C PARM(14)=FWU C PARM(15)=XWU C PARM(16)=FCO C PARM(17)=FHI C PARM(18)=FSP C PARM(19)=VWIND C PARM(20)=kapH C PARM(21)=tauLfin C-----------------------------------------------------------------------C SUBROUTINE DGLSO(DT,Y,DYDT,PARM) PARAMETER(IPDIM=21,IYDIM=IPDIM-4) IMPLICIT REAL*8(A-H,O-Z) 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/STERN/STARL,STARM,STARR,TEFF COMMON/CAW/AWSIO,AWMG,AWFE,AWH2O,AWC2H2,AWSI,AWSIS,AWH2S,AWAL, * AWFEO,AWMGO,AWCA 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/CCOR/V0COR,A0COR,EPSCOR,ALFCOR COMMON/CHIB/V0HIB,A0HIB,EPSHIB,ALFHIB COMMON/CSPI/V0SPI,A0SPI,EPSSPI,ALFSPI COMMON/CWUE/V0WUE,A0WUE,EPSWUE,ALFWUMG,ALFWUFE COMMON/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM COMMON/CVDR/VDROLI,VDRPYR,VDRQUA,VDRIRO,VDRWUE,VDRCOR,VDRHIB, * VDRSPI,VDRCAR,VDRSIC,VDRMGS DIMENSION Y(IYDIM),DYDT(IYDIM-1),PARM(IPDIM) VDROLI=0 VDRPYR=0 VDRQUA=0 VDRIRO=0 VDRWUE=0 VDRCOR=0 VDRHIB=0 VDRSPI=0 C -- wind model -- R=Y(1) V=Y(2) STARM=PARM(1) STARL=PARM(2) TEFF=PARM(3) STARR=DSQRT(STARL/(4D0*PI*SIG*TEFF**4)) IF(PARM(6).LE.0D0) THEN c TG=TEMP(R) CALL TEMPO(R,TG,TOLI,TPYR,TQUA,TIRO,TWUE,TCOR,THIB,TSPI) C ungenügende Daten TQUA=TG THIB=TG TSPI=TG TWUE=TG ELSE TAU=DMAX1(Y(3),0D0) TG=PARM(3)* * DSQRT(DSQRT(.5D0*(1D0-DSQRT(1D0-(STARR/R)**2)+1.5D0*TAU))) TOLI=TG TPYR=TG TQUA=TG TIRO=TG TWUE=TG TCOR=TG THIB=TG TSPI=TG END IF Y(NDUSPO+5)=TG RHO=PARM(5)/(4D0*PI*V*R**2) CGAS2=BK*TG/(2.333D0*PM) C -- degrees of condensation (spheres, only wuestite grains are cubes) -- FOL=4D0*PI*(Y(5)**3-A0OLI**3)*EPSOLI/(3D0*V0OLI*EPSSI) FPY=4D0*PI*(Y(7)**3-A0PYR**3)*EPSPYR/(3D0*V0PYR*EPSSI) FQU=4D0*PI*(Y(9)**3-A0QUA**3)*EPSQUA/(3D0*V0QUA*EPSSI) FIR=4D0*PI*(Y(10)**3-A0IRO**3)*EPSIRO/(3D0*V0IRO*EPSFE) FWU= (Y(11)**3-A0WUE**3)*EPSWUE/(V0WUE*EPSMG) FCO=4D0*PI*(Y(13)**3-A0COR**3)*EPSCOR/(3D0*V0COR*EPSAL) FHI=4D0*PI*(Y(14)**3-A0HIB**3)*EPSHIB/(3D0*V0HIB*EPSAL) FSP=4D0*PI*(Y(15)**3-A0SPI**3)*EPSSPI/(3D0*V0SPI*EPSAL) XOL=Y(6) XPY=Y(8) XWU=Y(12) C -- radius -- DYDT(1)=V C -- kappaH and Q/a -- IF(PARM(6).LE.0D0) THEN CALL EXTDO(R,DKH,QAOLI,QAPYR,QAQUA,QAIRO,QAWUE,QACOR,QAHIB, * QASPI) ELSE TEFF=PARM(3) TRAD=((TEFF**4-TG**4)*DEXP(-Y(4))+TG**4)**.25D0 DKH=DKSIL(TRAD)*(FOL+FPY)!DKAMOL(TRAD)*FOL+DKAMPY(TRAD)*FPY * +DKQUA(TRAD)*FQU+DKIR(TRAD)*FIR+DKGAS(TG,RHO) * +PKCOR(TRAD)*(FCO+FSP) QAOLI=0D0 QAPYR=0D0 QAQUA=0D0 QAIRO=0D0 QAWUE=0D0 QACOR=0D0 QAHIB=0D0 QASPI=0D0 END IF PARM(20)=DKH C -- velocity -- GAMMA=STARL*DKH/(4D0*PI*CL*GRAV*STARM) DYDT(2)=DMAX1(-GRAV*STARM*(1D0-GAMMA)/R**2,1D-30) Y(NDUSPO+6)=GAMMA C -- optical depths -- DYDT(3)=-V*RHO*DKH*(STARR/R)**2 DYDT(4)=V*RHO*DKH C -- molecules -- PHG=RHO*BK*TG/((1D0+4D0*EPSHE)*PM) PH2=.5D0*PHG PH2O=(EPSO-EPSC-(1D0+3D0*FOL+2D0*FPY+FQU)*EPSSI- * (1.5D0*FCO+2D0*FSP+(19D0/6D0)*FHI)*EPSAL-FWU*EPSMG)*PHG PSIO=DMAX1((1D0-FOL-FPY-FQU)*EPSSI*PHG,0D0) PFE=DMAX1(((1D0-FIR)*EPSFE-(2D0*(1D0-XOL)*FOL+(1D0-XPY)*FPY)*EPSSI * -XWU*FWU*EPSMG)*PHG,0D0) PMG=DMAX1((EPSMG*(1D0-(1D0-XWU)*FWU)-(2D0*XOL*FOL+XPY*FPY)*EPSSI * -0.5D0*FSP*EPSAL)*PHG,0D0) PAL=DMAX1((1D0-FCO-FHI-FSP)*EPSAL*PHG,0D0) XGAS=PMG/(PMG+PFE) C -- drift velocities -- CALL VDUIP(R,VDROLI,VDRPYR,VDRQUA,VDRIRO,VDRWUE,VDRCOR,VDRHIB, * VDRSPI) C -- olivine -- DAGRSIO=ALFOLI*PSIO/DSQRT(2D0*PI*AWSIO*PM*BK*TG) DAGRMG=ALFOLI*PMG/DSQRT(2D0*PI*AWMG*PM*BK*TG) DAGRFE=ALFOLI*PFE/DSQRT(2D0*PI*AWFE*PM*BK*TG) DAGRH2O=ALFOLI*PH2O/(3D0*DSQRT(2D0*PI*AWH2O*PM*BK*TG))*.25D0 DAGROL=DMIN1(DAGRSIO,.5D0*(DAGRMG+DAGRFE),DAGRH2O)* * PHIMOLI(R,Y(5),TG,RHO,QAOLI,VDROLI) IF(TOLI.GT.650D0) THEN CALL KPOLVAP(TOLI,XOL,DKPOL) POLSIOEQ=(PH2*PAT/PH2O)**3/((PMG**XOL*PFE**(1D0-XOL))**2*DKPOL) DAVASIO=ALFOLI*POLSIOEQ/DSQRT(2D0*PI*AWSIO*PM*BK*TG) AOL=DAGROL/DAVASIO C growth for AOL>1 IF(AOL.GT.1D0) THEN DYDT(5)=V0OLI*(DAGROL-DAVASIO) C average x IF(FOL.GT.1D-7.AND.XOL.LT.0.999D0) THEN CALL OLIKP(TOLI,XOL,DKXOL) DYDT(6)=(3D0*V0OLI/Y(5))*((XGAS-XOL)*DAGROL+.5D0*ALFOLEX* * (PMG-PFE*DKXOL)/DSQRT(2D0*PI*AWMG*PM*BK*TG)) ELSE DYDT(6)=0D0 END IF ELSE C Olivine does not exist for AOL<1 DYDT(5)=0D0 DYDT(6)=0D0 END IF ELSE C growth at low temperature DYDT(5)=V0OLI*DAGROL DYDT(6)=(3D0*V0OLI/Y(5))*(XGAS-XOL)*DAGROL ENDIF C -- pyroxene -- DAGRSIO=ALFPYR*PSIO/DSQRT(2D0*PI*AWSIO*PM*BK*TG) DAGRMG=ALFPYR*PMG/DSQRT(2D0*PI*AWMG*PM*BK*TG) DAGRFE=ALFPYR*PFE/DSQRT(2D0*PI*AWFE*PM*BK*TG) DAGRH2O=ALFPYR*PH2O/(2D0*DSQRT(2D0*PI*AWH2O*PM*BK*TG))*(1D0/3D0) DAGRPY=DMIN1(DAGRSIO,DAGRMG+DAGRFE,DAGRH2O)* * PHIMPYR(R,Y(7),TG,RHO,QAPYR,VDRPYR) IF(TPYR.GT.650D0) THEN CALL KPPYVAP(TPYR,XPY,DKPPY) PPYSIOEQ=(PH2*PAT/PH2O)**2/(PMG**XPY*PFE**(1D0-XPY)*DKPPY) DAVASIO=ALFPYR*PPYSIOEQ/DSQRT(2D0*PI*AWSIO*PM*BK*TG) APY=DAGRPY/DAVASIO C growth for APY>1 IF(APY.GT.1D0) THEN DYDT(7)=V0PYR*(DAGRPY-DAVASIO) C average x IF(FPY.GT.1D-7.AND.XPY.LT.0.999D0) THEN CALL PYRKP(TPYR,XPY,DKXPY) DYDT(8)=(3D0*V0PYR/Y(7))*((XGAS-XPY)*DAGRPY * +ALFPYEX*(PMG-PFE*DKXPY)/DSQRT(2D0*PI*AWMG*PM*BK*TG)) ELSE DYDT(8)=0D0 END IF ELSE C Pyroxene does not exist for APY<1 DYDT(7)=0D0 DYDT(8)=0D0 ENDIF ELSE C growth at low temperature DYDT(7)=V0PYR*DAGRPY DYDT(8)=(3D0*V0PYR/Y(7))*(XGAS-XPY)*DAGRPY ENDIF C -- quartz -- GF=ALFQUA/DSQRT(2D0*PI*AWSIO*PM*BK*TG) PGRH2O=PH2O*DSQRT(AWSIO/AWH2O) PGRQU=DMIN1(PSIO,PGRH2O)*PHIMQUA(R,Y(9),TG,RHO,QUQUA,VDRQUA) IF(TG.GT.500D0) THEN PQUSIOEQ=PAT*PH2/(PH2O*DKPQUSIO(TG)) ELSE PQUSIOEQ=0D0 END IF DYDT(9)=V0QUA*GF*DMAX1(PGRQU-PQUSIOEQ,0D0) C -- iron -- GF=ALFIRO/DSQRT(2D0*PI*AWFE*PM*BK*TG) PGR=PFE*PHIMIRN(R,Y(10),TG,RHO,QAIRO,VDRIRO) IF(TIRO.GT.500D0) THEN PFEEQ=PVFE(TIRO) ELSE PFEEQ=0D0 ENDIF DYDT(10)=V0IRO*GF*DMAX1(PGR-PFEEQ,0D0) C -- Wuestite -- IF(TWUE.GT.500D0) THEN PWUFEEQ=PAT*(PH2/PH2O)*DKPWUFEO(TWUE,XWU) PWUMGEQ=PAT*(PH2/PH2O)*DKPWUMGO(TWUE,XWU) ELSE PWUFEEQ=0.D0 PWUMGEQ=0.D0 END IF F1=PMG*PHIMWUE(R,Y(11),TG,RHO,QAWUE,VDRWUE)-PWUMGEQ IF(F1.GT.0D0) THEN GRMGO=F1*ALFWUMG/DSQRT(2D0*PI*AWMG*PM*BK*TG) DYDT(11)=2D0*V0WUE*GRMGO F2=PFE*PHIMWUE(R,Y(11),TG,RHO,QAWUE,VDRWUE)-PWUFEEQ IF(F2.GT.0D0) THEN GRFEO=F2*ALFWUFE/DSQRT(2D0*PI*AWFE*PM*BK*TG) DYDT(11)=DYDT(11)+2D0*V0WUE*GRFEO DYDT(12)=6D0*V0WUE/Y(11)*((1D0-XWU)*GRFEO-XWU*GRMGO) ELSE DYDT(12)=0D0 END IF ELSE DYDT(11)=0D0 DYDT(12)=0D0 END IF C -- Corundum -- GF=ALFCOR/DSQRT(2D0*PI*AWAL*PM*BK*TG) DAGRAL=GF*PAL*.5D0 IF(TCOR.GT.650D0) THEN PCORALEQ=PAT*(PH2/PH2O)**1.5D0/DSQRT(DKPCOAL(TCOR)) ELSE PCORALEQ=0D0 END IF DAVAAL=PCORALEQ*GF DYDT(13)=V0COR*DMAX1(DAGRAL*PHIMCOR(R,Y(13),TG,RHO,QACOR,VDRCOR) * -DAVAAL,0D0) C -- Hibonite -- GF=ALFHIB/DSQRT(2D0*PI*AWAL*PM*BK*TG) DAGRAL=GF*PAL*(1D0/6D0) IF(THIB.GT.650D0) THEN PHIBALEQ=PAT*(PH2/PH2O)**1.5D0/DSQRT(DKPCOAL(THIB)) ELSE PHIBALEQ=0D0 END IF DAVAAL=PHIBALEQ*GF DYDT(14)=V0HIB*DMAX1(DAGRAL*PHIMHIB(R,Y(14),TG,RHO,QAHIB,VDRHIB) * -DAVAAL,0D0) C -- Spinel -- GF=ALFSPI/DSQRT(2D0*PI*AWAL*PM*BK*TG) DAGRAL=GF*PAL*.5D0 IF(TSPI.GT.650D0) THEN PSPIALEQ=PAT**1.5D0*(PH2/PH2O)**2/DSQRT(PMG*DKPSPAL(TSPI)) ELSE PSPIALEQ=0D0 END IF DAVAAL=PSPIALEQ*GF DYDT(15)=V0SPI*DMAX1(DAGRAL*PHIMSPI(R,Y(15),TG,RHO,QASPI,VDRSPI) * -DAVAAL,0D0) RETURN END C***********************************************************************C C Correction factors for dust drift C C***********************************************************************C DOUBLE PRECISION FUNCTION PHIMOLI(RR,A,TG,RHO,QAOLI,VDR) PARAMETER(KDIM=1000,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) DATA RHOC/3.3D0/,AWSIO/44.09D0/ C Q_ext/a Silicate dust, Draine, at 2 mu DATA QSIL/2.38D3/ IF(QAOLI.EQ.0D0) QAOLI=QSIL CALL VDUST(RHOC,QAOLI,RR,A,TG,RHO,VDR,PHI) PHIMOLI=PHI RETURN END DOUBLE PRECISION FUNCTION PHIMPYR(RR,A,TG,RHO,QAPYR,VDR) PARAMETER(KDIM=1000,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) DATA RHOC/3.3D0/,AWSIO/44.09D0/ C Q_ext/a Silicate dust, Draine, at 2 mu DATA QSIL/2.38D3/ IF(QAPYR.EQ.0D0) QAPYR=QSIL CALL VDUST(RHOC,QAPYR,RR,A,TG,RHO,VDR,PHI) PHIMPYR=PHI RETURN END DOUBLE PRECISION FUNCTION PHIMQUA(RR,A,TG,RHO,QAQUA,VDR) PARAMETER(KDIM=1000,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) DATA RHOC/3.3D0/,AWSIO/44.09D0/ C Q_ext/a Quartz, at 2 mu DATA QSIL/2.40D3/ IF(QAQUA.EQ.0D0) QAQUA=QSIL CALL VDUST(RHOC,QAQUA,RR,A,TG,RHO,VDR,PHI) PHIMQUA=PHI RETURN END DOUBLE PRECISION FUNCTION PHIMIRN(RR,A,TG,RHO,QAIRO,VDR) PARAMETER(KDIM=1000,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) DATA RHOC/7.8D0/,AWFE/55.847D0/ C Q_ext/a Iron dust, at 2 mu DATA QIRN/2.40D4/ IF(QAIRO.EQ.0D0) QAIRO=QIRN CALL VDUST(RHOC,QAIRO,RR,A,TG,RHO,VDR,PHI) PHIMIRN=PHI RETURN END DOUBLE PRECISION FUNCTION PHIMWUE(RR,A,TG,RHO,QAWUE,VDR) PARAMETER(KDIM=1000,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) DATA RHOC/6.0D0/,AWMG/24.312D0/ C Q_ext/a Wuestite, at 2 mu DATA QWUE/1.34D4/ IF(QAWUE.EQ.0D0) QAWUE=QWUE CALL VDUST(RHOC,QAWUE,RR,A,TG,RHO,VDR,PHI) PHIMWUE=PHI RETURN END DOUBLE PRECISION FUNCTION PHIMCOR(RR,A,TG,RHO,QACOR,VDR) PARAMETER(KDIM=1000,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) DATA RHOC/4.0D0/ C Q_ext/a Korundum, at 2 mu DATA QCOR/1.50D3/ IF(QACOR.EQ.0D0) QACOR=QCOR CALL VDUST(RHOC,QACOR,RR,A,TG,RHO,VDR,PHI) PHIMCOR=PHI RETURN END DOUBLE PRECISION FUNCTION PHIMHIB(RR,A,TG,RHO,QAHIB,VDR) PARAMETER(KDIM=1000,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) DATA RHOC/3.84D0/ C Q_ext/a Hibonite, at 2 mu !!!!!!!!!!!!!!! Offensichtlich Unsinn DATA QHIB/5.30D2/ IF(QAHIB.EQ.0D0) QAHIB=QHIB CALL VDUST(RHOC,QAHIB,RR,A,TG,RHO,VDR,PHI) PHIMHIB=PHI RETURN END DOUBLE PRECISION FUNCTION PHIMSPI(RR,A,TG,RHO,QASPI,VDR) PARAMETER(KDIM=1000,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) DATA RHOC/4.25D2/ C Q_ext/a Spinel, at 2 mu !!!!!!!!!!!!!!! Offensichtlich Unsinn DATA QSPI/3.43D0/ IF(QASPI.EQ.0D0) QASPI=QSPI CALL VDUST(RHOC,QASPI,RR,A,TG,RHO,VDR,PHI) PHIMSPI=PHI RETURN END C**********************************************************************C C Drift velocity of dust grains C C**********************************************************************C SUBROUTINE VDUST(RHOC,QDUST,R,A,TG,RHO,VDR,PHI) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT COMMON/STERN/STARL,STARM,STARR,TEFF C vfac=8*bk/(pi*m), sfac=2*bk/m, eddf=4*pi*cl*grav/3 DATA VFAC/9.023D7/,SFAC/7.087D7/,EDDF/3.353D4/ TAUBR=A*RHOC/(DSQRT(VFAC*TG)*RHO) VT=DSQRT(SFAC*TG) DLE=EDDF*STARM*RHOC/QDUST S0=TAUBR*GRAV*STARM/R**2*(STARL/DLE-1D0)/VT S2=(DSQRT(1D0+(9D0/16D0)*PI*S0**2)-1D0)*32D0/(9D0*PI) VDR=VT*DSQRT(S2) PHI=DSQRT(1D0+.25D0*PI*S2) RETURN END C***********************************************************************C C Search for grid-points for interpolation C C***********************************************************************C SUBROUTINE FIPOL(GRM,DM,IDIMG,IO,IL) DOUBLE PRECISION DM,GRM(IDIMG) IO=IDIMG IM=(IDIMG+1)/2 IL=1 100 CONTINUE IF(GRM(IM).LE.DM) THEN IL=IM ELSE IO=IM END IF IF(IO-IL.GE.2) THEN IM=(IO+IL)/2 GOTO 100 END IF RETURN END C**********************************************************************C C Interpolation of drift velocities C C**********************************************************************C SUBROUTINE VDUIP(R,VDOLI,VDPYR,VDQUA,VDIRO,VDWUE,VDCOR,VDHIB, * VDSPI) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(KDIM=1000,IFDIM=11,NSPDIM=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 CHARACTER*3 ADU LOGICAL FLDU COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) DIMENSION VI(NSPDIM) CALL FIPOL(RD,R,KK,IO,IL) IF(IL.EQ.KL.OR.IO.EQ.KK) THEN C At end-intervals RF=(R-RD(IL))/(RD(IO)-RD(IL)) DO L=1,8 VI(L)=VDR(IL,L)+(VDR(IO,L)-VDR(IL,L))*RF END DO ELSE C At middle-intervals IP=IL-1 IN=IO+1 DO L=1,8 D1L=(VDR(IL,L)-VDR(IP,L))/(RD(IL)-RD(IP)) D1M=(VDR(IO,L)-VDR(IL,L))/(RD(IO)-RD(IL)) D1R=(VDR(IN,L)-VDR(IO,L))/(RD(IN)-RD(IO)) D2M=(D1M-D1L)/(RD(IO)-RD(IP))+(D1R-D1M)/(RD(IN)-RD(IL)) VI(L)=VDR(IL,L)+X*(D1M+.5D0*X*D2M) END DO END IF VDOLI=VI(1) VDPYR=VI(2) VDQUA=VI(3) VDIRO=VI(4) VDWUE=VI(5) VDCOR=VI(6) VDHIB=VI(7) VDSPI=VI(8) RETURN END C***********************************************************************C C temperature interpolation C C***********************************************************************C DOUBLE PRECISION FUNCTION TEMP(R) PARAMETER(KDIM=1000,IFDIM=11,NSPDIM=11) IMPLICIT REAL*8(A-H,O-Z) 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/CTEMP/TS(KDIM,NSPDIM),TGH(KDIM) CALL FIPOL(RD,R,KK,IO,IL) TEMP=TG(IL)+(TG(IO)-TG(IL))*(R-RD(IL))/(RD(IO)-RD(IL)) RETURN END C***********************************************************************C C temperature interpolation C C***********************************************************************C SUBROUTINE TEMPO(R,TEMP,TOLI,TPYR,TQUA,TIRO,TWUE,TCOR,THIB,TSPI) PARAMETER(KDIM=1000,IFDIM=11,NSPDIM=11) IMPLICIT REAL*8(A-H,O-Z) 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/CTEMP/TS(KDIM,NSPDIM),TGH(KDIM) CHARACTER*3 ADU LOGICAL FLDU COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) DIMENSION TI(NSPDIM) CALL FIPOL(RD,R,KK,IO,IL) IF(IL.EQ.KL.OR.IO.EQ.KK) THEN C At end-intervals RF=(R-RD(IL))/(RD(IO)-RD(IL)) TEMP=TG(IL)+(TG(IO)-TG(IL))*RF DO L1=1,NDON L=NDONP(L1) TI(L1)=TS(IL,L)+(TS(IO,L)-TS(IL,L))*RF END DO ELSE C At middle-intervals IP=IL-1 IN=IO+1 D1L=(TG(IL)-TG(IP))/(RD(IL)-RD(IP)) D1M=(TG(IO)-TG(IL))/(RD(IO)-RD(IL)) D1R=(TG(IN)-TG(IO))/(RD(IN)-RD(IO)) D2M=(D1M-D1L)/(RD(IO)-RD(IP))+(D1R-D1M)/(RD(IN)-RD(IL)) X=R-RD(IL) TEMP=TG(IL)+X*(D1M+.5D0*X*D2M) DO L1=1,NDON L=NDONP(L1) D1L=(TS(IL,L)-TS(IP,L))/(RD(IL)-RD(IP)) D1M=(TS(IO,L)-TS(IL,L))/(RD(IO)-RD(IL)) D1R=(TS(IN,L)-TS(IO,L))/(RD(IN)-RD(IO)) D2M=(D1M-D1L)/(RD(IO)-RD(IP))+(D1R-D1M)/(RD(IN)-RD(IL)) TI(L1)=TS(IL,L)+X*(D1M+.5D0*X*D2M) END DO END IF TOLI=TI(1) TPYR=TI(2) TQUA=TI(3) TIRO=TI(4) TWUE=TI(5) TCOR=TI(6) THIB=TI(7) TSPI=TI(8) RETURN END C***********************************************************************C C extinction interpolation C C***********************************************************************C SUBROUTINE EXTDO(R,DKH,QAOLI,QAPYR,QAQUA,QAIRO,QAWUE,QACOR,QAHIB, * QASPI) PARAMETER(KDIM=1000,IFDIM=11,NSPDIM=11) IMPLICIT REAL*8(A-H,O-Z) 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/CRADM/VJ(KDIM),VH(KDIM),VK(KDIM),BV(KDIM,NSPDIM), * SPHF(KDIM),CHIBSP(KDIM,NSPDIM),CHIJSP(KDIM,NSPDIM), * CHIHSP(KDIM,NSPDIM) COMMON/CMOD/TAU(KDIM),TAL(KDIM),TAH(KDIM),CHIHT(KDIM),KLU CHARACTER*3 ADU LOGICAL FLDU COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) DIMENSION DKI(NSPDIM) CALL FIPOL(RD,R,KK,IO,IL) RF=(R-RD(IL))/(RD(IO)-RD(IL)) DKH=CHIHT(IL)+(CHIHT(IO)-CHIHT(IL))*RF DO L1=1,NDON L=NDONP(L1) DKI(L1)=CHIHSP(IL,L)+(CHIHSP(IO,L)-CHIHSP(IL,L))*RF END DO QAOLI=DKI(1) QAPYR=DKI(2) QAQUA=DKI(3) QAIRO=DKI(4) QAWUE=DKI(5) QACOR=DKI(6) QAHIB=DKI(7) QASPI=DKI(8) RETURN END C***********************************************************************C C flux averaged opacity C C***********************************************************************C DOUBLE PRECISION FUNCTION DKAPH(R) PARAMETER(KDIM=1000,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) 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 CALL FIPOL(RD,R,KK,IO,IL) DKL=CHIHT(IL) DKO=CHIHT(IO) DKAPH=DKL+(DKO-DKL)*(R-RD(IL))/(RD(IO)-RD(IL)) RETURN END C***********************************************************************C C temperature intialization C C***********************************************************************C SUBROUTINE INITEM(PARM) PARAMETER(KDIM=1000,IPDIM=21,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) 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 DIMENSION PARM(IPDIM) STARR=DSQRT(PARM(2)/(4D0*PI*SIG*PARM(3)**4)) RMAX=PARM(7)*STARR RMIN=2.0D0*STARR H=10D0**(DLOG10(RMAX/RMIN)/(KK-1D0)) RD(1)=RMIN DO K=2,KK RD(K)=RD(K-1)*H END DO TAU=0D0 V=1D6 DO K=1,KK TG(K)=PARM(3)* * DSQRT(DSQRT(.5D0*(1D0-DSQRT(1D0-(STARR/RD(K))**2)+1.5D0*TAU))) DKAP(K)=0D0 RHO(K)=PARM(5)/(4D0*PI*RD(K)**2*V) DO L=1,IFDIM VDR(K,L)=1D0 END DO END DO RETURN END C***********************************************************************C C storing data C C for M-star element mixture C C***********************************************************************C SUBROUTINE STOREO(YN,YA,PARM,KNEXT) PARAMETER(KDIM=1000,IPDIM=21,IYDIM=IPDIM-4,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) 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/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/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/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM COMMON/CVDR/VDROLI,VDRPYR,VDRQUA,VDRIRO,VDRWUE,VDRCOR,VDRHIB, * VDRSPI,VDRCAR,VDRSIC,VDRMGS DIMENSION YN(IYDIM),YA(IYDIM),YI(IYDIM),PARM(IPDIM) IF(KNEXT.LT.0) GOTO 200 C -- model -- IF(YN(1).LT.RD(KNEXT)) RETURN 100 F=(RD(KNEXT)-YA(1))/(YN(1)-YA(1)) DO I=2,NDUSPO+6 YI(I)=YN(I)+(YN(I)-YA(I))*F END DO VV(KNEXT)=YI(2) RHO(KNEXT)=PARM(5)/(4D0*PI*VV(KNEXT)*RD(KNEXT)**2) TG(KNEXT)=YI(NDUSPO+5) GAM(KNEXT)=YI(NDUSPO+6) C -- degrees of condensation -- C Y(5) : aOli C Y(6) : xOli C Y(7) : aPyr C Y(8) : xPyr C Y(9) : aQua C Y(10): aFe C Y(11): aWu C Y(12): xWu C Y(13): aCo C Y(14): aHi C Y(15): aSp FV(KNEXT,1)=4D0*PI*(YI(5)**3-A0OLI**3)*EPSOLI/(3D0*V0OLI*EPSSI) FV(KNEXT,2)=YI(6) FV(KNEXT,3)=4D0*PI*(YI(7)**3-A0PYR**3)*EPSPYR/(3D0*V0PYR*EPSSI) FV(KNEXT,4)=YI(8) FV(KNEXT,5)=4D0*PI*(YI(9)**3-A0QUA**3)*EPSQUA/(3D0*V0QUA*EPSSI) FV(KNEXT,6)=4D0*PI*(YI(10)**3-A0IRO**3)*EPSIRO/(3D0*V0IRO*EPSFE) FV(KNEXT,7)= (YI(11)**3-A0WUE**3)*EPSWUE/(V0WUE*EPSMG) FV(KNEXT,8)=YI(12) FV(KNEXT,9)=4D0*PI*(YI(13)**3-A0COR**3)*EPSCOR/(3D0*V0COR*EPSAL) FV(KNEXT,10)=4D0*PI*(YI(14)**3-A0HIB**3)*EPSHIB/(3D0*V0HIB*EPSAL) FV(KNEXT,11)=4D0*PI*(YI(15)**3-A0SPI**3)*EPSSPI/(3D0*V0SPI*EPSAL) AV(KNEXT,1)=YI(5) AV(KNEXT,2)=YI(7) AV(KNEXT,3)=YI(9) AV(KNEXT,4)=YI(13) AV(KNEXT,5)=YI(14) AV(KNEXT,6)=YI(15) AV(KNEXT,7)=YI(10) AV(KNEXT,8)=YI(11) AV(KNEXT,9)=0D0 AV(KNEXT,10)=0D0 AV(KNEXT,11)=0D0 VDR(KNEXT,1)=VDROLI VDR(KNEXT,2)=VDRPYR VDR(KNEXT,3)=VDRQUA VDR(KNEXT,4)=VDRIRO VDR(KNEXT,5)=VDRWUE VDR(KNEXT,6)=VDRCOR VDR(KNEXT,7)=VDRHIB VDR(KNEXT,8)=VDRSPI C -- extinction coefficients -- DKAP(KNEXT)=PARM(20) C -- increase pointer -- IF(KNEXT.LT.KK) THEN KNEXT=KNEXT+1 IF(YN(1).GE.RD(KNEXT)) GOTO 100 END IF RETURN C -- initialization -- 200 STARR=DSQRT(PARM(2)/(4D0*PI*SIG*PARM(3)**4)) DO K=1,-KNEXT VV(K)=PARM(4) RHO(K)=PARM(5)/(4D0*PI*VV(K)*RD(K)**2) TG(K)=PARM(3)*DSQRT(DSQRT(.5D0*(1D0-DSQRT(1D0-(STARR/RD(K))**2) * +1.5D0*PARM(6)))) DKAP(K)=0D0 DO I=1,IFDIM FV(K,I)=0D0 END DO END DO RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C S u b r o u t i n e s f o r C - s t a r s C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C***********************************************************************C C Integration of the wind equations (driver) C C for carbon rich element mixture C C***********************************************************************C SUBROUTINE MODELLC(PARM,NEWMO) PARAMETER(KDIM=1000,IPDIM=21,IYDIM=IPDIM-4,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) 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/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/CCAR/V0CAR,A0CAR,EPSCAR,ALFCAR COMMON/CSIC/V0SIC,A0SIC,EPSSIC,ALFSIC COMMON/CIRO/V0IRO,A0IRO,EPSIRO,ALFIRO COMMON/CMGS/V0MGS,A0MGS,EPSMGS,ALFMGS COMMON/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM DIMENSION YN(IYDIM),YA(IYDIM),DYDTN(IYDIM-1),DYDTA(IYDIM-1), * PARM(IPDIM) LOGICAL NEWMO 1 FORMAT(I8,1P20E10.3) C -- outer boundary -- STARR=DSQRT(PARM(2)/(4D0*PI*SIG*PARM(3)**4)) RMAX=PARM(7)*STARR C -- initial values -- RMIN=2.5D0*STARR CALL FIPOL(RD,RMIN,KK,KO,KL) YA(1)=RMIN YA(2)=PARM(4) YA(3)=PARM(6) YA(4)=0D0 YA(5)=A0CAR YA(6)=A0SIC YA(7)=A0IRO YA(8)=A0MGS IF(NEWMO) THEN CALL STOREC(YN,YA,PARM,-KL) END IF CALL DGLSC(DT,YA,DYDTA,PARM) C -- loop for time integration -- DT=1D5 DTA=DT KNEXT=KO ISTEP=0 1000 ISTEP=ISTEP+1 1010 CALL NEXTTC(DT,DTA,DTN,YN,YA,DYDTN,DYDTA,PARM) IF(MOD(ISTEP,10000).EQ.0) WRITE(6,1) ISTEP,DT,(YN(L),L=1,10) IF(DTN.LT.0.5D0*DT) THEN DT=DTN c PRINT *,'repeat last step:',ISTEP GOTO 1010 END IF CALL STOREC(YN,YA,PARM,KNEXT) C -- new time step-width -- DTA=DT IF(DTN.GT.1.2D0*DT) DT=1.2D0*DT IF(DTN.LT.DT) DT=DTN IF(DT.LT.1D0) DT=1D0 C -- save last step -- DO N=1,NDUSPC+6 YA(N)=YN(N) END DO DO N=1,NDUSPC+4 DYDTA(N)=DYDTN(N) END DO IF(YN(1).LT.RMAX) GOTO 1000 C -- store final values -- PARM(NDUAM+1)=4D0*PI*(YN(5)**3-A0CAR**3)*EPSCAR/(3D0*V0CAR*EPSC) PARM(NDUAM+2)=4D0*PI*(YN(6)**3-A0SIC**3)*EPSSIC/(3D0*V0SIC*EPSSI) PARM(NDUAM+3)=4D0*PI*(YN(7)**3-A0IRO**3)*EPSIRO/(3D0*V0IRO*EPSFE) PARM(NDUAM+4)=4D0*PI*(YN(8)**3-A0MGS**3)*EPSMGS/(3D0*V0MGS*EPSS) PARM(NDUSPC+8)=YN(2)*1D-5 PARM(NDUSPC+9)=YN(4) PARM(NDUSPC+10)=YN(3) RETURN END C***********************************************************************C C Integration of wind equations (Stepper) C C Method: Adams - Bashforth C C for carbon rich element mixture C C***********************************************************************C SUBROUTINE NEXTTC(DT,DTA,DTN,YN,YA,DYDTN,DYDTA,PARM) PARAMETER(IPDIM=21,IYDIM=IPDIM-4) IMPLICIT REAL*8(A-H,O-Z) DIMENSION YN(IYDIM),YA(IYDIM),DYDTN(IYDIM-1),DYDTA(IYDIM-1), * PARM(IPDIM) COMMON/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM COMMON/CTE/CHM,RCOS,RCOC CALL DGLSC(DT,YA,DYDTN,PARM) YN(NDUSPC+5)=YA(NDUSPC+5) YN(NDUSPC+6)=YA(NDUSPC+6) DO N=1,NDUSPC+4 YN(N)=YA(N)+DT*(DYDTN(N)+.5D0*(DYDTN(N)-DYDTA(N))*DT/DTA) END DO DTN=1D30 DO N=1,NDUSPC+4 X=CHM*DABS(YN(N)/DYDTN(N)) IF((N.EQ.3).OR.(N.EQ.4)) X=1D30 IF(X.LT.DTN) DTN=X END DO RETURN END C***********************************************************************C C Right hand sides of the DEQ-System C C for carbon rich element mixture C C***********************************************************************C C ---- Independent variables C Y(1) : r C Y(2) : v C Y(3) : tauL C Y(4) : tau* C Y(5) : aCar C Y(6) : aSiC C Y(7) : aFe C Y(8) : aMgS C Y(9) : Tg C Y(10) : Gamma C ---- Parameters C PARM(1) =STARM C PARM(2) =STARL C PARM(3) =TEFF C PARM(4) =VEXP C PARM(5) =DMP C PARM(6) =tauLini C PARM(7) =Rout C PARM(8) =FCA C PARM(9) =FSC C PARM(10)=FIR C PARM(11)=FMS C PARM(12)=VWIND C PARM(13)=kapH C PARM(14)=tauLfin C-----------------------------------------------------------------------C SUBROUTINE DGLSC(DT,Y,DYDT,PARM) PARAMETER(IPDIM=21,IYDIM=IPDIM-4) IMPLICIT REAL*8(A-H,O-Z) 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/STERN/STARL,STARM,STARR,TEFF COMMON/CAW/AWSIO,AWMG,AWFE,AWH2O,AWC2H2,AWSI,AWSIS,AWH2S,AWAL, * AWFEO,AWMGO,AWCA COMMON/CCAR/V0CAR,A0CAR,EPSCAR,ALFCAR COMMON/CSIC/V0SIC,A0SIC,EPSSIC,ALFSIC COMMON/CIRO/V0IRO,A0IRO,EPSIRO,ALFIRO COMMON/CMGS/V0MGS,A0MGS,EPSMGS,ALFMGS COMMON/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM COMMON/CVDR/VDROLI,VDRPYR,VDRQUA,VDRIRO,VDRWUE,VDRCOR,VDRHIB, * VDRSPI,VDRCAR,VDRSIC,VDRMGS DIMENSION Y(IYDIM),DYDT(IYDIM-1),PARM(IPDIM) C -- wind model -- R=Y(1) V=Y(2) STARM=PARM(1) STARL=PARM(2) TEFF=PARM(3) STARR=DSQRT(PARM(2)/(4D0*PI*SIG*PARM(3)**4)) IF(PARM(6).LE.0D0) THEN CALL TEMPC(R,TG,TCAR,TSIC,TIRO,TMGS) ELSE TAU=DMAX1(Y(3),0D0) TG=PARM(3)* * DSQRT(DSQRT(.5D0*(1D0-DSQRT(1D0-(STARR/R)**2)+1.5D0*TAU))) TCAR=TG TSIC=TG TIRO=TG TMGS=TG END IF Y(NDUSPC+5)=TG RHO=PARM(5)/(4D0*PI*V*R**2) CGAS2=BK*TG/(2.333D0*PM) C -- degrees of condensation -- FCA=4D0*PI*(Y(5)**3-A0CAR**3)*EPSCAR/(3D0*V0CAR*EPSC) FSC=4D0*PI*(Y(6)**3-A0SIC**3)*EPSSIC/(3D0*V0SIC*EPSSI) FIR=4D0*PI*(Y(7)**3-A0IRO**3)*EPSIRO/(3D0*V0IRO*EPSFE) FMS=4D0*PI*(Y(8)**3-A0MGS**3)*EPSMGS/(3D0*V0MGS*EPSS) C -- radius -- DYDT(1)=V C -- kappaH -- IF(PARM(6).LE.0D0) THEN CALL EXTDC(R,DKH,QACAR,QASIC,QAIRO,QAMGS) ELSE TRAD=(TEFF-TG)*DEXP(-Y(4))+TG DKH=DKCAR(TRAD)*FCA*EPSC/3.24D-4+DKSIC(TRAD)*FSC+DKIR(TRAD)*FIR * +DKGAS(TG,RHO) QACAR=0D0 QASIC=0D0 QAIRO=0D0 QAMGS=0D0 END IF PARM(13)=DKH C -- velocity -- GAMMA=PARM(2)/(4D0*PI*CL*GRAV*PARM(1))*DKH DYDT(2)=DMAX1(-GRAV*PARM(1)*(1D0-GAMMA)/R**2,1D-30) Y(NDUSPC+6)=GAMMA C -- optical depths -- DYDT(3)=-V*RHO*DKH*(STARR/R)**2 DYDT(4)=V*RHO*DKH C -- molecules -- PHG=RHO*BK*TG/((1D0+4D0*EPSHE)*PM) PH2=.5D0*PHG PC2H2=.5D0*(EPSC*(1D0-FCA)-EPSO-FSC*EPSSI)*PHG IF((1D0-FSC)*EPSSI.GT.EPSS) THEN PSIS=(1D0-FMS)*EPSS*PHG PSI=((1D0-FSC)*EPSSI-EPSS)*PHG PH2S=0D0 ELSE PH2S=(EPSS-(1D0-FSC)*EPSSI)*PHG PSIS=(1D0-FMS)*EPSS*PHG-PH2S PSI=0D0 END IF PFE=(1D0-FIR)*EPSFE*PHG PMG=(EPSMG-FMS*EPSS)*PHG C --- carbon --- G=ALFCAR/DSQRT(2D0*PI*AWC2H2*PM*BK*TG) IF(TG.LT.1100) THEN DYDT(5)=V0CAR*G*2D0*PC2H2*PHICCAR(R,Y(5),TG,RHO,QACAR,VDRCAR) ELSE DYDT(5)=0D0 ENDIF C -- Iron -- G=ALFIRO/DSQRT(2D0*PI*AWFE*PM*BK*TG) IF(TG.GT.500D0) THEN PFEEQ=PVFE(TIRO) DYDT(7)=V0IRO*G*DMAX1(PFE*PHIMIRN(R,Y(7),TG,RHO,QAIRO,VDRIRO) * -PFEEQ,0D0) ELSE DYDT(7)=V0IRO*G*PFE*PHIMIRN(R,Y(7),TG,RHO,QAIRO,VDRIRO) ENDIF C -- SiC -- G=ALFSIC/DSQRT(2D0*PI*AWSI*PM*BK*TG) G2=G*DSQRT(AWSI/AWC2H2) DJGRC2H2=G2*PC2H2 IF(TG.GT.650D0) THEN PSIEQ=PAT*DSQRT(PH2/(PC2H2*DKSICSIG(TSIC))) ELSE PSIEQ=0D0 END IF DJGRSI=G*DMAX1(PSI-PSIEQ,0D0) DJGRCA=DMAX1(.5D0*G2*PC2H2-G*PSIEQ,0D0) DJGR=DMIN1(DJGRSI,DJGRCA)*PHICSIC(R,Y(6),TG,RHO,QASIC,VDRSIC) DYDT(6)=V0SIC*DJGR C -- MgS -- G=ALFMGS/DSQRT(2D0*PI*AWH2S*PM*BK*TG) IF(TG.GT.500D0) THEN PH2SEQ=PAT*PH2/(PMG*DKMGSMG(TSIC)) ELSE PH2SEQ=0D0 END IF DYDT(8)= * V0MGS*G*DMAX1((PH2S+PSIS)*PHICMGS(R,Y(8),TG,RHO,QAMGS,VDRMGS) * -PH2SEQ,0D0) RETURN END C***********************************************************************C C Correction factors for dust drift C C***********************************************************************C DOUBLE PRECISION FUNCTION PHICCAR(RR,A,TG,RHO,QACAR,VDR) PARAMETER(KDIM=1000,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) DATA RHOC/2.26D0/ C Q_ext/a Carbon, Draine, at 2 mu DATA QCAR/2.65D4/ IF(QACAR.EQ.0D0) QACAR=QCAR CALL VDUST(RHOC,QACAR,RR,A,TG,RHO,VDR,PHI) PHICCAR=PHI RETURN END DOUBLE PRECISION FUNCTION PHICSIC(RR,A,TG,RHO,QASIC,VDR) PARAMETER(KDIM=1000,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) DATA RHOC/3.217D0/ C Q_ext/a Silicon carbide, Pegourie, at 2 mu DATA QSIC/1.86D3/ IF(QASIC.EQ.0D0) QASIC=QSIC CALL VDUST(RHOC,QASIC,RR,A,TG,RHO,VDR,PHI) PHICSIC=PHI RETURN END C Iron: as in O-rich case DOUBLE PRECISION FUNCTION PHICMGS(RR,A,TG,RHO,QAMGS,VDR) PARAMETER(KDIM=1000,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) DATA RHOC/2.84D0/ C Q_ext/a MgS, Biegeman, at 2 mu DATA QMGS/8.65D2/ IF(QAMGS.EQ.0D0) QAMGS=QMGS CALL VDUST(RHOC,QAMGS,RR,A,TG,RHO,VDR,PHI) PHICMGS=PHI RETURN END C***********************************************************************C C temperature interpolation C C***********************************************************************C SUBROUTINE TEMPC(R,TEMP,TCAR,TSIC,TIRO,TMGS) PARAMETER(KDIM=1000,IFDIM=11,NSPDIM=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/CTEMP/TS(KDIM,NSPDIM),TGH(KDIM) COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) DIMENSION TI(NSPDIM) CHARACTER*3 ADU CALL FIPOL(RD,R,KK,IO,IL) RF=(R-RD(IL))/(RD(IO)-RD(IL)) TEMP=TG(IL)+(TG(IO)-TG(IL))*RF DO L1=1,NDON L=NDONP(L1) TI(L1)=TS(IL,L)+(TS(IO,L)-TS(IL,L))*RF END DO TCAR=TI(1) TSIC=TI(2) TIRO=TI(3) TMGS=TI(4) RETURN END C***********************************************************************C C extinction interpolation C C***********************************************************************C SUBROUTINE EXTDC(R,DKH,QACAR,QASIC,QAIRO,QAMGS) PARAMETER(KDIM=1000,IFDIM=11,NSPDIM=11) IMPLICIT REAL*8(A-H,O-Z) 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/CRADM/VJ(KDIM),VH(KDIM),VK(KDIM),BV(KDIM,NSPDIM), * SPHF(KDIM),CHIBSP(KDIM,NSPDIM),CHIJSP(KDIM,NSPDIM), * CHIHSP(KDIM,NSPDIM) COMMON/CMOD/TAU(KDIM),TAL(KDIM),TAH(KDIM),CHIHT(KDIM),KLU CHARACTER*3 ADU LOGICAL FLDU COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) DIMENSION DKI(NSPDIM) CALL FIPOL(RD,R,KK,IO,IL) RF=(R-RD(IL))/(RD(IO)-RD(IL)) DKH=CHIHT(IL)+(CHIHT(IO)-CHIHT(IL))*RF DO L1=1,NDON L=NDONP(L1) DKI(L1)=CHIHSP(IL,L)+(CHIHSP(IO,L)-CHIHSP(IL,L))*RF END DO QACAR=DKI(1) QASIC=DKI(2) QAIRO=DKI(3) QAMGS=DKI(4) RETURN END C***********************************************************************C C storing data C C for C-star element mixture C C***********************************************************************C SUBROUTINE STOREC(YN,YA,PARM,KNEXT) PARAMETER(KDIM=1000,IPDIM=21,IYDIM=IPDIM-4,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) 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/CCAR/V0CAR,A0CAR,EPSCAR,ALFCAR COMMON/CSIC/V0SIC,A0SIC,EPSSIC,ALFSIC COMMON/CIRO/V0IRO,A0IRO,EPSIRO,ALFIRO COMMON/CMGS/V0MGS,A0MGS,EPSMGS,ALFMGS 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/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM COMMON/CVDR/VDROLI,VDRPYR,VDRQUA,VDRIRO,VDRWUE,VDRCOR,VDRHIB, * VDRSPI,VDRCAR,VDRSIC,VDRMGS DIMENSION YN(IYDIM),YA(IYDIM),YI(IYDIM),PARM(IPDIM) IF(KNEXT.LT.0) GOTO 200 C -- model -- 100 IF(YN(1).LT.RD(KNEXT)) RETURN F=(RD(KNEXT)-YA(1))/(YN(1)-YA(1)) DO I=2,NDUSPC+6 YI(I)=YN(I)+(YN(I)-YA(I))*F END DO VV(KNEXT)=YI(2) RHO(KNEXT)=PARM(5)/(4D0*PI*VV(KNEXT)*RD(KNEXT)**2) TG(KNEXT)=YI(NDUSPC+5) GAM(KNEXT)=YI(NDUSPC+6) C -- degrees of condensation -- FV(KNEXT,1)=4D0*PI*(YI(5)**3-A0CAR**3)*EPSCAR/(3D0*V0CAR*EPSC) FV(KNEXT,2)=4D0*PI*(YI(6)**3-A0SIC**3)*EPSSIC/(3D0*V0SIC*EPSSI) FV(KNEXT,3)=4D0*PI*(YI(7)**3-A0IRO**3)*EPSIRO/(3D0*V0IRO*EPSFE) FV(KNEXT,4)=4D0*PI*(YI(8)**3-A0MGS**3)*EPSMGS/(3D0*V0MGS*EPSS) AV(KNEXT,1)=0D0 AV(KNEXT,2)=0D0 AV(KNEXT,3)=0D0 AV(KNEXT,4)=0D0 AV(KNEXT,5)=0D0 AV(KNEXT,6)=0D0 AV(KNEXT,7)=YI(7) AV(KNEXT,8)=0D0 AV(KNEXT,9)=YI(5) AV(KNEXT,10)=YI(6) C mantle thickness MgS RMGS=(YI(8)**3-A0MGS**3+YI(6)**3)**(1D0/3D0) AV(KNEXT,11)=RMGS-YI(6) VDR(KNEXT,1)=VDRCAR VDR(KNEXT,2)=VDRSIC VDR(KNEXT,3)=VDRIRO VDR(KNEXT,4)=VDRMGS C -- extinction coefficients -- TGAS=YI(NDUSPC+5) TRAD=(PARM(3)-TGAS)*DEXP(-YI(4))+TGAS DKAP(KNEXT)=PARM(13) IF(KNEXT.LT.KK) THEN KNEXT=KNEXT+1 IF(YN(1).GE.RD(KNEXT)) GOTO 100 END IF RETURN C -- initialization -- 200 STARR=DSQRT(PARM(2)/(4D0*PI*SIG*PARM(3)**4)) DO K=1,-KNEXT VV(K)=PARM(4) RHO(K)=PARM(5)/(4D0*PI*VV(K)*RD(K)**2) TG(K)=PARM(3)*DSQRT(DSQRT(.5D0*(1D0-DSQRT(1D0-(STARR/RD(K))**2) * +1.5D0*PARM(6)))) DKAP(K)=0D0 DO I=1,IFDIM FV(K,I)=0D0 END DO END DO RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C S u b r o u t i n e s f o r S - s t a r s C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C***********************************************************************C C Iteration of temperature stratification (Lucy-approximation) C C for S-star element mixture C C***********************************************************************C SUBROUTINE ITTEMPLS(PARM) PARAMETER(IPDIM=21) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT COMMON/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM DIMENSION PARM(IPDIM) 10 FORMAT(1X,'it:',I3,3X,'tauLi:',1PE11.4,3X,'tauLf:',E11.4,3X, * 'L:',E11.4,3X,'Teff:',E11.4,3X,'M:',E11.4,3X,'Mp:',E11.4) C -- initial values -- TAULH=PARM(6) CALL INITEM(PARM) TU=TAULH TO=TAULH TM=TAULH C -- bisection method to determine TAUL -- IM=1 IT=0 100 IT=IT+1 IF(IT.GT.100) GOTO 200 PARM(6)=TM CALL MODELLS(PARM) TAUL=PARM(NDUSPS+10) IF(IT.GT.90) WRITE(6,10) IT,TM,TAUL,PARM(2)/SUNL, * PARM(3),PARM(1)/SUNM,PARM(5)*YEAR/SUNM IF(IM.NE.1) GOTO 110 C -- right boundary too small -- IF(TAUL.LT.0D0) THEN TU=TM TO=1.2D0*TO TM=TO ELSE IM=2 TM=TU END IF GOTO 100 110 IF(IM.NE.2) GOTO 120 C -- left boundary too big -- IF(TAUL.GE.0D0) THEN TU=.5D0*TU TM=TU ELSE IM=3 TM=.5D0*(TU+TO) END IF GOTO 100 120 IF((TO-TU)/TU.GT.1D-4.OR.TAUL.GT.1D-8) THEN C -- bisection -- IF(TAUL.GT.0D0) THEN TO=TM ELSE TU=TM END IF TM=.5D0*(TU+TO) GOTO 100 END IF C -- method successful -- PARM(NDUSPS+10)=TM PARM(6)=TAULH RETURN 200 IF(DABS(TAUL).LT.1D-10) RETURN PRINT *,'Iteration of TAUL in ITTEMPLS did not converge' STOP END C***********************************************************************C C Integration of the wind equations (driver) C C for S-star element mixture C C***********************************************************************C SUBROUTINE MODELLS(PARM) PARAMETER(KDIM=1000,IPDIM=21,IYDIM=IPDIM-4,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) 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/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/CIRO/V0IRO,A0IRO,EPSIRO,ALFIRO COMMON/CQUA/V0QUA,A0QUA,EPSQUA,ALFQUA COMMON/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM DIMENSION YN(IYDIM),YA(IYDIM),DYDTN(IYDIM-1),DYDTA(IYDIM-1), * PARM(IPDIM) 1 FORMAT(I8,1P20E10.3) C -- outer boundary -- STARR=DSQRT(PARM(2)/(4D0*PI*SIG*PARM(3)**4)) RMAX=PARM(7)*STARR C -- initial values -- RMIN=2.5D0*STARR CALL FIPOL(RD,RMIN,KK,KO,KL) YA(1)=RMIN YA(2)=PARM(4) YA(3)=PARM(6) YA(4)=0D0 YA(5)=A0IRO YA(6)=A0QUA CALL STORES(YN,YA,PARM,-KL) CALL DGLSS(DT,YA,DYDTA,PARM) C -- loop for time integration -- DT=1D5 DTA=DT KNEXT=KO ISTEP=0 1000 ISTEP=ISTEP+1 1010 CALL NEXTTS(DT,DTA,DTN,YN,YA,DYDTN,DYDTA,PARM) IF(DTN.LT.0.5D0*DT) THEN DT=DTN c PRINT *,'Repeat last step:',ISTEP GOTO 1010 END IF CALL STORES(YN,YA,PARM,KNEXT) C -- new time step-width -- DTA=DT IF(DTN.GT.1.2D0*DT) DT=1.2D0*DT IF(DTN.LT.DT) DT=DTN c -- save last step -- DO N=1,NDUSPS+6 YA(N)=YN(N) END DO DO N=1,NDUSPS+4 DYDTA(N)=DYDTN(N) END DO IF(YN(1).LT.RMAX) GOTO 1000 C -- store final values -- PARM(NDUAM+1)=4D0*PI*(YN(5)**3-A0IRO**3)*EPSIRO/(3D0*V0IRO*EPSFE) PARM(NDUAM+2)=4D0*PI*(YN(6)**3-A0QUA**3)*EPSQUA/(3D0*V0QUA*EPSSI) PARM(NDUSPS+8)=YN(2)*1D-5 PARM(NDUSPS+9)=YN(4) PARM(NDUSPS+10)=YN(3) RETURN END C***********************************************************************C C Integration of wind equations (Stepper) C C method: Adams - Bashforth C C for S-star element mixture C C***********************************************************************C SUBROUTINE NEXTTS(DT,DTA,DTN,YN,YA,DYDTN,DYDTA,PARM) PARAMETER(IPDIM=21,IYDIM=IPDIM-4) IMPLICIT REAL*8(A-H,O-Z) DIMENSION YN(IYDIM),YA(IYDIM),DYDTN(IYDIM-1),DYDTA(IYDIM-1), * PARM(IPDIM) COMMON/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM COMMON/CTE/CHM,RCOS,RCOC CALL DGLSS(DT,YA,DYDTN,PARM) YN(NDUSPS+5)=YA(NDUSPS+5) YN(NDUSPS+6)=YA(NDUSPS+6) DO N=1,NDUSPS+4 YN(N)=YA(N)+DT*(DYDTN(N)+.5D0*(DYDTN(N)-DYDTA(N))*DT/DTA) END DO DTN=1D30 DO N=1,NDUSPS+4 X=CHM*DABS(YN(N)/DYDTN(N)) IF((N.EQ.3).OR.(N.EQ.4)) X=1D30 IF(X.LT.DTN) DTN=X END DO RETURN END C**********************************************************************C C Right hand sides of the DEQ-System C C for S-star element mixture C C**********************************************************************C C ---- Independent variables C Y(1) : r C Y(2) : v C Y(3) : tauL C Y(4) : tau* C Y(5) : aIro C Y(6) : aQua C Y(7) : Tg C Y(8) : Gamma C ---- Parameters C PARM(1) =STARM C PARM(2) =STARL C PARM(3) =TEFF C PARM(4) =VEXP C PARM(5) =DMP C PARM(6) =tauLini C PARM(7) =Rout C PARM(8) =FIR C PARM(9) =FQU C PARM(10)= not used C PARM(11)= not used C PARM(12)=VWIND C PARM(13)=kapH C PARM(14)=tauLfin C-----------------------------------------------------------------------C SUBROUTINE DGLSS(DT,Y,DYDT,PARM) PARAMETER(IPDIM=21,IYDIM=IPDIM-4) IMPLICIT REAL*8(A-H,O-Z) 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/STERN/STARL,STARM,STARR,TEFF COMMON/CAW/AWSIO,AWMG,AWFE,AWH2O,AWC2H2,AWSI,AWSIS,AWH2S,AWAL, * AWFEO,AWMGO,AWCA COMMON/CIRO/V0IRO,A0IRO,EPSIRO,ALFIRO COMMON/CQUA/V0QUA,A0QUA,EPSQUA,ALFQUA COMMON/CVDR/VDROLI,VDRPYR,VDRQUA,VDRIRO,VDRWUE,VDRCOR,VDRHIB, * VDRSPI,VDRCAR,VDRSIC,VDRMGS COMMON/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM DIMENSION Y(IYDIM),DYDT(IYDIM-1),PARM(IPDIM) C -- wind model -- R=Y(1) V=Y(2) STARM=PARM(1) STARL=PARM(2) TEFF=PARM(3) STARR=DSQRT(STARL/(4D0*PI*SIG*TEFF**4)) IF(PARM(6).LE.0D0) THEN CALL TEMPS(R,TG,TIRO,TQUA) ELSE TAU=DMAX1(Y(3),0D0) TG=TEFF* * DSQRT(DSQRT(.5D0*(1D0-DSQRT(1D0-(STARR/R)**2)+1.5D0*TAU))) TIRO=TG TQUA=TV END IF Y(NDUSPS+5)=TG RHO=PARM(5)/(4D0*PI*V*R**2) CGAS2=BK*TG/(2.333D0*PM) C -- degrees of condensation -- FIR=4D0*PI*(Y(5)**3-A0IRO**3)*EPSIRO/(3D0*V0IRO*EPSFE) FQU=4D0*PI*(Y(6)**3-A0QUA**3)*EPSQUA/(3D0*V0QUA*EPSSI) C -- radius -- DYDT(1)=V C -- kappaH -- IF(PARM(6).LE.0D0) THEN CALL EXTDS(R,DKH,QAIRO,QAQUA) ELSE TRAD=(TEFF-TG)*DEXP(-Y(4))+TG DKH=DKIR(TRAD)*FIR+DKSIL(TRAD)*FQU+DKGAS(TG,RHO) QAIRO=0D0 QAQUA=0D0 END IF PARM(13)=DKH C -- velocity -- GAMMA=STARL/(4D0*PI*CL*GRAV*STARM)*DKH DYDT(2)=DMAX1(-GRAV*PARM(1)*(1D0-GAMMA)/R**2,1D-30) Y(NDUSPS+6)=GAMMA C -- optical depths -- DYDT(3)=-V*RHO*DKH*(STARR/R)**2 DKP=DKIR(TRAD)*FIR+DKSIL(TRAD)*FQU+DKGAS(TG,RHO) DYDT(4)=V*RHO*DKP C -- molecules -- PHG=RHO*BK*TG/((1D0+4D0*EPSHE)*PM) PH2O=DMAX1((EPSO-EPSC-(1D0+FQU)*EPSSI)*PHG,0D0) PSIO=((1D0-FQU)*EPSSI)*PHG PFE=(1D0-FIR)*EPSFE*PHG C -- iron -- GF=ALFIRO/DSQRT(2D0*PI*AWFE*PM*BK*TG) PGR=PFE*PHIMIRN(R,Y(5),TG,RHO,QAIRO,VDRIRO) IF(TIRO.GT.500D0) THEN PFEEQ=PVFE(TIRO) ELSE PFEEQ=0D0 ENDIF DYDT(5)=V0IRO*GF*DMAX1(PGR-PFEEQ,0D0) C -- quartz -- GF=ALFQUA/DSQRT(2D0*PI*AWSIO*PM*BK*TG) PGRH2O=PH2O*DSQRT(AWSIO/AWH2O) PGRQU=DMIN1(PSIO,PGRH2O)*PHIMQUA(R,Y(6),TG,RHO,QUQUA,VDRQUA) IF(TG.GT.500D0.AND.PH2O.GT.0D0) THEN PQUSIOEQ=PAT*PH2/(PH2O*DKPQUSIO(TG)) ELSE PQUSIOEQ=0D0 END IF DYDT(6)=V0QUA*GF*DMAX1(PGRQU-PQUSIOEQ,0D0) RETURN END C***********************************************************************C C temperature interpolation C C***********************************************************************C SUBROUTINE TEMPS(R,TEMP,TIRO,TQUA) PARAMETER(KDIM=1000,IFDIM=11,NSPDIM=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/CTEMP/TS(KDIM,NSPDIM),TGH(KDIM) COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) DIMENSION TI(NSPDIM) CHARACTER*3 ADU CALL FIPOL(RD,R,KK,IO,IL) RF=(R-RD(IL))/(RD(IO)-RD(IL)) TEMP=TG(IL)+(TG(IO)-TG(IL))*RF DO L1=1,NDON L=NDONP(L1) TI(L1)=TS(IL,L)+(TS(IO,L)-TS(IL,L))*RF END DO TIRO=TI(1) TQUA=TI(2) RETURN END C***********************************************************************C C extinction interpolation C C***********************************************************************C SUBROUTINE EXTDS(R,DKH,QAIRO,QAQUA) PARAMETER(KDIM=1000,IFDIM=11,NSPDIM=11) IMPLICIT REAL*8(A-H,O-Z) 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/CRADM/VJ(KDIM),VH(KDIM),VK(KDIM),BV(KDIM,NSPDIM), * SPHF(KDIM),CHIBSP(KDIM,NSPDIM),CHIJSP(KDIM,NSPDIM), * CHIHSP(KDIM,NSPDIM) COMMON/CMOD/TAU(KDIM),TAL(KDIM),TAH(KDIM),CHIHT(KDIM),KLU CHARACTER*3 ADU LOGICAL FLDU COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) DIMENSION DKI(NSPDIM) CALL FIPOL(RD,R,KK,IO,IL) RF=(R-RD(IL))/(RD(IO)-RD(IL)) DKH=CHIHT(IL)+(CHIHT(IO)-CHIHT(IL))*RF DO L1=1,NDON L=NDONP(L1) DKI(L1)=CHIHSP(IL,L)+(CHIHSP(IO,L)-CHIHSP(IL,L))*RF END DO QAIRO=DKI(1) QAQUA=DKI(2) RETURN END C**********************************************************************C C storing data C C for S-star element mixture C C**********************************************************************C SUBROUTINE STORES(YN,YA,PARM,KNEXT) PARAMETER(KDIM=1000,IPDIM=21,IYDIM=IPDIM-4,IFDIM=11) IMPLICIT REAL*8(A-H,O-Z) 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/CIRO/V0IRO,A0IRO,EPSIRO,ALFIRO COMMON/CQUA/V0QUA,A0QUA,EPSQUA,ALFQUA 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/CVDR/VDROLI,VDRPYR,VDRQUA,VDRIRO,VDRWUE,VDRCOR,VDRHIB, * VDRSPI,VDRCAR,VDRSIC,VDRMGS COMMON/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM DIMENSION YN(IYDIM),YA(IYDIM),YI(IYDIM),PARM(IPDIM) IF(KNEXT.LT.0) GOTO 200 C -- wind model -- IF(YN(1).LT.RD(KNEXT)) RETURN 100 F=(RD(KNEXT)-YA(1))/(YN(1)-YA(1)) DO I=2,NDUSPS+6 YI(I)=YN(I)+(YN(I)-YA(I))*F END DO VV(KNEXT)=YI(2) RHO(KNEXT)=PARM(5)/(4D0*PI*VV(KNEXT)*RD(KNEXT)**2) TG(KNEXT)=YI(NDUSPS+5) GAM(KNEXT)=YI(NDUSPS+6) C -- degrees of condensation -- FV(KNEXT,1)=4D0*PI*(YI(5)**3-A0IRO**3)*EPSIRO/(3D0*V0IRO*EPSFE) FV(KNEXT,2)=4D0*PI*(YI(6)**3-A0QUA**3)*EPSQUA/(3D0*V0QUA*EPSSI) VDR(KNEXT,1)=VDRIRO VDR(KNEXT,2)=VDRQUA C -- extinction coefficient -- TGAS=YI(NDUSPS+5) TRAD=(PARM(3)-TGAS)*DEXP(-YI(4))+TGAS DKAP(KNEXT)=DKIR(TRAD)*FV(KNEXT,1)+DKSIL(TRAD)*FV(KNEXT,2)+ * DKGAS(TG(KNEXT),RHO(KNEXT)) IF(KNEXT.LT.KK) THEN KNEXT=KNEXT+1 IF(YN(1).GE.RD(KNEXT)) GOTO 100 END IF RETURN C -- initialization -- 200 STARR=DSQRT(PARM(2)/(4D0*PI*SIG*PARM(3)**4)) DO K=1,-KNEXT VV(K)=PARM(4) RHO(K)=PARM(5)/(4D0*PI*VV(K)*RD(K)**2) TG(K)=PARM(3)*DSQRT(DSQRT(.5D0*(1D0-DSQRT(1D0-(STARR/RD(K))**2) * +1.5D0*PARM(6)))) DKAP(K)=0D0 DO I=1,IFDIM FV(K,I)=0D0 END DO END DO RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C S u b r o u t i n e s f o r t h e r m o d y n a m i c d a t a C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C************************** Olivine ***********************************C C mass-action constant for ol+3H_2-->SiO+2xMg+2(1-x)Fe+3H_2O C C (to use with pressures in bar) C C**********************************************************************C SUBROUTINE KPOLVAP(TG,XOL,DKPOL) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT GMG2SIO4=7.52334D+4/TG-9.38369D+5+(2.47581D+2-3.14980D-3*TG)*TG GFE2SIO4=6.84477D+4/TG-9.02146D+5+(2.51045D+2-4.44028D-3*TG)*TG GMIX=2*RGAS*TG*(XOL*DLOG(XOL)+(1D0-XOL)*DLOG(1D0-XOL)) GOL=XOL*GMG2SIO4+(1D0-XOL)*GFE2SIO4+GMIX GH2=4.25321D+5/TG-1.07123D+5+(2.69980D+1+(5.48280D-4- * 3.81498D-8*TG)*TG)*TG GH2O=8.66184D+5/TG-2.27857D+5+(5.61473D+1+(7.62548D-4- * 4.95254D-8*TG)*TG)*TG GSIO=1.82780D+5/TG-1.92839D+5+(3.03804D+1+(4.16079D-4- * 1.91035D-8*TG)*TG)*TG DKPOL=DEXP(-(GOL+3D0*(GH2-GH2O)-GSIO)/(RGAS*TG)) RETURN END C*************************** Olivine **********************************C C mass-action constant for calculating the composition x of olivine C C (to use with pressures in bar) C C**********************************************************************C SUBROUTINE OLIKP(TG,XOL,DKXOL) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT DGFOR=-1430140D0/TG-383328D0+(129.753D0+(-4.20868D-3+ * 5.32303D-8*TG)*TG)*TG DGFAY=-1436920D0/TG-347105D0+(133.217D0+(-5.49916D-3+ * 5.32303D-8*TG)*TG)*TG DGMIX=2*RGAS*TG*(XOL*DLOG(XOL)+(1D0-XOL)*DLOG(1D0-XOL)) DKXOL=DEXP(-(DGFAY-DGFOR+DGMIX/(1D0-XOL))/(2D0*RGAS*TG)) RETURN END C************************** Pyroxene **********************************C C mass-action constant for py+2H_2-->SiO+xMg+(1-x)Fe+2H_2O C C (to use with pressures in bar) C C**********************************************************************C SUBROUTINE KPPYVAP(TG,XPY,DKPPY) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT GMGSIO3=8.74400D+3/TG-6.92565D+5+(1.77877D+2-1.41412D-3*TG)*TG GFESIO3=-5.26506D+5/TG-6.68373D+05+(1.73783D+2+(-2.13971D-3+ * 4.94433D-7*TG)*TG)*TG GMIX=RGAS*TG*(XPY*DLOG(XPY)+(1D0-XPY)*DLOG(1D0-XPY)) GPY=XPY*GMGSIO3+(1D0-XPY)*GFESIO3+GMIX GH2=4.25321D+5/TG-1.07123D+5+(2.69980D+1+(5.48280D-4- * 3.81498D-8*TG)*TG)*TG GH2O=8.66184D+5/TG-2.27857D+5+(5.61473D+1+(7.62548D-4- * 4.95254D-8*TG)*TG)*TG GSIO=1.82780D+5/TG-1.92839D+5+(3.03804D+1+(4.16079D-4- * 1.91035D-8*TG)*TG)*TG DKPPY=DEXP(-(GPY+2D0*(GH2-GH2O)-GSIO)/(RGAS*TG)) RETURN END C************************** Pyroxene **********************************C C mass-action constant for calculating the composition x of pyroxene C C (to use with pressures in bar) C C**********************************************************************C SUBROUTINE PYRKP(TG,XPY,DKXPY) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT GMGSIO3=8.74400D+3/TG-6.92565D+5+(1.77877D+2-1.41412D-3*TG)*TG GFESIO3=-5.26506D+5/TG-6.68373D+05+(1.73783D+2+(-2.13971D-3+ * 4.94433D-7*TG)*TG)*TG DGMIX=RGAS*TG*(XPY*DLOG(XPY)+(1D0-XPY)*DLOG(1D0-XPY)) DKXPY=DEXP(-(GFESIO3-GMGSIO3+DGMIX/(1D0-XPY))/(RGAS*TG)) RETURN END C***************************** Quartz *********************************C C mass-action constant for SiO_2(s) + H_2 --> SiO + H_2O C C (to use with pressures in bar) C C**********************************************************************C DOUBLE PRECISION FUNCTION DKPQUSIO(TG) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT DGSIO=1.82780D+5/TG-1.92839D+5+(3.03804D+1+(4.16079D-4- * 1.91035D-8*TG)*TG)*TG DGH2=4.25321D+5/TG-1.07123D+5+(2.69980D+1+(5.48280D-4- * 3.81498D-8*TG)*TG)*TG DGH2O=8.66184D+5/TG-2.27857D+5+(5.61473D+1+(7.62548D-4- * 4.95254D-8*TG)*TG)*TG DGQUA=-4.44364D+5+(1.08531D+2-6.59213D-4*TG)*TG DKPQUSIO=DEXP(-(DGQUA+DGH2-DGSIO-DGH2O)/(RGAS*TG)) RETURN END C****************************** Iron **********************************C C vapor pressure iron (dyn/cm^2) C C**********************************************************************C DOUBLE PRECISION FUNCTION PVFE(TG) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT DGFE=-97638.6D0+(35.5605D0-.00112348D0*TG)*TG PVFE=PAT*DEXP(DGFE/(RGAS*TG)) RETURN END C**************************** Corundum ********************************C C mass-action constant for Al_2O_3(s) + 3H_2 --> 2Al + 3H_2O C C (to use with pressures in bar) C C**********************************************************************C DOUBLE PRECISION FUNCTION DKPCOAL(TG) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT DGH2=4.25321D+5/TG-1.07123D+5+(2.69980D+1+(5.48280D-4- * 3.81498D-8*TG)*TG)*TG DGH2O=8.66184D+5/TG-2.27857D+5+(5.61473D+1+(7.62548D-4- * 4.95254D-8*TG)*TG)*TG DGAL2O3=-7.32976D+5+(1.84782D+2-2.57313D-3*TG)*TG DKPCOAL=DEXP(-(DGAL2O3+3D0*(DGH2-DGH2O))/(RGAS*TG)) RETURN END C**************************** Spinel **********************************C C mass-action constant for MgAl_2O_4(s) + 4H_2 --> Mg + 2Al + 4H_2O C C (to use with pressures in bar) C C**********************************************************************C DOUBLE PRECISION FUNCTION DKPSPAL(TG) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT DGH2=4.25321D+5/TG-1.07123D+5+(2.69980D+1+(5.48280D-4- * 3.81498D-8*TG)*TG)*TG DGH2O=8.66184D+5/TG-2.27857D+5+(5.61473D+1+(7.62548D-4- * 4.95254D-8*TG)*TG)*TG DGMGAL2O4=2.50397D+5/TG-9.83584D+5+(2.54352D+2-4.24568D-3*TG)*TG DKPSPAL=DEXP(-(DGMGAL2O4+4D0*(DGH2-DGH2O))/(RGAS*TG)) RETURN END C*************************** Wuestite *********************************C C Equilibrium partial pressure of FeO over W"ustite C C**********************************************************************C DOUBLE PRECISION FUNCTION DKPWUFEO(TG,XWU) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT DGFEOs=-6.05591D+4/TG-2.22156D+5+6.78769D+1*TG-1.52278D-3*TG*TG GH2=4.25321D+5/TG-1.07123D+5+(2.69980D+1+(5.48280D-4- * 3.81498D-8*TG)*TG)*TG GH2O=8.66184D+5/TG-2.27857D+5+(5.61473D+1+(7.62548D-4- * 4.95254D-8*TG)*TG)*TG GAMFEO=DEXP((525D0/TG+0.48D0)*(1D0-XWU)**2) DKPWUFEO=XWU*GAMFEO*DEXP((DGFEOs-GH2O+GH2)/(RGAS*TG)) RETURN END C**************************** Wuestite ********************************C C Equilibrium partial pressure of MgO over W"ustite C C**********************************************************************C DOUBLE PRECISION FUNCTION DKPWUMGO(TG,XWU) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT DGMGOs=4.38728D+3/TG-2.38741D+5+(6.86582D+1+(-1.19852D-3- * 5.72304D-8*TG)*TG)*TG GH2=4.25321D+5/TG-1.07123D+5+(2.69980D+1+(5.48280D-4- * 3.81498D-8*TG)*TG)*TG GH2O=8.66184D+5/TG-2.27857D+5+(5.61473D+1+(7.62548D-4- * 4.95254D-8*TG)*TG)*TG GAMMGO=DEXP((386D0/TG+0.539D0)*XWU**2) DKPWUMGO=(1D0-XWU)*GAMMGO*DEXP((DGMGOs-GH2O+GH2)/(RGAS*TG)) RETURN END C***************************** Siliconcarbide *************************C C mass-action constant for 2Si + C2H2 -> 2SiC(s) + H2 C C (to use with pressures in bar) C C**********************************************************************C DOUBLE PRECISION FUNCTION DKSICSIG(TG) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT GSIC=6.68748D+4/TG-2.96740D+5+(7.645711D+1+(-8.17922D-4+ * 1.30155D-7*TG)*TG)*TG GH2=4.25321D+5/TG-1.07123D+5+(2.69980D+1+(5.48280D-4- * 3.81498D-8*TG)*TG)*TG GC2H2=8.11613D+5/TG-3.99041D+5+(9.15250D+1+(2.53889D-5- * 2.03949D-8*TG)*TG)*TG DKSICSIG=DEXP(-(2D0*GSIC+GH2-GC2H2)/(RGAS*TG)) RETURN END C************************** Siliconcarbide ****************************C C mass-action constant for 2SiS + C2H2 + H2 -> 2SiC(s) + 2H2S C C (to use with pressures in bar) C C**********************************************************************C DOUBLE PRECISION FUNCTION DKSICSIS(TG) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT GSIC=6.68748D+4/TG-2.96740D+5+(7.645711D+1+(-8.17922D-4+ * 1.30155D-7*TG)*TG)*TG GH2=4.25321D+5/TG-1.07123D+5+(2.69980D+1+(5.48280D-4- * 3.81498D-8*TG)*TG)*TG GC2H2=8.11613D+5/TG-3.99041D+5+(9.15250D+1+(2.53889D-5- * 2.03949D-8*TG)*TG)*TG GH2S=7.41911D+5/TG-1.81063D+5+(5.35114D+1+(4.94594D-4- * 3.39231D-8*TG)*TG)*TG GSIS=1.44507D+5/TG-1.49916D+5+(2.87381D+1+(3.96186D-4- * 2.18002D-8*TG)*TG)*TG DKSICSIS=DEXP(-(2D0*(GSIC+GH2S-GSIS)-GC2H2-GH2)/(RGAS*TG)) RETURN END C*************************** Magnesiumsulfide *************************C C mass-action constant for Mg + H2S -> MgS(s) + H2 C C (to use with pressures in bar) C C**********************************************************************C DOUBLE PRECISION FUNCTION DKMGSMG(TG) IMPLICIT REAL*8(A-H,O-Z) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT sGMGS=-7.93442D+4/TG-1.83584D+5+(6.33920D+1-8.54820D-4*TG)*TG GH2=4.25321D+5/TG-1.07123D+5+(2.69980D+1+(5.48280D-4- * 3.81498D-8*TG)*TG)*TG GH2S=7.41911D+5/TG-1.81063D+5+(5.35114D+1+(4.94594D-4- * 3.39231D-8*TG)*TG)*TG DKMGSMG=DEXP(-(sGMGS+GH2-GH2S)/(RGAS*TG)) RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C S u b r o u t i n e s f o r e x t i n c t i o n C C c o e f f i c i e n t s C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C**********************************************************************C C Extinction by dirty silicate (Ossenkopf) C C**********************************************************************C DOUBLE PRECISION FUNCTION DKSIL(TG) IMPLICIT REAL*8(A-H,O-Z) DATA a1r/1.721D-2/,e1r/1.102D0/,a2r/2.155e1/,e2r/0.238D0/, * a3r/9.119D-3/,e3r/0.975D0/,a4r/2.563D-1/,e4r/0.516D0/ f1r=a1r*TG**e1r f2r=a2r/TG**e2r f3r=a3r*TG**e3r f4r=a4r*TG**e4r fhr=-dlog(dexp(-f3r)+dexp(-f4r)) flr=1./(1./f1r**8+1./f2r**8)**.125 DKSIL=(fhr**4+flr**4)**.25 RETURN END C********************************************************************** C Extinction by iron * C********************************************************************** DOUBLE PRECISION FUNCTION DKIR(TG) IMPLICIT NONE REAL*8 TG DKIR=1.D0/DSQRT(DSQRT(1.D0/(3.341D-5*TG**1.632)**4+ * 1.D0/((6.405D-4*TG**.777)**4+(4.385D-7*TG**1.981)**4))) END FUNCTION DKIR C*********************************************************************** C Extinction by amorphous olivine * C*********************************************************************** DOUBLE PRECISION FUNCTION DKAMOL(TG) IMPLICIT NONE REAL*8 TG DKAMOL=1D0/Dsqrt((1D0/(3.240D-5*TG**2.061D0))**2+ * (1D0/(1.067D4/TG**1.713D0+2.022D-8*TG**2.383D0))**2) END FUNCTION DKAMOL C**********************************************************************C C Extinction by amorphous pyroxene C C**********************************************************************C DOUBLE PRECISION FUNCTION DKAMPY(TG) IMPLICIT REAL*8(A-H,O-Z) DKAMPY=1d0/dsqrt(1d0/(3.773d-5*TG**2.017)**2+1d0/( * (1.5d6/TG**2.679)**2+dsqrt((8.356d-4*TG**0.7336)**4+ * (1.261d-8*TG**2.272)**4))) RETURN END C**********************************************************************C C Extinction by quartz C C**********************************************************************C DOUBLE PRECISION FUNCTION DKQUA(TG) IMPLICIT REAL*8(A-H,O-Z) DKQUA=dsqrt(1D0/dsqrt(1d0/(3.898d-6*TG**2.045)**4+1d0/ * (3.100d5/TG**2.622)**4)+sqrt((2.023d-5*TG**1.074)**4+ * (9.394d-11*TG**2.701)**4)) RETURN END C**********************************************************************C C Extinction by amorphous corundum (Planck-mean) C C**********************************************************************C DOUBLE PRECISION FUNCTION PKCOR(TG) IMPLICIT REAL*8(A-H,O-Z) PKCOR=(1D0/(1D0/(2.201d-8*TG**3.021)+ * 1D0/(5.037d2*TG**(-1.366)))**4+(5.725d-8*TG**1.770)**4)**.25D0 RETURN END C**********************************************************************C C Extinction by amorphous corundum (Rosseland-mean) C C**********************************************************************C DOUBLE PRECISION FUNCTION RKCOR(TG) IMPLICIT REAL*8(A-H,O-Z) RKCOR=1D0/(1D0/(3.288d-8*TG**2.547)**2+ * 1D0/(1.192d3*TG**(-1.888))**2)**.5+ * DSQRT((8.300d-5*TG**0.735)**2+(3.808d-10*TG**2.376)**2) RETURN END C**********************************************************************C C Extinction by carbon C C**********************************************************************C DOUBLE PRECISION FUNCTION DKCAR(TG) IMPLICIT REAL*8(A-H,O-Z) DKCAR=1.d0/dsqrt(dsqrt(1.d0/(1.692d-4*TG**1.813)**4+1.d0/ * ((3.167d-2*TG**.612)**4+(2.938d-5*TG**1.750)**4))) RETURN END C**********************************************************************C C Extinction by silicon carbide C C**********************************************************************C DOUBLE PRECISION FUNCTION DKSIC(TG) IMPLICIT REAL*8(A-H,O-Z) DKSIC=1d0/dsqrt(dsqrt(1d0/(9.560d-6*TG**1.710d0)**4+1d0/ * ((7.791d-4*TG**.8808d0)**4+(1.171d-6*TG**1.879d0)**4))) RETURN END C**********************************************************************C C Extinction by the gas C C**********************************************************************C DOUBLE PRECISION FUNCTION DKGAS(TG,RHO) IMPLICIT REAL*8(A-H,O-Z) DKGAS=1D-8*RHO**0.66667D0*TG**3 RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C B a s i c D a t a C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C**********************************************************************C C Initializing basic data C C**********************************************************************C SUBROUTINE INIWIND(ZSCALE,CORATIO) IMPLICIT REAL*8(A-H,O-Z) PARAMETER(KDIM=1000,IFDIM=11,IPDIM=21) 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/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/CTE/CHM,RCOS,RCOC COMMON/CNDUSP/NDUSPO,NDUSPC,NDUSPS,NDUAM DIMENSION ZEL(30),AEL(30) C******************************* Constants ***************************** HP=6.62618D-27 ! Planck constant CL=2.99793D+10 ! speed of light BK=1.38066D-16 ! Boltzmann constant PM=1.67265D-24 ! proton mass SIG=5.6703D-5 ! radiation constant GRAV=6.672D-8 ! gravitational constant PI=3.141592653589793D0 SUNM=1.991D+33 ! solar mass SUNR=6.9595D+10 ! solar radius SUNL=3.84D+33 ! solar luminosity YEAR=3.1558D+7 ! seconds/year RGAS=1.987165D0 ! gas constant in cal/mole PAT=1013250D0 ! conversion pressure in atm to dyn/cm^2 C**************************** Element abundances ********************** C -- Grevesse, N., Sauval, A.J. 1998, Space Sci. Rev., 85, 161 -- C -- Allende Prieto et al. ApJL 556, L63 ; 573, L137 EPSHE=8.51D-2 EPSC =2.45D-4 EPSN =8.32D-5 EPSO =4.90D-4 EPSNE=1.20D-4 EPSNA=2.14D-6 EPSMG=3.80D-5 EPSAL=2.95D-6 EPSSI=3.55D-5 EPSS =2.14D-5 EPSAR=3.63D-6 EPSCA=2.29D-6 EPSFE=3.16D-5 EPSNI=1.78D-6 OPEN(10,FILE='abun.txt',STATUS='OLD',ERR=20) N=1 10 READ(10,*,END=11) ZEL(N),AEL(N) N=N+1 GOTO 10 11 CLOSE(10) DO I=1,N IF(ZEL(I).EQ.2) EPSHE=AEL(I) IF(ZEL(I).EQ.6) EPSC =AEL(I) IF(ZEL(I).EQ.7) EPSN =AEL(I) IF(ZEL(I).EQ.8) EPSO =AEL(I) IF(ZEL(I).EQ.10) EPSNE=AEL(I) IF(ZEL(I).EQ.11) EPSNA=AEL(I) IF(ZEL(I).EQ.12) EPSMG=AEL(I) IF(ZEL(I).EQ.13) EPSAL=AEL(I) IF(ZEL(I).EQ.14) EPSSI=AEL(I) IF(ZEL(I).EQ.16) EPSS =AEL(I) IF(ZEL(I).EQ.18) EPSAR=AEL(I) IF(ZEL(I).EQ.20) EPSCA=AEL(I) IF(ZEL(I).EQ.26) EPSFE=AEL(I) IF(ZEL(I).EQ.28) EPSNI=AEL(I) END DO 20 CONTINUE c Scaling EPSC =ZSCALE*EPSC EPSN =ZSCALE*EPSN EPSO =ZSCALE*EPSO EPSNE=ZSCALE*EPSNE EPSNA=ZSCALE*EPSNA EPSMG=ZSCALE*EPSMG EPSAL=ZSCALE*EPSAL EPSSI=ZSCALE*EPSSi EPSS =ZSCALE*EPSS EPSAR=ZSCALE*EPSAR EPSCA=ZSCALE*EPSCA EPSFE=ZSCALE*EPSFE EPSNI=ZSCALE*EPSNI IF(CORATIO.GT.0D0) EPSC=CORATIO*EPSO C*************************** technical parameters ********************** C maximum change of variables per time step CHM=0.01D0 C C/O ratio at the M-S-transition RCOS=0.9D0 C C/O ratio at the S-C-transition RCOC=1D0 C*********************** Parameters for dust growth ******************** CALL DEFST C*************** Number of parameters for M, S, C stars **************** C Oxygen rich mixture NDUSPO=11 C Carbon rich mixture NDUSPC=4 C Mixture for S-stars NDUSPS=2 C Check I=10+MAX(NDUSPO,NDUSPC,NDUSPS) IF(IPDIM.LT.I) THEN PRINT*,'IPDIM too small, at least',I,' is required' STOP END IF NDUAM=7 END SUBROUTINE INIWIND C**********************************************************************C C Definition of parameters for dust formation C C**********************************************************************C SUBROUTINE DEFST IMPLICIT REAL*8(A-H,O-Z) PARAMETER(NSPDIM=11,IPDIM=21) COMMON/CONST/HP,CL,BK,PM,SIG,GRAV,PI,SUNM,SUNR,SUNL,YEAR,RGAS,PAT COMMON/CAW/AWSIO,AWMG,AWFE,AWH2O,AWC2H2,AWSI,AWSIS,AWH2S,AWAL, * AWFEO,AWMGO,AWCA 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/CASOL/AFOR,AFAY,AENS,AFRS,AQUA,ASIC,ASIL CHARACTER*3 ADU LOGICAL FLDU COMMON/CCON/NDONP(NSPDIM),NDON,IFIRSTC,FLDU(NSPDIM),ADU(NSPDIM) C --------- MgFeSiO4 (olivine) -------- ADU(1)='Oli' AFOR=140.69D0 ! molecular weight forsterite AFAY=203.77D0 ! molecular weight fayalite AOLI=172.23D0 ! molecular weight of olivine with x=0.5 RHOOLI=3.75D0 ! mass density of olivine V0OLI=AOLI*PM/RHOOLI ! volume of one formula unit EPSOLI=3D-13 ! number of seed particles per H nucleus A0OLI=1D-7 ! radius of seed particles AWSIO=44.09D0 ! molecular weight SiO AWMG=24.312D0 ! atomic weight Mg AWFE=55.847D0 ! atomic weight Fe AWH2O=18.051D0 ! molecular weight H_2O ALFOLI=.15D0 ! growth coefficient for olivine ALFOLEX=1D0*ALFOLI ! exchange coefficient olivine C --------- (Mg,Fe)SiO3 (pyroxene) ---- ADU(2)='Pyr' AENS=100.39D0 ! molecular weight enstatite AFRS=131.93D0 ! molecular weight ferrosilite APYR=116.12D0 ! molecular weight of pyroxene with x=0.5 RHOPYR=3.58D0 ! mass density of pyroxene V0PYR=APYR*PM/RHOPYR ! volume of one formula unit EPSPYR=3D-13 ! number of seed particles per H nucleus A0PYR=1D-7 ! radius of seed particles ALFPYR=.15D0 ! growth coefficient for pyroxene ALFPYEX=1.D0*ALFPYR ! exchange coefficient pyroxene C -------- Quartz -------- ADU(3)='Qtz' AQUA=60.08D0 ! molecular weight quartz RHOQUA=2.2D0 ! mass density of quartz V0QUA=AQUA*PM/RHOQUA ! volume of one formula unit EPSQUA=1D-13 ! number of seed particles per H nucleus A0QUA=1D-7 ! radius of seed particles ALFQUA=.05D0 ! growth coefficient for quartz C -------- Corundum -------- ADU(4)='Cor' ACOR=101.96D0 ! molecular weight corundum RHOCOR=3.97D0 ! mass density of corundum V0COR=ACOR*PM/RHOCOR ! volume of one formula unit EPSCOR=3D-13 ! number of seed particles per H nucleus A0COR=1D-7 ! radius of seed particles AWAL=26.982D0 ! atomic weight of Al ALFCOR=.10D0 ! growth coefficient for corundum C -------- Hibonite -------- ADU(5)='Hib' AHIB=667.8D0 ! molecular weight hibonite RHOHIB=3.84D0 ! mass density of hibonite V0HIB=AHIB*PM/RHOHIB ! volume of one formula unit EPSHIB=3D-13 ! number of seed particles per H nucleus A0HIB=1D-7 ! radius of seed particles AWCA=40.078D0 ! atomic weight of Ca ALFHIB=0.10D0 ! growth coefficient for hibonite C --------- Spinel --------- ADU(6)='Spi' ASPI=142.27D0 ! molecular weight spinel RHOSPI=3.58D0 ! mass density of spinel V0SPI=ASPI*PM/RHOSPI ! volume of one formula unit EPSSPI=3D-13 ! number of seed particles per H nucleus A0SPI=1D-7 ! radius of seed particles ALFSPI=0.10D0 ! growth coefficient for spinel C --------- Iron ----------- ADU(7)='Irn' AIRO=55.847D0 ! molecular weight iron RHOIRO=7.86D0 ! mass density of iron V0IRO=AIRO*PM/RHOIRO ! volume of one formula unit EPSIRO=1D-13 ! number of seed particles per H nucleus A0IRO=1D-7 ! radius of seed particles ALFIRO=0.8D0 ! growth coefficient for iron C --------- Wuestite ------- ADU(8)='Wue' AL=2.65D-8 ! Length of cube edge V0WUE=AL**3 ! volume of one formula unit RHOWUE=6.0D0 EPSWUE=3D-13 ! number of seed particles per H nucleus A0WUE=1D-7 ! size of seed particles ALFWUMG=0.00D0 ! growth coefficient for Mg ALFWUFE=ALFWUMG*DSQRT(AWFE/AWMG) ! growth coefficient for Fe AWFEO=71.844D0 ! molecular weight FeO AWMGO=40.304D0 ! molecular weight MgO C --------- Carbon --------- ADU(9)='Car' ACAR=12D0 ! molecular weight carbon RHOCAR=2.26D0 ! mass density of carbon V0CAR=ACAR*PM/RHOCAR ! volume of one formula unit EPSCAR=3D-13 ! number of seed particles per H nucleus A0CAR=1D-7 ! radius of seed particles AWC2H2=26D0 ! molecular weight C_2H_2 ALFCAR=.3D0 ! growth coefficient for carbon C --------- SiC ------------ ADU(10)='SiC' ASIC=40.1D0 ! molecular weight silicon carbide RHOSIC=3.217D0 ! mass density of silicon carbide V0SIC=ASIC*PM/RHOSIC ! volume of one formula unit EPSSIC=3D-13 ! number of seed particles per H nucleus A0SIC=1D-7 ! radius of seed particles AWSI=28D0 ! atomic weight Si AWSIS=60D0 ! molecular weight SiS ALFSIC=0.3D0 ! growth coefficient for silicon carbide C --------- MgS ------------ ADU(11)='MgS' AMGS=56.38D0 ! molecular weight magnesium sulphide RHOMGS=2.84D0 ! mass density of magnesium sulphide V0MGS=AMGS*PM/RHOMGS ! volume of one formula unit EPSMGS=1D-13 ! number of seed particles per H nucleus A0MGS=5D-6 ! radius of seed particles AWH2S=34D0 ! molecular weight H_2S ALFMGS=0.3D0 ! growth coefficient for magnesium sulphide RETURN END