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

Last change on this file since 224 was 210, checked in by raasch, 16 years ago

updates in dvr routines for new dvr version

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