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

Last change on this file since 4338 was 4331, checked in by suehring, 5 years ago

New diagnostic output for 10-m wind speed; Diagnostic output of 2-m potential temperature moved to diagnostic output

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