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

Last change on this file since 2232 was 2232, checked in by suehring, 7 years ago

Adjustments according new topography and surface-modelling concept implemented

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