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

Last change on this file since 4409 was 4360, checked in by suehring, 5 years ago

Bugfix in output of time-averaged plant-canopy quanities; Output of plant-canopy data only where tall canopy is defined; land-surface model: fix wrong location strings; tests: update urban test case; all source code files: copyright update

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