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

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

New:
---

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

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

For unknown variables (CASE DEFAULT) call new subroutine user_data_output_dvrp

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

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

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

Changed:


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

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

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

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

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

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

non-zero wall_heatflux,

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

humidity/scalar/salinity in a future release.

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

Errors:


Bugfix: summation of sums_l_l in diffusivities.

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

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

Bugfix: model_string needed a default value (combine_plot_fields)

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

Bugfix: assignment of fluxes at walls

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

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

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