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

Last change on this file since 3946 was 3943, checked in by maronga, 5 years ago

bugfixes in urban surface model; output of greenz roof transpiration added/corrected; minor formatting improvements

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