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

Last change on this file since 2233 was 2233, checked in by suehring, 4 years ago

last commit documented

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