source: palm/trunk/UTIL/analyze_particle_netcdf_data.f90 @ 589

Last change on this file since 589 was 1, checked in by raasch, 18 years ago

Initial repository layout and content

File size: 8.0 KB
RevLine 
[1]1 PROGRAM analyze_particle_netcdf_data
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! Initial revision 12/07/05
11!
12!
13! Description:
14! ------------
15! This is an EXAMPLE routine how to read NetCDF particle data output by PALM.
16! As an example, the mean heigth and radius of all particles are calculated
17! and output for each output time available on the file
18!
19!
20! Any additonal analyzation requested has to be added by the user by following
21! the steps given in this example!
22!
23!
24!
25! This routine must be compiled with:
26! decalpha:
27!    f95 -fast -r8 -I/usr/local/netcdf-3.5.1/include
28!    -L/usr/local/netcdf-3.5.1/lib -lnetcdf
29! IBM-Regatta:
30!    xlf95 -qsuffix=cpp=f90 -qrealsize=8 -q64 -qmaxmem=-1 -Q
31!    -I /aws/dataformats/netcdf-3.5.0/netcdf-64-32-3.5.0/include_F90_64
32!    -L/aws/dataformats/netcdf-3.5.0/netcdf-64-32-3.5.0/lib -lnetcdf -O3
33! IBM-Regatta KISTI:
34!    xlf95 -qsuffix=cpp=f90 -qrealsize=8 -q64 -qmaxmem=-1 -Q
35!    -I /applic/netcdf64/src/f90
36!    -L/applic/lib/NETCDF64 -lnetcdf -O3
37! IMUK:
38!    ifort analyze_particle_netcdf_data
39!    -I /muksoft/packages/netcdf/linux/include -axW -r8 -nbs -Vaxlib
40!    -L /muksoft/packages/netcdf/linux/lib -lnetcdf
41! NEC-SX6:
42!    sxf90 analyze_particle_netcdf_data.f90
43!    -I /pool/SX-6/netcdf/netcdf-3.6.0-p1/include  -C hopt -Wf '-A idbl4'
44!    -D__netcdf -L/pool/SX-6/netcdf/netcdf-3.6.0-p1/lib -lnetcdf
45!------------------------------------------------------------------------------!
46
47    USE netcdf
48
49    IMPLICIT NONE
50
51!
52!-- Local variables
53    CHARACTER (LEN=7)    ::  id_char
54    CHARACTER (LEN=2000) ::  string
55
56
57    INTEGER ::  f, fn, i, id_dim_time, id_var_rnop, id_var_r, id_var_time, &
58                id_var_x, id_var_y, id_var_z, n, nc_stat, ntl, start
59
60    INTEGER, DIMENSION(1000) ::  id_set
61    INTEGER, DIMENSION(1)    ::  id_dim_time_old
62
63    INTEGER, DIMENSION(:), ALLOCATABLE   ::  total_nop
64    INTEGER, DIMENSION(:,:), ALLOCATABLE ::  nop
65
66    LOGICAL ::  found
67
68    REAL ::  mean_r, mean_z
69    REAL, DIMENSION(:), ALLOCATABLE ::  prt_r, prt_x, prt_y, prt_z, tl
70
71!
72!-- Check, if file from PE0 exists. If it does not exist, PALM did not
73!-- create any output for this cross-section.
74    fn = 0
75    WRITE (id_char,'(''_'',I4.4)')  fn
76    INQUIRE ( FILE=id_char, EXIST=found )
77
78!
79!-- Find out the number of files (equal to the number of PEs which
80!-- have been used in PALM) and open them
81    IF ( .NOT. found )  THEN
82       PRINT*, '+++ no file _0000 found in current working directory'
83       STOP
84    ENDIF
85
86    DO  WHILE ( found )
87
88!
89!--    Open NetCDF dataset
90       nc_stat = NF90_OPEN( id_char, NF90_NOWRITE, id_set(fn) )
91       IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 1 )
92       fn = fn + 1
93       WRITE (id_char,'(''_'',I4.4)')  fn
94       INQUIRE ( FILE=id_char, EXIST=found )
95
96    ENDDO
97    fn = fn - 1
98
99    PRINT*, '*** ', fn+1, ' files found'
100
101!
102!-- Get the run description header and print it out
103    string = ' '
104    nc_stat = NF90_GET_ATT( id_set(0), NF90_GLOBAL, 'title', string )
105    IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 2 )
106    PRINT*, '*** run: ', TRIM( string )
107
108!
109!-- Get the available time levels
110    nc_stat = NF90_INQ_VARID( id_set(0), 'time', id_var_time )
111    IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 3 )
112
113    nc_stat = NF90_INQUIRE_VARIABLE( id_set(0), id_var_time, &
114                                     dimids = id_dim_time_old )
115    IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 4 )
116    id_dim_time = id_dim_time_old(1)
117
118    nc_stat = NF90_INQUIRE_DIMENSION( id_set(0), id_dim_time, len = ntl )
119    IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 5 )
120    ALLOCATE( tl(1:ntl) )
121    print*, 'ntl=',ntl
122
123    nc_stat = NF90_GET_VAR( id_set(0), id_var_time, tl )
124    IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 6 )
125
126    DO  n = 1, ntl
127       print*, '*** time_level(', n, ') =', tl(n)
128    ENDDO
129
130!
131!-- Get the number of particles used
132    nc_stat = NF90_INQ_VARID( id_set(0), 'real_num_of_prt', id_var_rnop )
133    IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 7 )
134
135    ALLOCATE( nop(1:ntl,0:fn), total_nop(1:ntl) )
136
137    DO  f = 0, fn
138
139       nc_stat = NF90_GET_VAR( id_set(f), id_var_rnop, nop(1:ntl,f) )
140       IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 8 )
141
142    ENDDO
143
144    total_nop = 0
145    DO  n = 1, ntl
146
147       DO  f = 0, fn
148          total_nop(n) = total_nop(n) + nop(n,f)
149       ENDDO
150
151       PRINT*, '*** time = ', tl(n), ' total # of particles: ', total_nop(n)
152
153    ENDDO
154
155!
156!-- Get the particle x and y coordinates
157    nc_stat = NF90_INQ_VARID( id_set(0), 'pt_x', id_var_x )
158    IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 9 )
159
160    nc_stat = NF90_INQ_VARID( id_set(0), 'pt_y', id_var_y )
161    IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 10 )
162
163    nc_stat = NF90_INQ_VARID( id_set(0), 'pt_z', id_var_z )
164    IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 11 )
165
166    nc_stat = NF90_INQ_VARID( id_set(0), 'pt_radius', id_var_r )
167    IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 12 )
168
169    PRINT*, ' '
170!
171!-- Loop over all timelevels
172    DO  n = 1, ntl
173
174       ALLOCATE( prt_x(total_nop(n)), prt_y(total_nop(n)), &
175                 prt_z(total_nop(n)), prt_r(total_nop(n)) )
176
177       start = 1
178
179!
180!--    Read the data from the files (one file per processor)
181       DO  f = 0, fn
182
183          nc_stat = NF90_GET_VAR( id_set(f), id_var_x,           &
184                                  prt_x(start:start+nop(n,f)-1), &
185                                  start = (/ 1, n /),            &
186                                  count = (/ nop(n,f), 1 /) )
187          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 13 )
188
189          nc_stat = NF90_GET_VAR( id_set(f), id_var_y,           &
190                                  prt_y(start:start+nop(n,f)-1), &
191                                  start = (/ 1, n /),            &
192                                  count = (/ nop(n,f), 1 /) )
193          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 14 )
194
195          nc_stat = NF90_GET_VAR( id_set(f), id_var_z,           &
196                                  prt_z(start:start+nop(n,f)-1), &
197                                  start = (/ 1, n /),            &
198                                  count = (/ nop(n,f), 1 /) )
199          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 15 )
200
201          nc_stat = NF90_GET_VAR( id_set(f), id_var_r,           &
202                                  prt_r(start:start+nop(n,f)-1), &
203                                  start = (/ 1, n /),            &
204                                  count = (/ nop(n,f), 1 /) )
205          IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( nc_stat, 16 )
206
207          start = start + nop(n,f)
208
209       ENDDO
210
211       mean_z = 0.0
212       mean_r = 0.0
213       DO  i = 1, total_nop(n)
214          mean_z = mean_z + prt_z(i)
215          mean_r = mean_r + prt_r(i)
216       ENDDO
217       mean_z = mean_z / total_nop(n)
218       mean_r = mean_r / total_nop(n)
219
220       PRINT*, '*** time = ', tl(n), ' mean height = ', mean_z, &
221                                     ' mean radius = ', mean_r
222
223!
224!--    prt_x, prt_y, prt_z, and prt_r contain the particle coordinates and
225!--    radii, respectively. Please output these data or carry out the
226!--    requested analyzing in this program before you deallocate the arrays.
227
228       DEALLOCATE( prt_x, prt_y, prt_z, prt_r )
229
230    ENDDO
231
232
233 END PROGRAM analyze_particle_netcdf_data
234
235
236
237 SUBROUTINE handle_netcdf_error( nc_stat, position )
238!
239!-- Prints out a text message corresponding to the current NetCDF status
240
241    USE netcdf
242
243    IMPLICIT NONE
244
245    INTEGER, INTENT(IN) ::  nc_stat, position
246
247    IF ( nc_stat /= NF90_NOERR )  THEN
248       PRINT*, '+++ analyze_particle_netcdf_data'
249       PRINT*, '    netcdf error: ', TRIM( nf90_strerror( nc_stat ) )
250       PRINT*, '    position = ', position
251       STOP
252    ENDIF
253
254 END SUBROUTINE handle_netcdf_error
Note: See TracBrowser for help on using the repository browser.