**==NEUTRON.F C Sorts out where n's get used up - I hope! C NEEDS TO BE OVERHAULED! SUBROUTINE NEUTRON(ISTAR) IMPLICIT REAL*8(A-H, O-Z) PARAMETER (MAXMSH = 2000) COMMON /NEUTRO/ FRAC(50,MAXMSH) COMMON H(60,MAXMSH),DH(60,MAXMSH),EPS,V(2),NMESH,JIN,ID(100),IE(100) COMMON /NUCMAT/ HNUC(100,MAXMSH), DHNUC(100,MAXMSH) COMMON /SODDS / W5(2), CMG, CSI, CFE, W8(11), DT, W6(23), IW(50), TRB COMMON /STAT2 / W1(19), RPP, R33, R34, RBE, RBP, : RPC, RPN, RPO, R3A, RAC, RAN, RAO, RANE, RCC, RCO, : ROO, RGNE, RGMG, RCCG, RPNG, WR(21) COMMON /RATES / W7(1),RPD,RPBE,RPLI,REBE,RPC13,RPN15A,RPO17a, : RAC13n, RAN15g, RAO17n, RLIA, RPB11, RPC14, RAC14, : RPO18,RPO18a,RAO18,RAO18n,RPF19,RPF19A,RAF19,RPNE21, : RANE21,RANE21n,RPNE22,RANE22,RANE22n,RnNA22p,RnNA22a,RPNA22, : RPNA23,RPNA23n,RPNA23a,RPMG24,RAMG24,RPMG25T,RPMG25M,RPMG25G, : RAMG25,RAMG25n,RAMG25p,RPMG26,RPMG26Tn,RPMG26Mn,RPMG26Gn, : RAMG26,RAMG26n,RAMG26p,Rg26Tp,Rn26Tp,Rn26Ta,Rp26T,Rg26Mp, : Rn26Mp,Rn26Ma,Rp26M,Rg26Gp,Rn26Gp,Rn26Ga,Rp26G,RPAL27, : RPAL27a,RAAL27n,RANA23nT,RANA23nM,RANA23nG,RPSI28,RPSI29, : RPSI30,RPN15,RPO17, RPNE20 COMMON /ABUND / XA(10), WW(14) COMMON /ABUND2/ XA2(50), WW2(50) COMMON /OP / ZS, VLEDD, VM, GR, GRAD, ETH, RLF, EGR, R, QQ DIMENSION BARYN(50) DIMENSION PROTONS(MAXMSH) DATA BARYN/62.0, 1.0, 2.0, 3.0, 7.0, 7.0, 11.0, 13.0, 14.0, 15.0, : 17.0, 18.0, 19.0, 21.0, 22.0, 22.0, 23.0, 24.0, 25.0, 26.0, : 26.0,26.0,27.0,28.0,29.0,30.0,31.0,32.0,33.0,34.0,56.0,57.0, : 58.0,59.0,60.0,59.0,58.0,59.0,60.0,61.0,1.0,4.0,12.0,14.0,16.0, : 20.0,1.0,1.0,1.0,1.0/ DO K=1,NMESH DO I=1,50 IF (I.NE.2.AND.I.NE.41) THEN C Don't change proton and neutron fractions - they are dealt with below DHNUC(I+50*(ISTAR-1),K) = DHNUC(I+50*(ISTAR-1),K) - : FRAC(I,K)*BARYN(I)*DHNUC(2+50*(ISTAR-1),K) END IF END DO PROTONS(K) = - FRAC(41,K)*BARYN(41)*DHNUC(2+50*(ISTAR-1),K) C DH(53,K) = DH(53,K) + FRAC(41,K)*BARYN(41)*DH(14,K) DHNUC(2+50*(ISTAR-1),K) = 0d0 END DO C Sort out remaining protons - RJS 27/8/04 DO K=1,NMESH C Call to statel to get electron density & no. abundances. AF = H(1+15*(ISTAR-1),K) + DH(1+15*(ISTAR-1),K) AT = H(2+15*(ISTAR-1),K) + DH(2+15*(ISTAR-1),K) CALL STATEL (0, AF, AT, ISTAR) C Blank any -ve abundances. There shouldn't be any anyway... C RJS 19/1/04 XA(1) = H(5+15*(ISTAR-1),K) + DH(5+15*(ISTAR-1),K) XA(2) = H(9+15*(ISTAR-1),K) + DH(9+15*(ISTAR-1),K) XA(3) = H(10+15*(ISTAR-1),K) + DH(10+15*(ISTAR-1),K) XA(4) = H(12+15*(ISTAR-1),K) + DH(12+15*(ISTAR-1),K) XA(5) = H(11+15*(ISTAR-1),K) + DH(11+15*(ISTAR-1),K) XA(6) = H(3+15*(ISTAR-1),K) + DH(3+15*(ISTAR-1),K) XA(7) = CMG*ZS XA(8) = CSI*ZS XA(9) = CFE*ZS DO I = 1, 10 IF (XA(I).LT.0d0) THEN XA(I) = 0d0 WW(I) = 0d0 END IF END DO DO I = 1, 50 XA2(I) = HNUC(I+50*(ISTAR-1),K) + DHNUC(I+50*(ISTAR-1),K) IF (XA2(I).LT.0d0) XA2(I) = 0d0 XA2(2) = 0d0 END DO C Pipe excess protons to nucrat as only source of them XA2(41) = PROTONS(K) C Work out abundances for nucrat2 - RJS DO I = 1, 50 WW2(I) = XA2(I)/BARYN(I) END DO CALL NUCRAT2(AT) RATETOT = 2d0*RPP + RPC + RPN + RPO + RPNG + RPD + RPBE + RPLI : + RPC13 + RPN15A + RPO17a + RPB11 + RPC14 + RPO18 + RPO18a : + RPF19 + RPF19A + RPNE21 + RPNE22 + RPNA22 + RPNA23 : + RPNA23n + RPNA23a + RPMG24 + RPMG25M + RPMG25G + RPMG26 : + RPMG26Mn + RPMG26Gn + Rp26M + Rp26G + RPAL27 + RPAL27a : + RPSI28 + RPSI29 + RPSI30 + RPN15 + RPO17 + RPNE20 IF (RATETOT.NE.0d0) THEN DO I = 1,50 FRAC(I,K) = 0d0 END DO FRAC(3 ,K) = - RPP + RPD FRAC(4 ,K) = - RPD FRAC(5 ,K) = RPLI FRAC(6 ,K) = RPBE FRAC(7 ,K) = RPB11 FRAC(8 ,K) = RPC13 - RPC FRAC(9 ,K) = RPC14 FRAC(10,K) = RPN15A - RPN - RPO18a + RPN15 FRAC(11,K) = - RPO + RPO17a + RPO17 FRAC(12,K) = RPO18 + RPO18a - RPO17 FRAC(13,K) = RPF19 + RPF19A - RPO18 FRAC(14,K) = RPNE21 - RPNE20 FRAC(15,K) = RPNE22 FRAC(16,K) = RPNA22 - RPNE21 FRAC(17,K) = RPNA23 + RPNA23n + RPNA23a - RPNE22 FRAC(18,K) = RPMG24 - RPNA23 - RPAL27a FRAC(19,K) = RPMG25M + RPMG25G FRAC(20,K) = RPMG26 + RPMG26Mn + RPMG26Gn FRAC(21,K) = Rp26M - RPMG25M - RPMG26Mn FRAC(22,K) = Rp26G - RPMG25G - RPMG26Gn FRAC(23,K) = RPAL27 + RPAL27a - RPMG26 FRAC(24,K) = RPSI28 - RPAL27 FRAC(25,K) = RPSI29 - RPSI28 FRAC(26,K) = RPSI30 - RPSI29 FRAC(27,K) = - RPSI30 C Why didn't I do the major isotopes as well? DO I=1,50 FRAC(I,K) = FRAC(I,K)/RATETOT DHNUC(I+50*(ISTAR-1),K) = DHNUC(I+50*(ISTAR-1),K) : - FRAC(I,K)*BARYN(I)*PROTONS(K) END DO C DHNUC(41+50*(ISTAR-1),K) = 0d0 END IF END DO RETURN END