PROGRAM PNDIO1 IMPLICIT REAL*8 (A-H,O-Z) REAL*4 CPUT, CPUTT C PARAMETER(MJ=10,MJ1=20,MJM=40,MJ2=20,MJ3=10, PARAMETER(MJM=40,MJ1=MJM/2,MJ2=MJ1,MJ=MJM/4,MJ3=MJ, + LY0=1,LY1=LY0+3,LY2=LY1+2,LY3=LY2+3,LYM=10, C + MK=20,MK1=20,MKM=20, + MKM=MJ1,MK1=MKM,MK=MKM, + KS0=1,KS1=3,KSM=3, + M1=LY0*MJ+(LY1-LY0)*MJ1+(LY2-LY1)*MJM+(LY3-LY2)*MJ2+ + (LYM-LY3)*MJ3+1, + M2=2*(KS0*MK+KS1*MK1)+KSM*MKM-1, CC M1=MJ*10+1,M2=MJ*11-1, + N=M1*M2, + NL=N, ML1=M1, N1DB=100, + M0Q=LY0*MJ+(LY1-LY0)*MJ1+MJM+1) C C DIMENSION IXS(M1),IYS(M2) +, AAG(N), B1(N) +, F(N+M1), R(N),WK(N+M1+1), BW(-M1:N+M1) +, AA(N),AB(-M1:N),AC(-M1:N),DI(-M1:N),P(-M1:N+M1) +, V(-M1:N+M1), W(-M1:N+M1),X(-M1:N+M1),AP(N),UB(-M1:N),AE(-M1:N) +, AB1(-M1:N+M1),AC1(-M1:N+M1),UB1(-M1:N+M1),AE1(-M1:N+M1) +, AF(-M1:N),AF1(-M1:N+M1) +, RRA(-M1:N+M1),PA(-M1:N+M1),RR(-M1:N+M1), R0(N) +, DXK(M2+1),DMJMKV(M2+1),DMJPKV(M2+1) +, PSI(-ML1:NL+ML1), DYM(NL) +, DP(N), QNA(N),QND(N), DN(N), DNL(-ML1:NL+ML1) +, DPL(-ML1:NL+ML1) +, XID(M2+M1), QNDMN(N+M1),QNAMN(N+M1) +, WKR(N+M1), XIDP(M2+M1) +, DHX(M2+1),DHY(M1) +, DYI(M1),DIKL(M1),DIKR(M1),DKIL(M1),DKIR(M1) C DIMENSION VBB(N1DB),VDB(N1DB) +, DMOBN(NL+ML1), EN(NL+ML1), ENY(NL+ML1), DMOBP(NL+ML1) +, EJN(NL+ML1), EJP(NL+ML1) +, ENX(NL+ML1), AENX(NL+ML1), AENY(NL+ML1) +, GR(NL+ML1), GRWK(NL+ML1), GRU(NL+ML1) +, DJN(NL+ML1), DJP(NL+ML1) +, DJNH(NL+ML1) +, DVN(NL+ML1) +, DVNX(NL+ML1) +, YMK(M1) C +, DFN(NL+ML1), DFP(NL+ML1), EINSN(NL+ML1),EINSP(NL+ML1) C CCC CCC* C DATA EPSN/0.500D-9 /, EPSNL/0.50D-9 /, EPSNP/1.00D-5 / C DATA EPSN/0.500D-9 /, EPSNL/0.50D-4 /, EPSNP/1.00D-5 / C1 DATA EPSN/0.500D-10 /, EPSNL/0.50D-11 /, EPSNP/1.00D-11 / DATA EPSN/0.500D-9 /, EPSNL/0.50D-9 /, EPSNP/1.00D-9 / C +, NKAI/ 300 / +, NKAI/ 300 / CC CALL CLOCK0 C**2A EPSDIV= 2.0D+1 C1 EPSDIV= 1.0D+3 C1 EPSDIL= 1.0D+3 EPSDIL= 1.0D+1 VB1= 0.00D0 VB2= 0.00D0 VBSP= 0.00D0 EPSNO= 2.00D-3 C1 EPSNO= 1.00D-4 C? EPSNO= 1.00D-2 C EPSNO= 0.40D-2 VD1= 0.000D0 VD2= -0.6000D0 VDSP= -0.6000D0 EPSDIP= 1.0D+1 C1 EPSDIP= 1.0D+3 CC VBACKK = 0.00D0 LWRT=0 C GAMAF=0.38D0 GAMAF=0.0D0 IGRST=1 * 1:ENX, -1:ENY, 0:zero * GAMA= 0.190D0 GAMAP= 0.190D0 C GAMA= 0.257D0 C GAMAP= 0.257D0 C GAMA=0.0D0 C GAMAP=0.0D0 C C------------------------------------------------------------- C DPND= 7.50D+19 C DPNDMI= 7.00D+19 C FDIEL= DPNDMI/ DPND C FDIEL=0.25D0 YWID=0.1D0 INDST=0 * 0:step, 1:slope, 2:exp *** INAST=0 C INAST=INDST C IRATI=0 * GAMA,GAMAP you 0: kotei*** C YREFCM= 1.0D-6 BET0MB= 0.5D0 FVNMB= 1.0D0 FVNMBP= 1.0D0 TAWNFI=0.4D-12 TWN10J=TAWNFI*1.0D10 FCTUTN=1.0D0 C C IGUM=1 IDBG=0 IDBGU=0 ISTB=1 * ninput etc IPIN=1 ININ=1 * 0:CALC C CALL BIASP(VB1,VB2,VBSP, VD1,VD2,VDSP, VBB,VDB,N1DB, NOB) C USHIRP=0.95D0 CC USHIROB=0.0D0; USHIRP=0.98D0(BEST) UTOP=1.01D0 UTOPN=1.01D0 C SI=12.0D0 C SIO2=SI/3.0D0 CC ZWW= 20.0D0 CC ZWW= 90.0D0 ZWW= 380.0D0 C EINS= 1.0D0/2.58510D-2 C/ TRMC= 1.0D-15 C C------------------------------------------------------------- C CC DPNDH= 7.50D+19 CC DPNDHM= DPNDH*TRMC C DPND= 7.50D+19 C CC DPNAZ= 0.58D+18 DPNA= 0.58D+18 C SINI= 1.5D+10 DPNDM= DPND*TRMC DPNAM= DPNA*TRMC CC DPNAZM= DPNAZ*TRMC SINIM= SINI*TRMC DNDOCM= (DPNDM+DSQRT(DPNDM*DPNDM + +4.0D0*SINIM*SINIM))*0.5D0 DPOCM= (SINIM*SINIM)/DNDOCM C C DNAOCM= (DPNAM+DSQRT(DPNAM*DPNAM + +4.0D0*SINIM*SINIM))*0.5D0 DNBACK= (SINIM*SINIM)/DNAOCM DPBACK= DNAOCM C C VOMCN= DLOG(DNDOCM/SINIM)/EINS VOMCP= DLOG(DNAOCM/SINIM)/EINS C WRITE(6,*) ' % VOMCN, VOMCP, DNDOCM, DNAOCM= ' + , VOMCN, VOMCP, DNDOCM, DNAOCM C C CALL DXKLST(M2, MK,MK1,MKM, KXA,KXB,KXC,KX1, KS0,KS1,KSM + , DXK,DHX, KXBL,KXCR) C CALL DYILST(M1, MJ,MJ1,MJM, MJ2,MJ3, LY0,LY1,LY2,LY3 + , JY0,JY1,JY2,JY3, MSU,MSM,MSD, DYI,DHY, M0Q) CC CC J1=M1*(KXB+MKM-1) +MSD J2=M1*(KXB+2*MKM-1) +MSD J3=M1*(KXB+MKM-1) +MSM J4=M1*(KXB+2*MKM-1) +MSM J5=M1*(KXB+MKM-1) +MSU J6=M1*(KXB+2*MKM-1) +MSU C WRITE(6,*) ' % M1= ',M1,' M2= ',M2,' N= ',N WRITE(6,*) ' % MJ,MJ1,MJM,MJ2,MJ3= ',MJ,MJ1,MJM,MJ2,MJ3 WRITE(6,*) ' % MK,MK1,MKM= ',MK,MK1,MKM WRITE(6,*) ' % EPSNO= ',EPSNO WRITE(6,*) ' % EINS, IGRST= ',EINS, IGRST C WRITE(6,*) ' % J1,J2,J3,J4,J5,J6, M0Q= ',J1,J2,J3,J4,J5,J6,M0Q WRITE(6,*) ' % DXK(KXB),DXK(KXC), DYI(JY1),DYI(JY2)= ' + , DXK(KXB),DXK(KXC), DYI(JY1),DYI(JY2) WRITE(6,*) ' % DHX(KXB),DHX(KXC), DHY(JY1),DHY(JY2)= ' + , DHX(KXB),DHX(KXC), DHY(JY1),DHY(JY2) C C CALL NDDOPE( N,M1,M2,M0Q, QND, DPNDM, WK,QNDMN + ,YMK,MJM,MJ1, MSD,INDST, FDIEL,YWID, JNCUT) WRITE(6,*) ' % QNDMN(J1),QNDMN(J3),QNDMN(J5)= ' + , QNDMN(J1),QNDMN(J3),QNDMN(J5) C CALL NADOPE( N,M1,M2,M0Q, QNA, DPNAM, WK,QNAMN + ,YMK,MJM,MJ2, MSU,INAST, FDIEL,YWID, JPCUT) WRITE(6,*) ' % QNAMN(J1),QNAMN(J3),QNAMN(J5)= ' + , QNAMN(J1),QNAMN(J3),QNAMN(J5) C CALL SYOKIB(N,M2,M1,N, DNL,SINIM, PSI,EINS, QND,QNA, DPL +, DN,DP) WRITE(6,*) ' % PSI(J1),PSI(J3),PSI(J5)= ' + , PSI(J1),PSI(J3),PSI(J5) C/ CALL MKEINS(N,M1,M2, EINS, EINSN,EINSP) CC C C------------------------------------------------------------- DO 682 I=0, M2-1 DO 682 J=1,M1 K=M1*I+J 682 BW(K)=PSI(K) C CC * *** BIAS STEP *** * INB=0 KCTB=0 INLB=0 KCTLB=0 INPB=0 KCTPB=0 CC VDPRE=0.0D0 VBPRE=0.0D0 CC CC DO 5111 IB=1, NOB VBACKK=VBB(IB) VDRAIK=VDB(IB) CC/ VBACKG= VBACKK-VOMCP VDRAIN= VDRAIK+VOMCN CC/ CALL VDGETA(N,M1,M2, VDRAIK, VDPRE,QND, BW, VOMCN) VDPRE=VDRAIK CC CALL VBGETA(N,M1,M2, VBACKK, VBPRE,QNA, BW, VOMCP) VBPRE=VBACKK C * *** START *** * INO=0 KCTO=0 INLO=0 KCTLO=0 INPO=0 KCTPO=0 CC KCTOO=0 KCTLOO=0 KCTPOO=0 CC INDWRITE=0 CPUTT=0.0 CC CC DO 1230 IN= 1, NKAI * CALL POISON(N,M1,B1,F,R, AA,AB,AC,DI,P,W,X,AP,UB,AE +, AB1,AC1,UB1,AE1,RRA,PA,RR, AAG, IN +, M2, WK, DXK,SI +, DN,DP,M1, IGUM, VDRAIN +, VBACKG, DMJMKV,DMJPKV, J1,J2,J3,J4,J5,J6, IDBG +, BW, EPSN,EPSDIV, KCTO,INO +, QNDMN,QNAMN, WKR, DHX,DHY +, DYI,DIKL,DIKR,DKIL,DKIR, CNRP, USHIRP +, EINSN,EINSP, N) CC ** DO 280 I=0,M2-1 DO 280 J=1,M1 K= M1*I+J 280 PSI(K)= BW(K) * * CC CALL ENCAL(N,M1,M2,PSI, EN, DXK,DHY,M1, ENY + ,DNL,DPL, DFN,DFP + ,EINSN,EINSP, GAMAF, ENX) C/ CC CC IF(IRATI.EQ.1) THEN RATIO=VDRAIN**8/(1.0D0+VDRAIN**8) * RATIO=VGATE**8/(1.0D0+VGATE**8) ELSE RATIO=1.0D0 ENDIF GAMAW= GAMA*RATIO GAMAPW=GAMAP*RATIO C WRITE(6,*) ' % VDRAIK,VDRAIN, VBACKK,VBACKG= ' C + ,VDRAIK,VDRAIN, VBACKK,VBACKG C WRITE(6,*) ' % GAMA,GAMAP, RATIO, GAMAW,GAMAPW= ' C + ,GAMA,GAMAP, RATIO, GAMAW,GAMAPW CC CC C IMOBTN= 0 IUTNM=1 C CALL MOBCAL(N,M1,M2, DMOBN, QNAMN,QNDMN,N,M1, ENY +, DFN, EINSN, GAMAW, IUTNM, EINS +, YREFCM, BET0MB, FVNMB, IMOBTN,TWN10J +, FCTUTN) CC IUTN=1 C CALL MOBCAP(N,M1,M2, DMOBP, QNAMN,QNDMN,N,M1, ENY +, DFP, EINSP, GAMAPW, IUTN, EINS +, YREFCM, BET0MB, FVNMBP) CC CCC============================ kokomade ========================= CC CC IF(ININ.EQ.0) THEN CALL NINPUC(N,M1,M2, DNL,PSI +, EINS,SINIM, VBACKK,VDRAIK,QND, QNA, DNDOCM, DNBACK) ELSE CALL NINPUT(N,M1,M2, DNL,PSI, DNDOCM +, B1,F,R,WK, AA,AB,AC,DI,P,W,X,AP,UB,AE +, AB1,AC1,UB1,AE1,RRA,PA,RR,DXK +, EPSNL, IN,AAG +, EPSDIL, KCTLO,INLO,CNRLO, DNBACK +, UTOPN, IDBG, V,R0,ISTB,AF,AF1 +, DYI,DIKL,DIKR,DKIL,DKIR, M1, DMOBN +, GR, EINSN, J1,J2,J3,J4,J5,J6) C CC ENDIF CC/ DO 510 I= 0, M2-1 DO 510 J= 1, M1 K= M1*I+J 510 DN(K)= DNL(K) ** * IF(IPIN.EQ.0) THEN CALL PINPUT(N,M1,M2, DPL,PSI +, EINS,SINIM, VBACKK,VDRAIK,QND, QNA, DPOCM, DPBACK) ELSE C CALL PINPU1(N,M1,M2, DPL,PSI, DPOCM +, B1,F,R,WK, AA,AB,AC,DI,P,W,X,AP,UB,AE +, AB1,AC1,UB1,AE1,RRA,PA,RR,DXK +, EPSNP, IN,AAG +, EPSDIP, KCTPO,INPO,CNRPO +, DPBACK, UTOP, IDBG, V,R0,ISTB,AF,AF1 +, DYI,DIKL,DIKR,DKIL,DKIR, M1, DMOBP +, GR, EINSP, J1,J2,J3,J4,J5,J6) C ENDIF CC/ DO 504 I= 0, M2-1 DO 504 J= 1, M1 K= M1*I+J 504 DP(K)= DPL(K) ** ** DO 333 I=1,N 333 WK(I)=DP(I)*DN(I) ZAN=0.0D0 DO 334 I=0,M2-1 DO 334 J=1,M1 KP=M1*I+J 334 ZAN=ZAN+(WK(KP)-SINIM*SINIM)**2 ZAN=DSQRT(ZAN) C# C# CCC==================================================C CCC==================================================C CC */ CALL JNCAL(N,M1,M2, WK,DXK, DMOBN +, PSI,DNL +, DJN, DYI,M1, EINSN, EJN, ENX,AENX, DJNH,DVN,DVNX +, ENY,AENY) C CALL JPCAL(N,M1,M2, WK,DXK, DMOBP +, PSI,DPL +, DJP, DYI,M1, EINSP, EJP) CC CC IF(IGRST.EQ.1) THEN CALL GRSET(N,M1, M2, M1 +, GR,QNAMN,QNDMN,N, SINIM,DP,DN,DHX,DHY, GRWK +, DJN,DJP,AENX, AENX , GRU) * ENX *** C ELSE IF(IGRST.EQ.-1) THEN CALL GRSET(N,M1, M2, M1 +, GR,QNAMN,QNDMN,N, SINIM,DP,DN,DHX,DHY, GRWK +, DJN,DJP,AENY, AENY , GRU) * ENY *** C ELSE DO 11 I=1,N+M1 11 GR(I)=0.0D0 C ENDIF CCC CCC */ CALL ZEROCR(N,M1,B1,F,R, AA,AB,AC,DI,P,W,X,AP,UB,AE +, AB1,AC1,UB1,AE1,RRA,PA,RR, AAG) ** IGUMD=2 CCC CALL REVMTR(N,M1,M2,AA,AB,AC,WK, DXK,SI +, AAG,DN,DP,M1, IGUMD, DHX,DHY +, DYI,DIKL,DIKR,DKIL,DKIR, EINSN,EINSP, N) * CALL FINPUT(N,M1,F,WK,VDRAIN, M2 +, SI,VBACKG,DXK,DMJMKV,DMJPKV, IGUMD, DP,DN, M1, BW +, AC,AB,AAG, SNRM, QNDMN,QNAMN, WKR, DHX,DHY +, DYI, EINSN,EINSP, N) * CNRO=0.0D0 DO 136 I=1,N 136 CNRO=CNRO+F(I)*F(I) ERRO=DSQRT(CNRO)/SNRM C CCC CCC CCC CALL IDINTB(N,M1,M2, WK,DXK, GID,DMOBN * --- +, DYM,PSI,DNL, XID +, DYI,M1, EINSN, ZWW) C C CALL IDINTP(N,M1,M2, WK,DXK, GIDP,DMOBP +, DYM,PSI,DPL, XIDP +, DYI,M1, EINSP, ZWW) CC CALL IDVGTB(N,M1,M2, WK,DXK, GIDVG,DMOBN +, DYM,PSI,DNL +, DYI,M1, EINSN, ZWW) CC CALL IDVGTP(N,M1,M2, WK,DXK, GIDVGP,DMOBP +, DYM,PSI,DPL +, DYI,M1, EINSP, ZWW) CC * INW=IN * IF(LWRT.EQ.1) THEN CC * IF((ERRO.LT.EPSNO*4.0).AND.(INDWRITE.EQ.0) ) THEN KCTOI=KCTO-KCTOO KCTLOI=KCTLO-KCTLOO KCTPOI=KCTPO-KCTPOO INI=INO-INOO C CALL OUTWR(M1,N, ZAN, INW,J1,J2,J3,J4,J5,J6 +, CNRP,INO,KCTO,ERRO, CNRLO,INLO,KCTLO +, VDRAIK,GID, BW,DN +, KXC,KX1,KXA,KXB, KXBL,KXCR +, IB, VDRAIN, INB,INLB,KCTB,KCTLB +, DP,CNRPO,INPO,KCTPO,INPB,KCTPB, GIDP +, GIDVG,GIDVGP,KCTOI,KCTLOI,KCTPOI,INI +, VBACKK,VBACKG, MSM) CC * KCTOO=KCTO KCTLOO=KCTLO KCTPOO=KCTPO INOO=INO GIDO=GID INDWRITE=INDWRITE+1 * CALL CLOCK(CPUT) CPUTT=CPUTT+CPUT WRITE(6,*) ' *4 EPS CPUTIME = ', CPUT,CPUTT,GID,INO CALL CLOCK0 END IF CC CC IF((ERRO.LT.EPSNO*1.0 ).AND.(INDWRITE.EQ.1) ) THEN KCTOI=KCTO-KCTOO KCTLOI=KCTLO-KCTLOO KCTPOI=KCTPO-KCTPOO INI=INO-INOO C CALL OUTWR(M1,N, ZAN, INW,J1,J2,J3,J4,J5,J6 +, CNRP,INO,KCTO,ERRO, CNRLO,INLO,KCTLO +, VDRAIK,GID, BW,DN +, KXC,KX1,KXA,KXB, KXBL,KXCR +, IB, VDRAIN, INB,INLB,KCTB,KCTLB +, DP,CNRPO,INPO,KCTPO,INPB,KCTPB, GIDP +, GIDVG,GIDVGP,KCTOI,KCTLOI,KCTPOI,INI +, VBACKK,VBACKG, MSM) C * KCTOO=KCTO KCTLOO=KCTLO KCTPOO=KCTPO INOO=INO GIDE=GID+(GID-GIDO) GIDO=GID INDWRITE=INDWRITE+1 * CALL CLOCK(CPUT) CPUTT=CPUTT+CPUT WRITE(6,*) ' *1 EPS CPUTIME = ', CPUT,CPUTT,GID,INO,GIDE CALL CLOCK0 END IF END IF * CC IF(ERRO.LT.EPSNO) GO TO 1250 * 1230 CONTINUE 1250 CONTINUE * VN=0.0D0 NNEG=0 VP=0.0D0 NEGP=0 * DO 521 J= 0, M2-1 DO 521 K= 1, ML1 I= M1*J+K IL=ML1*J+K IF(DN(I).LT.VN) THEN VN= DN(I) ENDIF IF(DN(I).LT.1.0D-56) THEN DN(I)=1.0D-56 DNL(IL)=1.0D-56 NNEG=NNEG+1 ENDIF CC IF(DP(I).LT.VP) THEN VP= DP(I) ENDIF IF(DP(I).LT.1.0D-56) THEN DP(I)=1.0D-56 DPL(IL)=1.0D-56 NEGP=NEGP+1 ENDIF 521 CONTINUE IF(VN.LT.0.0D0 ) THEN WRITE(*,*) '*** DN IS NEGATIVE !! NNEG,VN= ***',NNEG,VN,'***' ELSE WRITE(*,*) ' NNEG,VN= ',NNEG,VN,'***' END IF CC IF(VP.LT.0.0D0 ) THEN WRITE(*,*) '*** DP IS NEGATIVE !! NEGP,VP= ***',NEGP,VP,'***' ELSE WRITE(*,*) ' NEGP,VP= ',NEGP,VP,'***' END IF C INB=INB+INO KCTB=KCTB+KCTO INLB=INLB+INLO KCTLB=KCTLB+KCTLO INPB=INPB+INPO KCTPB=KCTPB+KCTPO C# CCCC IF(IB.EQ.NOB) THEN WRITE(6,*) ' ' CALL BPNWR(M1,N, JY0,MJM,JY1,JY2,JY3 +, KXC,KX1,KXA,KXB, KXBL,KXCR +, BW,DN,DP) CC CC CALL GRWKWR(M1,N, JY0,MJM,JY1,JY2,JY3 +, KXC,KX1,KXA,KXB, KXBL,KXCR +, EINSN,EINSP, GRWK, NL,ML1, M2) CC ENDIF C/ C * Debug You *** IF(IDBGU .EQ. 2) THEN TTT=100.0D0 DO 894 I=2,M2 DO 894 J=2,ML1 K=ML1*(I-1)+J B1KP=11610.145D0/EINSN(K) IF (B1KP .LT. TTT) THEN WRITE(6,*) '<MJM:', MJ2,MJM WRITE(6,*) ' ** BW,DN,DP,TN,TP,ENY,DMOBN,DMOBP,DFN **' IX=KXBL DO 111 J=JY2+MJ2, JY1+MJM, -1 I=M1*(IX-1)+J WRITE(6,2021) J,BW(I),DN(I),DP(I) +, (11610.145D0/EINSN(I)) +, (11610.145D0/EINSP(I)), ENY(I),DMOBN(I),DMOBP(I),DFN(I) +, DN(I)*DP(I) +, SINIM*SINIM 111 CONTINUE WRITE(6,*) ' ---' WRITE(6,*) ' MJM->MJ1:', MJM,MJ1 WRITE(6,*) ' ** BW,DN,DP,TN,TP,ENY,DMOBN,DMOBP,DFN **' C DO 113 J=JY1+MJM, JY0+2*MJ1, -1 I=M1*(IX-1)+J WRITE(6,2021) J,BW(I),DN(I),DP(I) +, (11610.145D0/EINSN(I)) +, (11610.145D0/EINSP(I)), ENY(I),DMOBN(I),DMOBP(I),DFN(I) +, DN(I)*DP(I) +, SINIM*SINIM 113 CONTINUE C WRITE(6,*) ' %%%' WRITE(6,*) ' %%% IX=2,M2' WRITE(6,*) ' ** BW,DN,DP,TN,TP,ENY,DMOBN,DMOBP,DFN,ENX **' C DO 3331 IX=2,M2 I=M1*(IX-1)+MSM+1 WRITE(6,2030) IX,BW(I),DN(I),DP(I) C +, (11610.145D0/EINSN(I)) C +, (11610.145D0/EINSP(I)) +, (EINSN(I)) +, (EINSP(I)) +, ENY(I),DMOBN(I),DMOBP(I), DFN(I),ENX(I) 3331 CONTINUE C C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 2029 FORMAT(1H ,I4,2F12.5,4D14.6) WRITE(6,*) ' ++++++++++++++++++++++++++' WRITE(6,*) ' MSU+MJ2->M0Q,JPCUT:', MSU+MJ2,M0Q,JPCUT WRITE(6,*) ' MJ2->MJM:', MJ2,MJM WRITE(6,*) ' FDIEL,YWID,INDST,INAST:', FDIEL,YWID,INDST,INAST WRITE(6,*) ' ** J,YMK,BW,QNDMN,QNAMN, DN,DP **' IX=KXBL DO 161 J=JY2+MJ2, JY1+MJM, -1 I=M1*(IX-1)+J WRITE(6,2029) J,YWID-YMK(J),BW(I),QNDMN(I),QNAMN(I) + ,DN(I),DP(I) 161 CONTINUE WRITE(6,*) ' ++++' WRITE(6,*) ' MSD-MJ1->M0Q,JNCUT:', MSD-MJ1,M0Q,JNCUT WRITE(6,*) ' MJM->MJ1:', MJM,MJ1 WRITE(6,*) ' ** J,YMK,BW,QNDMN,QNAMN, DN,DP **' C DO 163 J=JY1+MJM, JY0+2*MJ1, -1 I=M1*(IX-1)+J WRITE(6,2029) J,YWID-YMK(J),BW(I),QNDMN(I),QNAMN(I) + ,DN(I),DP(I) 163 CONTINUE CCC C C 1900 FORMAT(3D26.18) 9100 FORMAT(1H ,9D14.6) 9101 FORMAT(1H ,I4,9D14.6) 5820 FORMAT(1H,I4,D12.4,D12.4,D12.4,D12.4,D12.4,D12.4,D11.2,I5,I6,I5 +, I5,F7.1, D12.3 ) CC STOP END C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CC CC SUBROUTINE GRTWMX(M1,N,M2 +, GR ,NL,ML1 +, EINSN,EINSP, GRMX,GRMN,TNMX,TNMN,TPMX,TPMN +, EN,ENX,ENY, DFN,DFP, EJN,EJP) IMPLICIT REAL*8(A-H,O-Z) DIMENSION GR(NL+ML1) +, EINSN(NL+ML1), EINSP(NL+ML1) +, EN(NL+ML1), ENX(NL+ML1),ENY(NL+ML1) +, DFN(NL+ML1), DFP(NL+ML1) +, EJN(NL+ML1), EJP(NL+ML1) * * GRMX=GR(3) GRMN=GR(3) IGRMX=1 JGRMX=3 IGRMN=1 JGRMN=3 DO 511 I= 0, M2-1 DO 511 J= 1, ML1 * - K= ML1*I+J IF(GR(K).GT.GRMX) THEN GRMX=GR(K) IGRMX=I+1 JGRMX=J ENDIF IF(GR(K).LT.GRMN) THEN GRMN=GR(K) IGRMN=I+1 JGRMN=J ENDIF 511 CONTINUE * KMX=ML1*(IGRMX-1)+JGRMX CC IF(EINSN(KMX).NE.0.0D0 .AND. EINSP(KMX).NE.0.0D0) THEN WRITE(6,1010) '*GR:MAX,MIN ',IGRMX,JGRMX,GRMX, IGRMN,JGRMN,GRMN * , ENX(KMX) * , 11610.145D0/EINSN(KMX), 11610.145D0/EINSP(KMX) ELSE WRITE(6,1010) '*GR:MAX,MIN ',IGRMX,JGRMX,GRMX, IGRMN,JGRMN,GRMN * , ENX(KMX) C * , 11610.145D0/EINSN(KMX), 11610.145D0/EINSP(KMX) ENDIF C 1010 FORMAT(1H ,A12,I4,I4,D14.6, I4,I4,D14.6, D14.6, 2D14.6) CC * TNMX= 11610.145D0/EINSN(ML1+3) TNMN= 11610.145D0/EINSN(ML1+3) ITNMX=2 JTNMX=3 ITNMN=2 JTNMN=3 DO 371 I= 1, M2-1 DO 371 J= 2, ML1 K= ML1*I+J WKTN= 11610.145D0/EINSN(K) IF(WKTN.GT.TNMX) THEN TNMX=WKTN ITNMX=I+1 JTNMX=J ENDIF IF(WKTN.LT.TNMN) THEN TNMN=WKTN ITNMN=I+1 JTNMN=J ENDIF 371 CONTINUE * KMX=ML1*(ITNMX-1)+JTNMX WRITE(6,1010) '*TN:MAX,MIN ',ITNMX,JTNMX,TNMX, ITNMN,JTNMN,TNMN * , ENX(KMX) CC CC * TPMX= 11610.145D0/EINSP(ML1+3) TPMN= 11610.145D0/EINSP(ML1+3) ITPMX=2 JTPMX=3 ITPMN=2 JTPMN=3 DO 381 I= 1, M2-1 DO 381 J= 2, ML1 K= ML1*I+J WKTP= 11610.145D0/EINSP(K) IF(WKTP.GT.TPMX) THEN TPMX=WKTP ITPMX=I+1 JTPMX=J ENDIF IF(WKTP.LT.TPMN) THEN TPMN=WKTP ITPMN=I+1 JTPMN=J ENDIF 381 CONTINUE * KMX=ML1*(ITPMX-1)+JTPMX WRITE(6,1010) '*TP:MAX,MIN ',ITPMX,JTPMX,TPMX, ITPMN,JTPMN,TPMN * , ENX(KMX) CCC CC -------- C*EN TNMX= EN(ML1+3) TNMN= EN(ML1+3) ITNMX=2 JTNMX=3 ITNMN=2 JTNMN=3 DO 671 I= 1, M2-1 DO 671 J= 2, ML1 K= ML1*I+J WKTN= EN(K) IF(WKTN.GT.TNMX) THEN TNMX=WKTN ITNMX=I+1 JTNMX=J ENDIF IF(WKTN.LT.TNMN) THEN TNMN=WKTN ITNMN=I+1 JTNMN=J ENDIF 671 CONTINUE * WRITE(6,1010) '*EN:MAX,MIN ',ITNMX,JTNMX,TNMX, ITNMN,JTNMN,TNMN CC C*ENX TNMX= ENX(ML1+3) TNMN= ENX(ML1+3) ITNMX=2 JTNMX=3 ITNMN=2 JTNMN=3 DO 681 I= 1, M2-1 DO 681 J= 1, ML1 * - K= ML1*I+J WKTN= ENX(K) IF(WKTN.GT.TNMX) THEN TNMX=WKTN ITNMX=I+1 JTNMX=J ENDIF IF(WKTN.LT.TNMN) THEN TNMN=WKTN ITNMN=I+1 JTNMN=J ENDIF 681 CONTINUE * KMX=ML1*(ITNMX-1)+JTNMX KMN=ML1*(ITNMN-1)+JTNMN WRITE(6,1010) '*ENX:MAX,MIN ',ITNMX,JTNMX,TNMX, ITNMN,JTNMN,TNMN WRITE(6,2010) '-----TN(MX),TN(MN), TP(MX),TP(MN):' * , 11610.145D0/EINSN(KMX), 11610.145D0/EINSN(KMN) * , 11610.145D0/EINSP(KMX), 11610.145D0/EINSP(KMN) 2010 FORMAT(1H ,A34,2D14.6,2D14.6) CC C*ENY TNMX= ENY(ML1+3) TNMN= ENY(ML1+3) ITNMX=2 JTNMX=3 ITNMN=2 JTNMN=3 DO 691 I= 1, M2-1 DO 691 J= 2, ML1 K= ML1*I+J WKTN= ENY(K) IF(WKTN.GT.TNMX) THEN TNMX=WKTN ITNMX=I+1 JTNMX=J ENDIF IF(WKTN.LT.TNMN) THEN TNMN=WKTN ITNMN=I+1 JTNMN=J ENDIF 691 CONTINUE * WRITE(6,1010) '*ENY:MAX,MIN ',ITNMX,JTNMX,TNMX, ITNMN,JTNMN,TNMN CC C*DFN+ TNMX= DFN(ML1+3) TNMN= DFN(ML1+3) ITNMX=2 JTNMX=3 ITNMN=2 JTNMN=3 DO 1691 I= 1, M2-1 DO 1691 J= 2, ML1 K= ML1*I+J WKTN= DFN(K) IF(WKTN.GT.TNMX) THEN TNMX=WKTN ITNMX=I+1 JTNMX=J ENDIF IF(WKTN.LT.TNMN) THEN TNMN=WKTN ITNMN=I+1 JTNMN=J ENDIF 1691 CONTINUE * WRITE(6,1010) '*DFN:MAX,MIN ',ITNMX,JTNMX,TNMX, ITNMN,JTNMN,TNMN CC C*DFP+ TNMX= DFP(ML1+3) TNMN= DFP(ML1+3) ITNMX=2 JTNMX=3 ITNMN=2 JTNMN=3 DO 1091 I= 1, M2-1 DO 1091 J= 2, ML1 K= ML1*I+J WKTN= DFP(K) IF(WKTN.GT.TNMX) THEN TNMX=WKTN ITNMX=I+1 JTNMX=J ENDIF IF(WKTN.LT.TNMN) THEN TNMN=WKTN ITNMN=I+1 JTNMN=J ENDIF 1091 CONTINUE * WRITE(6,1010) '*DFP:MAX,MIN ',ITNMX,JTNMX,TNMX, ITNMN,JTNMN,TNMN C// C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CC C*EJN TNMX= EJN(ML1+3) TNMN= EJN(ML1+3) ITNMX=2 JTNMX=3 ITNMN=2 JTNMN=3 DO 1681 I= 1, M2-1 DO 1681 J= 2, ML1 K= ML1*I+J WKTN= EJN(K) IF(WKTN.GT.TNMX) THEN TNMX=WKTN ITNMX=I+1 JTNMX=J ENDIF IF(WKTN.LT.TNMN) THEN TNMN=WKTN ITNMN=I+1 JTNMN=J ENDIF 1681 CONTINUE * WRITE(6,1010) '*EJN:MAX,MIN ',ITNMX,JTNMX,TNMX, ITNMN,JTNMN,TNMN CC C*EJP TNMX= EJP(ML1+3) TNMN= EJP(ML1+3) ITNMX=2 JTNMX=3 ITNMN=2 JTNMN=3 DO 2681 I= 1, M2-1 DO 2681 J= 2, ML1 K= ML1*I+J WKTN= EJP(K) IF(WKTN.GT.TNMX) THEN TNMX=WKTN ITNMX=I+1 JTNMX=J ENDIF IF(WKTN.LT.TNMN) THEN TNMN=WKTN ITNMN=I+1 JTNMN=J ENDIF 2681 CONTINUE * WRITE(6,1010) '*EJP:MAX,MIN ',ITNMX,JTNMX,TNMX, ITNMN,JTNMN,TNMN CC CC RETURN END * C C SUBROUTINE DYILST(M1, MJ,MJ1,MJM, MJ2,MJ3, LY0,LY1,LY2,LY3 + , JY0,JY1,JY2,JY3, MSU,MSM,MSD, DYI,DHY, M0Q) IMPLICIT REAL*8(A-H,O-Z) DIMENSION DYI(M1),DHY(M1) C JY0=LY0*MJ +1 JY1=JY0+(LY1-LY0)*MJ1 JY2=JY1+(LY2-LY1)*MJM JY3=JY2+(LY3-LY2)*MJ2 MSD=JY1 C MSM=JY1+MJM MSM=M0Q MSU=JY2 CC DYI(1)=0.0D0 DO 15 I=2,JY0 15 DYI(I)=DFLOAT(MJ) DO 16 I=JY0+1,JY1 16 DYI(I)=DFLOAT(MJ1) DO 17 I=JY1+1,JY2 17 DYI(I)=DFLOAT(MJM) DO 18 I=JY2+1,JY3 18 DYI(I)=DFLOAT(MJ2) DO 19 I=JY3+1,M1 19 DYI(I)=DFLOAT(MJ3) C DO 102 I=2,M1 102 DHY(I)=1.0D0/DYI(I) C RETURN END CCC CCC SUBROUTINE DXKLST(M2, MK,MK1,MKM, KXA,KXB,KXC,KX1, KS0,KS1,KSM + , DXK,DHX, KXBL,KXCR) IMPLICIT REAL*8(A-H,O-Z) DIMENSION DXK(M2+1),DHX(M2+1) C KXA=KS0*MK KXB=KXA+KS1*MK1 KXC=KXB+KSM*MKM KX1=KXC+KS1*MK1 C KXBL=KXB+MKM KXCR=KXC-MKM CC DO 5 I=1,KXA 5 DXK(I)=DFLOAT(MK) DO 6 I=KXA+1,KXB 6 DXK(I)=DFLOAT(MK1) DO 7 I=KXB+1,KXC 7 DXK(I)=DFLOAT(MKM) DO 8 I=KXC+1,KX1 8 DXK(I)=DFLOAT(MK1) C DO 9 I=KX1+1,M2+1 9 DXK(I)=DFLOAT(MK) C DO 101 I=1,M2+1 101 DHX(I)=1.0D0/DXK(I) C RETURN END CCC CCC SUBROUTINE OUTWR(M1,N, ZAN, IN,J1,J2,J3,J4,J5,J6 +, CNRP,INO,KCTO,ERRO, CNRLO,INLO,KCTLO +, VDRAIK,GID, BW,DN +, KXC,KX1,KXA,KXB, KXBL,KXCR +, IB, VDRAIN, INB,INLB,KCTB,KCTLB +, DP,CNRPO,INPO,KCTPO,INPB,KCTPB, GIDP +, GIDVG,GIDVGP,KCTOI,KCTLOI,KCTPOI,INI +, VBACKK,VBACKG, MSM) IMPLICIT REAL*8(A-H,O-Z) DIMENSION BW(-M1:N+M1),DN(N), DP(N) * WRITE(6,*) ' ' WRITE(6,*) ' -------------------- OUT LOOP --------------------' WRITE(6,12) IB,' *** VBACKK, VDRAIK = *** ',VBACKK,VDRAIK WRITE(6,*) ' VBACKG, VDRAIN = ',VBACKG,VDRAIN CC WRITE(6,11) ' ** ZAN = ** ',ZAN * IF(INI.NE.0) THEN CTM=FLOAT(KCTOI)/INI CTLM=FLOAT(KCTLOI)/INI CTPM=FLOAT(KCTPOI)/INI ELSE CTM=FLOAT(KCTOI)/INO CTLM=FLOAT(KCTLOI)/INO CTPM=FLOAT(KCTPOI)/INO END IF WRITE(6,5810) IN , BW(J1),BW(J2),BW(J3),BW(J4),BW(J5),BW(J6) + , CNRP,INO,KCTO, KCTOI,INI, CTM, ERRO WRITE(6,5820) IN , DN(J1),DN(J2),DN(J3),DN(J4),DN(J5),DN(J6) + , CNRLO,INLO,KCTLO,KCTLOI,INI,CTLM WRITE(6,5820) IN , DP(J1),DP(J2),DP(J3),DP(J4),DP(J5),DP(J6) + , CNRPO,INPO,KCTPO,KCTPOI,INI,CTPM C# C WRITE(6,*) ' *** GID = ***' +, GID, ' (AMPERE)' WRITE(6,*) ' *** GIDP = ***', GIDP, ' (AMPERE)' WRITE(6,*) ' GIDVG,GIDVGP= ', GIDVG,GIDVGP, ' (AMPERE)' CC WRITE(6,9100) BW(M1*(KXA-1)+MSM),BW(M1*(KXB-1)+MSM) +, BW(M1*(KXBL-1)+MSM),BW(M1*(KXCR-1)+MSM) +, BW(M1*(KXC-1)+MSM),BW(M1*(KX1-1)+MSM) WRITE(6,9100) DN(M1*(KXA-1)+MSM),DN(M1*(KXB-1)+MSM) +, DN(M1*(KXBL-1)+MSM),DN(M1*(KXCR-1)+MSM) +, DN(M1*(KXC-1)+MSM),DN(M1*(KX1-1)+MSM) WRITE(6,9100) DP(M1*(KXA-1)+MSM),DP(M1*(KXB-1)+MSM) +, DP(M1*(KXBL-1)+MSM),DP(M1*(KXCR-1)+MSM) +, DP(M1*(KXC-1)+MSM),DP(M1*(KX1-1)+MSM) CC WRITE(6,15) ' ** INB,KCTB// INLB,KCTLB **' +, INB,KCTB,INLB,KCTLB WRITE(6,15) ' ** INPB,KCTPB// **' +, INPB,KCTPB C# WRITE(6,*) ' -------------------OUT LOOP END ------------------' * 11 FORMAT(1H ,A ,D20.7,F10.5) 12 FORMAT(1H ,I4,A ,2F10.5) 15 FORMAT(1H ,A ,4I15) 9100 FORMAT(1H ,9D14.6) 5810 FORMAT(1H,I4,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,D11.2,I5,I6,I5 +, I5,F7.1, D12.3 ) 5820 FORMAT(1H,I4,D12.4,D12.4,D12.4,D12.4,D12.4,D12.4,D11.2,I5,I6,I5 +, I5,F7.1, D12.3 ) RETURN END C C C SUBROUTINE TOKOLB(N,M1,M2,IXS,NOXS, IYS,NOYS, F,MJM) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IXS(M1), IYS(M2), F(N) C OPEN(61,FILE='./DIOBW1.DATA') C WRITE(61,*) ' % BW WRITE ****' CC C DO 82 KK= NOYS,1,-1 K=IYS(KK) WRITE(6,2000) (F(M1*(IXS(J)-1)+K),J=1,NOXS),' ; ' 82 CONTINUE C CC --------------------------------------------- CC WRITE(61,*) ' z=[ ' C *C MATLAB You *** DO 80 KK= 1,NOYS K=IYS(KK) WRITE(61,2000) (F(M1*(IXS(J)-1)+K),J=1,NOXS),' ; ' 80 CONTINUE CCC C WRITE(61,*) ' ]; ' IF(MJM.EQ.1) THEN WRITE(61,*) ' % contour(0:1:9, 0:1:9, z,10);' WRITE(61,*) ' contour(0:1:9, 0:1:10, z,10);' ELSE WRITE(61,*) ' contour(0:1:9, 0:0.5:10, z,10);' ENDIF C 2000 FORMAT(1H ,10F7.3,A3) RETURN END CC CC SUBROUTINE IXSLIST(M1,MK,MK1,MKM,IXS,NOXS) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IXS(M1) C NOXS=10 IXS(1)=MK IXS(2)= IXS(1)+MK1 IXS(3)= IXS(2)+MK1 IXS(4)= IXS(3)+MK1 IXS(5)= IXS(4)+MKM IXS(6)= IXS(5)+MKM IXS(7)= IXS(6)+MKM IXS(8)= IXS(7)+MK1 IXS(9)= IXS(8)+MK1 IXS(10)=IXS(9)+MK1 C RETURN END CCC CCC SUBROUTINE IYSLIST(M2,MJ,MJ1,MJM,MJ2,MJ3,IYS,NOYS) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION IYS(M2) C IF(MJM.EQ.1) THEN NOYS=11 IYS(1)=1 IYS(2)= IYS(1)+MJ IYS(3)= IYS(2)+MJ1 IYS(4)= IYS(3)+MJ1 IYS(5)= IYS(4)+MJ1 IYS(6)= IYS(5)+MJM IYS(7)= IYS(6)+MJM IYS(8)= IYS(7)+MJ2 IYS(9)= IYS(8)+MJ2 IYS(10)= IYS(9)+MJ2 IYS(11)=IYS(10)+MJ3 ELSE CCCCCC NOYS=21 IYS(1)=1 IYS(2)= IYS(1)+MJ/2 IYS(3)= IYS(2)+MJ/2 IYS(4)= IYS(3)+MJ1/2 IYS(5)= IYS(4)+MJ1/2 IYS(6)= IYS(5)+MJ1/2 IYS(7)= IYS(6)+MJ1/2 IYS(8)= IYS(7)+MJ1/2 IYS(9)= IYS(8)+MJ1/2 IYS(10)= IYS(9)+MJM/2 IYS(11)= IYS(10)+MJM/2 IYS(12)= IYS(11)+MJM/2 IYS(13)= IYS(12)+MJM/2 IYS(14)= IYS(13)+MJ2/2 IYS(15)= IYS(14)+MJ2/2 IYS(16)= IYS(15)+MJ2/2 IYS(17)= IYS(16)+MJ2/2 IYS(18)= IYS(17)+MJ2/2 IYS(19)= IYS(18)+MJ2/2 IYS(20)=IYS(19)+MJ3/2 IYS(21)=IYS(20)+MJ3/2 ENDIF C C RETURN END CCC CCC *PROCESS NOOPLIST,NOSOURCE C////////////////////////////////////////////////// C ** **CC ** ** ** SUBROUTINE UPILU4(N,M,AA,AB,AC,AE,AF,UB,DI,U,AB1,AC1,AE1,AF1,UB1) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M:N),AC(-M:N),DI(-M:N), + AE(-M:N),AF(-M:N),UB(-M:N), + AB1(-M:N+M), AC1(-M:N+M),AE1(-M:N+M),AF1(-M:N+M),UB1(-M:N+M) *** DO 100 I = 1, N *** CC W= ABS(AC1(I))+ABS(AB1(I))+ABS(AB(I))+ABS(AC(I)) CC UUI = MAX(W,AA(I)) *** AE1(I)= -AC1(I)*UB(I-M)*DI(I-M) AF1(I)= -AE1(I)*UB(I-M+1)*DI(I-M+1) UB1(I)= AB1(I)-AC1(I)*AE(I-M)*DI(I-M)-AE1(I)*AF(I-M+1)*DI(I-M+1) * C/ U: UTOP DIV=U*AA(I)-UB1(I)*UB(I-1)*DI(I-1)-AC1(I)*AC(I-M)*DI(I-M) + -AE1(I)*AE(I-M+1)*DI(I-M+1)-AF1(I)*AF(I-M+2)*DI(I-M+2) C DIV= AA(I)-UB1(I)*UB(I-1)*DI(I-1)-AC1(I)*AC(I-M)*DI(I-M) C + -AE1(I)*AE(I-M+1)*DI(I-M+1)-AF1(I)*AF(I-M+2)*DI(I-M+2) DI(I)=1.0D0/DIV * UB(I)= AB(I)-AE1(I)*AC(I-M+1)*DI(I-M+1) + -AF1(I)*AE(I-M+2)*DI(I-M+2) AF(I)= -UB1(I)*AE(I-1)*DI(I-1) AE(I)= -UB1(I)*AC(I-1)*DI(I-1) * 100 CONTINUE RETURN END * CC CC&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& CC CCC SUBROUTINE JNCAL(NL,ML1,M2, WKL,DXK, DMOBN +, PSI,DNL +, DJN, DYI,M1F, EINSN, EJN, ENX,AENX, DJNH,DVN,DVNX +, ENY,AENY) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION WKL(NL+ML1),DXK(M2+1), DNL(-ML1:NL+ML1) +, PSI(-ML1:NL+ML1) +, DMOBN(NL+ML1), ENX(NL+ML1), AENX(NL+ML1) +, DJN(NL+ML1), DYI(M1F), EINSN(NL+ML1), EJN(NL+ML1) +, DJNH(NL+ML1) +, DVN(NL+ML1) +, DVNX(NL+ML1) +, ENY(NL+ML1), AENY(NL+ML1) CCC M2CN= (M2+1)/2 WEE= 1.6022D-19 CCC C// DO 99 I=1,NL+ML1 99 WKL(I)=0.0D0 ** DO 222 I= 2, M2 DO 223 J= 2,ML1 K= ML1*(I-1)+J * WKL(K) = DMOBN(K)/EINSN(K) * 223 CONTINUE 222 CONTINUE C// * * C// DO 200 J=2, M2 DO 200 J=1, M2+1 DO 210 K=2, ML1 I= ML1*(J-1)+K CCC IDH= I IF (J.GT.M2CN) THEN IDH= I-ML1 ENDIF CCC DMKLIH = 0.5D0/DXK(J) WPSI= PSI(IDH)-PSI(IDH-1) DYMI= WPSI*DYI(K) CCC IF (J.LE.M2CN) THEN CCC PT=-(DMOBN(I)*DYMI*DNL(I-1)-BER4(EINSN(I)*WPSI) + *WKL(I)*(DNL(I)-DNL(I-1))*DYI(K)) ELSE TQ=-(DMOBN(IDH+ML1)*DYMI*DNL(IDH-1)-BER4(EINSN(IDH+ML1) + *WPSI) + *WKL(IDH+ML1)*(DNL(IDH)-DNL(IDH-1))*DYI(K)) PT= TQ ENDIF CCC C* PT:Jy ** PTH=PT*(DMKLIH) CCC DHYMH= 0.5D0/DYI(K) FPSI= PSI(I)-PSI(I-ML1) DXMI= FPSI*DXK(J) SW=-(DMOBN(I)*DXMI*DNL(I-ML1)-BER4(EINSN(I)*FPSI) + *WKL(I)*(DNL(I)-DNL(I-ML1))*DXK(J)) C* SW:Jx ** SWH=SW*(DHYMH) C DJN(I)= DSQRT(PT*PT+SW*SW)*1.0D10 C// C* DVN: jn/e, DVNX: jnx(SW)/e *** DVN(I)=DJN(I)/DNL(I) DVNX(I)=SW*1.0D10 /DNL(I) CC * okashiibe (2004.06.02) DJNH(I)= DSQRT(PTH*PTH+SWH*SWH)*1.0D10 C## DJNH(I)=DSQRT(PT*PT+SW*SW)*1.0D10 *DMKLIH*DHYMH CC##################################### CC DJNI= DSQRT(PT*PT+SW*SW) IF(DJNI.GE.1.0D-56) THEN EJN(I)= DABS((-DXMI)*(-SW) + (-DYMI)*(-PT)) /DJNI */ ENX(I) ENY(I) *** * EJN=DABS(E dot Jn)/DABS(Jn) ELSE EJN(I)= 0.0D0 ENDIF CC/ AENX(I)= DABS(ENX(I)) AENY(I)= DABS(ENY(I)) * AENX:GR keisan ni tukau, AENY tsuika * CC CC##################################### CC 210 CONTINUE 200 CONTINUE C RETURN END * ** SUBROUTINE JPCAL(NL,ML1,M2, WKL,DXK, DMOBN +, PSI,DNL +, DJN, DYI,M1F, EINSN, EJN) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION WKL(NL+ML1),DXK(M2+1), DNL(-ML1:NL+ML1) +, PSI(-ML1:NL+ML1) +, DMOBN(NL+ML1) +, DJN(NL+ML1), DYI(M1F), EINSN(NL+ML1), EJN(NL+ML1) * CCC M2CN= (M2+1)/2 CCC C// DO 99 I=1,NL+ML1 99 WKL(I)=0.0D0 ** DO 222 I= 2, M2 DO 223 J= 2,ML1 K= ML1*(I-1)+J * WKL(K) = DMOBN(K)/EINSN(K) * 223 CONTINUE 222 CONTINUE C// * * C// DO 200 J=2, M2 DO 200 J=1, M2+1 DO 210 K=2, ML1 I= ML1*(J-1)+K CCC IDH= I IF (J.GT.M2CN) THEN IDH= I-ML1 ENDIF CCC CCCC DMKLIH = 0.5D0/DXK(J) WPSI= PSI(IDH)-PSI(IDH-1) DYMI= WPSI*DYI(K) CCC IF (J.LE.M2CN) THEN CCC PT=-( -DMOBN(I)*DYMI*DNL(I-1)-BER4(-EINSN(I)*WPSI) + *WKL(I)*(DNL(I)-DNL(I-1))*DYI(K)) ELSE TQ=-( -DMOBN(IDH+ML1)*DYMI*DNL(IDH-1)-BER4(-EINSN(IDH+ML1) + *WPSI) + *WKL(IDH+ML1)*(DNL(IDH)-DNL(IDH-1))*DYI(K)) PT=TQ ENDIF CCC FPSI= PSI(I)-PSI(I-ML1) DXMI= FPSI*DXK(J) SW=-( -DMOBN(I)*DXMI*DNL(I-ML1)-BER4(-EINSN(I)*FPSI) + *WKL(I)*(DNL(I)-DNL(I-ML1))*DXK(J)) C DJN(I)= DSQRT(PT*PT+SW*SW)*1.0D10 CC##################################### CC DJNI= DSQRT(PT*PT+SW*SW) IF(DJNI.GE.1.0D-56) THEN CC EJN(I)= DABS((-DXMI)*SW + (-DYMI)*PT) /DJNI EJN(I)= DABS((-DXMI)*(-SW) + (-DYMI)*(-PT)) /DJNI */ EJP ENX(I) ENY(I) *** * EJP=DABS(E dot Jp)/DABS(Jp) ELSE EJN(I)= 0.0D0 ENDIF CC CC##################################### CC 210 CONTINUE 200 CONTINUE C RETURN END * ** SUBROUTINE GRSET(NL,ML1, M2, M1 +, GR,QNAMN,QNDMN,N, SINIM,DP,DN,DHX,DHY, GRWK +, DJN,DJP,EN, ENP , GRU) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION GR(NL+ML1), QNAMN(N+M1),QNDMN(N+M1), DP(N),DN(N) +, DHX(M2+1),DHY(M1), GRWK(NL+ML1) CCCC +, GR1(NL+ML1),GR2(NL+ML1),GR3(NL+ML1) +, DJN(NL+ML1),DJP(NL+ML1), EN(NL+ML1), ENP(NL+ML1) +, GRU(NL+ML1) * CC DO 11 I=1,NL+ML1 GR(I)=0.0D0 GRWK(I)=0.0D0 11 GRU(I)=0.0D0 * C// CMM=1.0D-10 CCC * SRH TERM *** DO 301 I=0, M2-1 * DO 302 J=2, ML1-1 K=M1*I +J KL=ML1*I+J CC KX=I+1 KY=J C/ TNA=3.94D-4/(1.0D0+(QNAMN(K+1)+QNDMN(K+1))/7.1D0) TNB=3.94D-4/(1.0D0+(QNAMN(K+M1+1)+QNDMN(K+M1+1))/7.1D0) TNC=3.94D-4/(1.0D0+(QNAMN(K)+QNDMN(K))/7.1D0) TND=3.94D-4/(1.0D0+(QNAMN(K+M1)+QNDMN(K+M1))/7.1D0) TPA=3.94D-5/(1.0D0+(QNAMN(K+1)+QNDMN(K+1))/7.1D0) TPB=3.94D-5/(1.0D0+(QNAMN(K+M1+1)+QNDMN(K+M1+1))/7.1D0) TPC=3.94D-5/(1.0D0+(QNAMN(K)+QNDMN(K))/7.1D0) TPD=3.94D-5/(1.0D0+(QNAMN(K+M1)+QNDMN(K+M1))/7.1D0) C/ FA1=((SINIM*SINIM-DP(K)*DN(K))/ + (TNA*(DP(K)+SINIM)+TPA*(DN(K)+SINIM))) + *DHX(KX)*DHY(KY+1)*0.25D0 FB1=((SINIM*SINIM-DP(K)*DN(K))/ + (TNB*(DP(K)+SINIM)+TPB*(DN(K)+SINIM))) + *DHX(KX+1)*DHY(KY+1)*0.25D0 FC1=((SINIM*SINIM-DP(K)*DN(K))/ + (TNC*(DP(K)+SINIM)+TPC*(DN(K)+SINIM))) + *DHX(KX)*DHY(KY)*0.25D0 C# GRU(KL)=((SINIM*SINIM-DP(K)*DN(K))/ + (TNC*(DP(K)+SINIM)+TPC*(DN(K)+SINIM))) FD1=((SINIM*SINIM-DP(K)*DN(K))/ + (TND*(DP(K)+SINIM)+TPD*(DN(K)+SINIM))) + *DHX(KX+1)*DHY(KY)*0.25D0 C/ GR(KL)=CMM*(FA1+FB1+FC1+FD1) GR1KL=(FA1+FB1+FC1+FD1)/ + (0.25D0*(DHX(KX)+DHX(KX+1))*(DHY(KY)+DHY(KY+1))) GRWK(KL)=GR1KL 302 CONTINUE C# 301 CONTINUE CCC C// * AUGER TERM *** CND=0.28D0 CPD=0.099D0 CCC DO 401 I=0, M2-1 * DO 402 J=2, ML1-1 K=M1*I +J KL=ML1*I+J CC KX=I+1 KY=J C/ FF=(SINIM*SINIM-DP(K)*DN(K))*(CND*DN(K)+CPD*DP(K)) C# GRU(KL)= GRU(KL)+ FF C/ GR(KL)= GR(KL)+ CMM*FF* + (0.25D0*(DHX(KX)+DHX(KX+1))*(DHY(KY)+DHY(KY+1))) GRWK(KL)= GRWK(KL)+ FF 402 CONTINUE C# 401 CONTINUE CCC C// * IMPACT IONIZATION TERM *** AN=7.0D0 BN=12.3D0 AP=15.88D0 BP=20.36D0 * DO 501 I=0, M2-1 C DO 502 J=2, ML1-1 K=M1*I +J KL=ML1*I+J CC KX=I+1 KY=J C/ IF(I.EQ.0) THEN FA1=0.0D0 FC1=0.0D0 ELSE CCC IF(DABS(EN(KL+1)).GE.1.0D-56 .AND. + DABS(ENP(KL+1)).GE.1.0D-56) THEN FA1=(DJN(KL+1)*AN*DEXP(-BN/EN(KL+1)) + + DJP(KL+1)*AP*DEXP(-BP/ENP(KL+1))) + *DHX(KX)*DHY(KY+1)*0.25D0 ELSE FA1=0.0D0 ENDIF CC IF(DABS(EN(KL)).GE.1.0D-56 .AND. + DABS(ENP(KL)).GE.1.0D-56) THEN FC1=(DJN(KL)*AN*DEXP(-BN/EN(KL)) + + DJP(KL)*AP*DEXP(-BP/ENP(KL))) + *DHX(KX)*DHY(KY)*0.25D0 C# GRU(KL)= GRU(KL)+ + (DJN(KL)*AN*DEXP(-BN/EN(KL)) + + DJP(KL)*AP*DEXP(-BP/ENP(KL))) ELSE FC1=0.0D0 ENDIF CC ENDIF CCC IF(I.EQ.M2-1) THEN FB1=0.0D0 FD1=0.0D0 ELSE CC IF(DABS(EN(KL+ML1+1)).GE.1.0D-56 .AND. + DABS(ENP(KL+ML1+1)).GE.1.0D-56) THEN FB1=(DJN(KL+ML1+1)*AN*DEXP(-BN/EN(KL+ML1+1)) + + DJP(KL+ML1+1)*AP*DEXP(-BP/ENP(KL+ML1+1))) + *DHX(KX+1)*DHY(KY+1)*0.25D0 ELSE FB1=0.0D0 ENDIF CC IF(DABS(EN(KL+ML1)).GE.1.0D-56 .AND. + DABS(ENP(KL+ML1)).GE.1.0D-56) THEN FD1=(DJN(KL+ML1)*AN*DEXP(-BN/EN(KL+ML1)) + + DJP(KL+ML1)*AP*DEXP(-BP/ENP(KL+ML1))) + *DHX(KX+1)*DHY(KY)*0.25D0 ELSE FD1=0.0D0 ENDIF CCC ENDIF C/ GR(KL)=GR(KL)+CMM*(FA1+FB1+FC1+FD1) GR3KL=(FA1+FB1+FC1+FD1)/ + (0.25D0*(DHX(KX)+DHX(KX+1))*(DHY(KY)+DHY(KY+1))) GRWK(KL)=GRWK(KL)+GR3KL 502 CONTINUE C# * 501 CONTINUE C C------------------------------------------------- DO 444 II=1, M2 KKL=ML1*(II-1)+1 GR(KKL)=0.0D0 444 GRWK(KKL)=0.0D0 DO 445 II=1, M2 KKL=ML1*(II-1)+ML1 GR(KKL)=0.0D0 445 GRWK(KKL)=0.0D0 C------------------------------------------------- C# DO 32 I= 1, M2+1 32 GRU(ML1*(I-1)+1)= 0.0D0 CC DO 20 I = 1, ML1 20 GRU(I) = 0.0D0 CC * nakutemo yoi *** DO 21 I = NL+1, NL+ML1 21 GRU(I) = 0.0D0 * CCC RETURN END *PROCESS NOOPLIST,NOSOURCE * CC SUBROUTINE MKEINS(NL,ML1,M2, EINS, EINSN,EINSP) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION EINSN(NL+ML1), EINSP(NL+ML1) * DO 100 J=1, M2+1 DO 110 K=2, ML1 I= ML1*(J-1)+K CC EINSN(I)= EINS EINSP(I)= EINS C 110 CONTINUE 100 CONTINUE CC DO 32 I= 1, M2+1 EINSN(ML1*(I-1)+1)= 0.0D0 32 EINSP(ML1*(I-1)+1)= 0.0D0 * RETURN END * CC SUBROUTINE ENCAL(NL,ML1,M2,PSI, EN, DXK,DHY,M1, ENY + ,DNL,DPL, DFN,DFP + ,EINSN,EINSP, GAMA, ENX) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION PSI(-ML1:NL+ML1), EN(NL+ML1), DXK(M2+1),DHY(M1) + , ENY(NL+ML1), ENX(NL+ML1) + , DNL(-ML1:NL+ML1), DFN(NL+ML1) + , DPL(-ML1:NL+ML1), DFP(NL+ML1) + , EINSN(NL+ML1),EINSP(NL+ML1) * CCC M2CN= (M2+1)/2 CCC C// DO 100 J=2, M2 DO 100 J=1, M2+1 DO 110 K=2, ML1 I= ML1*(J-1)+K CCC IDH= I IF (J.GT.M2CN) THEN IDH= I-ML1 ENDIF CCC C GRADM1=PSI(I-1)-PSI(I) GRADM1=PSI(IDH-1)-PSI(IDH) GRADMM=PSI(I-ML1)-PSI(I) C// GRX=GRADMM/DHX(J) GRX=GRADMM*DXK(J) GRY=GRADM1/DHY(K) EN(I)= DSQRT(GRX*GRX+GRY*GRY) CC *** CCC ENY:-DPSIY/DY, ENX:-DPSIX/DX CCC CCC GRADR1=PSI(I-1)-PSI(I) GRYR=GRADR1/DHY(K) CCC C ENY(I)= (GRY) ENY(I)= (GRYR) ENX(I)= (GRX) C C** DRIVING FORCE SET (N) C/// DFNX=0.0D0 DFNY=0.0D0 IF(J.NE.1 .AND. J.NE.M2+1) THEN DDNX=(DNL(I)-DNL(I-ML1))*DXK(J) CCC DDNY=(DNL(I)-DNL(I-1))/DHY(K) DDNY=(DNL(IDH)-DNL(IDH-1))/DHY(K) C* DFNX=-GRX -GAMA*(2.0D0*DDNX)/((DNL(I)+DNL(I-ML1))*EINSN(I)) CCC DFNY=-GRY -GAMA*(2.0D0*DDNY)/((DNL(I)+DNL(I-1))*EINSN(I)) DFNY=-GRY -GAMA*(2.0D0*DDNY)/((DNL(IDH)+DNL(IDH-1))*EINSN(I)) ENDIF CC DFN(I)=DSQRT(DFNX*DFNX+ DFNY*DFNY) **// C** DRIVING FORCE SET (P) C/// DFNX=0.0D0 DFNY=0.0D0 IF(J.NE.1 .AND. J.NE.M2+1) THEN DDNX=(DPL(I)-DPL(I-ML1))*DXK(J) CCC DDNY=(DPL(I)-DPL(I-1))/DHY(K) DDNY=(DPL(IDH)-DPL(IDH-1))/DHY(K) C* DFNX=-GRX +GAMA *(2.0D0*DDNX)/((DPL(I)+DPL(I-ML1))*EINSP(I)) CCC DFNY=-GRY +GAMA*(2.0D0*DDNY)/((DPL(I)+DPL(I-1))*EINSP(I)) DFNY=-GRY+GAMA*(2.0D0*DDNY)/((DPL(IDH)+DPL(IDH-1))*EINSP(I)) ENDIF CC DFP(I)=DSQRT(DFNX*DFNX+ DFNY*DFNY) **// C 110 CONTINUE 100 CONTINUE C C------------------------------------------------------- DO 400 J=1, M2+1 K=1 I= ML1*(J-1)+K CC GRADM1=PSI(I-1)-PSI(I) GRADMM=PSI(I-ML1)-PSI(I) C// GRX=GRADMM/DHX(J) GRX=GRADMM*DXK(J) CC ENX(I)= (GRX) C 400 CONTINUE C C------------------------------------------------------- CC DO 20 I = 1, ML1 EN(I) = 0.0D0 ENY(I) = 0.0D0 ENX(I) = 0.0D0 DFN(I) = 0.0D0 20 DFP(I) = 0.0D0 DO 21 I = NL+1, NL+ML1 EN(I) = 0.0D0 ENY(I) = 0.0D0 ENX(I) = 0.0D0 DFN(I) = 0.0D0 21 DFP(I) = 0.0D0 C// DO 32 I= 1, M2+1 EN(ML1*(I-1)+1)= 0.0D0 ENY(ML1*(I-1)+1)= 0.0D0 DFN(ML1*(I-1)+1)= 0.0D0 32 DFP(ML1*(I-1)+1)= 0.0D0 * C/*** CC -------------- TOSIKOSI ---- * OHMIC CONTACT *** DO 403 I=1, M2-1 J=2 DFN(ML1*I+J)=EN(ML1*I+J) 403 DFP(ML1*I+J)=EN(ML1*I+J) * CC C ** BACK GATE **** DO 555 J=1, M2-1 C I= ML1*J +ML1 DFN(I)=EN(I) DFP(I)=EN(I) 555 CONTINUE C///////////////////////////////////// CCCC RETURN END * ** SUBROUTINE MOBCAL(NL,ML1,M2, DMOBN, QNAMN,QNDMN,N,M1, ENY +, DFN, EINSN,GAMA,IUTN,EINS +, YREFCM, BET0MB, FVNMB, IMOBTN,TWN10J +, FCTUTN) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DMOBN(NL+ML1), QNAMN(N+M1),QNDMN(N+M1) + , ENY(NL+ML1) + , DFN(NL+ML1), EINSN(NL+ML1) * CCCCCCCCCCCCCCCC C IMOBC= 1 C 1:universality , 0:mae *** CCCCCCCCCCCCCCCC C M2CN= (M2+1)/2 CCC VNSAT=90.0D0 C*** CC ECRIT=6.49D4 ECRIT=0.649D0 C*** C// DO 100 J=2, M2 DO 100 J=1, M2+1 DO 110 K=2, ML1 I= ML1*(J-1)+K CCC IDH= I IF (J.GT.M2CN) THEN IDH= I-ML1 ENDIF CCC KP= M1*(J-1)+K CI= QNAMN(KP)+QNDMN(KP) C** CI=CI*1.0D15 *** RMUNLI= 80.0D0+ 1350.0D0/(1.0D0+(CI/1.12D2)**0.72D0) CC C*** MU-LIS *** C##################### YREF= YREFCM*1.0D5 CC// C Y=DFLOAT(K-1)/DMJF0 C FY=2.0D0*DEXP(-(Y/YREF)**2)/(1.0D0+DEXP(-2.0D0*(Y/YREF)**2)) C FY=0.0D0 * FY=1: hyoumen CC// CCCCCCCCCCCCCCCCCCCCCCCCCCC C IF(IMOBC.EQ.1) THEN CCCCCCCCCCCCCCCCCCCCCCCCCCC C RMUREF=850.0D0 BET0= BET0MB C/ CC RMUNLI=(RMUREF+(RMUNLI-RMUREF)*(1.0D0-FY))/ * DSQRT(1.0D0+FY*BET0* DABS(ENY(IDH)) /ECRIT) CC C*** C IF(IMOBTN.EQ.1) THEN * IMOBTN:0 moto C Cmoto TWN=0.2D-2 C TWN=0.4D-2 TWN= TWN10J UUT0= 2.58510D-2 C IF(J.EQ.1 .OR. J.EQ.M2+1) THEN UUTN= UUT0 ELSE CC1 UUTN= 1.0D0/EINSN(IDH) UUTN= 1.0D0/EINSN(I) ENDIF C CCC DDUT=UUTN-UUT0 CCC CC WTN= 1.5D0*RMUNLI*(UUTN-UUT0)/ (TWN*VNSAT*VNSAT) WTN= 1.5D0*RMUNLI* DDUT/ (TWN*VNSAT*VNSAT) C*** CC DMLIE = RMUNLI/ (1.0D0+WTN) DMLIE = RMUNLI/ (1.0D0+FCTUTN*WTN) CC CC## CCCC ELSE W= 2.0D0*RMUNLI* DFN(I) /VNSAT C*** CC DMLIE = RMUNLI*2.0D0/(1.0D0+DSQRT(1.0D0+W*W)) DMLIE = RMUNLI*2.0D0/(1.0D0+DSQRT(1.0D0+ FVNMB*W*W)) ENDIF C/ DMLISY = DMLIE C CCCCCCCCCCCC C ENDIF CCCCCCCCCCCC C DMOBN(I) = DMLISY C// RMULIS=RMUNLI IF(IUTN.EQ.1) THEN C/ EINSN(I)=EINS/ (1.0D0+GAMA*(RMULIS/DMLIE -1.0D0)) EINSN(I)=EINS/ (1.0D0+GAMA*(RMULIS/DMLISY -1.0D0)) ENDIF C// 110 CONTINUE 100 CONTINUE CC DO 20 I = 1, ML1 DMOBN(I) = 0.0D0 20 EINSN(I) = 0.0D0 DO 21 I = NL+1, NL+ML1 DMOBN(I) = 0.0D0 21 EINSN(I) = 0.0D0 C// DO 32 I= 1, M2+1 DMOBN(ML1*(I-1)+1)= 0.0D0 32 EINSN(ML1*(I-1)+1)= 0.0D0 * C/*** CC -------------- TOSIKOSI ---- * OHMIC CONTACT *** DO 403 J=1, M2-1 403 EINSN(ML1*J+2)=EINS * CC C ** BACK GATE **** DO 505 J=1, M2-1 * - I= ML1*J+ ML1 EINSN(I)=EINS 505 CONTINUE C=================================== CCC * RETURN END * ** SUBROUTINE MOBCAP(NL,ML1,M2, DMOBN, QNAMN,QNDMN,N,M1, ENY +, DFP, EINSN,GAMA,IUTN,EINS +, YREFCM, BET0MB, FVNMB) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DMOBN(NL+ML1), QNAMN(N+M1),QNDMN(N+M1) + , ENY(NL+ML1) + , DFP(NL+ML1), EINSN(NL+ML1) * CCCCCCCCCCCCCCC C IMOBC= 1 C 1:universality , 0:mae *** CCCCCCCCCCCCCCC C M2CN= (M2+1)/2 CCC CC C ECRIT=1.87D4 ECRIT=0.187D0 C VPSAT=79.82D0 C C// DO 100 J=2, M2 DO 100 J=1, M2+1 DO 110 K=2, ML1 I= ML1*(J-1)+K CCC IDH= I IF (J.GT.M2CN) THEN IDH= I-ML1 ENDIF CCC KP= M1*(J-1)+K CI= QNAMN(KP)+QNDMN(KP) C** CI=CI*1.0D15 *** RMUNLI= 45.0D0+ 415.0D0/(1.0D0+(CI/2.23D2)**0.72D0) CC C*** MU-LIS *** C################################################# YREF= YREFCM*1.0D5 CC// C Y=DFLOAT(K-1)/DMJF0 C FY=2.0D0*DEXP(-(Y/YREF)**2)/(1.0D0+DEXP(-2.0D0*(Y/YREF)**2)) C FY=0.0D0 * FY=1: hyoumen CC// CCCCCCCCCCCCCCCCCCCCCCCCCC C IF(IMOBC.EQ.1) THEN CCCCCCCCCCCCCCCCCCCCCCCCCC C RMUREF=270.0D0 BET0= BET0MB C/ RMUNLI=(RMUREF+(RMUNLI-RMUREF)*(1.0D0-FY))/ * DSQRT(1.0D0+FY*BET0* DABS(ENY(IDH)) /ECRIT) C*** DMLIE = RMUNLI/(1.0D0+ FVNMB* RMUNLI* DFP(I) /VPSAT) C/ DMLISY = DMLIE C CCCCCCCCCCCC C ENDIF CCCCCCCCCCCC C DMOBN(I) = DMLISY C// RMULIS=RMUNLI IF(IUTN.EQ.1) THEN C/ EINSN(I)=EINS/ (1.0D0+GAMA*(RMULIS/DMLIE -1.0D0)) EINSNI =EINS/ (1.0D0+GAMA*(RMULIS/DMLISY -1.0D0)) EINSN(I)= MAX(2.9025D0, EINSNI) C**B EINSN(I)=EINS/ (1.0D0+GAMA*(RMULIS/DMLISY -1.0D0)) ENDIF C// 110 CONTINUE 100 CONTINUE CC DO 20 I = 1, ML1 DMOBN(I) = 0.0D0 20 EINSN(I) = 0.0D0 DO 21 I = NL+1, NL+ML1 DMOBN(I) = 0.0D0 21 EINSN(I) = 0.0D0 C// DO 32 I= 1, M2+1 DMOBN(ML1*(I-1)+1)= 0.0D0 32 EINSN(ML1*(I-1)+1)= 0.0D0 * C/*** CC -------------- TOSIKOSI ---- * OHMIC CONTACT *** DO 403 J=1, M2-1 403 EINSN(ML1*J+2)=EINS * CC C C ** BACK GATE **** DO 505 J=1, M2-1 * - I= ML1*J+ML1 EINSN(I)=EINS 505 CONTINUE C=================================== CCC * RETURN END * ** SUBROUTINE VBGETA(N,M1,M2, VBACKK, VBPRE, QNAH,BW, VOMCP) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION QNAH(N), BW(-M1:N+M1) *** VBSA=VBACKK-VBPRE CC DO 280 J=0, M2-1 DO 280 I=1, M1 KP= M1*J+ I IF(QNAH(KP).NE.0.0D0) THEN IF(DABS(VBSA).GT.0.0005D0) THEN BW(KP)=BW(KP)+VBSA ELSE BW(KP)=VBACKK-VOMCP ENDIF ENDIF 280 CONTINUE RETURN END * ** SUBROUTINE VDGETA(N,M1,M2, VDRAIK, VDPRE,QNDH,BW, VOMCN) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION QNDH(N), BW(-M1:N+M1) *** VDSA=VDRAIK-VDPRE CC DO 280 J=0, M2-1 DO 280 I=1, M1 KP= M1*J+ I IF(QNDH(KP).NE.0.0D0) THEN IF(DABS(VDSA).GT.0.0005D0) THEN BW(KP)=BW(KP)+VDSA ELSE BW(KP)=VOMCN+VDRAIK ENDIF ENDIF 280 CONTINUE RETURN END CCC CCC SUBROUTINE BIASP(VG1,VG2,VGSP, VD1,VD2,VDSP, VGB,VDB,N1DB +, NOB) IMPLICIT REAL*8(A-H,O-Z) DIMENSION VGB(N1DB), VDB(N1DB) * IF(VDSP.NE.0.0D0) THEN VGBW=VG1 NOB=0 DO 100 VI=VD1, VD2+0.000001D0, VDSP NOB=NOB+1 VDB(NOB)=VI C 100 VGB(NOB)=VGBW C/ IF(DABS(VD2-VDB(NOB)) .GT. 0.001D0) THEN NOB=NOB+1 VDB(NOB)=VD2 VGB(NOB)=VGBW ENDIF C/ ** ELSE IF(VGSP.NE.0.0D0) THEN NOB=0 DO 110 VI=VG1, VG2+0.000001D0, VGSP NOB=NOB+1 VGB(NOB)=VI 110 VDB(NOB)=VD1 C/ IF(DABS(VG2-VGB(NOB)) .GT. 0.001D0) THEN NOB=NOB+1 VGB(NOB)=VG2 VDB(NOB)=VD1 ENDIF C/ ** ELSE NOB=1 VDB(NOB)=VD1 C VGB(NOB)=VG1 ENDIF RETURN END CCC CCC SUBROUTINE SYOKIB(N,M2,ML1,NL, DNL +, SINIM, PSI,EINS, QND,QNA, DPL, DN,DP) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION DNL(-ML1:NL+ML1), DPL(-ML1:NL+ML1) +, PSI(-ML1:NL+ML1), QND(NL),QNA(NL), DN(N),DP(N) * SINIM2= SINIM*SINIM * DO 100 I=1, M2 DO 100 J=1,ML1 KL=ML1*(I-1)+ J CCC IF(QND(KL).NE.0.0D0 .AND. QNA(KL).NE.0.0D0) THEN DPL(KL)= SINIM2/QND(KL) DNL(KL)= QND(KL) ENDIF CCC IF(QNA(KL).NE.0.0D0 .AND. QND(KL).EQ.0.0D0) THEN DPL(KL)= QNA(KL) DNL(KL)= SINIM2/QNA(KL) ENDIF * IF(QND(KL).NE.0.0D0 .AND. QNA(KL).EQ.0.0D0) THEN DPL(KL)= SINIM2/QND(KL) DNL(KL)= QND(KL) ENDIF 100 CONTINUE C * PSI WO TUKURU *** DO 285 I=1, M2 DO 285 J=1,ML1 K= ML1*(I-1)+J 285 PSI(K)= DLOG(DNL(K)/SINIM)/EINS * C DO 511 I=1, M2 DO 511 J=1,ML1 K= ML1*(I-1)+J DN(K)=DNL(K) 511 DP(K)=DPL(K) CCC * RETURN END * *C SUBROUTINE NINPUC(NL,ML1,M2, DPL,PSI +, EINS,SINIM, VBACKK,VDRAIK,QNDH, QNAH, DNDOCM, DNBACK) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION QNDH(NL), DPL(-ML1:NL+ML1),PSI(-ML1:NL+ML1), QNAH(NL) *** CC DO 280 J=0, M2-1 DO 280 I=1, ML1 K= ML1*J+I IF(QNDH(K).NE.0.0D0 .AND. QNAH(K).EQ.0.0D0) THEN DPL(K )= SINIM*DEXP( EINS*(PSI(K)-VDRAIK)) * VD<0:Jun bias ENDIF IF(QNDH(K).EQ.0.0D0 .AND. QNAH(K).NE.0.0D0) THEN DPL(K )= SINIM*DEXP( EINS*(PSI(K)-VBACKK)) ENDIF CCC IF(QNDH(K).NE.0.0D0 .AND. QNAH(K).NE.0.0D0) THEN DPL(K )= SINIM*DEXP( EINS*(PSI(K)-VDRAIK)) * VD<0:Jun bias ENDIF CCC/ IF(DPL(K ).LT.1.0D-56) THEN DPL(K )=1.0D-56 ENDIF 280 CONTINUE C =========================================== C * DO 802 J=0, M2-1 K=ML1*J+1 802 DPL(K)= DNDOCM * DO 804 J=0, M2-1 K=ML1*J+ML1 804 DPL(K)= DNBACK CC RETURN END * *C SUBROUTINE PINPUT(NL,ML1,M2, DPL,PSI +, EINS,SINIM, VBACKK,VDRAIK,QNDH, QNAH, DNDOCM, DNBACK) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION QNDH(NL), DPL(-ML1:NL+ML1),PSI(-ML1:NL+ML1), QNAH(NL) *** CC DO 280 J=0, M2-1 DO 280 I=1, ML1 K= ML1*J+I IF(QNDH(K).NE.0.0D0 .AND. QNAH(K).EQ.0.0D0) THEN DPL(K )= SINIM*DEXP(EINS*(VDRAIK-PSI(K))) * VD<0:Jun bias ENDIF IF(QNDH(K).EQ.0.0D0 .AND. QNAH(K).NE.0.0D0) THEN DPL(K )= SINIM*DEXP(EINS*(VBACKK-PSI(K))) ENDIF CCC IF(QNDH(K).NE.0.0D0 .AND. QNAH(K).NE.0.0D0) THEN DPL(K )= SINIM*DEXP(EINS*(VDRAIK-PSI(K))) * VD<0:Jun bias ENDIF CCC/ IF(DPL(K ).LT.1.0D-56) THEN DPL(K )=1.0D-56 ENDIF 280 CONTINUE C =========================================== C * DO 802 J=0, M2-1 K=ML1*J+1 802 DPL(K)= DNDOCM * DO 804 J=0, M2-1 K=ML1*J+ML1 804 DPL(K)= DNBACK CC RETURN END * CCC ** SUBROUTINE IDINTB(NL,ML1,M2, WKL,DXK, GID,DMOBN +, DYM,PSI,DNL, XID +, DYI,M1F, EINSN, Z) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION WKL(NL+ML1),DXK(M2+1),DYM(NL), DNL(-ML1:NL+ML1) +, XID(M2), PSI(-ML1:NL+ML1) +, DMOBN(NL+ML1), DYI(M1F), EINSN(NL+ML1) * C// DO 99 I=1,NL+ML1 99 WKL(I)=0.0D0 ** DO 222 I= 2, M2 DO 223 J= 2,ML1 K= ML1*(I-1)+J * WKL(K) = DMOBN(K)/EINSN(K) * 223 CONTINUE 222 CONTINUE C// * * DMJ=DYI(2) DO 410 J=1, M2 K=ML1*(J-1)+2 410 DYM(K)= (PSI(K)- PSI(K-1))*DMJ C GID= 0.0D0 * SOURCE NO HIDARI KADO *** J=1 I= ML1*(J-1)+2 DMKRIH= 0.5D0/DXK(J+1) TQ= (DMOBN(I+ML1)*DYM(I)*DNL(I-1)-BER4(EINSN(I+ML1)*DYM(I)/DMJ) + *WKL(I+ML1)*(DNL(I)-DNL(I-1))*DMJ)*(-DMKRIH) GID= GID+TQ CC XID(J)=TQ * DO 200 J=2, M2-1 I= ML1*(J-1)+2 DMKLIH = 0.5D0/DXK(J) DMKRIH = 0.5D0/DXK(J+1) PT= (DMOBN(I)*DYM(I)*DNL(I-1)-BER4(EINSN(I)*DYM(I)/DMJ) + *WKL(I)*(DNL(I)-DNL(I-1))*DMJ)*(-DMKLIH) TQ= (DMOBN(I+ML1)*DYM(I)*DNL(I-1)-BER4(EINSN(I+ML1)*DYM(I)/DMJ) + *WKL(I+ML1)*(DNL(I)-DNL(I-1))*DMJ)*(-DMKRIH) GID= GID+PT+TQ CC XID(J)=PT+TQ 200 CONTINUE C CC * DRAIN NO MIGI KADO *** C C---------------------------------------------------------------- J= M2 I= ML1*(J-1)+2 DMKLIH= 0.5D0/DXK(J) PT= (DMOBN(I)*DYM(I)*DNL(I-1)-BER4(EINSN(I)*DYM(I)/DMJ) + *WKL(I)*(DNL(I)-DNL(I-1))*DMJ)*(-DMKLIH) GID= GID+PT CC XID(J)=PT C// C---------------------------------------------------------------- C// C * CCC * FAC= 1.6022D-19*1.0D10 *** FAC= 1.6022D-9 GID=-Z*GID*FAC C RETURN END * ** SUBROUTINE IDINTP(NL,ML1,M2, WKL,DXK, GID,DMOBN +, DYM,PSI,DNL, XID +, DYI,M1F, EINSN, Z) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION WKL(NL+ML1),DXK(M2+1),DYM(NL), DNL(-ML1:NL+ML1) +, XID(M2), PSI(-ML1:NL+ML1) +, DMOBN(NL+ML1), DYI(M1F), EINSN(NL+ML1) * C// DO 99 I=1,NL+ML1 99 WKL(I)=0.0D0 ** DO 222 I= 2, M2 DO 223 J= 2,ML1 K= ML1*(I-1)+J * WKL(K) = DMOBN(K)/EINSN(K) * 223 CONTINUE 222 CONTINUE C// * * DMJ=DYI(2) DO 410 J=1, M2 K=ML1*(J-1)+2 410 DYM(K)= (PSI(K)- PSI(K-1))*DMJ C GID= 0.0D0 * SOURCE NO HIDARI KADO *** J=1 I= ML1*(J-1)+2 DMKRIH= 0.5D0/DXK(J+1) TQ=(-DMOBN(I+ML1)*DYM(I)*DNL(I-1)-BER4(-EINSN(I+ML1)*DYM(I)/DMJ) + *WKL(I+ML1)*(DNL(I)-DNL(I-1))*DMJ)*(-DMKRIH) GID= GID+TQ CC XID(J)=TQ * DO 200 J=2, M2-1 I= ML1*(J-1)+2 DMKLIH = 0.5D0/DXK(J) DMKRIH = 0.5D0/DXK(J+1) PT=(-DMOBN(I)*DYM(I)*DNL(I-1)-BER4(-EINSN(I)*DYM(I)/DMJ) + *WKL(I)*(DNL(I)-DNL(I-1))*DMJ)*(-DMKLIH) TQ=(-DMOBN(I+ML1)*DYM(I)*DNL(I-1)-BER4(-EINSN(I+ML1)*DYM(I)/DMJ) + *WKL(I+ML1)*(DNL(I)-DNL(I-1))*DMJ)*(-DMKRIH) GID= GID+PT+TQ CC XID(J)=PT+TQ 200 CONTINUE C * DRAIN NO MIGI KADO *** CC C------------------------------------------------------------------ J= M2 I= ML1*(J-1)+2 DMKLIH= 0.5D0/DXK(J) PT= (-DMOBN(I)*DYM(I)*DNL(I-1)-BER4(-EINSN(I)*DYM(I)/DMJ) + *WKL(I)*(DNL(I)-DNL(I-1))*DMJ)*(-DMKLIH) GID= GID+PT CC XID(J)=PT C------------------------------------------------------------------ C// C * CCC * FAC= 1.6022D-19*1.0D10 *** FAC= 1.6022D-9 GID= Z*GID*FAC C RETURN END C ############### HAJIME ######################################### * ** * ** * ** SUBROUTINE IDVGTB(NL,ML1,M2, WKL,DXK, GID,DMOBN +, DYM,PSI,DNL +, DYI,M1F, EINSN, Z) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION WKL(NL+ML1),DXK(M2+1),DYM(NL), DNL(-ML1:NL+ML1) +, PSI(-ML1:NL+ML1) +, DMOBN(NL+ML1), DYI(M1F), EINSN(NL+ML1) * C// DO 99 I=1,NL+ML1 99 WKL(I)=0.0D0 ** DO 222 I= 2, M2 DO 223 J= 2,ML1 K= ML1*(I-1)+J * WKL(K) = DMOBN(K)/EINSN(K) * 223 CONTINUE 222 CONTINUE C// * * DMJ= DYI(ML1) DO 410 J=1, M2 K=ML1*(J-1)+ML1-1 410 DYM(K)= (PSI(K+1)- PSI(K))*DMJ C GID= 0.0D0 C * DO 200 J= 1, M2 I= ML1*(J-1)+ML1-1 DMKLIH = 0.5D0/DXK(J) DMKRIH = 0.5D0/DXK(J+1) PT=(DMOBN(I+1)*DYM(I)*DNL(I)-BER4(EINSN(I+1)*DYM(I)/DMJ) + *WKL(I+1)*(DNL(I+1)-DNL(I))*DMJ)*( DMKLIH) TQ=(DMOBN(I+ML1+1)*DYM(I)*DNL(I)-BER4(EINSN(I+ML1+1)*DYM(I)/DMJ) + *WKL(I+ML1+1)*(DNL(I+1)-DNL(I))*DMJ)*( DMKRIH) GID= GID+PT+TQ CC 200 CONTINUE C CC * CCC * FAC= 1.6022D-19*1.0D10 *** FAC= 1.6022D-9 GID=-Z*GID*FAC C RETURN END * * ** SUBROUTINE IDVGTP(NL,ML1,M2, WKL,DXK, GID,DMOBN +, DYM,PSI,DNL +, DYI,M1F, EINSN, Z) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION WKL(NL+ML1),DXK(M2+1),DYM(NL), DNL(-ML1:NL+ML1) +, PSI(-ML1:NL+ML1) +, DMOBN(NL+ML1), DYI(M1F), EINSN(NL+ML1) * C// DO 99 I=1,NL+ML1 99 WKL(I)=0.0D0 ** DO 222 I= 2, M2 DO 223 J= 2,ML1 K= ML1*(I-1)+J * WKL(K) = DMOBN(K)/EINSN(K) * 223 CONTINUE 222 CONTINUE C// * * DMJ= DYI(ML1) DO 410 J=1, M2 K=ML1*(J-1)+ML1-1 410 DYM(K)= (PSI(K+1)- PSI(K))*DMJ C GID= 0.0D0 C * DO 200 J= 1, M2 I= ML1*(J-1)+ML1-1 DMKLIH = 0.5D0/DXK(J) DMKRIH = 0.5D0/DXK(J+1) PT=(-DMOBN(I+1)*DYM(I)*DNL(I)-BER4(-EINSN(I+1)*DYM(I)/DMJ) + *WKL(I+1)*(DNL(I+1)-DNL(I))*DMJ)*( DMKLIH) TQ=(-DMOBN(I+ML1+1)*DYM(I)*DNL(I)-BER4(-EINSN(I+ML1+1)*DYM(I)/DMJ) + *WKL(I+ML1+1)*(DNL(I+1)-DNL(I))*DMJ)*( DMKRIH) GID= GID+PT+TQ CC 200 CONTINUE C C * CCC * FAC= 1.6022D-19*1.0D10 *** FAC= 1.6022D-9 GID= Z*GID*FAC C RETURN END * C ############### OWARI ######################################### * ** SUBROUTINE POISON(N,M1,B1,F,R, AA,AB,AC,DI,P,W,X,AP,UB,AE +, AB1,AC1,UB1,AE1,RRA,PA,RR, AAG, IN +, M2, WK, DXK,SI +, DN,DP,ML1, IGUM, VDRAIN +, VBACKG, DMJMKV,DMJPKV, J1,J2,J3,J4,J5,J6, IDBG +, BW, EPSN,EPSDIV, KCTO,INO +, QNDMN,QNAMN, WKR, DHX,DHY +, DYI,DIKL,DIKR,DKIL,DKIR, CNR, USHIRP +, EINSN,EINSP, NL) IMPLICIT REAL*8(A-H,O-Z) *** *** DIMENSION AAG(N), B1(N) +, F(N+M1), R(N),WK(N+M1+1), BW(-M1:N+M1) +, AA(N),AB(-M1:N),AC(-M1:N),DI(-M1:N),P(-M1:N+M1) +, W(-M1:N+M1),X(-M1:N+M1),AP(N),UB(-M1:N),AE(-M1:N) +, AB1(-M1:N+M1),AC1(-M1:N+M1),UB1(-M1:N+M1),AE1(-M1:N+M1) +, RRA(-M1:N+M1),PA(-M1:N+M1),RR(-M1:N+M1) +, DXK(M2+1),DMJMKV(M2+1),DMJPKV(M2+1) +, DP(N), DN(N) +, QNDMN(N+M1),QNAMN(N+M1) +, WKR(N+M1), DHX(M2+1),DHY(M1) +, DYI(M1),DIKL(M1),DIKR(M1),DKIL(M1),DKIR(M1) +, EINSN(NL+ML1),EINSP(NL+ML1) CC CC * DO 544 I=0,M2-1 544 BW(M1*I+1) = 0.0D0 * DO 545 I=1,M2 545 BW(M1*I) = 0.0D0 * CCC ** CALL ZEROCR(N,M1,B1,F,R, AA,AB,AC,DI,P,W,X,AP,UB,AE +, AB1,AC1,UB1,AE1,RRA,PA,RR, AAG) * * CALL REVMTR(N,M1,M2,AA,AB,AC,WK, DXK,SI +, AAG,DN,DP,ML1, IGUM, DHX,DHY +, DYI,DIKL,DIKR,DKIL,DKIR, EINSN,EINSP, NL) * CALL FINPUT(N,M1,F,WK,VDRAIN, M2 +, SI,VBACKG,DXK,DMJMKV,DMJPKV, IGUM, DP,DN, ML1, BW +, AC,AB,AAG, SNRM, QNDMN,QNAMN, WKR, DHX,DHY +, DYI, EINSN,EINSP, NL) * * BNRM=0.0 DO 110 I=1,N 110 BNRM=BNRM+F(I)*F(I) BNRM=DSQRT(BNRM) * CALL ICDCMP( N,M1,AA,AB,AC,AE,UB,DI, USHIRP) CC CALL ICDCMP3( N,M1,AA,AB,AC,AE,AF,UB,DI, USHIRP) * IF(IGUM.EQ.2) THEN DO 225 I=1,N 225 X(I)=0.0D0 ELSE DO 226 I=1,N 226 X(I)=BW(I) ENDIF * EPSICW=EPSN/EPSDIV KCT=0 IPMX= 4 DO 230 IP=1,IPMX * CALL ICCG12(AA,AB,AC,AE,UB,DI,P,R,W,F + ,X,AP,EPSICW,ERR,M1,N,KCOUNT,BNRM) CC CALL ICCG13(AA,AB,AC,AE,AF,UB,DI,P,R,W,F CC + ,X,AP,EPSICW,ERR,M1,N,KCOUNT,BNRM) * KCT=KCT+KCOUNT CNR=0.0D0 DO 116 I=1,N B1(I)=F(I)-(AC(I-M1)*X(I-M1)+AB(I-1)*X(I-1) + +AA(I)*X(I)+AB(I)*X(I+1)+AC(I)*X(I+M1)) 116 CNR=CNR+B1(I)*B1(I) CNR=DSQRT(CNR)/BNRM * IF(IDBG.GT.1) THEN WRITE(6,5800) IN , X(J1),X(J2),X(J3),X(J4),X(J5),X(J6) + , CNR,IP,KCT ENDIF IF(CNR.LT.EPSN) GO TO 250 EPSICW=EPSICW/2.0 * 230 CONTINUE 250 CONTINUE * IF(IGUM.EQ.2) THEN DO 219 I=1,N 219 BW(I)=BW(I)+X(I) ELSE DO 220 I=1,N 220 BW(I)=X(I) ENDIF KCTO=KCTO+KCT IF(IP.EQ.IPMX+1) IP=IP-1 INO=INO+IP C * *// DO 446 I=0,M2-1 446 BW(M1*I+1) = VDRAIN C DO 445 I=1,M2 445 BW(M1*I) = VBACKG * 5800 FORMAT(1H,I4,F12.6,F12.6,F12.6,F12.6,F12.6,F12.6,D20.7,I5,I6) RETURN END * ** SUBROUTINE REVMTR(N,M1,M2,A0,A1,A2,WK, DXK,SI +, AAG,DN,DP,ML1, IGUM, DHX,DHY +, DYI,DIKL,DIKR,DKIL,DKIR, EINSN,EINSP, NL) * IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A0(N),A1(-M1:N),A2(-M1:N),WK(N+M1+1) +, DXK(M2+1) +, AAG(N),DN(N),DP(N), DHX(M2+1),DHY(M1) +, DYI(M1),DIKL(M1),DIKR(M1),DKIL(M1),DKIR(M1) +, EINSN(NL+ML1),EINSP(NL+ML1) **** **** DO 10 I = 1, N+M1 10 WK(I) = SI * DO 51 I = 1, M1 51 WK(I) = 0.0D0 DO 52 I = N+1, N+M1 52 WK(I) = 0.0D0 CCC * * DO 100 J = 1, M2 * DO 40 I=2,M1 DIKL(I)=(DYI(I)/DXK(J))*0.5D0 DIKR(I)=(DYI(I)/DXK(J+1))*0.5D0 DKIL(I)=(DXK(J)/DYI(I))*0.5D0 40 DKIR(I)=(DXK(J+1)/DYI(I))*0.5D0 CC CC I=M1*(J-1)+1 A0(I)=1.00 A1(I)=0.00 A2(I)=0.00 * * DO 110 K= 2, M1-1 I=M1*(J-1)+K DJVL= DIKL(K) DJVR=DIKR(K) DKVL=DKIL(K) DKVR=DKIR(K) DJWL= DIKL(K+1) DJWR=DIKR(K+1) DKWL=DKIL(K+1) DKWR=DKIR(K+1) A0(I)= WK(I)*DJVL+WK(I+M1)*DJVR + + WK(I+1)*DJWL+WK(I+M1+1)*DJWR + + WK(I)*DKVL+WK(I+1)*DKWL + + WK(I+M1)*DKVR+WK(I+M1+1)*DKWR A1(I)=-(WK(I+1)*DJWL+WK(I+M1+1)*DJWR) 110 A2(I)=-(WK(I+M1)*DKVR+WK(I+M1+1)*DKWR) * * I=M1*J A0(I)=1.00 A1(I)=0.00 A2(I)=0.00 * 100 CONTINUE * CC DO 210 I=1, N 210 AAG(I)= A0(I) * CCC * IGUM=1 OR 2 START *** IF(IGUM.EQ.1 .OR. IGUM.EQ.2) THEN ** * * EPE0:= E/EPS0= 1.6D-19/8.85D-19 *** CCC EPE0= 0.1809D0 EPE0= 1.6022D0/8.854D0 C * * CCC * IROU=1 START *** CC DO 304 I=0, M2-1 * DO 305 J=2, ML1-1 K=M1*I +J C/ KIL=M1*I +J CC KX=I+1 KY=J FA1=(DP(K)*EINSP(KIL+1)+DN(K)*EINSN(KIL+1)) + *DHX(KX)*DHY(KY+1)*0.25D0 FB1=(DP(K)*EINSP(KIL+ML1+1)+DN(K)*EINSN(KIL+ML1+1)) + *DHX(KX+1)*DHY(KY+1)*0.25D0 FC1=(DP(K)*EINSP(KIL)+DN(K)*EINSN(KIL)) + *DHX(KX)*DHY(KY)*0.25D0 FD1=(DP(K)*EINSP(KIL+ML1)+DN(K)*EINSN(KIL+ML1)) + *DHX(KX+1)*DHY(KY)*0.25D0 305 A0(K)=A0(K)+EPE0* (FA1+FB1+FC1+FD1) * 304 CONTINUE * CCC * IROU=1 END *** * * *** ENDIF CCC * IGUM=1 OR 2 END *** ** A1(N)=0.0D0 DO 801 I = N -M1+1, N A2(I)=0.0D0 801 CONTINUE *** * VBACKG NO TAMENI *** DO 813 I= 0,M2-1 J=M1-1 813 A1(M1*I+J)=0.0D0 C// *** RETURN END * ** SUBROUTINE FINPUT(N,M1,F,WK,VDRAIN, M2 +, SI,VBACKG,DXK,DMJMKV,DMJPKV, IGUM, DP,DN, ML1, BW +, AC,AB,AAG, SNRM, QNDMN,QNAMN, WKR, DHX,DHY +, DYI, EINSN,EINSP, NL) ** IMPLICIT REAL*8 (A-H,O-Z) DIMENSION F(N+M1),WK(N+M1),DXK(M2+1),DMJMKV(M2+1) +, DMJPKV(M2+1) +, DP(N), DN(N), BW(-M1:N+M1) +, AC(-M1:N),AB(-M1:N),AAG(N) +, QNDMN(N+M1),QNAMN(N+M1) +, WKR(N+M1), DHX(M2+1),DHY(M1) +, DYI(M1) +, EINSN(NL+ML1),EINSP(NL+ML1) **** C DO 10 I=1,N+M1 10 WKR(I)=0.0D0 * * EPE0:= E/EPS0= 1.6D-19/8.85D-19 *** CCC EPE0= 0.1809D0 EPE0= 1.6022D0/8.854D0 C * * CCC * IROU=1 START *** CC DO 301 I=0, M2-1 * DO 302 J=2, ML1-1 K=M1*I +J CC KX=I+1 KY=J FA1=(DP(K)-DN(K)-QNAMN(K+1)+QNDMN(K+1)) + *DHX(KX)*DHY(KY+1)*0.25D0 FB1=(DP(K)-DN(K)-QNAMN(K+M1+1)+QNDMN(K+M1+1)) + *DHX(KX+1)*DHY(KY+1)*0.25D0 FC1=(DP(K)-DN(K)-QNAMN(K)+QNDMN(K)) + *DHX(KX)*DHY(KY)*0.25D0 FD1=(DP(K)-DN(K)-QNAMN(K+M1)+QNDMN(K+M1)) + *DHX(KX+1)*DHY(KY)*0.25D0 F(K)=EPE0*(FA1+FB1+FC1+FD1) 302 WKR(K)=(FA1+FB1+FC1+FD1)/ + (0.25D0*(DHX(KX)+DHX(KX+1))*(DHY(KY)+DHY(KY+1))) * 301 CONTINUE * CCC * IROU=1 END *** * * *** DO 620 I=0, M2-1 620 F(M1*I+1)= 0.0D0 DO 640 I=0, M2-1 640 F(M1*I+M1)= 0.0D0 * ** CCC SNRM= 0.0D0 DO 490 I=1,N 490 SNRM= SNRM+F(I)*F(I) SNRM= DSQRT(SNRM) * CCC * IGUM=1 START *** IF(IGUM.EQ.1) THEN * CCC * IROU=1 START *** CC DO 304 I=0, M2-1 * DO 305 J=2, ML1-1 K=M1*I +J C/ KIL=M1*I +J CC KX=I+1 KY=J FA1=(DP(K)*EINSP(KIL+1)+DN(K)*EINSN(KIL+1)) + *DHX(KX)*DHY(KY+1)*0.25D0 FB1=(DP(K)*EINSP(KIL+ML1+1)+DN(K)*EINSN(KIL+ML1+1)) + *DHX(KX+1)*DHY(KY+1)*0.25D0 FC1=(DP(K)*EINSP(KIL)+DN(K)*EINSN(KIL)) + *DHX(KX)*DHY(KY)*0.25D0 FD1=(DP(K)*EINSP(KIL+ML1)+DN(K)*EINSN(KIL+ML1)) + *DHX(KX+1)*DHY(KY)*0.25D0 305 F(K)=F(K)+BW(K)*EPE0* (FA1+FB1+FC1+FD1) * 304 CONTINUE * CCC * IROU=1 END *** * * DO 660 I=0, M2-1 660 F(M1*I+1)= 0.0D0 DO 670 I=0, M2-1 670 F(M1*I+M1)= 0.0D0 CC * ENDIF CCC * IGUM=1 END *** * CCC * IGUM=2 START *** IF(IGUM.EQ.2) THEN * DO 200 I=1,N+M1 200 WK(I)=0.0D0 * DO 460 I=1, N 460 WK(I)= -(AC(I-M1)*BW(I-M1)+AB(I-1)*BW(I-1)+AAG(I)*BW(I) + +AB(I)*BW(I+1)+AC(I)*BW(I+M1)) * C DO 802 I=0, M2-1 802 WK(M1*I+1)= 0.0D0 DO 813 I=0, M2-1 813 WK(M1*I+M1)= 0.0D0 CC * DO 450 I=1,N 450 F(I)= F(I)+WK(I) * ENDIF CCC * IGUM=2 END *** * CC **** DMDY2=DYI(2) DO 500 K=1,M2+1 500 DMJPKV(K)= DMDY2/DXK(K) * CCCC CCCC F(2)=F(2)+SI*0.5D0*DMJPKV(2)*VDRAIN DO 513 I=2, M2-1 513 F(M1*(I-1)+2)= + F(M1*(I-1)+2)+SI*0.5D0*(DMJPKV(I)+DMJPKV(I+1))*VDRAIN C// F(M1*(M2-1)+2)= + F(M1*(M2-1)+2)+SI*0.5D0*DMJPKV(M2)*VDRAIN C// CC **** *** BACKGATE NI - DEN I O ATAERU *** CC C// DMJM=DYI(M1) C// DO 600 I=1,M2+1 600 DMJMKV(I) = DMJM/DXK(I) * SIO=SI * C// C// F(M1*1-1)=F(M1*1-1)+0.5D0*SIO*DMJMKV(2)*VBACKG DO 601 I=2,M2-1 601 F(M1*I-1)=F(M1*I-1)+SIO*0.5D0*(DMJMKV(I)+DMJMKV(I+1))*VBACKG * F(M1*M2-1)=F(M1*M2-1)+0.5D0*SIO*DMJMKV(M2)*VBACKG C// * RETURN END * ** SUBROUTINE NINPUT(N,M1,M2, DNL,PSI, DNDOCM +, B1,F,R,WK, AA,AB,AC,DI,P,W,X,AP,UB,AE +, AB1,AC1,UB1,AE1,RRA,PA,RR,DXK +, EPSN, ITS,AAG +, EPSDIV, KCTLO,INLO,CNR, DNBACK +, USHIROB, IDBG, V,R0,ISTB,AF,AF1 +, DYJ,DJKL,DJKR,DKJL,DKJR, M1F, DMOBN +, GR, EINSN, J1,J2,J3,J4,J5,J6) IMPLICIT REAL*8 (A-H,O-Z) *** *** DIMENSION B1(N) +, F(N+M1), R(N),WK(N+M1+1), V(-M1:N+M1) +, AA(N),AB(-M1:N),AC(-M1:N),DI(-M1:N),P(-M1:N+M1) +, W(-M1:N+M1),X(-M1:N+M1),AP(N),UB(-M1:N),AE(-M1:N) +, AB1(-M1:N+M1),AC1(-M1:N+M1),UB1(-M1:N+M1),AE1(-M1:N+M1) +, RRA(-M1:N+M1),PA(-M1:N+M1),RR(-M1:N+M1) +, R0(N), DXK(M2+1), AF(-M1:N),AF1(-M1:N+M1) +, PSI(-M1:N+M1), DNL(-M1:N+M1), AAG(N) +, DYJ(M1F),DJKL(M1),DJKR(M1),DKJL(M1),DKJR(M1) +, DMOBN(N+M1), GR(N+M1) +, EINSN(N+M1) * C% DO 101 I=-M1, N+M1 V(I)= 0.0D0 AF1(I)= 0.0D0 101 CONTINUE DO 102 I=-M1, N 102 AF(I)= 0.0D0 DO 103 I=1, N 103 R0(I)= 0.0D0 C% C * EPSICW=EPSN/EPSDIV * IF(IDBG.GT.1) THEN WRITE(6,*) ' ' WRITE(6,*) ' <<<<<<<<<<<<<<<<<< NINPUT START >>>>>>>>>>>>>>>>>>' ENDIF C// C// * CALL ZEROCR(N,M1,B1,F,R, AA,AB,AC,DI,P,W,X,AP,UB,AE +, AB1,AC1,UB1,AE1,RRA,PA,RR, AAG) CC CALL FINMAB(N,M1, M2, WK,DXK, F +, PSI,DNDOCM , DNBACK, DMOBN +, DYJ,M1F , GR, EINSN) C * CALL MATRIB(N,M1,M2,AA,AB,AC,WK,DXK,AB1,AC1 +, PSI,DYJ,DJKL,DJKR,DKJL,DKJR, M1F +, DMOBN, EINSN) CC * CMX=-1.0D0 CMY=-1.0D0 DO 222 I= 2,M2 DO 223 J= 2,M1 K= M1*(I-1)+J * CPEKX= DABS(EINSN(K)*(PSI(K)-PSI(K-M1))) IF(CPEKX.GT.CMX) THEN CMX=CPEKX MXJ=J MXI=I ENDIF C// CPEKY= DABS(EINSN(K)*(PSI(K)-PSI(K-1))) IF(CPEKY.GT.CMY) THEN CMY=CPEKY MYJ=J MYI=I ENDIF * 223 CONTINUE 222 CONTINUE C// DO 522 I= 1,M2 DO 523 J= 1,M1 K= M1*(I-1)+J IF(AA(K).LE.0.0D0) THEN WRITE(6,*) ' .... AA<0(NINPUT)',J,I,AA(K) ENDIF * 523 CONTINUE 522 CONTINUE * CC BNRM=0.0D0 DO 110 I=1,N 110 BNRM=BNRM+F(I)*F(I) BNRM=DSQRT(BNRM) * * ISTB=1 IF(ISTB.EQ.0) THEN CALL ILUCM2(N,M1,AA,AB,AC,AE,UB,DI, USHIROB + ,AB1,AC1,AE1,UB1) ELSE *** CALL UPILU3(N,M1,AA,AB,AC,AE,AF,UB,DI, USHIROB<- UTOP CALL UPILU4(N,M1,AA,AB,AC,AE,AF,UB,DI, USHIROB + ,AB1,AC1,AE1,AF1,UB1) ENDIF CC DO 602 I=0, M2-1 K=M1*I+1 602 DNL(K)= 0.0D0 * DO 604 I=0, M2-1 K=M1*I+M1 604 DNL(K)= 0.0D0 CC CC DO 226 I=1,N 226 X(I)=DNL(I) * *** START *** KCT=0 INMX= 4 DO 230 IN=1,INMX * IF(ISTB.EQ.0) THEN CALL BCG12(AA,AB,AC,AE,UB,DI,AB1,AC1,AE1,UB1,P,R,W,F + ,X,AP,EPSICW,ERR,M1,N,KCOUNT,BNRM,RRA,PA,RR) ELSE CALL BCGSTB3(AA,AB,AC,AE,AF,UB,DI,AB1,AC1,AE1,AF1,UB1,P,V,W,F + ,R0 ,X,AP,EPSICW,ERR,M1,N,KCOUNT,BNRM,RRA,PA,RR) ENDIF CC * KCT=KCT+KCOUNT CNR=0.0D0 DO 116 I=1,N B1(I)=F(I)-(AC1(I)*X(I-M1)+AB1(I)*X(I-1) + +AA(I)*X(I)+AB(I)*X(I+1)+AC(I)*X(I+M1)) 116 CNR=CNR+B1(I)*B1(I) CNR=DSQRT(CNR)/BNRM * IF(IDBG.GT.1) THEN WRITE(6,5800) ITS, X(J1),X(J2),X(J3),X(J4),X(J5),X(J6) + , CNR,IN,KCT,BNRM ENDIF * IF(CNR.LT.EPSN) GO TO 250 EPSICW=EPSICW/2.0 * 230 CONTINUE 250 CONTINUE * DO 220 I=1,N 220 DNL(I)=X(I) C KCTLO=KCTLO+KCT IF(IN.EQ.INMX+1) IN=IN-1 INLO=INLO+IN * DO 802 I=0, M2-1 K=M1*I+1 802 DNL(K)= DNDOCM * DO 804 I=0, M2-1 K=M1*I+M1 804 DNL(K)= DNBACK CC C// C IF(IDBG.GT.1) THEN WRITE(6,*) ' <<<<<<<<<<<<<<<<<< END >>>>>>>>>>>>>>>>>>' WRITE(6,22) ' *** CMX,J,I ***',CMX,MXJ,MXI WRITE(6,22) ' *** CMY,J,I ***',CMY,MYJ,MYI 22 FORMAT(1H ,A16,F7.3,2I5) WRITE(6,*) ' ' ENDIF 5800 FORMAT(1H,I4,D12.4,D12.4,D12.4,D12.4,D12.4,D12.4,D20.7,I5,I6 + ,D15.7) RETURN END * ** C *** FUJITSU REAL*8 FUNCTION BER4(BEX) REAL*8 BEX, HX C *** HITACHI C REAL FUNCTION BER4*8(BEX) C REAL*8 BER4, BEX, HX C** HX=1.000D-4 HX=1.000D-2 * IF(BEX.LT.-HX) THEN BER4=BEX/(DEXP(BEX)-1.0D0) ELSE IF(BEX.GE.-HX .AND. BEX.LE.HX) THEN BER4=1.0D0-(0.5D0-0.0833333*BEX)*BEX * CCCC BER4=1.0D0-0.5D0*BEX ELSE IF(BEX.GT.HX) THEN BER4=BEX*DEXP(-BEX)/(1.0D0-DEXP(-BEX)) ENDIF * END * ** SUBROUTINE MATRIB(N,M,M2,A0,A1,A2,WK,DXK,AB1,AC1 +, PSI,DYJ,DJKL,DJKR,DKJL,DKJR, M1F +, DMOBN, EINSN) * IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A0(N),A1(-M:N),A2(-M:N),WK(N+M+1),AB1(-M:N+M) +, AC1(-M:N+M),DXK(M2+1),PSI(-M:N+M),DYJ(M1F),DJKL(M),DJKR(M) +, DKJL(M), DKJR(M) +, DMOBN(N+M), EINSN(N+M) * C// DO 99 I=1,N+M+1 99 WK(I)=0.0D0 ** DO 222 I= 2, M2 DO 223 J= 2,M K= M *(I-1)+J * WK(K) = DMOBN(K)/EINSN(K) * 223 CONTINUE 222 CONTINUE C// * * DO 100 J = 1, M2 DO 40 I=2,M IP= I CC DJKL(I)=(DYJ(IP)/DXK(J))*0.5D0 DJKR(I)=(DYJ(IP)/DXK(J+1))*0.5D0 DKJL(I)=(DXK(J)/DYJ(IP))*0.5D0 40 DKJR(I)=(DXK(J+1)/DYJ(IP))*0.5D0 * I=M*(J-1)+1 A0(I)= 1.0D0 A1(I)= 0.0D0 A2(I)= 0.0D0 AB1(I)=0.0D0 AC1(I)=0.0D0 CC CC * DO 110 K= 2, M-1 I=M*(J-1)+K DJVL= DJKL(K) DJVR=DJKR(K) DKVL=DKJL(K) DKVR=DKJR(K) DJWL= DJKL(K+1) DJWR=DJKR(K+1) DKWL=DKJL(K+1) DKWR=DKJR(K+1) GRADM1=PSI(I-1)-PSI(I) GRADP1=PSI(I+1)-PSI(I) GRADMM=PSI(I-M)-PSI(I) GRADPM=PSI(I+M)-PSI(I) A0(I)= (WK(I)*DJVL) *BER4(-EINSN(I)*GRADM1) + +(WK(I+M)*DJVR)*BER4(-EINSN(I+M)*GRADM1) + +(WK(I+1)*DJWL) *BER4(-EINSN(I+1)*GRADP1) + +(WK(I+M+1)*DJWR)*BER4(-EINSN(I+M+1)*GRADP1) + +(WK(I)*DKVL) *BER4(-EINSN(I)*GRADMM) + +(WK(I+1)*DKWL)*BER4(-EINSN(I+1)*GRADMM) + +(WK(I+M)*DKVR) *BER4(-EINSN(I+M)*GRADPM) + +(WK(I+M+1)*DKWR)*BER4(-EINSN(I+M+1)*GRADPM) AB1(I)= -(WK(I)*DJVL) *BER4(EINSN(I)*GRADM1) + -(WK(I+M)*DJVR)*BER4(EINSN(I+M)*GRADM1) A1(I)=-(WK(I+1)*DJWL) *BER4(EINSN(I+1)*GRADP1) + -(WK(I+M+1)*DJWR)*BER4(EINSN(I+M+1)*GRADP1) AC1(I)=-(WK(I)*DKVL) *BER4(EINSN(I)*GRADMM) + -(WK(I+1)*DKWL)*BER4(EINSN(I+1)*GRADMM) A2(I)=-(WK(I+M)*DKVR) *BER4(EINSN(I+M)*GRADPM) + -(WK(I+M+1)*DKWR)*BER4(EINSN(I+M+1)*GRADPM) 110 CONTINUE * CCC * I=M*(J-1)+ M A0(I)= 1.0D0 A1(I)= 0.0D0 A2(I)= 0.0D0 AB1(I)=0.0D0 AC1(I)=0.0D0 CCC * 100 CONTINUE * DO 801 I= N-M+1, N 801 A2(I)=0.0D0 A1(N)= 0.0D0 AB1(1)= 0.0D0 DO 780 I= 1, M 780 AC1(I)= 0.0D0 CC * M0C LINE NO KIRIKAKI O TSUKURU *** * DO 806 I=0, M2-1 806 AB1(M*I+2)= 0.0D0 CC DO 813 I=0,M2-1 J=M-1 813 A1(M*I+J)= 0.0D0 *C RETURN END * ** SUBROUTINE FINMAB(NL,ML1, M2, WKL,DXK, F +, PSI,DNDOCM , DNBACK, DMOBN +, DYJ,M1F , GR, EINSN) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION WKL(NL+ML1),DXK(M2+1), PSI(-ML1:NL+ML1), F(NL+ML1) +, DMOBN(NL+ML1), DYJ(M1F) +, GR(NL+ML1), EINSN(NL+ML1) * C// DO 99 I=1,NL+ML1 99 WKL(I)=0.0D0 ** DO 222 I= 2, M2 DO 223 J= 2,ML1 K= ML1*(I-1)+J * WKL(K) = DMOBN(K)/EINSN(K) * 223 CONTINUE 222 CONTINUE C// * CC// * G-R NO SET *** CCC DO 301 I=0, M2-1 C/ DO 302 J=2, ML1-1 C KL=ML1*I+J CC 302 F(KL)=GR(KL) C/ 301 CONTINUE * CCC * DO 620 I=0, M2-1 620 F(ML1*I+1)= 0.0D0 * DO 640 I=0, M2-1 640 F(ML1*I+ML1)= 0.0D0 ** CC DMJF= DYJ(2) DMJM= DYJ(ML1) ** ** DO 420 J= 1, M2 * DMKL=DXK(J) DMKR=DXK(J+1) * DJFL=(DMJF/DMKL)*0.5D0 DJFR=(DMJF/DMKR)*0.5D0 DKFL=(DMKL/DMJF)*0.5D0 DKFR=(DMKR/DMJF)*0.5D0 * I= ML1*(J-1)+2 GRADM1=PSI(I-1)-PSI(I) AB1 = -(WKL(I)*DJFL) *BER4(EINSN(I)*GRADM1) + -(WKL(I+ML1)*DJFR)*BER4(EINSN(I+ML1)*GRADM1) 420 F(I)=F(I)+ (-AB1*DNDOCM) * * CC DO 440 J= 1, M2 * DMKL=DXK(J) DMKR=DXK(J+1) * DJ4L=(DMJM/DMKL)*0.5D0 DJ4R=(DMJM/DMKR)*0.5D0 * I= ML1*(J-1)+ML1-1 GRADP1=PSI(I+1)-PSI(I) A1 = -(WKL(I+1)*DJ4L) *BER4(EINSN(I+1)*GRADP1) + -(WKL(I+ML1+1)*DJ4R)*BER4(EINSN(I+ML1+1)*GRADP1) C// F(I)=F(I) -A1*DNBACK 440 CONTINUE CC * *C RETURN END CC ****************************************************************** * ** SUBROUTINE PINPU1(N,M1,M2, DNL,PSI, DNDOCM +, B1,F,R,WK, AA,AB,AC,DI,P,W,X,AP,UB,AE +, AB1,AC1,UB1,AE1,RRA,PA,RR,DXK +, EPSN, ITS,AAG +, EPSDIP, KCTLO,INLO,CNR +, DNBACK, USHIROB, IDBG, V,R0,ISTB,AF,AF1 +, DYJ,DJKL,DJKR,DKJL,DKJR, M1F, DMOBN +, GR, EINSN, J1,J2,J3,J4,J5,J6) IMPLICIT REAL*8 (A-H,O-Z) *** *** DIMENSION B1(N) +, F(N+M1), R(N),WK(N+M1+1), V(-M1:N+M1) +, AA(N),AB(-M1:N),AC(-M1:N),DI(-M1:N),P(-M1:N+M1) +, W(-M1:N+M1),X(-M1:N+M1),AP(N),UB(-M1:N),AE(-M1:N) +, AB1(-M1:N+M1),AC1(-M1:N+M1),UB1(-M1:N+M1),AE1(-M1:N+M1) +, RRA(-M1:N+M1),PA(-M1:N+M1),RR(-M1:N+M1) +, R0(N), DXK(M2+1), AF(-M1:N),AF1(-M1:N+M1) +, PSI(-M1:N+M1), DNL(-M1:N+M1), AAG(N) +, DYJ(M1F),DJKL(M1),DJKR(M1),DKJL(M1),DKJR(M1) +, DMOBN(N+M1), GR(N+M1) +, EINSN(N+M1) * C DO 101 I=-M1, N+M1 V(I)= 0.0D0 AF1(I)= 0.0D0 101 CONTINUE DO 102 I=-M1, N 102 AF(I)= 0.0D0 DO 103 I=1, N 103 R0(I)= 0.0D0 CC CC * EPSICW=EPSN/EPSDIP * IF(IDBG.GT.1) THEN WRITE(6,*) ' ' WRITE(6,*) ' <<<<<<<<<<<<<<<<<< PINPU1 START >>>>>>>>>>>>>>>>>>' ENDIF C// C// * CALL ZEROCR(N,M1,B1,F,R, AA,AB,AC,DI,P,W,X,AP,UB,AE +, AB1,AC1,UB1,AE1,RRA,PA,RR, AAG) CC CALL FINMAP(N,M1, M2, WK,DXK, F +, PSI,DNDOCM , DNBACK, DMOBN +, DYJ,M1F , GR, EINSN) * CALL MATRIP(N,M1,M2,AA,AB,AC,WK,DXK,AB1,AC1 +, PSI,DYJ,DJKL,DJKR,DKJL,DKJR, M1F +, DMOBN, EINSN) CC * CMX=-1.0D0 CMY=-1.0D0 DO 222 I= 2, M2 DO 223 J= 2,M1 K= M1*(I-1)+J * CPEKX= DABS(EINSN(K)*(PSI(K)-PSI(K-M1))) IF(CPEKX.GT.CMX) THEN CMX=CPEKX MXJ=J MXI=I ENDIF C// CPEKY= DABS(EINSN(K)*(PSI(K)-PSI(K-1))) IF(CPEKY.GT.CMY) THEN CMY=CPEKY MYJ=J MYI=I ENDIF * 223 CONTINUE 222 CONTINUE C// DO 522 I= 1,M2 DO 523 J= 1,M1 K= M1*(I-1)+J IF(AA(K).LE.0.0D0) THEN WRITE(6,*) ' .... AA<0(PINPUT)',J,I,AA(K) ENDIF * 523 CONTINUE 522 CONTINUE * CC BNRM=0.0D0 DO 110 I=1,N 110 BNRM=BNRM+F(I)*F(I) BNRM=DSQRT(BNRM) * * ISTB=1 * IF(ISTB.EQ.0) THEN CALL ILUCM2(N,M1,AA,AB,AC,AE,UB,DI, USHIROB + ,AB1,AC1,AE1,UB1) ELSE *** CALL UPILU3(N,M1,AA,AB,AC,AE,AF,UB,DI, USHIROB<- UTOP CALL UPILU4(N,M1,AA,AB,AC,AE,AF,UB,DI, USHIROB + ,AB1,AC1,AE1,AF1,UB1) ENDIF CC DO 602 I=0, M2-1 K=M1*I+1 602 DNL(K)= 0.0D0 * DO 604 I=0, M2-1 K=M1*I+M1 604 DNL(K)= 0.0D0 CC CC DO 226 I=1,N 226 X(I)=DNL(I) * *** START *** KCT=0 INMX= 4 DO 230 IN=1,INMX * IF(ISTB.EQ.0) THEN CALL BCG12(AA,AB,AC,AE,UB,DI,AB1,AC1,AE1,UB1,P,R,W,F + ,X,AP,EPSICW,ERR,M1,N,KCOUNT,BNRM,RRA,PA,RR) ELSE CALL BCGSTB3(AA,AB,AC,AE,AF,UB,DI,AB1,AC1,AE1,AF1,UB1,P,V,W,F + ,R0 ,X,AP,EPSICW,ERR,M1,N,KCOUNT,BNRM,RRA,PA,RR) ENDIF CC * KCT=KCT+KCOUNT CNR=0.0D0 DO 116 I=1,N B1(I)=F(I)-(AC1(I)*X(I-M1)+AB1(I)*X(I-1) + +AA(I)*X(I)+AB(I)*X(I+1)+AC(I)*X(I+M1)) 116 CNR=CNR+B1(I)*B1(I) CNR=DSQRT(CNR)/BNRM * IF(IDBG.GT.1) THEN WRITE(6,5800) ITS, X(J1),X(J2),X(J3),X(J4),X(J5),X(J6) + , CNR,IN,KCT,BNRM ENDIF * IF(CNR.LT.EPSN) GO TO 250 EPSICW=EPSICW/2.0 * 230 CONTINUE 250 CONTINUE * DO 220 I=1,N 220 DNL(I)=X(I) C KCTLO=KCTLO+KCT IF(IN.EQ.INMX+1) IN=IN-1 INLO=INLO+IN * DO 802 I=0, M2-1 K=M1*I+1 802 DNL(K)= DNDOCM * DO 804 I=0, M2-1 K=M1*I+M1 804 DNL(K)= DNBACK CC C// C IF(IDBG.GT.1) THEN WRITE(6,*) ' <<<<<<<<<<<<<<<<<< END >>>>>>>>>>>>>>>>>>' WRITE(6,22) ' *** CMX,J,I ***',CMX,MXJ,MXI WRITE(6,22) ' *** CMY,J,I ***',CMY,MYJ,MYI 22 FORMAT(1H ,A16,F7.3,2I5) WRITE(6,*) ' ' ENDIF 5800 FORMAT(1H,I4,D12.4,D12.4,D12.4,D12.4,D12.4,D12.4,D20.7,I5,I6 + ,D15.7) RETURN END * ** SUBROUTINE MATRIP(N,M,M2,A0,A1,A2,WK,DXK,AB1,AC1 +, PSI,DYJ,DJKL,DJKR,DKJL,DKJR, M1F +, DMOBN, EINSN) * IMPLICIT REAL*8 (A-H,O-Z) DIMENSION A0(N),A1(-M:N),A2(-M:N),WK(N+M+1),AB1(-M:N+M) +, AC1(-M:N+M),DXK(M2+1),PSI(-M:N+M),DYJ(M1F),DJKL(M),DJKR(M) +, DKJL(M), DKJR(M) +, DMOBN(N+M), EINSN(N+M) * C// DO 99 I=1,N+M+1 99 WK(I)=0.0D0 ** DO 222 I= 2, M2 DO 223 J= 2,M K= M *(I-1)+J * WK(K) = DMOBN(K)/EINSN(K) * 223 CONTINUE 222 CONTINUE C// * * DO 100 J = 1, M2 DO 40 I=2,M IP=I CC DJKL(I)=(DYJ(IP)/DXK(J))*0.5D0 DJKR(I)=(DYJ(IP)/DXK(J+1))*0.5D0 DKJL(I)=(DXK(J)/DYJ(IP))*0.5D0 40 DKJR(I)=(DXK(J+1)/DYJ(IP))*0.5D0 * I=M*(J-1)+1 A0(I)= 1.0D0 A1(I)= 0.0D0 A2(I)= 0.0D0 AB1(I)=0.0D0 AC1(I)=0.0D0 C * DO 110 K= 2, M-1 I=M*(J-1)+K DJVL= DJKL(K) DJVR=DJKR(K) DKVL=DKJL(K) DKVR=DKJR(K) DJWL= DJKL(K+1) DJWR=DJKR(K+1) DKWL=DKJL(K+1) DKWR=DKJR(K+1) GRADM1=PSI(I-1)-PSI(I) GRADP1=PSI(I+1)-PSI(I) GRADMM=PSI(I-M)-PSI(I) GRADPM=PSI(I+M)-PSI(I) A0(I)= (WK(I)*DJVL) *BER4( EINSN(I)*GRADM1) + +(WK(I+M)*DJVR)*BER4( EINSN(I+M)*GRADM1) + +(WK(I+1)*DJWL) *BER4( EINSN(I+1)*GRADP1) + +(WK(I+M+1)*DJWR)*BER4( EINSN(I+M+1)*GRADP1) + +(WK(I)*DKVL) *BER4( EINSN(I)*GRADMM) + +(WK(I+1)*DKWL)*BER4( EINSN(I+1)*GRADMM) + +(WK(I+M)*DKVR) *BER4( EINSN(I+M)*GRADPM) + +(WK(I+M+1)*DKWR)*BER4( EINSN(I+M+1)*GRADPM) AB1(I)= -(WK(I)*DJVL) *BER4(-EINSN(I)*GRADM1) + -(WK(I+M)*DJVR)*BER4(-EINSN(I+M)*GRADM1) A1(I)=-(WK(I+1)*DJWL) *BER4(-EINSN(I+1)*GRADP1) + -(WK(I+M+1)*DJWR)*BER4(-EINSN(I+M+1)*GRADP1) AC1(I)=-(WK(I)*DKVL) *BER4(-EINSN(I)*GRADMM) + -(WK(I+1)*DKWL)*BER4(-EINSN(I+1)*GRADMM) A2(I)=-(WK(I+M)*DKVR) *BER4(-EINSN(I+M)*GRADPM) + -(WK(I+M+1)*DKWR)*BER4(-EINSN(I+M+1)*GRADPM) 110 CONTINUE * CCC * I=M*(J-1)+ M A0(I)= 1.0D0 A1(I)= 0.0D0 A2(I)= 0.0D0 AB1(I)=0.0D0 AC1(I)=0.0D0 CCC * 100 CONTINUE * DO 801 I= N-M+1, N 801 A2(I)=0.0D0 A1(N)= 0.0D0 AB1(1)= 0.0D0 DO 780 I= 1, M 780 AC1(I)= 0.0D0 CC * M0C LINE NO KIRIKAKI O TSUKURU *** * DO 806 I=0, M2-1 806 AB1(M*I+2)= 0.0D0 CC DO 813 I=0,M2-1 J=M-1 813 A1(M*I+J)= 0.0D0 * RETURN END * ** SUBROUTINE FINMAP(NL,ML1, M2, WKL,DXK, F +, PSI,DNDOCM , DNBACK, DMOBN +, DYJ,M1F , GR, EINSN) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION WKL(NL+ML1),DXK(M2+1), PSI(-ML1:NL+ML1), F(NL+ML1) +, DMOBN(NL+ML1), DYJ(M1F) +, GR(NL+ML1), EINSN(NL+ML1) * C// DO 99 I=1,NL+ML1 99 WKL(I)=0.0D0 ** DO 222 I= 2, M2 DO 223 J= 2,ML1 K= ML1*(I-1)+J * WKL(K) = DMOBN(K)/EINSN(K) * 223 CONTINUE 222 CONTINUE C// * ** CC// * G-R NO SET *** CCC DO 301 I=0, M2-1 C/ DO 302 J=2, ML1-1 C KL=ML1*I+J CC 302 F(KL)=GR(KL) C/ 301 CONTINUE * CCC * DO 620 I=0, M2-1 620 F(ML1*I+1)= 0.0D0 * DO 640 I=0, M2-1 640 F(ML1*I+ML1)= 0.0D0 ** CC ** DMJF= DYJ(2) DMJM= DYJ(ML1) ** DO 420 J= 1, M2 * DMKL=DXK(J) DMKR=DXK(J+1) * DJFL=(DMJF/DMKL)*0.5D0 DJFR=(DMJF/DMKR)*0.5D0 DKFL=(DMKL/DMJF)*0.5D0 DKFR=(DMKR/DMJF)*0.5D0 * I= ML1*(J-1)+2 GRADM1=PSI(I-1)-PSI(I) AB1 = -(WKL(I)*DJFL) *BER4(-EINSN(I)*GRADM1) + -(WKL(I+ML1)*DJFR)*BER4(-EINSN(I+ML1)*GRADM1) 420 F(I)=F(I)+ (-AB1*DNDOCM) * CC DO 440 J= 1, M2 * DMKL=DXK(J) DMKR=DXK(J+1) * DJ4L=(DMJM/DMKL)*0.5D0 DJ4R=(DMJM/DMKR)*0.5D0 * I= ML1*(J-1)+ML1-1 GRADP1=PSI(I+1)-PSI(I) A1 = -(WKL(I+1)*DJ4L) *BER4(-EINSN(I+1)*GRADP1) + -(WKL(I+ML1+1)*DJ4R)*BER4(-EINSN(I+ML1+1)*GRADP1) C// F(I)=F(I)+ (-A1*DNBACK) 440 CONTINUE CC * RETURN END CC CC C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SUBROUTINE NPLUST( N,M1,M2,M0Q, DIELDM, WK + ,YMK,MJM,MJ1, MSD,INDST, FDIEL,YWID, JNCUT) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION WK(N+M1), YMK(M1) C DO 21 I=1,M1 21 YMK(I)=0.0D0 C BMJM= 0.1D0/ DFLOAT(MJM) BMJ1= 0.1D0/ DFLOAT(MJ1) C DO 610 I=M0Q, MSD, -1 610 YMK(I)= -(I-M0Q)*BMJM C DO 620 I=MSD, MSD-MJ1, -1 620 YMK(I)= 0.1D0 -(I-MSD)*BMJ1 C J=M0Q CC do while( YMK(J).LT.YWID ) 90 CONTINUE IF(YMK(J).LT.YWID) THEN JNCUT=J J=J-1 GO TO 90 ENDIF CC enddo C IF(INDST.EQ.1) THEN AA=DIELDM*(FDIEL-1.0D0)/YWID BB=DIELDM C DO 220 J=2,M2 DO 220 K=JNCUT, M0Q I=M1*(J-1)+K YR=YWID- YMK(K) 220 WK(I)=AA*YR +BB ENDIF CC CC IF(INDST.EQ.2) THEN SIG6=YWID*DSQRT(-0.5D0/DLOG(FDIEL)) C DO 320 J=2,M2 DO 320 K=JNCUT, M0Q I=M1*(J-1)+K YR=YWID- YMK(K) 320 WK(I)=DIELDM*DEXP(-0.5D0*(YR**2/SIG6**2)) ENDIF C WRITE(6,*) ' ## NPLUST ##' WRITE(6,*) ' INDST,FDIEL,DIELDM,YWID,JNCUT=' +, INDST,FDIEL,DIELDM,YWID,JNCUT WRITE(6,*) ' AA,BB,SIG6=' +, AA,BB,SIG6 WRITE(6,*) ' ## NPLUST END ##' C RETURN END CC CC C SUBROUTINE PPLUST( N,M1,M2,M0Q, DIELDM, WK + ,YMK,MJM,MJ2, MSU,INDST, FDIEL,YWID, JPCUT) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION WK(N+M1), YMK(M1) C C DO 21 I=1,M1 C 21 YMK(I)=0.0D0 C BMJM= 0.1D0/ DFLOAT(MJM) BMJ2= 0.1D0/ DFLOAT(MJ2) C DO 610 I=M0Q, MSU 610 YMK(I)= (I-M0Q)*BMJM C DO 620 I=MSU, MSU+MJ2 620 YMK(I)= 0.1D0 +(I-MSU)*BMJ2 C J=M0Q CC do while( YMK(J).LT.YWID ) 90 CONTINUE IF(YMK(J).LT.YWID) THEN JPCUT=J J=J+1 GO TO 90 ENDIF CC enddo C IF(INDST.EQ.1) THEN AA=DIELDM*(FDIEL-1.0D0)/YWID BB=DIELDM C DO 220 J=2,M2 DO 220 K=JPCUT, M0Q+1, -1 I=M1*(J-1)+K YR=YWID- YMK(K) 220 WK(I)=AA*YR +BB ENDIF CC CC IF(INDST.EQ.2) THEN SIG6=YWID*DSQRT(-0.5D0/DLOG(FDIEL)) C DO 320 J=2,M2 DO 320 K=JPCUT, M0Q+1, -1 I=M1*(J-1)+K YR=YWID- YMK(K) 320 WK(I)=DIELDM*DEXP(-0.5D0*(YR**2/SIG6**2)) ENDIF C WRITE(6,*) ' ## PPLUST ##' WRITE(6,*) ' INDST,FDIEL,DIELDM,YWID,JPCUT=' +, INDST,FDIEL,DIELDM,YWID,JPCUT WRITE(6,*) ' AA,BB,SIG6=' +, AA,BB,SIG6 WRITE(6,*) ' ## PPLUST END ##' RETURN END CC C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CC CC SUBROUTINE NDDOPE( N,M1,M2,M0Q, QND,DIELDM, WK,QNDMN + ,YMK,MJM,MJ1, MSD,INDST, FDIEL,YWID, JNCUT) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION QND(N), WK(N+M1), QNDMN(N+M1) + , YMK(M1) C DO 110 I=1,N+M1 QNDMN(I)= 0.0D0 110 WK(I)= 0.0D0 C C DO 120 J=2,M2 DO 120 K=2,M0Q I= M1*(J-1)+K 120 WK(I)=DIELDM C CALL NPLUST( N,M1,M2,M0Q, DIELDM, WK + ,YMK,MJM,MJ1, MSD,INDST, FDIEL,YWID, JNCUT) C DO 410 J=2,M2 DO 410 K=2,M0Q I= M1*(J-1)+K 410 QNDMN(I)=WK(I) C C DO 500 J= 1, M2-2 DO 500 K= 2, M0Q-1 I= M1*J+K 500 QND(I)= (WK(I)+WK(I+M1)+WK(I+M1+1)+WK(I+1))*0.25D0 C DO 510 J= 1, M2-2 K= 1 I= M1*J+K 510 QND(I)= (WK(I+M1+1)+WK(I+1))*0.5D0 C DO 520 J= 1, M2-2 K= M0Q I= M1*J+K 520 QND(I)= (WK(I+M1)+WK(I))*0.5D0 C J= 0 DO 530 K= 2, M0Q-1 I= M1*J+K 530 QND(I)= (WK(I+M1)+WK(I+M1+1))*0.5D0 C J= M2-1 DO 540 K= 2, M0Q-1 I= M1*J+K 540 QND(I)= (WK(I)+WK(I+1))*0.5D0 C QND(1)= WK(2+M1) QND(M0Q)= WK(M0Q+M1) QND(M1*(M2-1)+1)= WK(M1*(M2-1)+2) QND(M1*(M2-1)+M0Q)= WK(M1*(M2-1)+M0Q) C RETURN END CCC CCC SUBROUTINE NADOPE( N,M1,M2,M0Q, QNA, DIELDM, WK,QNAMN + ,YMK,MJM,MJ2, MSU,INDST, FDIEL,YWID, JPCUT) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION QNA(N), WK(N+M1), QNAMN(N+M1) + , YMK(M1) C DO 110 I=1,N+M1 QNAMN(I)= 0.0D0 110 WK(I)= 0.0D0 C C DO 120 J=2,M2 DO 120 K=M0Q+1,M1 I= M1*(J-1)+K 120 WK(I)=DIELDM C CALL PPLUST( N,M1,M2,M0Q, DIELDM, WK + ,YMK,MJM,MJ2, MSU,INDST, FDIEL,YWID, JPCUT) C DO 410 J=2,M2 DO 410 K=M0Q+1,M1 I= M1*(J-1)+K 410 QNAMN(I)=WK(I) C DO 500 J= 1, M2-2 DO 500 K= M0Q+1, M1-1 I= M1*J+K 500 QNA(I)= (WK(I)+WK(I+M1)+WK(I+M1+1)+WK(I+1))*0.25D0 C DO 510 J= 1, M2-2 K= M0Q I= M1*J+K 510 QNA(I)= (WK(I+M1+1)+WK(I+1))*0.5D0 C DO 520 J= 1, M2-2 K= M1 I= M1*J+K 520 QNA(I)= (WK(I+M1)+WK(I))*0.5D0 C J= 0 DO 530 K= M0Q+1, M1-1 I= M1*J+K 530 QNA(I)= (WK(I+M1)+WK(I+M1+1))*0.5D0 C J= M2-1 DO 540 K= M0Q+1, M1-1 I= M1*J+K 540 QNA(I)= (WK(I)+WK(I+1))*0.5D0 C QNA(M0Q)= WK(M0Q+M1+1) QNA(M1)= WK(M1+M1) QNA(M1*(M2-1)+M0Q)= WK(M1*(M2-1)+M0Q+1) QNA(M1*(M2-1)+M1)= WK(M1*(M2-1)+M1) C RETURN END *PROCESS NOOPLIST,NOSOURCE CCC CCC SUBROUTINE ZEROCR(N,M1,B1,F,R, AA,AB,AC,DI,P,W,X,AP,UB,AE +, AB1,AC1,UB1,AE1,RRA,PA,RR, AAG) IMPLICIT REAL*8 (A-H,O-Z) *** *** DIMENSION B1(N), F(N+M1), R(N) +, AA(N),AB(-M1:N),AC(-M1:N),DI(-M1:N),P(-M1:N+M1) +, W(-M1:N+M1),X(-M1:N+M1),AP(N),UB(-M1:N),AE(-M1:N) +, AB1(-M1:N+M1),AC1(-M1:N+M1),UB1(-M1:N+M1),AE1(-M1:N+M1) +, RRA(-M1:N+M1),PA(-M1:N+M1),RR(-M1:N+M1), AAG(N) * DO 10 I=1, N+M1 10 F(I)= 0.0D0 * DO 20 I=1, N AAG(I)= 0.0D0 AP(I)= 0.0D0 B1(I)= 0.0D0 R(I)= 0.0D0 20 AA(I)= 0.0D0 * DO 30 I= -M1, N AE(I)= 0.0D0 UB(I)= 0.0D0 DI(I)= 0.0D0 AC(I)= 0.0D0 30 AB(I)= 0.0D0 * DO 40 I= -M1, N+M1 AB1(I)= 0.0D0 AC1(I)= 0.0D0 UB1(I)= 0.0D0 AE1(I)= 0.0D0 RR(I)= 0.0D0 PA(I)= 0.0D0 RRA(I)= 0.0D0 X(I)= 0.0D0 W(I)= 0.0D0 40 P(I)= 0.0D0 C RETURN END CCCC CCCC CC *PROCESS NOOPLIST,NOSOURCE * * SUBROUTINE ICDCMP3( N,M ,AA,AB,AC,AE,AF,UB,DI,USHI ) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M:N),AC(-M:N),DI(-M:N), + AE(-M:N), AF(-M:N),UB(-M:N) DO 100 I = 1, N DIV=AA(I)-UB(I-1)*UB(I-1)*DI(I-1)-AC(I-M)*AC(I-M)*DI(I-M) + -AE(I-M+1)*AE(I-M+1)*DI(I-M+1) + -AF(I-M+2)*AF(I-M+2)*DI(I-M+2) + -USHI*(AF(I-M+2)*AC(I-M+2)*DI(I-M+2) + +AF(I-M)*AC(I-M)*DI(I-M) + +UB(I-1)*AF(I-1)*DI(I-1)+UB(I-M+2)*AF(I-M+2)*DI(I-M+2)) DI(I)=1.0D0/DIV UB(I)=AB(I)-AC(I-M+1)*AE(I-M+1)*DI(I-M+1) + -AE(I-M+2)*AF(I-M+2)*DI(I-M+2) AE(I)=-UB(I-1)*AC(I-1)*DI(I-1) 100 AF(I)=-UB(I-1)*AE(I-1)*DI(I-1) RETURN END * * SUBROUTINE ICCG13(AA,AB,AC,AE,AF,UB,DI,P,R,W,B0 + ,X,AP,EPSLON,ERR,M1,N,KCOUNT,BNRM) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M1:N),AC(-M1:N),DI(-M1:N) +, P(-M1:N+M1),W(-M1:N+M1),B0(N),X(-M1:N+M1) +, AP(N),R(N),UB(-M1:N),AE(-M1:N),AF(-M1:N) DO 10 I = 1, N 10 R(I)=B0(I)-(AA(I)*X(I)+AB(I)*X(I+1)+AC(I)*X(I+M1) + +AB(I-1)*X(I-1)+AC(I-M1)*X(I-M1)) *** DO 30 I = 1, N 30 W(I)=(R(I)-AC(I-M1)*W(I-M1)-UB(I-1)*W(I-1) + -AE(I-M1+1)*W(I-M1+1)-AF(I-M1+2)*W(I-M1+2))*DI(I) DO 40 I = N, 1, -1 40 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1) + +AE(I)*W(I+M1-1)+AF(I)*W(I+M1-2) ) RR1 = 0.0D0 DO 60 I = 1, N P(I) = W(I) 60 RR1 = RR1+R(I)*W(I) *** ITERATION *** DO 1000 K = 1, 200 RAP = 0.0D0 DO 100 I = 1, N AP(I)=AA(I)*P(I)+AB(I)*P(I+1)+AC(I)*P(I+M1) + +AB(I-1)*P(I-1)+AC(I-M1)*P(I-M1) 100 RAP = RAP+P(I)*AP(I) ALPHA = RR1/RAP ERR = 0.0D0 DO 200 I = 1, N X(I) = X(I)+ALPHA*P(I) R(I) = R(I)-ALPHA*AP(I) 200 ERR = ERR+R(I)*R(I) ERR = DSQRT(ERR)/BNRM IF( ERR.LT.EPSLON ) GO TO 2000 *** DO 230 I = 1, N 230 W(I)=(R(I)-AC(I-M1)*W(I-M1)-UB(I-1)*W(I-1) + -AE(I-M1+1)*W(I-M1+1)-AF(I-M1+2)*W(I-M1+2))*DI(I) DO 240 I = N, 1, -1 240 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1) + +AE(I)*W(I+M1-1)+AF(I)*W(I+M1-2) ) RR2 = 0.0D0 DO 300 I = 1, N 300 RR2 = RR2+R(I)*W(I) BETA = RR2/RR1 RR1 = RR2 DO 400 I = 1, N 400 P(I) = W(I)+BETA*P(I) 1000 CONTINUE 2000 KCOUNT = K RETURN END * * SUBROUTINE UPILU3(N,M,AA,AB,AC,AE,AF,UB,DI,U,AB1,AC1,AE1,AF1,UB1) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M:N),AC(-M:N),DI(-M:N), + AE(-M:N),AF(-M:N),UB(-M:N), + AB1(-M:N+M), AC1(-M:N+M),AE1(-M:N+M),AF1(-M:N+M),UB1(-M:N+M) *** CC// ** U:USHIRO PARAMETER NOT USE ** U=0.0D0 CC DO 100 I = 1, N *** CC W= ABS(AC1(I))+ABS(AB1(I))+ABS(AB(I))+ABS(AC(I)) CC UUI = MAX(W,AA(I)) *** AE1(I)= -AC1(I)*UB(I-M)*DI(I-M) AF1(I)= -AE1(I)*UB(I-M+1)*DI(I-M+1) UB1(I)= AB1(I)-AC1(I)*AE(I-M)*DI(I-M)-AE1(I)*AF(I-M+1)*DI(I-M+1) * DIV= AA(I)-UB1(I)*UB(I-1)*DI(I-1)-AC1(I)*AC(I-M)*DI(I-M) + -AE1(I)*AE(I-M+1)*DI(I-M+1)-AF1(I)*AF(I-M+2)*DI(I-M+2) DI(I)=1.0D0/DIV * UB(I)= AB(I)-AE1(I)*AC(I-M+1)*DI(I-M+1) + -AF1(I)*AE(I-M+2)*DI(I-M+2) AF(I)= -UB1(I)*AE(I-1)*DI(I-1) AE(I)= -UB1(I)*AC(I-1)*DI(I-1) * 100 CONTINUE RETURN END * * * SUBROUTINE BCGSTB3(AA,AB,AC,AE,AF,UB,DI,AB1,AC1,AE1,AF1,UB1,P,V,W + ,B0 ,R0 ,X,AP,EPSLON,ERR,M1,N,KCOUNT,BNRM,AEE,EE,RR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M1:N),AC(-M1:N),DI(-M1:N) +, P(-M1:N+M1),W(-M1:N+M1),B0(N),X(-M1:N+M1) +, AP(N),V(-M1:N+M1),UB(-M1:N),AE(-M1:N),AF(-M1:N) +, AB1(-M1:N+M1),AC1(-M1:N+M1),UB1(-M1:N+M1),AE1(-M1:N+M1) +, AF1(-M1:N+M1),AEE(-M1:N+M1),EE(-M1:N+M1),RR(-M1:N+M1),R0(N) * DO 10 I = 1, N 10 V(I)=B0(I)-(AA(I)*X(I)+AB(I)*X(I+1)+AC(I)*X(I+M1) + +AB1(I)*X(I-1)+AC1(I)*X(I-M1)) *** KOKO DE WA V = B -AX DE ARU *** * DO 30 I = 1, N 30 W(I)=(V(I)-AC1(I)*W(I-M1)-UB1(I)*W(I-1)-AE1(I)*W(I-M1+1) + -AF1(I)*W(I-M1+2))*DI(I) * DO 40 I = N, 1, -1 40 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1)+AE(I)*W(I+M1-1) + +AF(I)*W(I+M1-2)) *** KOKO DE WA W = INV(LU)*(B-AX) DE ARU ****** CR= 0.0D0 DO 60 I = 1, N P(I) = W(I) RR(I)= W(I) R0(I)= W(I) 60 CR = CR + W(I)*W(I) *** ITERATION; KONO ATO, V O V,W O Q TO SHITE TSUKAU *** DO 1000 K = 1, 400 DO 100 I = 1, N AP(I)=AA(I)*P(I)+AB(I)*P(I+1)+AC(I)*P(I+M1) + +AB1(I)*P(I-1)+AC1(I)*P(I-M1) 100 CONTINUE *** DO 110 I = 1, N 110 W(I)=(AP(I)-AC1(I)*W(I-M1)-UB1(I)*W(I-1)-AE1(I)*W(I-M1+1) + -AF1(I)*W(I-M1+2))*DI(I) * DO 120 I = N, 1, -1 120 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1)+AE(I)*W(I+M1-1) + +AF(I)*W(I+M1-2)) *** W = INV(LU)*AP DE ARU *** SIG= 0.0D0 DO 130 I=1, N 130 SIG = SIG+R0(I)*W(I) ALFA = CR/SIG * DO 140 I=1,N 140 EE(I)= RR(I)-ALFA*W(I) * DO 150 I = 1, N 150 AEE(I)= AA(I)*EE(I)+AB(I)*EE(I+1)+AC(I)*EE(I+M1) + +AB1(I)*EE(I-1)+AC1(I)*EE(I-M1) *** DO 160 I = 1, N 160 V(I)=(AEE(I)-AC1(I)*V(I-M1)-UB1(I)*V(I-1)-AE1(I)*V(I-M1+1) + -AF1(I)*V(I-M1+2))*DI(I) * DO 170 I = N, 1, -1 170 V(I)=V(I)-DI(I)*(UB(I)*V(I+1)+AC(I)*V(I+M1)+AE(I)*V(I+M1-1) + +AF(I)*V(I+M1-2)) *** V = INV(LU)*A*E DE ARU *** EEV= 0.0D0 VV= 0.0D0 DO 180 I= 1,N EEV= EEV+EE(I)*V(I) 180 VV = VV +V(I)*V(I) OMEG =EEV/VV * DO 190 I= 1, N X(I)= X(I)+ALFA*P(I)+OMEG*EE(I) 190 RR(I)=EE(I)-OMEG*V(I) * CR=0.0D0 ERR = 0.0D0 DO 200 I = 1, N CR = CR+ R0(I)*RR(I) 200 ERR = ERR+RR(I)*RR(I) ERR = DSQRT(ERR)/BNRM IF( ERR.LT.EPSLON ) GO TO 2000 * BETA = CR/(OMEG*SIG) DO 300 I = 1, N 300 P(I)= RR(I)+BETA*(P(I)-OMEG*W(I)) * 1000 CONTINUE 2000 KCOUNT = K RETURN END * * SUBROUTINE BCGSTB2(AA,AB,AC,AE,UB,DI,AB1,AC1,AE1,UB1,P,V,W,B0 + ,R0 ,X,AP,EPSLON,ERR,M1,N,KCOUNT,BNRM,AEE,EE,RR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M1:N),AC(-M1:N),DI(-M1:N) +, P(-M1:N+M1),W(-M1:N+M1),B0(N),X(-M1:N+M1) +, AP(N),V(-M1:N+M1),UB(-M1:N),AE(-M1:N) +, AB1(-M1:N+M1),AC1(-M1:N+M1),UB1(-M1:N+M1),AE1(-M1:N+M1) +, AEE(-M1:N+M1),EE(-M1:N+M1),RR(-M1:N+M1),R0(N) * DO 10 I = 1, N 10 V(I)=B0(I)-(AA(I)*X(I)+AB(I)*X(I+1)+AC(I)*X(I+M1) + +AB1(I)*X(I-1)+AC1(I)*X(I-M1)) *** KOKO DE WA V = B -AX DE ARU *** DO 30 I = 1, N 30 W(I)=(V(I)-AC1(I)*W(I-M1)-UB1(I)*W(I-1)-AE1(I)*W(I-M1+1))*DI(I) DO 40 I = N, 1, -1 40 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1)+AE(I)*W(I+M1-1)) *** KOKO DE WA W = INV(LU)*(B-AX) DE ARU ****** CR= 0.0D0 DO 60 I = 1, N P(I) = W(I) RR(I)= W(I) R0(I)= W(I) 60 CR = CR + W(I)*W(I) *** ITERATION; KONO ATO, V O V,W O Q TO SHITE TSUKAU *** DO 1000 K = 1, 400 DO 100 I = 1, N AP(I)=AA(I)*P(I)+AB(I)*P(I+1)+AC(I)*P(I+M1) + +AB1(I)*P(I-1)+AC1(I)*P(I-M1) 100 CONTINUE DO 110 I = 1, N 110 W(I)=(AP(I)-AC1(I)*W(I-M1)-UB1(I)*W(I-1)-AE1(I)*W(I-M1+1))*DI(I) DO 120 I = N, 1, -1 120 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1)+AE(I)*W(I+M1-1)) *** W = INV(LU)*AP DE ARU *** SIG= 0.0D0 DO 130 I=1, N 130 SIG = SIG+R0(I)*W(I) ALFA = CR/SIG * DO 140 I=1,N 140 EE(I)= RR(I)-ALFA*W(I) * * DO 150 I = 1, N 150 AEE(I)= AA(I)*EE(I)+AB(I)*EE(I+1)+AC(I)*EE(I+M1) + +AB1(I)*EE(I-1)+AC1(I)*EE(I-M1) DO 160 I = 1, N 160 V(I)=(AEE(I)-AC1(I)*V(I-M1)-UB1(I)*V(I-1)-AE1(I)*V(I-M1+1))*DI(I) DO 170 I = N, 1, -1 170 V(I)=V(I)-DI(I)*(UB(I)*V(I+1)+AC(I)*V(I+M1)+AE(I)*V(I+M1-1)) *** V = INV(LU)*A*E DE ARU *** EEV= 0.0D0 VV= 0.0D0 DO 180 I= 1,N EEV= EEV+EE(I)*V(I) 180 VV = VV +V(I)*V(I) OMEG =EEV/VV * DO 190 I= 1, N X(I)= X(I)+ALFA*P(I)+OMEG*EE(I) 190 RR(I)=EE(I)-OMEG*V(I) * CR=0.0D0 ERR = 0.0D0 DO 200 I = 1, N CR = CR+ R0(I)*RR(I) 200 ERR = ERR+RR(I)*RR(I) ERR = DSQRT(ERR)/BNRM IF( ERR.LT.EPSLON ) GO TO 2000 **************************************** KOKOMADE ******************* BETA = CR/(OMEG*SIG) DO 300 I = 1, N 300 P(I)= RR(I)+BETA*(P(I)-OMEG*W(I)) * 1000 CONTINUE 2000 KCOUNT = K RETURN END * * SUBROUTINE ICDCMP( N,M ,AA,AB,AC,AE,UB,DI,U) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M:N),AC(-M:N),DI(-M:N), + AE(-M:N), UB(-M:N) DO 100 I = 1, N DIV=AA(I)-UB(I-1)*UB(I-1)*DI(I-1)-AC(I-M)*AC(I-M)*DI(I-M) + -AE(I-M+1)*AE(I-M+1)*DI(I-M+1) + -(UB(I-1)*AE(I-1)*DI(I-1)+UB(I-M+1)*AE(I-M+1)*DI(I-M+1))*U DI(I)=1.0D0/DIV UB(I)=AB(I)-AC(I-M+1)*AE(I-M+1)*DI(I-M+1) 100 AE(I)=-UB(I-1)*AC(I-1)*DI(I-1) RETURN END * * SUBROUTINE ICCG12(AA,AB,AC,AE,UB,DI,P,R,W,B0 + ,X,AP,EPSLON,ERR,M1,N,KCOUNT,BNRM) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M1:N),AC(-M1:N),DI(-M1:N) +, P(-M1:N+M1),W(-M1:N+M1),B0(N),X(-M1:N+M1) +, AP(N),R(N),UB(-M1:N),AE(-M1:N) DO 10 I = 1, N 10 R(I)=B0(I)-(AA(I)*X(I)+AB(I)*X(I+1)+AC(I)*X(I+M1) + +AB(I-1)*X(I-1)+AC(I-M1)*X(I-M1)) *** DO 30 I = 1, N 30 W(I)=(R(I)-AC(I-M1)*W(I-M1)-UB(I-1)*W(I-1) + -AE(I-M1+1)*W(I-M1+1))*DI(I) DO 40 I = N, 1, -1 40 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1) + +AE(I)*W(I+M1-1)) RR1 = 0.0D0 DO 60 I = 1, N P(I) = W(I) 60 RR1 = RR1+R(I)*W(I) *** ITERATION *** DO 1000 K = 1, 200 RAP = 0.0D0 DO 100 I = 1, N AP(I)=AA(I)*P(I)+AB(I)*P(I+1)+AC(I)*P(I+M1) + +AB(I-1)*P(I-1)+AC(I-M1)*P(I-M1) 100 RAP = RAP+P(I)*AP(I) ALPHA = RR1/RAP ERR = 0.0D0 DO 200 I = 1, N X(I) = X(I)+ALPHA*P(I) R(I) = R(I)-ALPHA*AP(I) 200 ERR = ERR+R(I)*R(I) ERR = DSQRT(ERR)/BNRM IF( ERR.LT.EPSLON ) GO TO 2000 *** DO 230 I = 1, N 230 W(I)=(R(I)-AC(I-M1)*W(I-M1)-UB(I-1)*W(I-1) + -AE(I-M1+1)*W(I-M1+1))*DI(I) DO 240 I = N, 1, -1 240 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1) + +AE(I)*W(I+M1-1)) RR2 = 0.0D0 DO 300 I = 1, N 300 RR2 = RR2+R(I)*W(I) BETA = RR2/RR1 RR1 = RR2 DO 400 I = 1, N 400 P(I) = W(I)+BETA*P(I) 1000 CONTINUE 2000 KCOUNT = K RETURN END * * SUBROUTINE ILUCM2( N,M ,AA,AB,AC,AE,UB,DI,U,AB1,AC1,AE1,UB1) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M:N),AC(-M:N),DI(-M:N), + AE(-M:N),UB(-M:N), + AB1(-M:N+M), AC1(-M:N+M),AE1(-M:N+M),UB1(-M:N+M) DO 100 I = 1, N * AE1(I)= -AC1(I)*UB(I-M)*DI(I-M) UB1(I)= AB1(I)-AC1(I)*AE(I-M)*DI(I-M) * DIV=AA(I)-UB1(I)*UB(I-1)*DI(I-1)-AC1(I)*AC(I-M)*DI(I-M) + -AE1(I)*AE(I-M+1)*DI(I-M+1) + -(UB1(I)*AE(I-1)*DI(I-1)+AE1(I)*UB(I-M+1)*DI(I-M+1))*U DI(I)=1.0D0/DIV * UB(I)= AB(I)-AE1(I)*AC(I-M+1)*DI(I-M+1) AE(I)= -UB1(I)*AC(I-1)*DI(I-1) * 100 CONTINUE RETURN END * SUBROUTINE ILUCM1( N,M ,AA,AB,AC, DI,U, AB1,AC1 ) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M:N),AC(-M:N),DI(-M:N), + AB1(-M:N+M), AC1(-M:N+M) DO 100 I = 1, N * * DIV=AA(I)-AB1(I)*AB(I-1)*DI(I-1)-AC1(I)*AC(I-M)*DI(I-M) + -(AB1(I)*AC(I-1)*DI(I-1)+AC1(I)*AB(I-M)*DI(I-M))*U DI(I)=1.0D0/DIV * 100 CONTINUE RETURN END * * SUBROUTINE BCG12(AA,AB,AC,AE,UB,DI,AB1,AC1,AE1,UB1,P,R,W,B0 + ,X,AP,EPSLON,ERR,M1,N,KCOUNT,BNRM,RRA,PA,RR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M1:N),AC(-M1:N),DI(-M1:N) +, P(-M1:N+M1),W(-M1:N+M1),B0(N),X(-M1:N+M1) +, AP(N),R(N),UB(-M1:N),AE(-M1:N) +, AB1(-M1:N+M1),AC1(-M1:N+M1),UB1(-M1:N+M1),AE1(-M1:N+M1) +, RRA(-M1:N+M1),PA(-M1:N+M1),RR(-M1:N+M1) * DO 10 I = 1, N 10 R(I)=B0(I)-(AA(I)*X(I)+AB(I)*X(I+1)+AC(I)*X(I+M1) + +AB1(I)*X(I-1)+AC1(I)*X(I-M1)) *** DO 30 I = 1, N 30 W(I)=(R(I)-AC1(I)*W(I-M1)-UB1(I)*W(I-1)-AE1(I)*W(I-M1+1))*DI(I) DO 40 I = N, 1, -1 40 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1)+AE(I)*W(I+M1-1)) RR1 = 0.0D0 DO 60 I = 1, N P(I) = W(I) PA(I)= W(I) RR(I)= W(I) RRA(I)= W(I) 60 RR1 = RR1+RR(I)*RRA(I) *** ITERATION *** DO 1000 K = 1, 200 DO 100 I = 1, N AP(I)=AA(I)*P(I)+AB(I)*P(I+1)+AC(I)*P(I+M1) + +AB1(I)*P(I-1)+AC1(I)*P(I-M1) 100 CONTINUE DO 110 I = 1, N 110 W(I)=(AP(I)-AC1(I)*W(I-M1)-UB1(I)*W(I-1)-AE1(I)*W(I-M1+1))*DI(I) DO 120 I = N, 1, -1 120 W(I)=W(I)-DI(I)*(UB(I)*W(I+1)+AC(I)*W(I+M1)+AE(I)*W(I+M1-1)) RAP= 0.0D0 DO 130 I=1, N 130 RAP = RAP+W(I)*PA(I) ALPHA = RR1/RAP ERR = 0.0D0 DO 200 I = 1, N X(I) = X(I)+ALPHA*P(I) RR(I) = RR(I)-ALPHA*W(I) 200 ERR = ERR+RR(I)*RR(I) ERR = DSQRT(ERR)/BNRM IF( ERR.LT.EPSLON ) GO TO 2000 *** DO 500 I = 1, N 500 W(I)=(PA(I)-AC(I-M1)*W(I-M1)-UB(I-1)*W(I-1) + -AE(I-M1+1)*W(I-M1+1))*DI(I) DO 510 I = N, 1, -1 510 W(I)=W(I)-DI(I)*(UB1(I+1)*W(I+1)+AC1(I+M1)*W(I+M1) + +AE1(I+M1-1)*W(I+M1-1)) DO 520 I= 1, N QA =AA(I)*W(I)+AB1(I+1)*W(I+1)+AC1(I+M1)*W(I+M1) + +AB(I-1)*W(I-1)+AC(I-M1)*W(I-M1) RRA(I)= RRA(I)-ALPHA*QA 520 CONTINUE RR2=RR1 RR1 = 0.0D0 DO 300 I = 1, N 300 RR1 = RR1+RR(I)*RRA(I) BETA = RR1/RR2 DO 400 I = 1, N P(I)= RR(I)+BETA*P(I) 400 PA(I) = RRA(I)+BETA*PA(I) 1000 CONTINUE 2000 KCOUNT = K RETURN END * * SUBROUTINE BCG11(AA,AB,AC, DI,AB1,AC1, P,R,W,B0 + ,X,AP,EPSLON,ERR,M1,N,KCOUNT,BNRM,RRA,PA,RR) IMPLICIT REAL*8 (A-H,O-Z) DIMENSION AA(N),AB(-M1:N),AC(-M1:N),DI(-M1:N) +, P(-M1:N+M1),W(-M1:N+M1),B0(N),X(-M1:N+M1) +, AP(N),R(N) +, AB1(-M1:N+M1),AC1(-M1:N+M1) +, RRA(-M1:N+M1),PA(-M1:N+M1),RR(-M1:N+M1) * DO 10 I = 1, N 10 R(I)=B0(I)-(AA(I)*X(I)+AB(I)*X(I+1)+AC(I)*X(I+M1) + +AB1(I)*X(I-1)+AC1(I)*X(I-M1)) *** DO 30 I = 1, N 30 W(I)= ( R(I)-AC1(I)*W(I-M1)-AB1(I)*W(I-1) )*DI(I) * DO 40 I = N, 1, -1 40 W(I)=W(I)-DI(I)*(AB(I)*W(I+1)+AC(I)*W(I+M1) ) RR1 = 0.0D0 DO 60 I = 1, N P(I) = W(I) PA(I)= W(I) RR(I)= W(I) RRA(I)= W(I) 60 RR1 = RR1+RR(I)*RRA(I) *** ITERATION *** DO 1000 K = 1, 800 * DO 100 I = 1, N AP(I)=AA(I)*P(I)+AB(I)*P(I+1)+AC(I)*P(I+M1) + +AB1(I)*P(I-1)+AC1(I)*P(I-M1) 100 CONTINUE * DO 110 I = 1, N 110 W(I)=(AP(I)-AC1(I)*W(I-M1)-AB1(I)*W(I-1) )*DI(I) DO 120 I = N, 1, -1 120 W(I)=W(I)-DI(I)*(AB(I)*W(I+1)+AC(I)*W(I+M1) ) * RAP= 0.0D0 DO 130 I=1, N 130 RAP = RAP+W(I)*PA(I) ALPHA = RR1/RAP ERR = 0.0D0 DO 200 I = 1, N X(I) = X(I)+ALPHA*P(I) RR(I) = RR(I)-ALPHA*W(I) 200 ERR = ERR+RR(I)*RR(I) ERR = DSQRT(ERR)/BNRM IF( ERR.LT.EPSLON ) GO TO 2000 *** DO 500 I = 1, N 500 W(I)=(PA(I)-AC(I-M1)*W(I-M1)-AB(I-1)*W(I-1))*DI(I) DO 510 I = N, 1, -1 510 W(I)=W(I)-DI(I)*(AB1(I+1)*W(I+1)+AC1(I+M1)*W(I+M1) ) * DO 520 I= 1, N QA =AA(I)*W(I)+AB1(I+1)*W(I+1)+AC1(I+M1)*W(I+M1) + +AB(I-1)*W(I-1)+AC(I-M1)*W(I-M1) RRA(I)= RRA(I)-ALPHA*QA 520 CONTINUE * RR2=RR1 * RR1 = 0.0D0 DO 300 I = 1, N 300 RR1 = RR1+RR(I)*RRA(I) * BETA = RR1/RR2 DO 400 I = 1, N P(I)= RR(I)+BETA*P(I) 400 PA(I) = RRA(I)+BETA*PA(I) 1000 CONTINUE 2000 KCOUNT = K RETURN END CCCC SUBROUTINE CLOCK0 DIMENSION IA(3) COMMON /CLO/ IB(3) C CALL ITIME(IA) DO 200 I=1,3 200 IB(I)=IA(I) CC WRITE(6,*) (IB(I),I=1,3) RETURN END CCCC SUBROUTINE CLOCK(CPUT) DIMENSION IA(3) COMMON /CLO/ IB(3) C CALL ITIME(IA) IJI= IA(1)-IB(1) IFUN=IA(2)-IB(2) IBYO=IA(3)-IB(3) CPUT=IJI*3600.0 +IFUN*60.0 +IBYO*1.0 CC WRITE(6,*) (IB(I),I=1,3),(IA(I),I=1,3),CPUT C CALL ITIME(IA) DO 200 I=1,3 200 IB(I)=IA(I) RETURN END CCC CCC CCC C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CC SUBROUTINE GRWKWR(M1,N, JY0,MJM,JY1,JY2,JY3 +, KXC,KX1,KXA,KXB, KXBL,KXCR +, EINSN,EINSP, GR, NL,ML1, M2) IMPLICIT REAL*8(A-H,O-Z) DIMENSION WK(N+M1), GR(NL+ML1) +, EINSN(NL+ML1), EINSP(NL+ML1) * DO 510 I=1,N+M1 510 WK(I)=0.0D0 * DO 511 I= 0, M2-1 DO 511 J= 1, ML1 K= ML1*I+J 511 WK(K)= GR(K) * WRITE(6,*) ' *** GRWK ***' CALL FWRT(M1,N,WK(1),JY0,MJM,JY1,JY2,JY3 +, KXC,KX1,KXA,KXB, KXBL,KXCR) WRITE(6,*) ' *** GRWK END ***' CC DO 370 I=1,N+M1 370 WK(I)=0.0D0 * DO 371 I= 1, M2-1 DO 371 J= 2, ML1 K= ML1*I+J C// IF(EINSN(K).NE.0.0D0) THEN WK(K)= 11610.145D0/EINSN(K) C// ENDIF 371 CONTINUE * WRITE(6,*) ' *** TN ***' CALL FWRT(M1,N,WK(1),JY0,MJM,JY1,JY2,JY3 +, KXC,KX1,KXA,KXB, KXBL,KXCR) WRITE(6,*) ' *** TN END ***' CC CC DO 380 I=1,N+M1 380 WK(I)=0.0D0 * DO 381 I= 1, M2-1 DO 381 J= 2, ML1 K= ML1*I+J C// IF(EINSP(K).NE.0.0D0) THEN WK(K)= 11610.145D0/EINSP(K) C// ENDIF 381 CONTINUE * WRITE(6,*) ' *** TP ***' CALL FWRT(M1,N,WK(1),JY0,MJM,JY1,JY2,JY3 +, KXC,KX1,KXA,KXB, KXBL,KXCR) WRITE(6,*) ' *** TP END ***' C// RETURN END C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% C CC SUBROUTINE BPNWR(M1,N, JY0,MJM,JY1,JY2,JY3 +, KXC,KX1,KXA,KXB, KXBL,KXCR +, BW,DN,DP) IMPLICIT REAL*8(A-H,O-Z) DIMENSION BW(-M1:N+M1),DN(N),DP(N) * WRITE(6,*) ' *** BW ***' CALL FWRT(M1,N,BW(1),JY0,MJM,JY1,JY2,JY3 +, KXC,KX1,KXA,KXB, KXBL,KXCR) WRITE(6,*) ' *** BW END ***' ** WRITE(6,*) ' ' CC WRITE(6,*) ' *** DN ***' CALL FWRT(M1,N,DN(1),JY0,MJM,JY1,JY2,JY3 +, KXC,KX1,KXA,KXB, KXBL,KXCR) WRITE(6,*) ' *** DN END ***' ** WRITE(6,*) ' ' CCCC WRITE(6,*) ' *** DP ***' CALL FWRT(M1,N,DP(1),JY0,MJM,JY1,JY2,JY3 +, KXC,KX1,KXA,KXB, KXBL,KXCR) WRITE(6,*) ' *** DP END ***' 9100 FORMAT(1H ,9D14.6) C RETURN END CC SUBROUTINE FWRT(M1,N,WW,JY0,MJM,JY1,JY2,JY3 +, KXC,KX1,KXA,KXB, KXBL,KXCR) IMPLICIT REAL*8(A-H,O-Z) DIMENSION WW(N), M1LST(30) * M1LST(1)=1 M1LST(2)=JY0 M1LST(3)=JY1 M1LST(4)=JY1+MJM M1LST(5)=JY2 M1LST(6)=JY3 M1LST(7)=M1 NDIV=7 ** DO 196 K= NDIV, 1, -1 I=M1LST(K) 196 WRITE(6,9100) WW(M1*(KXA-1)+I),WW(M1*(KXB-1)+I) +, WW(M1*(KXBL-1)+I),WW(M1*(KXCR-1)+I) +, WW(M1*(KXC-1)+I),WW(M1*(KX1-1)+I) 9100 FORMAT(1H ,9D14.6) RETURN END