**==FUNCS1.FOR SUBROUTINE FUNCS1(I, K, K1, K2, ISTAR) IMPLICIT REAL*8(A-H, L, M, O-Z) INTEGER MAXMSH PARAMETER (MAXMSH = 2000) COMMON H(60,MAXMSH),DH(60,MAXMSH),EPS,DEL,DH0,NMESH,JIN,IW(200) COMMON /TRANS / HT(28,MAXMSH,2) COMMON /AUXIN / ICL, ION, JW, IOP, INUC, IBC, ICN, IML(2), ISGTH, : IMO, IDIFF COMMON /SODDS / ALPHA, RML, CMG, CSI, CFE, CT(10), AGE, DT, M1, : EC, BM, ANG, CM, MTA, MTB, TM(2), T0, M0, TC(2), OS, AC, RCD, : RMG, RHL, XF, DR, AK1 ,RMT, AK2, ITH, IX, IY, IZ, IB, ISX(45), : TRB COMMON /ATDATA/ DH2(4), CCHI(26,9), OMG(27), AMASS(10), BN(10), IZZ(10) * extra common for mesh-spacing COMMON /PMESH / PMH(2), PME(2), IAGB COMMON /STAT2 / AP, ARHO, U, P, RHO, FK, T, SF, ST, ZT, GRADA, CP, : CHI, QP, PR, PG, PF, PT, EN, RPP, R33, R34, RBE, RBP, : RPC, RPN, RPO, R3A, RAC, RAN, RAO, RANE, RCC, RCO, : ROO, RGNE, RGMG, RCCG, RPNG, EX, ENX, EXX(3), AAWT, : AWT, WR(14) C NB - if new variables are added, will need to change the WI(n) numbers in the C rest of this subroutine. RJS 5/7/06 COMMON /INF / AF, AT, VX16, AM, VX1, VQK, AR, L, VX4, VX12, : VX20, VX14, HORB, HSPIN, VX3, WI(45) COMMON /DINF / DAF, DAT, DX16, DAM, DX1, V1(3), DX4, DX12, DX20, : DX14, DHORB, DHSPIN, DX3, DWI(45) COMMON /OUTF / BC1, BC2, BC3, BCHORB, BCHSPIN, VP, VPK, VR, VRK, : VT, VTK, VL, LK, : LQ, GTMA, MT, VM, VMK, QK, SG, WT, X1, X1T, X16, : X16T, X4, X4T, X12, X12T, X20, X20T, X14, X14T, : X24, SGTH, MU, X3, X3T, DLEV, DMIX(6), : D4, D12, D14, D16, D20, D3, WX(103) COMMON /ABUND / XH, XHE, XC, XN, XO, XNE, XMG, XSI, XFE, XHE3, XW(14) COMMON /OP / ZS, LEDD, M, DG, GRADT, ETH, RLF, EGR, R, Q COMMON /CNSTS / CPI, PI4, CLN10, CA, CB, CC, CD, CG, CR(2), CEVB, & CEN, CPL, CMEVMU, CSECYR, LSUN, MSUN, RSUN, TSUNYR COMMON /WT / PWT,SIG,MK,MT2,LQP,LKP COMMON /MIX / KBICZ,KICZ COMMON /MASLOS/ AIJ(6,5),baseN COMMON /ACCRET/ ACCOMPOS(7,31, 2) COMMON /EVMODE/ IMODE COMMON /INERTI/ VI(2) COMMON /MONTAG/ MIXV, MIXL COMMON /MIXFUD/ SGTHFAC, FACSGMIN, FACSG, ISGFAC COMMON /IONISE/ MSTORE(8,5), MSTORE2(5) COMMON /DIFFUS/ D(10), A12(10) DIMENSION XSPEC(6), DDMIX(10), MUXX(6) CBRT(VX) = DEXP(DLOG(VX)/3.0D0) PS(VX) = 0.5D0*(VX+DABS(VX)) RLOBE(VX) = 0.49D0*VX*VX/(0.6D0*VX*VX+DLOG(1.0D0+VX)) M = DEXP(AM) C Spin period to equns - this isn't the best way of doing this... BCHSPIN = HSPIN C Compute binary mass if in binary mode IF (IMODE.EQ.2) THEN BM = M + DEXP(WI(4)) END IF VM3 = CBRT(M) C Pierre's mass mod. MT = (M - DEXP(AM-DAM))/DT C MT = M*DAM/DT MT2 = MT C set up the composition variables XH = VX1 XHE = VX4 XC = VX12 XN = VX14 XO = VX16 XNE = VX20 XMG = CMG*ZS XSI = CSI*ZS XFE = CFE*ZS XHE3 = VX3 * Put everything else in Mg24. XMG = 1.0D0 - XH - XHE - XC - XO - XNE - XN - XSI - XFE - XHE3 IF (XMG .LT. 0.0D0) XMG = 0.0D0 * following if N14 is a variable rather than Ne20: C XNE = 1.0D0 - XH - XHE - XC - XN - XO - XMG - XSI - XFE C IF (XNE .LT. 0.0D0) XNE = 0.0D0 CALL STATEL (I, AF, AT, ISTAR) R = DEXP(AR) R2 = R*R R2M = 2.0D0/(PI4*RHO*R) C pressure gradient equation; gravity corrected for corotation C ANGMOM is the orbital angular momentum GRAV = 1.0D11*CG*M/R2 SILLYM = DMAX1(M,0.95D0*TM(ISTAR)) IF (IB.EQ.1) GO TO 1 FDT=AGE*3.1557D7+DT-T0 M1=M0-PS(-FDT)*MTA+PS(FDT)*MTB 1 SEP = BM*(ANG/(SILLYM*(BM-SILLYM)))**2 C Angular momentum in G=1 units, according to RPC C Replacing old separation with new one based on eigenvalue HORB IF (ISTAR.EQ.2) HORB = WI(13) IF (IMODE.EQ.2) THEN SEP = (TM(1)+TM(2))*(HORB/(TM(1)*TM(2)))**2.0 ELSE SEP = BM*(HORB/(TM(1)*(BM - TM(1))))**2.0 END IF * GRAV = GRAV - GRAV*BM*R*R2/(1.5D0*M*SEP**3) C gravity corrected for rotation OSPIN = (H(14+15*(ISTAR - 1),1) + DH(14+15*(ISTAR - 1),1))/VI(ISTAR) C Need to mess about with units to get this right OSPIN = OSPIN*SQRT(CG) ! now in s^-1 GRAV = GRAV - 1d11*2.0/3.0*OSPIN**2.0*R IF (1.0D11*CG*M/R2.LT.1d11*2.0/3.0*OSPIN**2.0*R.AND.I.EQ.-1) THEN WRITE (*,*) "break-up velocity reached!" write (*,*) ISTAR, 1.0D11*CG*M/R2, 1d11*2.0/3.0*OSPIN**2.0*R, VI(ISTAR) write (*,*) OSPIN, R, H(14+15*(ISTAR - 1),1) + DH(14+15*(ISTAR - 1),1) STOP END IF APM = -1.0D11*GRAV/(PI4*R2*P) C temperature gradient equation GRADR = 1.0D11*FK*P*L/(4.0D0*PI4*CC*GRAV*R2*PR) DG = GRADR - GRADA S2 = GRADA*CP*T HP = P/(GRAV*RHO) C mixing-length theory W0 = 0.5D0*S2*(ALPHA*ALPHA*HP/(9.0D0*CHI))**2 W2 = 546.75D0*W0*DMAX1(0.0D0,DG) + 73.0D0 W3 = CBRT(W2+DSQRT(W2*W2+12167.0D0)) W5 = DMAX1(W3-23.0D0/W3-2.0D0,0.0D0) WCV = W5*CHI/(3.0D0*ALPHA*HP) WL = ALPHA*HP*WCV MIXL = ALPHA*HP MIXV = WCV IF (DG.LT.0d0) MIXV = 0d0 C if GRADR < GRADA, we get WCV = 0 and GRADT = GRADR IF (KBICZ.GT.600.AND.K.GT.KBICZ.AND.IAGB.EQ.1) WCV = 0d0 GRADT = GRADR-4.0D0*HP*WCV**3/(ALPHA*S2*CHI) GTMA = GRADT - GRADA ATM = GRADT*APM C Old MSF 2000 CT1 = 1.0D1**(1.0D1*CT(1)) !!! should change input format !!! VP = CT(4)*AP + CT(5)*LOG((P+CT(9))/(P+CT1)) + & CT(2)*LOG((P+CT(10))/(P+CT1)) VPP = CT(4) + P/(P+CT1)*(CT(5)*(CT1-CT(9))/(P+CT(9)) & + CT(2)*(CT1-CT(10))/(P+CT(10))) CT10 = 2.0D4 !!! fixed, should change input format !!! VT = CT(7)*DLOG(T/(T+CT10)) VTT = CT(7)*CT10/(T+CT10) C mesh spacing equation, with modified pressure, temperature gradient equations IF (IAGB.EQ.1) THEN CP1 = CT(1)*PME(ISTAR) CP2 = CT(10)*PME(ISTAR) CP3 = CT(9)*PMH(ISTAR) VP = CT(4)*AP + CT(5)*LOG((P+CP3)/(P+CP1)) + & CT(2)*LOG((P+CP2)/(P+CP1)) VPP = CT(4) + P/(P+CP1)*(CT(5)*(CP1-CP3)/(P+CP3) & + CT(2)*(CP1-CP2)/(P+CP2)) CT10 = 2.0D4 !!! fixed, should change input format !!! VT = CT(7)*DLOG(T/(T+CT10)) VTT = CT(7)*CT10/(T+CT10) END IF C VR and VM must go to zero at centre like r**2, m**(2/3) VM = CT(6)*CBRT(TM(ISTAR)*TM(ISTAR)) VMM = VM + VM3*VM3 VM = DLOG(VM/VMM) VMM = -0.66666667D0/(VMM*VM3) VR = -CT(3)*DLOG(R2/CT(8)+1.0D0) R21 = (3D0*M/(PI4*RHO))**(2D0/3D0) VR1 = -CT(3)*DLOG(R21/CT(8)+1.0D0) VRR = -CT(3)/(R2+CT(8)) C Q is the quantity that the meshpoints should be at equal intervals of Q = VP + VT + VM + VR QK = VQK QM = VPP*APM + VTT*ATM + VMM + VRR*R2M MK = VQK/QM VMK = VMM*MK VPK = VPP*APM*MK VTK = VTT*ATM*MK VRK = VRR*R2M*MK WT = 1D-7*R2*CHI*DT/(VRK*VRK) PWT = WT C energy equation VL = L EN = EN + ENX *** assume thermal eq. in central region if timestep small * ITHC = ITH * IF(K.GT.19*NMESH/20 .AND. DT .LT. 1D-1*CSECYR) ITHC = 0 * LK = (EX + EN + EC - ITHC*T*(SF*DAF+ST*DAT)/DT)*MK * LQ = ITHC*CP*T*APM*MK *** weight by a function of WT that drops from 1 to 0 for WT between 10^10 *** and 10^9 IF (K.GT.39*NMESH/40) THEN AWT = 2D-16*R2*CHI*DT/(R2M*MK)**2 AWT = 2D-15*R2*CHI*DT/(R2M*MK)**2 AWT4 = AWT**4 AAWT = AWT4/(1D0 + AWT4) ELSE AAWT = 1.0D0 END IF LK = (EX + EN + EC - ITH*T*(SF*DAF+ST*DAT)/DT)*MK LKP = EX + EN + EC - ITH*T*(SF*DAF+ST*DAT)/DT SIG = SF*DAF LQ = ITH*CP*T*APM*MK LQP = CP*T*APM*GTMA L1 = M*(EX + EN + EC - ITH*T*(SF*DAF+ST*DAT)/DT) C radius equation R2K = MK/(0.5D0*PI4*RHO*R) C composition equations * simplified to exclude reactions involving Mg24 C hydrogen equation with MS baryon correction when ICN = 1 X1 = VX1 * X1T = (2.0*IX*((1-ICN)*RPP + RPC + RPNG + (1-2*ICN)*(RPN + RPO)) X1T = (IX*((1-ICN)*RPP*3 + RPC*2 + RPNG*2 -R33*2 +R34 +2* (1-2*ICN)*(RPN + RPO)) * X1T = (IX*((1-ICN)*RPP*3 + RPC*2 + RPNG*2 -R33*2 +R34 +2* (1-ICN)*(RPN + RPO)) : + DX1/DT)*MK * x1t = (-5d-7*(0.76d0-x1)/csecyr + dx1/dt)*mk C helium equation X4 = VX4 * X4T = (4.0*(-IX*(0.5*RPP + RPN + RPO)*(1-ICN) X4T = (4.0*(-IX*(RPN + RPO + R33 +R34)*(1-ICN) * : + IY*(3.0*R3A + RAC + 1.5*RAN + RAO + RANE) * : - IZ*(RCC + RCO + 2.0*ROO + RGNE + RGMG)) + DX4/DT)*MK : + IY*(3.0*R3A + RAC + 1.5*RAN + RAO) : - IZ*(RCC + RGNE)) + DX4/DT)*MK * x4t = (5d-7*(0.76-x1)/csecyr + dx4/dt)*mk C carbon equation X12 = VX12 X12T = (12.0*(IX*(RPC - RPN) - IY*(R3A - RAC) * : + IZ*(2.0*(RCC + RCCG) + RCO)) + DX12/DT)*MK : + IZ*2.0*RCC) + DX12/DT)*MK C nitrogen equation X14 = VX14 X14T = (14.0*(IX*(RPN + RPNG - RPC - RPO) + IY*RAN) + DX14/DT)*MK C oxygen equation X16 = VX16 X16T = (16.0*(IX*(RPO - RPNG) - IY*(RAC - RAO) * : + IZ*(RCO + 2.0*ROO - RGNE)) + DX16/DT)*MK : - IZ*RGNE) + DX16/DT)*MK C neon equation X20 = VX20 * X20T = (20.0*(IY*(RANE - RAN - RAO) + IZ*(RGNE - RGMG - RCC)) X20T = (20.0*(-IY*(RAN + RAO) + IZ*(RGNE - RCC)) : + DX20/DT)*MK C He3 equation X3 = VX3 X3T = (3.0*(-IX*RPP*(1-ICN)+IX*(2d0*R33+R34))+DX3/DT)*MK C fudged convective diffusion coefficient: OS > 0 gives overshooting B = PR/PG AMAP=0.1D0/DABS(APM*M) EG = DG+OS/(2.5D0+B*(2.0D1+1.6D1*B))/(AMAP+1.0D0) UG = PS(EG) SG = UG*UG*TC(ISTAR)/(QM*QK) SIG = SG C modified diffusion coefficient according to MLT SGMLT = (WL*1.0D-22*(PI4*RHO*R2)**2)/(3.0*MK) VG = UG/GRADR IF (IAGB.EQ.1) THEN SG = SGMLT IF (XH.GT.1.0D-6) THEN ! H-rich envelope... SG = SGMLT * AK1*VG/(1.0 - (1.0-AK1)*VG) ELSE ! interior... SG = SGMLT * AK2*VG/(1.0 - (1.0-AK2)*VG) END IF END IF IF (KBICZ.GT.600.AND.K.GT.KBICZ.AND.IAGB.EQ.1) SG = 0d0 C Mixing fudge SG = FACSG*SG C Compute mean molecular weight, inc. for partial ionisation C Should have IOP equal to 5 XH1 = H(5+15*(ISTAR - 1),K) XHE1 = H(9+15*(ISTAR - 1),K) XC1 = H(10+15*(ISTAR - 1),K) XN1 = H(12+15*(ISTAR - 1),K) XO1 = H(3+15*(ISTAR - 1),K) XNE1 = H(11+15*(ISTAR - 1),K) XHE31 = H(15+15*(ISTAR - 1),K) LKP = 0d0 XSPEC(1) = XH XSPEC(2) = XHE XSPEC(3) = XC XSPEC(4) = XN XSPEC(5) = XO XSPEC(6) = XHE3 DO II = 1,ION MUX = 1d0 DO JJ = 1, IZZ(II) MUX = MUX + JJ*MSTORE(JJ,II)/MSTORE2(II) END DO C Store mu of each species MUXX(II) = 1/(MUX/AMASS(II)) IF (II.EQ.2) MUXX(6) = 1/(MUX/AMASS(10)) LKP = LKP + MUX*XSPEC(II)/AMASS(II) C pretend He3 is like He4... IF (II.EQ.2) LKP = LKP + MUX*XSPEC(6)/AMASS(10) END DO LKP = 1/LKP MU = LKP MUREAL = MU C diffusion requires the use of the actual mu, not one based on full C ionisation. Note that it is the actual mu that is passed to printb C Need to test what this does to the TH boundary output stuff... IF (IDIFF.EQ.1) THEN C Paquette diffusion stuff CALL DIFFUSION(RHO,T) DO II = 1, 10 C can't diffuse through H if we don't have any... IF (XH.LT.0.1) D(II) = 0d0 DDMIX(II) = D(II)*(1.0D-11*(PI4*RHO*R2))**2d0/MK A12(II) = A12(II)*D(II)*ATM*(1.0D-11*(PI4*RHO*R2)) IF (II.LE.5) THEN D(II) = D(II)*GRAV/(CR(1)*T)*(MUXX(II) - MU) ELSE D(II) = D(II)*GRAV/(CR(1)*T)*(AMASS(II)/(1d0+IZZ(II)) - MU) IF (II.EQ.10) D(II) = D(II)*GRAV/(CR(1)*T)*(MUXX(6) - MU) END IF D(II) = D(II) - A12(II) D(II) = D(II)*(1.0D-11*(PI4*RHO*R2)) END DO D4 = D(2) D12 = D(3) D14 = D(4) D16 = D(5) D20 = D(6) D3 = D(10) DMIX(1) = DDMIX(2) DMIX(2) = DDMIX(3) DMIX(3) = DDMIX(4) DMIX(4) = DDMIX(5) DMIX(5) = DDMIX(6) DMIX(6) = DDMIX(10) ELSE DO II = 1,6 DMIX(II) = 0d0 END DO D4 = 0d0 D12 = 0d0 D14 = 0d0 D16 = 0d0 D20 = 0d0 D3 = 0d0 END IF IF (ISGTH.EQ.1) THEN MU = 1.0/(2.0*XH1 + 0.75*XHE1 + 7.0/12.0*XC1 + 8.0/14.0*XN1 + 9.0/16.0*XO1 : + 11.0/20.0*XNE1 + XHE31) ! + 13.0/24.0*XMG + 15.0/28.0*XSI + 27.0/56.0*XFE) C Full implicit mu -- is much less stable! C MU = 1.0/(2.0*XH + 0.75*XHE + 7.0/12.0*XC + 8.0/14.0*XN + 9.0/16.0*XO C : + 11.0/20.0*XNE + XHE3) ! + 13.0/24.0*XMG + 15.0/28.0*XSI + 27.0/56.0*XFE) CHI1 = 4.0D0*CC*PR/(FK*RHO*RHO*CP*T) DTH = 3.0D0*CHI1/(MU*(-GTMA)*APM*MK) DTH = DMAX1(0d0,DTH) SGTH = SGTHFAC*DTH*1.0D-22*(PI4*R2*RHO)**2/MK C IF (XHE31.LT.1d-12.AND.K.LT.400) SGTH = 0d0 IF (K.LE.4) THEN SGTH = 0d0 END IF ELSE SGTH = 0d0 END IF MT2 = SGTH C the outermost 25 meshpoints are always mixed, even if not convective; C IF ( K.LT.5 ) SG = TC(ISTAR)/(QM*QK) IF ( K.LT.0.03*NMESH ) SG = 1d2*TC(ISTAR)/(QM*QK) IF ( K.GT.NMESH-0.01*NMESH) SG = TC(ISTAR)/(QM*QK) IF ( I.LE.0 ) THEN C sundry numbers for PRINTB, FUNCS2 ETH = -T*(SF*DAF+ST*DAT)/DT + CP*T*APM*GTMA*MT EGR = 1.0D22*CG*M/R - U LEDD = PI4*CC*CG*M/FK HT(1, K, ISTAR) = RHO HT(2, K, ISTAR) = SG HT(3, K, ISTAR) = ZT HT(4, K, ISTAR) = MK DO J = 1,14 HT(4+J, K, ISTAR) = XW(J) END DO C Store thermohaline stuff HT(19, K, ISTAR) = SGTH HT(20, K, ISTAR) = MU C Store temp and ln f HT(21, K, ISTAR) = AF HT(22, K, ISTAR) = AT C HT(23-24,K) are for mass loss, computed in massloss.f C Store diffusion coeff related stuff HT(25, K, ISTAR) = GRAV/(CR(1)*T) HT(26, K, ISTAR) = 1.0D-11*(PI4*RHO*R2) HT(27, K, ISTAR) = MUREAL HT(28, K, ISTAR) = ATM END IF C Calc RLF for printb - only if I=-1 IF (K.EQ.1.AND.I.EQ.-1) THEN C surface boundary conditions IF (IMODE.EQ.2) THEN BM = TM(1) + TM(2) END IF RAT = M/(BM-M) RLF = AR - DLOG(SEP*RLOBE(CBRT(RAT))) IF (AGE.LT.1d4) THEN RLF = -1d-1 END IF END IF IF ( K.NE.K1 ) RETURN PRB = CA*TRB**4/3D0 BC2 = DLOG(FK/GRAV*(1.5D0*PG + 0.75D0*PR)*(PR-2D0*PRB)/PR) BC3 = 1.0D11*L - (0.75D0*PI4*CC*R2*(PR-2D0*PRB)) IF (I.GE.0) THEN IF (ISTAR.EQ.1) ISTAROTHER = 2 IF (ISTAR.EQ.2) ISTAROTHER = 1 ACCOMPOS(1,I+1,ISTAR) = 0.5*(H(5+15*(ISTAR-1),1) + H(5+15*(ISTAROTHER-1),1)) !X1 ACCOMPOS(2,I+1,ISTAR) = 0.5*(H(9+15*(ISTAR-1),1) + H(9+15*(ISTAROTHER-1),1)) !X4 ACCOMPOS(3,I+1,ISTAR) = 0.5*(H(10+15*(ISTAR-1),1) + H(10+15*(ISTAROTHER-1),1)) !X12 ACCOMPOS(4,I+1,ISTAR) = 0.5*(H(12+15*(ISTAR-1),1) + H(12+15*(ISTAROTHER-1),1)) !X14 ACCOMPOS(5,I+1,ISTAR) = 0.5*(H(3+15*(ISTAR-1),1) + H(3+15*(ISTAROTHER-1),1)) !X16 ACCOMPOS(6,I+1,ISTAR) = 0.5*(H(11+15*(ISTAR-1),1) + H(11+15*(ISTAROTHER-1),1)) !X20 ACCOMPOS(7,I+1,ISTAR) = 0.5*(H(15+15*(ISTAR-1),1) + H(15+15*(ISTAROTHER-1),1)) END IF RETURN END