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

Last change on this file since 1833 was 1823, checked in by hoffmann, 8 years ago

last commit documented

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