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

Last change on this file since 3448 was 3448, checked in by kanani, 3 years ago

Implementation of human thermal indices (from branch biomet_p2 at r3444)

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