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

Last change on this file since 1734 was 1704, checked in by raasch, 9 years ago

last commit documented

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