source: palm/trunk/SOURCE/data_output_2d.f90 @ 3295

Last change on this file since 3295 was 3294, checked in by raasch, 6 years ago

modularization of the ocean code

  • Property svn:keywords set to Id
File size: 94.2 KB
Line 
1!> @file data_output_2d.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: data_output_2d.f90 3294 2018-10-01 02:37:10Z raasch $
27! changes concerning modularization of ocean option
28!
29! 3274 2018-09-24 15:42:55Z knoop
30! Modularization of all bulk cloud physics code components
31!
32! 3241 2018-09-12 15:02:00Z raasch
33! unused variables removed
34!
35! 3176 2018-07-26 17:12:48Z suehring
36! Remove output of latent heat flux at urban-surfaces and set fill values
37! instead
38!
39! 3052 2018-05-31 06:11:20Z raasch
40! Do not open FORTRAN binary files in case of parallel netCDF I/O
41!
42! 3049 2018-05-29 13:52:36Z Giersch
43! Error messages revised
44!
45! 3045 2018-05-28 07:55:41Z Giersch
46! Error messages revised
47!
48! 3014 2018-05-09 08:42:38Z maronga
49! Added nzb_do and nzt_do for some modules for 2d output
50!
51! 3004 2018-04-27 12:33:25Z Giersch
52! precipitation_rate removed, case prr*_xy removed, to_be_resorted have to point
53! to ql_vp_av and not to ql_vp, allocation checks implemented (averaged data
54! will be assigned to fill values if no allocation happened so far)   
55!
56! 2963 2018-04-12 14:47:44Z suehring
57! Introduce index for vegetation/wall, pavement/green-wall and water/window
58! surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
59!
60! 2817 2018-02-19 16:32:21Z knoop
61! Preliminary gust module interface implemented
62!
63! 2805 2018-02-14 17:00:09Z suehring
64! Consider also default-type surfaces for surface temperature output.
65!
66! 2797 2018-02-08 13:24:35Z suehring
67! Enable output of ground-heat flux also at urban surfaces.
68!
69! 2743 2018-01-12 16:03:39Z suehring
70! In case of natural- and urban-type surfaces output surfaces fluxes in W/m2.
71!
72! 2742 2018-01-12 14:59:47Z suehring
73! Enable output of surface temperature
74!
75! 2735 2018-01-11 12:01:27Z suehring
76! output of r_a moved from land-surface to consider also urban-type surfaces
77!
78! 2718 2018-01-02 08:49:38Z maronga
79! Corrected "Former revisions" section
80!
81! 2696 2017-12-14 17:12:51Z kanani
82! Change in file header (GPL part)
83! Implementation of uv exposure model (FK)
84! Implementation of turbulence_closure_mod (TG)
85! Set fill values at topography grid points or e.g. non-natural-type surface
86! in case of LSM output (MS)
87!
88! 2512 2017-10-04 08:26:59Z raasch
89! upper bounds of cross section output changed from nx+1,ny+1 to nx,ny
90! no output of ghost layer data
91!
92! 2292 2017-06-20 09:51:42Z schwenkel
93! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
94! includes two more prognostic equations for cloud drop concentration (nc) 
95! and cloud water content (qc).
96!
97! 2277 2017-06-12 10:47:51Z kanani
98! Removed unused variables do2d_xy_n, do2d_xz_n, do2d_yz_n
99!
100! 2233 2017-05-30 18:08:54Z suehring
101!
102! 2232 2017-05-30 17:47:52Z suehring
103! Adjustments to new surface concept
104!
105!
106! 2190 2017-03-21 12:16:43Z raasch
107! bugfix for r2031: string rho replaced by rho_ocean
108!
109! 2031 2016-10-21 15:11:58Z knoop
110! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
111!
112! 2000 2016-08-20 18:09:15Z knoop
113! Forced header and separation lines into 80 columns
114!
115! 1980 2016-07-29 15:51:57Z suehring
116! Bugfix, in order to steer user-defined output, setting flag found explicitly
117! to .F.
118!
119! 1976 2016-07-27 13:28:04Z maronga
120! Output of radiation quantities is now done directly in the respective module
121!
122! 1972 2016-07-26 07:52:02Z maronga
123! Output of land surface quantities is now done directly in the respective
124! module
125!
126! 1960 2016-07-12 16:34:24Z suehring
127! Scalar surface flux added
128! Rename INTEGER variable s into s_ind, as s is already assigned to scalar
129!
130! 1849 2016-04-08 11:33:18Z hoffmann
131! precipitation_amount, precipitation_rate, prr moved to arrays_3d
132!
133! 1822 2016-04-07 07:49:42Z hoffmann
134! Output of bulk cloud physics simplified.
135!
136! 1788 2016-03-10 11:01:04Z maronga
137! Added output of z0q
138!
139! 1783 2016-03-06 18:36:17Z raasch
140! name change of netcdf routines and module + related changes
141!
142! 1745 2016-02-05 13:06:51Z gronemeier
143! Bugfix: test if time axis limit exceeds moved to point after call of check_open
144!
145! 1703 2015-11-02 12:38:44Z raasch
146! bugfix for output of single (*) xy-sections in case of parallel netcdf I/O
147!
148! 1701 2015-11-02 07:43:04Z maronga
149! Bugfix in output of RRTGM data
150!
151! 1691 2015-10-26 16:17:44Z maronga
152! Added output of Obukhov length (ol) and radiative heating rates  for RRTMG.
153! Formatting corrections.
154!
155! 1682 2015-10-07 23:56:08Z knoop
156! Code annotations made doxygen readable
157!
158! 1585 2015-04-30 07:05:52Z maronga
159! Added support for RRTMG
160!
161! 1555 2015-03-04 17:44:27Z maronga
162! Added output of r_a and r_s
163!
164! 1551 2015-03-03 14:18:16Z maronga
165! Added suppport for land surface model and radiation model output. In the course
166! of this action, the limits for vertical loops have been changed (from nzb and
167! nzt+1 to nzb_do and nzt_do, respectively in order to allow soil model output).
168! Moreover, a new vertical grid zs was introduced.
169!
170! 1359 2014-04-11 17:15:14Z hoffmann
171! New particle structure integrated.
172!
173! 1353 2014-04-08 15:21:23Z heinze
174! REAL constants provided with KIND-attribute
175!
176! 1327 2014-03-21 11:00:16Z raasch
177! parts concerning iso2d output removed,
178! -netcdf output queries
179!
180! 1320 2014-03-20 08:40:49Z raasch
181! ONLY-attribute added to USE-statements,
182! kind-parameters added to all INTEGER and REAL declaration statements,
183! kinds are defined in new module kinds,
184! revision history before 2012 removed,
185! comment fields (!:) to be used for variable explanations added to
186! all variable declaration statements
187!
188! 1318 2014-03-17 13:35:16Z raasch
189! barrier argument removed from cpu_log.
190! module interfaces removed
191!
192! 1311 2014-03-14 12:13:39Z heinze
193! bugfix: close #if defined( __netcdf )
194!
195! 1308 2014-03-13 14:58:42Z fricke
196! +local_2d_sections, local_2d_sections_l, ns
197! Check, if the limit of the time dimension is exceeded for parallel output
198! To increase the performance for parallel output, the following is done:
199! - Update of time axis is only done by PE0
200! - Cross sections are first stored on a local array and are written
201!   collectively to the output file by all PEs.
202!
203! 1115 2013-03-26 18:16:16Z hoffmann
204! ql is calculated by calc_liquid_water_content
205!
206! 1076 2012-12-05 08:30:18Z hoffmann
207! Bugfix in output of ql
208!
209! 1065 2012-11-22 17:42:36Z hoffmann
210! Bugfix: Output of cross sections of ql
211!
212! 1053 2012-11-13 17:11:03Z hoffmann
213! +qr, nr, qc and cross sections
214!
215! 1036 2012-10-22 13:43:42Z raasch
216! code put under GPL (PALM 3.9)
217!
218! 1031 2012-10-19 14:35:30Z raasch
219! netCDF4 without parallel file support implemented
220!
221! 1007 2012-09-19 14:30:36Z franke
222! Bugfix: missing calculation of ql_vp added
223!
224! 978 2012-08-09 08:28:32Z fricke
225! +z0h
226!
227! Revision 1.1  1997/08/11 06:24:09  raasch
228! Initial revision
229!
230!
231! Description:
232! ------------
233!> Data output of cross-sections in netCDF format or binary format
234!> to be later converted to NetCDF by helper routine combine_plot_fields.
235!> Attention: The position of the sectional planes is still not always computed
236!> ---------  correctly. (zu is used always)!
237!------------------------------------------------------------------------------!
238 SUBROUTINE data_output_2d( mode, av )
239 
240
241    USE arrays_3d,                                                             &
242        ONLY:  dzw, e, heatflux_output_conversion, nc, nr, p, pt,              &
243               precipitation_amount, prr, q, qc, ql, ql_c, ql_v, qr, s, tend,  &
244               u, v, vpt, w, zu, zw, waterflux_output_conversion, hyrho, d_exner
245
246    USE averaging
247
248    USE basic_constants_and_equations_mod,                                     &
249        ONLY:  c_p, lv_d_cp, l_v
250
251    USE bulk_cloud_model_mod,                                                  &
252        ONLY:  bulk_cloud_model, bcm_data_output_2d
253
254    USE control_parameters,                                                    &
255        ONLY:  data_output_2d_on_each_pe, data_output_xy,                      &
256               data_output_xz, data_output_yz, do2d,                           &
257               do2d_xy_last_time, do2d_xy_time_count,                          &
258               do2d_xz_last_time, do2d_xz_time_count,                          &
259               do2d_yz_last_time, do2d_yz_time_count,                          &
260               ibc_uv_b, io_blocks, io_group, land_surface, message_string,    &
261               ntdim_2d_xy, ntdim_2d_xz, ntdim_2d_yz,                          &
262               ocean_mode, psolver, section, simulated_time,                   &
263               time_since_reference_point, uv_exposure
264
265    USE cpulog,                                                                &
266        ONLY:  cpu_log, log_point 
267
268    USE gust_mod,                                                              &
269        ONLY:  gust_data_output_2d, gust_module_enabled
270
271    USE indices,                                                               &
272        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,       &
273               nzb, nzt, wall_flags_0
274
275    USE kinds
276
277    USE land_surface_model_mod,                                                &
278        ONLY:  lsm_data_output_2d, zs
279
280#if defined( __netcdf )
281    USE NETCDF
282#endif
283
284    USE netcdf_interface,                                                      &
285        ONLY:  fill_value, id_set_xy, id_set_xz, id_set_yz, id_var_do2d,       &
286               id_var_time_xy, id_var_time_xz, id_var_time_yz, nc_stat,        &
287               netcdf_data_format, netcdf_handle_error
288
289    USE ocean_mod,                                                             &
290        ONLY:  ocean_data_output_2d
291
292    USE particle_attributes,                                                   &
293        ONLY:  grid_particles, number_of_particles, particle_advection_start,  &
294               particles, prt_count
295   
296    USE pegrid
297
298    USE radiation_model_mod,                                                   &
299        ONLY:  radiation, radiation_data_output_2d
300
301    USE surface_mod,                                                           &
302        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_def_h,           &
303               surf_lsm_h, surf_usm_h
304
305    USE turbulence_closure_mod,                                                &
306        ONLY:  tcm_data_output_2d
307
308    USE uv_exposure_model_mod,                                                 &
309        ONLY:  uvem_data_output_2d
310
311
312    IMPLICIT NONE
313
314    CHARACTER (LEN=2)  ::  do2d_mode    !<
315    CHARACTER (LEN=2)  ::  mode         !<
316    CHARACTER (LEN=4)  ::  grid         !<
317    CHARACTER (LEN=50) ::  rtext        !<
318   
319    INTEGER(iwp) ::  av        !<
320    INTEGER(iwp) ::  ngp       !<
321    INTEGER(iwp) ::  file_id   !<
322    INTEGER(iwp) ::  flag_nr   !< number of masking flag
323    INTEGER(iwp) ::  i         !<
324    INTEGER(iwp) ::  if        !<
325    INTEGER(iwp) ::  is        !<
326    INTEGER(iwp) ::  iis       !<
327    INTEGER(iwp) ::  j         !<
328    INTEGER(iwp) ::  k         !<
329    INTEGER(iwp) ::  l         !<
330    INTEGER(iwp) ::  layer_xy  !<
331    INTEGER(iwp) ::  m         !<
332    INTEGER(iwp) ::  n         !<
333    INTEGER(iwp) ::  nis       !<
334    INTEGER(iwp) ::  ns        !<
335    INTEGER(iwp) ::  nzb_do    !< lower limit of the data field (usually nzb)
336    INTEGER(iwp) ::  nzt_do    !< upper limit of the data field (usually nzt+1)
337    INTEGER(iwp) ::  s_ind     !<
338    INTEGER(iwp) ::  sender    !<
339    INTEGER(iwp) ::  ind(4)    !<
340   
341    LOGICAL ::  found          !<
342    LOGICAL ::  resorted       !<
343    LOGICAL ::  two_d          !<
344   
345    REAL(wp) ::  mean_r        !<
346    REAL(wp) ::  s_r2          !<
347    REAL(wp) ::  s_r3          !<
348   
349    REAL(wp), DIMENSION(:), ALLOCATABLE     ::  level_z             !<
350    REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  local_2d            !<
351    REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  local_2d_l          !<
352
353    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  local_pf            !<
354    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  local_2d_sections   !<
355    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  local_2d_sections_l !<
356
357#if defined( __parallel )
358    REAL(wp), DIMENSION(:,:),   ALLOCATABLE ::  total_2d    !<
359#endif
360    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !<
361
362    NAMELIST /LOCAL/  rtext
363
364!
365!-- Immediate return, if no output is requested (no respective sections
366!-- found in parameter data_output)
367    IF ( mode == 'xy'  .AND.  .NOT. data_output_xy(av) )  RETURN
368    IF ( mode == 'xz'  .AND.  .NOT. data_output_xz(av) )  RETURN
369    IF ( mode == 'yz'  .AND.  .NOT. data_output_yz(av) )  RETURN
370
371    CALL cpu_log (log_point(3),'data_output_2d','start')
372
373    two_d = .FALSE.    ! local variable to distinguish between output of pure 2D
374                       ! arrays and cross-sections of 3D arrays.
375
376!
377!-- Depending on the orientation of the cross-section, the respective output
378!-- files have to be opened.
379    SELECT CASE ( mode )
380
381       CASE ( 'xy' )
382          s_ind = 1
383          ALLOCATE( level_z(nzb:nzt+1), local_2d(nxl:nxr,nys:nyn) )
384
385          IF ( netcdf_data_format > 4 )  THEN
386             ns = 1
387             DO WHILE ( section(ns,s_ind) /= -9999  .AND.  ns <= 100 )
388                ns = ns + 1
389             ENDDO
390             ns = ns - 1
391             ALLOCATE( local_2d_sections(nxl:nxr,nys:nyn,1:ns) )
392             local_2d_sections = 0.0_wp
393          ENDIF
394
395!
396!--       Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
397          IF ( myid == 0  .OR.  netcdf_data_format > 4 )  THEN
398             CALL check_open( 101+av*10 )
399          ENDIF
400          IF ( data_output_2d_on_each_pe  .AND.  netcdf_data_format < 5 )  THEN
401             CALL check_open( 21 )
402          ELSE
403             IF ( myid == 0 )  THEN
404#if defined( __parallel )
405                ALLOCATE( total_2d(0:nx,0:ny) )
406#endif
407             ENDIF
408          ENDIF
409
410       CASE ( 'xz' )
411          s_ind = 2
412          ALLOCATE( local_2d(nxl:nxr,nzb:nzt+1) )
413
414          IF ( netcdf_data_format > 4 )  THEN
415             ns = 1
416             DO WHILE ( section(ns,s_ind) /= -9999  .AND.  ns <= 100 )
417                ns = ns + 1
418             ENDDO
419             ns = ns - 1
420             ALLOCATE( local_2d_sections(nxl:nxr,1:ns,nzb:nzt+1) )
421             ALLOCATE( local_2d_sections_l(nxl:nxr,1:ns,nzb:nzt+1) )
422             local_2d_sections = 0.0_wp; local_2d_sections_l = 0.0_wp
423          ENDIF
424
425!
426!--       Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
427          IF ( myid == 0 .OR. netcdf_data_format > 4 )  THEN
428             CALL check_open( 102+av*10 )
429          ENDIF
430
431          IF ( data_output_2d_on_each_pe  .AND.  netcdf_data_format < 5 )  THEN
432             CALL check_open( 22 )
433          ELSE
434             IF ( myid == 0 )  THEN
435#if defined( __parallel )
436                ALLOCATE( total_2d(0:nx,nzb:nzt+1) )
437#endif
438             ENDIF
439          ENDIF
440
441       CASE ( 'yz' )
442          s_ind = 3
443          ALLOCATE( local_2d(nys:nyn,nzb:nzt+1) )
444
445          IF ( netcdf_data_format > 4 )  THEN
446             ns = 1
447             DO WHILE ( section(ns,s_ind) /= -9999  .AND.  ns <= 100 )
448                ns = ns + 1
449             ENDDO
450             ns = ns - 1
451             ALLOCATE( local_2d_sections(1:ns,nys:nyn,nzb:nzt+1) )
452             ALLOCATE( local_2d_sections_l(1:ns,nys:nyn,nzb:nzt+1) )
453             local_2d_sections = 0.0_wp; local_2d_sections_l = 0.0_wp
454          ENDIF
455
456!
457!--       Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
458          IF ( myid == 0 .OR. netcdf_data_format > 4 )  THEN
459             CALL check_open( 103+av*10 )
460          ENDIF
461
462          IF ( data_output_2d_on_each_pe  .AND.  netcdf_data_format < 5 )  THEN
463             CALL check_open( 23 )
464          ELSE
465             IF ( myid == 0 )  THEN
466#if defined( __parallel )
467                ALLOCATE( total_2d(0:ny,nzb:nzt+1) )
468#endif
469             ENDIF
470          ENDIF
471
472       CASE DEFAULT
473          message_string = 'unknown cross-section: ' // TRIM( mode )
474          CALL message( 'data_output_2d', 'PA0180', 1, 2, 0, 6, 0 )
475
476    END SELECT
477
478!
479!-- For parallel netcdf output the time axis must be limited. Return, if this
480!-- limit is exceeded. This could be the case, if the simulated time exceeds
481!-- the given end time by the length of the given output interval.
482    IF ( netcdf_data_format > 4 )  THEN
483       IF ( mode == 'xy'  .AND.  do2d_xy_time_count(av) + 1 >                  &
484            ntdim_2d_xy(av) )  THEN
485          WRITE ( message_string, * ) 'Output of xy cross-sections is not ',   &
486                          'given at t=', simulated_time, '&because the',       & 
487                          ' maximum number of output time levels is exceeded.'
488          CALL message( 'data_output_2d', 'PA0384', 0, 1, 0, 6, 0 )
489          CALL cpu_log( log_point(3), 'data_output_2d', 'stop' )
490          RETURN
491       ENDIF
492       IF ( mode == 'xz'  .AND.  do2d_xz_time_count(av) + 1 >                  &
493            ntdim_2d_xz(av) )  THEN
494          WRITE ( message_string, * ) 'Output of xz cross-sections is not ',   &
495                          'given at t=', simulated_time, '&because the',       & 
496                          ' maximum number of output time levels is exceeded.'
497          CALL message( 'data_output_2d', 'PA0385', 0, 1, 0, 6, 0 )
498          CALL cpu_log( log_point(3), 'data_output_2d', 'stop' )
499          RETURN
500       ENDIF
501       IF ( mode == 'yz'  .AND.  do2d_yz_time_count(av) + 1 >                  &
502            ntdim_2d_yz(av) )  THEN
503          WRITE ( message_string, * ) 'Output of yz cross-sections is not ',   &
504                          'given at t=', simulated_time, '&because the',       & 
505                          ' maximum number of output time levels is exceeded.'
506          CALL message( 'data_output_2d', 'PA0386', 0, 1, 0, 6, 0 )
507          CALL cpu_log( log_point(3), 'data_output_2d', 'stop' )
508          RETURN
509       ENDIF
510    ENDIF
511
512!
513!-- Allocate a temporary array for resorting (kji -> ijk).
514    ALLOCATE( local_pf(nxl:nxr,nys:nyn,nzb:nzt+1) )
515    local_pf = 0.0
516
517!
518!-- Loop of all variables to be written.
519!-- Output dimensions chosen
520    if = 1
521    l = MAX( 2, LEN_TRIM( do2d(av,if) ) )
522    do2d_mode = do2d(av,if)(l-1:l)
523
524    DO  WHILE ( do2d(av,if)(1:1) /= ' ' )
525
526       IF ( do2d_mode == mode )  THEN
527!
528!--       Set flag to steer output of radiation, land-surface, or user-defined
529!--       quantities
530          found = .FALSE.
531
532          nzb_do = nzb
533          nzt_do = nzt+1
534!
535!--       Before each output, set array local_pf to fill value
536          local_pf = fill_value
537!
538!--       Set masking flag for topography for not resorted arrays
539          flag_nr = 0
540         
541!
542!--       Store the array chosen on the temporary array.
543          resorted = .FALSE.
544          SELECT CASE ( TRIM( do2d(av,if) ) )
545             CASE ( 'e_xy', 'e_xz', 'e_yz' )
546                IF ( av == 0 )  THEN
547                   to_be_resorted => e
548                ELSE
549                   IF ( .NOT. ALLOCATED( e_av ) ) THEN
550                      ALLOCATE( e_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
551                      e_av = REAL( fill_value, KIND = wp )
552                   ENDIF
553                   to_be_resorted => e_av
554                ENDIF
555                IF ( mode == 'xy' )  level_z = zu
556
557             CASE ( 'lpt_xy', 'lpt_xz', 'lpt_yz' )
558                IF ( av == 0 )  THEN
559                   to_be_resorted => pt
560                ELSE
561                   IF ( .NOT. ALLOCATED( lpt_av ) ) THEN
562                      ALLOCATE( lpt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
563                      lpt_av = REAL( fill_value, KIND = wp )
564                   ENDIF
565                   to_be_resorted => lpt_av
566                ENDIF
567                IF ( mode == 'xy' )  level_z = zu
568
569             CASE ( 'lwp*_xy' )        ! 2d-array
570                IF ( av == 0 )  THEN
571                   DO  i = nxl, nxr
572                      DO  j = nys, nyn
573                         local_pf(i,j,nzb+1) = SUM( ql(nzb:nzt,j,i) *          &
574                                                    dzw(1:nzt+1) )
575                      ENDDO
576                   ENDDO
577                ELSE
578                   IF ( .NOT. ALLOCATED( lwp_av ) ) THEN
579                      ALLOCATE( lwp_av(nysg:nyng,nxlg:nxrg) )
580                      lwp_av = REAL( fill_value, KIND = wp )
581                   ENDIF
582                   DO  i = nxl, nxr
583                      DO  j = nys, nyn
584                         local_pf(i,j,nzb+1) = lwp_av(j,i)
585                      ENDDO
586                   ENDDO
587                ENDIF
588                resorted = .TRUE.
589                two_d = .TRUE.
590                level_z(nzb+1) = zu(nzb+1)
591
592             CASE ( 'ghf*_xy' )        ! 2d-array
593                IF ( av == 0 )  THEN
594                   DO  m = 1, surf_lsm_h%ns
595                      i                   = surf_lsm_h%i(m)           
596                      j                   = surf_lsm_h%j(m)
597                      local_pf(i,j,nzb+1) = surf_lsm_h%ghf(m)
598                   ENDDO
599                   DO  m = 1, surf_usm_h%ns
600                      i                   = surf_usm_h%i(m)           
601                      j                   = surf_usm_h%j(m)
602                      local_pf(i,j,nzb+1) = surf_usm_h%frac(ind_veg_wall,m)  *  &
603                                            surf_usm_h%wghf_eb(m)        +      &
604                                            surf_usm_h%frac(ind_pav_green,m) *  &
605                                            surf_usm_h%wghf_eb_green(m)  +      &
606                                            surf_usm_h%frac(ind_wat_win,m)   *  &
607                                            surf_usm_h%wghf_eb_window(m)
608                   ENDDO
609                ELSE
610                   IF ( .NOT. ALLOCATED( ghf_av ) ) THEN
611                      ALLOCATE( ghf_av(nysg:nyng,nxlg:nxrg) )
612                      ghf_av = REAL( fill_value, KIND = wp )
613                   ENDIF
614                   DO  i = nxl, nxr
615                      DO  j = nys, nyn
616                         local_pf(i,j,nzb+1) = ghf_av(j,i)
617                      ENDDO
618                   ENDDO
619                ENDIF
620
621                resorted = .TRUE.
622                two_d = .TRUE.
623                level_z(nzb+1) = zu(nzb+1)
624
625             CASE ( 'ol*_xy' )        ! 2d-array
626                IF ( av == 0 ) THEN
627                   DO  m = 1, surf_def_h(0)%ns
628                      i = surf_def_h(0)%i(m)
629                      j = surf_def_h(0)%j(m)
630                      local_pf(i,j,nzb+1) = surf_def_h(0)%ol(m)
631                   ENDDO
632                   DO  m = 1, surf_lsm_h%ns
633                      i = surf_lsm_h%i(m)
634                      j = surf_lsm_h%j(m)
635                      local_pf(i,j,nzb+1) = surf_lsm_h%ol(m)
636                   ENDDO
637                   DO  m = 1, surf_usm_h%ns
638                      i = surf_usm_h%i(m)
639                      j = surf_usm_h%j(m)
640                      local_pf(i,j,nzb+1) = surf_usm_h%ol(m)
641                   ENDDO
642                ELSE
643                   IF ( .NOT. ALLOCATED( ol_av ) ) THEN
644                      ALLOCATE( ol_av(nysg:nyng,nxlg:nxrg) )
645                      ol_av = REAL( fill_value, KIND = wp )
646                   ENDIF
647                   DO  i = nxl, nxr
648                      DO  j = nys, nyn
649                         local_pf(i,j,nzb+1) = ol_av(j,i)
650                      ENDDO
651                   ENDDO
652                ENDIF
653                resorted = .TRUE.
654                two_d = .TRUE.
655                level_z(nzb+1) = zu(nzb+1)
656
657             CASE ( 'p_xy', 'p_xz', 'p_yz' )
658                IF ( av == 0 )  THEN
659                   IF ( psolver /= 'sor' )  CALL exchange_horiz( p, nbgp )
660                   to_be_resorted => p
661                ELSE
662                   IF ( .NOT. ALLOCATED( p_av ) ) THEN
663                      ALLOCATE( p_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
664                      p_av = REAL( fill_value, KIND = wp )
665                   ENDIF
666                   IF ( psolver /= 'sor' )  CALL exchange_horiz( p_av, nbgp )
667                   to_be_resorted => p_av
668                ENDIF
669                IF ( mode == 'xy' )  level_z = zu
670
671             CASE ( 'pc_xy', 'pc_xz', 'pc_yz' )  ! particle concentration
672                IF ( av == 0 )  THEN
673                   IF ( simulated_time >= particle_advection_start )  THEN
674                      tend = prt_count
675!                      CALL exchange_horiz( tend, nbgp )
676                   ELSE
677                      tend = 0.0_wp
678                   ENDIF
679                   DO  i = nxl, nxr
680                      DO  j = nys, nyn
681                         DO  k = nzb, nzt+1
682                            local_pf(i,j,k) = tend(k,j,i)
683                         ENDDO
684                      ENDDO
685                   ENDDO
686                   resorted = .TRUE.
687                ELSE
688                   IF ( .NOT. ALLOCATED( pc_av ) ) THEN
689                      ALLOCATE( pc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
690                      pc_av = REAL( fill_value, KIND = wp )
691                   ENDIF
692!                   CALL exchange_horiz( pc_av, nbgp )
693                   to_be_resorted => pc_av
694                ENDIF
695
696             CASE ( 'pr_xy', 'pr_xz', 'pr_yz' )  ! mean particle radius (effective radius)
697                IF ( av == 0 )  THEN
698                   IF ( simulated_time >= particle_advection_start )  THEN
699                      DO  i = nxl, nxr
700                         DO  j = nys, nyn
701                            DO  k = nzb, nzt+1
702                               number_of_particles = prt_count(k,j,i)
703                               IF (number_of_particles <= 0)  CYCLE
704                               particles => grid_particles(k,j,i)%particles(1:number_of_particles)
705                               s_r2 = 0.0_wp
706                               s_r3 = 0.0_wp
707                               DO  n = 1, number_of_particles
708                                  IF ( particles(n)%particle_mask )  THEN
709                                     s_r2 = s_r2 + particles(n)%radius**2 * &
710                                            particles(n)%weight_factor
711                                     s_r3 = s_r3 + particles(n)%radius**3 * &
712                                            particles(n)%weight_factor
713                                  ENDIF
714                               ENDDO
715                               IF ( s_r2 > 0.0_wp )  THEN
716                                  mean_r = s_r3 / s_r2
717                               ELSE
718                                  mean_r = 0.0_wp
719                               ENDIF
720                               tend(k,j,i) = mean_r
721                            ENDDO
722                         ENDDO
723                      ENDDO
724!                      CALL exchange_horiz( tend, nbgp )
725                   ELSE
726                      tend = 0.0_wp
727                   ENDIF
728                   DO  i = nxl, nxr
729                      DO  j = nys, nyn
730                         DO  k = nzb, nzt+1
731                            local_pf(i,j,k) = tend(k,j,i)
732                         ENDDO
733                      ENDDO
734                   ENDDO
735                   resorted = .TRUE.
736                ELSE
737                   IF ( .NOT. ALLOCATED( pr_av ) ) THEN
738                      ALLOCATE( pr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
739                      pr_av = REAL( fill_value, KIND = wp )
740                   ENDIF
741!                   CALL exchange_horiz( pr_av, nbgp )
742                   to_be_resorted => pr_av
743                ENDIF
744
745             CASE ( 'pt_xy', 'pt_xz', 'pt_yz' )
746                IF ( av == 0 )  THEN
747                   IF ( .NOT. bulk_cloud_model ) THEN
748                      to_be_resorted => pt
749                   ELSE
750                   DO  i = nxl, nxr
751                      DO  j = nys, nyn
752                            DO  k = nzb, nzt+1
753                               local_pf(i,j,k) = pt(k,j,i) + lv_d_cp *         &
754                                                             d_exner(k) *      &
755                                                             ql(k,j,i)
756                            ENDDO
757                         ENDDO
758                      ENDDO
759                      resorted = .TRUE.
760                   ENDIF
761                ELSE
762                   IF ( .NOT. ALLOCATED( pt_av ) ) THEN
763                      ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
764                      pt_av = REAL( fill_value, KIND = wp )
765                   ENDIF
766                   to_be_resorted => pt_av
767                ENDIF
768                IF ( mode == 'xy' )  level_z = zu
769
770             CASE ( 'q_xy', 'q_xz', 'q_yz' )
771                IF ( av == 0 )  THEN
772                   to_be_resorted => q
773                ELSE
774                   IF ( .NOT. ALLOCATED( q_av ) ) THEN
775                      ALLOCATE( q_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
776                      q_av = REAL( fill_value, KIND = wp )
777                   ENDIF
778                   to_be_resorted => q_av
779                ENDIF
780                IF ( mode == 'xy' )  level_z = zu
781
782             CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
783                IF ( av == 0 )  THEN
784                   to_be_resorted => ql
785                ELSE
786                   IF ( .NOT. ALLOCATED( ql_av ) ) THEN
787                      ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
788                      ql_av = REAL( fill_value, KIND = wp )
789                   ENDIF
790                   to_be_resorted => ql_av
791                ENDIF
792                IF ( mode == 'xy' )  level_z = zu
793
794             CASE ( 'ql_c_xy', 'ql_c_xz', 'ql_c_yz' )
795                IF ( av == 0 )  THEN
796                   to_be_resorted => ql_c
797                ELSE
798                   IF ( .NOT. ALLOCATED( ql_c_av ) ) THEN
799                      ALLOCATE( ql_c_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
800                      ql_c_av = REAL( fill_value, KIND = wp )
801                   ENDIF
802                   to_be_resorted => ql_c_av
803                ENDIF
804                IF ( mode == 'xy' )  level_z = zu
805
806             CASE ( 'ql_v_xy', 'ql_v_xz', 'ql_v_yz' )
807                IF ( av == 0 )  THEN
808                   to_be_resorted => ql_v
809                ELSE
810                   IF ( .NOT. ALLOCATED( ql_v_av ) ) THEN
811                      ALLOCATE( ql_v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
812                      ql_v_av = REAL( fill_value, KIND = wp )
813                   ENDIF
814                   to_be_resorted => ql_v_av
815                ENDIF
816                IF ( mode == 'xy' )  level_z = zu
817
818             CASE ( 'ql_vp_xy', 'ql_vp_xz', 'ql_vp_yz' )
819                IF ( av == 0 )  THEN
820                   IF ( simulated_time >= particle_advection_start )  THEN
821                      DO  i = nxl, nxr
822                         DO  j = nys, nyn
823                            DO  k = nzb, nzt+1
824                               number_of_particles = prt_count(k,j,i)
825                               IF (number_of_particles <= 0)  CYCLE
826                               particles => grid_particles(k,j,i)%particles(1:number_of_particles)
827                               DO  n = 1, number_of_particles
828                                  IF ( particles(n)%particle_mask )  THEN
829                                     tend(k,j,i) =  tend(k,j,i) +                 &
830                                                    particles(n)%weight_factor /  &
831                                                    prt_count(k,j,i)
832                                  ENDIF
833                               ENDDO
834                            ENDDO
835                         ENDDO
836                      ENDDO
837!                      CALL exchange_horiz( tend, nbgp )
838                   ELSE
839                      tend = 0.0_wp
840                   ENDIF
841                   DO  i = nxl, nxr
842                      DO  j = nys, nyn
843                         DO  k = nzb, nzt+1
844                            local_pf(i,j,k) = tend(k,j,i)
845                         ENDDO
846                      ENDDO
847                   ENDDO
848                   resorted = .TRUE.
849                ELSE
850                   IF ( .NOT. ALLOCATED( ql_vp_av ) ) THEN
851                      ALLOCATE( ql_vp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
852                      ql_vp_av = REAL( fill_value, KIND = wp )
853                   ENDIF
854!                   CALL exchange_horiz( ql_vp_av, nbgp )
855                   to_be_resorted => ql_vp_av
856                ENDIF
857                IF ( mode == 'xy' )  level_z = zu
858
859             CASE ( 'qsws*_xy' )        ! 2d-array
860                IF ( av == 0 ) THEN
861                   local_pf(:,:,nzb+1) = REAL( fill_value, KIND = wp )
862!
863!--                In case of default surfaces, clean-up flux by density.
864!--                In case of land-surfaces, convert fluxes into
865!--                dynamic units
866                   DO  m = 1, surf_def_h(0)%ns
867                      i = surf_def_h(0)%i(m)
868                      j = surf_def_h(0)%j(m)
869                      k = surf_def_h(0)%k(m)
870                      local_pf(i,j,nzb+1) = surf_def_h(0)%qsws(m) *            &
871                                            waterflux_output_conversion(k)
872                   ENDDO
873                   DO  m = 1, surf_lsm_h%ns
874                      i = surf_lsm_h%i(m)
875                      j = surf_lsm_h%j(m)
876                      k = surf_lsm_h%k(m)
877                      local_pf(i,j,nzb+1) = surf_lsm_h%qsws(m) * l_v
878                   ENDDO
879                ELSE
880                   IF ( .NOT. ALLOCATED( qsws_av ) ) THEN
881                      ALLOCATE( qsws_av(nysg:nyng,nxlg:nxrg) )
882                      qsws_av = REAL( fill_value, KIND = wp )
883                   ENDIF
884                   DO  i = nxl, nxr
885                      DO  j = nys, nyn 
886                         local_pf(i,j,nzb+1) =  qsws_av(j,i)
887                      ENDDO
888                   ENDDO
889                ENDIF
890                resorted = .TRUE.
891                two_d = .TRUE.
892                level_z(nzb+1) = zu(nzb+1)
893
894             CASE ( 'qv_xy', 'qv_xz', 'qv_yz' )
895                IF ( av == 0 )  THEN
896                   DO  i = nxl, nxr
897                      DO  j = nys, nyn
898                         DO  k = nzb, nzt+1
899                            local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
900                         ENDDO
901                      ENDDO
902                   ENDDO
903                   resorted = .TRUE.
904                ELSE
905                   IF ( .NOT. ALLOCATED( qv_av ) ) THEN
906                      ALLOCATE( qv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
907                      qv_av = REAL( fill_value, KIND = wp )
908                   ENDIF
909                   to_be_resorted => qv_av
910                ENDIF
911                IF ( mode == 'xy' )  level_z = zu
912
913             CASE ( 'r_a*_xy' )        ! 2d-array
914                IF ( av == 0 )  THEN
915                   DO  m = 1, surf_lsm_h%ns
916                      i                   = surf_lsm_h%i(m)           
917                      j                   = surf_lsm_h%j(m)
918                      local_pf(i,j,nzb+1) = surf_lsm_h%r_a(m)
919                   ENDDO
920
921                   DO  m = 1, surf_usm_h%ns
922                      i   = surf_usm_h%i(m)           
923                      j   = surf_usm_h%j(m)
924                      local_pf(i,j,nzb+1) =                                          &
925                                 ( surf_usm_h%frac(ind_veg_wall,m)  *                &
926                                   surf_usm_h%r_a(m)       +                         & 
927                                   surf_usm_h%frac(ind_pav_green,m) *                &
928                                   surf_usm_h%r_a_green(m) +                         & 
929                                   surf_usm_h%frac(ind_wat_win,m)   *                &
930                                   surf_usm_h%r_a_window(m) )
931                   ENDDO
932                ELSE
933                   IF ( .NOT. ALLOCATED( r_a_av ) ) THEN
934                      ALLOCATE( r_a_av(nysg:nyng,nxlg:nxrg) )
935                      r_a_av = REAL( fill_value, KIND = wp )
936                   ENDIF
937                   DO  i = nxl, nxr
938                      DO  j = nys, nyn
939                         local_pf(i,j,nzb+1) = r_a_av(j,i)
940                      ENDDO
941                   ENDDO
942                ENDIF
943                resorted       = .TRUE.
944                two_d          = .TRUE.
945                level_z(nzb+1) = zu(nzb+1)
946
947             CASE ( 's_xy', 's_xz', 's_yz' )
948                IF ( av == 0 )  THEN
949                   to_be_resorted => s
950                ELSE
951                   IF ( .NOT. ALLOCATED( s_av ) ) THEN
952                      ALLOCATE( s_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
953                      s_av = REAL( fill_value, KIND = wp )
954                   ENDIF
955                   to_be_resorted => s_av
956                ENDIF
957
958             CASE ( 'shf*_xy' )        ! 2d-array
959                IF ( av == 0 ) THEN
960!
961!--                In case of default surfaces, clean-up flux by density.
962!--                In case of land- and urban-surfaces, convert fluxes into
963!--                dynamic units.
964                   DO  m = 1, surf_def_h(0)%ns
965                      i = surf_def_h(0)%i(m)
966                      j = surf_def_h(0)%j(m)
967                      k = surf_def_h(0)%k(m)
968                      local_pf(i,j,nzb+1) = surf_def_h(0)%shf(m) *             &
969                                            heatflux_output_conversion(k)
970                   ENDDO
971                   DO  m = 1, surf_lsm_h%ns
972                      i = surf_lsm_h%i(m)
973                      j = surf_lsm_h%j(m)
974                      k = surf_lsm_h%k(m)
975                      local_pf(i,j,nzb+1) = surf_lsm_h%shf(m) * c_p
976                   ENDDO
977                   DO  m = 1, surf_usm_h%ns
978                      i = surf_usm_h%i(m)
979                      j = surf_usm_h%j(m)
980                      k = surf_usm_h%k(m)
981                      local_pf(i,j,nzb+1) = surf_usm_h%shf(m) * c_p
982                   ENDDO
983                ELSE
984                   IF ( .NOT. ALLOCATED( shf_av ) ) THEN
985                      ALLOCATE( shf_av(nysg:nyng,nxlg:nxrg) )
986                      shf_av = REAL( fill_value, KIND = wp )
987                   ENDIF
988                   DO  i = nxl, nxr
989                      DO  j = nys, nyn
990                         local_pf(i,j,nzb+1) =  shf_av(j,i)
991                      ENDDO
992                   ENDDO
993                ENDIF
994                resorted = .TRUE.
995                two_d = .TRUE.
996                level_z(nzb+1) = zu(nzb+1)
997               
998             CASE ( 'ssws*_xy' )        ! 2d-array
999                IF ( av == 0 ) THEN
1000                   DO  m = 1, surf_def_h(0)%ns
1001                      i = surf_def_h(0)%i(m)
1002                      j = surf_def_h(0)%j(m)
1003                      local_pf(i,j,nzb+1) = surf_def_h(0)%ssws(m)
1004                   ENDDO
1005                   DO  m = 1, surf_lsm_h%ns
1006                      i = surf_lsm_h%i(m)
1007                      j = surf_lsm_h%j(m)
1008                      local_pf(i,j,nzb+1) = surf_lsm_h%ssws(m)
1009                   ENDDO
1010                   DO  m = 1, surf_usm_h%ns
1011                      i = surf_usm_h%i(m)
1012                      j = surf_usm_h%j(m)
1013                      local_pf(i,j,nzb+1) = surf_usm_h%ssws(m)
1014                   ENDDO
1015                ELSE
1016                   IF ( .NOT. ALLOCATED( ssws_av ) ) THEN
1017                      ALLOCATE( ssws_av(nysg:nyng,nxlg:nxrg) )
1018                      ssws_av = REAL( fill_value, KIND = wp )
1019                   ENDIF
1020                   DO  i = nxl, nxr
1021                      DO  j = nys, nyn 
1022                         local_pf(i,j,nzb+1) =  ssws_av(j,i)
1023                      ENDDO
1024                   ENDDO
1025                ENDIF
1026                resorted = .TRUE.
1027                two_d = .TRUE.
1028                level_z(nzb+1) = zu(nzb+1)               
1029
1030             CASE ( 't*_xy' )        ! 2d-array
1031                IF ( av == 0 )  THEN
1032                   DO  m = 1, surf_def_h(0)%ns
1033                      i = surf_def_h(0)%i(m)
1034                      j = surf_def_h(0)%j(m)
1035                      local_pf(i,j,nzb+1) = surf_def_h(0)%ts(m)
1036                   ENDDO
1037                   DO  m = 1, surf_lsm_h%ns
1038                      i = surf_lsm_h%i(m)
1039                      j = surf_lsm_h%j(m)
1040                      local_pf(i,j,nzb+1) = surf_lsm_h%ts(m)
1041                   ENDDO
1042                   DO  m = 1, surf_usm_h%ns
1043                      i = surf_usm_h%i(m)
1044                      j = surf_usm_h%j(m)
1045                      local_pf(i,j,nzb+1) = surf_usm_h%ts(m)
1046                   ENDDO
1047                ELSE
1048                   IF ( .NOT. ALLOCATED( ts_av ) ) THEN
1049                      ALLOCATE( ts_av(nysg:nyng,nxlg:nxrg) )
1050                      ts_av = REAL( fill_value, KIND = wp )
1051                   ENDIF
1052                   DO  i = nxl, nxr
1053                      DO  j = nys, nyn
1054                         local_pf(i,j,nzb+1) = ts_av(j,i)
1055                      ENDDO
1056                   ENDDO
1057                ENDIF
1058                resorted = .TRUE.
1059                two_d = .TRUE.
1060                level_z(nzb+1) = zu(nzb+1)
1061
1062             CASE ( 'tsurf*_xy' )        ! 2d-array
1063                IF ( av == 0 )  THEN
1064                   DO  m = 1, surf_def_h(0)%ns
1065                      i                   = surf_def_h(0)%i(m)           
1066                      j                   = surf_def_h(0)%j(m)
1067                      local_pf(i,j,nzb+1) = surf_def_h(0)%pt_surface(m)
1068                   ENDDO
1069
1070                   DO  m = 1, surf_lsm_h%ns
1071                      i                   = surf_lsm_h%i(m)           
1072                      j                   = surf_lsm_h%j(m)
1073                      local_pf(i,j,nzb+1) = surf_lsm_h%pt_surface(m)
1074                   ENDDO
1075
1076                   DO  m = 1, surf_usm_h%ns
1077                      i   = surf_usm_h%i(m)           
1078                      j   = surf_usm_h%j(m)
1079                      local_pf(i,j,nzb+1) = surf_usm_h%pt_surface(m)
1080                   ENDDO
1081
1082                ELSE
1083                   IF ( .NOT. ALLOCATED( tsurf_av ) ) THEN
1084                      ALLOCATE( tsurf_av(nysg:nyng,nxlg:nxrg) )
1085                      tsurf_av = REAL( fill_value, KIND = wp )
1086                   ENDIF
1087                   DO  i = nxl, nxr
1088                      DO  j = nys, nyn
1089                         local_pf(i,j,nzb+1) = tsurf_av(j,i)
1090                      ENDDO
1091                   ENDDO
1092                ENDIF
1093                resorted       = .TRUE.
1094                two_d          = .TRUE.
1095                level_z(nzb+1) = zu(nzb+1)
1096
1097             CASE ( 'u_xy', 'u_xz', 'u_yz' )
1098                flag_nr = 1
1099                IF ( av == 0 )  THEN
1100                   to_be_resorted => u
1101                ELSE
1102                   IF ( .NOT. ALLOCATED( u_av ) ) THEN
1103                      ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1104                      u_av = REAL( fill_value, KIND = wp )
1105                   ENDIF
1106                   to_be_resorted => u_av
1107                ENDIF
1108                IF ( mode == 'xy' )  level_z = zu
1109!
1110!--             Substitute the values generated by "mirror" boundary condition
1111!--             at the bottom boundary by the real surface values.
1112                IF ( do2d(av,if) == 'u_xz'  .OR.  do2d(av,if) == 'u_yz' )  THEN
1113                   IF ( ibc_uv_b == 0 )  local_pf(:,:,nzb) = 0.0_wp
1114                ENDIF
1115
1116             CASE ( 'u*_xy' )        ! 2d-array
1117                IF ( av == 0 )  THEN
1118                   DO  m = 1, surf_def_h(0)%ns
1119                      i = surf_def_h(0)%i(m)
1120                      j = surf_def_h(0)%j(m)
1121                      local_pf(i,j,nzb+1) = surf_def_h(0)%us(m)
1122                   ENDDO
1123                   DO  m = 1, surf_lsm_h%ns
1124                      i = surf_lsm_h%i(m)
1125                      j = surf_lsm_h%j(m)
1126                      local_pf(i,j,nzb+1) = surf_lsm_h%us(m)
1127                   ENDDO
1128                   DO  m = 1, surf_usm_h%ns
1129                      i = surf_usm_h%i(m)
1130                      j = surf_usm_h%j(m)
1131                      local_pf(i,j,nzb+1) = surf_usm_h%us(m)
1132                   ENDDO
1133                ELSE
1134                   IF ( .NOT. ALLOCATED( us_av ) ) THEN
1135                      ALLOCATE( us_av(nysg:nyng,nxlg:nxrg) )
1136                      us_av = REAL( fill_value, KIND = wp )
1137                   ENDIF
1138                   DO  i = nxl, nxr
1139                      DO  j = nys, nyn
1140                         local_pf(i,j,nzb+1) = us_av(j,i)
1141                      ENDDO
1142                   ENDDO
1143                ENDIF
1144                resorted = .TRUE.
1145                two_d = .TRUE.
1146                level_z(nzb+1) = zu(nzb+1)
1147
1148             CASE ( 'v_xy', 'v_xz', 'v_yz' )
1149                flag_nr = 2
1150                IF ( av == 0 )  THEN
1151                   to_be_resorted => v
1152                ELSE
1153                   IF ( .NOT. ALLOCATED( v_av ) ) THEN
1154                      ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1155                      v_av = REAL( fill_value, KIND = wp )
1156                   ENDIF
1157                   to_be_resorted => v_av
1158                ENDIF
1159                IF ( mode == 'xy' )  level_z = zu
1160!
1161!--             Substitute the values generated by "mirror" boundary condition
1162!--             at the bottom boundary by the real surface values.
1163                IF ( do2d(av,if) == 'v_xz'  .OR.  do2d(av,if) == 'v_yz' )  THEN
1164                   IF ( ibc_uv_b == 0 )  local_pf(:,:,nzb) = 0.0_wp
1165                ENDIF
1166
1167             CASE ( 'vpt_xy', 'vpt_xz', 'vpt_yz' )
1168                IF ( av == 0 )  THEN
1169                   to_be_resorted => vpt
1170                ELSE
1171                   IF ( .NOT. ALLOCATED( vpt_av ) ) THEN
1172                      ALLOCATE( vpt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1173                      vpt_av = REAL( fill_value, KIND = wp )
1174                   ENDIF
1175                   to_be_resorted => vpt_av
1176                ENDIF
1177                IF ( mode == 'xy' )  level_z = zu
1178
1179             CASE ( 'w_xy', 'w_xz', 'w_yz' )
1180                flag_nr = 3
1181                IF ( av == 0 )  THEN
1182                   to_be_resorted => w
1183                ELSE
1184                   IF ( .NOT. ALLOCATED( w_av ) ) THEN
1185                      ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1186                      w_av = REAL( fill_value, KIND = wp )
1187                   ENDIF
1188                   to_be_resorted => w_av
1189                ENDIF
1190                IF ( mode == 'xy' )  level_z = zw
1191
1192             CASE ( 'z0*_xy' )        ! 2d-array
1193                IF ( av == 0 ) THEN
1194                   DO  m = 1, surf_def_h(0)%ns
1195                      i = surf_def_h(0)%i(m)
1196                      j = surf_def_h(0)%j(m)
1197                      local_pf(i,j,nzb+1) = surf_def_h(0)%z0(m)
1198                   ENDDO
1199                   DO  m = 1, surf_lsm_h%ns
1200                      i = surf_lsm_h%i(m)
1201                      j = surf_lsm_h%j(m)
1202                      local_pf(i,j,nzb+1) = surf_lsm_h%z0(m)
1203                   ENDDO
1204                   DO  m = 1, surf_usm_h%ns
1205                      i = surf_usm_h%i(m)
1206                      j = surf_usm_h%j(m)
1207                      local_pf(i,j,nzb+1) = surf_usm_h%z0(m)
1208                   ENDDO
1209                ELSE
1210                   IF ( .NOT. ALLOCATED( z0_av ) ) THEN
1211                      ALLOCATE( z0_av(nysg:nyng,nxlg:nxrg) )
1212                      z0_av = REAL( fill_value, KIND = wp )
1213                   ENDIF
1214                   DO  i = nxl, nxr
1215                      DO  j = nys, nyn
1216                         local_pf(i,j,nzb+1) =  z0_av(j,i)
1217                      ENDDO
1218                   ENDDO
1219                ENDIF
1220                resorted = .TRUE.
1221                two_d = .TRUE.
1222                level_z(nzb+1) = zu(nzb+1)
1223
1224             CASE ( 'z0h*_xy' )        ! 2d-array
1225                IF ( av == 0 ) THEN
1226                   DO  m = 1, surf_def_h(0)%ns
1227                      i = surf_def_h(0)%i(m)
1228                      j = surf_def_h(0)%j(m)
1229                      local_pf(i,j,nzb+1) = surf_def_h(0)%z0h(m)
1230                   ENDDO
1231                   DO  m = 1, surf_lsm_h%ns
1232                      i = surf_lsm_h%i(m)
1233                      j = surf_lsm_h%j(m)
1234                      local_pf(i,j,nzb+1) = surf_lsm_h%z0h(m)
1235                   ENDDO
1236                   DO  m = 1, surf_usm_h%ns
1237                      i = surf_usm_h%i(m)
1238                      j = surf_usm_h%j(m)
1239                      local_pf(i,j,nzb+1) = surf_usm_h%z0h(m)
1240                   ENDDO
1241                ELSE
1242                   IF ( .NOT. ALLOCATED( z0h_av ) ) THEN
1243                      ALLOCATE( z0h_av(nysg:nyng,nxlg:nxrg) )
1244                      z0h_av = REAL( fill_value, KIND = wp )
1245                   ENDIF
1246                   DO  i = nxl, nxr
1247                      DO  j = nys, nyn
1248                         local_pf(i,j,nzb+1) =  z0h_av(j,i)
1249                      ENDDO
1250                   ENDDO
1251                ENDIF
1252                resorted = .TRUE.
1253                two_d = .TRUE.
1254                level_z(nzb+1) = zu(nzb+1)
1255
1256             CASE ( 'z0q*_xy' )        ! 2d-array
1257                IF ( av == 0 ) THEN
1258                   DO  m = 1, surf_def_h(0)%ns
1259                      i = surf_def_h(0)%i(m)
1260                      j = surf_def_h(0)%j(m)
1261                      local_pf(i,j,nzb+1) = surf_def_h(0)%z0q(m)
1262                   ENDDO
1263                   DO  m = 1, surf_lsm_h%ns
1264                      i = surf_lsm_h%i(m)
1265                      j = surf_lsm_h%j(m)
1266                      local_pf(i,j,nzb+1) = surf_lsm_h%z0q(m)
1267                   ENDDO
1268                   DO  m = 1, surf_usm_h%ns
1269                      i = surf_usm_h%i(m)
1270                      j = surf_usm_h%j(m)
1271                      local_pf(i,j,nzb+1) = surf_usm_h%z0q(m)
1272                   ENDDO
1273                ELSE
1274                   IF ( .NOT. ALLOCATED( z0q_av ) ) THEN
1275                      ALLOCATE( z0q_av(nysg:nyng,nxlg:nxrg) )
1276                      z0q_av = REAL( fill_value, KIND = wp )
1277                   ENDIF
1278                   DO  i = nxl, nxr
1279                      DO  j = nys, nyn
1280                         local_pf(i,j,nzb+1) =  z0q_av(j,i)
1281                      ENDDO
1282                   ENDDO
1283                ENDIF
1284                resorted = .TRUE.
1285                two_d = .TRUE.
1286                level_z(nzb+1) = zu(nzb+1)
1287
1288             CASE DEFAULT
1289
1290!
1291!--             Quantities of other modules
1292                IF ( .NOT. found  .AND.  bulk_cloud_model )  THEN
1293                   CALL bcm_data_output_2d( av, do2d(av,if), found, grid, mode,&
1294                                            local_pf, two_d, nzb_do, nzt_do )
1295                ENDIF
1296
1297                IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
1298                   CALL gust_data_output_2d( av, do2d(av,if), found, grid,     &
1299                                             local_pf, two_d, nzb_do, nzt_do )
1300                ENDIF
1301
1302                IF ( .NOT. found  .AND.  land_surface )  THEN
1303                   CALL lsm_data_output_2d( av, do2d(av,if), found, grid, mode,&
1304                                            local_pf, two_d, nzb_do, nzt_do )
1305                ENDIF
1306
1307                IF ( .NOT. found  .AND.  ocean_mode )  THEN
1308                   CALL ocean_data_output_2d( av, do2d(av,if), found, grid,    &
1309                                              mode, local_pf, nzb_do, nzt_do )
1310                ENDIF
1311
1312                IF ( .NOT. found  .AND.  radiation )  THEN
1313                   CALL radiation_data_output_2d( av, do2d(av,if), found, grid,&
1314                                                  mode, local_pf, two_d,       &
1315                                                  nzb_do, nzt_do  )
1316                ENDIF
1317
1318                IF ( .NOT. found )  THEN
1319                   CALL tcm_data_output_2d( av, do2d(av,if), found, grid, mode,&
1320                                             local_pf, nzb_do, nzt_do )
1321                ENDIF
1322
1323                IF ( .NOT. found  .AND.  uv_exposure )  THEN
1324                   CALL uvem_data_output_2d( av, do2d(av,if), found, grid,     &
1325                                             local_pf, two_d, nzb_do, nzt_do )
1326                ENDIF
1327
1328!
1329!--             User defined quantities
1330                IF ( .NOT. found )  THEN
1331                   CALL user_data_output_2d( av, do2d(av,if), found, grid,     &
1332                                             local_pf, two_d, nzb_do, nzt_do )
1333                ENDIF
1334
1335                resorted = .TRUE.
1336
1337                IF ( grid == 'zu' )  THEN
1338                   IF ( mode == 'xy' )  level_z = zu
1339                ELSEIF ( grid == 'zw' )  THEN
1340                   IF ( mode == 'xy' )  level_z = zw
1341                ELSEIF ( grid == 'zu1' ) THEN
1342                   IF ( mode == 'xy' )  level_z(nzb+1) = zu(nzb+1)
1343                ELSEIF ( grid == 'zs' ) THEN
1344                   IF ( mode == 'xy' )  level_z = zs
1345                ENDIF
1346
1347                IF ( .NOT. found )  THEN
1348                   message_string = 'no output provided for: ' //              &
1349                                    TRIM( do2d(av,if) )
1350                   CALL message( 'data_output_2d', 'PA0181', 0, 0, 0, 6, 0 )
1351                ENDIF
1352
1353          END SELECT
1354
1355!
1356!--       Resort the array to be output, if not done above. Flag topography
1357!--       grid points with fill values, using the corresponding maksing flag.
1358          IF ( .NOT. resorted )  THEN
1359             DO  i = nxl, nxr
1360                DO  j = nys, nyn
1361                   DO  k = nzb_do, nzt_do
1362                      local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),          &
1363                                               REAL( fill_value, KIND = wp ),  &
1364                                               BTEST( wall_flags_0(k,j,i),     &
1365                                                      flag_nr ) ) 
1366                   ENDDO
1367                ENDDO
1368             ENDDO
1369          ENDIF
1370
1371!
1372!--       Output of the individual cross-sections, depending on the cross-
1373!--       section mode chosen.
1374          is = 1
1375   loop1: DO WHILE ( section(is,s_ind) /= -9999  .OR.  two_d )
1376
1377             SELECT CASE ( mode )
1378
1379                CASE ( 'xy' )
1380!
1381!--                Determine the cross section index
1382                   IF ( two_d )  THEN
1383                      layer_xy = nzb+1
1384                   ELSE
1385                      layer_xy = section(is,s_ind)
1386                   ENDIF
1387
1388!
1389!--                Exit the loop for layers beyond the data output domain
1390!--                (used for soil model)
1391                   IF ( layer_xy > nzt_do )  THEN
1392                      EXIT loop1
1393                   ENDIF
1394
1395!
1396!--                Update the netCDF xy cross section time axis.
1397!--                In case of parallel output, this is only done by PE0
1398!--                to increase the performance.
1399                   IF ( simulated_time /= do2d_xy_last_time(av) )  THEN
1400                      do2d_xy_time_count(av) = do2d_xy_time_count(av) + 1
1401                      do2d_xy_last_time(av)  = simulated_time
1402                      IF ( myid == 0 )  THEN
1403                         IF ( .NOT. data_output_2d_on_each_pe  &
1404                              .OR.  netcdf_data_format > 4 )   &
1405                         THEN
1406#if defined( __netcdf )
1407                            nc_stat = NF90_PUT_VAR( id_set_xy(av),             &
1408                                                    id_var_time_xy(av),        &
1409                                             (/ time_since_reference_point /), &
1410                                         start = (/ do2d_xy_time_count(av) /), &
1411                                                    count = (/ 1 /) )
1412                            CALL netcdf_handle_error( 'data_output_2d', 53 )
1413#endif
1414                         ENDIF
1415                      ENDIF
1416                   ENDIF
1417!
1418!--                If required, carry out averaging along z
1419                   IF ( section(is,s_ind) == -1  .AND.  .NOT. two_d )  THEN
1420
1421                      local_2d = 0.0_wp
1422!
1423!--                   Carry out the averaging (all data are on the PE)
1424                      DO  k = nzb_do, nzt_do
1425                         DO  j = nys, nyn
1426                            DO  i = nxl, nxr
1427                               local_2d(i,j) = local_2d(i,j) + local_pf(i,j,k)
1428                            ENDDO
1429                         ENDDO
1430                      ENDDO
1431
1432                      local_2d = local_2d / ( nzt_do - nzb_do + 1.0_wp)
1433
1434                   ELSE
1435!
1436!--                   Just store the respective section on the local array
1437                      local_2d = local_pf(:,:,layer_xy)
1438
1439                   ENDIF
1440
1441#if defined( __parallel )
1442                   IF ( netcdf_data_format > 4 )  THEN
1443!
1444!--                   Parallel output in netCDF4/HDF5 format.
1445                      IF ( two_d ) THEN
1446                         iis = 1
1447                      ELSE
1448                         iis = is
1449                      ENDIF
1450
1451#if defined( __netcdf )
1452!
1453!--                   For parallel output, all cross sections are first stored
1454!--                   here on a local array and will be written to the output
1455!--                   file afterwards to increase the performance.
1456                      DO  i = nxl, nxr
1457                         DO  j = nys, nyn
1458                            local_2d_sections(i,j,iis) = local_2d(i,j)
1459                         ENDDO
1460                      ENDDO
1461#endif
1462                   ELSE
1463
1464                      IF ( data_output_2d_on_each_pe )  THEN
1465!
1466!--                      Output of partial arrays on each PE
1467#if defined( __netcdf )
1468                         IF ( myid == 0 )  THEN
1469                            WRITE ( 21 )  time_since_reference_point,          &
1470                                          do2d_xy_time_count(av), av
1471                         ENDIF
1472#endif
1473                         DO  i = 0, io_blocks-1
1474                            IF ( i == io_group )  THEN
1475                               WRITE ( 21 )  nxl, nxr, nys, nyn, nys, nyn
1476                               WRITE ( 21 )  local_2d
1477                            ENDIF
1478#if defined( __parallel )
1479                            CALL MPI_BARRIER( comm2d, ierr )
1480#endif
1481                         ENDDO
1482
1483                      ELSE
1484!
1485!--                      PE0 receives partial arrays from all processors and
1486!--                      then outputs them. Here a barrier has to be set,
1487!--                      because otherwise "-MPI- FATAL: Remote protocol queue
1488!--                      full" may occur.
1489                         CALL MPI_BARRIER( comm2d, ierr )
1490
1491                         ngp = ( nxr-nxl+1 ) * ( nyn-nys+1 )
1492                         IF ( myid == 0 )  THEN
1493!
1494!--                         Local array can be relocated directly.
1495                            total_2d(nxl:nxr,nys:nyn) = local_2d
1496!
1497!--                         Receive data from all other PEs.
1498                            DO  n = 1, numprocs-1
1499!
1500!--                            Receive index limits first, then array.
1501!--                            Index limits are received in arbitrary order from
1502!--                            the PEs.
1503                               CALL MPI_RECV( ind(1), 4, MPI_INTEGER,          &
1504                                              MPI_ANY_SOURCE, 0, comm2d,       &
1505                                              status, ierr )
1506                               sender = status(MPI_SOURCE)
1507                               DEALLOCATE( local_2d )
1508                               ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
1509                               CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp,    &
1510                                              MPI_REAL, sender, 1, comm2d,     &
1511                                              status, ierr )
1512                               total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
1513                            ENDDO
1514!
1515!--                         Relocate the local array for the next loop increment
1516                            DEALLOCATE( local_2d )
1517                            ALLOCATE( local_2d(nxl:nxr,nys:nyn) )
1518
1519#if defined( __netcdf )
1520                            IF ( two_d ) THEN
1521                               nc_stat = NF90_PUT_VAR( id_set_xy(av),       &
1522                                                       id_var_do2d(av,if),  &
1523                                                       total_2d(0:nx,0:ny), &
1524                             start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
1525                                             count = (/ nx+1, ny+1, 1, 1 /) )
1526                            ELSE
1527                               nc_stat = NF90_PUT_VAR( id_set_xy(av),       &
1528                                                       id_var_do2d(av,if),  &
1529                                                       total_2d(0:nx,0:ny), &
1530                            start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
1531                                             count = (/ nx+1, ny+1, 1, 1 /) )
1532                            ENDIF
1533                            CALL netcdf_handle_error( 'data_output_2d', 54 )
1534#endif
1535
1536                         ELSE
1537!
1538!--                         First send the local index limits to PE0
1539                            ind(1) = nxl; ind(2) = nxr
1540                            ind(3) = nys; ind(4) = nyn
1541                            CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0,       &
1542                                           comm2d, ierr )
1543!
1544!--                         Send data to PE0
1545                            CALL MPI_SEND( local_2d(nxl,nys), ngp,             &
1546                                           MPI_REAL, 0, 1, comm2d, ierr )
1547                         ENDIF
1548!
1549!--                      A barrier has to be set, because otherwise some PEs may
1550!--                      proceed too fast so that PE0 may receive wrong data on
1551!--                      tag 0
1552                         CALL MPI_BARRIER( comm2d, ierr )
1553                      ENDIF
1554
1555                   ENDIF
1556#else
1557#if defined( __netcdf )
1558                   IF ( two_d ) THEN
1559                      nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
1560                                              id_var_do2d(av,if),           &
1561                                              local_2d(nxl:nxr,nys:nyn),    &
1562                             start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
1563                                           count = (/ nx+1, ny+1, 1, 1 /) )
1564                   ELSE
1565                      nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
1566                                              id_var_do2d(av,if),           &
1567                                              local_2d(nxl:nxr,nys:nyn),    &
1568                            start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
1569                                           count = (/ nx+1, ny+1, 1, 1 /) )
1570                   ENDIF
1571                   CALL netcdf_handle_error( 'data_output_2d', 447 )
1572#endif
1573#endif
1574
1575!
1576!--                For 2D-arrays (e.g. u*) only one cross-section is available.
1577!--                Hence exit loop of output levels.
1578                   IF ( two_d )  THEN
1579                      IF ( netcdf_data_format < 5 )  two_d = .FALSE.
1580                      EXIT loop1
1581                   ENDIF
1582
1583                CASE ( 'xz' )
1584!
1585!--                Update the netCDF xz cross section time axis.
1586!--                In case of parallel output, this is only done by PE0
1587!--                to increase the performance.
1588                   IF ( simulated_time /= do2d_xz_last_time(av) )  THEN
1589                      do2d_xz_time_count(av) = do2d_xz_time_count(av) + 1
1590                      do2d_xz_last_time(av)  = simulated_time
1591                      IF ( myid == 0 )  THEN
1592                         IF ( .NOT. data_output_2d_on_each_pe  &
1593                              .OR.  netcdf_data_format > 4 )   &
1594                         THEN
1595#if defined( __netcdf )
1596                            nc_stat = NF90_PUT_VAR( id_set_xz(av),             &
1597                                                    id_var_time_xz(av),        &
1598                                             (/ time_since_reference_point /), &
1599                                         start = (/ do2d_xz_time_count(av) /), &
1600                                                    count = (/ 1 /) )
1601                            CALL netcdf_handle_error( 'data_output_2d', 56 )
1602#endif
1603                         ENDIF
1604                      ENDIF
1605                   ENDIF
1606
1607!
1608!--                If required, carry out averaging along y
1609                   IF ( section(is,s_ind) == -1 )  THEN
1610
1611                      ALLOCATE( local_2d_l(nxl:nxr,nzb_do:nzt_do) )
1612                      local_2d_l = 0.0_wp
1613                      ngp = ( nxr-nxl + 1 ) * ( nzt_do-nzb_do + 1 )
1614!
1615!--                   First local averaging on the PE
1616                      DO  k = nzb_do, nzt_do
1617                         DO  j = nys, nyn
1618                            DO  i = nxl, nxr
1619                               local_2d_l(i,k) = local_2d_l(i,k) +             &
1620                                                 local_pf(i,j,k)
1621                            ENDDO
1622                         ENDDO
1623                      ENDDO
1624#if defined( __parallel )
1625!
1626!--                   Now do the averaging over all PEs along y
1627                      IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1628                      CALL MPI_ALLREDUCE( local_2d_l(nxl,nzb_do),                &
1629                                          local_2d(nxl,nzb_do), ngp, MPI_REAL,   &
1630                                          MPI_SUM, comm1dy, ierr )
1631#else
1632                      local_2d = local_2d_l
1633#endif
1634                      local_2d = local_2d / ( ny + 1.0_wp )
1635
1636                      DEALLOCATE( local_2d_l )
1637
1638                   ELSE
1639!
1640!--                   Just store the respective section on the local array
1641!--                   (but only if it is available on this PE!)
1642                      IF ( section(is,s_ind) >= nys  .AND.  section(is,s_ind) <= nyn ) &
1643                      THEN
1644                         local_2d = local_pf(:,section(is,s_ind),nzb_do:nzt_do)
1645                      ENDIF
1646
1647                   ENDIF
1648
1649#if defined( __parallel )
1650                   IF ( netcdf_data_format > 4 )  THEN
1651!
1652!--                   Output in netCDF4/HDF5 format.
1653!--                   Output only on those PEs where the respective cross
1654!--                   sections reside. Cross sections averaged along y are
1655!--                   output on the respective first PE along y (myidy=0).
1656                      IF ( ( section(is,s_ind) >= nys  .AND.                   &
1657                             section(is,s_ind) <= nyn )  .OR.                  &
1658                           ( section(is,s_ind) == -1  .AND.  myidy == 0 ) )  THEN
1659#if defined( __netcdf )
1660!
1661!--                      For parallel output, all cross sections are first
1662!--                      stored here on a local array and will be written to the
1663!--                      output file afterwards to increase the performance.
1664                         DO  i = nxl, nxr
1665                            DO  k = nzb_do, nzt_do
1666                               local_2d_sections_l(i,is,k) = local_2d(i,k)
1667                            ENDDO
1668                         ENDDO
1669#endif
1670                      ENDIF
1671
1672                   ELSE
1673
1674                      IF ( data_output_2d_on_each_pe )  THEN
1675!
1676!--                      Output of partial arrays on each PE. If the cross
1677!--                      section does not reside on the PE, output special
1678!--                      index values.
1679#if defined( __netcdf )
1680                         IF ( myid == 0 )  THEN
1681                            WRITE ( 22 )  time_since_reference_point,          &
1682                                          do2d_xz_time_count(av), av
1683                         ENDIF
1684#endif
1685                         DO  i = 0, io_blocks-1
1686                            IF ( i == io_group )  THEN
1687                               IF ( ( section(is,s_ind) >= nys  .AND.          &
1688                                      section(is,s_ind) <= nyn )  .OR.         &
1689                                    ( section(is,s_ind) == -1  .AND.           &
1690                                      nys-1 == -1 ) )                          &
1691                               THEN
1692                                  WRITE (22)  nxl, nxr, nzb_do, nzt_do, nzb, nzt+1
1693                                  WRITE (22)  local_2d
1694                               ELSE
1695                                  WRITE (22)  -1, -1, -1, -1, -1, -1
1696                               ENDIF
1697                            ENDIF
1698#if defined( __parallel )
1699                            CALL MPI_BARRIER( comm2d, ierr )
1700#endif
1701                         ENDDO
1702
1703                      ELSE
1704!
1705!--                      PE0 receives partial arrays from all processors of the
1706!--                      respective cross section and outputs them. Here a
1707!--                      barrier has to be set, because otherwise
1708!--                      "-MPI- FATAL: Remote protocol queue full" may occur.
1709                         CALL MPI_BARRIER( comm2d, ierr )
1710
1711                         ngp = ( nxr-nxl + 1 ) * ( nzt_do-nzb_do + 1 )
1712                         IF ( myid == 0 )  THEN
1713!
1714!--                         Local array can be relocated directly.
1715                            IF ( ( section(is,s_ind) >= nys  .AND.              &
1716                                   section(is,s_ind) <= nyn )  .OR.             &
1717                                 ( section(is,s_ind) == -1  .AND.               &
1718                                   nys-1 == -1 ) )  THEN
1719                               total_2d(nxl:nxr,nzb_do:nzt_do) = local_2d
1720                            ENDIF
1721!
1722!--                         Receive data from all other PEs.
1723                            DO  n = 1, numprocs-1
1724!
1725!--                            Receive index limits first, then array.
1726!--                            Index limits are received in arbitrary order from
1727!--                            the PEs.
1728                               CALL MPI_RECV( ind(1), 4, MPI_INTEGER,          &
1729                                              MPI_ANY_SOURCE, 0, comm2d,       &
1730                                              status, ierr )
1731!
1732!--                            Not all PEs have data for XZ-cross-section.
1733                               IF ( ind(1) /= -9999 )  THEN
1734                                  sender = status(MPI_SOURCE)
1735                                  DEALLOCATE( local_2d )
1736                                  ALLOCATE( local_2d(ind(1):ind(2),            &
1737                                                     ind(3):ind(4)) )
1738                                  CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1739                                                 MPI_REAL, sender, 1, comm2d,  &
1740                                                 status, ierr )
1741                                  total_2d(ind(1):ind(2),ind(3):ind(4)) =      &
1742                                                                        local_2d
1743                               ENDIF
1744                            ENDDO
1745!
1746!--                         Relocate the local array for the next loop increment
1747                            DEALLOCATE( local_2d )
1748                            ALLOCATE( local_2d(nxl:nxr,nzb_do:nzt_do) )
1749
1750#if defined( __netcdf )
1751                            nc_stat = NF90_PUT_VAR( id_set_xz(av),             &
1752                                                 id_var_do2d(av,if),           &
1753                                                 total_2d(0:nx,nzb_do:nzt_do), &
1754                               start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
1755                                          count = (/ nx+1, 1, nzt_do-nzb_do+1, 1 /) )
1756                            CALL netcdf_handle_error( 'data_output_2d', 58 )
1757#endif
1758
1759                         ELSE
1760!
1761!--                         If the cross section resides on the PE, send the
1762!--                         local index limits, otherwise send -9999 to PE0.
1763                            IF ( ( section(is,s_ind) >= nys  .AND.              &
1764                                   section(is,s_ind) <= nyn )  .OR.             &
1765                                 ( section(is,s_ind) == -1  .AND.  nys-1 == -1 ) ) &
1766                            THEN
1767                               ind(1) = nxl; ind(2) = nxr
1768                               ind(3) = nzb_do;   ind(4) = nzt_do
1769                            ELSE
1770                               ind(1) = -9999; ind(2) = -9999
1771                               ind(3) = -9999; ind(4) = -9999
1772                            ENDIF
1773                            CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0,       &
1774                                           comm2d, ierr )
1775!
1776!--                         If applicable, send data to PE0.
1777                            IF ( ind(1) /= -9999 )  THEN
1778                               CALL MPI_SEND( local_2d(nxl,nzb_do), ngp,         &
1779                                              MPI_REAL, 0, 1, comm2d, ierr )
1780                            ENDIF
1781                         ENDIF
1782!
1783!--                      A barrier has to be set, because otherwise some PEs may
1784!--                      proceed too fast so that PE0 may receive wrong data on
1785!--                      tag 0
1786                         CALL MPI_BARRIER( comm2d, ierr )
1787                      ENDIF
1788
1789                   ENDIF
1790#else
1791#if defined( __netcdf )
1792                   nc_stat = NF90_PUT_VAR( id_set_xz(av),                   &
1793                                           id_var_do2d(av,if),              &
1794                                           local_2d(nxl:nxr,nzb_do:nzt_do), &
1795                            start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
1796                                       count = (/ nx+1, 1, nzt_do-nzb_do+1, 1 /) )
1797                   CALL netcdf_handle_error( 'data_output_2d', 451 )
1798#endif
1799#endif
1800
1801                CASE ( 'yz' )
1802!
1803!--                Update the netCDF yz cross section time axis.
1804!--                In case of parallel output, this is only done by PE0
1805!--                to increase the performance.
1806                   IF ( simulated_time /= do2d_yz_last_time(av) )  THEN
1807                      do2d_yz_time_count(av) = do2d_yz_time_count(av) + 1
1808                      do2d_yz_last_time(av)  = simulated_time
1809                      IF ( myid == 0 )  THEN
1810                         IF ( .NOT. data_output_2d_on_each_pe  &
1811                              .OR.  netcdf_data_format > 4 )   &
1812                         THEN
1813#if defined( __netcdf )
1814                            nc_stat = NF90_PUT_VAR( id_set_yz(av),             &
1815                                                    id_var_time_yz(av),        &
1816                                             (/ time_since_reference_point /), &
1817                                         start = (/ do2d_yz_time_count(av) /), &
1818                                                    count = (/ 1 /) )
1819                            CALL netcdf_handle_error( 'data_output_2d', 59 )
1820#endif
1821                         ENDIF
1822                      ENDIF
1823                   ENDIF
1824
1825!
1826!--                If required, carry out averaging along x
1827                   IF ( section(is,s_ind) == -1 )  THEN
1828
1829                      ALLOCATE( local_2d_l(nys:nyn,nzb_do:nzt_do) )
1830                      local_2d_l = 0.0_wp
1831                      ngp = ( nyn-nys+1 ) * ( nzt_do-nzb_do+1 )
1832!
1833!--                   First local averaging on the PE
1834                      DO  k = nzb_do, nzt_do
1835                         DO  j = nys, nyn
1836                            DO  i = nxl, nxr
1837                               local_2d_l(j,k) = local_2d_l(j,k) +             &
1838                                                 local_pf(i,j,k)
1839                            ENDDO
1840                         ENDDO
1841                      ENDDO
1842#if defined( __parallel )
1843!
1844!--                   Now do the averaging over all PEs along x
1845                      IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1846                      CALL MPI_ALLREDUCE( local_2d_l(nys,nzb_do),                &
1847                                          local_2d(nys,nzb_do), ngp, MPI_REAL,   &
1848                                          MPI_SUM, comm1dx, ierr )
1849#else
1850                      local_2d = local_2d_l
1851#endif
1852                      local_2d = local_2d / ( nx + 1.0_wp )
1853
1854                      DEALLOCATE( local_2d_l )
1855
1856                   ELSE
1857!
1858!--                   Just store the respective section on the local array
1859!--                   (but only if it is available on this PE!)
1860                      IF ( section(is,s_ind) >= nxl  .AND.  section(is,s_ind) <= nxr ) &
1861                      THEN
1862                         local_2d = local_pf(section(is,s_ind),:,nzb_do:nzt_do)
1863                      ENDIF
1864
1865                   ENDIF
1866
1867#if defined( __parallel )
1868                   IF ( netcdf_data_format > 4 )  THEN
1869!
1870!--                   Output in netCDF4/HDF5 format.
1871!--                   Output only on those PEs where the respective cross
1872!--                   sections reside. Cross sections averaged along x are
1873!--                   output on the respective first PE along x (myidx=0).
1874                      IF ( ( section(is,s_ind) >= nxl  .AND.                       &
1875                             section(is,s_ind) <= nxr )  .OR.                      &
1876                           ( section(is,s_ind) == -1  .AND.  myidx == 0 ) )  THEN
1877#if defined( __netcdf )
1878!
1879!--                      For parallel output, all cross sections are first
1880!--                      stored here on a local array and will be written to the
1881!--                      output file afterwards to increase the performance.
1882                         DO  j = nys, nyn
1883                            DO  k = nzb_do, nzt_do
1884                               local_2d_sections_l(is,j,k) = local_2d(j,k)
1885                            ENDDO
1886                         ENDDO
1887#endif
1888                      ENDIF
1889
1890                   ELSE
1891
1892                      IF ( data_output_2d_on_each_pe )  THEN
1893!
1894!--                      Output of partial arrays on each PE. If the cross
1895!--                      section does not reside on the PE, output special
1896!--                      index values.
1897#if defined( __netcdf )
1898                         IF ( myid == 0 )  THEN
1899                            WRITE ( 23 )  time_since_reference_point,          &
1900                                          do2d_yz_time_count(av), av
1901                         ENDIF
1902#endif
1903                         DO  i = 0, io_blocks-1
1904                            IF ( i == io_group )  THEN
1905                               IF ( ( section(is,s_ind) >= nxl  .AND.          &
1906                                      section(is,s_ind) <= nxr )  .OR.         &
1907                                    ( section(is,s_ind) == -1  .AND.           &
1908                                      nxl-1 == -1 ) )                          &
1909                               THEN
1910                                  WRITE (23)  nys, nyn, nzb_do, nzt_do, nzb, nzt+1
1911                                  WRITE (23)  local_2d
1912                               ELSE
1913                                  WRITE (23)  -1, -1, -1, -1, -1, -1
1914                               ENDIF
1915                            ENDIF
1916#if defined( __parallel )
1917                            CALL MPI_BARRIER( comm2d, ierr )
1918#endif
1919                         ENDDO
1920
1921                      ELSE
1922!
1923!--                      PE0 receives partial arrays from all processors of the
1924!--                      respective cross section and outputs them. Here a
1925!--                      barrier has to be set, because otherwise
1926!--                      "-MPI- FATAL: Remote protocol queue full" may occur.
1927                         CALL MPI_BARRIER( comm2d, ierr )
1928
1929                         ngp = ( nyn-nys+1 ) * ( nzt_do-nzb_do+1 )
1930                         IF ( myid == 0 )  THEN
1931!
1932!--                         Local array can be relocated directly.
1933                            IF ( ( section(is,s_ind) >= nxl  .AND.             &
1934                                   section(is,s_ind) <= nxr )   .OR.           &
1935                                 ( section(is,s_ind) == -1  .AND.  nxl-1 == -1 ) ) &
1936                            THEN
1937                               total_2d(nys:nyn,nzb_do:nzt_do) = local_2d
1938                            ENDIF
1939!
1940!--                         Receive data from all other PEs.
1941                            DO  n = 1, numprocs-1
1942!
1943!--                            Receive index limits first, then array.
1944!--                            Index limits are received in arbitrary order from
1945!--                            the PEs.
1946                               CALL MPI_RECV( ind(1), 4, MPI_INTEGER,          &
1947                                              MPI_ANY_SOURCE, 0, comm2d,       &
1948                                              status, ierr )
1949!
1950!--                            Not all PEs have data for YZ-cross-section.
1951                               IF ( ind(1) /= -9999 )  THEN
1952                                  sender = status(MPI_SOURCE)
1953                                  DEALLOCATE( local_2d )
1954                                  ALLOCATE( local_2d(ind(1):ind(2),            &
1955                                                     ind(3):ind(4)) )
1956                                  CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1957                                                 MPI_REAL, sender, 1, comm2d,  &
1958                                                 status, ierr )
1959                                  total_2d(ind(1):ind(2),ind(3):ind(4)) =      &
1960                                                                        local_2d
1961                               ENDIF
1962                            ENDDO
1963!
1964!--                         Relocate the local array for the next loop increment
1965                            DEALLOCATE( local_2d )
1966                            ALLOCATE( local_2d(nys:nyn,nzb_do:nzt_do) )
1967
1968#if defined( __netcdf )
1969                            nc_stat = NF90_PUT_VAR( id_set_yz(av),             &
1970                                                 id_var_do2d(av,if),           &
1971                                                 total_2d(0:ny,nzb_do:nzt_do), &
1972                            start = (/ is, 1, 1, do2d_yz_time_count(av) /),    &
1973                                       count = (/ 1, ny+1, nzt_do-nzb_do+1, 1 /) )
1974                            CALL netcdf_handle_error( 'data_output_2d', 61 )
1975#endif
1976
1977                         ELSE
1978!
1979!--                         If the cross section resides on the PE, send the
1980!--                         local index limits, otherwise send -9999 to PE0.
1981                            IF ( ( section(is,s_ind) >= nxl  .AND.              &
1982                                   section(is,s_ind) <= nxr )  .OR.             &
1983                                 ( section(is,s_ind) == -1  .AND.  nxl-1 == -1 ) ) &
1984                            THEN
1985                               ind(1) = nys; ind(2) = nyn
1986                               ind(3) = nzb_do;   ind(4) = nzt_do
1987                            ELSE
1988                               ind(1) = -9999; ind(2) = -9999
1989                               ind(3) = -9999; ind(4) = -9999
1990                            ENDIF
1991                            CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0,       &
1992                                           comm2d, ierr )
1993!
1994!--                         If applicable, send data to PE0.
1995                            IF ( ind(1) /= -9999 )  THEN
1996                               CALL MPI_SEND( local_2d(nys,nzb_do), ngp,         &
1997                                              MPI_REAL, 0, 1, comm2d, ierr )
1998                            ENDIF
1999                         ENDIF
2000!
2001!--                      A barrier has to be set, because otherwise some PEs may
2002!--                      proceed too fast so that PE0 may receive wrong data on
2003!--                      tag 0
2004                         CALL MPI_BARRIER( comm2d, ierr )
2005                      ENDIF
2006
2007                   ENDIF
2008#else
2009#if defined( __netcdf )
2010                   nc_stat = NF90_PUT_VAR( id_set_yz(av),                   &
2011                                           id_var_do2d(av,if),              &
2012                                           local_2d(nys:nyn,nzb_do:nzt_do), &
2013                            start = (/ is, 1, 1, do2d_xz_time_count(av) /), &
2014                                           count = (/ 1, ny+1, nzt_do-nzb_do+1, 1 /) )
2015                   CALL netcdf_handle_error( 'data_output_2d', 452 )
2016#endif
2017#endif
2018
2019             END SELECT
2020
2021             is = is + 1
2022          ENDDO loop1
2023
2024!
2025!--       For parallel output, all data were collected before on a local array
2026!--       and are written now to the netcdf file. This must be done to increase
2027!--       the performance of the parallel output.
2028#if defined( __netcdf )
2029          IF ( netcdf_data_format > 4 )  THEN
2030
2031                SELECT CASE ( mode )
2032
2033                   CASE ( 'xy' )
2034                      IF ( two_d ) THEN
2035                         nis = 1
2036                         two_d = .FALSE.
2037                      ELSE
2038                         nis = ns
2039                      ENDIF
2040!
2041!--                   Do not output redundant ghost point data except for the
2042!--                   boundaries of the total domain.
2043!                      IF ( nxr == nx  .AND.  nyn /= ny )  THEN
2044!                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
2045!                                                 id_var_do2d(av,if),           &
2046!                                                 local_2d_sections(nxl:nxr+1,  &
2047!                                                    nys:nyn,1:nis),            &
2048!                                                 start = (/ nxl+1, nys+1, 1,   &
2049!                                                    do2d_xy_time_count(av) /), &
2050!                                                 count = (/ nxr-nxl+2,         &
2051!                                                            nyn-nys+1, nis, 1  &
2052!                                                          /) )
2053!                      ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
2054!                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
2055!                                                 id_var_do2d(av,if),           &
2056!                                                 local_2d_sections(nxl:nxr,    &
2057!                                                    nys:nyn+1,1:nis),          &
2058!                                                 start = (/ nxl+1, nys+1, 1,   &
2059!                                                    do2d_xy_time_count(av) /), &
2060!                                                 count = (/ nxr-nxl+1,         &
2061!                                                            nyn-nys+2, nis, 1  &
2062!                                                          /) )
2063!                      ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
2064!                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
2065!                                                 id_var_do2d(av,if),           &
2066!                                                 local_2d_sections(nxl:nxr+1,  &
2067!                                                    nys:nyn+1,1:nis),          &
2068!                                                 start = (/ nxl+1, nys+1, 1,   &
2069!                                                    do2d_xy_time_count(av) /), &
2070!                                                 count = (/ nxr-nxl+2,         &
2071!                                                            nyn-nys+2, nis, 1  &
2072!                                                          /) )
2073!                      ELSE
2074                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
2075                                                 id_var_do2d(av,if),           &
2076                                                 local_2d_sections(nxl:nxr,    &
2077                                                    nys:nyn,1:nis),            &
2078                                                 start = (/ nxl+1, nys+1, 1,   &
2079                                                    do2d_xy_time_count(av) /), &
2080                                                 count = (/ nxr-nxl+1,         &
2081                                                            nyn-nys+1, nis, 1  &
2082                                                          /) )
2083!                      ENDIF   
2084
2085                      CALL netcdf_handle_error( 'data_output_2d', 55 )
2086
2087                   CASE ( 'xz' )
2088!
2089!--                   First, all PEs get the information of all cross-sections.
2090!--                   Then the data are written to the output file by all PEs
2091!--                   while NF90_COLLECTIVE is set in subroutine
2092!--                   define_netcdf_header. Although redundant information are
2093!--                   written to the output file in that case, the performance
2094!--                   is significantly better compared to the case where only
2095!--                   the first row of PEs in x-direction (myidx = 0) is given
2096!--                   the output while NF90_INDEPENDENT is set.
2097                      IF ( npey /= 1 )  THEN
2098                         
2099#if defined( __parallel )
2100!
2101!--                      Distribute data over all PEs along y
2102                         ngp = ( nxr-nxl+1 ) * ( nzt_do-nzb_do+1 ) * ns
2103                         IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2104                         CALL MPI_ALLREDUCE( local_2d_sections_l(nxl,1,nzb_do),  &
2105                                             local_2d_sections(nxl,1,nzb_do),    &
2106                                             ngp, MPI_REAL, MPI_SUM, comm1dy,  &
2107                                             ierr )
2108#else
2109                         local_2d_sections = local_2d_sections_l
2110#endif
2111                      ENDIF
2112!
2113!--                   Do not output redundant ghost point data except for the
2114!--                   boundaries of the total domain.
2115!                      IF ( nxr == nx )  THEN
2116!                         nc_stat = NF90_PUT_VAR( id_set_xz(av),                &
2117!                                             id_var_do2d(av,if),               &
2118!                                             local_2d_sections(nxl:nxr+1,1:ns, &
2119!                                                nzb_do:nzt_do),                &
2120!                                             start = (/ nxl+1, 1, 1,           &
2121!                                                do2d_xz_time_count(av) /),     &
2122!                                             count = (/ nxr-nxl+2, ns, nzt_do-nzb_do+1,  &
2123!                                                        1 /) )
2124!                      ELSE
2125                         nc_stat = NF90_PUT_VAR( id_set_xz(av),                &
2126                                             id_var_do2d(av,if),               &
2127                                             local_2d_sections(nxl:nxr,1:ns,   &
2128                                                nzb_do:nzt_do),                &
2129                                             start = (/ nxl+1, 1, 1,           &
2130                                                do2d_xz_time_count(av) /),     &
2131                                             count = (/ nxr-nxl+1, ns, nzt_do-nzb_do+1,  &
2132                                                1 /) )
2133!                      ENDIF
2134
2135                      CALL netcdf_handle_error( 'data_output_2d', 57 )
2136
2137                   CASE ( 'yz' )
2138!
2139!--                   First, all PEs get the information of all cross-sections.
2140!--                   Then the data are written to the output file by all PEs
2141!--                   while NF90_COLLECTIVE is set in subroutine
2142!--                   define_netcdf_header. Although redundant information are
2143!--                   written to the output file in that case, the performance
2144!--                   is significantly better compared to the case where only
2145!--                   the first row of PEs in y-direction (myidy = 0) is given
2146!--                   the output while NF90_INDEPENDENT is set.
2147                      IF ( npex /= 1 )  THEN
2148
2149#if defined( __parallel )
2150!
2151!--                      Distribute data over all PEs along x
2152                         ngp = ( nyn-nys+1 ) * ( nzt-nzb + 2 ) * ns
2153                         IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2154                         CALL MPI_ALLREDUCE( local_2d_sections_l(1,nys,nzb_do),  &
2155                                             local_2d_sections(1,nys,nzb_do),    &
2156                                             ngp, MPI_REAL, MPI_SUM, comm1dx,  &
2157                                             ierr )
2158#else
2159                         local_2d_sections = local_2d_sections_l
2160#endif
2161                      ENDIF
2162!
2163!--                   Do not output redundant ghost point data except for the
2164!--                   boundaries of the total domain.
2165!                      IF ( nyn == ny )  THEN
2166!                         nc_stat = NF90_PUT_VAR( id_set_yz(av),                &
2167!                                             id_var_do2d(av,if),               &
2168!                                             local_2d_sections(1:ns,           &
2169!                                                nys:nyn+1,nzb_do:nzt_do),      &
2170!                                             start = (/ 1, nys+1, 1,           &
2171!                                                do2d_yz_time_count(av) /),     &
2172!                                             count = (/ ns, nyn-nys+2,         &
2173!                                                        nzt_do-nzb_do+1, 1 /) )
2174!                      ELSE
2175                         nc_stat = NF90_PUT_VAR( id_set_yz(av),                &
2176                                             id_var_do2d(av,if),               &
2177                                             local_2d_sections(1:ns,nys:nyn,   &
2178                                                nzb_do:nzt_do),                &
2179                                             start = (/ 1, nys+1, 1,           &
2180                                                do2d_yz_time_count(av) /),     &
2181                                             count = (/ ns, nyn-nys+1,         &
2182                                                        nzt_do-nzb_do+1, 1 /) )
2183!                      ENDIF
2184
2185                      CALL netcdf_handle_error( 'data_output_2d', 60 )
2186
2187                   CASE DEFAULT
2188                      message_string = 'unknown cross-section: ' // TRIM( mode )
2189                      CALL message( 'data_output_2d', 'PA0180', 1, 2, 0, 6, 0 )
2190
2191                END SELECT                     
2192
2193          ENDIF
2194#endif
2195       ENDIF
2196
2197       if = if + 1
2198       l = MAX( 2, LEN_TRIM( do2d(av,if) ) )
2199       do2d_mode = do2d(av,if)(l-1:l)
2200
2201    ENDDO
2202
2203!
2204!-- Deallocate temporary arrays.
2205    IF ( ALLOCATED( level_z ) )  DEALLOCATE( level_z )
2206    IF ( netcdf_data_format > 4 )  THEN
2207       DEALLOCATE( local_pf, local_2d, local_2d_sections )
2208       IF( mode == 'xz' .OR. mode == 'yz' ) DEALLOCATE( local_2d_sections_l )
2209    ENDIF
2210#if defined( __parallel )
2211    IF ( .NOT.  data_output_2d_on_each_pe  .AND.  myid == 0 )  THEN
2212       DEALLOCATE( total_2d )
2213    ENDIF
2214#endif
2215
2216!
2217!-- Close plot output file.
2218    file_id = 20 + s_ind
2219
2220    IF ( data_output_2d_on_each_pe )  THEN
2221       DO  i = 0, io_blocks-1
2222          IF ( i == io_group )  THEN
2223             CALL close_file( file_id )
2224          ENDIF
2225#if defined( __parallel )
2226          CALL MPI_BARRIER( comm2d, ierr )
2227#endif
2228       ENDDO
2229    ELSE
2230       IF ( myid == 0 )  CALL close_file( file_id )
2231    ENDIF
2232
2233    CALL cpu_log( log_point(3), 'data_output_2d', 'stop' )
2234
2235 END SUBROUTINE data_output_2d
Note: See TracBrowser for help on using the repository browser.