**==REMESH.FOR SUBROUTINE REMESH(M, NCH, CH, CO, CC, CNE ) IMPLICIT REAL*8(A-H, O-Z) PARAMETER (MAXMSH = 2000) COMMON H(60,MAXMSH),DH(60,MAXMSH),EP(3),NMESH,JIN,ID(200) COMMON /NUCMAT/ HNUC(100,MAXMSH), DHNUC(100,MAXMSH) COMMON /SODDS / ALPHA, RML, CMG, CSI, CFE, CT(10), AGE, DT, AM1, : EC, BM, ANG, CM, AMTA, AMTB, TM(2), T0, AM0, TC(2), OS, AC, RCD, : RMG, RHL, XF, DR, AK1 ,RMT, AK2, IZ(4), IB, ISX(45), : TRB COMMON /INF / VAR(60) COMMON /DINF / DVAR(60) COMMON /OP / ZS, W(8), QQ COMMON /ATDATA/ DH2(4), CHI(26,9), OMG(27), AM(10), BN(10), JZ(10) COMMON /YUK / Q(MAXMSH) COMMON /EVMODE/ IMODE COMMON /CNSTS / CPI, PI4, CLN10, CA, CB, CLIGHT, CD, CG, CR(2), CEVB, & CEN, CPL, CMEVMU, CSECYR, CLSUN, CMSUN, RSUN, TSUNYR COMMON /INERTI/ VI(2) COMMON /ANGMOM/ VROT1, VROT2, FMAC, FAM, IRAM, IRS1, IRS2 write (32,*) "Remeshing model" C RJS 22/11/07 - Save NMESH, we need it for the second remesh pass NMESHORIG = NMESH DO ISTAR = 1,IMODE!2 NMESH = NMESHORIG C Calculate moment of inertia VI(ISTAR) = 0d0 DO K=1,NMESH - 1 VI(ISTAR) = VI(ISTAR) + 2.0/3.0*(DEXP(H(4+15*(ISTAR-1),K)) - : DEXP(H(4+15*(ISTAR-1),K+1)))*DEXP(H(7+15*(ISTAR-1),K))**2.0 END DO VI(ISTAR) = VI(ISTAR) + 2.0/5.0*DEXP(H(4+15*(ISTAR-1),NMESH))*DEXP(H(7+15*(ISTAR-1),NMESH))**2.0 IF ( NCH.LT.3 ) GO TO 41 * * Set initial composition. * The composition variables are NOT the actual mass fractions if we use * non-integer atomic masses, so we have to compute what they are * CHE = 1.0 - CH - ZS CN = 1D0 - CC - CO - CNE - CMG - CSI - CFE XH = CH*BN(1)/AM(1) XHE = CHE*BN(2)/AM(2) XC = CC*ZS*BN(3)/AM(3) XN = CN*ZS*BN(4)/AM(4) XO = CO*ZS*BN(5)/AM(5) XNE = CNE*ZS*BN(6)/AM(6) C should do He3 like this... XMG = CMG*ZS XSI = CSI*ZS XFE = CFE*ZS VMA = XH + XHE + XC + XN + XO + XNE + XMG + XSI + XFE C Minor element initial abundances from Anders & Grevesse 89 XD = 4.801d-5*(ZS/0.02) XHE3 = 2.929d-5*(ZS/0.02) XLI = 9.353d-9*(ZS/0.02) C No data available for Be-7, assume it's zero? XBE = 0d0 XB11 = 4.725d-9*(ZS/0.02) XC13 = 3.650d-5*(ZS/0.02) XC14 = 0d0 XN15 = 4.363d-6*(ZS/0.02) XO17 = 3.887d-6*(ZS/0.02) XO18 = 2.167d-5*(ZS/0.02) XF19 = 4.051d-7*(ZS/0.02) XNE21 = 4.127d-6*(ZS/0.02) XNE22 = 1.302d-4*(ZS/0.02) XNA22 = 0d0 XNA23 = 3.339d-5*(ZS/0.02) XMG24 = 5.148d-4*(ZS/0.02) XMG25 = 6.766d-5*(ZS/0.02) XMG26 = 7.760d-5*(ZS/0.02) X26M = 0d0 X26G = 0d0 XAL27 = 5.798d-5*(ZS/0.02) XSI28 = 6.530d-4*(ZS/0.02) XSI29 = 3.426d-5*(ZS/0.02) XSI30 = 2.352d-5*(ZS/0.02) XP31 = 8.155d-5*(ZS/0.02) XS32 = 3.958d-4*(ZS/0.02) XS33 = 3.222d-6*(ZS/0.02) XS34 = 1.866d-5*(ZS/0.02) XFE56 = 1.169d-3*(ZS/0.02) XFE57 = 2.855d-5*(ZS/0.02) XFE58 = 3.697d-6*(ZS/0.02) XFE59 = 0d0 XFE60 = 0d0 XCO59 = 3.358d-6*(ZS/0.02) XNI58 = 4.944d-5*(ZS/0.02) XNI59 = 0d0 XNI60 = 1.958d-5*(ZS/0.02) XNI61 = 8.594d-7*(ZS/0.02) DO 4 K = 1, NMESH H(5+15*(ISTAR-1), K) = XH/VMA H(3+15*(ISTAR-1), K) = XO/VMA H(9+15*(ISTAR-1), K) = XHE/VMA H(10+15*(ISTAR-1),K) = XC/VMA H(11+15*(ISTAR-1),K) = XNE/VMA H(12+15*(ISTAR-1),K) = XN/VMA C HE3 for structure code H(15+15*(ISTAR - 1),K) = XHE3/VMA C Beginning of minor variables, in mass order (except gallinoes) HNUC(1+50*(ISTAR-1),K) = 0d0 ! gallinoes HNUC(2+50*(ISTAR-1),K) = 0d0 ! neutrons HNUC(3+50*(ISTAR-1),K) = XD/VMA HNUC(4+50*(ISTAR-1),K) = XHE3/VMA HNUC(5+50*(ISTAR-1),K) = XLI/VMA HNUC(6+50*(ISTAR-1),K) = XBE/VMA HNUC(7+50*(ISTAR-1),K) = XB11/VMA HNUC(8+50*(ISTAR-1),K) = XC13/VMA HNUC(9+50*(ISTAR-1),K) = XC14/VMA HNUC(10+50*(ISTAR-1),K) = XN15/VMA HNUC(11+50*(ISTAR-1),K) = XO17/VMA HNUC(12+50*(ISTAR-1),K) = XO18/VMA HNUC(13+50*(ISTAR-1),K) = XF19/VMA HNUC(14+50*(ISTAR-1),K) = XNE21/VMA HNUC(15+50*(ISTAR-1),K) = XNE22/VMA HNUC(16+50*(ISTAR-1),K) = XNA22/VMA HNUC(17+50*(ISTAR-1),K) = XNA23/VMA HNUC(18+50*(ISTAR-1),K) = XMG24/VMA HNUC(19+50*(ISTAR-1),K) = XMG25/VMA HNUC(20+50*(ISTAR-1),K) = XMG26/VMA HNUC(21+50*(ISTAR-1),K) = X26M/VMA HNUC(22+50*(ISTAR-1),K) = X26G/VMA HNUC(23+50*(ISTAR-1),K) = XAL27/VMA HNUC(24+50*(ISTAR-1),K) = XSI28/VMA HNUC(25+50*(ISTAR-1),K) = XSI29/VMA HNUC(26+50*(ISTAR-1),K) = XSI30/VMA HNUC(27+50*(ISTAR-1),K) = XP31/VMA HNUC(28+50*(ISTAR-1),K) = XS32/VMA HNUC(29+50*(ISTAR-1),K) = XS33/VMA HNUC(30+50*(ISTAR-1),K) = XS34/VMA HNUC(31+50*(ISTAR-1),K) = XFE56/VMA HNUC(32+50*(ISTAR-1),K) = XFE57/VMA HNUC(33+50*(ISTAR-1),K) = XFE58/VMA HNUC(34+50*(ISTAR-1),K) = XFE59/VMA HNUC(35+50*(ISTAR-1),K) = XFE60/VMA HNUC(36+50*(ISTAR-1),K) = XCO59/VMA HNUC(37+50*(ISTAR-1),K) = XNI58/VMA HNUC(38+50*(ISTAR-1),K) = XNI59/VMA HNUC(39+50*(ISTAR-1),K) = XNI60/VMA HNUC(40+50*(ISTAR-1),K) = XNI61/VMA !X26T/VMA HNUC(41+50*(ISTAR-1),K) = XH/VMA HNUC(42+50*(ISTAR-1),K) = XHE/VMA HNUC(43+50*(ISTAR-1),K) = XC/VMA HNUC(44+50*(ISTAR-1),K) = XN/VMA HNUC(45+50*(ISTAR-1),K) = XO/VMA HNUC(46+50*(ISTAR-1),K) = XNE/VMA 4 CONTINUE 41 CONTINUE C Orbital stuff doesn't work properly - copy surface point to all others DO K=2,NMESH H(13,K) = H(13,1) H(14,K) = H(14,1) H(29,K) = H(29,1) END DO C Reset orbital AM using period in modin IF (IRAM.EQ.1.AND.ISTAR.EQ.1) THEN write (32,*) "Resetting orbital AM" DO K=1,NMESH H(13,K) = ANG END DO END IF IF (IRS1.EQ.1.AND.ISTAR.EQ.1) THEN write (32,*) "Resetting *1 spin" VROT = VROT1*1d3 OMEGA = VROT/(1d9*DEXP(H(7,1))) C Convert omega to code units OMEGA = OMEGA/DSQRT(CG) DO K=1,NMESH H(14,K) = VI(1)*OMEGA END DO END IF IF (IRS2.EQ.1.AND.ISTAR.EQ.2) THEN write (32,*) "Resetting *2 spin" VROT = VROT2*1d3 OMEGA = VROT/(1d9*DEXP(H(22,1))) OMEGA = OMEGA/DSQRT(CG) DO K=1,NMESH H(29,K) = VI(2)*OMEGA END DO END IF QE = 0D0 DO K = 1, NMESH DO J = 1, 15 !60 DVAR(J) = 0.0D0 VAR(J) = H(J+15*(ISTAR-1), K) END DO VAR(6) = 1.0D0 CALL FUNCS1(-1, K, 1, NMESH, ISTAR) Q(K) = QQ END DO PHI = Q(NMESH) - Q(1) Q1 = Q(1) Q2 = (M-1.0D0)/PHI IF ( NCH.GE.2 ) THEN C IF REQUIRED, REDISTRIBUTE THE MESH POINTS, EITHER BY CHANGING THEIR NUMBER C OR BY CHANGING THE MESH-SPACING FUNCTION. DO K = 1, NMESH VX = H(5+15*(ISTAR-1), K) + 1.0D-10 H(5+15*(ISTAR-1), K) = DLOG(VX) Q(K) = (Q(K)-Q1)*Q2 + 1.0D0 END DO KK = 1 DO K = 1, M DK = 0.0D0 IF ( K.EQ.M ) KK = NMESH IF ( K.NE.1 .AND. K.NE.M ) THEN DO I = 1, 50 IF ( K.GE.Q(KK+1) ) KK = KK + 1 IF ( K.LT.Q(KK+1) ) GO TO 10 END DO 10 CONTINUE DK = (K-Q(KK))/(Q(KK+1)-Q(KK)) END IF DO J = 1, 15 !60 DH(J+15*(ISTAR-1), K) = H(J+15*(ISTAR-1), KK) + : DK*(H(J+15*(ISTAR-1),KK+1)-H(J+15*(ISTAR-1),KK)) END DO C Remesh nucleosynthesis matrix DO J=1,50 DHNUC(J+50*(ISTAR-1),K) = HNUC(J+50*(ISTAR-1), KK) + : DK*(HNUC(J+50*(ISTAR-1), KK+1)-HNUC(J+50*(ISTAR-1), KK)) END DO END DO NMESH = M DO K = 1, NMESH DH(5+15*(ISTAR-1), K) = DEXP(DH(5+15*(ISTAR-1),K)) - 1.0D-10 C Is this line still necessary? IF ( DH(5+15*(ISTAR-1),K).LT.1.0D-5 ) DH(5+15*(ISTAR-1), K) = 0.0D0 DO J = 1, 15 !60 H(J+15*(ISTAR-1), K) = DH(J+15*(ISTAR-1), K) DH(J+15*(ISTAR-1), K) = 0.0D0 END DO DO J=1,50 HNUC(J+50*(ISTAR-1), K) = DHNUC(J+50*(ISTAR-1),K) DHNUC(J+50*(ISTAR-1),K) = 0d0 END DO END DO END IF DO K = 1, NMESH H(6+15*(ISTAR-1), K) = 1.0D0/Q2 END DO END DO RETURN END