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

Last change on this file since 2284 was 2277, checked in by kanani, 8 years ago

code documentation and cleanup

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