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

Last change on this file since 3163 was 3052, checked in by raasch, 6 years ago

Do not open FORTRAN binary files in case of parallel netCDF I/O

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