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

Last change on this file since 4617 was 4559, checked in by raasch, 4 years ago

files re-formatted to follow the PALM coding standard

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