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

Last change on this file since 1686 was 1683, checked in by knoop, 9 years ago

last commit documented

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