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

Last change on this file since 2809 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

  • Property svn:keywords set to Id
File size: 30.1 KB
RevLine 
[1]1 PROGRAM combine_plot_fields
2
3!------------------------------------------------------------------------------!
[2696]4! This file is part of the PALM model system.
[1046]5!
[2696]6! PALM is free software: you can redistribute it and/or modify it under the
7! terms of the GNU General Public License as published by the Free Software
8! Foundation, either version 3 of the License, or (at your option) any later
9! version.
[1046]10!
11! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
12! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
14!
15! You should have received a copy of the GNU General Public License along with
16! PALM. If not, see <http://www.gnu.org/licenses/>.
17!
[2718]18! Copyright 1997-2018  Leibniz Universitaet Hannover
[2696]19!------------------------------------------------------------------------------!
[1046]20!
[691]21! Current revisions:
[1]22! -----------------
[1469]23!
[2365]24!
[226]25! Former revisions:
26! -----------------
27! $Id: combine_plot_fields.f90 2718 2018-01-02 08:49:38Z schwenkel $
[2716]28! Corrected "Former revisions" section
29!
30! 2696 2017-12-14 17:12:51Z kanani
31! Change in file header (GPL part)
32!
33! 2669 2017-12-06 16:03:27Z raasch
[2669]34! data of 3d-nest runs are completely processed now
35!
36! 2523 2017-10-05 14:42:47Z kanani
[2523]37! Increased LEN for CHARACTER variable var_name, equal to the value in PALM
38!
39! 2512 2017-10-04 08:26:59Z raasch
[2512]40! PALM output does not contain ghost layer data any more
41! avs- and iso2d-related parts removed, handling of compressed data removed
42!
43! 2365 2017-08-21 14:59:59Z kanani
[2365]44! Vertical grid nesting implemented (SadiqHuq)
[226]45!
[2365]46! 1809 2016-04-05 20:13:28Z raasch
47!
[1809]48! 1808 2016-04-05 19:44:00Z raasch
49! cpu measurements are done with standard FORTRAN routine on every machine
50!
[1552]51! 1551 2015-03-03 14:18:16Z maronga
52! Adjustments for data output of soil model quantities
53!
[1469]54! 1468 2014-09-24 14:06:57Z maronga
55! Adapted for use on up to 6-digit processor cores (not tested)
56!
[1395]57! 1394 2014-05-06 10:17:31Z keck
58! KIND-attribute added to all INTEGER and REAL declaration statements
59!
[1047]60! 1046 2012-11-09 14:38:45Z maronga
61! code put under GPL (PALM 3.9)
62!
[692]63! 691 2011-03-04 08:45:30Z maronga
64! Bugfix: check for precursor ocean runs, removed typo
65!
[494]66! 493 2010-03-01 08:30:24Z raasch
67! Exit in case of already complete NetCDF data (due to parallel output in PALM)
68! cpu measurements included
69!
[226]70! 210 2008-11-06 08:54:02Z raasch
[210]71! Size of pf3d adjusted to the required output size (1 gridpoint less, along
72! all three dimensions), because output of a subset of the data
73! (pf3d(nxa:nxe...) in the NF90_PUT_VAR statement caused segmentation fault
74! with the INTEL compiler.
[191]75! Subdomain data are read into temporary arrays pf_tmp/pf3d_tmp in order to
76! avoid INTEL compiler warnings about (automatic) creation of temporary arrays
[179]77! Bugfix: three misplaced #endif directives
[1]78!
[139]79! 114 2007-10-10 00:03:15Z raasch
80! Bugfix: model_string needed a default value
81!
[114]82! Aug 07    Loop for processing of output by coupled runs, id_string does not
83!           contain modus any longer
84!
[1]85! 18/01/06  Output of time-averaged data
86!
87! 25/05/05  Errors removed
88!
89! 26/04/05  Output in NetCDF format, iso2d and avs output only if parameter
90!           file exists
91!
92! 31/10/01  All comments and messages translated into English
93!
94! 23/02/99  Keine Bearbeitung komprimierter 3D-Daten
95! Ursprungsversion vom 28/07/97
96!
97!
98! Description:
99! ------------
100! This routine combines data of the PALM-subdomains into one file. In PALM
101! every processor element opens its own file and writes 2D- or 3D-binary-data
102! into it (different files are opened for xy-, xz-, yz-cross-sections and
103! 3D-data). For plotting or analyzing these PE-data have to be collected and
104! to be put into single files, which is done by this routine.
105! Output format is NetCDF. Additionally, a data are output in a binary format
106! readable by ISO2D-software (cross-sections) and by AVS (3D-data).
107!------------------------------------------------------------------------------!
108
109#if defined( __netcdf )
110    USE netcdf
111#endif
112
113    IMPLICIT NONE
114
115!
116!-- Local variables
[2669]117    CHARACTER (LEN=2)    ::  modus
118    CHARACTER (LEN=4)    ::  model_string
[1468]119    CHARACTER (LEN=6)    ::  id_string
[2523]120    CHARACTER (LEN=30)   ::  dimname, var_name
[1]121    CHARACTER (LEN=40)   ::  filename
122
123    CHARACTER (LEN=2000), DIMENSION(0:1) ::  var_list
124
[1394]125    INTEGER, PARAMETER ::  iwp = 4, spk = SELECTED_REAL_KIND( 6 ), wp = 8
[1]126
[1394]127    INTEGER(iwp) ::  av, danz, i, id, j, k, model, models, nc_stat,            &
128                     nxa, nxag, nxe, nxeg, nya, nyag, nye, nyeg,               &
[2669]129                     nza, nze, pos, time_step, xa, xe, xxa, xxe, ya, ya_do,    &
130                     ya_tot, ye, ye_do, ye_tot, yya, yye, za, ze, zza, zze
[1]131
[1808]132    INTEGER(8)                        ::  count, count_rate
[493]133
[1394]134    INTEGER(iwp), DIMENSION(0:1)      ::  current_level, current_var,          &
135                                          fanz, id_set, id_var_time, num_var
136    INTEGER(iwp), DIMENSION(4)        ::  id_dims_loc
137    INTEGER(iwp), DIMENSION(0:1,4)    ::  id_dims
138    INTEGER(iwp), DIMENSION(0:1,1000) ::  id_var, levels
[1]139
[2669]140    LOGICAL  ::  found, nest3d, netcdf_output, netcdf_parallel, netcdf_0,      &
141                 netcdf_1, vnest
[1]142
[1394]143    REAL(wp) ::  cpu_start_time, cpu_end_time, dx, simulated_time
144    REAL(wp),  DIMENSION(:,:), ALLOCATABLE   ::  pf, pf_tmp
[191]145    REAL(spk), DIMENSION(:,:,:), ALLOCATABLE ::  pf3d, pf3d_tmp
[1]146
[2365]147
148
[1]149    PRINT*, ''
[114]150    PRINT*, ''
[108]151    PRINT*, '*** combine_plot_fields ***'
[114]152
153!
[2669]154!-- Find out if a coupled or nested run has been carried out
[108]155    INQUIRE( FILE='COUPLING_PORT_OPENED', EXIST=found )
[2365]156    INQUIRE( FILE='VNESTING_PORT_OPENED', EXIST=vnest )
[2669]157    INQUIRE( FILE='3DNESTING', EXIST=nest3d )
[108]158    IF ( found )  THEN
159       models = 2
160       PRINT*, '    coupled run'
[2365]161    ELSEIF ( vnest )  THEN
162       models = 2
163       PRINT*, '    Vertically nested grid coupling'
[2669]164    ELSEIF ( nest3d )  THEN
165       OPEN( 90, FILE='3DNESTING', FORM='FORMATTED' )
166       READ ( 90, '(I2)' )  models
167       CLOSE ( 90 )
168       PRINT*, '    3d-nest run'
169       PRINT*, '    number of nest domains = ', models
[108]170    ELSE
171       models = 1
172       PRINT*, '    uncoupled run'
173    ENDIF
[114]174
175!
[691]176!-- Find out if a precursor ocean run has been carried out
177    INQUIRE( FILE='PRECURSOR_OCEAN', EXIST=found )
178    IF ( found )  THEN
179       model_string = '_O'
180       PRINT*, '    precursor ocean run'
181    ELSE
182       model_string = ''
183    ENDIF
184
185!
[114]186!-- Do everything for each model
[108]187    DO model = 1, models
[114]188!
189!--    Set the model string used to identify the filenames
[2669]190       IF ( found  .OR.  vnest )  THEN
[108]191          PRINT*, ''
192          IF ( model == 2 )  THEN
[2365]193             IF ( vnest )  THEN
194                model_string = '_N'
195                PRINT*, '    now combining FINE data'
196                PRINT*, '    ========================'
197             ELSE
198                model_string = '_O'
199                PRINT*, '    now combining ocean data'
200                PRINT*, '    ========================'
201             ENDIF
[108]202          ELSE
[2365]203             IF ( vnest )  THEN
204                PRINT*, '    now combining COARSE data'
205                PRINT*, '    ============================='
206             ELSE
207                PRINT*, '    now combining atmosphere data'
208                PRINT*, '    ============================='
209             ENDIF
[108]210          ENDIF
[2669]211       ELSEIF ( nest3d )  THEN
212          PRINT*, ''
213          PRINT*, '--> processing nest id = ', model
214          IF ( model == 1 )  THEN
215             model_string = ''
216          ELSE
217             WRITE( model_string, '(A2,I2.2)' )  '_N', model
218          ENDIF
[108]219       ENDIF
[1]220!
[108]221!--    2D-arrays for ISO2D
222!--    Main loop for the three different cross-sections, starting with
223!--    xy-section
224       modus = 'XY'
225       PRINT*, ''
226       DO  WHILE ( modus == 'XY'  .OR.  modus == 'XZ'  .OR.  modus == 'YZ' )
[1]227!
[493]228!--       Take current time
229          CALL SYSTEM_CLOCK( count, count_rate )
230          cpu_start_time = REAL( count ) / REAL( count_rate )
231
232          netcdf_parallel = .FALSE.
233!
[108]234!--       Check, if file from PE0 exists. If it does not exist, PALM did not
235!--       create any output for this cross-section.
236          danz = 0
[1468]237          WRITE (id_string,'(I6.6)')  danz
[108]238          INQUIRE ( &
239               FILE='PLOT2D_'//modus//TRIM( model_string )//'_'//id_string, &
240               EXIST=found )
241!
242!--       Find out the number of files (equal to the number of PEs which
243!--       have been used in PALM) and open them
244          DO  WHILE ( found )
[1]245
[108]246             OPEN ( danz+110, &
247                  FILE='PLOT2D_'//modus//TRIM( model_string )//'_'//id_string, &
248                  FORM='UNFORMATTED' )
249             danz = danz + 1
[1468]250             WRITE (id_string,'(I6.6)')  danz
[108]251             INQUIRE ( &
252                  FILE='PLOT2D_'//modus//TRIM( model_string )//'_'//id_string, &
253                  EXIST=found )
[1]254
[108]255          ENDDO
[1]256
257!
[108]258!--       Inquire whether a NetCDF file exists
259          INQUIRE( FILE='DATA_2D_'//modus//'_NETCDF'//TRIM( model_string ), &
260               EXIST=netcdf_0 )
[1]261
262!
[108]263!--       Inquire whether a NetCDF file for time-averaged data exists
264          INQUIRE( FILE='DATA_2D_'//modus//'_AV_NETCDF'//TRIM( model_string ),&
265               EXIST=netcdf_1 )
[1]266
[108]267          IF ( netcdf_0  .OR.  netcdf_1 )  THEN
268             netcdf_output = .TRUE.
[493]269!
270!--          Inquire whether the NetCDF file is already complete (parallel
271!--          output)
272             INQUIRE( FILE='NO_COMBINE_PLOT_FIELDS_'//modus, &
273                      EXIST=netcdf_parallel )
274             IF ( netcdf_parallel )  THEN
275                netcdf_parallel = .TRUE.
276             ELSE
277                netcdf_parallel = .FALSE.
278             ENDIF
[108]279          ELSE
280             netcdf_output = .FALSE.
281          ENDIF
[1]282
283!
[108]284!--       Info-output
285          PRINT*, ''
[1]286#if defined( __netcdf )
[493]287          IF ( netcdf_output )  THEN
288             IF ( netcdf_parallel )  THEN
289             PRINT*, '    NetCDF ' // modus // '-data are in one file ', &
290                          '(NetCDF4-format) - merging not neccessary'
291             ELSE
292                PRINT*, '    NetCDF output enabled'
293             ENDIF
294          ENDIF
[1]295#else
[108]296          IF ( netcdf_output )  THEN
297             PRINT*, '--- Sorry, no NetCDF support on this host'
298             netcdf_output = .FALSE.
299          ENDIF
[1]300#endif
[493]301          IF ( .NOT. netcdf_parallel )  THEN
302             IF ( danz /= 0 )  THEN
303                PRINT*, '    ',modus,'-section:  ', danz, ' file(s) found'
304             ELSE
305                PRINT*, '    no ', modus, '-section data available'
306             ENDIF
[108]307          ENDIF
[1]308
[493]309          IF ( netcdf_output  .AND. .NOT. netcdf_parallel  .AND.  danz /= 0 ) &
310          THEN
[1]311#if defined( __netcdf )
[108]312             DO  av = 0, 1
[1]313
[108]314                IF ( av == 0  .AND.  .NOT.  netcdf_0 )  CYCLE
315                IF ( av == 1  .AND.  .NOT.  netcdf_1 )  CYCLE
[1]316
317!
[108]318!--             Open NetCDF dataset
319                IF ( av == 0 )  THEN
320                   filename = 'DATA_2D_'//modus//'_NETCDF' &
321                        //TRIM( model_string )
322                ELSE
323                   filename = 'DATA_2D_'//modus//'_AV_NETCDF' &
324                        //TRIM( model_string )
325                ENDIF
326                nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set(av) )
327                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 1 )
[1]328
329!
[108]330!--             Get the list of variables (order of variables corresponds with
331!--             the order of data on the binary file)
332                var_list(av) = ' '    ! GET_ATT does not assign trailing blanks
333                nc_stat = NF90_GET_ATT( id_set(av), NF90_GLOBAL, 'VAR_LIST', &
334                     var_list(av) )
335                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 2 )
[1]336
337!
[108]338!--             Inquire id of the time coordinate variable
339                nc_stat = NF90_INQ_VARID( id_set(av), 'time', id_var_time(av) )
340                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 3 )
[1]341
342!
[108]343!--             Count number of variables; there is one more semicolon in the
344!--             string than variable names
345                num_var(av) = -1
346                DO  i = 1, LEN( var_list(av) )
347                   IF ( var_list(av)(i:i) == ';' )  num_var(av) = num_var(av) +1
348                ENDDO
[1]349
350!
[108]351!--             Extract the variable names from the list and inquire their
352!--             NetCDF IDs
353                pos = INDEX( var_list(av), ';' )
[1]354!
[108]355!--             Loop over all variables
356                DO  i = 1, num_var(av)
[1]357
358!
[108]359!--                Extract variable name from list
360                   var_list(av) = var_list(av)(pos+1:)
361                   pos = INDEX( var_list(av), ';' )
362                   var_name = var_list(av)(1:pos-1)
[1]363
364!
[108]365!--                Get variable ID from name
366                   nc_stat = NF90_INQ_VARID( id_set(av), TRIM( var_name ), &
367                        id_var(av,i) )
368                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 4 )
[1]369
370!
[108]371!--                Get number of x/y/z levels for that variable
372                   nc_stat = NF90_INQUIRE_VARIABLE( id_set(av), id_var(av,i), &
373                        dimids = id_dims_loc )
374                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 5 )
375                   id_dims(av,:) = id_dims_loc
[1]376
377!
[108]378!--                Inquire dimension ID
379                   DO  j = 1, 4
380                      nc_stat = NF90_INQUIRE_DIMENSION( id_set(av), &
381                           id_dims(av,j), dimname, levels(av,i) )
382                      IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 6 )
[1]383
[108]384                      IF ( modus == 'XY' .AND. INDEX(dimname, 'z') /= 0 )  EXIT
385                      IF ( modus == 'XZ' .AND. INDEX(dimname, 'y') /= 0 )  EXIT
386                      IF ( modus == 'YZ' .AND. INDEX(dimname, 'x') /= 0 )  EXIT
387                   ENDDO
388
[1]389                ENDDO
390
[108]391             ENDDO   ! av = 0, 1
[1]392
[179]393#endif
[108]394          ENDIF
[1]395
396!
[108]397!--       Read the arrays, as long as the end of the file is reached
[493]398          IF ( .NOT. netcdf_parallel )  THEN
[1]399
[493]400             fanz          =         0
401             current_level =         1
402             current_var   = 999999999
[1]403
[493]404             DO  WHILE ( danz /= 0 )
405
[1]406!
[493]407!--             Loop over all files (reading data of the subdomains)
408                DO  id = 0, danz-1
[1]409!
[493]410!--                File from PE0 contains special information at the beginning,
411!--                concerning the lower and upper indices of the total-domain
[2512]412!--                used in PALM (nxa, nxe, nya, nye).
[493]413!--                Allocate necessary arrays, open the output file and write
414!--                the coordinate informations needed by ISO2D.
415                   IF ( id == 0  .AND.  fanz(0) == 0  .AND.  fanz(1) == 0 ) THEN
[2512]416
[493]417                      READ ( id+110 )  nxa, nxe, nya, nye
[2512]418                      ALLOCATE ( pf(nxa:nxe,nya:nye) )
[1551]419!
420!--                   Set actual domain bounds to total domain
421                      ya_do = nya
422                      ye_do = nye
423
[108]424                   ENDIF
[1]425!
[493]426!--                Read output time
427                   IF ( netcdf_output  .AND.  id == 0 )  THEN
[2512]428                      READ ( id+110, END=998 )  simulated_time, time_step, av
[108]429                   ENDIF
[1]430!
[493]431!--                Read subdomain indices
[1551]432                   READ ( id+110, END=998 )  xa, xe, ya, ye, ya_tot, ye_tot
433
[1]434!
[493]435!--                IF the PE made no output (in case that no part of the
436!--                cross-section is situated on this PE), indices have the
437!--                value -1
438                   IF ( .NOT. ( xa == -1  .AND.  xe == -1  .AND. &
439                                ya == -1  .AND.  ye == -1 ) )  THEN
[1551]440
441
[1]442!
[493]443!--                   Read the subdomain grid-point values
444                      ALLOCATE( pf_tmp(xa:xe,ya:ye) )
445                      READ ( id+110 )  pf_tmp
[1551]446!
447!--                   Calculate indices on atmospheric grid (e.g. for soil model
448!--                   quantities)
449                      IF ( ya /= ya_tot .OR. ye /= ye_tot )  THEN
450                         ye_do = ye - ya
451                         ya_do = ya
452                         pf(xa:xe,ya_do:ye_do) = pf_tmp
453                      ELSE
454                         ye_do = nye
455                         ya_do = nya
456                         pf(xa:xe,ya:ye) = pf_tmp
457                      ENDIF
458
[493]459                      DEALLOCATE( pf_tmp )
460                   ENDIF
461                   IF ( id == 0 )  fanz(av) = fanz(av) + 1
[1]462
[493]463                ENDDO
[1]464!
[2512]465!--             Write data in NetCDF format
[493]466                IF ( netcdf_output )  THEN
[1]467#if defined( __netcdf )
468!
[493]469!--                Check if a new time step has begun; if yes write data to
470!--                time axis
471                   IF ( current_var(av) > num_var(av) )  THEN
472                      current_var(av) = 1
473                      nc_stat = NF90_PUT_VAR( id_set(av), id_var_time(av), &
474                                              (/ simulated_time /),        &
475                                              start = (/ time_step /),     &
476                                              count = (/ 1 /) )
477                      IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 7 )
478                   ENDIF
[1]479
480!
[493]481!--                Now write the data; this is mode dependent
482                   SELECT CASE ( modus )
[1]483
[493]484                      CASE ( 'XY' )
485                         nc_stat = NF90_PUT_VAR( id_set(av),                   &
[1]486                                           id_var(av,current_var(av)),         &
487                                           pf(nxa:nxe,nya:nye),                &
488                             start = (/ 1, 1, current_level(av), time_step /), &
489                                      count = (/ nxe-nxa+1, nye-nya+1, 1, 1 /) )
[493]490                         IF ( nc_stat /= NF90_NOERR )  THEN
491                            CALL handle_netcdf_error( 8 )
492                         ENDIF
[1]493                 
[493]494                      CASE ( 'XZ' )
495                         nc_stat = NF90_PUT_VAR( id_set(av),                   &
[1]496                                           id_var(av,current_var(av)),         &
[1551]497                                           pf(nxa:nxe,ya_do:ye_do),            &
[1]498                             start = (/ 1, current_level(av), 1, time_step /), &
[1551]499                                      count = (/ nxe-nxa+1, 1, ye_do-ya_do+1, 1 /) )
[493]500                         IF ( nc_stat /= NF90_NOERR )  THEN
501                            CALL handle_netcdf_error( 9 )
502                         ENDIF
[1]503
[493]504                      CASE ( 'YZ' )
505                         nc_stat = NF90_PUT_VAR( id_set(av),                   &
[1]506                                           id_var(av,current_var(av)),         &
[1551]507                                           pf(nxa:nxe,ya_do:ye_do),            &
[1]508                             start = (/ current_level(av), 1, 1, time_step /), &
[1551]509                                      count = (/ 1, nxe-nxa+1, ye_do-ya_do+1, 1 /) )
[493]510                         IF ( nc_stat /= NF90_NOERR )  THEN
511                            CALL handle_netcdf_error( 10 )
512                         ENDIF
[1]513
[493]514                   END SELECT
[1]515
516!
[493]517!--                Data is written, check if max level is reached
518                   current_level(av) = current_level(av) + 1
519                   IF ( current_level(av) > levels(av,current_var(av)) )  THEN
520                      current_level(av) = 1
521                      current_var(av)   = current_var(av) + 1
522                   ENDIF
[108]523
[179]524#endif
[493]525                ENDIF
[1]526
[493]527             ENDDO
[1]528
[493]529          ENDIF
530
531998       IF ( danz /= 0  .AND.  .NOT. netcdf_parallel )  THEN
[1]532!
[108]533!--          Print the number of the arrays processed
534             WRITE (*,'(16X,I4,A)')  fanz(0)+fanz(1), ' array(s) processed'
535             IF ( fanz(1) /= 0 )  THEN
536                WRITE (*,'(16X,I4,A)')  fanz(1), ' array(s) are time-averaged'
537             ENDIF
[1]538
539!
[108]540!--          Close all files and deallocate arrays
541             DO  id = 0, danz-1
542                CLOSE ( id+110 )
543             ENDDO
544             CLOSE ( 2 )
[2512]545             DEALLOCATE ( pf )
[1]546
547!
[493]548!--          Close the NetCDF file
549             IF ( netcdf_output )  THEN
[1]550#if defined( __netcdf )
[493]551                IF ( netcdf_0 )  THEN
552                   nc_stat = NF90_CLOSE( id_set(0) )
553                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 11 )
554                ENDIF
555                IF ( netcdf_1 )  THEN
556                   nc_stat = NF90_CLOSE( id_set(1) )
557                   IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 12 )
558                ENDIF
559#endif
[108]560             ENDIF
[493]561          ENDIF
562
563!
564!--       Output required cpu time
565          IF ( danz /= 0  .AND.  .NOT. netcdf_parallel )  THEN
566             CALL SYSTEM_CLOCK( count, count_rate )
567             cpu_end_time = REAL( count ) / REAL( count_rate )
568             WRITE (*,'(5X,A,F9.3,A)')  'Required cpu-time: ', &
569                                        cpu_end_time-cpu_start_time, ' sec'
[1]570          ENDIF
571
572!
[108]573!--       Choose the next cross-section
574          SELECT CASE ( modus )
575             CASE ( 'XY' )
576                modus = 'XZ'
577             CASE ( 'XZ' )
578                modus = 'YZ'
579             CASE ( 'YZ' )
580                modus = 'no'
581          END SELECT
[1]582
[108]583       ENDDO
[1]584
585
586!
[108]587!--    Combine the 3D-arrays
[493]588       netcdf_parallel = .FALSE.
[1]589
590!
[493]591!--    Info-output
592       PRINT*, ' '
593
594!
595!--    Take current time
596       CALL SYSTEM_CLOCK( count, count_rate )
597       cpu_start_time = REAL( count ) / REAL( count_rate )
598
599!
[108]600!--    Inquire whether a NetCDF file exists
601       INQUIRE( FILE='DATA_3D_NETCDF'//TRIM( model_string ), EXIST=netcdf_0 )
[1]602
603!
[108]604!--    Inquire whether a NetCDF file for time-averaged data exists
605       INQUIRE( FILE='DATA_3D_AV_NETCDF'//TRIM( model_string ), EXIST=netcdf_1 )
[1]606
[108]607       IF ( netcdf_0  .OR.  netcdf_1 )  THEN
608          netcdf_output = .TRUE.
[493]609!
610!--       Inquire whether the NetCDF file is already complete (parallel output)
611          INQUIRE( FILE='NO_COMBINE_PLOT_FIELDS_3D', EXIST=netcdf_parallel )
612          IF ( netcdf_parallel )  THEN
613             netcdf_parallel = .TRUE.
614          ELSE
615             netcdf_parallel = .FALSE.
616          ENDIF
[108]617       ELSE
618          netcdf_output = .FALSE.
619       ENDIF
[1]620
621!
[493]622!--    Check, if file from PE0 exists; not neccessary in case of parallel
623!--    PALM output
624       IF ( .NOT. netcdf_parallel )  THEN
625          danz = 0
[1468]626          WRITE (id_string,'(I6.6)')  danz
[493]627          INQUIRE ( &
628               FILE='PLOT3D_DATA'//TRIM( model_string )//'_'//TRIM( id_string ),  &
629               EXIST=found )
630       ELSE
631          found = .FALSE.
632       ENDIF
[1]633
634!
[108]635!--    Find out the number of files and open them
[2512]636       DO  WHILE ( found )
[1]637
[108]638          OPEN ( danz+110, &
639               FILE='PLOT3D_DATA'//TRIM( model_string )//'_'//TRIM(id_string), &
640               FORM='UNFORMATTED')
641          danz = danz + 1
[1468]642          WRITE (id_string,'(I6.6)')  danz
[108]643          INQUIRE ( &
644               FILE='PLOT3D_DATA'//TRIM( model_string )//'_'//TRIM(id_string), &
645               EXIST=found )
[1]646
[108]647       ENDDO
[1]648
649#if defined( __netcdf )
[493]650       IF ( netcdf_output )  THEN
651          IF ( netcdf_parallel )  THEN
652             PRINT*, '    NetCDF data are in one file (NetCDF4-format)', &
653                          ' - merging not neccessary'
654          ELSE
655             PRINT*, '    NetCDF output enabled'
656          ENDIF
657       ENDIF
[1]658#else
659       IF ( netcdf_output )  THEN
660          PRINT*, '--- Sorry, no NetCDF support on this host'
661          netcdf_output = .FALSE.
662       ENDIF
663#endif
[493]664       IF ( .NOT. netcdf_parallel )  THEN
665          IF ( danz /= 0 )  THEN
666             PRINT*, '    3D-data:     ', danz, ' file(s) found'
[108]667          ELSE
[2512]668             PRINT*, '    no 3D-data file available'
[108]669          ENDIF
[1]670       ENDIF
671
[493]672       IF ( netcdf_output  .AND. .NOT. netcdf_parallel  .AND.  danz /= 0 )  THEN
[1]673#if defined( __netcdf )
[108]674          DO  av = 0, 1
[1]675
[108]676             IF ( av == 0  .AND.  .NOT.  netcdf_0 )  CYCLE
677             IF ( av == 1  .AND.  .NOT.  netcdf_1 )  CYCLE
[1]678
679!
[108]680!--          Open NetCDF dataset
681             IF ( av == 0 )  THEN
682                filename = 'DATA_3D_NETCDF'//TRIM( model_string )
683             ELSE
684                filename = 'DATA_3D_AV_NETCDF'//TRIM( model_string )
685             ENDIF
686             nc_stat = NF90_OPEN( filename, NF90_WRITE, id_set(av) )
687             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 13 )
[1]688
689
690!
[108]691!--          Get the list of variables (order of variables corresponds with the
692!--          order of data on the binary file)
693             var_list(av) = ' '    ! GET_ATT does not assign trailing blanks
694             nc_stat = NF90_GET_ATT( id_set(av), NF90_GLOBAL, 'VAR_LIST', &
695                  var_list(av) )
696             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 14 )
[1]697
698!
[108]699!--          Inquire id of the time coordinate variable
700             nc_stat = NF90_INQ_VARID( id_set(av), 'time', id_var_time(av) )
701             IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 15 )
[1]702
703!
[108]704!--          Count number of variables; there is one more semicolon in the
705!--          string than variable names
706             num_var(av) = -1
707             DO  i = 1, LEN( var_list(av) )
708                IF ( var_list(av)(i:i) == ';' )  num_var(av) = num_var(av) + 1
709             ENDDO
[1]710
711!
[108]712!--          Extract the variable names from the list and inquire their NetCDF
713!--          IDs
714             pos = INDEX( var_list(av), ';' )
[1]715!
[108]716!--          Loop over all variables
717             DO  i = 1, num_var(av)
[1]718
719!
[108]720!--             Extract variable name from list
721                var_list(av) = var_list(av)(pos+1:)
722                pos = INDEX( var_list(av), ';' )
723                var_name = var_list(av)(1:pos-1)
[1]724
725!
[108]726!--             Get variable ID from name
727                nc_stat = NF90_INQ_VARID( id_set(av), TRIM( var_name ), &
728                     id_var(av,i) )
729                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 16 )
[1]730
[108]731             ENDDO
[1]732
[108]733          ENDDO    ! av=0,1
[1]734
[179]735#endif
[108]736       ENDIF
[1]737
738!
[108]739!--    Read arrays, until the end of the file is reached
[493]740       IF ( .NOT. netcdf_parallel )  THEN
[1]741
[493]742          current_var = 999999999
743          fanz = 0
744          DO  WHILE ( danz /= 0 )
745
[1]746!
[493]747!--          Loop over all files
748             DO  id = 0, danz-1
[1]749!
[2512]750!--             File from PE0 contains at the beginning the index bounds
751!--             of PALM's total domain.
752!--             Allocate the array for storing the total domain data
[493]753                IF ( id == 0  .AND.  fanz(0) == 0  .AND.  fanz(1) == 0 )  THEN
[2512]754!                   READ ( id+110 )  nxag, nxeg, nyag, nyeg
[493]755                   READ ( id+110 )  nxa, nxe, nya, nye, nza, nze
756                   ALLOCATE ( pf3d(nxa:nxe,nya:nye,nza:nze) )
[108]757                ENDIF
[1]758
759!
[493]760!--             Read output time
761                IF ( netcdf_output  .AND.  id == 0 )  THEN
762                   IF ( netcdf_1 )  THEN
763                      READ ( id+110, END=999 )  simulated_time, time_step, av
764                   ELSE
[1]765!
[493]766!--                   For compatibility with earlier PALM versions
767                      READ ( id+110, END=999 )  simulated_time, time_step
768                      av = 0
769                   ENDIF
[108]770                ENDIF
[1]771
772!
[493]773!--             Read subdomain indices and grid point values
774                READ ( id+110, END=999 )  xa, xe, ya, ye, za, ze
775                ALLOCATE( pf3d_tmp(xa:xe,ya:ye,za:ze) )
776                READ ( id+110 )  pf3d_tmp
[210]777
[493]778                xxa = MAX( nxa, xa )
779                xxe = MIN( nxe, xe )
780                yya = MAX( nya, ya )
781                yye = MIN( nye, ye )
[1551]782                DO  k = za, ze
[493]783                   DO  j = yya, yye
784                      DO  i = xxa, xxe
785                         pf3d(i,j,k) = pf3d_tmp(i,j,k)
786                      ENDDO
[210]787                   ENDDO
788                ENDDO
789
[493]790                DEALLOCATE( pf3d_tmp )
791                IF ( id == 0 )  fanz(av) = fanz(av) + 1
[1]792
[493]793             ENDDO
[1]794
795!
[2512]796!--          Write data of the total domain in NetCDF format
[493]797             IF ( netcdf_output )  THEN
[1]798#if defined( __netcdf )
799!
[493]800!--             Check if a new time step has begun; if yes write data to time
801!--             axis
802                IF ( current_var(av) > num_var(av) )  THEN
803                   current_var(av) = 1
804                   nc_stat = NF90_PUT_VAR( id_set(av), id_var_time(av), &
805                                      (/ simulated_time /),&
806                                      start = (/ time_step /), count = (/ 1 /) )
807                   IF ( nc_stat /= NF90_NOERR ) CALL handle_netcdf_error( 17 )
808                ENDIF
[1]809
810!
[493]811!--             Now write the data
812                nc_stat = NF90_PUT_VAR( id_set(av), id_var(av,current_var(av)),&
[1551]813                                        pf3d(nxa:nxe,nya:nye,za:ze), start = (/ 1, 1, 1, time_step /),&
814                              count = (/ nxe-nxa+1, nye-nya+1, ze-za+1, 1 /) )
[493]815                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 18 )
[1]816
[493]817                current_var(av) = current_var(av) + 1
[1]818
819#endif
[493]820             ENDIF
[1]821
[493]822          ENDDO
[1]823
[493]824       ENDIF
825
826999    IF ( danz /= 0  .AND.  .NOT. netcdf_parallel )  THEN
[1]827!
[108]828!--       Print the number of arrays processed
829          WRITE (*,'(16X,I4,A)')  fanz(0)+fanz(1), ' array(s) processed'
830          IF ( fanz(1) /= 0 )  THEN
831             WRITE (*,'(16X,I4,A)')  fanz(1), ' array(s) are time-averaged'
832          ENDIF
[1]833!
[108]834!--       Close all files and deallocate array
835          DO  id = 0, danz-1
836             CLOSE ( id+110 )
837          ENDDO
838          CLOSE ( 2 )
839          DEALLOCATE ( pf3d )
[1]840!
[108]841!--       Close the NetCDF file
842          IF ( netcdf_output )  THEN
[1]843#if defined( __netcdf )
[108]844             IF ( netcdf_0 )  THEN
845                nc_stat = NF90_CLOSE( id_set(0) )
846                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 19 )
847             ENDIF
848             IF ( netcdf_1 )  THEN
849                nc_stat = NF90_CLOSE( id_set(1) )
850                IF ( nc_stat /= NF90_NOERR )  CALL handle_netcdf_error( 20 )
851             ENDIF
852#endif
[1]853          ENDIF
[493]854
855!
856!--       Output required cpu time
857          CALL SYSTEM_CLOCK( count, count_rate )
858          cpu_end_time = REAL( count ) / REAL( count_rate )
859          WRITE (*,'(5X,A,F9.3,A)')  'Required cpu-time: ', &
860                                     cpu_end_time-cpu_start_time, ' sec'
861
[1]862       ENDIF
863
[108]864    ENDDO  ! models
[1]865
[108]866
[1]867 CONTAINS
868
869
870    SUBROUTINE handle_netcdf_error( errno )
871!
872!--    Prints out a text message corresponding to the current NetCDF status
873
874       IMPLICIT NONE
875
876       INTEGER, INTENT(IN) ::  errno
877
[179]878#if defined( __netcdf )
[1]879       IF ( nc_stat /= NF90_NOERR )  THEN
880          PRINT*, '+++ combine_plot_fields  netcdf: ', av, errno, &
881                  TRIM( nf90_strerror( nc_stat ) )
882       ENDIF
[179]883#endif
[1]884
885    END SUBROUTINE handle_netcdf_error
886
887
888 END PROGRAM combine_plot_fields
889
890
891
Note: See TracBrowser for help on using the repository browser.