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

Last change on this file since 3046 was 3046, checked in by Giersch, 6 years ago

Remaining error messages revised, comments extended

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