HPC tutorial heading HPC tutorial heading
  •                          
 
NCAR

 

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.