diff --git a/fortran/trimspNL.F b/fortran/trimspNL.F index 6ec41bc..9d00cdc 100644 --- a/fortran/trimspNL.F +++ b/fortran/trimspNL.F @@ -334,6 +334,9 @@ C CHARACTER Variables DATA PL6SUM/0.D0/,PLSUM/0.D0/ DATA ENL2B/0.D0/,ENUCLB/0.D0/,EINELI/0.D0/,EIL2I/0.D0/ DATA ENUCLI/0.D0/,ENL2I/0.D0/ +C BUGFIX: transmitted projectile loss accumulators were read before +C guaranteed assignment when at least one projectile transmits. + DATA ENUCLT/0.D0/,ENL2T/0.D0/,EINELT/0.D0/,EIL2T/0.D0/ DATA PLSB/0.D0/,PL2SB/0.D0/,PL3SB/0.D0/,PL4SB/0.D0/ DATA PL5SB/0.D0/,PL6SB/0.D0/ DATA EELWC/0.D0/,EELWC2/0.D0/,EELWC3/0.D0/,EELWC4/0.D0/ @@ -359,6 +362,9 @@ C CHARACTER Variables DATA ELE/MAXDNL5*0.D0/,ELI/MAXDNL5*0.D0/ DATA ICDIRI/MAXNL5p2*0/ DATA ICSUM/0/,ICSUMS/0/,ICDI/0/,ISPA/0/,ISPAT/0/ +C BUGFIX: recoil counters are used in IRL output and incremented +C later; they must not start from compiler-dependent garbage. + DATA ICDIR/0/,ICSBR/0/,ICSUMR/0/ DATA Z2/MAXNL*0.D0/,M2/MAXNL*0.D0/ DATA KLM1/MAXNL*0.D0/,CHM1/MAXNL*0.D0/ DATA SB/MAXNL*0.D0/,KLM/MAXNLp25*0.D0/ @@ -520,7 +526,9 @@ C 44 CONTINUE NJJ=0 C Loop over all defined layers DO I=1,L - if(.not.EQUAL(DX(K)/CW-DBLE(IDINT(DX(K)/CW)),0.D0)) then +C BUGFIX: this loop iterates with I; using K here read an +C uninitialised/stale index and could test the wrong layer. + if(.not.EQUAL(DX(I)/CW-DBLE(IDINT(DX(I)/CW)),0.D0)) then depth_interval_flag = 0 endif if (I.ge.2) XX(I)=XX(I-1)+DX(I) @@ -1729,6 +1737,14 @@ C IAG=IDINT(EXIRT*20.D0+1.D0) KADST(IAG)=KADST(IAG)+1 KDSTJ(IAG,JJR(IREC1,2))=KDSTJ(IAG,JJR(IREC1,2))+1 +C BUGFIX: transmission-sputter matrices below use IG. +C Compute it here instead of reusing a stale value. + IG=2 + SXR(IREC1)=DMAX1(SXR(IREC1),1.0D-12) + U=CYR(IREC1)/SXR(IREC1) + IF(DABS(U).GT.1.D0) U = SIGN(1.D0,U) + ACS=DACOS(U) + IG=IDINT(DGW*ACS+2.D0) C C 4 GROUPS:ION IN , PKA ;ION IN , SKA ;ION OUT, PKA ;ION OUT, SKA C @@ -1820,7 +1836,12 @@ C TEST(IV)=E(IV).LE.EF.OR.X(IV).LT.-1.D0*SU.OR.X(IV).GT.SUT ENDDO IVMIN=1+ILLZ(IH1,TEST,1) - IF(IVMIN.GT.IH1) GO TO 90 + IF(IVMIN.GT.IH1) THEN +C BUGFIX: no projectile needs special handling; ensure IVMAX +C is defined before the shared label-90 bookkeeping. + IVMAX=0 + GO TO 90 + ENDIF IVMAX=IH1-ILLZ(IH1,TEST,-1) DO 120 IV=IVMIN,IVMAX 160 IF(IV.GT.IH1) GO TO 90 @@ -2012,7 +2033,9 @@ C ET5SUM=ET5SUM+ET3*ETQ ET6SUM=ET6SUM+ET3*ET3 IPT = MAX0( MIN0( IDINT(PL(IV)/CW+1.D0), 100), 1) - IPLT(IP)=IPLT(IP)+1 +C BUGFIX: IPT is computed for transmitted projectiles here. +C IP belongs to different logic and may be uninitialised/stale. + IPLT(IPT)=IPLT(IPT)+1 PLQT=PL(IV)*PL(IV) PL3T=PLQT*PL(IV) PLST=PLST+PL(IV) @@ -2069,9 +2092,11 @@ C MEPT(IE,IPT) = MEPT(IE,IPT)+1 MEPT(NE,IPT) = MEPT(NE,IPT)+1 MEPT(IE,102) = MEPT(IE,102)+1 - EMAT(IG,IAGB) = EMAT(IG,IAGB)+ES - EMAT(IG,22) = EMAT(IG,22)+ES - EMAT(NG,IAGB) = EMAT(NG,IAGB)+ES +C BUGFIX: this is the transmission branch; use EST, not ES +C from the backscattering branch. + EMAT(IG,IAGB) = EMAT(IG,IAGB)+EST + EMAT(IG,22) = EMAT(IG,22)+EST + EMAT(NG,IAGB) = EMAT(NG,IAGB)+EST GO TO 110 C C PROJECTILE IS REFLECTED BACK INTO THE TARGET BY THE SURF. BARRIER @@ -2231,6 +2256,10 @@ C ENDDO 134 CONTINUE C +C BUGFIX: if no projectile satisfied TEST, IVMIN is IH1+1 and +C IVMAX was not set in this block. All particles were already +C advanced by the loop above, so skip the trailing range. + IF(IVMIN.GT.IH1) GO TO 132 IF(IVMAX.LT.IVMIN) GO TO 132 DO IV=IVMAX+1,IH1 LLL(IV) = MIN0(ISRCHFGT(L,XX(1),1,X(IV)),L) @@ -2266,7 +2295,9 @@ C E0keV=E0/1.D3 EsigkeV=Esig/1.D3 - IF(ZT(1,2).LT.1.0D-3) THEN +C BUGFIX: single-element target should be decided from NJ(1). +C ZT(1,2) is not read for NJ(1)=1 and may be uninitialised. + IF(NJ(1).EQ.1) THEN epsilon = 32.55D0*(MT(1,1)/M1)/(1.D0+(MT(1,1)/M1))* 1.D0/(Z1 & *ZT(1,1)*DSQRT(Z1**(2.D0/3.D0)+ZT(1,1)**(2.D0/3.D0))) & * E0keV @@ -3272,8 +3303,10 @@ C KDSTL(I,1)=KDSTL(I,1)+KDSTJ(I,J) ENDDO ENDDO - DO J=NJ(1)+1,JT(3) - KDSTL(I,2)=KDSTL(I,2)+KDSTJ(I,J) + DO I=1,20 + DO J=NJ(1)+1,JT(3) + KDSTL(I,2)=KDSTL(I,2)+KDSTJ(I,J) + ENDDO ENDDO 1766 CONTINUE DO J=1,2 @@ -4013,6 +4046,10 @@ C IMPLICIT NONE INTEGER INIV1,INIV3 REAL*8 FG(128),FFG(64) +C BUGFIX: VELOC reuses buffered Gaussian deviates between calls. +C The buffer indices and arrays must be saved and initialized. + SAVE INIV1,INIV3,FG,FFG + DATA INIV1/0/,INIV3/0/ REAL*8 COSX,COSY,COSZ,SINE REAL*8 M1,VELC,ZARG REAL*8 VELX,VELY,VELZ,VELQ,VEL,E