[Fortran] 纯文本查看 复制代码
CHARACTER*16 FILENAME
CHARACTER*1 S,S1,S2
CHARACTER*16 CDATE,CTIME
COMMON /COM1/ENDA,A1,A2,A3,Z,C1,C2,HM0
DATA MK,PAI,Z/6,3.14159265,0.68/,S1,S2/1HY,1Hy/,
& N,PH,E1,EDA0,RX,US,X0,XE,C1,C2/33,0.5E9,2.21E11,0.05,
& 0.02,1.0,-2.5,1.5,0.31,0.31/,
& CDATE,CTIME/'The date is','The time is'/
OPEN(4,FILE='OUT.DAT',STATUS='UNKNOWN')
OPEN(8,FILE='FILM.DAT',STATUS='UNKNOWN')
OPEN(10,FILE='PRESS.DAT',STATUS='UNKNOWN')
1 FORMAT(20X,A12,I2.2,':',I2.2,':',I4.4)
2 FORMAT(20X,A12,I2.2,':',I2.2,':',I2.2,'.',I2.2)
WRITE(*,*)'Show the example or not (Y or N)?'
READ(*,'(A)')S
IF(S.EQ.S1.OR.S.EQ.S2)GOTO 5
WRITE(*,*)'N,PH,E1,EDA0,RX,US,X0,XE='
READ(*,*)N,PH,E1,EDA0,RX,US,X0,XE
WRITE(*,*)' Change iteration factors or not (Y or N) ?'
READ(*,'(A)')S
IF(S.EQ.S1.OR.S.EQ.S2)THEN
WRITE(*,*)'C1,C2='
READ(*,*)C1,C2
ENDIF
5 CALL GETDAT(IYR,IMON,IDAY)
WRITE(4,1)CDATE,IMON,IDAY,IYR
WRITE(*,1)CDATE,IMON,IDAY,IYR
WRITE(*,20)N,X0,XE,PH,E1,EDA0,RX,US,C1,C2
WRITE(4,20)N,X0,XE,PH,E1,EDA0,RX,US,C1,C2
20 FORMAT(/////,
& 15X,'+--------+--------+--------+--------+--------+',/,
& 15X,'| N | X0 | XE | W(N/M) | E(Pa) |',/,
& 15X,'+--------+--------+--------+--------+--------+',/,
& 15X,'| ',I4,' | ',F6.3,' | ',F6.3,' |',E8.3,'|',E8.3,'|',/,
& 15X,'+--------+--------+--------+--------+--------+',/,
& 15X,'|Eda(PaS)| R (M) | Us(M/S)| C1 | C2 |',/,
& 15X,'+--------+--------+--------+--------+--------+',/,
& 15X,'| ',F6.3,' | ',F6.3,' | ',F6.3,' | ',F6.3,' | ',F6.3,'|',/,
& 15X,'+--------+--------+--------+--------+--------+',/)
H00=0.0
MM=N-1
LMIN=ALOG(N-1.)/ALOG(2.)-1.99
LMAX=LMIN
U=EDA0*US/(2.*E1*RX)
A1=ALOG(EDA0)+9.67
A2=5.1E-9*PH
A3=0.59/(PH*1.E-9)
B=PAI*PH*RX/E1
W=2.*PAI*PH/(3.*E1)*(B/RX)**2
ALFA=Z*5.1E-9*A1
G=ALFA*E1
AM=W*(2.*U)**(-0.75)
AL=G*(2.*U)**(0.25)
HM0=3.63*(RX/B)**2*G**0.49*U**0.68*W**(-0.073)
ENDA=12.*U*(E1/PH)*(RX/B)**3
WRITE(*,*)' Wait please '
CALL SUBAK(MM)
CALL MULTI(N,LMIN,LMAX,MK,X0,XE,H00)
CALL GETDAT(IYR,IMON,IDAY)
WRITE(4,1)CDATE,IMON,IDAY,IYR
WRITE(*,1)CDATE,IMON,IDAY,IYR
STOP
END
C**************************************************************
SUBROUTINE MULTI(N,LMIN,LMAX,MK,X0,XE,H00)
DIMENSION X(100),Y(100),H(5000),RO(5000),R0(5000),
& EPS(5000),P(5000),F(5000),R(5000),G(20)
COMMON /COMK/KL
DATA NMAX,G0/100,2.0943951/
NN=(N+1)/2
KL=LMIN
CALL KNDX(KL,N,NS,NMAX,DX,X0,XE,X,Y)
CALL INITI(N,X,Y,P(NS))
12 CALL HREE(N,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),0)
14 DO 100 L=LMIN,LMAX
KL=L
G(KL)=G0
20 KK=2
CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),R(NS),0)
KK=1
CALL CHAN(N,R0,R(NS))
CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),R0,1)
CALL GCAL(N,G1,P(NS))
G(KL-1)=G(KL)-DX*DX*G1
N2=N
NS0=NS
KL=KL-1
CALL KNDX(KL,N,NS,NMAX,DX,X0,XE,X,Y)
CALL TRANS(N,N2,N2*N2,H,RO,EPS,P(NS),P(NS0),R(NS),R0)
CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),R(NS),2)
CALL GCAL(N,G1,P(NS))
G(KL)=G(KL)+DX*DX*G1
G0=G(KL)
CALL HREE(N,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),1)
IF(KL.NE.1)GOTO 20
KK=15
CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),R(NS),0)
40 NS0=NS
KL=KL+1
N1=N
CALL KNDX(KL,N,NS,NMAX,DX,X0,XE,X,Y)
G0=G(KL)
CALL TRAP(N1,N,P(NS0),P(NS))
CALL HREE(N,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),0)
KK=1
CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),R(NS),0)
IF(KL.LT.L)GOTO 40
CALL CHAN(N,R0,R(NS))
CALL ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,P(NS),F(NS),R0,1)
100 CONTINUE
M=M+1
IF(M.LT.MK)GOTO 14
120 CALL OUPT(N,X,Y,H,P(NS))
RETURN
END
C**************************************************************
SUBROUTINE OUPT(N,X,Y,H,P)
DIMENSION X(N),Y(N),H(N,N),P(N,N)
NN=(N+1)/2
A=0.0
WRITE(8,*)N,NN
WRITE(8,110)A,(Y(I),I=1,N)
DO I=1,N
WRITE(8,110)X(I),(H(I,J),J=1,N)
END DO
WRITE(10,*)N,NN
WRITE(10,110)A,(Y(I),I=1,N)
DO I=1,N
WRITE(10,110)X(I),(P(I,J),J=1,N)
END DO
110 FORMAT(34(1X,E12.6))
RETURN
END
C**************************************************************
SUBROUTINE TRAP(N1,N2,P1,P2)
DIMENSION P1(N1,N1),P2(N2,N2)
DO 10 I=1,N1-1
II=2*I
DO 10 J=1,N1-1
JJ=2*J
10 P2(II,JJ)=P2(II,JJ)+0.25*(P1(I,J)+P1(I+1,J)+
& P1(I+1,J+1)+P1(I,J+1)-P2(II-1,JJ-1)-P2(II-1,JJ+1)-
& P2(II+1,JJ+1)-P2(II+1,JJ-1))
DO 20 I=1,N1-1
II=2*I
DO 20 J=1,N1
JJ=2*J-1
P2(II,JJ)=P2(II,JJ)+0.5*(P1(I,J)+P1(I+1,J)-
& P2(II-1,JJ)-P2(II+1,JJ))
20 P2(JJ,II)=P2(JJ,II)+0.5*(P1(J,I)+P1(J,I+1)-
& P2(JJ,II-1)-P2(JJ,II+1))
DO 30 I=1,N1
II=2*I-1
DO 30 J=1,N1
JJ=2*J-1
30 P2(II,JJ)=P1(I,J)
RETURN
END
C**************************************************************
SUBROUTINE INITI(N,X,Y,P)
DIMENSION X(N),Y(N),P(N,N)
NN=(N+1)/2
DO 10 I=1,N
D=1.-X(I)*X(I)
DO 10 J=1,NN
C=D-Y(J)*Y(J)
IF(C.LE.0.0)P(I,J)=0.0
10 IF(C.GT.0.0)P(I,J)=SQRT(C)
DO 20 I=1,N
DO 20 J=NN+1,N
JJ=N-J+1
20 P(I,J)=P(I,JJ)
RETURN
END
C**************************************************************
SUBROUTINE CHAN(N,R1,R2)
DIMENSION R1(N,N),R2(N,N)
DO 10 I=1,N
DO 10 J=1,N
10 R1(I,J)=R2(I,J)
RETURN
END
C**************************************************************
SUBROUTINE GCAL(N,G,P)
DIMENSION P(N,N)
G=0.
DO 10 I=1,N
DO 10 J=1,N
10 G=G+P(I,J)
RETURN
END
C**************************************************************
SUBROUTINE KNDX(K,N,NS,NMAX,DX,X0,XE,X,Y)
DIMENSION INDEX(6),NNDEX(6),X(NMAX),Y(NMAX)
DATA INDEX,NNDEX/1,82,371,1460,5685,16616,
& 9,17,33,65,129,257/
N=NNDEX(K)
NS=INDEX(K)
DX=(XE-X0)/(N-1.)
Y0=-0.5*(XE-X0)
DO 10 I=1,N
X(I)=X0+(I-1)*DX
10 Y(I)=Y0+(I-1)*DX
RETURN
END
C**************************************************************
SUBROUTINE HREE(N,DX,H00,G0,X,Y,H,RO,EPS,P,F,KG)
DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),EPS(N,N),F(N,N)
DIMENSION W(150,150),P0(150,150)
COMMON /COM1/ENDA,A1,A2,A3,Z,C1,C2,HM0/COMK/KL
& /COMAK/AK(0:100,0:100)
DATA KK,NW,PAI1/0,150,0.2026423/
NN=(N+1)/2
CALL VI(N,DX,P,W,NW)
HMIN=1.E3
DO 30 I=1,N
W1=0.5*X(I)*X(I)
DO 30 J=1,NN
H0=W1+0.5*Y(J)*Y(J)+W(I,J)
IF(KG.EQ.1)GOTO 20
IF(H0+F(I,J).LT.HMIN)HMIN=H0+F(I,J)
H(I,J)=H0
GOTO 30
20 F(I,J)=H(I,J)-H00-H0
F(I,N-J+1)=F(I,J)
30 CONTINUE
IF(KK.EQ.0)H00=HM0-HMIN
KK=1
IF(KG.EQ.1)RETURN
H0=H00+HMIN
IF(H0.LE.0.0)GOTO 40
IF(KL.NE.1)GOTO 50
W1=0.0
DO 32 I=1,N
DO 32 J=1,N
32 W1=W1+P(I,J)
W1=DX*DX*W1/G0
DW=1.-W1
HM0=H0*W1
40 H00=-HMIN+HM0
50 DO 60 I=1,N
DO 60 J=1,NN
IF(P(I,J).LT.0.0)P(I,J)=0.0
EDA=EXP(A1*(-1.+(1.+A2*P(I,J))**Z))
H(I,J)=H00+H(I,J)+F(I,J)
RO(I,J)=(A3+1.34*P(I,J))/(A3+P(I,J))
60 EPS(I,J)=RO(I,J)*H(I,J)**3/(ENDA*EDA)
DO 70 J=NN+1,N
JJ=N-J+1
DO 70 I=1,N
H(I,J)=H(I,JJ)
RO(I,J)=RO(I,JJ)
70 EPS(I,J)=EPS(I,JJ)
RETURN
END
[Fortran] 纯文本查看 复制代码
C**************************************************************
SUBROUTINE ITER(N,KK,DX,H00,G0,X,Y,H,RO,EPS,P,F,R,KG)
DIMENSION X(N),Y(N),P(N,N),H(N,N),RO(N,N),
& EPS(N,N),F(N,N),R(N,N)
DIMENSION D(200),A(1000),B(600),ID(200)
COMMON /COM1/ENDA,A1,A2,A3,Z,C1,C2,C3
& /COMAK/AK(0:100,0:100)
DATA KG1,PAI1/0,0.2026423/
IF(KG1.NE.0)GOTO 2
KG1=1
AK00=AK(0,0)
AK10=AK(1,0)
AK20=AK(2,0)
BK00=AK00-AK10
BK10=AK10-0.25*(AK00+2.*AK(1,1)+AK(2,0))
BK20=AK20-0.25*(AK10+2.*AK(2,1)+AK(3,0))
2 NN=(N+1)/2
MM=N-1
DX1=1./DX
DX2=DX*DX
DX3=1./DX2
DX4=0.3*DX2
DO 100 K=1,KK
DO 70 J=2,NN
J0=J-1
J1=J+1
JJ=N-J+1
IF(KG.NE.0)GOTO 20
IA=1
8 MM=N-IA
IF(P(MM,J0).GT.1.E-6)GOTO 20
IF(P(MM,J).GT.1.E-6)GOTO 20
IF(P(MM,J1).GT.1.E-6)GOTO 20
IA=IA+1
IF(IA.LT.N)GOTO 8
GOTO 70
20 IF(MM.LT.N-1)MM=MM+1
D2=0.5*(EPS(1,J)+EPS(2,J))
DO 50 I=2,MM
I0=I-1
I1=I+1
II=5*I0
D1=D2
D2=0.5*(EPS(I1,J)+EPS(I,J))
D4=0.5*(EPS(I,J0)+EPS(I,J))
D5=0.5*(EPS(I,J1)+EPS(I,J))
P1=P(I0,JJ)
P2=P(I1,JJ)
P3=P(I,JJ)
P4=P(I,JJ+1)
P5=P(I,JJ-1)
D3=D1+D2+D4+D5
IF(KG.NE.0)GOTO 32
IF(J.EQ.NN.AND.ID(I).EQ.1)P(I,J)=P(I,J)-0.5*C2*D(I)
IF(D1.GE.DX4)GOTO 30
IF(D2.GE.DX4)GOTO 30
IF(D4.GE.DX4)GOTO 30
IF(D5.GE.DX4)GOTO 30
ID(I)=1
IF(J.EQ.NN)P5=P4
A(II+1)=PAI1*(RO(I0,J)*BK10-RO(I,J)*BK20)
A(II+2)=DX3*(D1+0.25*D3)+PAI1*(RO(I0,J)*BK00-RO(I,J)*BK10)
A(II+3)=-1.25*DX3*D3+PAI1*(RO(I0,J)*BK10-RO(I,J)*BK00)
A(II+4)=DX3*(D2+0.25*D3)+PAI1*(RO(I0,J)*BK20-RO(I,J)*BK10)
A(II+5)=-DX3*(D1*P1+D2*P2+D4*P4+D5*P5-D3*P3)+
& DX1*(RO(I,J)*H(I,J)-RO(I0,J)*H(I0,J))+R(I,J)
GOTO 50
30 ID(I)=0
P4=P(I,J0)
IF(J.EQ.NN)P5=P4
A(II+1)=PAI1*(RO(I0,J)*AK10-RO(I,J)*AK20)
A(II+2)=DX3*D1+PAI1*(RO(I0,J)*AK00-RO(I,J)*AK10)
A(II+3)=-DX3*D3+PAI1*(RO(I0,J)*AK10-RO(I,J)*AK00)
A(II+4)=DX3*D2+PAI1*(RO(I0,J)*AK20-RO(I,J)*AK10)
A(II+5)=-DX3*(D1*P1+D2*P2+D4*P4+D5*P5-D3*P3)+
& DX1*(RO(I,J)*H(I,J)-RO(I0,J)*H(I0,J))+R(I,J)
GOTO 50
32 IF(KG.EQ.2)GOTO 40
R(I,J)=-DX3*(D1*P1+D2*P2+D4*P4+D5*P5-D3*P3)+
& DX1*(RO(I,J)*H(I,J)-RO(I0,J)*H(I0,J))+R(I,J)
GOTO 50
40 R(I,J)=DX3*(D1*P1+D2*P2+D4*P4+D5*P5-D3*P3)-
& DX1*(RO(I,J)*H(I,J)-RO(I0,J)*H(I0,J))+R(I,J)
50 CONTINUE
IF(KG.NE.0)GOTO 70
CALL TRA4(MM,D,A,B)
DO 54 I=2,MM
IF(ID(I).EQ.0)GOTO 52
DD=D(I+1)
IF(I.EQ.MM)DD=0
P(I,J)=P(I,J)+C2*(D(I)-0.25*(D(I-1)+DD))
IF(J0.NE.1)P(I,J0)=P(I,J0)-0.25*C2*D(I)
IF(P(I,J0).LT.0.)P(I,J0)=0.0
IF(J1.GE.NN)GOTO 54
P(I,J1)=P(I,J1)-0.25*C2*D(I)
GOTO 54
52 P(I,J)=P(I,J)+C1*D(I)
54 IF(P(I,J).LT.0.0)P(I,J)=0.0
70 CONTINUE
IF(KG.NE.0)GOTO 100
DO 80 J=1,NN
JJ=N+1-J
DO 80 I=1,N
80 P(I,JJ)=P(I,J)
CALL HREE(N,DX,H00,G0,X,Y,H,RO,EPS,P,F,0)
100 CONTINUE
RETURN
END
C**************************************************************
SUBROUTINE TRANS(N1,N2,NMAX,H,RO,EPS,P1,P2,R1,R2)
DIMENSION H(NMAX),RO(NMAX),EPS(NMAX),
& P1(N1,N1),P2(N2,N2),R1(N1,N1),R2(N2,N2)
DO 10 I=1,N1
II=2*I-1
NI1=(I-1)*N1
NI2=2*(I-1)*N2
DO 10 J=1,N1
JJ=2*J-1
NJ1=NI1+J
NJ2=NI2+JJ
H(NJ1)=H(NJ2)
RO(NJ1)=RO(NJ2)
EPS(NJ1)=EPS(NJ2)
P1(I,J)=P2(II,JJ)
10 R1(I,J)=R2(II,JJ)
RETURN
END
C**************************************************************
SUBROUTINE TRA4(N,D,A,B)
DIMENSION D(N),A(5,N),B(3,N)
C=1./A(3,N)
B(1,N)=-A(1,N)*C
B(2,N)=-A(2,N)*C
B(3,N)=A(5,N)*C
DO 10 I=1,N-2
IN=N-I
IN1=IN+1
C=1./(A(3,IN)+A(4,IN)*B(2,IN1))
B(1,IN)=-A(1,IN)*C
B(2,IN)=-(A(2,IN)+A(4,IN)*B(1,IN1))*C
10 B(3,IN)=(A(5,IN)-A(4,IN)*B(3,IN1))*C
D(1)=0.0
D(2)=B(3,2)
DO 20 I=3,N
20 D(I)=B(1,I)*D(I-2)+B(2,I)*D(I-1)+B(3,I)
RETURN
END
C**************************************************************
SUBROUTINE VI(N,DX,P,V,NW)
DIMENSION P(N,N),V(NW,NW)
COMMON /COMAK/AK(0:100,0:100)
PAI1=0.2026423
DO 40 I=1,N
DO 40 J=1,N
H0=0.0
DO 30 K=1,N
IK=IABS(I-K)
DO 30 L=1,N
JL=IABS(J-L)
30 H0=H0+AK(IK,JL)*P(K,L)
40 V(I,J)=H0*DX*PAI1
RETURN
END
C**************************************************************
SUBROUTINE SUBAK(MM)
COMMON /COMAK/AK(0:100,0:100)
S(X,Y)=X+SQRT(X**2+Y**2)
DO 10 I=0,MM
XP=I+0.5
XM=I-0.5
DO 10 J=0,I
YP=J+0.5
YM=J-0.5
A1=S(YP,XP)/S(YM,XP)
A2=S(XM,YM)/S(XP,YM)
A3=S(YM,XM)/S(YP,XM)
A4=S(XP,YP)/S(XM,YP)
AK(I,J)=XP*ALOG(A1)+YM*ALOG(A2)+
& XM*ALOG(A3)+YP*ALOG(A4)
10 AK(J,I)=AK(I,J)
RETURN
END
C*********************************************************