diff --git a/fortran/trimspNL.F b/fortran/trimspNL.F index 428199e..6126c7f 100644 --- a/fortran/trimspNL.F +++ b/fortran/trimspNL.F @@ -387,97 +387,52 @@ C part. refl. coeff. from Thomas et al. 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 +C elements per layer. JMAX=5 - if (OldNew(innam).eq.1) then C This part reads the input file (new format) - OPEN(UNIT=99,file=errnam,STATUS='replace') - OPEN(UNIT=11,file=innam,STATUS='unknown',ERR=1359) + OPEN(UNIT=99,file=errnam,STATUS='replace') + OPEN(UNIT=11,file=innam,STATUS='unknown',ERR=1359) C First line: properties of projectile C Ordered as: Z, mass number (amu), imp. energy, dist. of imp. energy, angle, dist. angles, - READ(11,*) Z1,M1,E0,Esig,ALPHA,ALPHASIG,EF,ESB,SHEATH,ERC + READ(11,*) Z1,M1,E0,Esig,ALPHA,ALPHASIG,EF,ESB,SHEATH,ERC C Second line: simulation related parameters C Ordered as: Number of particles, seed, seed, seed, initial depth, RD, depth increment, CA, KK0, KDEE1,KDEE2,IPOT - READ(11,*) NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,KK0R,KDEE1,KDEE2,IPOT - & ,IPOTR,IRL + READ(11,*) NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,KK0R,KDEE1,KDEE2,IPOT + & ,IPOTR,IRL C Third line: Number of layers - read(11,66) NLayers - 66 format(9x,I4) - LMAX=NLayers + read(11,66) NLayers + 66 format(9x,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 -C always 1.0??), chemical formula (chem) - READ(11,*) chem(I),DX(I),RHO(I),CK(I) + DO I=1,NLayers +C Chemical formula (chem), Thickness (DX), density (RHO), and correction factor +C (CK, it is always 1.0??) +C Possibly add the number of elements in the layer + READ(11,*) chem(I),DX(I),RHO(I),CK(I) C Atomic numbers - 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) C Mass numbers (amu) - 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) C Concentration (stoichiometry) - 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) C Surface binding energy - READ(11,*) SBE(I,1),SBE(I,2),SBE(I,3),SBE(I,4),SBE(I,5) -C Displacement energy - READ(11,*) ED(I,1),ED(I,2),ED(I,3),ED(I,4),ED(I,5) -C Bulk binding energy - READ(11,*) BE(I,1),BE(I,2),BE(I,3),BE(I,4),BE(I,5) + READ(11,*) SBE(I,1),SBE(I,2),SBE(I,3),SBE(I,4),SBE(I,5) +C Displacement energy, always 30 eV?? + READ(11,*) ED(I,1),ED(I,2),ED(I,3),ED(I,4),ED(I,5) +C Bulk binding energy, usually zero + READ(11,*) BE(I,1),BE(I,2),BE(I,3),BE(I,4),BE(I,5) C value A-1 of the ziegler tables - READ(11,*) CH1(I,1),CH1(I,2),CH1(I,3),CH1(I,4),CH1(I,5) + READ(11,*) CH1(I,1),CH1(I,2),CH1(I,3),CH1(I,4),CH1(I,5) C value A-2 of the ziegler tables - READ(11,*) CH2(I,1),CH2(I,2),CH2(I,3),CH2(I,4),CH2(I,5) + READ(11,*) CH2(I,1),CH2(I,2),CH2(I,3),CH2(I,4),CH2(I,5) C value A-3 of the ziegler tables - 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) C value A-4 of the ziegler tables - 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) 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 -C This part reads the input file (old format, 7 layers) -C To be phased out soon (aim for beginning of 2023). - OPEN(UNIT=99,file=errnam,STATUS='replace') - OPEN(UNIT=11,file=innam,STATUS='unknown',ERR=1359) - -C First line: properties of projectile - READ(11,*) Z1,M1,E0,Esig,ALPHA,ALPHASIG,EF,ESB,SHEATH,ERC -C Second line: simulation related parameters - READ(11,*) NH,RI,RI2,RI3,X0,RD,CW,CA,KK0,KK0R,KDEE1,KDEE2,IPOT - & ,IPOTR,IRL -C Third line: layer structure. To be replaced by number of layers -C and then each layer with its properties: Thickness (DX), density -C (RHO), and correction factor (CK, it is always 1.0??) - 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) -C Here we read the 7 layer structure - DO I=1,7 -C Atomic numbers - READ(11,*) ZT(I,1),ZT(I,2),ZT(I,3),ZT(I,4),ZT(I,5) -C Mass numbers (amu) - READ(11,*) MT(I,1),MT(I,2),MT(I,3),MT(I,4),MT(I,5) -C Concentration - READ(11,*) CO(I,1),CO(I,2),CO(I,3),CO(I,4),CO(I,5) -C Surface binding energy - READ(11,*) SBE(I,1),SBE(I,2),SBE(I,3),SBE(I,4),SBE(I,5) -C Displacement energy - READ(11,*) ED(I,1),ED(I,2),ED(I,3),ED(I,4),ED(I,5) -C Bulk binding energy - READ(11,*) BE(I,1),BE(I,2),BE(I,3),BE(I,4),BE(I,5) -C value A-1 of the ziegler tables - READ(11,*) CH1(I,1),CH1(I,2),CH1(I,3),CH1(I,4),CH1(I,5) -C value A-2 of the ziegler tables - READ(11,*) CH2(I,1),CH2(I,2),CH2(I,3),CH2(I,4),CH2(I,5) -C value A-3 of the ziegler tables - READ(11,*) CH3(I,1),CH3(I,2),CH3(I,3),CH3(I,4),CH3(I,5) -C value A-4 of the ziegler tables - READ(11,*) CH4(I,1),CH4(I,2),CH4(I,3),CH4(I,4),CH4(I,5) -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 - endif + READ(11,*) CH5(I,1),CH5(I,2),CH5(I,3),CH5(I,4),CH5(I,5) + ENDDO 1359 CLOSE(UNIT=11) @@ -548,11 +503,13 @@ C If the thickness of layer K is not an integer multiples of depth increment 44 CONTINUE DO I=1,L +C This loop finds out how many elements are in the layer. +C I can be removed in the number is know from the input file. DO J=1,JMAX 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 +C I am guessing NJ(I) is the number of elements in layer I 156 NJ(I)=J-1 ENDDO JT(1) = 0 @@ -4397,37 +4354,6 @@ C in seconds from beginning of year END - INTEGER FUNCTION OldNew(filename) -C This funnction returns 0 for old input format and 1 for new input -C format - CHARACTER filename*12 - REAL dummy(21) - INTEGER idummy -C Assume old format - OldNew = 0 - -C Read first three lines, the third line differs between the two -C format - OPEN(UNIT=11,file=filename,STATUS='unknown') - READ(11,*) - READ(11,*) - READ(11,*,ERR=10) dummy(1),dummy(2),dummy(3),dummy(4),dummy(5) - & ,dummy(6),dummy(7),dummy(8),dummy(9),dummy(10),dummy(11) - & ,dummy(12),dummy(13),dummy(14),dummy(15),dummy(16),dummy(17) - & ,dummy(18),dummy(19),dummy(20),dummy(21) - GOTO 20 -C If you got here it means that the file in in the old format - 10 OldNew = 1 - 20 CLOSE(11) -C if (OldNew.eq.1) then -C write(*,*) "The input file is in the new format" -C else -C write(*,*) "The input file is in the old format" -C endif - - RETURN - END - SUBROUTINE RANLUX(RVEC,LENV) C Subtract-and-borrow random number generator proposed by C Marsaglia and Zaman, implemented by F. James with the name