source: palm/trunk/UTIL/combine_plot_fields.f90 @ 158

Last change on this file since 158 was 139, checked in by raasch, 17 years ago

New:
---

Plant canopy model of Watanabe (2004,BLM 112,307-341) added.
It can be switched on by the inipar parameter plant_canopy.
The inipar parameter canopy_mode can be used to prescribe a
plant canopy type. The default case is a homogeneous plant
canopy. Heterogeneous distributions of the leaf area
density and the canopy drag coefficient can be defined in the
new routine user_init_plant_canopy (user_interface).
The inipar parameters lad_surface, lad_vertical_gradient and
lad_vertical_gradient_level can be used in order to
prescribe the vertical profile of leaf area density. The
inipar parameter drag_coefficient determines the canopy
drag coefficient.
Finally, the inipar parameter pch_index determines the
index of the upper boundary of the plant canopy.

Allow new case bc_uv_t = 'dirichlet_0' for channel flow.

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

Pressure boundary conditions for vertical walls added to the multigrid solver.
They are applied using new wall flag arrays (wall_flags_..) which are defined
for each grid level. New argument gls added to routine user_init_grid
(user_interface).

Frequence of sorting particles can be controlled with new particles_par
parameter dt_sort_particles. Sorting is moved from the SGS timestep loop in
advec_particles after the end of this loop.

advec_particles, check_parameters, data_output_dvrp, header, init_3d_model, init_grid, init_particles, init_pegrid, modules, package_parin, parin, plant_canopy_model, read_var_list, read_3d_binary, user_interface, write_var_list, write_3d_binary

Changed:


Redefine initial nzb_local as the actual total size of topography (later the
extent of topography in nzb_local is reduced by 1dx at the E topography walls
and by 1dy at the N topography walls to form the basis for nzb_s_inner);
for consistency redefine 'single_building' case.

Vertical profiles now based on nzb_s_inner; they are divided by
ngp_2dh_s_inner (scalars, procucts of scalars) and ngp_2dh (staggered velocity
components and their products, procucts of scalars and velocity components),
respectively.

Allow two instead of one digit to specify isosurface and slicer variables.

Status of 3D-volume NetCDF data file only depends on switch netcdf_64bit_3d (check_open)

prognostic_equations include the respective wall_*flux in the parameter list of
calls of diffusion_s. Same as before, only the values of wall_heatflux(0:4)
can be assigned. At present, wall_humidityflux, wall_qflux, wall_salinityflux,
and wall_scalarflux are kept zero. diffusion_s uses the respective wall_*flux
instead of wall_heatflux. This update serves two purposes:

  • it avoids errors in calculations with humidity/scalar/salinity and prescribed

non-zero wall_heatflux,

  • it prepares PALM for a possible assignment of wall fluxes of

humidity/scalar/salinity in a future release.

buoyancy, check_open, data_output_dvrp, diffusion_s, diffusivities, flow_statistics, header, init_3d_model, init_dvrp, init_grid, modules, prognostic_equations

Errors:


Bugfix: summation of sums_l_l in diffusivities.

Several bugfixes in the ocean part: Initial density rho is calculated
(init_ocean). Error in initializing u_init and v_init removed
(check_parameters). Calculation of density flux now starts from
nzb+1 (production_e).

Bugfix: pleft/pright changed to pnorth/psouth in sendrecv of particle tail
numbers along y, small bugfixes in the SGS part (advec_particles)

Bugfix: model_string needed a default value (combine_plot_fields)

Bugfix: wavenumber calculation for even nx in routines maketri (poisfft)

Bugfix: assignment of fluxes at walls

Bugfix: absolute value of f must be used when calculating the Blackadar mixing length (init_1d_model)

advec_particles, check_parameters, combine_plot_fields, diffusion_s, diffusivities, init_ocean, init_1d_model, poisfft, production_e

  • Property svn:keywords set to Id
File size: 25.7 KB
Line 
1 PROGRAM combine_plot_fields
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: combine_plot_fields.f90 139 2007-11-29 09:37:41Z steinfeld $
11!
12! 114 2007-10-10 00:03:15Z raasch
13! Bugfix: model_string needed a default value
14!
15! Aug 07    Loop for processing of output by coupled runs, id_string does not
16!           contain modus any longer
17!
18! 18/01/06  Output of time-averaged data
19!
20! 25/05/05  Errors removed
21!
22! 26/04/05  Output in NetCDF format, iso2d and avs output only if parameter
23!           file exists
24!
25! 31/10/01  All comments and messages translated into English
26!
27! 23/02/99  Keine Bearbeitung komprimierter 3D-Daten
28! Ursprungsversion vom 28/07/97
29!
30!
31! Description:
32! ------------
33! This routine combines data of the PALM-subdomains into one file. In PALM
34! every processor element opens its own file and writes 2D- or 3D-binary-data
35! into it (different files are opened for xy-, xz-, yz-cross-sections and
36! 3D-data). For plotting or analyzing these PE-data have to be collected and
37! to be put into single files, which is done by this routine.
38! Output format is NetCDF. Additionally, a data are output in a binary format
39! readable by ISO2D-software (cross-sections) and by AVS (3D-data).
40!
41! This routine must be compiled with:
42! decalpha:
43!    f95 -cpp -D__netcdf -fast -r8 -I/usr/local/netcdf-3.5.1/include
44!    -L/usr/local/netcdf-3.5.1/lib -lnetcdf
45! IBM-Regatta:
46!    xlf95 -qsuffix=cpp=f90 -WF,-D__netcdf -qrealsize=8 -q64 -qmaxmem=-1 -Q
47!    -I /aws/dataformats/netcdf-3.6.0-p1/64-32/include-64
48!    -L/aws/dataformats/netcdf-3.6.0-p1/64-32/lib -lnetcdf -O3
49! IBM-Regatta KISTI:
50!    xlf95 -qsuffix=cpp=f90 -WF,-D__netcdf -qrealsize=8 -q64 -qmaxmem=-1 -Q
51!    -I /applic/netcdf64/src/f90
52!    -L/applic/lib/NETCDF64 -lnetcdf -O3
53! IBM-Regatta Yonsei (gfdl5):
54!    xlf95 -qsuffix=cpp=f90 -WF,-D__netcdf -qrealsize=8 -q64 -qmaxmem=-1 -Q
55!    -I /usr1/users/raasch/pub/netcdf-3.6.0-p1/include
56!    -L/usr1/users/raasch/pub/netcdf-3.6.0-p1/lib -lnetcdf -O3
57! IMUK:
58!    ifort combine...f90 -o combine...x
59!    -cpp -D__netcdf -I /muksoft/packages/netcdf/linux/include -axW -r8 -nbs
60!    -Vaxlib -L /muksoft/packages/netcdf/linux/lib -lnetcdf
61! NEC-SX6:
62!    sxf90 combine...f90 -o combine...x
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! Sun Fire X4600
66!    pgf95 combine...f90 -o combine...x
67!    -Mpreprocess -D__netcdf -I /home/usr5/mkanda/netcdf-3.6.1/src/f90 -r8
68!    -fast -L/home/usr5/mkanda/netcdf-3.6.1/src/libsrc -lnetcdf
69! FIMM:
70!    ifort combine...f90 -o combine...x
71!    -axW -cpp -openmp -r8 -nbs -convert little_endian -D__netcdf
72!    -I /local/netcdf/include -Vaxlib -L/local/netcdf/lib -lnetcdf
73!------------------------------------------------------------------------------!
74
75#if defined( __netcdf )
76    USE netcdf
77#endif
78
79    IMPLICIT NONE
80
81!
82!-- Local variables
83    CHARACTER (LEN=2)    ::  modus, model_string
84    CHARACTER (LEN=4)    ::  id_string
85    CHARACTER (LEN=10)   ::  dimname, var_name
86    CHARACTER (LEN=40)   ::  filename
87
88    CHARACTER (LEN=2000), DIMENSION(0:1) ::  var_list
89
90    INTEGER, PARAMETER ::  spk = SELECTED_REAL_KIND( 6 )
91
92    INTEGER ::  av, danz, i, id,             &
93                j, model, models, nc_stat, nxa, nxag, nxe, nxeg, nya,   &
94                nyag, nye, nyeg, nza, nzag, nze, nzeg, pos, time_step, xa, xe, &
95                ya, ye, za, ze
96
97    INTEGER, DIMENSION(0:1) ::  current_level, current_var, fanz, id_set, &
98         id_var_time, num_var
99
100    INTEGER, DIMENSION(4) ::  id_dims_loc
101
102    INTEGER, DIMENSION(0:1,4) ::  id_dims
103
104    INTEGER, DIMENSION(0:1,1000) ::  id_var, levels
105
106    LOGICAL ::  avs_output, compressed, found, iso2d_output, netcdf_output, &
107                netcdf_0, netcdf_1
108
109    REAL ::  dx, simulated_time
110    REAL, DIMENSION(:),   ALLOCATABLE   ::  eta, ho, hu
111    REAL, DIMENSION(:,:), ALLOCATABLE   ::  pf
112    REAL(spk), DIMENSION(:,:,:), ALLOCATABLE ::  pf3d
113
114    PRINT*, ''
115    PRINT*, ''
116    PRINT*, '*** combine_plot_fields ***'
117
118!
119!-- Find out if a coupled run has been carried out
120    INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
121    IF ( found )  THEN
122       models = 2
123       PRINT*, '    coupled run'
124    ELSE
125       models = 1
126       PRINT*, '    uncoupled run'
127    ENDIF
128
129!
130!-- Do everything for each model
131    DO model = 1, models
132!
133!--    Set the model string used to identify the filenames
134       model_string = ''
135       IF ( models == 2 )  THEN
136          PRINT*, ''
137          PRINT*, '*** combine_plot_fields ***'
138          IF ( model == 2 )  THEN
139             model_string = '_O'
140             PRINT*, '    now combining ocean data'
141             PRINT*, '    ========================'
142          ELSE
143             PRINT*, '    now combining atmosphere data'
144             PRINT*, '    ============================='
145          ENDIF
146       ENDIF
147!
148!--    2D-arrays for ISO2D
149!--    Main loop for the three different cross-sections, starting with
150!--    xy-section
151       modus = 'XY'
152       PRINT*, ''
153       DO  WHILE ( modus == 'XY'  .OR.  modus == 'XZ'  .OR.  modus == 'YZ' )
154!
155!--       Check, if file from PE0 exists. If it does not exist, PALM did not
156!--       create any output for this cross-section.
157          danz = 0
158          WRITE (id_string,'(I4.4)')  danz
159          INQUIRE ( &
160               FILE='PLOT2D_'//modus//TRIM( model_string )//'_'//id_string, &
161               EXIST=found )
162!
163!--       Find out the number of files (equal to the number of PEs which
164!--       have been used in PALM) and open them
165          DO  WHILE ( found )
166
167             OPEN ( danz+110, &
168                  FILE='PLOT2D_'//modus//TRIM( model_string )//'_'//id_string, &
169                  FORM='UNFORMATTED' )
170             danz = danz + 1
171             WRITE (id_string,'(I4.4)')  danz
172             INQUIRE ( &
173                  FILE='PLOT2D_'//modus//TRIM( model_string )//'_'//id_string, &
174                  EXIST=found )
175
176          ENDDO
177
178!
179!--       Inquire whether an iso2d parameter file exists
180          INQUIRE( FILE='PLOT2D_'//modus//'_GLOBAL'//TRIM( model_string ), &
181               EXIST=iso2d_output )
182
183!
184!--       Inquire whether a NetCDF file exists
185          INQUIRE( FILE='DATA_2D_'//modus//'_NETCDF'//TRIM( model_string ), &
186               EXIST=netcdf_0 )
187
188!
189!--       Inquire whether a NetCDF file for time-averaged data exists
190          INQUIRE( FILE='DATA_2D_'//modus//'_AV_NETCDF'//TRIM( model_string ),&
191               EXIST=netcdf_1 )
192
193          IF ( netcdf_0  .OR.  netcdf_1 )  THEN
194             netcdf_output = .TRUE.
195          ELSE
196             netcdf_output = .FALSE.
197          ENDIF
198
199!
200!--       Info-output
201          PRINT*, ''
202          PRINT*, '*** combine_plot_fields ***'
203#if defined( __netcdf )
204          IF ( netcdf_output )  PRINT*, '    NetCDF output enabled'
205#else
206          IF ( netcdf_output )  THEN
207             PRINT*, '--- Sorry, no NetCDF support on this host'
208             netcdf_output = .FALSE.
209          ENDIF
210#endif
211          IF ( danz /= 0 )  THEN
212             PRINT*, '    ',modus,'-section:  ', danz, ' file(s) found'
213          ELSE
214             PRINT*, '    no ', modus, '-section data available'
215          ENDIF
216
217          IF ( netcdf_output  .AND.  danz /= 0 )  THEN
218#if defined( __netcdf )
219             DO  av = 0, 1
220
221                IF ( av == 0  .AND.  .NOT.  netcdf_0 )  CYCLE
222                IF ( av == 1  .AND.  .NOT.  netcdf_1 )  CYCLE
223
224!
225!--             Open NetCDF dataset
226                IF ( av == 0 )  THEN
227                   filename = 'DATA_2D_'//modus//'_NETCDF' &
228                        //TRIM( model_string )
229                ELSE
230                   filename = 'DATA_2D_'//modus//'_AV_NETCDF' &
231                        //TRIM( model_string )
232                ENDIF
233                nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set(av) )
234                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 1 )
235
236!
237!--             Get the list of variables (order of variables corresponds with
238!--             the order of data on the binary file)
239                var_list(av) = ' '    ! GET_ATT does not assign trailing blanks
240                nc_stat = NF90_GET_ATT( id_set(av), NF90_GLOBAL, 'VAR_LIST', &
241                     var_list(av) )
242                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 2 )
243
244!
245!--             Inquire id of the time coordinate variable
246                nc_stat = NF90_INQ_VARID( id_set(av), 'time', id_var_time(av) )
247                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 3 )
248
249!
250!--             Count number of variables; there is one more semicolon in the
251!--             string than variable names
252                num_var(av) = -1
253                DO  i = 1, LEN( var_list(av) )
254                   IF ( var_list(av)(i:i) == ';' )  num_var(av) = num_var(av) +1
255                ENDDO
256
257!
258!--             Extract the variable names from the list and inquire their
259!--             NetCDF IDs
260                pos = INDEX( var_list(av), ';' )
261!
262!--             Loop over all variables
263                DO  i = 1, num_var(av)
264
265!
266!--                Extract variable name from list
267                   var_list(av) = var_list(av)(pos+1:)
268                   pos = INDEX( var_list(av), ';' )
269                   var_name = var_list(av)(1:pos-1)
270
271!
272!--                Get variable ID from name
273                   nc_stat = NF90_INQ_VARID( id_set(av), TRIM( var_name ), &
274                        id_var(av,i) )
275                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 4 )
276
277!
278!--                Get number of x/y/z levels for that variable
279                   nc_stat = NF90_INQUIRE_VARIABLE( id_set(av), id_var(av,i), &
280                        dimids = id_dims_loc )
281                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 5 )
282                   id_dims(av,:) = id_dims_loc
283
284!
285!--                Inquire dimension ID
286                   DO  j = 1, 4
287                      nc_stat = NF90_INQUIRE_DIMENSION( id_set(av), &
288                           id_dims(av,j), dimname, levels(av,i) )
289                      IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 6 )
290
291                      IF ( modus == 'XY' .AND. INDEX(dimname, 'z') /= 0 )  EXIT
292                      IF ( modus == 'XZ' .AND. INDEX(dimname, 'y') /= 0 )  EXIT
293                      IF ( modus == 'YZ' .AND. INDEX(dimname, 'x') /= 0 )  EXIT
294                   ENDDO
295
296                ENDDO
297
298             ENDDO   ! av = 0, 1
299
300          ENDIF
301#endif
302
303!
304!--       Read the arrays, as long as the end of the file is reached
305          fanz          =         0
306          current_level =         1
307          current_var   = 999999999
308
309          DO  WHILE ( danz /= 0 )
310
311!
312!--          Loop over all files (reading data of the subdomains)
313             DO  id = 0, danz-1
314!
315!--             File from PE0 contains special information at the beginning,
316!--             concerning the lower and upper indices of the total-domain used
317!--             in PALM (nxag, nxeg, nyag, nyeg) and the lower and upper indices
318!--             of the array to be writte by this routine (nxa, nxe, nya,
319!--             nye). Usually in the horizontal directions nxag=-1 and nxa=0
320!--             while all other variables have the same value (i.e. nxeg=nxe).
321!--             Allocate necessary arrays, open the output file and write
322!--             the coordinate informations needed by ISO2D.
323                IF ( id == 0  .AND.  fanz(0) == 0  .AND.  fanz(1) == 0 )  THEN
324                   READ ( id+110 )  nxag, nxeg, nyag, nyeg
325                   READ ( id+110 )  nxa, nxe, nya, nye
326                   ALLOCATE ( eta(nya:nye), ho(nxa:nxe), hu(nxa:nxe), &
327                        pf(nxag:nxeg,nyag:nyeg) )
328                   READ ( id+110 )  dx, eta, hu, ho
329
330                   IF ( iso2d_output )  THEN
331                      OPEN ( 2, FILE='PLOT2D_'//modus//TRIM( model_string ), &
332                           FORM='UNFORMATTED' )
333                      WRITE ( 2 )  dx, eta, hu, ho
334                   ENDIF
335                ENDIF
336!
337!--             Read output time
338                IF ( netcdf_output  .AND.  id == 0 )  THEN
339                   IF ( netcdf_1 )  THEN
340                      READ ( id+110, END=998 )  simulated_time, time_step, av
341                   ELSE
342!
343!--                   For compatibility with earlier PALM versions
344                      READ ( id+110, END=998 )  simulated_time, time_step
345                      av = 0
346                   ENDIF
347                ENDIF
348!
349!--             Read subdomain indices
350                READ ( id+110, END=998 )  xa, xe, ya, ye
351!
352!--             IF the PE made no output (in case that no part of the
353!--             cross-section is situated on this PE), indices have the
354!--             value -1
355                IF ( .NOT. ( xa == -1  .AND.  xe == -1  .AND. &
356                             ya == -1  .AND.  ye == -1 ) )  THEN
357!
358!--                Read the subdomain grid-point values
359                   READ ( id+110 )  pf(xa:xe,ya:ye)
360                ENDIF
361                IF ( id == 0 )  fanz(av) = fanz(av) + 1
362
363             ENDDO
364!
365!--          Write the data of the total domain cross-section
366             IF ( iso2d_output )  WRITE ( 2 )  pf(nxa:nxe,nya:nye)
367       
368!
369!--          Write same data in NetCDF format
370             IF ( netcdf_output )  THEN
371#if defined( __netcdf )
372!
373!--             Check if a new time step has begun; if yes write data to time
374!--             axis
375                IF ( current_var(av) > num_var(av) )  THEN
376                   current_var(av) = 1
377                   nc_stat = NF90_PUT_VAR( id_set(av), id_var_time(av), &
378                        (/ simulated_time /),        &
379                        start = (/ time_step /),     &
380                        count = (/ 1 /) )
381                   IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 7 )
382                ENDIF
383
384!
385!--             Now write the data; this is mode dependent
386                SELECT CASE ( modus )
387
388                   CASE ( 'XY' )
389                      nc_stat = NF90_PUT_VAR( id_set(av),                      &
390                                           id_var(av,current_var(av)),         &
391                                           pf(nxa:nxe,nya:nye),                &
392                             start = (/ 1, 1, current_level(av), time_step /), &
393                                      count = (/ nxe-nxa+1, nye-nya+1, 1, 1 /) )
394                      IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 8)
395                 
396                   CASE ( 'XZ' )
397                      nc_stat = NF90_PUT_VAR( id_set(av),                      &
398                                           id_var(av,current_var(av)),         &
399                                           pf(nxa:nxe,nya:nye),                &
400                             start = (/ 1, current_level(av), 1, time_step /), &
401                                      count = (/ nxe-nxa+1, 1, nye-nya+1, 1 /) )
402                      IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 9)
403
404                   CASE ( 'YZ' )
405                      nc_stat = NF90_PUT_VAR( id_set(av),                      &
406                                           id_var(av,current_var(av)),         &
407                                           pf(nxa:nxe,nya:nye),                &
408                             start = (/ current_level(av), 1, 1, time_step /), &
409                                      count = (/ 1, nxe-nxa+1, nye-nya+1, 1 /) )
410                      IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error(10)
411
412                END SELECT
413
414!
415!--             Data is written, check if max level is reached
416                current_level(av) = current_level(av) + 1
417                IF ( current_level(av) > levels(av,current_var(av)) )  THEN
418                   current_level(av) = 1
419                   current_var(av)   = current_var(av) + 1
420                ENDIF
421
422             ENDIF
423#endif
424
425          ENDDO
426
427998       IF ( danz /= 0 )  THEN
428!
429!--          Print the number of the arrays processed
430             WRITE (*,'(16X,I4,A)')  fanz(0)+fanz(1), ' array(s) processed'
431             IF ( fanz(1) /= 0 )  THEN
432                WRITE (*,'(16X,I4,A)')  fanz(1), ' array(s) are time-averaged'
433             ENDIF
434
435!
436!--          Close all files and deallocate arrays
437             DO  id = 0, danz-1
438                CLOSE ( id+110 )
439             ENDDO
440             CLOSE ( 2 )
441             DEALLOCATE ( eta, ho, hu, pf )
442          ENDIF
443
444!
445!--       Close the NetCDF file
446          IF ( netcdf_output  .AND.  danz /= 0 )  THEN
447#if defined( __netcdf )
448             IF ( netcdf_0 )  THEN
449                nc_stat = NF90_CLOSE( id_set(0) )
450                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 11 )
451             ENDIF
452             IF ( netcdf_1 )  THEN
453                nc_stat = NF90_CLOSE( id_set(1) )
454                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 12 )
455             ENDIF
456#endif
457          ENDIF
458
459!
460!--       Choose the next cross-section
461          SELECT CASE ( modus )
462             CASE ( 'XY' )
463                modus = 'XZ'
464             CASE ( 'XZ' )
465                modus = 'YZ'
466             CASE ( 'YZ' )
467                modus = 'no'
468          END SELECT
469
470       ENDDO
471
472
473!
474!--    Combine the 3D-arrays
475
476!
477!--    Inquire whether an avs fld file exists
478       INQUIRE( FILE='PLOT3D_FLD'//TRIM( model_string ), EXIST=avs_output )
479
480!
481!--    Inquire whether a NetCDF file exists
482       INQUIRE( FILE='DATA_3D_NETCDF'//TRIM( model_string ), EXIST=netcdf_0 )
483
484!
485!--    Inquire whether a NetCDF file for time-averaged data exists
486       INQUIRE( FILE='DATA_3D_AV_NETCDF'//TRIM( model_string ), EXIST=netcdf_1 )
487
488       IF ( netcdf_0  .OR.  netcdf_1 )  THEN
489          netcdf_output = .TRUE.
490       ELSE
491          netcdf_output = .FALSE.
492       ENDIF
493
494!
495!--    Check, if file from PE0 exists
496       danz = 0
497       WRITE (id_string,'(I4.4)')  danz
498       INQUIRE ( &
499            FILE='PLOT3D_DATA'//TRIM( model_string )//'_'//TRIM( id_string ),  &
500            EXIST=found )
501
502!
503!--    Combination only works, if data are not compressed. In that case,
504!--    PALM created a flag file (PLOT3D_COMPRESSED)
505       INQUIRE ( FILE='PLOT3D_COMPRESSED'//TRIM( model_string ), &
506            EXIST=compressed )
507
508!
509!--    Find out the number of files and open them
510       DO  WHILE ( found  .AND.  .NOT. compressed )
511
512          OPEN ( danz+110, &
513               FILE='PLOT3D_DATA'//TRIM( model_string )//'_'//TRIM(id_string), &
514               FORM='UNFORMATTED')
515          danz = danz + 1
516          WRITE (id_string,'(I4.4)')  danz
517          INQUIRE ( &
518               FILE='PLOT3D_DATA'//TRIM( model_string )//'_'//TRIM(id_string), &
519               EXIST=found )
520
521       ENDDO
522
523!
524!--    Info-output
525       PRINT*, ' '
526       PRINT*, '*** combine_plot_fields ***'
527#if defined( __netcdf )
528       IF ( netcdf_output )  PRINT*, '    NetCDF output enabled'
529#else
530       IF ( netcdf_output )  THEN
531          PRINT*, '--- Sorry, no NetCDF support on this host'
532          netcdf_output = .FALSE.
533       ENDIF
534#endif
535       IF ( danz /= 0 )  THEN
536          PRINT*, '    3D-data:     ', danz, ' file(s) found'
537       ELSE
538          IF ( found .AND. compressed )  THEN
539             PRINT*, '+++ no 3D-data processing, since data are compressed'
540          ELSE
541             PRINT*, '    no 3D-data file available'
542          ENDIF
543       ENDIF
544
545       IF ( netcdf_output  .AND.  danz /= 0 )  THEN
546#if defined( __netcdf )
547          DO  av = 0, 1
548
549             IF ( av == 0  .AND.  .NOT.  netcdf_0 )  CYCLE
550             IF ( av == 1  .AND.  .NOT.  netcdf_1 )  CYCLE
551
552!
553!--          Open NetCDF dataset
554             IF ( av == 0 )  THEN
555                filename = 'DATA_3D_NETCDF'//TRIM( model_string )
556             ELSE
557                filename = 'DATA_3D_AV_NETCDF'//TRIM( model_string )
558             ENDIF
559             nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set(av) )
560             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 13 )
561
562
563!
564!--          Get the list of variables (order of variables corresponds with the
565!--          order of data on the binary file)
566             var_list(av) = ' '    ! GET_ATT does not assign trailing blanks
567             nc_stat = NF90_GET_ATT( id_set(av), NF90_GLOBAL, 'VAR_LIST', &
568                  var_list(av) )
569             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 14 )
570
571!
572!--          Inquire id of the time coordinate variable
573             nc_stat = NF90_INQ_VARID( id_set(av), 'time', id_var_time(av) )
574             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 15 )
575
576!
577!--          Count number of variables; there is one more semicolon in the
578!--          string than variable names
579             num_var(av) = -1
580             DO  i = 1, LEN( var_list(av) )
581                IF ( var_list(av)(i:i) == ';' )  num_var(av) = num_var(av) + 1
582             ENDDO
583
584!
585!--          Extract the variable names from the list and inquire their NetCDF
586!--          IDs
587             pos = INDEX( var_list(av), ';' )
588!
589!--          Loop over all variables
590             DO  i = 1, num_var(av)
591
592!
593!--             Extract variable name from list
594                var_list(av) = var_list(av)(pos+1:)
595                pos = INDEX( var_list(av), ';' )
596                var_name = var_list(av)(1:pos-1)
597
598!
599!--             Get variable ID from name
600!                print*, '*** find id for "',TRIM( var_name ),'" begin'
601                nc_stat = NF90_INQ_VARID( id_set(av), TRIM( var_name ), &
602                     id_var(av,i) )
603                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 16 )
604!                print*, '*** find id for "',TRIM( var_name ),'" end'
605
606             ENDDO
607
608          ENDDO    ! av=0,1
609
610       ENDIF
611#endif
612
613!
614!--    Read arrays, until the end of the file is reached
615       current_var = 999999999
616       fanz = 0
617       DO  WHILE ( danz /= 0 )
618
619!
620!--       Loop over all files
621          DO  id = 0, danz-1
622!
623!--          File from PE0 contains special information at the beginning,
624!--          concerning the lower and upper indices of the total-domain used in
625!--          PALM (nxag, nxeg, nyag, nyeg, nzag, nzeg) and the lower and upper
626!--          indices of the array to be written by this routine (nxa, nxe, nya,
627!--          nye, nza, nze). Usually nxag=-1 and nxa=0, nyag=-1 and nya=0,
628!--          nzeg=nz and nze=nz_plot3d.
629!--          Allocate necessary array and open the output file.
630             IF ( id == 0  .AND.  fanz(0) == 0  .AND.  fanz(1) == 0 )  THEN
631                READ ( id+110 )  nxag, nxeg, nyag, nyeg, nzag, nzeg
632                READ ( id+110 )  nxa, nxe, nya, nye, nza, nze
633                ALLOCATE ( pf3d(nxag:nxeg,nyag:nyeg,nzag:nzeg) )
634                IF ( avs_output )  THEN
635                   OPEN ( 2, FILE='PLOT3D_DATA'//TRIM( model_string ), &
636                        FORM='UNFORMATTED' )
637                ENDIF
638             ENDIF
639
640!
641!--          Read output time
642             IF ( netcdf_output  .AND.  id == 0 )  THEN
643                IF ( netcdf_1 )  THEN
644                   READ ( id+110, END=999 )  simulated_time, time_step, av
645                ELSE
646!
647!--                For compatibility with earlier PALM versions
648                   READ ( id+110, END=999 )  simulated_time, time_step
649                   av = 0
650                ENDIF
651             ENDIF
652
653!
654!--          Read subdomain indices and grid point values
655             READ ( id+110, END=999 )  xa, xe, ya, ye, za, ze
656             READ ( id+110 )  pf3d(xa:xe,ya:ye,za:ze)
657             IF ( id == 0 )  fanz(av) = fanz(av) + 1
658
659          ENDDO
660
661!
662!--       Write data of the total domain
663          IF ( avs_output )  WRITE ( 2 )  pf3d(nxa:nxe,nya:nye,nza:nze)
664       
665!
666!--       Write same data in NetCDF format
667          IF ( netcdf_output )  THEN
668#if defined( __netcdf )
669!
670!--          Check if a new time step has begun; if yes write data to time axis
671             IF ( current_var(av) > num_var(av) )  THEN
672                current_var(av) = 1
673                nc_stat = NF90_PUT_VAR( id_set(av), id_var_time(av), &
674                                     (/ simulated_time /),&
675                                     start = (/ time_step /), count = (/ 1 /) )
676                IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 17 )
677             ENDIF
678
679!
680!--          Now write the data
681             nc_stat = NF90_PUT_VAR( id_set(av), id_var(av,current_var(av)), &
682                                  pf3d(nxa:nxe,nya:nye,nza:nze),      &
683                                  start = (/ 1, 1, 1, time_step /),   &
684                              count = (/ nxe-nxa+1, nye-nya+1, nze-nza+1, 1 /) )
685             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 18 )
686
687             current_var(av) = current_var(av) + 1
688
689#endif
690          ENDIF
691
692       ENDDO
693
694999    IF ( danz /= 0 )  THEN
695!
696!--       Print the number of arrays processed
697          WRITE (*,'(16X,I4,A)')  fanz(0)+fanz(1), ' array(s) processed'
698          IF ( fanz(1) /= 0 )  THEN
699             WRITE (*,'(16X,I4,A)')  fanz(1), ' array(s) are time-averaged'
700          ENDIF
701!
702!--       Close all files and deallocate array
703          DO  id = 0, danz-1
704             CLOSE ( id+110 )
705          ENDDO
706          CLOSE ( 2 )
707          DEALLOCATE ( pf3d )
708!
709!--       Close the NetCDF file
710          IF ( netcdf_output )  THEN
711#if defined( __netcdf )
712             IF ( netcdf_0 )  THEN
713                nc_stat = NF90_CLOSE( id_set(0) )
714                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 19 )
715             ENDIF
716             IF ( netcdf_1 )  THEN
717                nc_stat = NF90_CLOSE( id_set(1) )
718                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 20 )
719             ENDIF
720#endif
721          ENDIF
722       ENDIF
723
724    ENDDO  ! models
725
726
727 CONTAINS
728
729
730    SUBROUTINE handle_netcdf_error( errno )
731!
732!--    Prints out a text message corresponding to the current NetCDF status
733
734       IMPLICIT NONE
735
736       INTEGER, INTENT(IN) ::  errno
737
738       IF ( nc_stat /= NF90_NOERR )  THEN
739          PRINT*, '+++ combine_plot_fields  netcdf: ', av, errno, &
740                  TRIM( nf90_strerror( nc_stat ) )
741       ENDIF
742
743    END SUBROUTINE handle_netcdf_error
744
745
746 END PROGRAM combine_plot_fields
747
748
749
Note: See TracBrowser for help on using the repository browser.