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

Last change on this file since 4869 was 4828, checked in by Giersch, 4 years ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

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