Suspected bug fixes.

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