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

The new file format gives same results for a structure up to 7 layers.

Browse Source
This commit is contained in:
salman 2013-03-28 12:50:32 +00:00
parent 280b770a6c
commit 1cdc02a760
1 changed files with 135 additions and 86 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

221
trimsp/src/trimspNL.F
View File

@ -184,12 +184,15 @@ c Included the source of ranlux for random numbers generation
c into trimsp7l source. No need for cern libraries to be installed.
c
c Feb 2013: Zaher Salman PSI
c Started cleaning up the code.
c Cleaned up the code.
c When possible remove loop numbres and use do-enddo instead.
c Using proper fortran line indentation (emacs style).
c Created TimeStamp subroutine to generate time stamp and removed
c original code from main (two places)
c
c Mar 2013: Zaher Salman PSI
c Started implementation of a more flexible input file format, with
c flexible number of layer, not limited to max 7.
c-------------------------------------------
c check OS
c
@ -200,6 +203,13 @@ c
#endif
IMPLICIT NONE
C These parameters are related to the maximum number of layers
INTEGER MAXNL,MAXNL5,MAXNLp25,MAXNL5p2,MAXNLm15
PARAMETER (MAXNL=100)
PARAMETER (MAXNL5=MAXNL*5)
PARAMETER (MAXNLp25=MAXNL*MAXNL5)
PARAMETER (MAXNL5p2=MAXNL5*MAXNL5*100)
PARAMETER (MAXNLm15=(MAXNL-1)*5)
LOGICAL TEST(64),TESTR(2000),TEST1(2000)
LOGICAL EQUAL
INTEGER*4 ISRCHEQ,ISRCHFGT,ISRCHFGE,ILLZ
@ -214,29 +224,31 @@ c
INTEGER*4 ISEED,ISEED2,ISEED3
INTEGER*4 JJR(2000,2),INOUT(2000,2),LRR(2000,2)
INTEGER*4 IDMAX(2000),IKR(2000)
INTEGER*4 number_in_layer(7),laufzahl
INTEGER*4 number_in_layer(MAXNL),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 ITSP(35),ISPALT(7)
& ,ISPIPT(35),ISPIST(35),ISPOPT(35),ISPOST(35)
INTEGER*4 KO(600,35,2)
INTEGER*4 ICD(100,MAXNL5),ICDT(100),ICDJT(MAXNL5) ,ICDIRJ(MAXNL5
& ,MAXNL5),ICDR(100,MAXNL5),ICDTR(100),ICDJTR(MAXNL5)
& ,ICDIRI(100,MAXNL5,MAXNL5) ,ICDIRN(100,MAXNL5),ICDITR(MAXNL5)
INTEGER*4 KADB(20),KADT(20),KADS(20),KADST(20) ,KADRIP(20,MAXNLm15
& ),KADRIS(20,MAXNLm15),KADROP(20,MAXNLm15),KADROS(20,MAXNLm15)
& ,KADSJ(20,MAXNLm15),KADSL(20,6),KDSTJ(20,MAXNLm15),KDSTL(20,6
& )
INTEGER*4 IBSP(MAXNL5),ISPAL(MAXNL),IBSPL(MAXNL5),ISPIP(MAXNL5)
& ,ISPIS(MAXNL5),ISPOP(MAXNL5),ISPOS(MAXNL5)
INTEGER*4 ITSP(MAXNL5),ISPALT(MAXNL)
& ,ISPIPT(MAXNL5),ISPIST(MAXNL5),ISPOPT(MAXNL5),ISPOST(MAXNL5)
INTEGER*4 KO(600,MAXNL5,2)
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)
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)
INTEGER*4 NJ(7),JT(7),ILD(7)
INTEGER*4 MEAS(102,22,MAXNLm15),MAGS(62,22,MAXNLm15),MAGSA(62,32
& ,MAXNLm15),MEAGS(102,12,22,MAXNLm15) ,MEASL(75,21,MAXNLm15)
INTEGER*4 MEAST(102,22,MAXNLm15),MAGST(62,22,MAXNLm15),MEASTL(75
& ,21,MAXNLm15)
INTEGER*4 NJ(MAXNL),JT(MAXNL),ILD(MAXNL)
INTEGER*4 LLL(64),JJJ(64),IK(64)
INTEGER*4 me(5000),nli(0:7),irpl(7)
INTEGER*4 me(5000),nli(MAXNL),irpl(MAXNL)
INTEGER*4 IT,NPROJ
INTEGER*4 IB,IBL
INTEGER*4 IIRP,IIPL,ICDTT,ICDTTR
@ -274,42 +286,58 @@ C REAL Variables
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 ELE(100,MAXNL5),ELI(100,MAXNL5),ELP(100,MAXNL5),ELD(100
& ,MAXNL5) ,ELET(MAXNL5),ELIT(MAXNL5),ELPT(MAXNL5),ELDT(MAXNL5)
& ,ELER(100,MAXNL5),ELIR(100,MAXNL5),ELPR(100,MAXNL5),ELDR(100
& ,MAXNL5) ,ELETR(MAXNL5),ELITR(MAXNL5),ELPTR(MAXNL5)
& ,ELDTR(MAXNL5)
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)
& ,RKADSJ(20,MAXNLm15),RKADSL(20,MAXNL),RKDSTJ(20,MAXNLm15)
& ,RKDSTL(20,MAXNL)
REAL*8 EBSP(MAXNL5),ESPAL(MAXNL) ,SPY(MAXNL5),SPE(MAXNL5)
& ,REY(MAXNL5),EMSP(MAXNL5),ESPIP(MAXNL5),ESPIS(MAXNL5)
& ,ESPOP(MAXNL5),ESPOS(MAXNL5) ,RIP(MAXNL5),RIS(MAXNL5)
& ,ROP(MAXNL5),ROS(MAXNL5) ,REIP(MAXNL5),REIS(MAXNL5)
& ,REOP(MAXNL5),REOS(MAXNL5),ESPMIP(MAXNL5),ESPMIS(MAXNL5)
& ,ESPMOP(MAXNL5),ESPMOS(MAXNL5) ,RIPJ(MAXNL5),RISJ(MAXNL5)
& ,ROPJ(MAXNL5),ROSJ(MAXNL5) ,REIPJ(MAXNL5),REISJ(MAXNL5)
& ,REOPJ(MAXNL5),REOSJ(MAXNL5)
REAL*8 ETSP(MAXNL5),ESPALT(MAXNL) ,SPYT(MAXNL5),SPET(MAXNL5)
& ,REYT(MAXNL5),EMSPT(MAXNL5),ESPIPT(MAXNL5),ESPIST(MAXNL5)
& ,ESPOPT(MAXNL5),ESPOST(MAXNL5),RIPT(MAXNL5),RIST(MAXNL5)
& ,ROPT(MAXNL5),ROST(MAXNL5) ,REIPT(MAXNL5),REIST(MAXNL5)
& ,REOPT(MAXNL5),REOST(MAXNL5) ,ESPMIPT(MAXNL5),ESPMIST(MAXNL5)
& ,ESPMOPT(MAXNL5),ESPMOST(MAXNL5)
REAL*8 SPEM(MAXNL5),SPE2S(MAXNL5),SPE3S(MAXNL5),SPE4S(MAXNL5)
& ,SPE5S(MAXNL5) ,SPE6S(MAXNL5),VSPE(MAXNL5),SSPE(MAXNL5)
& ,GSPE(MAXNL5),BSPE(MAXNL5)
REAL*8 SPE1SL(MAXNL5),SPE2SL(MAXNL5),SPE3SL(MAXNL5),SPE4SL(MAXNL5)
& ,SPE5SL(MAXNL5),SPE6SL(MAXNL5)
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,MAXNLm15),EASTL(75,MAXNLm15)
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(MAXNL),DX(MAXNL),RHO(MAXNL)
REAL*8 Z2(MAXNL),M2(MAXNL),LM(MAXNL)
REAL*8 PDMAX(MAXNL),ARHO(MAXNL),AM(MAXNL)
REAL*8 FM(MAXNL),EPS0(MAXNL),ASIG(MAXNL)
REAL*8 K2(MAXNL),CK(MAXNL),KLM1(MAXNL)
REAL*8 SB(MAXNL),DLI(MAXNL)
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 CH1(7,5),CH2(7,5),CH3(7,5),CH4(7,5),CH5(7,5)
REAL*8 CHM1(7)
REAL*8 ZT(MAXNL,5),MT(MAXNL,5),CO(MAXNL,5)
& ,SBE(MAXNL,5),ED(MAXNL,5),BE(MAXNL,5),
& COM(5,MAXNL)
REAL*8 MU(MAXNL5,MAXNL5),EC(MAXNL5,MAXNL5),A(MAXNL5,MAXNL5)
& ,F(MAXNL5,MAXNL5) ,KL(MAXNL5,MAXNL5),KOR(MAXNL5,MAXNL5)
& ,KLM(MAXNL,MAXNL5)
REAL*8 MU1(MAXNL5),EC1(MAXNL5),A1(MAXNL5),F1(MAXNL5),KL1(MAXNL5)
& ,KOR1(MAXNL5) ,DI(MAXNL5),EP(MAXNL5),ZZ(MAXNL5),TM(MAXNL5)
REAL*8 CH1(MAXNL,5),CH2(MAXNL,5),CH3(MAXNL,5)
& ,CH4(MAXNL,5),CH5(MAXNL,5)
REAL*8 CHM1(MAXNL)
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(MAXNLm15),SEESB(MAXNLm15),THESB(MAXNLm15)
& ,FOESB(MAXNLm15) ,SGMESB(MAXNLm15),DFIESB(MAXNLm15)
& ,DSEESB(MAXNLm15),DTHESB(MAXNLm15)
REAL*8 pi,c,E0,de,alfa,z1,rtheta,cu,enot,esb,est,esp
REAL*8 E2,AB,FP,AN
REAL*8 Esig,Epar
@ -411,7 +439,8 @@ C CHARACTER Variables
DATA ET/0.D0/,PLST/0.D0/,PL2ST/0.D0/,PL3ST/0.D0/
DATA PL4ST/0.D0/,PL5ST/0.D0/,PL6ST/0.D0/
DATA SEM/0.D0/,IT/0/,NPROJ/0/,NREC1/0/,NREC2/0/
DATA NH/0/,IB/0/,IBL/0/,NJ/7*0/,NLI/8*0/,DLI/7*0.D0/
DATA NH/0/,IB/0/,IBL/0/,NJ/MAXNL*0/,NLI/MAXNL*0/
DATA DLI/MAXNL*0.D0/
DATA tryE/0/,negE/0/
DATA IIRP/0/,IIPL/0/,ICDTT/0/,ICDTTR/0/,TDMENR/0.D0/
DATA TION/0.D0/,TIONR/0.D0/,TDENT/0.D0/,TDENTR/0.D0/
@ -456,40 +485,50 @@ C CHARACTER Variables
DATA PHONR/100*0.D0/
DATA ELGD/100*0.D0/,ELGDR/100*0.D0/
DATA ICDT/100*0/,ICDTR/100*0/
DATA ICDR/3500*0/,ICDIRN/3500*0/,IONR/100*0.D0/
DATA ICDR/MAXNLp25*0/,ICDIRN/MAXNLp25*0/,IONR/100*0.D0/
DATA DENTR/100*0.D0/,DMGNR/100*0.D0/
DATA IPL/100*0/,IPLB/100*0/,IPLT/100*0/
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/
DATA ELE/3500*0.D0/,ELI/3500*0.D0/
DATA ICDIRI/122500*0/
DATA IRPL/MAXNL*0/
DATA ICDJT/MAXNL5*0/,ICDJTR/MAXNL5*0/,ICDITR/MAXNL5*0/
DATA ICD/MAXNLp25*0/,ELP/MAXNLp25*0.D0/,ELD/MAXNLp25*0.D0/
DATA ELE/MAXNLp25*0.D0/,ELI/MAXNLp25*0.D0/
DATA ICDIRI/MAXNL5p2*0/
DATA ICSUM/0/,ICSUMS/0/,ICDI/0/,ISPA/0/,ISPAT/0/
DATA Z2/7*0.D0/,M2/7*0.D0/
DATA KLM1/7*0.D0/,CHM1/7*0.D0/,SB/7*0.D0/,KLM/245*0.D0/
DATA Z2/MAXNL*0.D0/,M2/MAXNL*0.D0/
DATA KLM1/MAXNL*0.D0/,CHM1/MAXNL*0.D0/
DATA SB/MAXNL*0.D0/,KLM/MAXNLp25*0.D0/
DATA ME/5000*0/,EMX/0.D0/,ESPAT/0.D0/,ESPA/0.D0/
DATA IBSP/35*0/,IBSPL/35*0/,EBSP/35*0.D0/,ISPAL/7*0/
DATA ITSP/35*0/,ETSP/35*0.D0/
DATA ESPAL/7*0.D0/,ESPALT/7*0.D0/,ISPALT/7*0/
DATA SPE2S/35*0.D0/,SPE3S/35*0.D0/,SPE4S/35*0.D0/
DATA SPE5S/35*0.D0/,SPE6S/35*0.D0/
DATA SPE1SL/35*0.D0/,SPE2SL/35*0.D0/,SPE3SL/35*0.D0/
DATA SPE4SL/35*0.D0/,SPE5SL/35*0.D0/,SPE6SL/35*0.D0/
DATA ELET/35*0.D0/,ELPT/35*0.D0/,ELDT/35*0.D0/,ELIT/35*0.D0/
DATA ELETR/35*0.D0/,ELITR/35*0.D0/,ELPTR/35*0.D0/
DATA ELDTR/35*0.D0/
DATA IBSP/MAXNL5*0/,IBSPL/MAXNL5*0/,EBSP/MAXNL5*0.D0/
DATA ISPAL/MAXNL*0/
DATA ITSP/MAXNL5*0/,ETSP/MAXNL5*0.D0/
DATA ESPAL/MAXNL*0.D0/,ESPALT/MAXNL*0.D0/
DATA ISPALT/MAXNL*0/
DATA SPE2S/MAXNL5*0.D0/,SPE3S/MAXNL5*0.D0/,SPE4S/MAXNL5*0.D0/
DATA SPE5S/MAXNL5*0.D0/,SPE6S/MAXNL5*0.D0/
DATA SPE1SL/MAXNL5*0.D0/,SPE2SL/MAXNL5*0.D0/,SPE3SL/MAXNL5*0.D0/
DATA SPE4SL/MAXNL5*0.D0/,SPE5SL/MAXNL5*0.D0/,SPE6SL/MAXNL5*0.D0/
DATA ELET/MAXNL5*0.D0/,ELPT/MAXNL5*0.D0/,ELDT/MAXNL5*0.D0/
DATA ELIT/MAXNL5*0.D0/
DATA ELETR/MAXNL5*0.D0/,ELITR/MAXNL5*0.D0/,ELPTR/MAXNL5*0.D0/
DATA ELDTR/MAXNL5*0.D0/
DATA Epar/0.D0/
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/
DATA number_in_layer /MAXNL*0/
innam=filein//inext
outnam=fileout//outext
rgenam=fileout//rgeext
errnam=fileout//errext
C LMAX is maximum number of layers and JMAX is maximum number of
C elements per layer. Assume these values for old files and change
C LMAX as needed for the new format.
LMAX=7
JMAX=5
if (OldNew(innam).eq.1) then
OPEN(UNIT=99,file=errnam)
OPEN(UNIT=99,file=errnam,STATUS='replace')
OPEN(UNIT=11,file=innam,STATUS='unknown',ERR=1359)
C This part reads the input file (new format)
@ -501,6 +540,7 @@ C Second line: simulation related parameters
C Third line: Number of layers
read(11,66) NLayers
66 format(8x,I4)
LMAX=NLayers
C Here we read the NLayers structure
DO I=1,NLayers
C Thickness (DX), density (RHO), and correction factor (CK, it is
@ -530,7 +570,7 @@ C value A-5 of the ziegler tables
READ(11,*) CH5(I,1),CH5(I,2),CH5(I,3),CH5(I,4),CH5(I,5)
ENDDO
else
OPEN(UNIT=99,file=errnam)
OPEN(UNIT=99,file=errnam,STATUS='replace')
OPEN(UNIT=11,file=innam,STATUS='unknown',ERR=1359)
C This part reads the input file (old format, 7 layers)
@ -576,25 +616,25 @@ C value A-5 of the ziegler tables
C open statement for output files, removed from line 2449 ff to here
OPEN(UNIT=21,FILE=outnam,STATUS='new',ERR=6000)
OPEN(UNIT=21,FILE=outnam)
GOTO 6001
6000 WRITE(*,*)' File schon vorhanden, Gib neue Ausgabedatei an (A8)'
READ(*,'(A8)') fileout
outnam=fileout//outext
rgenam=fileout//rgeext
OPEN(UNIT=21,FILE=outnam,STATUS='new',ERR=6000)
OPEN(UNIT=21,FILE=outnam,STATUS='new')
6001 OPEN(UNIT=22,FILE=rgenam,STATUS='new')
WRITE(21,1000)
1000 FORMAT(1H1/,6X,'* PROGRAM TRVMC95 - V TrimSP7L 17.Oct.02 TP *')
C Get simulation start time
CALL TimeStamp(day_start,month_start,year_start,hour_start
& ,min_start,sec_start,seconds_start_total)
WRITE(21,*)
WRITE(21,10050)day_start,month_start,year_start,hour_start
WRITE(21,1050)day_start,month_start,year_start,hour_start
& ,min_start,sec_start
WRITE(*,10050)day_start,month_start,year_start,hour_start
WRITE(*,1050)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
1050 FORMAT(1x,' TrimSp simulation started at: ',A2,'.',A4,1x,A4,1x,A2
& ,':',A2,':',A2)
C SET INTERVAL CONSTANTS FOR OUTPUT
@ -622,13 +662,13 @@ C SET INTERVAL CONSTANTS FOR OUTPUT
C CALCULATION OF CHARGE AND MASS DEPENDENT CONSTANTS
PI2=2.D0*PI
ABC=AB*FP
LMAX=7
JMAX=5
C This function is used only once to get the number of defined
C layers. It should be replaced!
L=ISRCHEQ(LMAX,DX(1),1,0.D0)-1
C Checks wether depth interval is an integer denominator of layer thickness or not
C If not, calculated implantation profile is not correct.
depth_interval_flag = 1
DO K=1,L-1
IF(.NOT.EQUAL(DX(K)/CW-DBLE(IDINT(DX(K)/CW)),0.D0)) THEN
@ -643,16 +683,25 @@ C If not, calculated implantation profile is not correct.
IF(EQUAL(CO(I,J),0.D0)) GOTO 156
ENDDO
J=JMAX+1
C I am guessing NJ(I) is the number of elements in layer I
156 NJ(I)=J-1
ENDDO
JT(1) = 0
JT(2) = NJ(1)
JT(3) = NJ(1)+NJ(2)
JT(4) = JT(3)+NJ(3)
JT(5) = JT(4)+NJ(4)
JT(6) = JT(5)+NJ(5)
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)
do I=3,L
JT(I)=JT(I-1)+NJ(I-1)
LJ=LJ+NJ(I)
enddo
C JT(3) = JT(2)+NJ(2)
C JT(4) = JT(3)+NJ(3)
C JT(5) = JT(4)+NJ(4)
C JT(6) = JT(5)+NJ(5)
C JT(7) = JT(6)+NJ(6)
C LJ = NJ(1)+NJ(2)+NJ(3)+NJ(4)+NJ(5)+NJ(6)+NJ(7)
XX(1)=DX(1)
DO I=2,L
XX(I)=XX(I-1)+DX(I)