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

Last change on this file since 1916 was 1852, checked in by hoffmann, 9 years ago

last commit documented

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