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

Last change on this file since 4011 was 3994, checked in by suehring, 6 years ago

new module for diagnostic output quantities added + output of turbulence intensity

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