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

Last change on this file since 218 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
RevLine 
[1]1 PROGRAM combine_plot_fields
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
[210]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.
[191]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
[179]12! Bugfix: three misplaced #endif directives
[1]13!
14! Former revisions:
15! -----------------
[114]16! $Id: combine_plot_fields.f90 210 2008-11-06 08:54:02Z letzel $
17!
[139]18! 114 2007-10-10 00:03:15Z raasch
19! Bugfix: model_string needed a default value
20!
[114]21! Aug 07    Loop for processing of output by coupled runs, id_string does not
22!           contain modus any longer
23!
[1]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
[108]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
[1]79!------------------------------------------------------------------------------!
80
81#if defined( __netcdf )
82    USE netcdf
83#endif
84
85    IMPLICIT NONE
86
87!
88!-- Local variables
[108]89    CHARACTER (LEN=2)    ::  modus, model_string
90    CHARACTER (LEN=4)    ::  id_string
[1]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,             &
[210]99                j, k, model, models, nc_stat, nxa, nxag, nxe, nxeg, nya,   &
[1]100                nyag, nye, nyeg, nza, nzag, nze, nzeg, pos, time_step, xa, xe, &
[210]101                xxa, xxe, ya, ye, yya, yye, za, ze, zza, zze
[1]102
[108]103    INTEGER, DIMENSION(0:1) ::  current_level, current_var, fanz, id_set, &
104         id_var_time, num_var
[1]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
[191]117    REAL, DIMENSION(:,:), ALLOCATABLE   ::  pf, pf_tmp
118    REAL(spk), DIMENSION(:,:,:), ALLOCATABLE ::  pf3d, pf3d_tmp
[1]119
120    PRINT*, ''
[114]121    PRINT*, ''
[108]122    PRINT*, '*** combine_plot_fields ***'
[114]123
124!
125!-- Find out if a coupled run has been carried out
[108]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
[114]134
135!
136!-- Do everything for each model
[108]137    DO model = 1, models
[114]138!
139!--    Set the model string used to identify the filenames
140       model_string = ''
[108]141       IF ( models == 2 )  THEN
142          PRINT*, ''
143          PRINT*, '*** combine_plot_fields ***'
144          IF ( model == 2 )  THEN
[114]145             model_string = '_O'
[108]146             PRINT*, '    now combining ocean data'
147             PRINT*, '    ========================'
148          ELSE
149             PRINT*, '    now combining atmosphere data'
150             PRINT*, '    ============================='
151          ENDIF
152       ENDIF
[1]153!
[108]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' )
[1]160!
[108]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 )
[1]172
[108]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 )
[1]181
[108]182          ENDDO
[1]183
184!
[108]185!--       Inquire whether an iso2d parameter file exists
186          INQUIRE( FILE='PLOT2D_'//modus//'_GLOBAL'//TRIM( model_string ), &
187               EXIST=iso2d_output )
[1]188
189!
[108]190!--       Inquire whether a NetCDF file exists
191          INQUIRE( FILE='DATA_2D_'//modus//'_NETCDF'//TRIM( model_string ), &
192               EXIST=netcdf_0 )
[1]193
194!
[108]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 )
[1]198
[108]199          IF ( netcdf_0  .OR.  netcdf_1 )  THEN
200             netcdf_output = .TRUE.
201          ELSE
202             netcdf_output = .FALSE.
203          ENDIF
[1]204
205!
[108]206!--       Info-output
207          PRINT*, ''
208          PRINT*, '*** combine_plot_fields ***'
[1]209#if defined( __netcdf )
[108]210          IF ( netcdf_output )  PRINT*, '    NetCDF output enabled'
[1]211#else
[108]212          IF ( netcdf_output )  THEN
213             PRINT*, '--- Sorry, no NetCDF support on this host'
214             netcdf_output = .FALSE.
215          ENDIF
[1]216#endif
[108]217          IF ( danz /= 0 )  THEN
218             PRINT*, '    ',modus,'-section:  ', danz, ' file(s) found'
219          ELSE
220             PRINT*, '    no ', modus, '-section data available'
221          ENDIF
[1]222
[108]223          IF ( netcdf_output  .AND.  danz /= 0 )  THEN
[1]224#if defined( __netcdf )
[108]225             DO  av = 0, 1
[1]226
[108]227                IF ( av == 0  .AND.  .NOT.  netcdf_0 )  CYCLE
228                IF ( av == 1  .AND.  .NOT.  netcdf_1 )  CYCLE
[1]229
230!
[108]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 )
[1]241
242!
[108]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 )
[1]249
250!
[108]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 )
[1]254
255!
[108]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
[1]262
263!
[108]264!--             Extract the variable names from the list and inquire their
265!--             NetCDF IDs
266                pos = INDEX( var_list(av), ';' )
[1]267!
[108]268!--             Loop over all variables
269                DO  i = 1, num_var(av)
[1]270
271!
[108]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)
[1]276
277!
[108]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 )
[1]282
283!
[108]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
[1]289
290!
[108]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 )
[1]296
[108]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
[1]302                ENDDO
303
[108]304             ENDDO   ! av = 0, 1
[1]305
[179]306#endif
[108]307          ENDIF
[1]308
309!
[108]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
[1]314
[108]315          DO  WHILE ( danz /= 0 )
[1]316
317!
[108]318!--          Loop over all files (reading data of the subdomains)
319             DO  id = 0, danz-1
[1]320!
[108]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
[1]335
[108]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
[1]341                ENDIF
342!
[108]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
[1]348!
[108]349!--                   For compatibility with earlier PALM versions
350                      READ ( id+110, END=998 )  simulated_time, time_step
351                      av = 0
352                   ENDIF
[1]353                ENDIF
354!
[108]355!--             Read subdomain indices
356                READ ( id+110, END=998 )  xa, xe, ya, ye
[1]357!
[108]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
[1]363!
[108]364!--                Read the subdomain grid-point values
[191]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 )
[108]369                ENDIF
370                IF ( id == 0 )  fanz(av) = fanz(av) + 1
[1]371
[108]372             ENDDO
[1]373!
[108]374!--          Write the data of the total domain cross-section
375             IF ( iso2d_output )  WRITE ( 2 )  pf(nxa:nxe,nya:nye)
[1]376       
377!
[108]378!--          Write same data in NetCDF format
379             IF ( netcdf_output )  THEN
[1]380#if defined( __netcdf )
381!
[108]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
[1]392
393!
[108]394!--             Now write the data; this is mode dependent
395                SELECT CASE ( modus )
[1]396
[108]397                   CASE ( 'XY' )
398                      nc_stat = NF90_PUT_VAR( id_set(av),                      &
[1]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 /) )
[108]403                      IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 8)
[1]404                 
[108]405                   CASE ( 'XZ' )
406                      nc_stat = NF90_PUT_VAR( id_set(av),                      &
[1]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 /) )
[108]411                      IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 9)
[1]412
[108]413                   CASE ( 'YZ' )
414                      nc_stat = NF90_PUT_VAR( id_set(av),                      &
[1]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 /) )
[108]419                      IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error(10)
[1]420
[108]421                END SELECT
[1]422
423!
[108]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
[179]431#endif
[1]432             ENDIF
433
[108]434          ENDDO
[1]435
[108]436998       IF ( danz /= 0 )  THEN
[1]437!
[108]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
[1]443
444!
[108]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
[1]452
453!
[108]454!--       Close the NetCDF file
455          IF ( netcdf_output  .AND.  danz /= 0 )  THEN
[1]456#if defined( __netcdf )
[108]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
[1]466          ENDIF
467
468!
[108]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
[1]478
[108]479       ENDDO
[1]480
481
482!
[108]483!--    Combine the 3D-arrays
[1]484
485!
[108]486!--    Inquire whether an avs fld file exists
487       INQUIRE( FILE='PLOT3D_FLD'//TRIM( model_string ), EXIST=avs_output )
[1]488
489!
[108]490!--    Inquire whether a NetCDF file exists
491       INQUIRE( FILE='DATA_3D_NETCDF'//TRIM( model_string ), EXIST=netcdf_0 )
[1]492
493!
[108]494!--    Inquire whether a NetCDF file for time-averaged data exists
495       INQUIRE( FILE='DATA_3D_AV_NETCDF'//TRIM( model_string ), EXIST=netcdf_1 )
[1]496
[108]497       IF ( netcdf_0  .OR.  netcdf_1 )  THEN
498          netcdf_output = .TRUE.
499       ELSE
500          netcdf_output = .FALSE.
501       ENDIF
[1]502
503!
[108]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 )
[1]510
511!
[108]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 )
[1]516
517!
[108]518!--    Find out the number of files and open them
519       DO  WHILE ( found  .AND.  .NOT. compressed )
[1]520
[108]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 )
[1]529
[108]530       ENDDO
[1]531
532!
[108]533!--    Info-output
534       PRINT*, ' '
535       PRINT*, '*** combine_plot_fields ***'
[1]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
[108]544       IF ( danz /= 0 )  THEN
545          PRINT*, '    3D-data:     ', danz, ' file(s) found'
[1]546       ELSE
[108]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
[1]552       ENDIF
553
[108]554       IF ( netcdf_output  .AND.  danz /= 0 )  THEN
[1]555#if defined( __netcdf )
[108]556          DO  av = 0, 1
[1]557
[108]558             IF ( av == 0  .AND.  .NOT.  netcdf_0 )  CYCLE
559             IF ( av == 1  .AND.  .NOT.  netcdf_1 )  CYCLE
[1]560
561!
[108]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 )
[1]570
571
572!
[108]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 )
[1]579
580!
[108]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 )
[1]584
585!
[108]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
[1]592
593!
[108]594!--          Extract the variable names from the list and inquire their NetCDF
595!--          IDs
596             pos = INDEX( var_list(av), ';' )
[1]597!
[108]598!--          Loop over all variables
599             DO  i = 1, num_var(av)
[1]600
601!
[108]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)
[1]606
607!
[108]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 )
[1]612
[108]613             ENDDO
[1]614
[108]615          ENDDO    ! av=0,1
[1]616
[179]617#endif
[108]618       ENDIF
[1]619
620!
[108]621!--    Read arrays, until the end of the file is reached
622       current_var = 999999999
623       fanz = 0
624       DO  WHILE ( danz /= 0 )
[1]625
626!
[108]627!--       Loop over all files
628          DO  id = 0, danz-1
[1]629!
[108]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
[210]640                ALLOCATE ( pf3d(nxa:nxe,nya:nye,nza:nze) )
[108]641                IF ( avs_output )  THEN
642                   OPEN ( 2, FILE='PLOT3D_DATA'//TRIM( model_string ), &
643                        FORM='UNFORMATTED' )
644                ENDIF
[1]645             ENDIF
646
647!
[108]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
[1]653!
[108]654!--                For compatibility with earlier PALM versions
655                   READ ( id+110, END=999 )  simulated_time, time_step
656                   av = 0
657                ENDIF
[1]658             ENDIF
659
660!
[108]661!--          Read subdomain indices and grid point values
662             READ ( id+110, END=999 )  xa, xe, ya, ye, za, ze
[191]663             ALLOCATE( pf3d_tmp(xa:xe,ya:ye,za:ze) )
664             READ ( id+110 )  pf3d_tmp
[210]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
[191]680             DEALLOCATE( pf3d_tmp )
[108]681             IF ( id == 0 )  fanz(av) = fanz(av) + 1
[1]682
[108]683          ENDDO
[1]684
685!
[108]686!--       Write data of the total domain
687          IF ( avs_output )  WRITE ( 2 )  pf3d(nxa:nxe,nya:nye,nza:nze)
[1]688       
689!
[108]690!--       Write same data in NetCDF format
691          IF ( netcdf_output )  THEN
[1]692#if defined( __netcdf )
693!
[108]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), &
[1]698                                     (/ simulated_time /),&
699                                     start = (/ time_step /), count = (/ 1 /) )
[108]700                IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 17 )
701             ENDIF
[1]702
703!
[108]704!--          Now write the data
705             nc_stat = NF90_PUT_VAR( id_set(av), id_var(av,current_var(av)), &
[210]706                                     pf3d, start = (/ 1, 1, 1, time_step /), &
[1]707                              count = (/ nxe-nxa+1, nye-nya+1, nze-nza+1, 1 /) )
[108]708             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 18 )
[1]709
[108]710             current_var(av) = current_var(av) + 1
[1]711
712#endif
[108]713          ENDIF
[1]714
[108]715       ENDDO
[1]716
[108]717999    IF ( danz /= 0 )  THEN
[1]718!
[108]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
[1]724!
[108]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 )
[1]731!
[108]732!--       Close the NetCDF file
733          IF ( netcdf_output )  THEN
[1]734#if defined( __netcdf )
[108]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
[1]744          ENDIF
745       ENDIF
746
[108]747    ENDDO  ! models
[1]748
[108]749
[1]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
[179]761#if defined( __netcdf )
[1]762       IF ( nc_stat /= NF90_NOERR )  THEN
763          PRINT*, '+++ combine_plot_fields  netcdf: ', av, errno, &
764                  TRIM( nf90_strerror( nc_stat ) )
765       ENDIF
[179]766#endif
[1]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.