**==SOLVER.FOR SUBROUTINE SOLVER (NOC, ITER, JTER, ERR, IE, NWRT5) IMPLICIT REAL*8(A-H, O-Z) SAVE PARAMETER (MAXMSH = 2000) COMMON H(60,MAXMSH),DH(60,MAXMSH),EPS,DEL,DH0,NMESH,JIN,IW(200) COMMON /NUCMAT/ HNUC(100,MAXMSH), DHNUC(100,MAXMSH) COMMON /SOLV / C(MAXMSH+1,51,51),S(50,151),ER(60),D,KH,NE,N2,N5, : N6, N10, N12, N14, MM, IA(1), NE1, NE2, NE3, NB, : NEV, NF, J1, J2, IH, JH, ID(90), K1, K2, NV, NW COMMON /DTCONT/ VLHP(2), VLEP(2), VLCP(2), RLFP(2), TOTMP(2), VLHC(2), : VLEC(2), VLCC(2), RLFC(2), TOTMC(2) COMMON /MIXFUD/ SGTHFAC, FACSGMIN, FACSG, ISGFAC DIMENSION MK(60), ERT(60), IE(100) LOGICAL REDUCE REDUCE = .FALSE. C EG-style mixing convergence trick IF (NOC.EQ.1.AND.ISGFAC.EQ.1) THEN FACSG = FACSGMIN ELSE FACSG = 1d0 END IF C 1/9/09 RJS -- the following section of code has been removed. It was originally there C to stop problems that occurred during matrix inversion due to the use of some elements C from previous iterations (i.e. the nucleosynthesis). This massively slowed the code in C non-nucleosynthesis mode, and no longer seems to be a problem. If you find you need it, C add it back in again :) C Clear C and S C DO J = 1, 51 C DO I = 1,51 C DO K = 1,NMESH+1 !MAXMSH+1 C C(K,I,J) = 0d0 C END DO C END DO C END DO C DO J = 1,151 C DO I = 1,50 C S(I,J) = 0d0 C END DO C END DO C Replace loop with written out stuff NE1 = IE(1) NE2 = IE(2) NE3 = IE(3) NB = IE(4) NEV = IE(5) NF = IE(6) J1 = IE(7) J2 = IE(8) IH = IE(9) JH = IE(10) DO I=1,90 !45 ID(I) = IE(10+I) END DO NE = NE1 + NE2 ERR = 0.0D0 IF ( NE.EQ.0 ) RETURN NV = NE + NEV NW = NV + NV N1 = NB + 1 N2 = NE + 1 N3 = NE + NB N4 = N3 + 1 N5 = NE + NE N6 = N5 + 1 N7 = N5 + NB N8 = N7 + 1 N9 = N5 + NE N10 = N9 + 1 N11 = N9 + NEV N12 = N10 + NEV N13 = NE - NB N15 = N13 + 1 N16 = N13 + NEV N14 = N16 + 1 NE4 = NE3 + NE2 C SURFACE IS AT K1'TH MESHPOINT, CENTRE AT K2'TH K1 = 1 + 0.01D0*J1*(NMESH-1) K2 = NMESH - 0.01D0*J2*(NMESH-1) KMESH = K2 - K1 + 1 C DETERMINE 'TYPICAL' VALUES FOR EACH VARIABLE IF (NOC.LE.1) THEN DO J = 1, NV L = ID(J) ER(L) = 1.0D0 IF (L.EQ.14.OR.L.EQ.29) ER(L) = DMAX1(1d-4,H(L,1)) DO K = K1, K2 ER(L) = DMAX1(ER(L), DABS(H(L,K))) END DO END DO ELSE ISTAR = NOC - 1 IF (ISTAR.EQ.1) write (32,*) "Star 1 nucleosynthesis" IF (ISTAR.EQ.2) write (32,*) "Star 2 nucleosynthesis" C Scale minor elements relative to most abundant DO J = 1,NV IF (J.LE.30) THEN L = ID(J) ELSE L = ID(30) + (J-30) END IF ER(L) = 1d-8 ER(1) = 1d0 ER(2) = 1d0 C ER(28) = 1d0 DO K = K1, K2 ER(L) = DMAX1(ER(L), DABS(HNUC(L+50*(ISTAR-1),K))) END DO END DO END IF MM = 2 - MIN0(1, NE4) C BEGIN ITERATIVE LOOP 600 CONTINUE DO KTER = 1, ITER D = 0.0D0 K = K1 KH = IH C EVALUATE FUNCTIONS AT SURFACE MESHPOINT CALL DIFRNS(K, NOC, ITER, N15, NE, 60+NEV) CALL DIVIDE(N15, N6, NE, N7, N8, K) K = K1 + 1 IF ( K.GT.K2 ) GO TO 100 C DITTO NEXT-TO-SURFACE, ELIMINATING SOME UNKNOWNS CALL DIFRNS(K, NOC, ITER, 1, NE, 30) CALL ELIMN8(N2, NE, N3, N15, N4, N5, K-1) CALL DIVIDE(1, N4, NE, N7, N8, K) 50 CONTINUE K = K + 1 IF ( K.LE.K2 ) THEN IF ( K.EQ.K1+3 .OR. K.EQ.JH-1 .OR. K.EQ.JH+2 .OR. : K.EQ.K2-1 ) KH = IH - KH C DITTO REMAINING MESHPOINTS CALL DIFRNS(K, NOC, ITER, 1, NE, 30) CALL ELIMN8(1, NE4, NB, N15, N1, NE, K-2) CALL ELIMN8(N1, NE4, NE, 1, N4, N5, K-1) CALL ELIMN8(N2, NE, N3, N15, N4, N5, K-1) CALL DIVIDE(1, N4, NE, N7, N8, K) GO TO 50 END IF CALL DIFRNS(K, NOC, ITER, 1, N16, 60) CALL ELIMN8(N2, N16, N3, N15, N4, N5, K-2) CALL ELIMN8(N4, N16, N5, 1, N8, N9, K-1) CALL ELIMN8(N6, N16, N7, N15, N8, N9, K-1) C SOLVE FOR CORRECTIONS AT CENTRAL MESHPOINT CALL DIVIDE(1, N8, N16, N11, N12, K) C BY BACK-SUBSTITUTION FIND CORRECTIONS THROUGHOUT NN = 1 60 CONTINUE K = K - 1 IF ( K.NE.K1 .OR. NB.NE.0 ) THEN IF ( K.EQ.K1 ) NN = N15 KK = K + 1 DO J = 1, N16 IF ( J.GE.N15 ) KK = K2 + 1 DO I = NN, NE C(K, I, N14) = C(K, I, N14) - C(K, I, J) : *C(KK, J, N14) END DO END DO IF ( K.GT.K1 ) GO TO 60 END IF DO K = K1, K2 C Put corrections into first column IF ( NEV.NE.0 ) THEN DO I = N2, NV C(K, I, 1) = C(K2+1, I-NB, N14) END DO END IF IF ( NB.NE.0 ) THEN DO I = 1, NB C(K, I, 1) = C(K, I+N13, N14) END DO END IF DO I = 1, N13 C(K, I+NB, 1) = C(K+1, I, N14) END DO END DO 100 CONTINUE IF ( IH.NE.0 ) THEN DO K = K1, K2 WRITE(24, 99001) (C(K,I,1), I=1, NV) END DO 99001 FORMAT (10E12.4) END IF ERR = 0.0D0 ERMAX = 0.0D0 DO J = 1, NV C ESTIMATE ACCURACY OF ITERATION VX = 0.0D0 DO K = K1, K2 VZ = C(K, J, 1) IF ( DABS(VZ).GT.DABS(VX) ) VX = VZ IF ( VX.EQ.VZ ) MK(J) = K ERR = ERR + DABS(VZ) END DO ERT(J) = VX IF (ABS(ERT(J)).GT.ERMAX) JMAX = J IF (ABS(ERT(J)).GT.ERMAX) KMAX = MK(J) ERMAX = MAX(ERMAX, ABS(ERT(J))) END DO ERR = ERR/(NV*KMESH) * * Extra code to reduce FAC when in difficulty. * IF (KTER .GT. 3 .AND. ERR .GE. 0.99D0*ERRPR) REDUCE = .TRUE. IF (REDUCE) THEN IF (ERRPR .GT. ERR) THEN FAC = 1.2D0*FAC IF (FAC .GE. 1D0) THEN FAC = 1D0 REDUCE = .FALSE. END IF ELSE FAC = 0.8D0*FAC END IF ELSE FAC = DEL/DMAX1(ERR, DEL) END IF FRR = DLOG(ERR)*0.43429D0 IF (NOC.GE.2.AND.FRR.GT.-5.85) FAC = 0.9d0 IF (NV.LE.12) THEN IF (KTER.GT.NWRT5 .OR. ERR.GE.1D-2 .OR. ERMAX.GE.1D1) : WRITE(32, 99002) KTER, FRR, FAC, (MK(J), ERT(J), J=1, NV) ELSE IF (KTER.GT.NWRT5 .OR. ERR.GE.1D-2 .OR. ERMAX.GE.1D1) : WRITE(32, 99003) KTER, FRR, FAC, (MK(J), ERT(J), J=1, NV) END IF 99002 FORMAT (1X, I2, 2F7.2, 12(I4,F7.4)) 99003 FORMAT (1X, I2, 2F7.2, 12(I4,F7.4),3(/,17X,12(I4,F7.4))) C RJS 10/11/02 - Added lines to print out reason for failure IF ((KTER.LT.6).OR.ITER.GT.30) THEN !NOC.GE.2.AND. IF (ERR.GE.1d-1) ERR = 1d-2 IF (ERMAX.GT.1d1) ERMAX = 1d0 END IF IF (ERR.GE.1d-1) THEN WRITE(32,*) "ERR tolerance exceeded:",ERR END IF IF (ERMAX.GE.1D1) THEN WRITE(32,*) "ERMAX tolerance exceeded in variable ", JMAX WRITE(32,*) "at meshpoint ", KMAX, " ERMAX=", ERMAX END IF C NaN trap IF (ERR.NE.ERR) THEN WRITE(32,*) "SNAFU, restarting" ERR = 1d1 END IF CALL FLUSH (32) IF ( IH.GT.0 ) IH = IH - 1 IE(9) = IH C APPLY CORRECTIONS, SCALED DOWN IF TOO LARGE IF (NOC.LE.1) THEN DO I = 1, NV J = ID(I) DO K = K1, K2 DH(J, K) = DH(J, K) - C(K, I, 1)*FAC*ER(J) END DO END DO ELSE DO I = 1, NV IF (I.LE.30) THEN J = ID(I) ELSE J = ID(30) + (I-30) END IF DO K = K1, K2 DHNUC(J+50*(ISTAR-1),K) = DHNUC(J+50*(ISTAR-1),K) : - C(K, I, 1)*FAC*ER(J) END DO END DO END IF * * Don't let the error get too large. * JTER = KTER IF (ERR.LT.EPS.AND.FACSG.NE.1d0) THEN FACSG = 10.0*FACSG FACSG = DMIN1(1d0,FACSG) write(32,*) "Boosting mixing", FACSG ERR = 1d-3 C Yes, I know I shouldn't do this, but haven't worked out another way yet GOTO 600 END IF C Sort out neutrons once convergence ok IF (ERR.LT.EPS.AND.NOC.GE.2) CALL NEUTRON(ISTAR) IF (ERR .LT. EPS .OR. ERR.GE.1d-1 .OR. ERMAX.GE.1D1) RETURN 200 ERRPR = ERR C CONTINUE ITERATING IF NOT YET ACCURATE ENOUGH END DO RETURN END