- 积分
- 111
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2011-7-26
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
麻烦高手们帮我看看,为什么老出现Segmentation fault错误
1 PROGRAM chazhi
2 IMPLICIT NONE
3 INTEGER,PARAMETER::IMT=360,JMT=180,KMT=29
4 real(kind=4),dimension(IMT,JMT,KMT)::ts,ss
5 !integer,dimension(JMT)::XX,O
6 real,dimension(JMT)::XX,O
7 real(kind=4),dimension(JMT)::S,DS,DDS,P
8 !integer,allocatable::X(:)
9 real,allocatable::X(:)
10 real,allocatable::Y(:),DY(:),DDY(:),H(:)
11 integer::i,j,k,count,M,N,u,v,w
12 real::DY1,DYN,T
13
14
15 open(30,file="s00mn1_30",form="formatted",status="old")
16 open(60,file="chazhi_30",form="formatted",status="replace")
17 100 format(10f8.4)
18 do k=1,29
19 read(30,100) ((ts(i,j,k),i=1,IMT),j=1,JMT)
20 end do
21
22 do j=1,JMT
23 XX(j)=j
24 end do
25
26 outter: do w=1,KMT
27 inner: do u=1,IMT
28 count=0
29 do v=1,JMT
30 if(ts(u,v,1)/=-99.9999)then
31 count=count+1
32 O(count)=v
33 P(count)=ts(u,v,1)
34 end if
35 end do
36 N=count
37 allocate(X(N))
38 allocate(Y(N))
39 allocate(DY(N))
40 allocate(DDY(N))
41 allocate(H(N))
42
43 X=O(1:N)
44 Y=P(1:N)
45
46 DY1=0
47 DYN=0
48 N=count
49 M=JMT
50 ! print*,X
51 ! print*,Y
52 ! print*,XX
53 ! print*,N,M
54
55 call ESPL1(X,Y,N,DY1,DYN,XX,M,DY,DDY,S,DS,DDS,T,H)
56 ! write(*,100),S
57 print*,"AAA"
58 do j=1,M
59 if(S(j)<0)then
60 S(j)=0.0
61 else if(s(j)>41.4029)then
62 S(j)=41.500
63 end if
64 end do
65 ss(u,1:JMT,1)=S(1:M)
66 ! print*,S
67
68 deallocate(X)
69 deallocate(Y)
70 deallocate(DY)
71 deallocate(DDY)
72 deallocate(H)
73 print*,"ddd"
74 end do inner
75 end do outter
76 print*,"222"
77 close(30)
78 close(60)
79 stop
80 end
100 SUBROUTINE ESPL1(X,Y,N,DY1,DYN,XX,M,DY,DDY,S,DS,DDS,T,H)
101
102
103 IMPLICIT NONE
104 real,DIMENSION(N)::X,Y,DY,DDY,H,T
105 real,DIMENSION(M)::XX,S,DS,DDS
106 INTEGER::DY1,DYN,I,N,M,J
107 real::H0,H1,BETA,ALPHA
108
109 DY(1)=0.0
110 H(1)=DY1
111 H0=X(2)-X(1)
112 DO 10 J=2,N-1
113 H1=X(J+1)-X(J)
114 ALPHA=H0/(H0+H1)
115 BETA=(1.0-ALPHA)*(Y(J)-Y(J-1))/H0
116 BETA=3.0*(BETA+ALPHA*(Y(J+1)-Y(J))/H1)
117 DY(J)=-ALPHA/(2.0+(1.0-ALPHA)*DY(J-1))
118 H(J)=(BETA-(1.0-ALPHA)*H(J-1))
119 H(J)=H(J)/(2.0+(1.0-ALPHA)*DY(J-1))
120 H0=H1
121 10 CONTINUE
122 DY(N)=DYN
123 DO 20 J=N-1,1,-1
124 20 DY(J)=DY(J)*DY(J+1)+H(J)
125 DO 30 J=1,N-1
126 30 H(J)=X(J+1)-X(J)
127 DO 40 J=1,N-1
128 H1=H(J)*H(J)
129 DDY(J)=6.0*(Y(J+1)-Y(J))/H1-2.0*(2.0*DY(J)+DY(J+1))/H(J)
130 40 CONTINUE
131 H1=H(N-1)*H(N-1)
132 DDY(N)=6.0*(Y(N-1)-Y(N))/H1+2.0*(2.0*DY(N)+DY(N-1))/H(N-1)
133 T=0.0
134 DO 50 I=1,N-1
135 H1=0.5*H(I)*(Y(I)+Y(I+1))
136 H1=H1-H(I)*H(I)*H(I)*(DDY(I)+DDY(I+1))/24.0
137 T=T+H1
138 50 CONTINUE
139 DO 70 J=1,M
140 IF (XX(J).GE.X(N)) THEN
141 I=N-1
142 ELSE
143 I=1
144 60 IF (XX(J).GT.X(I+1)) THEN
145 I=I+1
146 GOTO 60
147 END IF
148 END IF
149 H1=(X(I+1)-XX(J))/H(I)
150 S(J)=(3.0*H1*H1-2.0*H1*H1*H1)*Y(I)
151 S(J)=S(J)+H(I)*(H1*H1-H1*H1*H1)*DY(I)
152 DS(J)=6.0*(H1*H1-H1)*Y(I)/H(I)
153 DS(J)=DS(J)+(3.0*H1*H1-2.0*H1)*DY(I)
154 DDS(J)=(6.0-12.0*H1)*Y(I)/(H(I)*H(I))
155 DDS(J)=DDS(J)+(2.0-6.0*H1)*DY(I)/H(I)
156 H1=(XX(J)-X(I))/H(I)
157 S(J)=S(J)+(3.0*H1*H1-2.0*H1*H1*H1)*Y(I+1)
158 S(J)=S(J)-H(I)*(H1*H1-H1*H1*H1)*DY(I+1)
159 DS(J)=DS(J)-6.0*(H1*H1-H1)*Y(I+1)/H(I)
160 DS(J)=DS(J)+(3.0*H1*H1-2.0*H1)*DY(I+1)
161 DDS(J)=DDS(J)+(6.0-12.0*H1)*Y(I+1)/(H(I)*H(I))
162 DDS(J)=DDS(J)-(2.0-6.0*H1)*DY(I+1)/H(I)
163 70 CONTINUE
164 RETURN
165 END
|
|