source: palm/tags/release-5.0/UTIL/analyze_particle_netcdf_data.f90 @ 3979

Last change on this file since 3979 was 2696, checked in by kanani, 6 years ago

Merge of branch palm4u into trunk

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