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

Last change on this file since 1962 was 1961, checked in by suehring, 8 years ago

last commit documented

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