|
Example 1 - shallow.f - Basic serial code
Download shallow.f code.
1. PROGRAM SHALOW
2. C
3. C BENCHMARK WEATHER PREDICTION PROGRAM FOR COMPARING THE
4. C PREFORMANCE OF CURRENT SUPERCOMPUTERS. THE MODEL IS
5. C BASED OF THE PAPER - THE DYNAMICS OF FINITE-DIFFERENCE
6. C MODELS OF THE SHALLOW-WATER EQUATIONS, BY ROBERT SADOURNY
7. C J. ATM. SCIENCES, VOL 32, NO 4, APRIL 1975.
8. C
9. C CODE BY PAUL N. SWARZTRAUBER, NATIONAL CENTER FOR
10. C ATMOSPHERIC RESEARCH, BOULDER, CO, OCTOBER 1984.
11. C
12. DIMENSION U(65,65),V(65,65),P(65,65),UNEW(65,65),VNEW(65,65),
13. 1 PNEW(65,65),UOLD(65,65),VOLD(65,65),POLD(65,65),
14. 2 CU(65,65),CV(65,65),Z(65,65),H(65,65),PSI(65,65)
15. REAL MFS100,MFS200,MFS300
16. C
17. C NOTE BELOW THAT TWO DELTA T (TDT) IS SET TO DT ON THE FIRST
18. C CYCLE AFTER WHICH IT IS RESET TO DT+DT.
19. C
20. DT = 90.
21. TDT = DT
22. C
23. DX = 1.E5
24. DY = 1.E5
25. A = 1.E6
26. ALPHA = .001
27. ITMAX = 400
28. MPRINT = 100
29. M = 64
30. N = 64
31. MP1 = M+1
32. NP1 = N+1
33. EL = FLOAT(N)*DX
34. PI = 4.*ATAN(1.)
35. TPI = PI+PI
36. DI = TPI/FLOAT(M)
37. DJ = TPI/FLOAT(N)
38. PCF = PI*PI*A*A/(EL*EL)
39. C
40. C INITIAL VALUES OF THE STREAM FUNCTION AND P
41. C
42. DO 50 J=1,NP1
43. DO 50 I=1,MP1
44. PSI(I,J) = A*SIN((FLOAT(I)-.5)*DI)*SIN((FLOAT(J)-.5)*DJ)
45. P(I,J) = PCF*(COS(2.*FLOAT(I-1)*DI)
46. 1 +COS(2.*FLOAT(J-1)*DJ))+50000.
47. 50 CONTINUE
48. C
49. C INITIALIZE VELOCITIES
50. C
51. DO 60 J=1,N
52. DO 60 I=1,M
53. U(I+1,J) = -(PSI(I+1,J+1)-PSI(I+1,J))/DY
54. V(I,J+1) = (PSI(I+1,J+1)-PSI(I,J+1))/DX
55. 60 CONTINUE
56. C
57. C PERIODIC CONTINUATION
58. C
59. DO 70 J=1,N
60. U(1,J) = U(M+1,J)
61. V(M+1,J+1) = V(1,J+1)
62. 70 CONTINUE
63. DO 75 I=1,M
64. U(I+1,N+1) = U(I+1,1)
65. V(I,1) = V(I,N+1)
66. 75 CONTINUE
67. U(1,N+1) = U(M+1,1)
68. V(M+1,1) = V(1,N+1)
69. DO 86 J=1,NP1
70. DO 86 I=1,MP1
71. UOLD(I,J) = U(I,J)
72. VOLD(I,J) = V(I,J)
73. POLD(I,J) = P(I,J)
74. 86 CONTINUE
75. C
76. C PRINT INITIAL VALUES
77. C
78. WRITE(6,390) N,M,DX,DY,DT,ALPHA
79. 390 FORMAT("1NUMBER OF POINTS IN THE X DIRECTION",I8,/
80. 1 " NUMBER OF POINTS IN THE Y DIRECTION",I8,/
81. 2 " GRID SPACING IN THE X DIRECTION ",F8.0,/
82. 3 " GRID SPACING IN THE Y DIRECTION ",F8.0,/
83. 4 " TIME STEP ",F8.0,/
84. 5 " TIME FILTER PARAMETER ",F8.3)
85. MNMIN = MIN0(M,N)
86. WRITE(6,391) (P(I,I),I=1,MNMIN)
87. 391 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF P ' //,(8E15.6))
88. WRITE(6,392) (U(I,I),I=1,MNMIN)
89. 392 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF U ' //,(8E15.6))
90. WRITE(6,393) (V(I,I),I=1,MNMIN)
91. 393 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF V ' //,(8E15.6))
92. TSTART = SECOND(DUM)
93. T300 = 1.
94. TIME = 0.
95. NCYCLE = 0
96. 90 NCYCLE = NCYCLE + 1
97. C
98. C COMPUTE CAPITAL U, CAPITAL V, Z AND H
99. C
100. FSDX = 4./DX
101. FSDY = 4./DY
102. T100 = SECOND(DUM)
103. DO 100 J=1,N
104. DO 100 I=1,M
105. CU(I+1,J) = .5*(P(I+1,J)+P(I,J))*U(I+1,J)
106. CV(I,J+1) = .5*(P(I,J+1)+P(I,J))*V(I,J+1)
107. Z(I+1,J+1) = (FSDX*(V(I+1,J+1)-V(I,J+1))-FSDY*(U(I+1,J+1)
108. 1 -U(I+1,J)))/(P(I,J)+P(I+1,J)+P(I+1,J+1)+P(I,J+1))
109. H(I,J) = P(I,J)+.25*(U(I+1,J)*U(I+1,J)+U(I,J)*U(I,J)
110. 1 +V(I,J+1)*V(I,J+1)+V(I,J)*V(I,J))
111. 100 CONTINUE
112. T100 = SECOND(DUM)-T100
113. C
114. C PERIODIC CONTINUATION
115. C
116. DO 110 J=1,N
117. CU(1,J) = CU(M+1,J)
118. CV(M+1,J+1) = CV(1,J+1)
119. Z(1,J+1) = Z(M+1,J+1)
120 H(M+1,J) = H(1,J)
121. 110 CONTINUE
122. DO 115 I=1,M
123. CU(I+1,N+1) = CU(I+1,1)
124. CV(I,1) = CV(I,N+1)
125. Z(I+1,1) = Z(I+1,N+1)
126. H(I,N+1) = H(I,1)
127. 115 CONTINUE
128. CU(1,N+1) = CU(M+1,1)
129. CV(M+1,1) = CV(1,N+1)
130. Z(1,1) = Z(M+1,N+1)
131. H(M+1,N+1) = H(1,1)
132. C
133. C COMPUTE NEW VALUES U,V AND P
134. C
135. TDTS8 = TDT/8.
136. TDTSDX = TDT/DX
137. TDTSDY = TDT/DY
138. T200 = SECOND(DUM)
139. DO 200 J=1,N
140. DO 200 I=1,M
141. UNEW(I+1,J) = UOLD(I+1,J)+
142. 1 TDTS8*(Z(I+1,J+1)+Z(I+1,J))*(CV(I+1,J+1)+CV(I,J+1)+CV(I,J)
143. 2 +CV(I+1,J))-TDTSDX*(H(I+1,J)-H(I,J))
144. VNEW(I,J+1) = VOLD(I,J+1)-TDTS8*(Z(I+1,J+1)+Z(I,J+1))
145. 1 *(CU(I+1,J+1)+CU(I,J+1)+CU(I,J)+CU(I+1,J))
146. 2 -TDTSDY*(H(I,J+1)-H(I,J))
147. PNEW(I,J) = POLD(I,J)-TDTSDX*(CU(I+1,J)-CU(I,J))
148. 1 -TDTSDY*(CV(I,J+1)-CV(I,J))
149. 200 CONTINUE
150. T200 = SECOND(DUM)-T200
151. C
152. C PERIODIC CONTINUATION
153. C
154. DO 210 J=1,N
155. UNEW(1,J) = UNEW(M+1,J)
156. VNEW(M+1,J+1) = VNEW(1,J+1)
157. PNEW(M+1,J) = PNEW(1,J)
158. 210 CONTINUE
159. DO 215 I=1,M
160. UNEW(I+1,N+1) = UNEW(I+1,1)
161. VNEW(I,1) = VNEW(I,N+1)
162. PNEW(I,N+1) = PNEW(I,1)
163. 215 CONTINUE
164. UNEW(1,N+1) = UNEW(M+1,1)
165. VNEW(M+1,1) = VNEW(1,N+1)
166. PNEW(M+1,N+1) = PNEW(1,1)
167. IF(NCYCLE .GT. ITMAX) CALL EXIT
168. TIME = TIME + DT
169. IF(MOD(NCYCLE,MPRINT) .NE. 0) GO TO 370
170. PTIME = TIME/3600.
171. WRITE(6,350) NCYCLE,PTIME
172. 350 FORMAT(//' CYCLE NUMBER',I5,' MODEL TIME IN HOURS', F6.2)
173. WRITE(6,355) (PNEW(I,I),I=1,MNMIN)
174. 355 FORMAT(/' DIAGONAL ELEMENTS OF P ' //(8E15.6))
175. WRITE(6,360) (UNEW(I,I),I=1,MNMIN)
176. 360 FORMAT(/' DIAGONAL ELEMENTS OF U ' //(8E15.6))
177. WRITE(6,365) (VNEW(I,I),I=1,MNMIN)
178. 365 FORMAT(/' DIAGONAL ELEMENTS OF V ' //(8E15.6))
179. MFS100 = 24.*M*N/T100/1.E6
180. MFS200 = 26.*M*N/T200/1.E6
181. MFS300 = 15.*M*N/T300/1.E6
182. CTIME = SECOND(DUM)-TSTART
183. TCYC = CTIME/FLOAT(NCYCLE)
184. WRITE(6,375) NCYCLE,CTIME,TCYC,T310,MFS310,T200,MFS200,T300,MFS300
185. 375 FORMAT(' CYCLE NUMBER',I5,' TOTAL COMPUTER TIME', E15.6,
186. 1 ' TIME PER CYCLE', E15.6, /
187. 2 ' TIME AND MEGAFLOPS FOR LOOP 310 ', E15.6,F6.1/
188. 3 ' TIME AND MEGAFLOPS FOR LOOP 200 ', E15.6,F6.1/
189. 4 ' TIME AND MEGAFLOPS FOR LOOP 300 ', E15.6,F6.1/ )
190.
191. C
192. C TIME SMOOTHING AND UPDATE FOR NEXT CYCLE
193. C
194. 370 IF(NCYCLE .LE. 1) GO TO 310
195. T300 = SECOND(DUM)
196. DO 300 J=1,N
197. DO 300 I=1,M
198. UOLD(I,J) = U(I,J)+ALPHA*(UNEW(I,J)-2.*U(I,J)+UOLD(I,J))
199. VOLD(I,J) = V(I,J)+ALPHA*(VNEW(I,J)-2.*V(I,J)+VOLD(I,J))
200. POLD(I,J) = P(I,J)+ALPHA*(PNEW(I,J)-2.*P(I,J)+POLD(I,J))
201. U(I,J) = UNEW(I,J)
202. V(I,J) = VNEW(I,J)
203. P(I,J) = PNEW(I,J)
204. 300 CONTINUE
205. T300 = SECOND(DUM)-T300
206. C
207. C PERIODIC CONTINUATION
208. C
209. DO 320 J=1,N
210. UOLD(M+1,J) = UOLD(1,J)
211. VOLD(M+1,J) = VOLD(1,J)
212. POLD(M+1,J) = POLD(1,J)
213. U(M+1,J) = U(1,J)
214. V(M+1,J) = V(1,J)
215. P(M+1,J) = P(1,J)
216. 320 CONTINUE
217. DO 325 I=1,M
218. UOLD(I,N+1) = UOLD(I,1)
219. VOLD(I,N+1) = VOLD(I,1)
220. POLD(I,N+1) = POLD(I,1)
221. U(I,N+1) = U(I,1)
222. V(I,N+1) = V(I,1)
223. P(I,N+1) = P(I,1)
224. 325 CONTINUE
225. UOLD(M+1,N+1) = UOLD(1,1)
226. VOLD(M+1,N+1) = VOLD(1,1)
227. POLD(M+1,N+1) = POLD(1,1)
228. U(M+1,N+1) = U(1,1)
229. V(M+1,N+1) = V(1,1)
230. P(M+1,N+1) = P(1,1)
231. GO TO 90
232. 310 TDT = TDT+TDT
233. DO 400 J=1,NP1
234. DO 400 I=1,MP1
235. UOLD(I,J) = U(I,J)
236. VOLD(I,J) = V(I,J)
237. POLD(I,J) = P(I,J)
238. U(I,J) = UNEW(I,J)
239. V(I,J) = VNEW(I,J)
240. P(I,J) = PNEW(I,J)
241. 400 CONTINUE
242. GO TO 90
243. END
244. function second(dum)
245. real time
246. call cpu_time(time)
247. second=time
248. end
Line 27: ITMAX - the number of time steps
Line 57: Periodic continuation - The last column of points is the same
as the first column, and the top row of points is the same as the
bottom row. Thus, the values "wrap around" back to the beginning, like
on a sphere. The first column is not set initially - it gets set after
all of the other columns are calculated, then the last value gets
duplicated into the first column. Similarly for the rows. This is an
effort to mimic a sphere on what's actually a rectangle.
If you have questions about this document, please contact
SCD Customer Support.
You can also reach us by telephone 24 hours a day, seven days a week at
303-497-1278.
Additional contact methods:
consult1@ucar.edu
and during
business hours
in NCAR Mesa Lab Suite 39.
© Copyright 2006. University Corporation for Atmospheric
Research (UCAR). All Rights Reserved.
|
