From aaeb7f969a8199936004a68c6dcf1b597077fafd Mon Sep 17 00:00:00 2001 From: salman Date: Thu, 18 Jun 2026 10:55:37 +0200 Subject: [PATCH] Extract the implanted depth/profile table writer --- fortran/trimspNL.F | 376 +++++++++++++++++++++++++-------------------- 1 file changed, 210 insertions(+), 166 deletions(-) diff --git a/fortran/trimspNL.F b/fortran/trimspNL.F index 348470b..1bc3684 100644 --- a/fortran/trimspNL.F +++ b/fortran/trimspNL.F @@ -2572,172 +2572,13 @@ C & ,76X,'MEAN NR OF EL.COLL.(E > EDISPL.) 2-1 = ',F10.3/ & ,76X,'MEAN NR OF EL.COLL.(E > EDISPL.) 2-2 = ',F10.3) 1453 CONTINUE - 590 WRITE(21,600) - 600 FORMAT(1H1,///,5X,8HDEPTH(A),2X,9HPARTICLES,2X,10HNORM.DEPTH,1X, - & 10HPATHLENGTH,3X,10HINLOSS(EV),2X,10HTELOSS(EV),2X - & ,10HELLOSS(EV), 2X,10HDAMAGE(EV),2X,10HPHONON(EV),2X - & ,10HCASCAD(EV),5X,3HDPA/) - - IF(depth_interval_flag.EQ.0) THEN -C WRITE(22,*)' CALCULATED IMPLANTATION PROFILE DID NOT ', -C & 'AGREE WITH LAYER THICKNESS' - WRITE(21,*)' CALCULATED IMPLANTATION PROFILE DID NOT ', - & 'AGREE WITH LAYER THICKNESS' - ENDIF - -#if defined (OS_WIN) - WRITE(22,6002) - 6002 FORMAT(1H1,///,5X,8HDEPTH(A),2X,9HPARTICLES) -#else - write(22,'(a)') ' DEPTH PARTICLES' -#endif - IF(YH.LT.1.D0) GO TO 603 - DO I=0,MAXD1 -C This normalization of the implantation profile is wrong if the -C simulation does not account for all implanted projectiles, -C e.g. when the number of points in the profile is not enough to -C reach the maximum implantation depth. - RIRP(I) = DBLE(IRP(I))/YH -C To resolve this issue we need to check whether we reached maximum -C implantation or not. If not rerun calculation with larger step -C size! - ENDDO - 603 D1=0. - D2=CW - WRITE(21,601) D1,IRP(0),RIRP(0) - 601 FORMAT(4X,3H-SU,1H-,F6.0,I10,E12.4) - DO J=1,LJ - DO I=1,MAXD - ICDT(I)=ICDT(I)+ICD(I,J) - ICDTR(I)=ICDTR(I)+ICDR(I,J) - ENDDO - ENDDO - DO K=1,NJ(1) - DO J=1,NJ(1) - DO I=1,MAXD - ICDIRN(I,J)=ICDIRN(I,J)+ICDIRI(I,K,J) - ENDDO - ENDDO - ENDDO - DO I=0,MAXD1 - IIRP=IIRP+IRP(I) - TRIRP=TRIRP+RIRP(I) - ENDDO - DO I=1,MAXD - IIPL=IIPL+IPL(I) - TION=TION+ION(I) - TDENT=TDENT+DENT(I) - TDMGN=TDMGN+DMGN(I) - TELGD=TELGD+ELGD(I) - TPHON=TPHON+PHON(I) - TCASMO=TCASMO+CASMOT(I) - ICDTT=ICDTT+ICDT(I) - TIONR=TIONR+IONR(I) - TDENTR=TDENTR+DENTR(I) - TELGDR=TELGDR+ELGDR(I) - TDMGNR=TDMGNR+DMGNR(I) - TPHONR=TPHONR+PHONR(I) - ICDTTR=ICDTTR+ICDTR(I) - ENDDO - do im1=MAXD,1,-1 - if(ipl(im1).ne.0.or.(.NOT.EQUAL(ion(im1),0.D0))) goto 20 - enddo - im1=1 - 20 im1=min0(im1+2,MAXD) - DO I=1,im1 - WRITE(21,700) D1,D2,IRP(I),RIRP(I),IPL(I),ION(I),DENT(I), - & DMGN(I),ELGD(I),PHON(I),CASMOT(I),ICDT(I) - Dmid=(D2-D1)/2+D1 - WRITE(22,701) Dmid,IRP(I) - 700 FORMAT(1X,F6.0,1H-,F6.0,I10,E12.4,I10,1P1E14.4,5E12.4,I8) - 701 FORMAT(1X, F16.4, 2X, I10) - D1=D2 - D2=D2+CW - ENDDO - WRITE(21,604) D2-CW,IRP(MAXD1),RIRP(MAXD1) - 604 FORMAT(1X,F6.0,1H-,3X,3HSUT,I10,E12.4) - WRITE(21,710) IIRP,TRIRP,IIPL,TION,TDENT,TDMGN,TELGD,TPHON,TCASMO - & ,ICDTT - 710 FORMAT(/14X,I10,1P1E12.4,I10,1E14.4,5E12.4,I8) - DO J=1,NJ(1) - DO I=1,MAXD - ELET(J)=ELET(J)+ELE(I,J) - ELIT(J)=ELIT(J)+ELI(I,J) - ELPT(J)=ELPT(J)+ELP(I,J) - ELDT(J)=ELDT(J)+ELD(I,J) - ICDJT(J)=ICDJT(J)+ICD(I,J) - ICDJTR(J)=ICDJTR(J)+ICDR(I,J) - ICDITR(J)=ICDITR(J)+ICDIRN(I,J) - ENDDO - ENDDO - - WRITE(21,1521) - 1521 FORMAT(1H1,4X,'DEPTH(A)' - & ,3X,' INLOSS(1)',3X,'ELLOSS(1)',3X,'DAMAGE(1)',3X,'PHONON(1)' - & ,2X,' INLOSS(2)',3X,'ELLOSS(2)',3X,'DAMAGE(2)',3X,'PHONON(2)' - & ,2X,'DPA(1)',2X,'DPA(2)'/) - D1=0. - D2=CW - do im2=MAXD,1,-1 - if(.NOT.EQUAL(eli(im2,1),0.D0).or.(.NOT.EQUAL(eli(im2,2),0.D0)) - & ) goto 30 - - enddo - im2=1 - 30 im2=MIN0(im2+2,MAXD) - DO 1525 I=1,im2 - WRITE(21,1523) D1,D2,ELI(I,1),ELE(I,1),ELD(I,1),ELP(I,1), - & ELI(I,2),ELE(I,2),ELD(I,2),ELP(I,2),ICD(I,1),ICD(I,2) - 1523 FORMAT(1X,F6.0,1H-,F6.0,1P8E12.4,2I8) - D1=D2 - D2=D2+CW - 1525 CONTINUE - WRITE(21,1533) ELIT(1),ELET(1),ELDT(1),ELPT(1),ELIT(2),ELET(2) - & ,ELDT(2),ELPT(2),ICDJT(1),ICDJT(2) - 1533 FORMAT(/14X,1P8E12.4,2I8///) - DO I=1,L-1 - ILD(I)=IDINT(XX(I)/CW+0.01D0) - IF(ILD(I).GT.MAXD) ILD(I)=MAXD - DO J=1,ILD(I) - DLI(I)=DLI(I)+DMGN(J) - ENDDO - ENDDO - DLI(L)=TDMGN - DO 1493 I=L,2,-1 - DLI(I)=DLI(I)-DLI(I-1) - 1493 CONTINUE - DO 1494 I=1,L - WRITE(21,1495) I,DLI(I) - 1495 FORMAT(/5X,'DAMAGE IN LAYER ',I3,' : ',1P1E12.4) - 1494 CONTINUE - 1455 CONTINUE - if(irl.eq.0) goto 1497 - WRITE(21,1496) - 1496 FORMAT(1H1,/,5X,'RECOILS') - WRITE(21,1597) - 1597 FORMAT(///,5X,8HDEPTH(A), 5X,10HINLOSS(EV),3X,10HTELOSS(EV),3X - & ,10HELLOSS(EV), 3X,10HDAMAGE(EV),3X,10HPHONON(EV),5X,3HDPA, - & 2X,6HDPA(1),2X,6HDPA(2), 1X,5H(1-1),1X,5H(1-2),1X,5H(2-1),1X - & ,5H(2-2)/) - D1=0.D0 - D2=CW - do im3=MAXD,1,-1 - if (.not.equal(ionr(im3),0.D0)) go to 31 - enddo - im3=1 - 31 im3=MIN0(im3+2,MAXD) - DO I=1,im3 - WRITE(21,1595) D1,D2,IONR(I),DENTR(I),DMGNR(I),ELGDR(I), - & PHONR(I),ICDTR(I),ICDIRN(I,1),ICDIRN(I,2) - & ,ICDIRI(I,1,1),ICDIRI(I,1,2),ICDIRI(I,2,1),ICDIRI(I,2,2) - 1595 FORMAT(1X,F6.0,1H-,F6.0,1P1E14.4,4E13.4,3I8,4I6) - D1=D2 - D2=D2+CW - ENDDO - WRITE(21,1596) TIONR,TDENTR,TDMGNR,TELGDR,TPHONR - & ,ICDTTR,ICDITR(1),ICDITR(2) - 1596 FORMAT(/14X,1P1E14.4,4E13.4,3I8) - 1497 continue + call writeImplantationProfile(depth_interval_flag,YH,MAXD, + & MAXD1,LJ,NJ,L,XX,CW,IRL,IRP,RIRP,IPL,ION,DENT,DMGN, + & ELGD,PHON,CASMOT,ICD,ICDT,ICDR,ICDTR,ICDIRI,ICDIRN, + & IIRP,TRIRP,IIPL,TION,TDENT,TDMGN,TELGD,TPHON,TCASMO, + & ICDTT,TIONR,TDENTR,TELGDR,TDMGNR,TPHONR,ICDTTR,ELE, + & ELI,ELP,ELD,ELET,ELIT,ELPT,ELDT,ICDJT,ICDJTR,ICDITR, + & ILD,DLI,IONR,DENTR,DMGNR,ELGDR,PHONR,MAXNL,MAXNL5) C C BACKSCATTERING C @@ -3809,6 +3650,209 @@ C======================================================================= C OUTPUT WRITER SUBROUTINES C======================================================================= +C======================================================================= +C WRITE IMPLANTATION PROFILE AND RECOIL DEPTH TABLES +C======================================================================= + subroutine writeImplantationProfile(depth_interval_flag,YH, + & MAXD,MAXD1,LJ,NJ,L,XX,CW,IRL,IRP,RIRP,IPL,ION,DENT, + & DMGN,ELGD,PHON,CASMOT,ICD,ICDT,ICDR,ICDTR,ICDIRI, + & ICDIRN,IIRP,TRIRP,IIPL,TION,TDENT,TDMGN,TELGD,TPHON, + & TCASMO,ICDTT,TIONR,TDENTR,TELGDR,TDMGNR,TPHONR,ICDTTR, + & ELE,ELI,ELP,ELD,ELET,ELIT,ELPT,ELDT,ICDJT,ICDJTR, + & ICDITR,ILD,DLI,IONR,DENTR,DMGNR,ELGDR,PHONR,MAXNL, + & MAXNL5) + implicit none + integer MAXD,MAXD1,LJ,L,IRL,MAXNL,MAXNL5 + integer depth_interval_flag + integer NJ(MAXNL),IRP(0:MAXD1),IPL(MAXD) + integer ICD(MAXD,MAXNL5),ICDT(MAXD),ICDR(MAXD,MAXNL5) + integer ICDTR(MAXD),ICDIRI(MAXD,MAXNL5,MAXNL5) + integer ICDIRN(MAXD,MAXNL5),ICDJT(MAXNL5),ICDJTR(MAXNL5) + integer ICDITR(MAXNL5),ILD(MAXNL) + integer IIRP,IIPL,ICDTT,ICDTTR + integer I,J,K,im1,im2,im3 + real*8 YH,CW,XX(MAXNL) + real*8 RIRP(0:MAXD1),ION(MAXD),DENT(MAXD),DMGN(MAXD) + real*8 ELGD(MAXD),PHON(MAXD),CASMOT(MAXD) + real*8 TION,TDENT,TDMGN,TELGD,TPHON,TCASMO,TRIRP + real*8 IONR(MAXD),DENTR(MAXD),DMGNR(MAXD),ELGDR(MAXD) + real*8 PHONR(MAXD),TIONR,TDENTR,TDMGNR,TELGDR,TPHONR + real*8 ELE(MAXD,MAXNL5),ELI(MAXD,MAXNL5) + real*8 ELP(MAXD,MAXNL5),ELD(MAXD,MAXNL5) + real*8 ELET(MAXNL5),ELIT(MAXNL5),ELPT(MAXNL5),ELDT(MAXNL5) + real*8 DLI(MAXNL) + real*8 D1,D2,Dmid + logical EQUAL + + 590 WRITE(21,600) + 600 FORMAT(1H1,///,5X,8HDEPTH(A),2X,9HPARTICLES,2X,10HNORM.DEPTH,1X, + & 10HPATHLENGTH,3X,10HINLOSS(EV),2X,10HTELOSS(EV),2X + & ,10HELLOSS(EV), 2X,10HDAMAGE(EV),2X,10HPHONON(EV),2X + & ,10HCASCAD(EV),5X,3HDPA/) + + IF(depth_interval_flag.EQ.0) THEN +C WRITE(22,*)' CALCULATED IMPLANTATION PROFILE DID NOT ', +C & 'AGREE WITH LAYER THICKNESS' + WRITE(21,*)' CALCULATED IMPLANTATION PROFILE DID NOT ', + & 'AGREE WITH LAYER THICKNESS' + ENDIF + +#if defined (OS_WIN) + WRITE(22,6002) + 6002 FORMAT(1H1,///,5X,8HDEPTH(A),2X,9HPARTICLES) +#else + write(22,'(a)') ' DEPTH PARTICLES' +#endif + IF(YH.LT.1.D0) GO TO 603 + DO I=0,MAXD1 +C This normalization of the implantation profile is wrong if the +C simulation does not account for all implanted projectiles, +C e.g. when the number of points in the profile is not enough to +C reach the maximum implantation depth. + RIRP(I) = DBLE(IRP(I))/YH +C To resolve this issue we need to check whether we reached maximum +C implantation or not. If not rerun calculation with larger step +C size! + ENDDO + 603 D1=0. + D2=CW + WRITE(21,601) D1,IRP(0),RIRP(0) + 601 FORMAT(4X,3H-SU,1H-,F6.0,I10,E12.4) + DO J=1,LJ + DO I=1,MAXD + ICDT(I)=ICDT(I)+ICD(I,J) + ICDTR(I)=ICDTR(I)+ICDR(I,J) + ENDDO + ENDDO + DO K=1,NJ(1) + DO J=1,NJ(1) + DO I=1,MAXD + ICDIRN(I,J)=ICDIRN(I,J)+ICDIRI(I,K,J) + ENDDO + ENDDO + ENDDO + DO I=0,MAXD1 + IIRP=IIRP+IRP(I) + TRIRP=TRIRP+RIRP(I) + ENDDO + DO I=1,MAXD + IIPL=IIPL+IPL(I) + TION=TION+ION(I) + TDENT=TDENT+DENT(I) + TDMGN=TDMGN+DMGN(I) + TELGD=TELGD+ELGD(I) + TPHON=TPHON+PHON(I) + TCASMO=TCASMO+CASMOT(I) + ICDTT=ICDTT+ICDT(I) + TIONR=TIONR+IONR(I) + TDENTR=TDENTR+DENTR(I) + TELGDR=TELGDR+ELGDR(I) + TDMGNR=TDMGNR+DMGNR(I) + TPHONR=TPHONR+PHONR(I) + ICDTTR=ICDTTR+ICDTR(I) + ENDDO + do im1=MAXD,1,-1 + if(ipl(im1).ne.0.or.(.NOT.EQUAL(ion(im1),0.D0))) goto 20 + enddo + im1=1 + 20 im1=min0(im1+2,MAXD) + DO I=1,im1 + WRITE(21,700) D1,D2,IRP(I),RIRP(I),IPL(I),ION(I),DENT(I), + & DMGN(I),ELGD(I),PHON(I),CASMOT(I),ICDT(I) + Dmid=(D2-D1)/2+D1 + WRITE(22,701) Dmid,IRP(I) + 700 FORMAT(1X,F6.0,1H-,F6.0,I10,E12.4,I10,1P1E14.4,5E12.4,I8) + 701 FORMAT(1X, F16.4, 2X, I10) + D1=D2 + D2=D2+CW + ENDDO + WRITE(21,604) D2-CW,IRP(MAXD1),RIRP(MAXD1) + 604 FORMAT(1X,F6.0,1H-,3X,3HSUT,I10,E12.4) + WRITE(21,710) IIRP,TRIRP,IIPL,TION,TDENT,TDMGN,TELGD,TPHON,TCASMO + & ,ICDTT + 710 FORMAT(/14X,I10,1P1E12.4,I10,1E14.4,5E12.4,I8) + DO J=1,NJ(1) + DO I=1,MAXD + ELET(J)=ELET(J)+ELE(I,J) + ELIT(J)=ELIT(J)+ELI(I,J) + ELPT(J)=ELPT(J)+ELP(I,J) + ELDT(J)=ELDT(J)+ELD(I,J) + ICDJT(J)=ICDJT(J)+ICD(I,J) + ICDJTR(J)=ICDJTR(J)+ICDR(I,J) + ICDITR(J)=ICDITR(J)+ICDIRN(I,J) + ENDDO + ENDDO + + WRITE(21,1521) + 1521 FORMAT(1H1,4X,'DEPTH(A)' + & ,3X,' INLOSS(1)',3X,'ELLOSS(1)',3X,'DAMAGE(1)',3X,'PHONON(1)' + & ,2X,' INLOSS(2)',3X,'ELLOSS(2)',3X,'DAMAGE(2)',3X,'PHONON(2)' + & ,2X,'DPA(1)',2X,'DPA(2)'/) + D1=0. + D2=CW + do im2=MAXD,1,-1 + if(.NOT.EQUAL(eli(im2,1),0.D0).or.(.NOT.EQUAL(eli(im2,2),0.D0)) + & ) goto 30 + + enddo + im2=1 + 30 im2=MIN0(im2+2,MAXD) + DO 1525 I=1,im2 + WRITE(21,1523) D1,D2,ELI(I,1),ELE(I,1),ELD(I,1),ELP(I,1), + & ELI(I,2),ELE(I,2),ELD(I,2),ELP(I,2),ICD(I,1),ICD(I,2) + 1523 FORMAT(1X,F6.0,1H-,F6.0,1P8E12.4,2I8) + D1=D2 + D2=D2+CW + 1525 CONTINUE + WRITE(21,1533) ELIT(1),ELET(1),ELDT(1),ELPT(1),ELIT(2),ELET(2) + & ,ELDT(2),ELPT(2),ICDJT(1),ICDJT(2) + 1533 FORMAT(/14X,1P8E12.4,2I8///) + DO I=1,L-1 + ILD(I)=IDINT(XX(I)/CW+0.01D0) + IF(ILD(I).GT.MAXD) ILD(I)=MAXD + DO J=1,ILD(I) + DLI(I)=DLI(I)+DMGN(J) + ENDDO + ENDDO + DLI(L)=TDMGN + DO 1493 I=L,2,-1 + DLI(I)=DLI(I)-DLI(I-1) + 1493 CONTINUE + DO 1494 I=1,L + WRITE(21,1495) I,DLI(I) + 1495 FORMAT(/5X,'DAMAGE IN LAYER ',I3,' : ',1P1E12.4) + 1494 CONTINUE + 1455 CONTINUE + if(irl.eq.0) goto 1497 + WRITE(21,1496) + 1496 FORMAT(1H1,/,5X,'RECOILS') + WRITE(21,1597) + 1597 FORMAT(///,5X,8HDEPTH(A), 5X,10HINLOSS(EV),3X,10HTELOSS(EV),3X + & ,10HELLOSS(EV), 3X,10HDAMAGE(EV),3X,10HPHONON(EV),5X,3HDPA, + & 2X,6HDPA(1),2X,6HDPA(2), 1X,5H(1-1),1X,5H(1-2),1X,5H(2-1),1X + & ,5H(2-2)/) + D1=0.D0 + D2=CW + do im3=MAXD,1,-1 + if (.not.equal(ionr(im3),0.D0)) go to 31 + enddo + im3=1 + 31 im3=MIN0(im3+2,MAXD) + DO I=1,im3 + WRITE(21,1595) D1,D2,IONR(I),DENTR(I),DMGNR(I),ELGDR(I), + & PHONR(I),ICDTR(I),ICDIRN(I,1),ICDIRN(I,2) + & ,ICDIRI(I,1,1),ICDIRI(I,1,2),ICDIRI(I,2,1),ICDIRI(I,2,2) + 1595 FORMAT(1X,F6.0,1H-,F6.0,1P1E14.4,4E13.4,3I8,4I6) + D1=D2 + D2=D2+CW + ENDDO + WRITE(21,1596) TIONR,TDENTR,TDMGNR,TELGDR,TPHONR + & ,ICDTTR,ICDITR(1),ICDITR(2) + 1596 FORMAT(/14X,1P1E14.4,4E13.4,3I8) + 1497 continue + return + end + C======================================================================= C WRITE INPUT AND BEAM SUMMARY C=======================================================================