LMU/musrsim
LMU/musrsim
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Cleanup code in progress.

Browse Source
This commit is contained in:
salman 2013-02-20 14:53:56 +00:00
parent f9269571f5
commit abca9a346e
1 changed files with 276 additions and 319 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

163
trimsp/src/trimsp7l.F
View File

@ -149,7 +149,7 @@ C in the 100 depth intervals agrees wit
C number of particles stopped in the different layers
C IF 0 not, message will be written in the range file
C UNIT 21 and UNIT22
c November 2000. Tanya Riseman
c Nov 2000: Tanya Riseman
c Starts porting program to Open VMS and DEC UNIX.
c Compile options: f90/list/warn/ext
c Minor problems with continuation characters in wrong column and
@ -162,7 +162,7 @@ c column 6. Straddling can be non-portable!
c Comment out !DEC$REAL:8 because all
c reals are already declared REAL*8. Add implicit none.
c
c June 2002. Thomas Prokscha PSI
c Jun 2002: Thomas Prokscha PSI
c Stopped porting to f90, use f77 only.
c replaced all non-standard f77 function by standard functions.
c use ranlux random number generator from the CERN library (libmathlib.a must
@ -170,26 +170,19 @@ c be installed) to get rid of machine specific random number generators.
c Add pre-compiler instructions for making different output for the .rge files
c on Windows and Unix.
c
c October 2002 Thomas Prokscha PSI
c Oct 2002: Thomas Prokscha PSI
c corrected error in the calculation of the Thomas-Fermi reduced energy:
c it was 1/(Z1 Z2 * sqrt( Z1**2/3 + Z2**2/3))
c it must be:
c 1/(Z1 Z2 * sqrt( Z1**(2/3) + Z2**(2/3)))
c
c August 2009 Zaher Salman PSI
c changed input file reading to non-formatted. This still works
c Aug 2009: Zaher Salman PSI
c Changed input file reading to non-formatted. This still works
c with old input files, but makes it much easier to create new
c input files without the hassle of formatting.
c Included the source of ranlux for random numbers generation
c into trimsp7l source. No need for cern libraries to be installed.
c
CDIR$ NOLIST
C
cTR !DEC$REAL:8
C
C IMPLICIT INTEGER (i-j)
C IMPLICIT REAL*8 (a-h,k-z)
C
c-------------------------------------------
c check OS
c
@ -237,9 +230,7 @@ c
# ,MEAGS(102,12,22,30)
# ,MEASL(75,21,30)
INTEGER*4 MEAST(102,22,30),MAGST(62,22,30)
C # ,MEAGST(102,36,22,10) von Eckstein herauskommentiert
# ,MEASTL(75,21,30)
C REAL*8 MEAGSL(75,36,21),EAGSL(75) von Eckstein herauskommentiert
INTEGER*4 NJ(7),JT(7),ILD(7)
INTEGER*4 LLL(64),JJJ(64),IK(64)
INTEGER*4 me(5000),nli(0:7),irpl(7)
@ -257,7 +248,7 @@ C REAL*8 MEAGSL(75,36,21),EAGSL(75) von Eckstein herauskommentiert
INTEGER*4 IPB1,KIB,IPT,IE,IERLOG,IAGB,KIT,IMA,IIM
INTEGER*4 im1,im2,im3,IG2,ies,ias
INTEGER*4 JE,JA,JG,JTJ,JTK,JTL
C REAL Variablen
C REAL Variables
REAL*8 CVMGT
REAL*8 X(64),Y(64),Z(64),E(64),PL(64)
# ,COSX(64),COSY(64),COSZ(64),SINE(64)
@ -320,13 +311,11 @@ C REAL Variablen
# ,KL(35,35),KOR(35,35),KLM(7,35)
REAL*8 MU1(35),EC1(35),A1(35),F1(35),KL1(35),KOR1(35)
# ,DI(35),EP(35),ZZ(35),TM(35)
C REAL*8 SL(64,5),SM(64,5),SH(64,5),CH(92,12),CHE(92,9)
REAL*8 CH1(7,5),CH2(7,5),CH3(7,5),CH4(7,5),CH5(7,5)
REAL*8 CHM1(7)
REAL*8 SM(64),SH(64,5)
REAL*8 FIESB(30),SEESB(30),THESB(30),FOESB(30)
# ,SGMESB(30),DFIESB(30),DSEESB(30),DTHESB(30)
C REAL*8 ESVDL(2000)
REAL*8 pi,c,E0,de,alfa,z1,rtheta,cu,enot,esb,est,esp
REAL*8 E2,AB,FP,AN
REAL*8 Esig,Epar
@ -406,7 +395,7 @@ c fuer part. reflec. coeff. Berechnung
REAL*8 PL1S,PL2S,PL3S,PL4S,PL5S,PL6S
REAL*8 YH,HN,CST,BI,BIL,YSP,YSPL,EEE
REAL*4 random,ran2(2)
C CHARACTER Variablen
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
@ -416,10 +405,10 @@ C CHARACTER Variablen
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
CHARACTER min_start*2,min_stop*2,sec_start*2,sec_stop*2
C
COMMON /A/ M1,VELC,ZARG
COMMON /B/ TI,SHEATH,CALFA
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'/
@ -477,7 +466,6 @@ C
DATA ICDR/3500*0/,ICDIRN/3500*0/,IONR/100*0.D0/
DATA DENTR/100*0.D0/,DMGNR/100*0.D0/
DATA IPL/100*0/,IPLB/100*0/,IPLT/100*0/
c DATA IRP/102*0/ gibt witzigweise einen Fehler, aber warum????
DATA IRPL/7*0/
DATA ICDJT/35*0/,ICDJTR/35*0/,ICDITR/35*0/
DATA ICD/3500*0/,ELP/3500*0.D0/,ELD/3500*0.D0/
@ -498,21 +486,18 @@ c DATA IRP/102*0/ gibt witzigweise einen Fehler, aber warum????
DATA ELETR/35*0.D0/,ELITR/35*0.D0/,ELPTR/35*0.D0/
DATA ELDTR/35*0.D0/
DATA Epar/0.D0/
c part. refl. coeff. from Thomas et al.
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 /7*0/
C
C EXTERNAL CVMGT,ILLZ,SCOPY,ISRCHEQ,ISRCHFGE,ISRCHFGT
C
innam=filein//inext
outnam=fileout//outext
rgenam=fileout//rgeext
errnam=fileout//errext
C
OPEN(UNIT=99,file=errnam,STATUS='new')
OPEN(UNIT=11,file=innam,STATUS='unknown',ERR=13591)
C OPEN(UNIT=31,NAME='energ.dat',STATUS='new')
C
READ(11,*) Z1,M1,E0,Esig,ALPHA,ALPHASIG,EF,ESB,SHEATH,ERC
READ(11,*) NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,KK0R,KDEE1,KDEE2
# ,IPOT,IPOTR,IRL
@ -546,24 +531,15 @@ C
WRITE(*,*) ZT(I,1),ZT(I,2),ZT(I,3),ZT(I,4),ZT(I,5)
WRITE(*,*) MT(I,1),MT(I,2),MT(I,3),MT(I,4),MT(I,5)
WRITE(*,*) CO(I,1),CO(I,2),CO(I,3),CO(I,4),CO(I,5)
c WRITE(*,103) SBE(I,1),SBE(I,2),SBE(I,3),SBE(I,4),SBE(I,5)
c WRITE(*,103) ED(I,1),ED(I,2),ED(I,3),ED(I,4),ED(I,5)
c WRITE(*,103) BE(I,1),BE(I,2),BE(I,3),BE(I,4),BE(I,5)
c WRITE(*,107) CH1(I,1),CH1(I,2),CH1(I,3),CH1(I,4),CH1(I,5)
c WRITE(*,107) CH2(I,1),CH2(I,2),CH2(I,3),CH2(I,4),CH2(I,5)
c WRITE(*,107) CH3(I,1),CH3(I,2),CH3(I,3),CH3(I,4),CH3(I,5)
c WRITE(*,107) CH4(I,1),CH4(I,2),CH4(I,3),CH4(I,4),CH4(I,5)
c WRITE(*,107) CH5(I,1),CH5(I,2),CH5(I,3),CH5(I,4),CH5(I,5)
1359 CONTINUE
C
100 FORMAT(2F7.2,1F12.2,7F9.2)
101 FORMAT(I9,3F8.0,1F7.2,1F7.0,2F7.2,6I4,I3)
102 FORMAT(7F9.2,14F7.2)
103 FORMAT(5F9.4)
107 FORMAT(5F12.6)
C
C open statement for output files, removed from line 2449 ff to here
C
OPEN(UNIT=21,FILE=outnam,STATUS='new',ERR=6000)
GOTO 6001
6000 WRITE(*,*)' File schon vorhanden, Gib neue Ausgabedatei an (A8)'
@ -574,13 +550,11 @@ C
6001 OPEN(UNIT=22,FILE=rgenam,STATUS='new')
WRITE(21,1000)
1000 FORMAT(1H1/,6X,'* PROGRAM TRVMC95 - V TrimSP7L 17.Oct.02 TP *')
c TP 1000 FORMAT(1H1/,6X,'* PROGRAM TRVMC95 - Vers. TrimSP7L 22.Nov.00 *')
c TR 1000 FORMAT(1H1/,6X,'** PROGRAM TRVMC95 - V TrimSP7L 26.06.00 **')
C
C 1st CALL DATE_AND_TIME
CALL DATE_AND_TIME(Real_Clock(1),Real_Clock(2),Real_Clock(3),
# Date_Time)
C
IF(Date_Time(2).EQ.1) THEN
month_start='Jan.'
days_start_total=Date_Time(3)
@ -621,13 +595,13 @@ C
C in seconds from beginning of year
seconds_start_total=Date_Time(7)+(Date_Time(6)*60)+
# (Date_Time(5)*3600)+(days_start_total-1)*86400
C
READ(Real_Clock(1)(1:4),'(A4)')year_start
READ(Real_Clock(1)(7:8),'(A2)')day_start
READ(Real_Clock(2)(1:2),'(A2)')hour_start
READ(Real_Clock(2)(3:4),'(A2)')min_start
READ(Real_Clock(2)(5:6),'(A2)')sec_start
C
WRITE(21,*)
WRITE(21,10050)day_start,month_start,year_start,
# hour_start,min_start,sec_start
@ -635,7 +609,6 @@ C
# 1x,A2,':',A2,':',A2)
C SET INTERVAL CONSTANTS FOR OUTPUT
C
DE = 1.D0
DA = 3.D0
DG = 3.D0
@ -656,18 +629,17 @@ C
DAW = BW/DA
DGW = BW/DG
DGIK = BW/DGI
C
C CALCULATION OF CHARGE AND MASS DEPENDENT CONSTANTS
C
PI2=2.D0*PI
ABC=AB*FP
LMAX=7
JMAX=5
L=ISRCHEQ(LMAX,DX(1),1,0.D0)-1
C
C Checks wether depth interval is an integer denominator of layer thickness or not
C If not, calculated implantation profile is not correct.
C
depth_interval_flag = 1
LOOP_Check_layer_thick : DO K=1,L-1
IF(.NOT.EQUAL(DX(K)/CW-DBLE(IDINT(DX(K)/CW)),0.D0)) THEN
@ -679,7 +651,6 @@ C
DO 165 I=1,L
DO 155 J=1,JMAX
IF(EQUAL(CO(I,J),0.D0)) GOTO 156
C IF(CO(I,J).D0EQ.D00.D000) GO TO 156
155 CONTINUE
J=JMAX+1
156 NJ(I)=J-1
@ -819,9 +790,7 @@ C ZBL POTENTIAL (IPOTR=3)
C KLM(LL,I) = KLM(LL,I)+CK(LL)*CO(LL,J)*KL(I,J+JT(LL))
KLM(LL,I) = KLM(LL,I)+CO(LL,J)*KL(I,J+JT(LL))
193 CONTINUE
C
C ALPHA = CVMGT( .001, ALPHA, ALPHA.EQ.0. )
C ALPHA = CVMGT( 179.999, ALPHA, ALPHA.EQ.180.)
ALPHA = CVMGT( .001D0, ALPHA, EQUAL(ALPHA,0.D0))
ALPHA = CVMGT( 179.999D0, ALPHA, EQUAL(ALPHA,180.D0))
@ -831,9 +800,9 @@ C ALPHA = CVMGT( 179.999, ALPHA, ALPHA.EQ.180.)
8883 FORMAT(1X,'ERROR : IF ALPHA.GE.90. THEN IT MUST BE X0.LE.0.')
GO TO 8000
8882 CONTINUE
C
C SET CONSTANT DISTANCES
C
TT = XX(L)
INEL = 0
HLM = CVMGT( 0.D0, -.5D0*LM(1), INEL.EQ.0)
@ -848,7 +817,7 @@ C
SUT = DMAX1(SUTR,DMAX1(SUT1,SUT2))
XC = CVMGT( X0, -SU, X0.GT.0.D0)
RT = TT-RD
C
IF(E0.GE.0.D0) GO TO 51
C
C SET CONSTANTS FOR MAXWELLIAN DISTRIBUTION
@ -868,13 +837,10 @@ C
IY = IDINT(RI)
IY2 = IDINT(RI2)
IY3 = IDINT(RI3)
CC ANFANG = RANSET(IY)
CC ANFANG = SRAND48(IY)
ISEED = IY
ISEED2 = IY2
ISEED3 = IY3
C WRITE(*,*)ISEED
C
IF ( E0.GT.0.D0 ) GO TO 47
IF ( ALPHA.GE.0.D0 ) THEN
C
@ -908,7 +874,6 @@ C
C FIXED PROJECTILE ENERGY
DO IV=1,NUM
E(IV) = E0
C WRITE(*,*)' Da Esig=0 ist E=E0'
ENDDO
ELSE
C GAUSSIAN ENERGY DISTRIBUTION
@ -920,16 +885,13 @@ C GAUSSIAN ENERGY DISTRIBUTION
GO TO 5200
ENDIF
E(IV)=Epar
C WRITE(*,*)E(IV),Epar,E0
ENDDO
ENDIF
C
C die nachfolgende Zeile wurden von Linie 633 hier hin verschoben
C
SFE = DMIN1(SB(1),SB(L))
C
IF ( ALPHA.GE.0.D0 ) THEN
C
IF(EQUAL(ALPHASIG,0.D0))THEN
C
C FIXED PROJECTILE ANGLE
@ -967,22 +929,16 @@ C
C COSINE ANGLE DISTRIBUTION (THREE-DIMENSIONAL)
C
CDIR$ IVDEP
DO 53 IV=1,NUM
CC RPHI = PI2*RANF()
CC RPHI = PI2*DRAND48()
CC RPHI = PI2*DBLE(RAN(ISEED))
DO IV=1,NUM
call ranlux(ran2,2)
RPHI = PI2*DBLE(ran2(1))
CC RTHETA = RANF()
CC RTHETA = DRAND48()
CC RTHETA = DBLE(RAN(ISEED))
RTHETA = DBLE(ran2(2))
COSX(IV) = DSQRT(RTHETA)
SINE(IV) = DSQRT(1.D0-RTHETA)
COSY(IV) = SINE(IV)*DCOS(RPHI)
COSZ(IV) = SINE(IV)*DSIN(RPHI)
53 CONTINUE
C
ENDDO
ELSEIF (EQUAL(ALPHA,-1.D0).AND.X0.GT.0.D0) THEN
C ELSEIF ( ALPHA.EQ.-1. AND. X0.GT.0. ) THEN
C
@ -1331,45 +1287,46 @@ C
DO 73 J=1,NJ(LLL(IV))
SM(IV)=SM(IV)+SH(IV,J)*CO(LLL(IV),J)
73 CONTINUE
DO 78 IV=1,IH1
DEE(IV)=CVMGT(CHM1(LLL(IV))*DSQRT(EM(IV)),SM(IV),EM(IV).LE.10.D0)
78 CONTINUE
DO 69 IV=1,IH1
DEE(IV)=10.D0*ASIGT(IV)*
# CVMGT(0.D0,DEE(IV),X(IV).LT.HLM.OR.X(IV).GT.HLMT)
69 CONTINUE
DO IV=1,IH1
DEE(IV)=CVMGT(CHM1(LLL(IV))*DSQRT(EM(IV)),SM(IV),EM(IV).LE.10
& .D0)
ENDDO
DO IV=1,IH1
DEE(IV)=10.D0*ASIGT(IV)*CVMGT(0.D0,DEE(IV),X(IV).LT.HLM.OR
& .X(IV).GT.HLMT)
ENDDO
GO TO 40
19 FHE=CVMGT(1.3333D0,1.D0,M1.LT.4.00D0)
DO 25 IV=1,IH1
DO IV=1,IH1
SM(IV)=0.D0
EM(IV)=E(IV)*0.001D0*FHE
25 CONTINUE
DO 74 IV=1,IH1
DO 74 J=1,NJ(LLL(IV))
SH(IV,J)=CH1(LLL(IV),J)*EM(IV)**CH2(LLL(IV),J)*
# (CH3(LLL(IV),J)/(EM(IV)*0.001D0))
# *DLOG(DABS(1.D0+(CH4(LLL(IV),J)/(EM(IV)*0.001D0))+
# CH5(LLL(IV),J)*EM(IV)*0.001D0))
# /(CH1(LLL(IV),J)*EM(IV)**CH2(LLL(IV),J)+
# (CH3(LLL(IV),J)/(EM(IV)*0.001D0))
# *DLOG(DABS(1.D0+(CH4(LLL(IV),J)/(EM(IV)*0.001D0))+
# CH5(LLL(IV),J)*EM(IV)*0.001D0)))
74 CONTINUE
DO 92 IV=1,IH1
DO 92 J=1,NJ(LLL(IV))
ENDDO
DO IV=1,IH1
DO J=1,NJ(LLL(IV))
SH(IV,J)=CH1(LLL(IV),J)*EM(IV)**CH2(LLL(IV),J)* (CH3(LLL(IV)
& ,J)/(EM(IV)*0.001D0)) *DLOG(DABS(1.D0+(CH4(LLL(IV),J)
& /(EM(IV)*0.001D0))+ CH5(LLL(IV),J)*EM(IV)*0.001D0))
& /(CH1(LLL(IV),J)*EM(IV)**CH2(LLL(IV),J)+ (CH3(LLL(IV),J
& )/(EM(IV)*0.001D0)) *DLOG(DABS(1.D0+(CH4(LLL(IV),J)
& /(EM(IV)*0.001D0))+ CH5(LLL(IV),J)*EM(IV)*0.001D0)))
ENDDO
ENDDO
DO IV=1,IH1
DO J=1,NJ(LLL(IV))
SM(IV)=SM(IV)+SH(IV,J)*CO(LLL(IV),J)
92 CONTINUE
DO 79 IV=1,IH1
DEE(IV)=10.D0*ASIGT(IV)*
# CVMGT(0.D0,SM(IV),X(IV).LT.HLM.OR.X(IV).GT.HLMT)
79 CONTINUE
ENDDO
ENDDO
DO IV=1,IH1
DEE(IV)=10.D0*ASIGT(IV)* CVMGT(0.D0,SM(IV),X(IV).LT.HLM.OR.X(IV
& ).GT.HLMT)
ENDDO
40 CONTINUE
C
DO 44 IV=1,IH1
DO IV=1,IH1
DEL=DMAX1(1.0D-20,DENS(IV)+DEE(IV))
DENS(IV)=CVMGT(E(IV)*DENS(IV)/DEL,DENS(IV),DEL.GT.E(IV))
DEE(IV)=CVMGT(E(IV)*DEE(IV)/DEL,DEE(IV),DEL.GT.E(IV))
44 CONTINUE
ENDDO
C
C INCREMENT OF DAMAGE, CASCADE AND PHONON ENERGY
C