 |
 |
|
||
|
|
|
Example 4 - shallow4.f - Code with netCDF I/O added to facilitate data portability and visualizationDownload shallow4.f code.
1. PROGRAM SHALOW 2. * In this version, shallow4.f, initial and calculated values 3. * of U, V, and P are written to a netCDF file 4. * for later use in visualizing the results. The netCDF data 5. * management library is freely available from 6. * http://www.unidata.ucar.edu/software/netcdf 7. * This code is still serial but has been brought up to modern 8. * Fortran constructs and uses portable intrinsic Fortran 90 timing routines 9. * This can be compiled on the IBM SP using: 10. * xlf90 -qmaxmem=-1 -g -o shallow4 -qfixed=132 -qsclk=micro \ 11. * -I/usr/local/include shallow4.f -L/usr/local/lib32/r4i4 -l netcdf 12. * where the -L and -I point to local installation of netCDF 13. 14. include 'netcdf.inc' 15. C 16. C BENCHMARK WEATHER PREDICTION PROGRAM FOR COMPARING THE 17. C PREFORMANCE OF CURRENT SUPERCOMPUTERS. THE MODEL IS 18. C BASED OF THE PAPER - THE DYNAMICS OF FINITE-DIFFERENCE 19. C MODELS OF THE SHALLOW-WATER EQUATIONS, BY ROBERT SADOURNY 20. C J. ATM. SCIENCES, VOL 32, NO 4, APRIL 1975. 21. C 22. C CODE BY PAUL N. SWARZTRAUBER, NATIONAL CENTER FOR 23. C ATMOSPHERIC RESEARCH, BOULDER, CO, OCTOBER 1984. 24. C Modified by Juliana Rew, NCAR, January 2006 25. 26. integer mlen,nlen 27. 28. parameter (m_len = 65) 29. parameter (n_len = 65) 30. parameter (m = 64) 31. parameter (n = 64) 32. DIMENSION U(m_len,n_len),V(m_len,n_len),P(m_len,n_len),UNEW(m_len,n_len),VNEW(m_len,n_len), 33. 1 PNEW(m_len,n_len),UOLD(m_len,n_len),VOLD(m_len,n_len),POLD(m_len,n_len), 34. 2 CU(m_len,n_len),CV(m_len,n_len),Z(m_len,n_len),H(m_len,n_len),PSI(m_len,n_len) 35. 36. * REAL MFS100,MFS200,MFS300 37. double precision MFS100,MFS200,MFS300 38. double precision T100,T200,T300,TSTART,CTIME,TCYC 39. integer c1,c2,r,max 40. integer ncid,p_id,u_id,v_id 41. integer istart(3) 42. integer icount(3) 43. * prepare netCDF file to receive model output data 44. call netcdf_setup(m_len,n_len,ncid,p_id,u_id,v_id,istart,icount) 45. C 46. C NOTE BELOW THAT TWO DELTA T (TDT) IS SET TO DT ON THE FIRST 47. C CYCLE AFTER WHICH IT IS RESET TO DT+DT. 48. C 49. 50. DT = 90. 51. TDT = DT 52. C 53. DX = 1.E5 54. DY = 1.E5 55. A = 1.E6 56. ALPHA = .001 57. ITMAX = 400 58. MPRINT = 100 59. * M and N could be set as constants above to make it easier 60. * to change the size of the problem above 61. * M = 64 62. * N = 64 63. MP1 = M+1 64. NP1 = N+1 65. EL = FLOAT(N)*DX 66. PI = 4.*ATAN(1.) 67. TPI = PI+PI 68. DI = TPI/FLOAT(M) 69. DJ = TPI/FLOAT(N) 70. PCF = PI*PI*A*A/(EL*EL) 71. C 72. C INITIAL VALUES OF THE STREAM FUNCTION AND P 73. C 74. DO J=1,NP1 75. DO I=1,MP1 76. PSI(I,J) = A*SIN((FLOAT(I)-.5)*DI)*SIN((FLOAT(J)-.5)*DJ) 77. * PSI(I,J) = 0 78. P(I,J) = PCF*(COS(2.*FLOAT(I-1)*DI) 79. 1 +COS(2.*FLOAT(J-1)*DJ))+50000. 80. END DO 81. 50 END DO 82. C 83. C INITIALIZE VELOCITIES 84. C 85. DO J=1,N 86. DO I=1,M 87. U(I+1,J) = -(PSI(I+1,J+1)-PSI(I+1,J))/DY 88. V(I,J+1) = (PSI(I+1,J+1)-PSI(I,J+1))/DX 89. END DO 90. 60 END DO 91. C 92. C PERIODIC CONTINUATION 93. C 94. DO J=1,N 95. U(1,J) = U(M+1,J) 96. V(M+1,J+1) = V(1,J+1) 97. 70 END DO 98. DO I=1,M 99. U(I+1,N+1) = U(I+1,1) 100. V(I,1) = V(I,N+1) 101. 75 END DO 102. U(1,N+1) = U(M+1,1) 103. V(M+1,1) = V(1,N+1) 104. DO J=1,NP1 105. DO I=1,MP1 106. UOLD(I,J) = U(I,J) 107. VOLD(I,J) = V(I,J) 108. POLD(I,J) = P(I,J) 109. END DO 110. 86 END DO 111. C 112. C PRINT INITIAL VALUES 113. C 114. WRITE(6,390) N,M,DX,DY,DT,ALPHA 115. 390 FORMAT("1NUMBER OF POINTS IN THE X DIRECTION",I8,/ 116. 1 " NUMBER OF POINTS IN THE Y DIRECTION",I8,/ 117. 2 " GRID SPACING IN THE X DIRECTION ",F8.0,/ 118. 3 " GRID SPACING IN THE Y DIRECTION ",F8.0,/ 119. 4 " TIME STEP ",F8.0,/ 120 5 " TIME FILTER PARAMETER ",F8.3) 121. MNMIN = MIN0(M,N) 122. WRITE(6,391) (P(I,I),I=1,MNMIN) 123. 391 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF P ' //,(8E15.6)) 124. WRITE(6,392) (U(I,I),I=1,MNMIN) 125. 392 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF U ' //,(8E15.6)) 126. WRITE(6,393) (V(I,I),I=1,MNMIN) 127. 393 FORMAT(/' INITIAL DIAGONAL ELEMENTS OF V ' //,(8E15.6)) 128. * write intial values of p, u, and v into a netCDF file 129. 130. call my_ncwrite(ncid,p_id,istart,icount,p,m_len,n_len) 131. call my_ncwrite(ncid,u_id,istart,icount,u,m_len,n_len) 132. call my_ncwrite(ncid,v_id,istart,icount,v,m_len,n_len) 133. * TSTART = SECOND(DUM) 134. call system_clock (count=c1, count_rate=r, count_max=max) 135. TSTART = c1 136. T300 = 1. 137. TIME = 0. 138. NCYCLE = 0 139. * 90 NCYCLE = NCYCLE + 1 140. DO ncycle=1,itmax 141. C 142. C COMPUTE CAPITAL U, CAPITAL V, Z AND H 143. C 144. FSDX = 4./DX 145. FSDY = 4./DY 146. * T100 = SECOND(DUM) 147. call system_clock(count=c1, count_rate=r,count_max=max) 148. T100 = c1 149. DO J=1,N 150. DO I=1,M 151. CU(I+1,J) = .5*(P(I+1,J)+P(I,J))*U(I+1,J) 152. CV(I,J+1) = .5*(P(I,J+1)+P(I,J))*V(I,J+1) 153. Z(I+1,J+1) =(FSDX*(V(I+1,J+1)-V(I,J+1))-FSDY*(U(I+1,J+1) 154. 1 -U(I+1,J)))/(P(I,J)+P(I+1,J)+P(I+1,J+1)+P(I,J+1)) 155. H(I,J) = P(I,J)+.25*(U(I+1,J)*U(I+1,J)+U(I,J)*U(I,J) 156. 1 +V(I,J+1)*V(I,J+1)+V(I,J)*V(I,J)) 157. END DO 158. 100 END DO 159. * T100 = SECOND(DUM)-T100 160. call system_clock(count=c2,count_rate=r,count_max=max) 161. T100 = dble(c2-T100)/dble(r) 162. C 163. C PERIODIC CONTINUATION 164. C 165. DO J=1,N 166. CU(1,J) = CU(M+1,J) 167. CV(M+1,J+1) = CV(1,J+1) 168. Z(1,J+1) = Z(M+1,J+1) 169. H(M+1,J) = H(1,J) 170. 110 END DO 171. DO I=1,M 172. CU(I+1,N+1) = CU(I+1,1) 173. CV(I,1) = CV(I,N+1) 174. Z(I+1,1) = Z(I+1,N+1) 175. H(I,N+1) = H(I,1) 176. 115 END DO 177. CU(1,N+1) = CU(M+1,1) 178. CV(M+1,1) = CV(1,N+1) 179. Z(1,1) = Z(M+1,N+1) 180. H(M+1,N+1) = H(1,1) 181. C 182. C COMPUTE NEW VALUES U,V AND P 183. C 184. TDTS8 = TDT/8. 185. TDTSDX = TDT/DX 186. TDTSDY = TDT/DY 187. * T200 = SECOND(DUM) 188. call system_clock(count=c1, count_rate=r, count_max=max) 189. T200 = c1 190. DO J=1,N 191. DO I=1,M 192. UNEW(I+1,J) = UOLD(I+1,J)+ 193. 1 TDTS8*(Z(I+1,J+1)+Z(I+1,J))*(CV(I+1,J+1)+CV(I,J+1)+CV(I,J) 194. 2 +CV(I+1,J))-TDTSDX*(H(I+1,J)-H(I,J)) 195. VNEW(I,J+1) = VOLD(I,J+1)-TDTS8*(Z(I+1,J+1)+Z(I,J+1)) 196. 1 *(CU(I+1,J+1)+CU(I,J+1)+CU(I,J)+CU(I+1,J)) 197. 2 -TDTSDY*(H(I,J+1)-H(I,J)) 198. PNEW(I,J) = POLD(I,J)-TDTSDX*(CU(I+1,J)-CU(I,J)) 199. 1 -TDTSDY*(CV(I,J+1)-CV(I,J)) 200. END DO 201. 200 END DO 202. * T200 = SECOND(DUM)-T200 203. call system_clock(count=c2, count_rate=r, count_max=max) 204. T200 = dble(c2 -T200)/dble(r) 205. C 206. C PERIODIC CONTINUATION 207. C 208. DO J=1,N 209. UNEW(1,J) = UNEW(M+1,J) 210. VNEW(M+1,J+1) = VNEW(1,J+1) 211. PNEW(M+1,J) = PNEW(1,J) 212. 210 END DO 213. DO I=1,M 214. UNEW(I+1,N+1) = UNEW(I+1,1) 215. VNEW(I,1) = VNEW(I,N+1) 216. PNEW(I,N+1) = PNEW(I,1) 217. 215 END DO 218. UNEW(1,N+1) = UNEW(M+1,1) 219. VNEW(M+1,1) = VNEW(1,N+1) 220. PNEW(M+1,N+1) = PNEW(1,1) 221. * jr IF(NCYCLE .GT. ITMAX) CALL EXIT 222. 223. TIME = TIME + DT 224. * jr IF(MOD(NCYCLE,MPRINT) .NE. 0) GO TO 370 225. IF(MOD(NCYCLE,MPRINT) .EQ. 0) then 226. * jr 370 IF(NCYCLE .LE. 1) GO TO 310 227. 228. PTIME = TIME/3600. 229. WRITE(6,350) NCYCLE,PTIME 230. 350 FORMAT(//' CYCLE NUMBER',I5,' MODEL TIME IN HOURS', F6.2) 231. WRITE(6,355) (PNEW(I,I),I=1,MNMIN) 232. 355 FORMAT(/' DIAGONAL ELEMENTS OF P ' //(8E15.6)) 233. WRITE(6,360) (UNEW(I,I),I=1,MNMIN) 234. 360 FORMAT(/' DIAGONAL ELEMENTS OF U ' //(8E15.6)) 235. WRITE(6,365) (VNEW(I,I),I=1,MNMIN) 236. 365 FORMAT(/' DIAGONAL ELEMENTS OF V ' //(8E15.6)) 237. * jr added MFS310--don't know what actual mult factor should be 238. * jr changed divide by 1 million to divide by 100K since system_clock 239. * jr resolution is millisec rather than cpu_time's 10 millisec 240. MFS310 = 24.*M*N/T310/1.D5 241. * MFS100 = 24.*M*N/T100/1.E6 242. MFS100 = 24.*M*N/T100/1.D5 243. * MFS200 = 26.*M*N/T200/1.E6 244. MFS200 = 26.*M*N/T200/1.D5 245. * MFS300 = 15.*M*N/T300/1.E6 246. MFS300 = 15.*M*N/T300/1.D5 247. * CTIME = SECOND(DUM)-TSTART 248. call system_clock(count=c2, count_rate=r,count_max=max) 249. CTIME = dble(c2 - TSTART)/dble(r) 250. TCYC = CTIME/FLOAT(NCYCLE) 251. WRITE(6,375) NCYCLE,CTIME,TCYC,T310,MFS310,T200,MFS200,T300,MFS300 252. * 375 FORMAT(' CYCLE NUMBER',I5,' TOTAL COMPUTER TIME', E15.6, 253. * 1 ' TIME PER CYCLE', E15.6, / 254. * 2 ' TIME AND MEGAFLOPS FOR LOOP 310 ', E15.6,F6.1/ 255. * 3 ' TIME AND MEGAFLOPS FOR LOOP 200 ', E15.6,F6.1/ 256. * 4 ' TIME AND MEGAFLOPS FOR LOOP 300 ', E15.6,F6.1/ ) 257. 375 FORMAT(' CYCLE NUMBER',I5,' TOTAL COMPUTER TIME', D15.6, 258. 1 ' TIME PER CYCLE', D15.6, / 259. 2 ' TIME AND MEGAFLOPS FOR LOOP 310 ', D15.6,2x,D6.1/ 260. 3 ' TIME AND MEGAFLOPS FOR LOOP 200 ', D15.6,2x,D6.1/ 261. 4 ' TIME AND MEGAFLOPS FOR LOOP 300 ', D15.6,2x,D6.1/ ) 262. endif 263. * append calculated values of p, u, and v to netCDF file 264. istart(3) = ncycle + 1 265. * shape of record to be written (one ncycle at a time) 266. call my_ncwrite(ncid,p_id,istart,icount,p,m_len,n_len) 267. call my_ncwrite(ncid,u_id,istart,icount,u,m_len,n_len) 268. call my_ncwrite(ncid,v_id,istart,icount,v,m_len,n_len) 269. C 270. C TIME SMOOTHING AND UPDATE FOR NEXT CYCLE 271. C 272. * 370 IF(NCYCLE .LE. 1) GO TO 310 273. 370 IF(NCYCLE .GT. 1) then 274. * T300 = SECOND(DUM) 275. call system_clock(count=c1,count_rate=r,count_max=max) 276. T300 = c1 277. DO J=1,N 278. DO I=1,M 279. UOLD(I,J) = U(I,J)+ALPHA*(UNEW(I,J)-2.*U(I,J)+UOLD(I,J)) 280. VOLD(I,J) = V(I,J)+ALPHA*(VNEW(I,J)-2.*V(I,J)+VOLD(I,J)) 281. POLD(I,J) = P(I,J)+ALPHA*(PNEW(I,J)-2.*P(I,J)+POLD(I,J)) 282. U(I,J) = UNEW(I,J) 283. V(I,J) = VNEW(I,J) 284. P(I,J) = PNEW(I,J) 285. END DO 286. 300 END DO 287. * T300 = SECOND(DUM)-T300 288. call system_clock(count=c2,count_rate=r, count_max=max) 289. T300 = dble(c2 - T300)/dble(r) 290. C 291. C PERIODIC CONTINUATION 292. C 293. DO J=1,N 294. UOLD(M+1,J) = UOLD(1,J) 295. VOLD(M+1,J) = VOLD(1,J) 296. POLD(M+1,J) = POLD(1,J) 297. U(M+1,J) = U(1,J) 298. V(M+1,J) = V(1,J) 299. P(M+1,J) = P(1,J) 300. 320 END DO 301. DO I=1,M 302. UOLD(I,N+1) = UOLD(I,1) 303. VOLD(I,N+1) = VOLD(I,1) 304. POLD(I,N+1) = POLD(I,1) 305. U(I,N+1) = U(I,1) 306. V(I,N+1) = V(I,1) 307. P(I,N+1) = P(I,1) 308. 325 END DO 309. UOLD(M+1,N+1) = UOLD(1,1) 310. VOLD(M+1,N+1) = VOLD(1,1) 311. POLD(M+1,N+1) = POLD(1,1) 312. U(M+1,N+1) = U(1,1) 313. V(M+1,N+1) = V(1,1) 314. P(M+1,N+1) = P(1,1) 315. else 316. 310 TDT = TDT+TDT 317. call system_clock(count=c1, count_rate=r,count_max=max) 318. T310 = c1 319. DO J=1,NP1 320. DO I=1,MP1 321. UOLD(I,J) = U(I,J) 322. VOLD(I,J) = V(I,J) 323. POLD(I,J) = P(I,J) 324. U(I,J) = UNEW(I,J) 325. V(I,J) = VNEW(I,J) 326. P(I,J) = PNEW(I,J) 327. END DO 328. 400 END DO 329. call system_clock(count=c2, count_rate=r, count_max=max) 330. T310 = dble(c2 - T310)/dble(r) 331. endif 332. 333. * end big loop 90 334. End do 335. * close the netCDF file 336. iret = nf_close(ncid) 337. call check_err(iret) 338. END 339. 340. 341. subroutine check_err(iret) 342. integer iret 343. include 'netcdf.inc' 344. if(iret .ne. NF_NOERR) then 345. print *, nf_strerror(iret) 346. stop 347. endif 348. end 349. 350. subroutine netcdf_setup(m_len,n_len,ncid,p_id,u_id,v_id,istart,icount) 351. * Input args: m_len, n_len 352. * Output args: ncid, p_id,u_id,v_id,istart,icount) 353. integer m_len,n_len 354. * declarations for netCDF library 355. include 'netcdf.inc' 356. * error status return 357. integer iret 358. * netCDF id 359. integer ncid 360. * dimension ids 361. integer m_dim 362. integer n_dim 363. integer time_dim 364. * variable ids 365. integer p_id 366. integer u_id 367. integer v_id 368. * rank (number of dimensions) for each variable 369. integer p_rank, u_rank, v_rank 370. parameter (p_rank = 3) 371. parameter (u_rank = 3) 372. parameter (v_rank = 3) 373. * variable shapes 374. integer p_dims(p_rank) 375. integer u_dims(u_rank) 376. integer v_dims(v_rank) 377. integer istart(p_rank) 378. integer icount(p_rank) 379. 380. * enter define mode 381. iret = nf_create('shallowdat.nc', NF_CLOBBER, 382. + ncid) 383. call check_err(iret) 384. * define dimensions 385. iret = nf_def_dim(ncid, 'm', m_len, m_dim) 386. call check_err(iret) 387. iret = nf_def_dim(ncid, 'n', n_len, n_dim) 388. call check_err(iret) 389. * time is an unlimited dimension so that any number of 390. * records can be added 391. iret = nf_def_dim(ncid, 'time', NF_UNLIMITED, time_dim) 392. call check_err(iret) 393. * define variables 394. p_dims(1) = m_dim 395. p_dims(2) = n_dim 396. p_dims(3) = time_dim 397. iret = nf_def_var(ncid, 'p', NF_REAL, p_rank, p_dims, p_id) 398. call check_err(iret) 399. u_dims(1) = m_dim 400. u_dims(2) = n_dim 401. u_dims(3) = time_dim 402. iret = nf_def_var(ncid, 'u', NF_REAL, u_rank, u_dims, u_id) 403. call check_err(iret) 404. v_dims(1) = m_dim 405. v_dims(2) = n_dim 406. v_dims(3) = time_dim 407. iret = nf_def_var(ncid, 'v', NF_REAL, v_rank, v_dims, v_id) 408. call check_err(iret) 409. * start netCDF write at the first (1,1,1) position of the array 410. istart(1) = 1 411. istart(2) = 1 412. istart(3) = 1 413. * shape of record to be written (one ncycle at a time) 414. icount(1) = m_len 415. icount(2) = n_len 416. icount(3) = 1 417. 418. * leave define mode 419. iret = nf_enddef(ncid) 420. call check_err(iret) 421. 422. * end of netCDF definitions 423. end subroutine netcdf_setup 424. 425. subroutine my_ncwrite(id,varid,istart,icount,var,m_len,n_len) 426. * Input args: id, varid,istart,icount,var,m_len_n_len 427. * Write a whole array out to the netCDF file 428. integer id,varid,iret 429. integer icount(3) 430. integer istart(3) 431. real var (m_len,n_len) 432. iret = nf_put_vara_real(id,varid,istart,icount,var) 433. call check_err(iret) 434. 435. end subroutine my_ncwrite 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. |
|
|
|
|