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

Last change on this file since 4578 was 4481, checked in by maronga, 5 years ago

Bugfix for copyright updates in document_changes; copyright update applied to all files

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