Cleanup code in progress.

This commit is contained in:
salman 2013-02-21 13:06:27 +00:00
parent 6b949a2841
commit a1b7d7365a

View File

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