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

More fortran code optimizations and cleanup.

Browse Source
This commit is contained in:
salman
2023-01-25 14:38:51 +01:00
parent 6cc2ccf353
commit 687537eb45
2 changed files with 40 additions and 55 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fortran/trimspNL.F
+40 -55
View File
@@ -206,7 +206,7 @@ C CH1,2,3,4,5 are values of A-1,2,3,4,5 of the Ziegler tables
REAL*8 E2,AB,FP,AN
REAL*8 Esig,Epar
REAL*8 E0keV,EsigkeV
c fuer part. reflec. coeff. Berechnung
c for part. reflec. coeff. calculation
REAL*8 epsilon, prcoeff,PRC(6)
REAL*8 cossin,rphi,vanlb,vailb,vanlt,vailt,phip,ta,ta2
REAL*8 exir,phipr,u,rq,es,vanli,vaili
@@ -490,11 +490,12 @@ C SET INTERVAL CONSTANTS FOR OUTPUT
NE1 = NE -1
NA1 = NA -1
NG1 = NG -1
IF(E0.LT.0.) GO TO 2
E0DE = 100.0D0/(E0*DE)
GO TO 3
2 E0DE = 10.0D0/(DABS(E0)*DE)
3 BW = 180.D0/PI
IF(E0.LT.0.) then
E0DE = 10.0D0/(DABS(E0)*DE)
else
E0DE = 100.0D0/(E0*DE)
endif
BW = 180.D0/PI
DAW = BW/DA
DGW = BW/DG
DGIK = BW/DGI
@@ -511,25 +512,26 @@ 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 DO K=1,L-1
C 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
C depth_interval_flag = 0
C GO TO 44
C ENDIF
C ENDDO
C 44 CONTINUE
JT(1) = 0
JT(2) = NJ(1)
LJ= NJ(1)+NJ(2)
XX(1)=DX(1)
C Loop over all defined layers
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
depth_interval_flag = 0
endif
if (I.ge.2) XX(I)=XX(I-1)+DX(I)
if (I.ge.3) then
JT(I)=JT(I-1)+NJ(I-1)
@@ -668,6 +670,7 @@ C SET CONSTANT DISTANCES
XC = CVMGT( X0, -SU, X0.GT.0.D0)
RT = TT-RD
IF(E0.GE.0.D0) GO TO 51
C
C SET CONSTANTS FOR MAXWELLIAN DISTRIBUTION
@@ -693,21 +696,16 @@ C
IF ( E0.GT.0.D0 ) GO TO 47
IF ( ALPHA.GE.0.D0 ) THEN
C
C MAXWELLIAN VELOCITY DISTRIBUTION
C
CALL VELOCV(FG,FFG,E,COSX,COSY,COSZ,SINE,NUM)
DO IV=1,NUM
EMX = EMX+E(IV)
ENDDO
DO iv=1,num
ne = IDINT(DMIN1(5000.D0,e(iv)+1.D0))
C me(ne) is undefined here!
write (*,*) 'me=',me(ne)
me(ne) = me(ne)+1
ENDDO
GO TO 56
C
C MAXWELLIAN ENERGY DISTRIBUTION
C
ELSE
CALL ENERGV(FE,E,COSX,COSY,COSZ,SINE,NUM)
DO IV=1,NUM
@@ -823,8 +821,6 @@ C
Y(IV) = LM(JL) *RA *COSY(IV)
Z(IV) = LM(JL) *RA *COSZ(IV)
PL(IV) = CVMGT(0.D0,LM(JL)*RA,X0.LE.0.0)
ENDDO
DO IV=1,NUM
LLL(IV) = JL
ENDDO
C
@@ -852,10 +848,7 @@ C
DO IV=1,IH1
call ranlux(random, 1)
JJJ(IV) = ISRCHFGE(NJ(LLL(IV)),COM(1,LLL(IV)),1,
& DBLE(random(1)))
& +JT(LLL(IV))
ENDDO
DO IV=1,IH1
& DBLE(random(1)))+JT(LLL(IV))
EPS(IV)=E(IV)*F1(JJJ(IV))
ENDDO
DO IV=1,IH1
@@ -1149,8 +1142,6 @@ C
ICDI=ICDI+IDINT(CVMGT(1.D0,0.D0,DEN(IV).GT.DI(JJJ(IV))))
ICSUMS=ICSUMS+IDINT(CVMGT(1.D0,0.D0,DEN(IV).GT.SB(1)))
ICSUM=ICSUM+IDINT(CVMGT(1.D0,0.D0,DENS(IV).GT.0.D0))
ENDDO
DO IV=1,IH1
DEN2=DEN(IV)*DEN(IV)
DEN3=DEN2*DEN(IV)
EEL=EEL+DEN(IV)
@@ -1281,9 +1272,6 @@ C
& *CSYR(IREC1,2)
ZR(IREC1,2)=ZR(IREC1,2)+(LM(LRR(IREC1,2)) +TAUPSR(IREC1,2))
& *CSZR(IREC1,2)
ENDDO
DO IREC1=1,NREC2
CXR(IREC1)=CSXR(IREC1,2)
CYR(IREC1)=CSYR(IREC1,2)
CZR(IREC1)=CSZR(IREC1,2)
@@ -1564,8 +1552,6 @@ C
ICDIRJ(JJR(IREC1,2),JJR(IREC1,1))= ICDIRJ(JJR(IREC1,2)
& ,JJR(IREC1,1)) +IDINT(CVMGT(1.D0,0.D0,T(IREC1).GT
& .DI(JJR(IREC1,1))))
ENDDO
DO IREC1=1,NREC2
TEL=TEL+TS(IREC1)
TR=TR+T(IREC1)
EI=EI+DEER(IREC1)
@@ -1853,7 +1839,7 @@ C
I1 = MAX0( MIN0( IDINT(X(IV)/CW+1.D0), MAXD1), 0)
IRP(I1)=IRP(I1)+1
C
C Calculated the stopped particles in each layer
C Calculate the stopped particles in each layer
C
LowTiefe = 0.D0
UpTiefe = DX(1)
@@ -2250,8 +2236,6 @@ C
IF(IVMAX.LT.IVMIN) GO TO 132
DO IV=IVMAX+1,IH1
LLL(IV) = MIN0(ISRCHFGT(L,XX(1),1,X(IV)),L)
ENDDO
DO IV=IVMAX+1,IH1
X(IV)=X(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSX(IV)
Y(IV)=Y(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSY(IV)
Z(IV)=Z(IV)+(LM(LLL(IV))+TAUPSI(IV))*COSZ(IV)
@@ -2272,14 +2256,14 @@ C
ima = 1
1010 ima = MIN0(ima+2,5000)
C I am not sure what is this file for?!
open(20,file='edist')
do ne=1,ima
write(20,1020) ne,me(ne)
enddo
1020 format(1x,2i6)
close(20)
C open(20,file='edist')
C do ne=1,ima
C write(20,1020) ne,me(ne)
C enddo
C 1020 format(1x,2i6)
C close(20)
C
C Berechnung der part. reflec. coeff. nach Thomas et al.
C Calculate the part. reflec. coeff. according to Thomas et al.
C
E0keV=E0/1.D3
EsigkeV=Esig/1.D3
@@ -2494,8 +2478,6 @@ C
DO J=1,LJ
ISPA = ISPA+IBSP(J)
ESPA = ESPA+EBSP(J)
ENDDO
DO J=1,LJ
ISPAT = ISPAT+ITSP(J)
ESPAT = ESPAT+ETSP(J)
ENDDO
@@ -2515,7 +2497,9 @@ C
nli(i)=IDINT(xx(i)/cw+0.01)
enddo
do i=1,l
C we have a problem here, what is nli(i-1) for i=1 (assumed 0??)
do j=nli(i-1)+1,nli(i)
write(*,*) 'i=',i,'j=',j,nli(i-1)+1,nli(i)
irpl(i)=irpl(i)+irp(j)
enddo
enddo
@@ -3388,9 +3372,11 @@ C
& ,(RKDSTJ(I,J),I=1,20)
ENDDO
1800 CONTINUE
c
c The file for33 is created here
c
C
C The file fort.33 is created here
C This is the summary file writing procedure
C
7802 FORMAT(6x,'Energy',4x,'SigmaE',5x,'Alpha',2x,'SigAlpha',4x,'ntot',
& 5x,'imp',2x,'backsc',3x,'trans',3x,'tried',4x,'negE',3x,
& 'range',6x,'straggeling',2x, 'Eback',7x,'sigEback',4x
@@ -3403,8 +3389,6 @@ C & I0,3x))
inquire(FILE='fort.33',EXIST=FORT33)
if (.not.FORT33) then
open(33)
C This is where we can insert the chemical formula
C WRITE(33,7802) ('impL',k,k=1,NLayers)
WRITE(33,7802) (chem(k),k=1,NLayers)
else
open(33,access='append')
@@ -3413,8 +3397,9 @@ C WRITE(33,7802) ('impL',k,k=1,NLayers)
& ,negE,FIX0,SIGMAX,FIB0,SIGMAB,FIT0,SIGMAT,epsilon,prcoeff
& ,(number_in_layer(k),k=1,NLayers)
CLOSE(33)
c
c End of file for33
C
C End of file fort.33
C
C TOP AND FRONT LINES FOR MATRICES
C
trimspNL
BIN
View File
Binary file not shown.