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

Last change on this file since 3272 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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