**==STATEF.FOR SUBROUTINE STATEF(FL, TL) IMPLICIT REAL*8(A-H, N-Z) PARAMETER (EPS=1.0D-30) COMMON /AUXIN / ICL, ION, JW, IOP, INUC, IBC, ICN, IML(2), ISGTH, : IMO, IDIFF COMMON /STAT1 / CSX(10), CS(90,127,10), CNIU(60,41,2), W(18400), : JCSX COMMON /STAT2 / PL, RL, U, P, RHO, FK, T, SF, ST, ZT, GRADA, CP, : CH, S, PR, PG, PF, PT, EN, WR(41) COMMON /STATFD/ RE, PE, QE, RET, PET, QET, REF, PEF, QEF, FDD(7) COMMON /ABUND / XA(10), NA(10), NEO, NIO, NZZ, AVM COMMON /ATDATA/ DH2, D1, D2, D3, CHI(26,9), OMG(27), AM(10), BN(10), & IZ(10) COMMON /CNSTS / CPI, PI4, CLN10, CA, CB, CC, CD, CG, CR, CT, CEVB, & CEN, CPL, CW(6) COMMON /IONISE/ NSTORE(8,5), NSTORE2(5) COMMON /OPDAT / cbase,obase,opT(141),opR(31),ZS FXP(VX) = DEXP(MAX(-50.0D0,MIN(50.0D0,VX))) DIMENSION HA(26), VA(26), tkappa(2,2),Xcompos(5),COcompos(8) data Xcompos /0.0d0,0.03d0,0.1d0,0.35d0,0.7d0/ data COcompos /0.0d0,0.01d0,0.03d0,0.1d0,0.2d0,0.4d0,0.6d0,1.0d0/ DATA JT, JX /2, 2/ * Evaluate Fermi-Dirac integrals according to Eggleton, Faulkner & * Flannery (1973): RE, PE, SE and UE correspond to rho*, P*, S* and U*. * PSI is the usual degeneracy parameter F = EXP(FL) T = EXP(TL) UF = F/(1.0+F) WF = SQRT(1.0+F) PSI = FL + 2.0*(WF-LOG(1.0+WF)) G = CT*T*WF CALL FDIRAC(F, G) PE = G*PE PET = PET + 1.0 PEF = PEF + 0.5*UF QE = QE/(RE*WF) SE = QE + 2.0*WF - PSI SEF = QE*(QEF-REF-0.5*UF) - 1.0/WF SET = QE*(QET-RET) UE = SE + PSI - PE/(RE*CT*T) * Some quantities that do not depend on the state of ionization: * the NA are the element number densities (per baryon); AVM is the average * mass per baryon (in amu); NEO and NIO are the numbers of electrons and * ions assuming complete ionization. NEO = 0.0 NIO = 0.0 NZZ = 0.0 AVM = 0.0 DO I = 1, 10 NA(I) = XA(I)/BN(I) AVM = AVM + AM(I)*NA(I) NIO = NIO + NA(I) NEO = NEO + IZ(I)*NA(I) NZZ = NZZ + IZ(I)*IZ(I)*NA(I) END DO * PRESSI gives a crude model for pressure ionization and (if ICL=1) * a model for Coulomb interactions, returning corrections to the electron * chemical potential, pressure, entropy and internal energy. * TI is 1eV/kT, DE is 1 amu * number of electrons/cm3 TI = CEVB/T DE = RE*CD CALL PRESSI(ICL, TI, DE, REF, RET, F, DC, DVT, DVF, DPA, DPAT, & DPAF, DSA, DSAT, DSAF, DUA) DV = DC - PSI DVF = DVF - WF * Contributions of the completely ionized species NE = 0.0 NEF = 0.0 NET = 0.0 SI = 0.0 SIF = 0.0 SIT = 0.0 UI = 0.0 DO I = ION+1, 9 NE = NE + IZ(I)*NA(I) VM = AM(I)*SQRT(AM(I)) SI = SI - NA(I)*LOG(NA(I)/VM + EPS) END DO * Calculate ionization of the first ION elements. DO I = ION, 1, -1 SHA = 1.0 SJHA = 0.0 SCHA = 0.0 * compute potentials VA and number ratios HA of ionization state J * relative to the ground state DO J = 1, IZ(I) VA(J) = -CHI(J,I)*TI + J*DV IF(J.EQ.1) THEN HA(J) = FXP(VA(J))*OMG(IZ(I))/OMG(IZ(I)+1) ELSE HA(J) = HA(J-1)*FXP(VA(J) - VA(J-1))*OMG(IZ(I)+1-J) & /OMG(IZ(I)+2-J) END IF SHA = SHA + HA(J) SJHA = SJHA + J*HA(J) SCHA = SCHA + CHI(J,I)*TI*HA(J) NSTORE(J,I) = HA(J) END DO NSTORE2(I) = SHA VM = AM(I)*SQRT(AM(I)) SI = SI + NA(I)*LOG(VM) IF(I.GT.1) THEN * contributions to electron number density NE, entropy SI and * internal energy UI for Helium and heavier VX = NA(I)/SHA SI = SI - VX*LOG(VX/OMG(IZ(I)+1) + EPS) DO J = 1, IZ(I) NX = HA(J)*VX NXF = NX*DVF*(J - SJHA/SHA) NXT = NXF*DVT/DVF + NX*(CHI(J,I)*TI - SCHA/SHA) NE = NE + J*NX NEF = NEF + J*NXF NET = NET + J*NXT SIF = SIF - VA(J)*NXF SIT = SIT - VA(J)*NXT SI = SI - NX*LOG(NX/OMG(IZ(I)+1-J) + EPS) UI = UI + CHI(J,I)*NX END DO END IF END DO * Ionization and molecular dissciation of Hydrogen. * partition function for H2 from Vardya (1960), Webbink (1975) DH2TI = DH2*TI D1TI = D1*TI D2TI = (D2*TI)**2 D3TI = (D3*TI)**3 ZET = 1.0 - (1.0 + DH2TI)*EXP(-DH2TI) DZH2T = -DH2TI**2*EXP(-DH2TI)/ZET DZH2TT = (DH2TI - 2.0 - DZH2T)*DZH2T ZH2 = 6608.8*ZET*DH2TI**(-2.5)*EXP(-D1TI - D2TI - D3TI) ZH2T = 2.5 + D1TI + 2.0*D2TI + 3.0*D3TI + DZH2T ZH2TT = -D1TI - 4.0*D2TI - 9.0*D3TI + DZH2TT ZH2S = ZH2*SQRT(8D0)/VM H2A = CEN*(ZH2S/4.0)*DE/(T*SQRT(T))*EXP(DH2TI) H2BT = DH2TI + 1.5 - ZH2T H2AT = RET - H2BT * solve for densities of H+, H, and H2 AA = 2*H2A + HA(1)*(1.0 + HA(1)) AB = NE + HA(1)*(NE - NA(1)) AC = NA(1)*NE HG = 2*AC/(SQRT(AB*AB + 4*AA*AC) + AB) HI = HA(1)*HG NE = NE + HI EN = 1.0/NE H2 = H2A*HG*HG*EN NI = NIO - H2 * derivatives w.r.t. F and T AA = NE + 4*HG*H2A AB = HA(1)*(NE - 2*H2) AD = 1/(AA + AB) AC = 2*H2*AD AF = (NEF - NE*REF)*AC AT = (NET - NE*H2AT)*AC AE = HG*AD BF = DVF*AE BT = (CHI(1,1)*TI + DVT)*AE HGF = AF - AB*BF HGT = AT - AB*BT HIF = HA(1)*(AF + AA*BF) HIT = HA(1)*(AT + AA*BT) NEF = NEF + HIF NET = NET + HIT H2F = H2*REF + EN*(2*H2A*HG*HGF - H2*NEF) H2T = H2*H2AT + EN*(2*H2A*HG*HGT - H2*NET) * hydrogen contribution to entropy, internal energy SIF = SIF - VA(1)*HIF - H2BT*H2F SIT = SIT - VA(1)*HIT - H2BT*H2T + H2*(ZH2T + ZH2TT) SI = SI - HI*LOG(HI/OMG(1) + EPS) - HG*LOG(HG/OMG(2) + EPS) & - H2*(LOG(H2/ZH2S + EPS) + ZH2T) UI = UI + CHI(1,1)*HI + 0.5*DH2*(HI + HG) * DB is 1 amu * number of baryons/cm3; RHO is the mass density in g/cm3 DB = EN*DE DL = LOG(DB) RHO = DB*AVM RL = LOG(RHO) RT = RET - EN*NET RF = REF - EN*NEF * second call to PRESSI compensates for spurious pressure and entropy terms DE = DB*NEO CALL PRESSI(0, TI, DE, RF, RT, F, DC, DVT, DVF, DPB, DPBT, DPBF, : DSB, DSBT, DSBF, DUB) * pressure terms PE = CB*PE TCR = T*CR P0 = TCR*DB PI = NI*P0 T4 = T*T*T*T PR = CA*T4/3.0 B = 4.0*PR/P0 PG = PE + PI + TCR*(DPA-DPB) P = PG + PR PF = (PE*PEF + PI*RF - H2F*P0 + TCR*(DPAF-DPBF))/P PT = (PE*PET + PI + PI*RT - H2T*P0 + TCR*(DPAT-DPBT) + PR*4.0)/P PL = LOG(P) * entropy, in erg/g/K DSF = NEF*DSA + NE*DSAF - NEO*DSBF - RF*B DST = NET*DSA + NE*DSAT - NEO*DSBT - (RT - 3.0)*B SF = CR*(-NI*RF + NEF*SE + NE*SEF + SIF + DSF)/AVM ST = CR*( NI*(1.5-RT) + NET*SE + NE*SET + SIT + DST)/AVM S = CR*(SE*NE + DSA*NE - DSB*NEO + B + (1.5*TL - DL + 2.5 : - LOG(CEN))*NI + SI)/AVM * internal energy, in erg/g U = TCR*(UE*NE + DUA*NE - DUB*NEO + 1.5*NI + ZH2T*H2 + 0.75*B : + TI*UI)/AVM * other thermodynamic quantities Q = PT*SF - PF*ST CP = -Q/PF GRADA = SF/Q GAMMA = Q/(RT*SF-RF*ST) ZT = SQRT(ABS((NE*REF/WF + NZZ)/NI)) * TST ought to be zero, if all the above programming is correct TST = SF/CR - P*(RT*PF-RF*PT)/(TCR*RHO) *** END OF THERMODYNAMIC CALCULATION. BEGINNING OF TABLE-BASED CALCULATIONS FR = RL/CLN10 TF = TL/CLN10 * Opacity tables from Alexander & Ferguson (1994; molecular), Itoh (1983; * electron conduction) and Iglesias & Rogers (1992; the rest) C XF = 2.0*X + Y + 1.0 C RJS -- Attempt at estimating molecular opacities C RJS 14/7/07 - Compute pressures for HG, HI and H2 IF (T.LT.1d4.AND.IMO.EQ.1) THEN C Low T's cause trouble... TFAKE = DMAX1(1d3,T) PHG = RHO*HG*(CR*TFAKE) PHI = RHO*HI*(CR*TFAKE) ! Note that P(ressure)HI spells PHI!!!! PH2 = RHO*H2*(CR*TFAKE) PC = RHO*NA(3)*(CR*TFAKE) PN = RHO*NA(4)*(CR*TFAKE) PO = RHO*NA(5)*(CR*TFAKE) C Compute dissociation coefficients - from Rossi & Maciel '83 THETA = 5040d0/TFAKE TH2 = THETA*THETA TH3 = TH2*THETA TH4 = TH3*THETA DKCN = 1d1**(1.2699d1-8.0228*THETA+1.8920d-2*TH2-6.1584e-4*TH3 : -1.1931e-6*TH4) DKCO = 1d1**(1.3590d1-1.1523d1*THETA+6.6519e-2*TH2-6.3343d-3*TH3 : +2.3781e-4*TH4) DKC2 = 1d1**(1.2654d1-6.4287*THETA+3.4061d-2*TH2-3.0069d-3*TH3 : +2.0995d-4*TH4) DKN2 = 1d1**(1.3236d1-1.0177d1*THETA+6.3664d-2*TH2-5.9114d-3*TH3 : +2.1845d-4*TH4) DKOH = 1d1**(1.2198d1-4.8641*THETA+6.6674d-2*TH2-5.6697d-3*TH3 : +1.9546d-4*TH4) DKH2O = 1d1**(2.5315d1-1.0343d1*THETA+1.0807d-1*TH2-8.7408d-3*TH3 : +2.9465d-4*TH4) C Iterate computation of these things PsubN = PN PsubC = PC PsubO = PO IF (NA(3).GT.NA(5)) THEN PsubO = 0d0 PsubC = PC ELSE PsubC = 0d0 PsubO = PO END IF PClast = 0d0 PNlast = 0d0 POlast = 0d0 ILOOP = 0 10 continue ILOOP = ILOOP + 1 PsubO = PO/(1+PsubC/DKCO+PHG/DKOH+PHG*PHG/(DKOH*DKH2O)) PsubC = PC/(1+PsubO/DKCO+PsubC/DKC2+PsubN/DKCN) PsubN = -(1+PsubC/DKCN) PsubN = PsubN + SQRT(PsubN*PsubN+4*PN/DKN2) PsubN = DKN2/2d0*PsubN deltaC = dabs((PsubC-PClast)/PsubC) deltaN = dabs((PsubN-PNlast)/PsubN) deltaO = dabs((PsubO-POlast)/PsubO) PClast = PsubC PNlast = PsubN POlast = PsubO IF (deltaC.lt.1d-3.and.deltaN.lt.1d-3.and.deltaO.lt.1d-3) THEN CONTINUE ELSE IF (ILOOP.GT.5000) THEN C write (*,*) "Convergence failed", T, TFAKE, deltaC, deltaN, deltaO ELSE GOTO 10 END IF END IF C Compute molecular partial pressures PCO = PsubC*PsubO/DKCO PC2 = PsubC*PsubC/DKC2 PN2 = PsubN*PsubN/DKN2 PCN = PsubC*PsubN/DKCN POH = PsubO*PHG/DKOH PH2O = PHG*PHG*PsubO/(DKOH*DKH2O) C Compute number abundances from pressures - note these are number fractions NCO = PCO/(CR*TFAKE*RHO) NC2 = PC2/(CR*TFAKE*RHO) NN2 = PN2/(CR*TFAKE*RHO) NCN = PCN/(CR*TFAKE*RHO) NOH = POH/(CR*TFAKE*RHO) NH2O = PH2O/(CR*TFAKE*RHO) C Fits to opacity as a function of temp: C CN per molecule - Scalo & Ulrich 1975 CAMU = 1.6605402D-24 FKCN = 1d1**(-19.12+2.22479*THETA-2.8069*TH2+0.76000*TH3-0.078384*TH4) C Assume C2 as opaque as CN FKC2 = FKCN*NC2/(CAMU) FKCN = FKCN*NCN/(CAMU) C CO, OH and H2O from Keeley 1970 FKCO = 2.75d-26*NCO/CAMU FKOH = 17*(NOH/CAMU)*(1.4d-21*(T/1d4)**6d0)/(0.1+(T/1d4)**6d0) C H2O, second term modified as in Marigo 2002 FKH2O = (2.6d-27/(4.23d-4 + (T/1d4)**4d0)) + : (9.72d-18*DEXP(-3.2552/(T/1d4))/(1+3.78d3*(T/1d4)**1d1)) FKH2O = 18*NH2O*FKH2O/CAMU C Note that I have been unable to reproduce Marigo's plots with this, CN seems C about right though and that's the important one FKTOT = FKCO+FKC2+FKCN+FKOH+FKH2O END IF if(IOP.LE.1) THEN XH = NA(1)*AM(1)/AVM XHE = NA(2)*AM(2)/AVM XF = 2.0*XH + XHE + 1.0 IF ( JX.LT.2 ) JX = 2 IF ( JX.GT.JCSX ) JX = JCSX JX = JX - 2 1 CONTINUE JX = JX + 1 IF ( XF.LT.CSX(JX+1) .AND. JX.LT.JCSX-1 ) GO TO 1 IF ( XF.GE.CSX(JX) .AND. JX.GT.1 ) THEN JX = JX - 2 GO TO 1 ELSE XT = (XF-CSX(JX))/(CSX(JX+1)-CSX(JX)) XU = 1.0 - XT IF (IOP.EQ.0) THEN * linear interpolation in opacity table TT = 20.0*(TF - 2.95) JT = MAX(1, MIN(126, INT(TT))) TT = TT - JT TU = 1.0 - TT RT = 4.0*(FR + 12.25) JR = MAX(1, MIN(89, INT(RT))) RT = MAX(0.0, RT - JR) RU = 1.0 - RT DO J = 1, 2 K = JX + J-1 IF (J .EQ. 1) FKL = TT*(RT*CS(JR+1,JT+1,K) & + RU*CS(JR,JT+1,K)) + TU*(RT*CS(JR+1,JT,K) & + RU*CS(JR,JT,K)) IF (J .EQ. 2) FKH = TT*(RT*CS(JR+1,JT+1,K) & + RU*CS(JR,JT+1,K)) + TU*(RT*CS(JR+1,JT,K) & + RU*CS(JR,JT,K)) END DO ELSE * bicubic spline interpolation CALL OPACTY(JX,TF,FR,FKL,FKH) END IF FK = XT*10D0**FKH + XU*10D0**FKL CH = 4D0*CC*PR/(FK*RHO*RHO*CP*T) * Neutrino loss rates from Itoh et al (1983-1992) EN = 0.0 IF ( TF .GE. 7.0 .AND. FR .GE. 0.0 ) THEN TT = 20.0*(TF - 6.95) IT = MAX(1, MIN(59, INT(TT))) TT = TT - IT TU = 1.0 - TT RT = 4.0*(FR + 0.25) IR = MAX(1, MIN(40, INT(RT))) RT = RT - IR RU = 1.0 - RT ENP = TT*(RT*CNIU(IT+1,IR+1,1) + RU*CNIU(IT+1,IR,1)) : + TU*(RT*CNIU(IT,IR+1,1) + RU*CNIU(IT,IR,1)) ENB = TT*(RT*CNIU(IT+1,IR+1,2) + RU*CNIU(IT+1,IR,2)) : + TU*(RT*CNIU(IT,IR+1,2) + RU*CNIU(IT,IR,2)) EN = -10.0D0**ENP - NZZ*10.0D0**ENB END IF END IF ELSE !! DO MY NEW OPACITY STUFF - JJ 5/11/02 fT=log10(exp(TL)) C RJS - note fR is now log rho/(T6^3) not log rho as it was above fR=log10(RHO)-3d0*(fT-6d0) fX = NA(1)*AM(1)/AVM fY = NA(2)*AM(2)/AVM if(IOP.eq.2) then fC=0d0 fO=0d0 else fC = (NA(3)*AM(3)/AVM-cbase)/(1d0-ZS-fX) fO = (NA(5)*AM(5)/AVM-obase)/(1d0-ZS-fX) endif if(fX.lt.1d-10.and.IOP.ge.5) then fO=1d0-fX-fY-ZS-fC endif ! add in N,Ne,Mg,Si if(fC.le.1d-6.and.fO.le.1d-6) then !Pure X only Spline Opacity FKL=0d0 FKH=0d0 if(fX.ge.0.1) then if(fX.ge.0.35) then call Xopacty(1+3*61,fT,fR,FKL,FKH) FK=((0.7-fX)*10d0**FKL+(fX-0.35)*10d0**FKH)/0.35 else call Xopacty(1+2*61,fT,fR,FKL,FKH) FK=((0.35-fX)*10d0**FKL+(fX-0.1)*10d0**FKH)/0.25 endif else if(fX.ge.0.03) then call Xopacty(1+61,fT,fR,FKL,FKH) FK=((0.1-fX)*10d0**FKL+(fX-0.03)*10d0**FKH)/0.07 else call Xopacty(1,fT,fR,FKL,FKH) FK=((0.03-fX)*10d0**FKL+(fX)*10d0**FKH)/0.03 endif endif else if(fC.lt.1d-20) fC=0d0 if(fO.lt.1d-20) fO=0d0 !Full 5D interpolate if(fX.ge.0.1) then if(fX.ge.0.35) then iX=3 else iX=2 endif else if(fX.ge.0.03) then iX=1 else iX=0 endif endif do K2=1,8 if(COcompos(K2).le.fC) iC=K2 if(COcompos(K2).le.fO) iO=K2 enddo if(iC.le.7) then JX=iO+(iC-1)*8+iX*61 else JX=6*8+7+iX*61+iO endif do im=0,1 do in=0,1 if(iO.le.6) then JK=JX+im*8+in else JK=JX+im*7+in endif call Xopacty(JK,fT,fR,FKL,FKH) FK=((Xcompos(iX+2)-fX)*10d0**FKL+(fX-Xcompos(iX+1)) : *10d0**FKH)/(Xcompos(iX+2)-Xcompos(iX+1)) tkappa(im+1,in+1)=FK enddo enddo C RJS 19/09/05 C Fudge to keep things inside the array - I have no idea what effect this has... IF (iC.GT.7) iC = 7 IF (iO.GT.7) iO = 7 et=(fC-COcompos(iC))/(COcompos(iC+1)-COcompos(iC)) eu=(fO-COcompos(iO))/(COcompos(iO+1)-COcompos(iO)) FK=(1d0-et)*(1d0-eu)*tkappa(1,1)+et*(1d0-eu)*tkappa(2,1) : +et*eu*tkappa(2,2)+(1d0-et)*eu*tkappa(1,2) endif CH = 4D0*CC*PR/(FK*RHO*RHO*CP*T) C Restore FR to log rho FR = RL/CLN10 TF = TL/CLN10 * Neutrino loss rates from Itoh et al (1983-1992) EN = 0.0 IF ( TF .GE. 7.0 .AND. FR .GE. 0.0 ) THEN TT = 20.0*(TF - 6.95) IT = MAX(1, MIN(59, INT(TT))) TT = TT - IT TU = 1.0 - TT RT = 4.0*(FR + 0.25) IR = MAX(1, MIN(40, INT(RT))) RT = RT - IR RU = 1.0 - RT ENP = TT*(RT*CNIU(IT+1,IR+1,1) + RU*CNIU(IT+1,IR,1)) : + TU*(RT*CNIU(IT,IR+1,1) + RU*CNIU(IT,IR,1)) ENB = TT*(RT*CNIU(IT+1,IR+1,2) + RU*CNIU(IT+1,IR,2)) : + TU*(RT*CNIU(IT,IR+1,2) + RU*CNIU(IT,IR,2)) EN = -10.0D0**ENP - NZZ*10.0D0**ENB END IF END IF C Add molecular opacities to tabulated ones IF (T.LT.1d4.AND.IMO.EQ.1) THEN FK = FKTOT + FK END IF RETURN END