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

Optimize fortran code and make it more readable. Omit multiple loops on the layer and fold everything in one loop.

Browse Source
This commit is contained in:
salman
2023-01-25 10:35:21 +01:00
parent 3c7d8a0e65
commit e5a0a8048a
2 changed files with 50 additions and 56 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
fortran/Makefile
View File

@@ -14,7 +14,8 @@ WARN = # -Wall -Wextra
DIALECT = -std=gnu
prefix = /usr/local
# OPS = -c $(DIALECT) $(WARN) $(DEBUG)
FCFLAGS = $(DIALECT) $(WARN) $(DEBUG) -O3 -mcmodel=large
#FCFLAGS = $(DIALECT) $(WARN) $(DEBUG) -O3 -mcmodel=large
FCFLAGS = $(DIALECT) $(WARN) $(DEBUG) -O3 -mcmodel=medium
all : trimspNL

103
fortran/trimspNL.F
View File

@@ -58,7 +58,7 @@ C accordingly.
PARAMETER (MAXD=1000)
PARAMETER (MAXNL=100)
C Maximum number of elements in each layer, was limited to 5.
PARAMETER (MAXEL=20)
PARAMETER (MAXEL=5)
PARAMETER (MAXD1=MAXD+1)
PARAMETER (MAXD2=MAXD+2)
PARAMETER (MAXD5=MAXD*MAXEL)
@@ -72,7 +72,7 @@ C Maximum number of elements in each layer, was limited to 5.
PARAMETER (MAXNLm15=(MAXNL-1)*MAXEL)
LOGICAL TEST(64),TESTR(2000),TEST1(2000)
LOGICAL EQUAL,FORT33
INTEGER*4 ISRCHEQ,ISRCHFGT,ISRCHFGE,ILLZ
INTEGER*4 ISRCHFGT,ISRCHFGE,ILLZ
INTEGER N,L,LL,NH,NUM,KK
INTEGER I,J,IV
INTEGER tryE,negE
@@ -283,8 +283,9 @@ c fuer part. reflec. coeff. Berechnung
REAL*4 random(1),ran2(2)
C CHARACTER Variables
CHARACTER*18 DPOT,DPOTR,DKDEE1,DKDEE2
CHARACTER filein*8,inext*4,fileout*8,outext*4,innam*12,outnam*12
CHARACTER rgenam*12,rgeext*4,errnam*12,errext*4
CHARACTER*60 filein,fileout
CHARACTER*64 innam,outnam,rgenam,errnam
CHARACTER*4 inext,outext,rgeext,errext
CHARACTER errcom*72
CHARACTER month_start*4,month_stop*4,day_start*2,day_stop*2
CHARACTER year_start*4,year_stop*4,hour_start*2,hour_stop*2
@@ -298,7 +299,6 @@ C CHARACTER Variables
DATA AB/0.52917725D0/, FP/0.885341377D0/, AN/0.60221367D0/
DATA inext/'.inp'/,outext/'.out'/,rgeext/'.rge'/
DATA errext/'.err'/
DATA filein/'eingabe1'/,fileout/'ausgabe1'/
DATA ET/0.D0/,PLST/0.D0/,PL2ST/0.D0/,PL3ST/0.D0/
DATA PL4ST/0.D0/,PL5ST/0.D0/,PL6ST/0.D0/
@@ -381,10 +381,21 @@ C part. refl. coeff. from Thomas et al.
DATA PRC/0.825D0,21.41D0,8.606D0,0.6425D0,1.907D0,1.927D0/
DATA number_in_layer/MAXNL*0/
innam=filein//inext
outnam=fileout//outext
rgenam=fileout//rgeext
errnam=fileout//errext
CALL getarg(1, filein)
C if there is no input argument take defailt file names
if (filein.eq.'') then
filein = 'eingabe1'
fileout = 'ausgabe1'
else
fileout = filein
endif
innam=trim(filein)//inext
outnam=trim(fileout)//outext
rgenam=trim(fileout)//rgeext
errnam=trim(fileout)//errext
write (*,*) innam
C LMAX is maximum number of layers and JMAX is maximum number of
C elements per layer.
@@ -450,14 +461,9 @@ C READ(11,*) CH5(I,1),CH5(I,2),CH5(I,3),CH5(I,4),CH5(I,5)
C open statement for output files, removed from line 2449 ff to here
OPEN(UNIT=21,FILE=outnam)
GOTO 6001
READ(*,'(A8)') fileout
outnam=fileout//outext
rgenam=fileout//rgeext
OPEN(UNIT=21,FILE=outnam,STATUS='replace')
6001 OPEN(UNIT=22,FILE=rgenam,STATUS='replace')
WRITE(21,1000)
1000 FORMAT(1H1/,6X,'* TrimSPNL 02.Apr.13 *')
1000 FORMAT(1H1/,6X,'* TrimSPNL v1.1.0 *')
C Get simulation start time
CALL TimeStamp(day_start,month_start,year_start,hour_start
@@ -465,7 +471,7 @@ C Get simulation start time
WRITE(21,*)
WRITE(21,1050)day_start,month_start,year_start,hour_start
& ,min_start,sec_start
WRITE(*,*) 'Version 1.0.1'
WRITE(*,*) 'Version 1.1.0'
WRITE(*,1050)day_start,month_start,year_start,hour_start
& ,min_start,sec_start
1050 FORMAT(1x,' Start: ',A2,'.',A4,1x,A4,1x,A2
@@ -499,62 +505,49 @@ C CALCULATION OF CHARGE AND MASS DEPENDENT CONSTANTS
L=NLayers
C Checks wether depth interval is an integer denominator of layer thickness or not
C If not, calculated implantation profile is not correct.
C Checks wether depth interval is an integer denominator of layer thickness or not.
C If not, the calculated number of implanted particles in each layer from bin width
C is wrong. Calculation form layer thickness should be correct.
C The rge files (stopping profiles) and summary files should be fine.
C CW is depth increment and DX(K) is the thickness of layer K
depth_interval_flag = 1
DO K=1,L-1
write(*,*) DX(K),CW
IF(.NOT.EQUAL(DX(K)/CW-DBLE(IDINT(DX(K)/CW)),0.D0)) THEN
C If the thickness of layer K is not an integer multiples of depth increment
depth_interval_flag = 0
write (*,*) 'I am here'
GO TO 44
ENDIF
ENDDO
44 CONTINUE
DO I=1,L
C This loop finds out how many elements are in the layer.
C I can be removed in the number is know from the input file.
C DO J=1,JMAX
C IF(EQUAL(CO(I,J),0.D0)) GOTO 156
C ENDDO
C J=JMAX+1
C I am guessing NJ(I) is the number of elements in layer I
C 156 NJ(I)=J-1
ENDDO
JT(1) = 0
JT(2) = NJ(1)
LJ= NJ(1)+NJ(2)
do I=3,L
JT(I)=JT(I-1)+NJ(I-1)
LJ=LJ+NJ(I)
enddo
XX(1)=DX(1)
DO I=2,L
XX(I)=XX(I-1)+DX(I)
ENDDO
C Loop over all defined layers
NJJ=0
DO I=1,L
if (I.ge.2) XX(I)=XX(I-1)+DX(I)
if (I.ge.3) then
JT(I)=JT(I-1)+NJ(I-1)
LJ=LJ+NJ(I)
endif
DO J=1,NJ(I)
Z2(I)=Z2(I)+CO(I,J)*ZT(I,J)
M2(I)=M2(I)+CO(I,J)*MT(I,J)
ENDDO
ENDDO
DO LL=1,L
ARHO(LL) = RHO(LL)*AN/M2(LL)
LM(LL) = ARHO(LL)**(-1.D0/3.D0)
ASIG(LL) = LM(LL)*ARHO(LL)
PDMAX(LL) = LM(LL)/DSQRT(PI)
K2(LL) = .133743D0*Z2(LL)**(2.D0/3.D0)/DSQRT(M2(LL))
AM(LL) = CA*ABC*(Z2(LL)**(-1.D0/3.D0))
FM(LL) = AM(LL)*M2(LL)/(Z1*Z2(LL)*E2*(M1+M2(LL)))
EPS0(LL) = FM(LL)*E0
ENDDO
ARHO(I) = RHO(I)*AN/M2(I)
LM(I) = ARHO(I)**(-1.D0/3.D0)
ASIG(I) = LM(I)*ARHO(I)
PDMAX(I) = LM(I)/DSQRT(PI)
K2(I) = .133743D0*Z2(I)**(2.D0/3.D0)/DSQRT(M2(I))
AM(I) = CA*ABC*(Z2(I)**(-1.D0/3.D0))
FM(I) = AM(I)*M2(I)/(Z1*Z2(I)*E2*(M1+M2(I)))
EPS0(I) = FM(I)*E0
C Loop over all defined layers
NJJ=0
do I=1,L
C For each layer calculate the following
do J=1,NJ(I)
ZZ(NJJ+J) = ZT(I,J)
@@ -563,14 +556,14 @@ C For each layer calculate the following
EP(NJJ+J) = BE(I,J)
enddo
NJJ=NJJ+NJ(I)
enddo
DO I=1,L
COM(1,I) = CO(I,1)
DO J=1,NJ(I)-1
COM(J+1,I) = COM(J,I)+CO(I,J+1)
ENDDO
ENDDO
DO J = 1,LJ
MU1(J) = M1/TM(J)
EC1(J) = 4.D0*MU1(J)/((1.D0+MU1(J))*(1.D0+MU1(J)))
@@ -632,6 +625,7 @@ C ZBL POTENTIAL (IPOTR=3)
ENDDO
ENDDO
ENDIF
DO LL=1,L
DO J=1,NJ(LL)
KLM1(LL) = KLM1(LL)+CO(LL,J)*KL1(J+JT(LL))*CK(LL)
@@ -1859,7 +1853,7 @@ C
I1 = MAX0( MIN0( IDINT(X(IV)/CW+1.D0), MAXD1), 0)
IRP(I1)=IRP(I1)+1
C
C Berechnung der gestoppten Teilchen im jeweiligen Layer
C Calculated the stopped particles in each layer
C
LowTiefe = 0.D0
UpTiefe = DX(1)
@@ -4220,8 +4214,7 @@ C why is it needed??
INTEGER FUNCTION ISRCHFGE(N,ARRAY,INC,TARGT)
IMPLICIT NONE
INTEGER I,N,J,INC
REAL*8 ARRAY(N)
REAL*8 TARGT
REAL*8 ARRAY(N),TARGT
J=1
IF(INC.LT.0) J=N*(-INC)