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

Last change on this file since 1560 was 1556, checked in by maronga, 10 years ago

last commit documented

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