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

Last change on this file since 1736 was 1310, checked in by raasch, 11 years ago

update of GPL copyright

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