From a1b7d7365a7c951f89352df2efe046075786dd89 Mon Sep 17 00:00:00 2001 From: Zaher Salman Date: Thu, 21 Feb 2013 13:06:27 +0000 Subject: [PATCH] Cleanup code in progress. --- trimsp/src/trimsp7l.F | 893 ++++++++++++++++++++---------------------- 1 file changed, 423 insertions(+), 470 deletions(-) diff --git a/trimsp/src/trimsp7l.F b/trimsp/src/trimsp7l.F index 38c1fe5..2c03b9b 100644 --- a/trimsp/src/trimsp7l.F +++ b/trimsp/src/trimsp7l.F @@ -191,7 +191,7 @@ c #else #define OS_UNIX #endif -c + IMPLICIT NONE LOGICAL TEST(64),TESTR(2000),TEST1(2000) LOGICAL EQUAL @@ -211,26 +211,24 @@ c INTEGER*4 IDMAX(2000),IKR(2000) INTEGER*4 number_in_layer(7),laufzahl INTEGER*4 IRP(0:101),IPL(100),IPLB(100),IPLT(100) - INTEGER*4 ICD(100,35),ICDT(100),ICDJT(35),ICDIRJ(35,35) - # ,ICDR(100,35),ICDTR(100),ICDJTR(35) - # ,ICDIRI(100,35,35),ICDIRN(100,35),ICDITR(35) - INTEGER*4 KADB(20),KADT(20),KADS(20),KADST(20) - # ,KADRIP(20,30),KADRIS(20,30),KADROP(20,30),KADROS(20,30) - # ,KADSJ(20,30),KADSL(20,6),KDSTJ(20,30),KDSTL(20,6) - INTEGER*4 IBSP(35),ISPAL(7),IBSPL(35) - # ,ISPIP(35),ISPIS(35),ISPOP(35),ISPOS(35) + INTEGER*4 ICD(100,35),ICDT(100),ICDJT(35),ICDIRJ(35,35),ICDR(100 + & ,35),ICDTR(100),ICDJTR(35) ,ICDIRI(100,35,35),ICDIRN(100,35) + & ,ICDITR(35) + INTEGER*4 KADB(20),KADT(20),KADS(20),KADST(20) ,KADRIP(20,30) + & ,KADRIS(20,30),KADROP(20,30),KADROS(20,30) ,KADSJ(20,30) + & ,KADSL(20,6),KDSTJ(20,30),KDSTL(20,6) + INTEGER*4 IBSP(35),ISPAL(7),IBSPL(35) ,ISPIP(35),ISPIS(35) + & ,ISPOP(35),ISPOS(35) INTEGER*4 ITSP(35),ISPALT(7) - # ,ISPIPT(35),ISPIST(35),ISPOPT(35),ISPOST(35) + & ,ISPIPT(35),ISPIST(35),ISPOPT(35),ISPOST(35) INTEGER*4 KO(600,35,2) - INTEGER*4 MEAB(102,22),MAGB(62,22),MEAGB(102,36,22) - # ,MEABL(75,21),MEPB(102,102) + INTEGER*4 MEAB(102,22),MAGB(62,22),MEAGB(102,36,22) ,MEABL(75,21) + & ,MEPB(102,102) INTEGER*4 MEAT(102,22),MAGT(62,22),MEAGT(102,36,22), - # MEATL(75,21),MEPT(102,102) + & MEATL(75,21),MEPT(102,102) INTEGER*4 MEAS(102,22,30),MAGS(62,22,30),MAGSA(62,32,30) - # ,MEAGS(102,12,22,30) - # ,MEASL(75,21,30) - INTEGER*4 MEAST(102,22,30),MAGST(62,22,30) - # ,MEASTL(75,21,30) + & ,MEAGS(102,12,22,30) ,MEASL(75,21,30) + INTEGER*4 MEAST(102,22,30),MAGST(62,22,30) ,MEASTL(75,21,30) 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) @@ -250,72 +248,62 @@ c INTEGER*4 JE,JA,JG,JTJ,JTK,JTL 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) - REAL*8 EPS(64),DEN(64),DEE(64),DENS(64),DEES(64) - # ,CX(64),CY(64),CZ(64),SX(64),X1(64),ASIGT(64),EM(64) - REAL*8 EX1(64),EX2(64),EX3(64),P(64),TAU(64),EX4(64) - # ,B(64),R(64),C2(64),S2(64),CT(64),ST(64),V(64),V1(64) - # ,CPHI(64),SPHI(64),CPSI(64),SPSI(64),TAUPSI(64) - # ,ENUCL(64),EINEL(64),ENUCL2(64),EINEL2(64) - REAL*8 ER(2000,2),XR(2000,2),YR(2000,2),ZR(2000,2) - # ,CSXR(2000,2),CSYR(2000,2),CSZR(2000,2),TAUR(2000) - # ,SNXR(2000,2),CPSIR(2000,2),SPSIR(2000,2),CPHIR(2000,2) - # ,SPHIR(2000,2),TAUPSR(2000,2) + REAL*8 X(64),Y(64),Z(64),E(64),PL(64) ,COSX(64),COSY(64),COSZ(64) + & ,SINE(64) + REAL*8 EPS(64),DEN(64),DEE(64),DENS(64),DEES(64) ,CX(64),CY(64) + & ,CZ(64),SX(64),X1(64),ASIGT(64),EM(64) + REAL*8 EX1(64),EX2(64),EX3(64),P(64),TAU(64),EX4(64) ,B(64),R(64) + & ,C2(64),S2(64),CT(64),ST(64),V(64),V1(64) ,CPHI(64),SPHI(64) + & ,CPSI(64),SPSI(64),TAUPSI(64) ,ENUCL(64),EINEL(64),ENUCL2(64) + & ,EINEL2(64) + REAL*8 ER(2000,2),XR(2000,2),YR(2000,2),ZR(2000,2) ,CSXR(2000,2) + & ,CSYR(2000,2),CSZR(2000,2),TAUR(2000) ,SNXR(2000,2) + & ,CPSIR(2000,2),SPSIR(2000,2),CPHIR(2000,2) ,SPHIR(2000,2) + & ,TAUPSR(2000,2) REAL*8 EPSR(2000),T(2000),TS(2000),DEER(2000),DEERS(2000) - # ,PR(2000),BR(2000),EX1R(2000),EX2R(2000),EX3R(2000) - # ,CTR(2000),STR(2000),ASIGTR(2000),EX4R(2000) - # ,X2(2000),RR(2000),VR(2000) - # ,V1R(2000),CXR(2000),CYR(2000),CZR(2000) - # ,SXR(2000),C2R(2000),S2R(2000),CUR(2000) - REAL*8 RIRP(0:101) - # ,CASMOT(100),PHON(100),DENT(100),ION(100),DMGN(100) - # ,CASMOTR(100),PHONR(100),DENTR(100),IONR(100),DMGNR(100) - # ,ELGD(100),ELGDR(100) - REAL*8 ELE(100,35),ELI(100,35),ELP(100,35),ELD(100,35) - # ,ELET(35),ELIT(35),ELPT(35),ELDT(35) - # ,ELER(100,35),ELIR(100,35),ELPR(100,35),ELDR(100,35) - # ,ELETR(35),ELITR(35),ELPTR(35),ELDTR(35) - REAL*8 AI(20),RKADB(20),RKADT(20) - # ,RKADS(20),RKADST(20) - # ,RKADSJ(20,30),RKADSL(20,7) - # ,RKDSTJ(20,30),RKDSTL(20,7) - REAL*8 EBSP(35),ESPAL(7) - # ,SPY(35),SPE(35),REY(35),EMSP(35) - # ,ESPIP(35),ESPIS(35),ESPOP(35),ESPOS(35) - # ,RIP(35),RIS(35),ROP(35),ROS(35) - # ,REIP(35),REIS(35),REOP(35),REOS(35) - # ,ESPMIP(35),ESPMIS(35),ESPMOP(35),ESPMOS(35) - # ,RIPJ(35),RISJ(35),ROPJ(35),ROSJ(35) - # ,REIPJ(35),REISJ(35),REOPJ(35),REOSJ(35) - REAL*8 ETSP(35),ESPALT(7) - # ,SPYT(35),SPET(35),REYT(35),EMSPT(35) - # ,ESPIPT(35),ESPIST(35),ESPOPT(35),ESPOST(35) - # ,RIPT(35),RIST(35),ROPT(35),ROST(35) - # ,REIPT(35),REIST(35),REOPT(35),REOST(35) - # ,ESPMIPT(35),ESPMIST(35),ESPMOPT(35),ESPMOST(35) - REAL*8 SPEM(35),SPE2S(35),SPE3S(35),SPE4S(35),SPE5S(35) - # ,SPE6S(35),VSPE(35),SSPE(35),GSPE(35),BSPE(35) + & ,PR(2000),BR(2000),EX1R(2000),EX2R(2000),EX3R(2000),CTR(2000) + & ,STR(2000),ASIGTR(2000),EX4R(2000) ,X2(2000),RR(2000) + & ,VR(2000) ,V1R(2000),CXR(2000),CYR(2000),CZR(2000) ,SXR(2000) + & ,C2R(2000),S2R(2000),CUR(2000) + REAL*8 RIRP(0:101) ,CASMOT(100),PHON(100),DENT(100),ION(100) + & ,DMGN(100) ,CASMOTR(100),PHONR(100),DENTR(100),IONR(100) + & ,DMGNR(100) ,ELGD(100),ELGDR(100) + REAL*8 ELE(100,35),ELI(100,35),ELP(100,35),ELD(100,35) ,ELET(35) + & ,ELIT(35),ELPT(35),ELDT(35) ,ELER(100,35),ELIR(100,35) + & ,ELPR(100,35),ELDR(100,35) ,ELETR(35),ELITR(35),ELPTR(35) + & ,ELDTR(35) + REAL*8 AI(20),RKADB(20),RKADT(20) ,RKADS(20),RKADST(20) + & ,RKADSJ(20,30),RKADSL(20,7) ,RKDSTJ(20,30),RKDSTL(20,7) + REAL*8 EBSP(35),ESPAL(7) ,SPY(35),SPE(35),REY(35),EMSP(35) + & ,ESPIP(35),ESPIS(35),ESPOP(35),ESPOS(35) ,RIP(35),RIS(35) + & ,ROP(35),ROS(35) ,REIP(35),REIS(35),REOP(35),REOS(35) + & ,ESPMIP(35),ESPMIS(35),ESPMOP(35),ESPMOS(35) ,RIPJ(35) + & ,RISJ(35),ROPJ(35),ROSJ(35) ,REIPJ(35),REISJ(35),REOPJ(35) + & ,REOSJ(35) + REAL*8 ETSP(35),ESPALT(7) ,SPYT(35),SPET(35),REYT(35),EMSPT(35) + & ,ESPIPT(35),ESPIST(35),ESPOPT(35),ESPOST(35) ,RIPT(35) + & ,RIST(35),ROPT(35),ROST(35) ,REIPT(35),REIST(35),REOPT(35) + & ,REOST(35) ,ESPMIPT(35),ESPMIST(35),ESPMOPT(35),ESPMOST(35) + REAL*8 SPEM(35),SPE2S(35),SPE3S(35),SPE4S(35),SPE5S(35) ,SPE6S(35) + & ,VSPE(35),SSPE(35),GSPE(35),BSPE(35) REAL*8 SPE1SL(35),SPE2SL(35),SPE3SL(35),SPE4SL(35),SPE5SL(35) - # ,SPE6SL(35) + & ,SPE6SL(35) REAL*8 ELOG(75),EMA(62,22),EABL(75) REAL*8 EMAT(62,22),EATL(75),EASL(75,30),EASTL(75,30) REAL*8 FG(128),FFG(64) - REAL*8 XX(7),DX(7),RHO(7),Z2(7),M2(7),LM(7),PDMAX(7) - # ,ARHO(7),AM(7),FM(7),EPS0(7),ASIG(7),K2(7),CK(7) - # ,KLM1(7),SB(7),DLI(7) + REAL*8 XX(7),DX(7),RHO(7),Z2(7),M2(7),LM(7),PDMAX(7),ARHO(7),AM(7) + & ,FM(7),EPS0(7),ASIG(7),K2(7),CK(7) ,KLM1(7),SB(7),DLI(7) REAL*8 UpTiefe,LowTiefe - REAL*8 ZT(7,5),MT(7,5),CO(7,5),SBE(7,5),ED(7,5),BE(7,5) - # ,COM(5,7) - REAL*8 MU(35,35),EC(35,35),A(35,35),F(35,35) - # ,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) + REAL*8 ZT(7,5),MT(7,5),CO(7,5),SBE(7,5),ED(7,5),BE(7,5) ,COM(5,7) + REAL*8 MU(35,35),EC(35,35),A(35,35),F(35,35) ,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) 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) + REAL*8 FIESB(30),SEESB(30),THESB(30),FOESB(30) ,SGMESB(30) + & ,DFIESB(30),DSEESB(30),DTHESB(30) REAL*8 pi,c,E0,de,alfa,z1,rtheta,cu,enot,esb,est,esp REAL*8 E2,AB,FP,AN REAL*8 Esig,Epar @@ -381,8 +369,8 @@ c fuer part. reflec. coeff. Berechnung REAL*8 FIT0,SET,THT,FOT,FIT,SIT,SIGMAT,DFIT0,DSET,DTHT REAL*8 FIPT0,SEPT,THPT,FOPT,FIPT,SIPT,SIGMPT,DFIPT0,DSEPT,DTHPT REAL*8 FIES0,SEES,THES,FOES,FIES,SIES,SIGMES,DFIES0,DSEES,DTHES - REAL*8 FIES0L,SEESL,THESL,FOESL,FIESL,SIESL,SIGMSL, - # DFIESL,DSEESL,DTHESL + REAL*8 FIES0L,SEESL,THESL,FOESL,FIESL,SIESL,SIGMSL,DFIESL,DSEESL + & ,DTHESL REAL*8 X1SD,X2SD,X3SD,X4SD,X5SD,X6SD REAL*8 ACSUMR,ACDISR,ACSBER,ACSUR,ACDIR,ACSBR REAL*8 ACDR11,ACDR12,ACDR21,ACDR22 @@ -499,12 +487,12 @@ C part. refl. coeff. from Thomas et al. OPEN(UNIT=11,file=innam,STATUS='unknown',ERR=13591) 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 - READ(11,*) DX(1),DX(2),DX(3),DX(4),DX(5),DX(6),DX(7), - # RHO(1),RHO(2),RHO(3),RHO(4),RHO(5),RHO(6),RHO(7), - # CK(1),CK(2),CK(3),CK(4),CK(5),CK(6),CK(7) - DO 135 I=1,7 + READ(11,*) NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,KK0R,KDEE1,KDEE2,IPOT + & ,IPOTR,IRL + READ(11,*) DX(1),DX(2),DX(3),DX(4),DX(5),DX(6),DX(7),RHO(1),RHO(2) + & ,RHO(3),RHO(4),RHO(5),RHO(6),RHO(7), CK(1),CK(2),CK(3),CK(4) + & ,CK(5),CK(6),CK(7) + DO I=1,7 READ(11,*) ZT(I,1),ZT(I,2),ZT(I,3),ZT(I,4),ZT(I,5) READ(11,*) MT(I,1),MT(I,2),MT(I,3),MT(I,4),MT(I,5) READ(11,*) CO(I,1),CO(I,2),CO(I,3),CO(I,4),CO(I,5) @@ -516,22 +504,22 @@ C part. refl. coeff. from Thomas et al. READ(11,*) CH3(I,1),CH3(I,2),CH3(I,3),CH3(I,4),CH3(I,5) READ(11,*) CH4(I,1),CH4(I,2),CH4(I,3),CH4(I,4),CH4(I,5) READ(11,*) CH5(I,1),CH5(I,2),CH5(I,3),CH5(I,4),CH5(I,5) - 135 CONTINUE + ENDDO 13591 CLOSE(UNIT=21) WRITE(*,*) Z1,M1,E0,Esig,ALPHA,ALPHASIG,EF,ESB,SHEATH,ERC - WRITE(*,*) NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,KK0R,KDEE1,KDEE2 - # ,IPOT,IPOTR,IRL - WRITE(*,*) DX(1),DX(2),DX(3),DX(4),DX(5),DX(6),DX(7), - # RHO(1),RHO(2),RHO(3),RHO(4),RHO(5),RHO(6),RHO(7), - # CK(1),CK(2),CK(3),CK(4),CK(5),CK(6),CK(7) - DO 1359 I=1,7 + WRITE(*,*) NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,KK0R,KDEE1,KDEE2 ,IPOT + & ,IPOTR,IRL + WRITE(*,*) DX(1),DX(2),DX(3),DX(4),DX(5),DX(6),DX(7),RHO(1),RHO(2) + & ,RHO(3),RHO(4),RHO(5),RHO(6),RHO(7), CK(1),CK(2),CK(3),CK(4) + & ,CK(5),CK(6),CK(7) + DO I=1,7 WRITE(*,'(1x,I2,A7)')i,'. Layer' 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) - 1359 CONTINUE + ENDDO 100 FORMAT(2F7.2,1F12.2,7F9.2) 101 FORMAT(I9,3F8.0,1F7.2,1F7.0,2F7.2,6I4,I3) @@ -553,7 +541,7 @@ C open statement for output files, removed from line 2449 ff to here C 1st CALL DATE_AND_TIME CALL DATE_AND_TIME(Real_Clock(1),Real_Clock(2),Real_Clock(3), - # Date_Time) + & Date_Time) IF(Date_Time(2).EQ.1) THEN month_start='Jan.' @@ -593,8 +581,8 @@ C 1st CALL DATE_AND_TIME days_start_total=31+28+31+30+31+30+31+31+30+31+30+Date_Time(3) ENDIF 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 + seconds_start_total=Date_Time(7)+(Date_Time(6)*60)+(Date_Time(5) + & *3600)+(days_start_total-1)*86400 READ(Real_Clock(1)(1:4),'(A4)')year_start READ(Real_Clock(1)(7:8),'(A2)')day_start @@ -603,10 +591,10 @@ C in seconds from beginning of year READ(Real_Clock(2)(5:6),'(A2)')sec_start WRITE(21,*) - WRITE(21,10050)day_start,month_start,year_start, - # hour_start,min_start,sec_start -10050 FORMAT(1x,' TrimSp simulation started at: ',A2,'.',A4,1x,A4, - # 1x,A2,':',A2,':',A2) + WRITE(21,10050)day_start,month_start,year_start,hour_start + & ,min_start,sec_start +10050 FORMAT(1x,' TrimSp simulation started at: ',A2,'.',A4,1x,A4,1x,A2 + & ,':',A2,':',A2) C SET INTERVAL CONSTANTS FOR OUTPUT DE = 1.D0 @@ -647,14 +635,14 @@ C If not, calculated implantation profile is not correct. EXIT LOOP_Check_layer_thick ENDIF ENDDO LOOP_Check_layer_thick -C - DO 165 I=1,L - DO 155 J=1,JMAX + + DO I=1,L + DO J=1,JMAX IF(EQUAL(CO(I,J),0.D0)) GOTO 156 - 155 CONTINUE + ENDDO J=JMAX+1 156 NJ(I)=J-1 - 165 CONTINUE + ENDDO JT(1) = 0 JT(2) = NJ(1) JT(3) = NJ(1)+NJ(2) @@ -664,231 +652,243 @@ C JT(7) = JT(6)+NJ(6) LJ = NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+NJ(7) XX(1)=DX(1) - DO 170 I=2,L - 170 XX(I)=XX(I-1)+DX(I) - DO 180 I=1,L - DO 180 J=1,NJ(I) + DO I=2,L + XX(I)=XX(I-1)+DX(I) + ENDDO + DO I=1,L + 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) - 180 CONTINUE - DO 185 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 - 185 CONTINUE - DO 186 J = 1,NJ(1) - ZZ(J) = ZT(1,J) - TM(J) = MT(1,J) - DI(J) = ED(1,J) - 186 EP(J) = BE(1,J) - DO 187 J = 1,NJ(2) - ZZ(NJ(1)+J) = ZT(2,J) - TM(NJ(1)+J) = MT(2,J) - DI(NJ(1)+J) = ED(2,J) - 187 EP(NJ(1)+J) = BE(2,J) - DO 188 J = 1,NJ(3) - ZZ(NJ(1)+NJ(2)+J) = ZT(3,J) - TM(NJ(1)+NJ(2)+J) = MT(3,J) - DI(NJ(1)+NJ(2)+J) = ED(3,J) - 188 EP(NJ(1)+NJ(2)+J) = BE(3,J) - DO 1880 J = 1,NJ(4) - ZZ(NJ(1)+NJ(2)+NJ(3)+J) = ZT(4,J) - TM(NJ(1)+NJ(2)+NJ(3)+J) = MT(4,J) - DI(NJ(1)+NJ(2)+NJ(3)+J) = ED(4,J) - 1880 EP(NJ(1)+NJ(2)+NJ(3)+J) = BE(4,J) - DO 1881 J = 1,NJ(5) - ZZ(NJ(1)+NJ(2)+NJ(3)+NJ(4)+J) = ZT(5,J) - TM(NJ(1)+NJ(2)+NJ(3)+NJ(4)+J) = MT(5,J) - DI(NJ(1)+NJ(2)+NJ(3)+NJ(4)+J) = ED(5,J) - 1881 EP(NJ(1)+NJ(2)+NJ(3)+NJ(4)+J) = BE(5,J) - DO 1882 J = 1,NJ(6) - ZZ(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+J) = ZT(6,J) - TM(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+J) = MT(6,J) - DI(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+J) = ED(6,J) - 1882 EP(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+J) = BE(6,J) - DO 18803 J = 1,NJ(7) - ZZ(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+J) = ZT(7,J) - TM(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+J) = MT(7,J) - DI(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+J) = ED(7,J) -18803 EP(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+J) = BE(7,J) - DO 189 I=1,L - COM(1,I) = CO(I,1) - DO 189 J=1,NJ(I)-1 - COM(J+1,I) = COM(J,I)+CO(I,J+1) - 189 CONTINUE - DO 190 J = 1,LJ - MU1(J) = M1/TM(J) - EC1(J) = 4.D0*MU1(J)/((1.D0+MU1(J))*(1.D0+MU1(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 + DO J = 1,NJ(1) + ZZ(J) = ZT(1,J) + TM(J) = MT(1,J) + DI(J) = ED(1,J) + EP(J) = BE(1,J) + ENDDO + DO J = 1,NJ(2) + ZZ(NJ(1)+J) = ZT(2,J) + TM(NJ(1)+J) = MT(2,J) + DI(NJ(1)+J) = ED(2,J) + EP(NJ(1)+J) = BE(2,J) + ENDDO + DO J = 1,NJ(3) + ZZ(NJ(1)+NJ(2)+J) = ZT(3,J) + TM(NJ(1)+NJ(2)+J) = MT(3,J) + DI(NJ(1)+NJ(2)+J) = ED(3,J) + EP(NJ(1)+NJ(2)+J) = BE(3,J) + ENDDO + DO J = 1,NJ(4) + ZZ(NJ(1)+NJ(2)+NJ(3)+J) = ZT(4,J) + TM(NJ(1)+NJ(2)+NJ(3)+J) = MT(4,J) + DI(NJ(1)+NJ(2)+NJ(3)+J) = ED(4,J) + EP(NJ(1)+NJ(2)+NJ(3)+J) = BE(4,J) + ENDDO + DO J = 1,NJ(5) + ZZ(NJ(1)+NJ(2)+NJ(3)+NJ(4)+J) = ZT(5,J) + TM(NJ(1)+NJ(2)+NJ(3)+NJ(4)+J) = MT(5,J) + DI(NJ(1)+NJ(2)+NJ(3)+NJ(4)+J) = ED(5,J) + EP(NJ(1)+NJ(2)+NJ(3)+NJ(4)+J) = BE(5,J) + ENDDO + DO J = 1,NJ(6) + ZZ(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+J) = ZT(6,J) + TM(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+J) = MT(6,J) + DI(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+J) = ED(6,J) + EP(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+J) = BE(6,J) + ENDDO + DO J = 1,NJ(7) + ZZ(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+J) = ZT(7,J) + TM(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+J) = MT(7,J) + DI(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+J) = ED(7,J) + EP(NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+J) = BE(7,J) + 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))) 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) - 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)) - 190 CONTINUE - IF(IPOT.EQ.1) THEN -C KR-C POTENTIAL (IPOT=1) - DO 194 J=1,LJ - KOR1(J) = 0.0389205D0*KL1(J)/(PI*A1(J)*A1(J)) - 194 CONTINUE - ELSEIF (IPOT.EQ.2) THEN + A1(J) = CVMGT(CA*ABC*(ZZ(J)**(-1.D0/3.D0)),CA*ABC/(Z1**0.23D0 + & +ZZ(J)**0.23D0),IPOT.LT.3) + 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 +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 C MOLIERE POTENTIAL (IPOT=2) - DO 195 J=1,LJ - KOR1(J) = 0.045D0*KL1(J)/(PI*A1(J)*A1(J)) - 195 CONTINUE - ELSEIF (IPOT.EQ.3) THEN + DO J=1,LJ + KOR1(J) = 0.045D0*KL1(J)/(PI*A1(J)*A1(J)) + ENDDO + ELSEIF (IPOT.EQ.3) THEN C ZBL POTENTIAL - DO 196 J=1,LJ - KOR1(J) = 0.0203253D0*KL1(J)/(PI*A1(J)*A1(J)) - 196 CONTINUE - ENDIF - DO 191 I = 1,LJ - DO 191 J = 1,LJ - MU(I,J) = TM(I)/TM(J) - 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) -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))) - 191 CONTINUE - IF (IPOTR.EQ.1) THEN -C KR-C POTENTIAL (IPOTR=1) - DO 197 I = 1,LJ - DO 197 J = 1,LJ - KOR(I,J) = 0.0389205D0*KL(I,J)/(PI*A(I,J)*A(I,J)) - 197 CONTINUE - ELSEIF (IPOTR.EQ.2) THEN -C MOLIERE POTENTIAL (IPOTR=2) - DO 198 I = 1,LJ - DO 198 J = 1,LJ - KOR(I,J) = 0.045D0*KL(I,J)/(PI*A(I,J)*A(I,J)) - 198 CONTINUE - ELSEIF (IPOTR.EQ.3) THEN -C ZBL POTENTIAL (IPOTR=3) - DO 184 I = 1,LJ - DO 184 J = 1,LJ - KOR(I,J) = 0.0203253D0*KL(I,J)/(PI*A(I,J)*A(I,J)) - 184 CONTINUE - ENDIF - DO 192 LL=1,L - DO 192 J=1,NJ(LL) - KLM1(LL) = KLM1(LL)+CO(LL,J)*KL1(J+JT(LL))*CK(LL) - CHM1(LL) = CHM1(LL)+CO(LL,J)*CH1(LL,J) - SB(LL) = SB(LL)+CO(LL,J)*SBE(LL,J) - 192 CONTINUE - DO 193 I=1,LJ - DO 193 LL = 1,L - DO 193 J=1,NJ(LL) -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 + DO J=1,LJ + KOR1(J) = 0.0203253D0*KL1(J)/(PI*A1(J)*A1(J)) + ENDDO + ENDIF + DO I = 1,LJ + DO J = 1,LJ + MU(I,J) = TM(I)/TM(J) + 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) +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 +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 +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 +C ZBL POTENTIAL (IPOTR=3) + DO I = 1,LJ + DO J = 1,LJ + KOR(I,J) = 0.0203253D0*KL(I,J)/(PI*A(I,J)*A(I,J)) + 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) + CHM1(LL) = CHM1(LL)+CO(LL,J)*CH1(LL,J) + SB(LL) = SB(LL)+CO(LL,J)*SBE(LL,J) + ENDDO + ENDDO + DO I=1,LJ + DO LL = 1,L + DO J=1,NJ(LL) + KLM(LL,I) = KLM(LL,I)+CO(LL,J)*KL(I,J+JT(LL)) + ENDDO + ENDDO + ENDDO - ALPHA = CVMGT( .001D0, ALPHA, EQUAL(ALPHA,0.D0)) - ALPHA = CVMGT( 179.999D0, ALPHA, EQUAL(ALPHA,180.D0)) + ALPHA = CVMGT( .001D0, ALPHA, EQUAL(ALPHA,0.D0)) + ALPHA = CVMGT( 179.999D0, ALPHA, EQUAL(ALPHA,180.D0)) - IF(ALPHA.GE.90.0.AND.X0.LE.0.0) GO TO 8881 - GO TO 8882 - 8881 WRITE(6,8883) - 8883 FORMAT(1X,'ERROR : IF ALPHA.GE.90. THEN IT MUST BE X0.LE.0.') - GO TO 8000 - 8882 CONTINUE + IF(ALPHA.GE.90.0.AND.X0.LE.0.0) GO TO 8881 + GO TO 8882 + 8881 WRITE(6,8883) + 8883 FORMAT(1X,'ERROR : IF ALPHA.GE.90. THEN IT MUST BE X0.LE.0.') + GO TO 8000 + 8882 CONTINUE C SET CONSTANT DISTANCES + + TT = XX(L) + INEL = 0 + HLM = CVMGT( 0.D0, -.5D0*LM(1), INEL.EQ.0) + HLMT = CVMGT( TT, TT+.5D0*LM(L), INEL.EQ.0) + SU1 = PDMAX(1) + PDMAX(1) + SU2 = PDMAX(1)*(1.D0+KK0) + SUR = PDMAX(1)*(1.D0+KK0R) + SU = DMAX1(SUR,DMAX1(SU1,SU2)) + SUT1 = TT + PDMAX(L) + PDMAX(L) + SUT2 = TT + PDMAX(L)*(1.D0+KK0) + SUTR = TT + PDMAX(L)*(1.D0+KK0R) + SUT = DMAX1(SUTR,DMAX1(SUT1,SUT2)) + XC = CVMGT( X0, -SU, X0.GT.0.D0) + RT = TT-RD - TT = XX(L) - INEL = 0 - HLM = CVMGT( 0.D0, -.5D0*LM(1), INEL.EQ.0) - HLMT = CVMGT( TT, TT+.5D0*LM(L), INEL.EQ.0) - SU1 = PDMAX(1) + PDMAX(1) - SU2 = PDMAX(1)*(1.D0+KK0) - SUR = PDMAX(1)*(1.D0+KK0R) - SU = DMAX1(SUR,DMAX1(SU1,SU2)) - SUT1 = TT + PDMAX(L) + PDMAX(L) - SUT2 = TT + PDMAX(L)*(1.D0+KK0) - SUTR = TT + PDMAX(L)*(1.D0+KK0R) - SUT = DMAX1(SUTR,DMAX1(SUT1,SUT2)) - XC = CVMGT( X0, -SU, X0.GT.0.D0) - RT = TT-RD - - IF(E0.GE.0.D0) GO TO 51 -C + IF(E0.GE.0.D0) GO TO 51 +C C SET CONSTANTS FOR MAXWELLIAN DISTRIBUTION C - TI = -1.D0*E0 - ZARG = DSQRT(TI/(M1*2.D0)) - VELC = SHEATH/M1 -C + TI = -1.D0*E0 + ZARG = DSQRT(TI/(M1*2.D0)) + VELC = SHEATH/M1 +C C NUMBERS FOR VECTORIZED LOOPS C - 51 NUM = MIN( 64, NH) - IH = NUM - IH1 = NUM + 51 NUM = MIN( 64, NH) + IH = NUM + IH1 = NUM C C RANDOM START CONDITIONS -C - IY = IDINT(RI) - IY2 = IDINT(RI2) - IY3 = IDINT(RI3) - ISEED = IY - ISEED2 = IY2 - ISEED3 = IY3 - - IF ( E0.GT.0.D0 ) GO TO 47 - IF ( ALPHA.GE.0.D0 ) THEN +C + IY = IDINT(RI) + IY2 = IDINT(RI2) + IY3 = IDINT(RI3) + ISEED = IY + ISEED2 = IY2 + ISEED3 = IY3 + + 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 49 IV=1,NUM + CALL VELOCV(FG,FFG,E,COSX,COSY,COSZ,SINE,NUM) + DO IV=1,NUM EMX = EMX+E(IV) - 49 CONTINUE - DO iv=1,num - ne = IDINT(DMIN1(5000.D0,e(iv)+1.D0)) - me(ne) = me(ne)+1 - ENDDO -c - GO TO 56 + ENDDO + DO iv=1,num + ne = IDINT(DMIN1(5000.D0,e(iv)+1.D0)) + me(ne) = me(ne)+1 + ENDDO + GO TO 56 C C MAXWELLIAN ENERGY DISTRIBUTION C - ELSE -C - CALL ENERGV(FE,E,COSX,COSY,COSZ,SINE,NUM) - DO 48 IV=1,NUM + ELSE + CALL ENERGV(FE,E,COSX,COSY,COSZ,SINE,NUM) + DO IV=1,NUM EMX = EMX+E(IV) - 48 CONTINUE - GO TO 56 - ENDIF -C - 47 CONTINUE -C - IF(EQUAL(Esig,0.D0)) THEN + ENDDO + GO TO 56 + ENDIF + 47 CONTINUE + IF(EQUAL(Esig,0.D0)) THEN C FIXED PROJECTILE ENERGY - DO IV=1,NUM - E(IV) = E0 - ENDDO - ELSE + DO IV=1,NUM + E(IV) = E0 + ENDDO + ELSE C GAUSSIAN ENERGY DISTRIBUTION - DO IV=1,NUM - 5200 CALL ENERGGAUSS(ISEED2,Esig,Epar,E0) - tryE = tryE+1 - IF(Epar.LE.0.0D0) THEN - negE = negE+1 - GO TO 5200 - ENDIF - E(IV)=Epar - ENDDO - ENDIF + DO IV=1,NUM + 5200 CALL ENERGGAUSS(ISEED2,Esig,Epar,E0) + tryE = tryE+1 + IF(Epar.LE.0.0D0) THEN + negE = negE+1 + GO TO 5200 + ENDIF + E(IV)=Epar + ENDDO + ENDIF C -C die nachfolgende Zeile wurden von Linie 633 hier hin verschoben +C die nachfolgende Zeile wurden von Linie 633 hier hin verschoben C SFE = DMIN1(SB(1),SB(L)) IF ( ALPHA.GE.0.D0 ) THEN @@ -896,134 +896,100 @@ C C C FIXED PROJECTILE ANGLE C -C nachfolgende drei Zeilen waren vorher bei LINIE 633 - ALFA = ALPHA /BW - CALFA = DCOS(ALFA) - SALFA = DSIN(ALFA) - DO IV=1,NUM - COSX(IV) = CALFA - COSY(IV) = SALFA - COSZ(IV) = 0.D0 - SINE(IV) = COSY(IV) -C WRITE(88,*)ALPHA,ALPHASIG,CALFA,SALFA - ENDDO - ELSE +C nachfolgende drei Zeilen waren vorher bei LINIE 633 + ALFA = ALPHA /BW + CALFA = DCOS(ALFA) + SALFA = DSIN(ALFA) + DO IV=1,NUM + COSX(IV) = CALFA + COSY(IV) = SALFA + COSZ(IV) = 0.D0 + SINE(IV) = COSY(IV) + ENDDO + ELSE C C 1-D GAUSSIAN DISTRIBUTION OF ANGLE C - DO IV=1,NUM - CALL ALPHAGAUSS(ISEED3,ALPHASIG,ALPHA,ALFA,ALPHApar, - + CALFA,SALFA,BW) - COSX(IV) = CALFA - COSY(IV) = SALFA - COSZ(IV) = 0.D0 - SINE(IV) = COSY(IV) -C WRITE(88,'(5(F12.5))')ALPHA,ALPHASIG, -C + ALPHApar,CALFA,SALFA - ENDDO - ENDIF -C - ELSEIF (EQUAL(ALPHA,-2.D0)) THEN -C ELSEIF ( ALPHA.EQ.-2. ) THEN + DO IV=1,NUM + CALL ALPHAGAUSS(ISEED3,ALPHASIG,ALPHA,ALFA,ALPHApar, + & CALFA,SALFA,BW) + COSX(IV) = CALFA + COSY(IV) = SALFA + COSZ(IV) = 0.D0 + SINE(IV) = COSY(IV) + ENDDO + ENDIF + ELSEIF (EQUAL(ALPHA,-2.D0)) THEN C C COSINE ANGLE DISTRIBUTION (THREE-DIMENSIONAL) C -CDIR$ IVDEP - DO IV=1,NUM - call ranlux(ran2,2) - RPHI = PI2*DBLE(ran2(1)) - 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) - ENDDO + DO IV=1,NUM + call ranlux(ran2,2) + RPHI = PI2*DBLE(ran2(1)) + 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) + ENDDO - ELSEIF (EQUAL(ALPHA,-1.D0).AND.X0.GT.0.D0) THEN -C ELSEIF ( ALPHA.EQ.-1. AND. X0.GT.0. ) THEN + ELSEIF (EQUAL(ALPHA,-1.D0).AND.X0.GT.0.D0) THEN C C RANDOM DISTRIBUTION C -CDIR$ IVDEP - DO 50 IV=1,NUM -CC RPHI = PI2*RANF() -CC RPHI = PI2*DRAND48() -CC RPHI = PI2*DBLE(RAN(ISEED)) - 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) = 1.D0 -2.D0*RTHETA - SINE(IV) = DSQRT( 1.D0 -COSX(IV)*COSX(IV)) - COSY(IV) = SINE(IV) *DSIN(RPHI) - COSZ(IV) = SINE(IV) *DCOS(RPHI) - 50 CONTINUE -C -C ELSEIF ( ALPHA.EQ.-1. AND. X0.LE.0. ) THEN - ELSEIF (EQUAL(ALPHA,-1.D0).AND.X0.LE.0.D0) THEN -C -CDIR$ IVDEP - DO 55 IV=1,NUM -CC RPHI = PI2*RANF() -CC RPHI = PI2*DRAND48() -CC RPHI = PI2*DBLE(RAN(ISEED)) - 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) = 1.D0 -RTHETA - SINE(IV) = DSQRT( 1.D0 -COSX(IV)*COSX(IV)) - COSY(IV) = SINE(IV) *DSIN(RPHI) - COSZ(IV) = SINE(IV) *DCOS(RPHI) - 55 CONTINUE -C - ENDIF -C - 56 IF ( X0.GT.0.D0 ) GO TO 59 + DO IV=1,NUM + call ranlux(ran2,2) + RPHI = PI2*DBLE(ran2(1)) + RTHETA = DBLE(ran2(2)) + COSX(IV) = 1.D0 -2.D0*RTHETA + SINE(IV) = DSQRT( 1.D0 -COSX(IV)*COSX(IV)) + COSY(IV) = SINE(IV) *DSIN(RPHI) + COSZ(IV) = SINE(IV) *DCOS(RPHI) + ENDDO + ELSEIF (EQUAL(ALPHA,-1.D0).AND.X0.LE.0.D0) THEN + DO IV=1,NUM + call ranlux(ran2,2) + RPHI = PI2*DBLE(ran2(1)) + RTHETA = DBLE(ran2(2)) + COSX(IV) = 1.D0 -RTHETA + SINE(IV) = DSQRT( 1.D0 -COSX(IV)*COSX(IV)) + COSY(IV) = SINE(IV) *DSIN(RPHI) + COSZ(IV) = SINE(IV) *DCOS(RPHI) + ENDDO + ENDIF + 56 IF ( X0.GT.0.D0 ) GO TO 59 C C EXTERNAL START C - DO 57 IV=1,NUM - SINA = SINE(IV) - COSX(IV) = DSQRT( ( E(IV)*COSX(IV)*COSX(IV) +ESB) - & /( E(IV) +ESB)) - SINE(IV) = DSQRT( 1.D0 -COSX(IV)*COSX(IV)) - COSY(IV) = COSY(IV) *SINE(IV) /SINA - COSZ(IV) = COSZ(IV) *SINE(IV) /SINA - E(IV) = E(IV) + ESB - 57 CONTINUE + DO IV=1,NUM + SINA = SINE(IV) + COSX(IV) = DSQRT( ( E(IV)*COSX(IV)*COSX(IV)+ESB)/(E(IV)+ESB)) + SINE(IV) = DSQRT( 1.D0 -COSX(IV)*COSX(IV)) + COSY(IV) = COSY(IV) *SINE(IV) /SINA + COSZ(IV) = COSZ(IV) *SINE(IV) /SINA + E(IV) = E(IV) + ESB + ENDDO C C LOCUS OF FIRST COLLISION C - 59 JL = ISRCHFGT(L,XX(1),1,X0) -C WRITE(*,*)X0 - DO 58 IV=1,NUM -CC RA = CVMGT(RANF(),1.0,X0.LE.0.0) -CC RA = CVMGT(DRAND48(),1.0,X0.LE.0.0) -CC RA = CVMGT(DBLE(RAN(ISEED)),1.D0,X0.LE.0.0D0) - call ranlux(random, 1) - RA = CVMGT(DBLE(random),1.D0,X0.LE.0.0D0) - X(IV) = XC + LM(JL) *RA *COSX(IV) - 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) - 58 CONTINUE -C - DO 199 IV=1,NUM - LLL(IV) = JL - 199 CONTINUE + 59 JL = ISRCHFGT(L,XX(1),1,X0) + DO IV=1,NUM + call ranlux(random, 1) + RA = CVMGT(DBLE(random),1.D0,X0.LE.0.0D0) + X(IV) = XC + LM(JL) *RA *COSX(IV) + 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 C PROJECTILE LOOP C - 1 CONTINUE -C + 1 CONTINUE NPROJ=NPROJ+1 -C - DO 63 IV=1,IH1 + DO IV=1,IH1 CX(IV)=COSX(IV) CY(IV)=COSY(IV) CZ(IV)=COSZ(IV) @@ -1031,7 +997,7 @@ C DEES(IV)=0.D0 DENS(IV)=0.D0 DEN(IV)=0.D0 - 63 CONTINUE + ENDDO KK1=KK0 C C COLLISION LOOP (INCLUDES WEAK SIMULTANEOUS COLL. FOR KK1.LT.4) @@ -1040,45 +1006,31 @@ C C C CHOICE OF COLLISION PARTNERS C - DO 298 IV=1,IH1 + DO IV=1,IH1 call ranlux(random, 1) - JJJ(IV) = ISRCHFGE(NJ(LLL(IV)),COM(1,LLL(IV)),1 -CC # ,RANF())+JT(LLL(IV)) -CC # ,DRAND48())+JT(LLL(IV)) -CC # ,DBLE(RAN(ISEED)))+JT(LLL(IV)) - # ,DBLE(random))+JT(LLL(IV)) - 298 CONTINUE - DO 67 IV=1,IH1 + JJJ(IV) = ISRCHFGE(NJ(LLL(IV)),COM(1,LLL(IV)),1,DBLE(random)) + & +JT(LLL(IV)) + ENDDO + DO IV=1,IH1 EPS(IV)=E(IV)*F1(JJJ(IV)) - 67 CONTINUE -C -CDIR$ IVDEP - DO 64 IV=1,IH1 + ENDDO + DO IV=1,IH1 C C RANDOM AZIMUTHAL ANGLE AND IMPACT PARAMETER C -CC PHIP=PI2*RANF() -CC PHIP=PI2*DRAND48() -CC PHIP=PI2*DBLE(RAN(ISEED)) call ranlux(ran2, 2) PHIP=PI2*DBLE(ran2(1)) CPHI(IV)=DCOS(PHIP) SPHI(IV)=DSIN(PHIP) -CC P(IV)=PDMAX(LLL(IV))*DSQRT(RANF()+KK) -CC P(IV)=PDMAX(LLL(IV))*DSQRT(DRAND48()+KK) -CC P(IV)=PDMAX(LLL(IV))*DSQRT(DBLE(RAN(ISEED))+KK) P(IV)=PDMAX(LLL(IV))*DSQRT(DBLE(ran2(2))+KK) C C POSITION OF TARGET ATOM C X1(IV)=X(IV)-P(IV)*CPHI(IV)*SX(IV) P(IV)=CVMGT(1.D10,P(IV),X1(IV).LT.0.D0.OR.X1(IV).GT.TT) -C IF(A1(JJJ(IV)).EQ.0.) WRITE(99,'(A50)')' A1 vor Label 64 ' B(IV)=P(IV)/A1(JJJ(IV)) - 64 CONTINUE + ENDDO CALL SCOPY(IH1,B,1,R,1) -C WRITE(99,*)IH1,B(IV),R(IV) -C CALL MAGICKRC(C2(1),S2(1),B(1),R(1),EPS(1),IH1) IVMIN=1 IVMAX=IH1 C @@ -1086,35 +1038,35 @@ C MAGIC (DETERMINATION OF SCATTERING ANGLE : KRYPTON-CARBON POT.) C IF(IPOT.NE.1) GO TO 4101 C KRYPTON-CARBON POTENTIAL -C CALL MAGICKRC(C2(1),S2(1),B(1),R(1),EPS(1),IH1) - 104 DO 105 IV=IVMIN,IVMAX + 104 DO IV=IVMIN,IVMAX IF(R(IV).LT.1.D-20)THEN - WRITE(99,'(A70)')'in DO 104 R(IV)<1.D-20 -> 0.00001D0 gesetzt' - R(IV)=0.00001D0 - ENDIF - EX1(IV)=DEXP(-.278544D0*R(IV)) - EX2(IV)=DEXP(-.637174D0*R(IV)) - EX3(IV)=DEXP(-1.919249D0*R(IV)) + WRITE(99,'(A70)' + & )'in DO 104 R(IV)<1.D-20 -> 0.00001D0 gesetzt' + R(IV)=0.00001D0 + ENDIF + EX1(IV)=DEXP(-.278544D0*R(IV)) + EX2(IV)=DEXP(-.637174D0*R(IV)) + EX3(IV)=DEXP(-1.919249D0*R(IV)) - RR1=1.D0/R(IV) -C IF(R(IV).EQ.0.D0)WRITE(99,'(1x,F15.7,1x,A50)')R(IV),' Label 104 ' - V(IV)=(.190945D0*EX1(IV)+.473674D0*EX2(IV)+.335381D0*EX3(IV))*RR1 - FR=B(IV)*B(IV)*RR1+V(IV)*R(IV)/EPS(IV)-R(IV) - V1(IV)=-(V(IV)+.05318658408D0*EX1(IV)+.301812757276D0*EX2(IV)+ - 1 .643679648869D0*EX3(IV))*RR1 - FR1=-B(IV)*B(IV)*RR1*RR1+(V(IV)+V1(IV)*R(IV))/EPS(IV)-1.D0 - Q=FR/FR1 - R(IV)=R(IV)-Q - TEST(IV)=DABS(Q/R(IV)).GT.0.001D0 - 105 CONTINUE + RR1=1.D0/R(IV) + V(IV)=(.190945D0*EX1(IV)+.473674D0*EX2(IV)+.335381D0*EX3(IV)) + & *RR1 + FR=B(IV)*B(IV)*RR1+V(IV)*R(IV)/EPS(IV)-R(IV) + V1(IV)=-(V(IV)+.05318658408D0*EX1(IV)+.301812757276D0*EX2(IV)+ + & .643679648869D0*EX3(IV))*RR1 + FR1=-B(IV)*B(IV)*RR1*RR1+(V(IV)+V1(IV)*R(IV))/EPS(IV)-1.D0 + Q=FR/FR1 + R(IV)=R(IV)-Q + TEST(IV)=DABS(Q/R(IV)).GT.0.001D0 + ENDDO C GET MAX AND MIN INDEX OF TEST FAILURES IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),1) IF(IVMIN.GT.IVMAX) GO TO 106 IVMAX=IVMAX-ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),-1) IF(IVMIN.GT.IVMAX) GO TO 106 GO TO 104 - 106 DO 108 IV=1,IH1 + 106 DO IV=1,IH1 ROCINV=-0.5D0*V1(IV)/(EPS(IV)-V(IV)) SQE=DSQRT(DABS(EPS(IV))) CC=(.235809D0+SQE)/(.126000D0+SQE) @@ -1123,30 +1075,31 @@ C GET MAX AND MIN INDEX OF TEST FAILURES DELTA=(R(IV)-B(IV))*AA*FF/(FF+1.D0) C=(ROCINV*(B(IV)+DELTA)+1.D0)/(ROCINV*R(IV)+1.D0) C2(IV)=DMIN1(1.0D0,C*C) - 108 S2(IV)=1.D0-(1.D0*C2(IV)) + S2(IV)=1.D0-(1.D0*C2(IV)) + ENDDO GO TO 4103 -C + 4101 IF(IPOT.NE.2) GO TO 4102 C MOLIERE POTENTIAL C CALL MAGICMOL(C2(1),S2(1),B(1),R(1),EPS(1),IH1) - 4104 DO 4105 IV=IVMIN,IVMAX + 4104 DO IV=IVMIN,IVMAX IF(R(IV).LT.1.D-20)THEN - WRITE(99,'(A70)')'in DO 4104 R(IV)<1.D-20 -> 0.00001D0 gesetzt' - R(IV)=0.00001D0 - ENDIF - EX1(IV)=DEXP(-.3D0*R(IV)) - EX2(IV)=DEXP(-1.2D0*R(IV)) - EX3(IV)=DEXP(-6.0D0*R(IV)) -C IF(R(IV).EQ.0.D0)WRITE(99,'(A50)')' R nach Label 4104 ' - RR1=1.D0/R(IV) - V(IV)=(.35D0*EX1(IV)+.55D0*EX2(IV)+.10D0*EX3(IV))*RR1 - FR=B(IV)*B(IV)*RR1+V(IV)*R(IV)/EPS(IV)-R(IV) - V1(IV)=-(V(IV)+.105D0*EX1(IV)+.66D0*EX2(IV)+.6D0*EX3(IV))*RR1 - FR1=-B(IV)*B(IV)*RR1*RR1+(V(IV)+V1(IV)*R(IV))/EPS(IV)-1.D0 - Q=FR/FR1 - R(IV)=R(IV)-Q - TEST(IV)=DABS(Q/R(IV)).GT.0.001D0 - 4105 CONTINUE + WRITE(99,'(A70)' + & )'in DO 4104 R(IV)<1.D-20 -> 0.00001D0 gesetzt' + R(IV)=0.00001D0 + ENDIF + EX1(IV)=DEXP(-.3D0*R(IV)) + EX2(IV)=DEXP(-1.2D0*R(IV)) + EX3(IV)=DEXP(-6.0D0*R(IV)) + RR1=1.D0/R(IV) + V(IV)=(.35D0*EX1(IV)+.55D0*EX2(IV)+.10D0*EX3(IV))*RR1 + FR=B(IV)*B(IV)*RR1+V(IV)*R(IV)/EPS(IV)-R(IV) + V1(IV)=-(V(IV)+.105D0*EX1(IV)+.66D0*EX2(IV)+.6D0*EX3(IV))*RR1 + FR1=-B(IV)*B(IV)*RR1*RR1+(V(IV)+V1(IV)*R(IV))/EPS(IV)-1.D0 + Q=FR/FR1 + R(IV)=R(IV)-Q + TEST(IV)=DABS(Q/R(IV)).GT.0.001D0 + ENDDO C GET MAX AND MIN INDEX OF TEST FAILURES IVMIN=IVMIN+ILLZ(IVMAX-IVMIN+1,TEST(IVMIN),1) IF(IVMIN.GT.IVMAX) GO TO 4106