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

Last change on this file since 4329 was 4329, checked in by motisi, 18 months ago

Renamed wall_flags_0 to wall_flags_static_0

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