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

Last change on this file since 1745 was 1745, checked in by gronemeier, 9 years ago

Bugfix:calculation of time levels for parallel NetCDF output

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