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

Last change on this file since 3537 was 3525, checked in by kanani, 6 years ago

Changes related to clean-up of biometeorology (by dom_dwd_user)

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