HPC tutorial heading HPC tutorial heading
  •                          
 
NCAR

 

Example 4 - shallow4.f - Code with netCDF I/O added to facilitate data portability and visualization

Download 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

Line 14. include 'netcdf.inc' netCDF provides an interface for array-oriented data access and a library that provides an implementation of the interface. The netCDF library also defines a machine-independent format for representing scientific data. Together the interface, library, and format support the creation, access, and sharing of scientific data. Although netCDF files are binary, they are portable and self-documenting, making it easier for the user to examine and use the data. The package is freely available from: http://www.unidata.ucar.edu/software/netcdf

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.