LMU/TRIMSP
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Split the file into sections with clear banners

Browse Source
This commit is contained in:
salman
2026-06-18 10:08:13 +02:00
parent d2d9f786e9
commit 91dd6f96a1
1 changed files with 99 additions and 32 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fortran/trimspNL.F
+99 -32
View File
@@ -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