From e5a0a8048abab1a620496bc5773e6ce7bdfa8447 Mon Sep 17 00:00:00 2001 From: salman Date: Wed, 25 Jan 2023 10:35:21 +0100 Subject: [PATCH] Optimize fortran code and make it more readable. Omit multiple loops on the layer and fold everything in one loop. --- fortran/Makefile | 3 +- fortran/trimspNL.F | 103 +++++++++++++++++++++------------------------ 2 files changed, 50 insertions(+), 56 deletions(-) diff --git a/fortran/Makefile b/fortran/Makefile index c967ccd..5667dee 100644 --- a/fortran/Makefile +++ b/fortran/Makefile @@ -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 diff --git a/fortran/trimspNL.F b/fortran/trimspNL.F index 50f5b48..f0b145f 100644 --- a/fortran/trimspNL.F +++ b/fortran/trimspNL.F @@ -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)