|
Example 3 - shallow3.f - Code with portable timing instrumentation
added. Serial.
Download shallow3.f code.
1. PROGRAM SHALOW
2. * In this version, shallow3.f, the code is still serial
3. * but has been brought up to modern Fortran constructs and
4. * uses the Fortran 90 intrinsic timing routine system_clock,
5. * which is portable
6. * This can be compiled on the IBM SP using
7. * xlf90 -g -o shallow3 -qfixed=132 -qsclk=micro shallow3.f
8.
9. C
10. C BENCHMARK WEATHER PREDICTION PROGRAM FOR COMPARING THE
11. C PREFORMANCE OF CURRENT SUPERCOMPUTERS. THE MODEL IS
12. C BASED OF THE PAPER - THE DYNAMICS OF FINITE-DIFFERENCE
13. C MODELS OF THE SHALLOW-WATER EQUATIONS, BY ROBERT SADOURNY
14. C J. ATM. SCIENCES, VOL 32, NO 4, APRIL 1975.
15. C
16. C CODE BY PAUL N. SWARZTRAUBER, NATIONAL CENTER FOR
17. C ATMOSPHERIC RESEARCH, BOULDER, CO, OCTOBER 1984.
18. C Modifed by Juliana Rew, NCAR, January 2006
19.
20. DIMENSION U(65,65),V(65,65),P(65,65),UNEW(65,65),VNEW(65,65),
21. 1 PNEW(65,65),UOLD(65,65),VOLD(65,65),POLD(65,65),
22. 2 CU(65,65),CV(65,65),Z(65,65),H(65,65),PSI(65,65)
23. * REAL MFS100,MFS200,MFS300
24. double precision MFS100,MFS200,MFS300
25. C
26. C NOTE BELOW THAT TWO DELTA T (TDT) IS SET TO DT ON THE FIRST
27. C CYCLE AFTER WHICH IT IS RESET TO DT+DT.
28. C
29. * external second
30. integer c1,c2,r,max
31. double precision T100,T200,T300,TSTART,CTIME,TCYC
32. DT = 90.
33. TDT = DT
34. C
35. DX = 1.E5
36. DY = 1.E5
37. A = 1.E6
38. ALPHA = .001
39. ITMAX = 400
40. MPRINT = 100
41. M = 64
42. N = 64
43. MP1 = M+1
44. NP1 = N+1
45. EL = FLOAT(N)*DX
46. PI = 4.*ATAN(1.)
47. TPI = PI+PI
48. DI = TPI/FLOAT(M)
49. DJ = TPI/FLOAT(N)
50. PCF = PI*PI*A*A/(EL*EL)
51. C
52. C INITIAL VALUES OF THE STREAM FUNCTION AND P
53. C
54. DO J=1,NP1
55. DO I=1,MP1
56. PSI(I,J) = A*SIN((FLOAT(I)-.5)*DI)*SIN((FLOAT(J)-.5)*DJ)
57. P(I,J) = PCF*(COS(2.*FLOAT(I-1)*DI)
58. 1 +COS(2.*FLOAT(J-1)*DJ))+50000.
59. END DO
60. 50 END DO
61. C
62. C INITIALIZE VELOCITIES
63. C
64. DO J=1,N
65. DO I=1,M
66. U(I+1,J) = -(PSI(I+1,J+1)-PSI(I+1,J))/DY
67. V(I,J+1) = (PSI(I+1,J+1)-PSI(I,J+1))/DX
68. END DO
69. 60 END DO
70. C
71. C PERIODIC CONTINUATION
72. C
73. DO J=1,N
74. U(1,J) = U(M+1,J)
75. V(M+1,J+1) = V(1,J+1)
76. 70 END DO
77. DO I=1,M
78. U(I+1,N+1) = U(I+1,1)
79. V(I,1) = V(I,N+1)
80. 75 END DO
81. U(1,N+1) = U(M+1,1)
82. V(M+1,1) = V(1,N+1)
83. DO J=1,NP1
84. DO I=1,MP1
85. UOLD(I,J) = U(I,J)
86. VOLD(I,J) = V(I,J)
87. POLD(I,J) = P(I,J)
88. END DO
89. 86 END DO
90. C
91. C PRINT INITIAL VALUES
92. C
93. WRITE(6,390) N,M,DX,DY,DT,ALPHA
94. 390 FORMAT("1NUMBER OF POINTS IN THE X DIRECTION",I8,/
95. 1 " NUMBER OF POINTS IN THE Y DIRECTION",I8,/
96. 2 " GRID SPACING IN THE X DIRECTION ",F8.0,/
97. 3 " GRID SPACING IN THE Y DIRECTION ",F8.0,/
98. 4 " TIME STEP ",F8.0,/
99. 5 " TIME FILTER PARAMETER ",F8.3)
100. MNMIN = MIN0(M,N)
101. WRITE(6,391) (P(I,I),I=1,MNMIN)
102. 391 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF P ' //,(8E15.6))
103. WRITE(6,392) (U(I,I),I=1,MNMIN)
104. 392 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF U ' //,(8E15.6))
105. WRITE(6,393) (V(I,I),I=1,MNMIN)
106. 393 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF V ' //,(8E15.6))
107. * TSTART = SECOND(DUM)
108. call system_clock (count=c1, count_rate=r, count_max=max)
109. TSTART = c1
110. T300 = 1.
111. TIME = 0.
112. NCYCLE = 0
113. * 90 NCYCLE = NCYCLE + 1
114. DO ncycle=1,itmax
115. C
116. C COMPUTE CAPITAL U, CAPITAL V, Z AND H
117. C
118. FSDX = 4./DX
119. FSDY = 4./DY
120 * T100 = SECOND(DUM)
121. call system_clock(count=c1, count_rate=r,count_max=max)
122. T100 = c1
123. DO J=1,N
124. DO I=1,M
125. CU(I+1,J) = .5*(P(I+1,J)+P(I,J))*U(I+1,J)
126. CV(I,J+1) = .5*(P(I,J+1)+P(I,J))*V(I,J+1)
127. Z(I+1,J+1) =(FSDX*(V(I+1,J+1)-V(I,J+1))-FSDY*(U(I+1,J+1)
128. 1 -U(I+1,J)))/(P(I,J)+P(I+1,J)+P(I+1,J+1)+P(I,J+1))
129. H(I,J) = P(I,J)+.25*(U(I+1,J)*U(I+1,J)+U(I,J)*U(I,J)
130. 1 +V(I,J+1)*V(I,J+1)+V(I,J)*V(I,J))
131. END DO
132. 100 END DO
133. * T100 = SECOND(DUM)-T100
134. call system_clock(count=c2,count_rate=r,count_max=max)
135. T100 = dble(c2-T100)/dble(r)
136. C
137. C PERIODIC CONTINUATION
138. C
139. DO J=1,N
140. CU(1,J) = CU(M+1,J)
141. CV(M+1,J+1) = CV(1,J+1)
142. Z(1,J+1) = Z(M+1,J+1)
143. H(M+1,J) = H(1,J)
144. 110 END DO
145. DO I=1,M
146. CU(I+1,N+1) = CU(I+1,1)
147. CV(I,1) = CV(I,N+1)
148. Z(I+1,1) = Z(I+1,N+1)
149. H(I,N+1) = H(I,1)
150. 115 END DO
151. CU(1,N+1) = CU(M+1,1)
152. CV(M+1,1) = CV(1,N+1)
153. Z(1,1) = Z(M+1,N+1)
154. H(M+1,N+1) = H(1,1)
155. C
156. C COMPUTE NEW VALUES U,V AND P
157. C
158. TDTS8 = TDT/8.
159. TDTSDX = TDT/DX
160. TDTSDY = TDT/DY
161. * T200 = SECOND(DUM)
162. call system_clock(count=c1, count_rate=r, count_max=max)
163. T200 = c1
164. DO J=1,N
165. DO I=1,M
166. UNEW(I+1,J) = UOLD(I+1,J)+
167. 1 TDTS8*(Z(I+1,J+1)+Z(I+1,J))*(CV(I+1,J+1)+CV(I,J+1)+CV(I,J)
168. 2 +CV(I+1,J))-TDTSDX*(H(I+1,J)-H(I,J))
169. VNEW(I,J+1) = VOLD(I,J+1)-TDTS8*(Z(I+1,J+1)+Z(I,J+1))
170. 1 *(CU(I+1,J+1)+CU(I,J+1)+CU(I,J)+CU(I+1,J))
171. 2 -TDTSDY*(H(I,J+1)-H(I,J))
172. PNEW(I,J) = POLD(I,J)-TDTSDX*(CU(I+1,J)-CU(I,J))
173. 1 -TDTSDY*(CV(I,J+1)-CV(I,J))
174. END DO
175. 200 END DO
176. * T200 = SECOND(DUM)-T200
177. call system_clock(count=c2, count_rate=r, count_max=max)
178. T200 = dble(c2 -T200)/dble(r)
179. C
180. C PERIODIC CONTINUATION
181. C
182. DO J=1,N
183. UNEW(1,J) = UNEW(M+1,J)
184. VNEW(M+1,J+1) = VNEW(1,J+1)
185. PNEW(M+1,J) = PNEW(1,J)
186. 210 END DO
187. DO I=1,M
188. UNEW(I+1,N+1) = UNEW(I+1,1)
189. VNEW(I,1) = VNEW(I,N+1)
190. PNEW(I,N+1) = PNEW(I,1)
191. 215 END DO
192. UNEW(1,N+1) = UNEW(M+1,1)
193. VNEW(M+1,1) = VNEW(1,N+1)
194. PNEW(M+1,N+1) = PNEW(1,1)
195. * jr IF(NCYCLE .GT. ITMAX) CALL EXIT
196.
197. TIME = TIME + DT
198. * jr IF(MOD(NCYCLE,MPRINT) .NE. 0) GO TO 370
199. IF(MOD(NCYCLE,MPRINT) .EQ. 0) then
200. PTIME = TIME/3600.
201. WRITE(6,350) NCYCLE,PTIME
202. 350 FORMAT(//' CYCLE NUMBER',I5,' MODEL TIME IN HOURS', F6.2)
203. WRITE(6,355) (PNEW(I,I),I=1,MNMIN)
204. 355 FORMAT(/' DIAGONAL ELEMENTS OF P ' //(8E15.6))
205. WRITE(6,360) (UNEW(I,I),I=1,MNMIN)
206. 360 FORMAT(/' DIAGONAL ELEMENTS OF U ' //(8E15.6))
207. WRITE(6,365) (VNEW(I,I),I=1,MNMIN)
208. 365 FORMAT(/' DIAGONAL ELEMENTS OF V ' //(8E15.6))
209. * jr added MFS310--don't know what actual mult factor should be
210. * jr changed divide by 1 million to divide by 100K since system_clock
211. * jr resolution is millisec rather than cpu_time's 10 millisec
212. MFS310 = 24.*M*N/T310/1.D5
213. * MFS100 = 24.*M*N/T100/1.E6
214. MFS100 = 24.*M*N/T100/1.D5
215. * MFS200 = 26.*M*N/T200/1.E6
216. MFS200 = 26.*M*N/T200/1.D5
217. * MFS300 = 15.*M*N/T300/1.E6
218. MFS300 = 15.*M*N/T300/1.D5
219. * CTIME = SECOND(DUM)-TSTART
220. call system_clock(count=c2, count_rate=r,count_max=max)
221. CTIME = dble(c2 - TSTART)/dble(r)
222. TCYC = CTIME/FLOAT(NCYCLE)
223. WRITE(6,375) NCYCLE,CTIME,TCYC,T310,MFS310,T200,MFS200,T300,MFS300
224. * 375 FORMAT(' CYCLE NUMBER',I5,' TOTAL COMPUTER TIME', E15.6,
225. * 1 ' TIME PER CYCLE', E15.6, /
226. * 2 ' TIME AND MEGAFLOPS FOR LOOP 310 ', E15.6,F6.1/
227. * 3 ' TIME AND MEGAFLOPS FOR LOOP 200 ', E15.6,F6.1/
228. * 4 ' TIME AND MEGAFLOPS FOR LOOP 300 ', E15.6,F6.1/ )
229. 375 FORMAT(' CYCLE NUMBER',I5,' TOTAL COMPUTER TIME', D15.6,
230. 1 ' TIME PER CYCLE', D15.6, /
231. 2 ' TIME AND MEGAFLOPS FOR LOOP 310 ', D15.6,2x,D6.1/
232. 3 ' TIME AND MEGAFLOPS FOR LOOP 200 ', D15.6,2x,D6.1/
233. 4 ' TIME AND MEGAFLOPS FOR LOOP 300 ', D15.6,2x,D6.1/ )
234. endif
235. C
236. C TIME SMOOTHING AND UPDATE FOR NEXT CYCLE
237. C
238. * 370 IF(NCYCLE .LE. 1) GO TO 310
239. 370 IF(NCYCLE .GT. 1) then
240. * T300 = SECOND(DUM)
241. call system_clock(count=c1,count_rate=r,count_max=max)
242. T300 = c1
243.
244. DO J=1,N
245. DO I=1,M
246. UOLD(I,J) = U(I,J)+ALPHA*(UNEW(I,J)-2.*U(I,J)+UOLD(I,J))
247. VOLD(I,J) = V(I,J)+ALPHA*(VNEW(I,J)-2.*V(I,J)+VOLD(I,J))
248. POLD(I,J) = P(I,J)+ALPHA*(PNEW(I,J)-2.*P(I,J)+POLD(I,J))
249. U(I,J) = UNEW(I,J)
250. V(I,J) = VNEW(I,J)
251. P(I,J) = PNEW(I,J)
252. END DO
253. 300 END DO
254. * T300 = SECOND(DUM)-T300
255. call system_clock(count=c2,count_rate=r,count_max=max)
256. T300 = dble(c2 - T300)/dble(r)
257. C
258. C PERIODIC CONTINUATION
259. C
260. DO J=1,N
261. UOLD(M+1,J) = UOLD(1,J)
262. VOLD(M+1,J) = VOLD(1,J)
263. POLD(M+1,J) = POLD(1,J)
264. U(M+1,J) = U(1,J)
265. V(M+1,J) = V(1,J)
266. P(M+1,J) = P(1,J)
267. 320 END DO
268. DO I=1,M
269. UOLD(I,N+1) = UOLD(I,1)
270. VOLD(I,N+1) = VOLD(I,1)
271. POLD(I,N+1) = POLD(I,1)
272. U(I,N+1) = U(I,1)
273. V(I,N+1) = V(I,1)
274. P(I,N+1) = P(I,1)
275. 325 END DO
276. UOLD(M+1,N+1) = UOLD(1,1)
277. VOLD(M+1,N+1) = VOLD(1,1)
278. POLD(M+1,N+1) = POLD(1,1)
279. U(M+1,N+1) = U(1,1)
280. V(M+1,N+1) = V(1,1)
281. P(M+1,N+1) = P(1,1)
282. else
283. 310 TDT = TDT+TDT
284. call system_clock(count=c1, count_rate=r, count_max=max)
285. T310 = c1
286. DO J=1,NP1
287. DO I=1,MP1
288. UOLD(I,J) = U(I,J)
289. VOLD(I,J) = V(I,J)
290. POLD(I,J) = P(I,J)
291. U(I,J) = UNEW(I,J)
292. V(I,J) = VNEW(I,J)
293. P(I,J) = PNEW(I,J)
294. END DO
295. 400 END DO
296. call system_clock(count=c2,count_rate=r,count_max=max)
297. T310 = dble(c2 - T310)/dble(r)
298. endif
299. * end big loop 90
300. End do
301.
302. END
303.
304. * function second(dum)
305. * real time
306. * call cpu_time(time)
307. * second=time
308. * end
Line 304. The second function was written to give the elapsed time in
seconds. It calls the IBM XLF routine cpu_time, so is not portable to
other platforms. It was replaced throughout the code with the portable
Fortran 90 intrinsic, system_clock.
Lines 108, 121, 134, 162, 177, 220, 241, 255, 284, and 296 - calls to system_clock were placed before and after each loop that we wished to time.
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.
|
