From 91dd6f96a1ad48c9402e3214e8b16e7eca19ac8c Mon Sep 17 00:00:00 2001 From: salman Date: Thu, 18 Jun 2026 10:08:13 +0200 Subject: [PATCH] Split the file into sections with clear banners --- fortran/trimspNL.F | 131 ++++++++++++++++++++++++++++++++++----------- 1 file changed, 99 insertions(+), 32 deletions(-) diff --git a/fortran/trimspNL.F b/fortran/trimspNL.F index ec456ce..5e32a7e 100644 --- a/fortran/trimspNL.F +++ b/fortran/trimspNL.F @@ -45,6 +45,9 @@ c #endif IMPLICIT NONE +C======================================================================= +C DECLARATIONS AND GLOBAL STORAGE +C======================================================================= CHARACTER*16 TRIMSP_VERSION C These parameters are related to the maximum number of layers MAXNL C and maximum number of points in the depth distribution MAXD @@ -70,6 +73,18 @@ C Maximum number of elements in each layer, was limited to 5. PARAMETER (MAXNL5p2=MAXNL5*MAXNL5*MAXD) PARAMETER (MAXNLm15=(MAXNL-1)*MAXEL) PARAMETER (TRIMSP_VERSION='1.3.2') +C +C Named integer options used by the old input format. The input +C file still contains the numeric values; these PARAMETER names make +C comparisons and output descriptions easier to read. +C + INTEGER POT_KRC,POT_MOLIERE,POT_ZBL + INTEGER STOP_LS,STOP_OR,STOP_MIXED,STOP_ICRU49,STOP_ZIEGLER + INTEGER IRL_OFF + PARAMETER (POT_KRC=1,POT_MOLIERE=2,POT_ZBL=3) + PARAMETER (STOP_LS=1,STOP_OR=2,STOP_MIXED=3) + PARAMETER (STOP_ICRU49=4,STOP_ZIEGLER=5) + PARAMETER (IRL_OFF=0) LOGICAL TEST(64),TESTR(2000),TEST1(2000) LOGICAL EQUAL INTEGER*4 ISRCHFGT,ISRCHFGE,ILLZ @@ -295,6 +310,9 @@ C CHARACTER Variables COMMON /A/ M1,VELC,ZARG COMMON /B/ TI,SHEATH,CALFA +C======================================================================= +C STATIC INITIALISATION +C======================================================================= DATA PI/3.14159265358979D0/, ICW/100/, E2/14.399651D0/ DATA AB/0.52917725D0/, FP/0.885341377D0/, AN/0.60221367D0/ DATA inext/'.inp'/,outext/'.out'/,rgeext/'.rge'/ @@ -507,6 +525,9 @@ C RGENAM range output file name, FILEIN//'.rge' C ERRNAM error/log file name, FILEIN//'.err' C======================================================================= +C======================================================================= +C COMMAND LINE AND FILE NAMES +C======================================================================= CALL getarg(1, filein) C Require an explicit run basename for input/output naming. if (filein.eq.'') then @@ -527,6 +548,9 @@ C LMAX is maximum number of layers and JMAX is maximum number of C elements per layer. JMAX=MAXEL +C======================================================================= +C READ GUI-GENERATED INPUT FILE +C======================================================================= C This part reads the input file (new format). C The JavaScript/Electron frontend writes this sequential block layout, C so the READ order below must stay in sync with CreateInpFile(). @@ -589,6 +613,9 @@ C value A-5 of the ziegler tables 1359 CLOSE(UNIT=11) +C======================================================================= +C OPEN OUTPUT FILES AND START CLOCK +C======================================================================= C open statement for output files, removed from line 2449 ff to here OPEN(UNIT=21,FILE=outnam) 6001 OPEN(UNIT=22,FILE=rgenam,STATUS='replace') @@ -608,6 +635,9 @@ C Get simulation start time 1050 FORMAT(1x,' Start: ',A2,'.',A4,1x,A4,1x,A2 & ,':',A2,':',A2) +C======================================================================= +C OUTPUT BINNING CONSTANTS +C======================================================================= C SET INTERVAL CONSTANTS FOR OUTPUT DE = 1.D0 DA = 3.D0 @@ -631,6 +661,9 @@ C SET INTERVAL CONSTANTS FOR OUTPUT DGW = BW/DG DGIK = BW/DGI +C======================================================================= +C TARGET SETUP AND PRECOMPUTED TABLES +C======================================================================= C CALCULATION OF CHARGE AND MASS DEPENDENT CONSTANTS PI2=2.D0*PI ABC=AB*FP @@ -704,22 +737,22 @@ C For each layer calculate the following EC1(J) = 4.D0*MU1(J)/((1.D0+MU1(J))*(1.D0+MU1(J))) C KR-C (IPOT=1), MOLIERE (IPOT=2), ZBL POTENTIAL (IPOT=3) A1(J) = CVMGT(CA*ABC*(ZZ(J)**(-1.D0/3.D0)),CA*ABC/(Z1**0.23D0 - & +ZZ(J)**0.23D0),IPOT.LT.3) + & +ZZ(J)**0.23D0),IPOT.LT.POT_ZBL) F1(J) = A1(J)*TM(J)/(Z1*ZZ(J)*E2*(M1+TM(J))) KL1(J) = 1.212D0*Z1**(7.D0/6.D0)*ZZ(J)/ ((Z1**(2.D0/3.D0)+ZZ(J) & **(2.D0/3.D0))**1.5D0*DSQRT(M1)) ENDDO - IF(IPOT.EQ.1) THEN + IF(IPOT.EQ.POT_KRC) THEN C KR-C POTENTIAL (IPOT=1) DO J=1,LJ KOR1(J) = 0.0389205D0*KL1(J)/(PI*A1(J)*A1(J)) ENDDO - ELSEIF (IPOT.EQ.2) THEN + ELSEIF (IPOT.EQ.POT_MOLIERE) THEN C MOLIERE POTENTIAL (IPOT=2) DO J=1,LJ KOR1(J) = 0.045D0*KL1(J)/(PI*A1(J)*A1(J)) ENDDO - ELSEIF (IPOT.EQ.3) THEN + ELSEIF (IPOT.EQ.POT_ZBL) THEN C ZBL POTENTIAL DO J=1,LJ KOR1(J) = 0.0203253D0*KL1(J)/(PI*A1(J)*A1(J)) @@ -731,28 +764,29 @@ C ZBL POTENTIAL EC(I,J) = 4.D0*MU(I,J)/((1.D0+MU(I,J))*(1.D0+MU(I,J))) C KR-C , MOLIERE , ZBL POTENTIAL A(I,J)= CVMGT(CA*ABC/(DSQRT(ZZ(I))+DSQRT(ZZ(J)))**(2.D0/3.D0 - & ),CA*ABC/(ZZ(I)**0.23D0+ZZ(J)**0.23D0),IPOTR.LT.3) + & ),CA*ABC/(ZZ(I)**0.23D0+ZZ(J)**0.23D0), + & IPOTR.LT.POT_ZBL) C ZBL POTENTIAL F(I,J) = A(I,J)*TM(J)/(ZZ(I)*ZZ(J)*E2*(TM(I)+TM(J))) KL(I,J) = 1.212D0*ZZ(I)**(7.D0/6.D0)*ZZ(J)/ ((ZZ(I)**(2.D0 & /3.D0)+ZZ(J)**(2.D0/3.D0))**1.5D0*DSQRT(TM(I))) ENDDO ENDDO - IF (IPOTR.EQ.1) THEN + IF (IPOTR.EQ.POT_KRC) THEN C KR-C POTENTIAL (IPOTR=1) DO I = 1,LJ DO J = 1,LJ KOR(I,J) = 0.0389205D0*KL(I,J)/(PI*A(I,J)*A(I,J)) ENDDO ENDDO - ELSEIF (IPOTR.EQ.2) THEN + ELSEIF (IPOTR.EQ.POT_MOLIERE) THEN C MOLIERE POTENTIAL (IPOTR=2) DO I = 1,LJ DO J = 1,LJ KOR(I,J) = 0.045D0*KL(I,J)/(PI*A(I,J)*A(I,J)) ENDDO ENDDO - ELSEIF (IPOTR.EQ.3) THEN + ELSEIF (IPOTR.EQ.POT_ZBL) THEN C ZBL POTENTIAL (IPOTR=3) DO I = 1,LJ DO J = 1,LJ @@ -806,6 +840,9 @@ C SET CONSTANT DISTANCES IF(E0.GE.0.D0) GO TO 51 C +C======================================================================= +C INITIAL PROJECTILE BATCH +C======================================================================= C SET CONSTANTS FOR MAXWELLIAN DISTRIBUTION C TI = -1.D0*E0 @@ -957,6 +994,9 @@ C LLL(IV) = JL ENDDO C +C======================================================================= +C PRIMARY PROJECTILE TRANSPORT LOOP +C======================================================================= C PROJECTILE LOOP C Each pass transports one projectile history. The state vectors C X/Y/Z, E, direction cosines, current layer index, etc. are updated @@ -975,6 +1015,9 @@ C ENDDO KK1=KK0 C +C======================================================================= +C PROJECTILE COLLISION CALCULATION +C======================================================================= C COLLISION LOOP (INCLUDES WEAK SIMULTANEOUS COLL. FOR KK1.LT.4) C KK controls the treatment of weak simultaneous collisions. C For every active projectile we choose a collision partner, sample an @@ -1013,7 +1056,7 @@ C C C MAGIC (DETERMINATION OF SCATTERING ANGLE : KRYPTON-CARBON POT.) C - IF(IPOT.NE.1) GO TO 4101 + IF(IPOT.NE.POT_KRC) GO TO 4101 C KRYPTON-CARBON POTENTIAL 104 DO IV=IVMIN,IVMAX @@ -1056,7 +1099,7 @@ C GET MAX AND MIN INDEX OF TEST FAILURES ENDDO GO TO 4103 - 4101 IF(IPOT.NE.2) GO TO 4102 + 4101 IF(IPOT.NE.POT_MOLIERE) GO TO 4102 C MOLIERE POTENTIAL C CALL MAGICMOL(C2(1),S2(1),B(1),R(1),EPS(1),IH1) 4104 DO IV=IVMIN,IVMAX @@ -1096,7 +1139,7 @@ C GET MAX AND MIN INDEX OF TEST FAILURES ENDDO GO TO 4103 - 4102 IF(IPOT.NE.3) GO TO 4103 + 4102 IF(IPOT.NE.POT_ZBL) GO TO 4103 C ZBL POTENTIAL C CALL MAGICZBL(C2(1),S2(1),B(1),R(1),EPS(1),IH1) 5104 DO IV=IVMIN,IVMAX @@ -1163,7 +1206,7 @@ C CPSI(IV)=CVMGT(TA2,-TA2,CU.GT.0.D0) SPSI(IV)=DABS(TA)*TA2 DEEOR=CVMGT(KOR1(JJJ(IV))*DSQRT(DABS(E(IV)))*EX1(IV),0.D0, - & KDEE1.EQ.2.OR.KDEE1.EQ.3) + & KDEE1.EQ.STOP_OR.OR.KDEE1.EQ.STOP_MIXED) DENS(IV)=DENS(IV)+DEN(IV) DEES(IV)=DEES(IV)+DEEOR ENDDO @@ -1176,12 +1219,17 @@ C C C END OF COLLISION LOOP C +C======================================================================= +C PROJECTILE ENERGY LOSS AND DAMAGE ACCUMULATION +C======================================================================= C INELASTIC ENERGY LOSS( 5 POSSIBILITIES) C DO IV=1,IH1 ASIGT(IV)=(LM(LLL(IV))-TAU(IV)+TAUPSI(IV))*ARHO(LLL(IV)) TAUPSI(IV)=TAU(IV)*DABS(CPSI(IV)) ENDDO +C KDEE1 uses STOP_LS, STOP_OR, STOP_MIXED, STOP_ICRU49, +C and STOP_ZIEGLER in this order. GO TO(15,16,17,18,19),KDEE1 15 DO IV=1,IH1 DEE(IV)=CVMGT(0.D0,KLM1(LLL(IV))*ASIGT(IV)*DSQRT(E(IV)),X(IV) @@ -1321,8 +1369,11 @@ C ENDDO 89 CONTINUE - IF(IRL.EQ.0) GO TO 27 + IF(IRL.EQ.IRL_OFF) GO TO 27 C +C======================================================================= +C RECOIL GENERATION AND TRANSPORT +C======================================================================= C VECTORIZED RECOIL LOOP C C TARGET RECOIL ATOM SECTION @@ -1450,7 +1501,7 @@ C C C MAGIC (DETERMINATION OF SCATTERING ANGLE : KRYPTON-CARBON POT.) C - IF(IPOTR.NE.1) GO TO 4201 + IF(IPOTR.NE.POT_KRC) GO TO 4201 C KR-C POTENTIAL C CALL MAGICKRC(C2R(1),S2R(1),BR(1),RR(1),EPSR(1),NREC2) 205 DO IV=IVMIN,IVMAX @@ -1494,7 +1545,7 @@ C GET MAX AND MIN INDEX OF TEST FAILURES ENDDO GO TO 4203 - 4201 IF(IPOTR.NE.2) GO TO 4202 + 4201 IF(IPOTR.NE.POT_MOLIERE) GO TO 4202 C MOLIERE POTENTIAL C CALL MAGICMOL(C2R(1),S2R(1),BR(1),RR(1),EPSR(1),NREC2) 4205 DO IV=IVMIN,IVMAX @@ -1537,7 +1588,7 @@ C GET MAX AND MIN INDEX OF TEST FAILURES ENDDO GO TO 4203 - 4202 IF(IPOTR.NE.3) GO TO 4203 + 4202 IF(IPOTR.NE.POT_ZBL) GO TO 4203 C ZBL POTENTIAL C CALL MAGICZBL(C2R(1),S2R(1),BR(1),RR(1),EPSR(1),NREC2) 5205 DO 5206 IV=IVMIN,IVMAX @@ -1589,7 +1640,7 @@ C T1=CVMGT(T(IREC1),0.D0,KKR.EQ.3) TR1=TR1+T1 DEEORR=CVMGT(0.D0,KOR(JJR(IREC1,2),JJR(IREC1,1)) - & * DSQRT(ER(IREC1,2))*EX1R(IREC1),KDEE2.EQ.1) + & * DSQRT(ER(IREC1,2))*EX1R(IREC1),KDEE2.EQ.STOP_LS) DEERS(IREC1)=DEERS(IREC1)+DEEORR TAUR(IREC1)=CVMGT(PR(IREC1)*DSQRT(S2R(IREC1)/C2R(IREC1)),0 & .D0,KKR.EQ.0) @@ -1645,6 +1696,7 @@ C & *ARHO(LRR(IREC1,2)) TAUPSR(IREC1,2)=TAUR(IREC1)*DABS(CPSIR(IREC1,2)) ENDDO +C KDEE2 uses STOP_LS, STOP_OR, and STOP_MIXED in this order. GO TO(115,116,117),KDEE2 115 DO IREC1=1,NREC2 DEER(IREC1)=CVMGT(0.D0,KLM(LRR(IREC1,2), JJR(IREC1,2)) @@ -1959,6 +2011,9 @@ C 27 CONTINUE IF(IH1.EQ.0.AND.IH.EQ.NH) GO TO 140 C +C======================================================================= +C PROJECTILE STOP, BACKSCATTER AND TRANSMISSION TESTS +C======================================================================= C PROJECTILE CANDIDATE FOR REFLECTION C DO IV=1,IH1 @@ -2408,6 +2463,9 @@ C advanced by the loop above, so skip the trailing range. 140 IF(NREC1+NREC2.GT.0) GO TO 83 C C +C======================================================================= +C WRITE MAIN OUTPUT REPORT +C======================================================================= C PRINTOUT C C @@ -2506,7 +2564,7 @@ C ENDIF 1418 FORMAT(/1X,I3,6H.LAYER,17X,5F6.2,3X,5F7.2,3X,5F6.2) 1416 CONTINUE - IF(KDEE1.LT.4) GO TO 1421 + IF(KDEE1.LT.STOP_ICRU49) GO TO 1421 WRITE(21,1419) 1419 FORMAT(//30X,'CH1',10X,'CH2',10X,'CH3',10X,'CH4',10X,'CH5') DO 1417 I=1,L @@ -2523,23 +2581,23 @@ C 1415 FORMAT(/1X,I3,6H.LAYER,17X,5F13.6) 1423 FORMAT(/25X,5F13.6) 1421 CONTINUE - IF(IPOT.EQ.1) DPOT='KR-C POTENTIAL' - IF(IPOT.EQ.2) DPOT='mod. MOLIERE ' - IF(IPOT.EQ.3) DPOT='ZBL POTENTIAL' - IF(IPOTR.EQ.1) DPOTR='KR-C POTENTIAL' - IF(IPOTR.EQ.2) DPOTR='MOLIERE POTENTIAL' - IF(IPOTR.EQ.3) DPOTR='ZBL POTENTIAL' + IF(IPOT.EQ.POT_KRC) DPOT='KR-C POTENTIAL' + IF(IPOT.EQ.POT_MOLIERE) DPOT='mod. MOLIERE ' + IF(IPOT.EQ.POT_ZBL) DPOT='ZBL POTENTIAL' + IF(IPOTR.EQ.POT_KRC) DPOTR='KR-C POTENTIAL' + IF(IPOTR.EQ.POT_MOLIERE) DPOTR='MOLIERE POTENTIAL' + IF(IPOTR.EQ.POT_ZBL) DPOTR='ZBL POTENTIAL' WRITE(21,1411) DPOT,DPOTR 1411 FORMAT(//7X,'INTERACTION POTENTIAL : PROJECTILE-TARGET : ',A18 & ,' TARGET-TARGET : ',A18) - IF(KDEE1.EQ.1) DKDEE1='LINDHARD-SCHARFF' - IF(KDEE1.EQ.2) DKDEE1='OEN-ROBINSON' - IF(KDEE1.EQ.3) DKDEE1='50% LS 50% OR' - IF(KDEE1.EQ.4) DKDEE1='AZ nach ICRU49' - IF(KDEE1.EQ.5) DKDEE1='ZIEGLER' - IF(KDEE2.EQ.1) DKDEE2='LINDHARD-SCHARFF' - IF(KDEE2.EQ.2) DKDEE2='OEN-ROBINSON' - IF(KDEE2.EQ.3) DKDEE2='50% LS 50% OR' + IF(KDEE1.EQ.STOP_LS) DKDEE1='LINDHARD-SCHARFF' + IF(KDEE1.EQ.STOP_OR) DKDEE1='OEN-ROBINSON' + IF(KDEE1.EQ.STOP_MIXED) DKDEE1='50% LS 50% OR' + IF(KDEE1.EQ.STOP_ICRU49) DKDEE1='AZ nach ICRU49' + IF(KDEE1.EQ.STOP_ZIEGLER) DKDEE1='ZIEGLER' + IF(KDEE2.EQ.STOP_LS) DKDEE2='LINDHARD-SCHARFF' + IF(KDEE2.EQ.STOP_OR) DKDEE2='OEN-ROBINSON' + IF(KDEE2.EQ.STOP_MIXED) DKDEE2='50% LS 50% OR' WRITE(21,1413) DKDEE1,DKDEE2 1413 FORMAT(//7X,'INELASTIC LOSS MODEL : PROJECTILE-TARGET : ',A18 & ,' TARGET-TARGET : ',A18) @@ -3569,6 +3627,9 @@ C & I0,3x)) C The *.seq file is the compact per-run summary consumed by the GUI for C scan plots (fractions stopped in each layer, backscattering, C transmission, mean depth, etc.). Each run writes its own summary file. +C======================================================================= +C WRITE SEQUENCE SUMMARY FILE +C======================================================================= OPEN(UNIT=33,FILE=seqnam,STATUS='replace') WRITE(33,7802) (chem(k),k=1,NLayers) WRITE(33,7801)E0keV,EsigkeV,ALPHA,ALPHASIG,NH,IIM,IB,IT,tryE @@ -4027,6 +4088,9 @@ C +C======================================================================= +C UTILITY SUBROUTINES AND FUNCTIONS +C======================================================================= SUBROUTINE MOMENTS(FIM0,SEM,THM,FOM,FIM,SIM,SIGMA,DFIM0,DSEM,DTHM, # X1S,X2S,X3S,X4S,X5S,X6S,Y) IMPLICIT NONE @@ -4544,6 +4608,9 @@ C in seconds from beginning of year END +C======================================================================= +C RANDOM NUMBER GENERATOR +C======================================================================= SUBROUTINE RANLUX(RVEC,LENV) C Subtract-and-borrow random number generator proposed by C Marsaglia and Zaman, implemented by F. James with the name