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

Last change on this file since 4441 was 4441, checked in by suehring, 19 months ago

Change order of dimension in surface arrays %frac, %emissivity and %albedo to allow for better vectorization in the radiation interactions; Set back turbulent length scale to 8 x grid spacing in the parametrized mode for the synthetic turbulence generator (was accidentally changed in last commit)

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