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

Last change on this file since 3516 was 3467, checked in by suehring, 5 years ago

Branch salsa @3446 re-integrated into trunk

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