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

Last change on this file since 4467 was 4457, checked in by raasch, 5 years ago

ghost point exchange modularized, bugfix for wrong 2d-exchange

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