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

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

upper bounds of cross section and 3d output changed from nx+1,ny+1 to nx,ny; no output if redundant ghost layer data to NetCDF files

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