|
Example 2 - shallow2.f - code modernized to Fortran90 standard.
Serial.
Download shallow2.f code.
1. PROGRAM SHALOW
2. * In this version, shallow2.f, the code has been brought to
3. * more modern Fortran90 constructs. Gotos have been removed.
4. * This is still a serial code. Timing routines used are for
5. * IBM, so are not portable.
6. * This can be compiled on the IBM SP using:
7. * xlf90 -g -o shallow2 -qfixed=132 shallow2.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 Modified 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.
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. DT = 90.
31. TDT = DT
32. C
33. DX = 1.E5
34. DY = 1.E5
35. A = 1.E6
36. ALPHA = .001
37. ITMAX = 400
38. MPRINT = 100
39. M = 64
40. N = 64
41. MP1 = M+1
42. NP1 = N+1
43. EL = FLOAT(N)*DX
44. PI = 4.*ATAN(1.)
45. TPI = PI+PI
46. DI = TPI/FLOAT(M)
47. DJ = TPI/FLOAT(N)
48. PCF = PI*PI*A*A/(EL*EL)
49. C
50. C INITIAL VALUES OF THE STREAM FUNCTION AND P
51. C
52. DO J=1,NP1
53. DO I=1,MP1
54. PSI(I,J) = A*SIN((FLOAT(I)-.5)*DI)*SIN((FLOAT(J)-.5)*DJ)
55. P(I,J) = PCF*(COS(2.*FLOAT(I-1)*DI)
56. 1 +COS(2.*FLOAT(J-1)*DJ))+50000.
57. END DO
58. 50 END DO
59. C
60. C INITIALIZE VELOCITIES
61. C
62. DO J=1,N
63. DO I=1,M
64. U(I+1,J) = -(PSI(I+1,J+1)-PSI(I+1,J))/DY
65. V(I,J+1) = (PSI(I+1,J+1)-PSI(I,J+1))/DX
66. END DO
67. 60 END DO
68. C
69. C PERIODIC CONTINUATION
70. C
71. DO J=1,N
72. U(1,J) = U(M+1,J)
73. V(M+1,J+1) = V(1,J+1)
74. 70 END DO
75. DO I=1,M
76. U(I+1,N+1) = U(I+1,1)
77. V(I,1) = V(I,N+1)
78. 75 END DO
79. U(1,N+1) = U(M+1,1)
80. V(M+1,1) = V(1,N+1)
81. DO J=1,NP1
82. DO I=1,MP1
83. UOLD(I,J) = U(I,J)
84. VOLD(I,J) = V(I,J)
85. POLD(I,J) = P(I,J)
86. END DO
87. 86 END DO
88. C
89. C PRINT INITIAL VALUES
90. C
91. WRITE(6,390) N,M,DX,DY,DT,ALPHA
92. 390 FORMAT("1NUMBER OF POINTS IN THE X DIRECTION",I8,/
93. 1 " NUMBER OF POINTS IN THE Y DIRECTION",I8,/
94. 2 " GRID SPACING IN THE X DIRECTION ",F8.0,/
95. 3 " GRID SPACING IN THE Y DIRECTION ",F8.0,/
96. 4 " TIME STEP ",F8.0,/
97. 5 " TIME FILTER PARAMETER ",F8.3)
98. MNMIN = MIN0(M,N)
99. WRITE(6,391) (P(I,I),I=1,MNMIN)
100. 391 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF P ' //,(8E15.6))
101. WRITE(6,392) (U(I,I),I=1,MNMIN)
102. 392 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF U ' //,(8E15.6))
103. WRITE(6,393) (V(I,I),I=1,MNMIN)
104. 393 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF V ' //,(8E15.6))
105. TSTART = SECOND(DUM)
106. T300 = 1.
107. TIME = 0.
108. NCYCLE = 0
109. * 90 NCYCLE = NCYCLE + 1
110. DO ncycle=1,itmax
111. C
112. C COMPUTE CAPITAL U, CAPITAL V, Z AND H
113. C
114. FSDX = 4./DX
115. FSDY = 4./DY
116. T100 = SECOND(DUM)
117. DO J=1,N
118. DO I=1,M
119. CU(I+1,J) = .5*(P(I+1,J)+P(I,J))*U(I+1,J)
120 CV(I,J+1) = .5*(P(I,J+1)+P(I,J))*V(I,J+1)
121. Z(I+1,J+1) = (FSDX*(V(I+1,J+1)-V(I,J+1))-FSDY*(U(I+1,J+1)
122. 1 -U(I+1,J)))/(P(I,J)+P(I+1,J)+P(I+1,J+1)+P(I,J+1))
123. H(I,J) = P(I,J)+.25*(U(I+1,J)*U(I+1,J)+U(I,J)*U(I,J)
124. 1 +V(I,J+1)*V(I,J+1)+V(I,J)*V(I,J))
125. END DO
126. 100 END DO
127. T100 = SECOND(DUM)-T100
128. C
129. C PERIODIC CONTINUATION
130. C
131. DO J=1,N
132. CU(1,J) = CU(M+1,J)
133. CV(M+1,J+1) = CV(1,J+1)
134. Z(1,J+1) = Z(M+1,J+1)
135. H(M+1,J) = H(1,J)
136. 110 END DO
137. DO I=1,M
138. CU(I+1,N+1) = CU(I+1,1)
139. CV(I,1) = CV(I,N+1)
140. Z(I+1,1) = Z(I+1,N+1)
141. H(I,N+1) = H(I,1)
142. 115 END DO
143. CU(1,N+1) = CU(M+1,1)
144. CV(M+1,1) = CV(1,N+1)
145. Z(1,1) = Z(M+1,N+1)
146. H(M+1,N+1) = H(1,1)
147. C
148. C COMPUTE NEW VALUES U,V AND P
149. C
150. TDTS8 = TDT/8.
151. TDTSDX = TDT/DX
152. TDTSDY = TDT/DY
153. T200 = SECOND(DUM)
154. DO J=1,N
155. DO I=1,M
156. UNEW(I+1,J) = UOLD(I+1,J)+
157. 1 TDTS8*(Z(I+1,J+1)+Z(I+1,J))*(CV(I+1,J+1)+CV(I,J+1)+CV(I,J)
158. 2 +CV(I+1,J))-TDTSDX*(H(I+1,J)-H(I,J))
159. VNEW(I,J+1) = VOLD(I,J+1)-TDTS8*(Z(I+1,J+1)+Z(I,J+1))
160. 1 *(CU(I+1,J+1)+CU(I,J+1)+CU(I,J)+CU(I+1,J))
161. 2 -TDTSDY*(H(I,J+1)-H(I,J))
162. PNEW(I,J) = POLD(I,J)-TDTSDX*(CU(I+1,J)-CU(I,J))
163. 1 -TDTSDY*(CV(I,J+1)-CV(I,J))
164. END DO
165. 200 END DO
166. T200 = SECOND(DUM)-T200
167. C
168. C PERIODIC CONTINUATION
169. C
170. DO J=1,N
171. UNEW(1,J) = UNEW(M+1,J)
172. VNEW(M+1,J+1) = VNEW(1,J+1)
173. PNEW(M+1,J) = PNEW(1,J)
174. 210 END DO
175. DO I=1,M
176. UNEW(I+1,N+1) = UNEW(I+1,1)
177. VNEW(I,1) = VNEW(I,N+1)
178. PNEW(I,N+1) = PNEW(I,1)
179. 215 END DO
180. UNEW(1,N+1) = UNEW(M+1,1)
181. VNEW(M+1,1) = VNEW(1,N+1)
182. PNEW(M+1,N+1) = PNEW(1,1)
183. * jr IF(NCYCLE .GT. ITMAX) CALL EXIT
184.
185. TIME = TIME + DT
186. * jr IF(MOD(NCYCLE,MPRINT) .NE. 0) GO TO 370
187. IF(MOD(NCYCLE,MPRINT) .EQ. 0) then
188.
189. PTIME = TIME/3600.
190. WRITE(6,350) NCYCLE,PTIME
191. 350 FORMAT(//' CYCLE NUMBER',I5,' MODEL TIME IN HOURS', F6.2)
192. WRITE(6,355) (PNEW(I,I),I=1,MNMIN)
193. 355 FORMAT(/' DIAGONAL ELEMENTS OF P ' //(8E15.6))
194. WRITE(6,360) (UNEW(I,I),I=1,MNMIN)
195. 360 FORMAT(/' DIAGONAL ELEMENTS OF U ' //(8E15.6))
196. WRITE(6,365) (VNEW(I,I),I=1,MNMIN)
197. 365 FORMAT(/' DIAGONAL ELEMENTS OF V ' //(8E15.6))
198. MFS100 = 24.*M*N/T100/1.D6
199. MFS200 = 26.*M*N/T200/1.D6
200. MFS300 = 15.*M*N/T300/1.D6
201. CTIME = SECOND(DUM)-TSTART
202. TCYC = CTIME/FLOAT(NCYCLE)
203. WRITE(6,375) NCYCLE,CTIME,TCYC,T310,MFS310,T200,MFS200,T300,MFS300
204. 375 FORMAT(' CYCLE NUMBER',I5,' TOTAL COMPUTER TIME', E15.6,
205. 1 ' TIME PER CYCLE', E15.6, /
206. 2 ' TIME AND MEGAFLOPS FOR LOOP 310 ', E15.6,F6.1/
207. 3 ' TIME AND MEGAFLOPS FOR LOOP 200 ', E15.6,F6.1/
208. 4 ' TIME AND MEGAFLOPS FOR LOOP 300 ', E15.6,F6.1/ )
209. endif
210. C
211. C TIME SMOOTHING AND UPDATE FOR NEXT CYCLE
212. C
213. * 370 IF(NCYCLE .LE. 1) GO TO 310
214. 370 IF(NCYCLE .GT. 1) then
215. T300 = SECOND(DUM)
216. DO J=1,N
217. DO I=1,M
218. UOLD(I,J) = U(I,J)+ALPHA*(UNEW(I,J)-2.*U(I,J)+UOLD(I,J))
219. VOLD(I,J) = V(I,J)+ALPHA*(VNEW(I,J)-2.*V(I,J)+VOLD(I,J))
220. POLD(I,J) = P(I,J)+ALPHA*(PNEW(I,J)-2.*P(I,J)+POLD(I,J))
221. U(I,J) = UNEW(I,J)
222. V(I,J) = VNEW(I,J)
223. P(I,J) = PNEW(I,J)
224. END DO
225. 300 END DO
226. T300 = SECOND(DUM)-T300
227. C
228. C PERIODIC CONTINUATION
229. C
230. DO J=1,N
231. UOLD(M+1,J) = UOLD(1,J)
232. VOLD(M+1,J) = VOLD(1,J)
233. POLD(M+1,J) = POLD(1,J)
234. U(M+1,J) = U(1,J)
235. V(M+1,J) = V(1,J)
236. P(M+1,J) = P(1,J)
237. 320 END DO
238. DO I=1,M
239. UOLD(I,N+1) = UOLD(I,1)
240. VOLD(I,N+1) = VOLD(I,1)
241. POLD(I,N+1) = POLD(I,1)
242. U(I,N+1) = U(I,1)
243. V(I,N+1) = V(I,1)
244. P(I,N+1) = P(I,1)
245. 325 END DO
246. UOLD(M+1,N+1) = UOLD(1,1)
247. VOLD(M+1,N+1) = VOLD(1,1)
248. POLD(M+1,N+1) = POLD(1,1)
249. U(M+1,N+1) = U(1,1)
250. V(M+1,N+1) = V(1,1)
251. P(M+1,N+1) = P(1,1)
252. else
253.
243. 310 TDT = TDT+TDT
244. DO J=1,NP1
245. DO I=1,MP1
246. UOLD(I,J) = U(I,J)
247. VOLD(I,J) = V(I,J)
248. POLD(I,J) = P(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. 400 END DO
254. endif
255. *end big loop 90
256. End do
257.
258. END
259.
260. function second(dum)
261. real time
262. call cpu_time(time)
263. second=time
264. end
Line 186 - *jr IF(MOD(NCYCLE,MPRINT) .NE. 0) GO TO 370
Goto removed and replaced by do loop.
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.
|
