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

Last change on this file since 3014 was 3014, checked in by maronga, 3 years ago

series of bugfixes

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