HPC tutorial heading HPC tutorial heading
  •                          
 
NCAR

 

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.