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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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