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

Last change on this file since 3049 was 3049, checked in by Giersch, 3 years ago

Revision history corrected

  • 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 3049 2018-05-29 13:52:36Z Giersch $
27! Error messages revised
28!
29! 3045 2018-05-28 07:55:41Z Giersch
30! Error messages revised
31!
32! 3014 2018-05-09 08:42:38Z maronga
33! Added nzb_do and nzt_do for some modules for 2d output
34!
35! 3004 2018-04-27 12:33:25Z Giersch
36! precipitation_rate removed, case prr*_xy removed, to_be_resorted have to point
37! to ql_vp_av and not to ql_vp, allocation checks implemented (averaged data
38! will be assigned to fill values if no allocation happened so far)   
39!
40! 2963 2018-04-12 14:47:44Z suehring
41! Introduce index for vegetation/wall, pavement/green-wall and water/window
42! surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
43!
44! 2817 2018-02-19 16:32:21Z knoop
45! Preliminary gust module interface implemented
46!
47! 2805 2018-02-14 17:00:09Z suehring
48! Consider also default-type surfaces for surface temperature output.
49!
50! 2797 2018-02-08 13:24:35Z suehring
51! Enable output of ground-heat flux also at urban surfaces.
52!
53! 2743 2018-01-12 16:03:39Z suehring
54! In case of natural- and urban-type surfaces output surfaces fluxes in W/m2.
55!
56! 2742 2018-01-12 14:59:47Z suehring
57! Enable output of surface temperature
58!
59! 2735 2018-01-11 12:01:27Z suehring
60! output of r_a moved from land-surface to consider also urban-type surfaces
61!
62! 2718 2018-01-02 08:49:38Z maronga
63! Corrected "Former revisions" section
64!
65! 2696 2017-12-14 17:12:51Z kanani
66! Change in file header (GPL part)
67! Implementation of uv exposure model (FK)
68! Implementation of turbulence_closure_mod (TG)
69! Set fill values at topography grid points or e.g. non-natural-type surface
70! in case of LSM output (MS)
71!
72! 2512 2017-10-04 08:26:59Z raasch
73! upper bounds of cross section output changed from nx+1,ny+1 to nx,ny
74! no output of ghost layer data
75!
76! 2292 2017-06-20 09:51:42Z schwenkel
77! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
78! includes two more prognostic equations for cloud drop concentration (nc) 
79! and cloud water content (qc).
80!
81! 2277 2017-06-12 10:47:51Z kanani
82! Removed unused variables do2d_xy_n, do2d_xz_n, do2d_yz_n
83!
84! 2233 2017-05-30 18:08:54Z suehring
85!
86! 2232 2017-05-30 17:47:52Z suehring
87! Adjustments to new surface concept
88!
89!
90! 2190 2017-03-21 12:16:43Z raasch
91! bugfix for r2031: string rho replaced by rho_ocean
92!
93! 2031 2016-10-21 15:11:58Z knoop
94! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
95!
96! 2000 2016-08-20 18:09:15Z knoop
97! Forced header and separation lines into 80 columns
98!
99! 1980 2016-07-29 15:51:57Z suehring
100! Bugfix, in order to steer user-defined output, setting flag found explicitly
101! to .F.
102!
103! 1976 2016-07-27 13:28:04Z maronga
104! Output of radiation quantities is now done directly in the respective module
105!
106! 1972 2016-07-26 07:52:02Z maronga
107! Output of land surface quantities is now done directly in the respective
108! module
109!
110! 1960 2016-07-12 16:34:24Z suehring
111! Scalar surface flux added
112! Rename INTEGER variable s into s_ind, as s is already assigned to scalar
113!
114! 1849 2016-04-08 11:33:18Z hoffmann
115! precipitation_amount, precipitation_rate, prr moved to arrays_3d
116!
117! 1822 2016-04-07 07:49:42Z hoffmann
118! Output of bulk cloud physics simplified.
119!
120! 1788 2016-03-10 11:01:04Z maronga
121! Added output of z0q
122!
123! 1783 2016-03-06 18:36:17Z raasch
124! name change of netcdf routines and module + related changes
125!
126! 1745 2016-02-05 13:06:51Z gronemeier
127! Bugfix: test if time axis limit exceeds moved to point after call of check_open
128!
129! 1703 2015-11-02 12:38:44Z raasch
130! bugfix for output of single (*) xy-sections in case of parallel netcdf I/O
131!
132! 1701 2015-11-02 07:43:04Z maronga
133! Bugfix in output of RRTGM data
134!
135! 1691 2015-10-26 16:17:44Z maronga
136! Added output of Obukhov length (ol) and radiative heating rates  for RRTMG.
137! Formatting corrections.
138!
139! 1682 2015-10-07 23:56:08Z knoop
140! Code annotations made doxygen readable
141!
142! 1585 2015-04-30 07:05:52Z maronga
143! Added support for RRTMG
144!
145! 1555 2015-03-04 17:44:27Z maronga
146! Added output of r_a and r_s
147!
148! 1551 2015-03-03 14:18:16Z maronga
149! Added suppport for land surface model and radiation model output. In the course
150! of this action, the limits for vertical loops have been changed (from nzb and
151! nzt+1 to nzb_do and nzt_do, respectively in order to allow soil model output).
152! Moreover, a new vertical grid zs was introduced.
153!
154! 1359 2014-04-11 17:15:14Z hoffmann
155! New particle structure integrated.
156!
157! 1353 2014-04-08 15:21:23Z heinze
158! REAL constants provided with KIND-attribute
159!
160! 1327 2014-03-21 11:00:16Z raasch
161! parts concerning iso2d output removed,
162! -netcdf output queries
163!
164! 1320 2014-03-20 08:40:49Z raasch
165! ONLY-attribute added to USE-statements,
166! kind-parameters added to all INTEGER and REAL declaration statements,
167! kinds are defined in new module kinds,
168! revision history before 2012 removed,
169! comment fields (!:) to be used for variable explanations added to
170! all variable declaration statements
171!
172! 1318 2014-03-17 13:35:16Z raasch
173! barrier argument removed from cpu_log.
174! module interfaces removed
175!
176! 1311 2014-03-14 12:13:39Z heinze
177! bugfix: close #if defined( __netcdf )
178!
179! 1308 2014-03-13 14:58:42Z fricke
180! +local_2d_sections, local_2d_sections_l, ns
181! Check, if the limit of the time dimension is exceeded for parallel output
182! To increase the performance for parallel output, the following is done:
183! - Update of time axis is only done by PE0
184! - Cross sections are first stored on a local array and are written
185!   collectively to the output file by all PEs.
186!
187! 1115 2013-03-26 18:16:16Z hoffmann
188! ql is calculated by calc_liquid_water_content
189!
190! 1076 2012-12-05 08:30:18Z hoffmann
191! Bugfix in output of ql
192!
193! 1065 2012-11-22 17:42:36Z hoffmann
194! Bugfix: Output of cross sections of ql
195!
196! 1053 2012-11-13 17:11:03Z hoffmann
197! +qr, nr, qc and cross sections
198!
199! 1036 2012-10-22 13:43:42Z raasch
200! code put under GPL (PALM 3.9)
201!
202! 1031 2012-10-19 14:35:30Z raasch
203! netCDF4 without parallel file support implemented
204!
205! 1007 2012-09-19 14:30:36Z franke
206! Bugfix: missing calculation of ql_vp added
207!
208! 978 2012-08-09 08:28:32Z fricke
209! +z0h
210!
211! Revision 1.1  1997/08/11 06:24:09  raasch
212! Initial revision
213!
214!
215! Description:
216! ------------
217!> Data output of cross-sections in netCDF format or binary format
218!> to be later converted to NetCDF by helper routine combine_plot_fields.
219!> Attention: The position of the sectional planes is still not always computed
220!> ---------  correctly. (zu is used always)!
221!------------------------------------------------------------------------------!
222 SUBROUTINE data_output_2d( mode, av )
223 
224
225    USE arrays_3d,                                                             &
226        ONLY:  dzw, e, heatflux_output_conversion, nc, nr, p, pt,              &
227               precipitation_amount, prr, q, qc, ql, ql_c, ql_v, ql_vp, qr,    &
228               rho_ocean, s, sa, tend, u, v, vpt, w, zu, zw,                   &
229               waterflux_output_conversion
230       
231    USE averaging
232       
233    USE cloud_parameters,                                                      &
234        ONLY:  cp, hyrho, l_d_cp, l_v, pt_d_t
235               
236    USE control_parameters,                                                    &
237        ONLY:  cloud_physics, data_output_2d_on_each_pe, data_output_xy,       &
238               data_output_xz, data_output_yz, do2d,                           &
239               do2d_xy_last_time, do2d_xy_time_count,                          &
240               do2d_xz_last_time, do2d_xz_time_count,                          &
241               do2d_yz_last_time, do2d_yz_time_count,                          &
242               ibc_uv_b, io_blocks, io_group, land_surface, message_string,    &
243               ntdim_2d_xy, ntdim_2d_xz, ntdim_2d_yz,                          &
244               psolver, section, simulated_time, simulated_time_chr,           &
245               time_since_reference_point, uv_exposure
246       
247    USE cpulog,                                                                &
248        ONLY:  cpu_log, log_point
249       
250    USE grid_variables,                                                        &
251        ONLY:  dx, dy
252
253    USE gust_mod,                                                              &
254        ONLY:  gust_data_output_2d, gust_module_enabled
255       
256    USE indices,                                                               &
257        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nz,   &
258               nzb, nzt, wall_flags_0
259               
260    USE kinds
261   
262    USE land_surface_model_mod,                                                &
263        ONLY:  lsm_data_output_2d, zs
264   
265#if defined( __netcdf )
266    USE NETCDF
267#endif
268
269    USE netcdf_interface,                                                      &
270        ONLY:  fill_value, id_set_xy, id_set_xz, id_set_yz, id_var_do2d,       &
271               id_var_time_xy, id_var_time_xz, id_var_time_yz, nc_stat,        &
272               netcdf_data_format, netcdf_handle_error
273
274    USE particle_attributes,                                                   &
275        ONLY:  grid_particles, number_of_particles, particle_advection_start,  &
276               particles, prt_count
277   
278    USE pegrid
279
280    USE radiation_model_mod,                                                   &
281        ONLY:  radiation, radiation_data_output_2d
282
283    USE surface_mod,                                                           &
284        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_def_h,           &
285               surf_lsm_h, surf_usm_h
286
287    USE turbulence_closure_mod,                                                &
288        ONLY:  tcm_data_output_2d
289
290    USE uv_exposure_model_mod,                                                 &
291        ONLY:  uvem_data_output_2d
292
293
294    IMPLICIT NONE
295
296    CHARACTER (LEN=2)  ::  do2d_mode    !<
297    CHARACTER (LEN=2)  ::  mode         !<
298    CHARACTER (LEN=4)  ::  grid         !<
299    CHARACTER (LEN=25) ::  section_chr  !<
300    CHARACTER (LEN=50) ::  rtext        !<
301   
302    INTEGER(iwp) ::  av        !<
303    INTEGER(iwp) ::  ngp       !<
304    INTEGER(iwp) ::  file_id   !<
305    INTEGER(iwp) ::  flag_nr   !< number of masking flag
306    INTEGER(iwp) ::  i         !<
307    INTEGER(iwp) ::  if        !<
308    INTEGER(iwp) ::  is        !<
309    INTEGER(iwp) ::  iis       !<
310    INTEGER(iwp) ::  j         !<
311    INTEGER(iwp) ::  k         !<
312    INTEGER(iwp) ::  l         !<
313    INTEGER(iwp) ::  layer_xy  !<
314    INTEGER(iwp) ::  m         !<
315    INTEGER(iwp) ::  n         !<
316    INTEGER(iwp) ::  nis       !<
317    INTEGER(iwp) ::  ns        !<
318    INTEGER(iwp) ::  nzb_do    !< lower limit of the data field (usually nzb)
319    INTEGER(iwp) ::  nzt_do    !< upper limit of the data field (usually nzt+1)
320    INTEGER(iwp) ::  psi       !<
321    INTEGER(iwp) ::  s_ind     !<
322    INTEGER(iwp) ::  sender    !<
323    INTEGER(iwp) ::  ind(4)    !<
324   
325    LOGICAL ::  found          !<
326    LOGICAL ::  resorted       !<
327    LOGICAL ::  two_d          !<
328   
329    REAL(wp) ::  mean_r        !<
330    REAL(wp) ::  s_r2          !<
331    REAL(wp) ::  s_r3          !<
332   
333    REAL(wp), DIMENSION(:), ALLOCATABLE     ::  level_z             !<
334    REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  local_2d            !<
335    REAL(wp), DIMENSION(:,:), ALLOCATABLE   ::  local_2d_l          !<
336
337    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  local_pf            !<
338    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  local_2d_sections   !<
339    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  local_2d_sections_l !<
340
341#if defined( __parallel )
342    REAL(wp), DIMENSION(:,:),   ALLOCATABLE ::  total_2d    !<
343#endif
344    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !<
345
346    NAMELIST /LOCAL/  rtext
347
348!
349!-- Immediate return, if no output is requested (no respective sections
350!-- found in parameter data_output)
351    IF ( mode == 'xy'  .AND.  .NOT. data_output_xy(av) )  RETURN
352    IF ( mode == 'xz'  .AND.  .NOT. data_output_xz(av) )  RETURN
353    IF ( mode == 'yz'  .AND.  .NOT. data_output_yz(av) )  RETURN
354
355    CALL cpu_log (log_point(3),'data_output_2d','start')
356
357    two_d = .FALSE.    ! local variable to distinguish between output of pure 2D
358                       ! arrays and cross-sections of 3D arrays.
359
360!
361!-- Depending on the orientation of the cross-section, the respective output
362!-- files have to be opened.
363    SELECT CASE ( mode )
364
365       CASE ( 'xy' )
366          s_ind = 1
367          ALLOCATE( level_z(nzb:nzt+1), local_2d(nxl:nxr,nys:nyn) )
368
369          IF ( netcdf_data_format > 4 )  THEN
370             ns = 1
371             DO WHILE ( section(ns,s_ind) /= -9999  .AND.  ns <= 100 )
372                ns = ns + 1
373             ENDDO
374             ns = ns - 1
375             ALLOCATE( local_2d_sections(nxl:nxr,nys:nyn,1:ns) )
376             local_2d_sections = 0.0_wp
377          ENDIF
378
379!
380!--       Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
381          IF ( myid == 0  .OR.  netcdf_data_format > 4 )  THEN
382             CALL check_open( 101+av*10 )
383          ENDIF
384          IF ( data_output_2d_on_each_pe )  THEN
385             CALL check_open( 21 )
386          ELSE
387             IF ( myid == 0 )  THEN
388#if defined( __parallel )
389                ALLOCATE( total_2d(0:nx,0:ny) )
390#endif
391             ENDIF
392          ENDIF
393
394       CASE ( 'xz' )
395          s_ind = 2
396          ALLOCATE( local_2d(nxl:nxr,nzb:nzt+1) )
397
398          IF ( netcdf_data_format > 4 )  THEN
399             ns = 1
400             DO WHILE ( section(ns,s_ind) /= -9999  .AND.  ns <= 100 )
401                ns = ns + 1
402             ENDDO
403             ns = ns - 1
404             ALLOCATE( local_2d_sections(nxl:nxr,1:ns,nzb:nzt+1) )
405             ALLOCATE( local_2d_sections_l(nxl:nxr,1:ns,nzb:nzt+1) )
406             local_2d_sections = 0.0_wp; local_2d_sections_l = 0.0_wp
407          ENDIF
408
409!
410!--       Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
411          IF ( myid == 0 .OR. netcdf_data_format > 4 )  THEN
412             CALL check_open( 102+av*10 )
413          ENDIF
414
415          IF ( data_output_2d_on_each_pe )  THEN
416             CALL check_open( 22 )
417          ELSE
418             IF ( myid == 0 )  THEN
419#if defined( __parallel )
420                ALLOCATE( total_2d(0:nx,nzb:nzt+1) )
421#endif
422             ENDIF
423          ENDIF
424
425       CASE ( 'yz' )
426          s_ind = 3
427          ALLOCATE( local_2d(nys:nyn,nzb:nzt+1) )
428
429          IF ( netcdf_data_format > 4 )  THEN
430             ns = 1
431             DO WHILE ( section(ns,s_ind) /= -9999  .AND.  ns <= 100 )
432                ns = ns + 1
433             ENDDO
434             ns = ns - 1
435             ALLOCATE( local_2d_sections(1:ns,nys:nyn,nzb:nzt+1) )
436             ALLOCATE( local_2d_sections_l(1:ns,nys:nyn,nzb:nzt+1) )
437             local_2d_sections = 0.0_wp; local_2d_sections_l = 0.0_wp
438          ENDIF
439
440!
441!--       Parallel netCDF4/HDF5 output is done on all PEs, all other on PE0 only
442          IF ( myid == 0 .OR. netcdf_data_format > 4 )  THEN
443             CALL check_open( 103+av*10 )
444          ENDIF
445
446          IF ( data_output_2d_on_each_pe )  THEN
447             CALL check_open( 23 )
448          ELSE
449             IF ( myid == 0 )  THEN
450#if defined( __parallel )
451                ALLOCATE( total_2d(0:ny,nzb:nzt+1) )
452#endif
453             ENDIF
454          ENDIF
455
456       CASE DEFAULT
457          message_string = 'unknown cross-section: ' // TRIM( mode )
458          CALL message( 'data_output_2d', 'PA0180', 1, 2, 0, 6, 0 )
459
460    END SELECT
461
462!
463!-- For parallel netcdf output the time axis must be limited. Return, if this
464!-- limit is exceeded. This could be the case, if the simulated time exceeds
465!-- the given end time by the length of the given output interval.
466    IF ( netcdf_data_format > 4 )  THEN
467       IF ( mode == 'xy'  .AND.  do2d_xy_time_count(av) + 1 >                  &
468            ntdim_2d_xy(av) )  THEN
469          WRITE ( message_string, * ) 'Output of xy cross-sections is not ',   &
470                          'given at t=', simulated_time, '&because the',       & 
471                          ' maximum number of output time levels is exceeded.'
472          CALL message( 'data_output_2d', 'PA0384', 0, 1, 0, 6, 0 )
473          CALL cpu_log( log_point(3), 'data_output_2d', 'stop' )
474          RETURN
475       ENDIF
476       IF ( mode == 'xz'  .AND.  do2d_xz_time_count(av) + 1 >                  &
477            ntdim_2d_xz(av) )  THEN
478          WRITE ( message_string, * ) 'Output of xz cross-sections is not ',   &
479                          'given at t=', simulated_time, '&because the',       & 
480                          ' maximum number of output time levels is exceeded.'
481          CALL message( 'data_output_2d', 'PA0385', 0, 1, 0, 6, 0 )
482          CALL cpu_log( log_point(3), 'data_output_2d', 'stop' )
483          RETURN
484       ENDIF
485       IF ( mode == 'yz'  .AND.  do2d_yz_time_count(av) + 1 >                  &
486            ntdim_2d_yz(av) )  THEN
487          WRITE ( message_string, * ) 'Output of yz cross-sections is not ',   &
488                          'given at t=', simulated_time, '&because the',       & 
489                          ' maximum number of output time levels is exceeded.'
490          CALL message( 'data_output_2d', 'PA0386', 0, 1, 0, 6, 0 )
491          CALL cpu_log( log_point(3), 'data_output_2d', 'stop' )
492          RETURN
493       ENDIF
494    ENDIF
495
496!
497!-- Allocate a temporary array for resorting (kji -> ijk).
498    ALLOCATE( local_pf(nxl:nxr,nys:nyn,nzb:nzt+1) )
499    local_pf = 0.0
500
501!
502!-- Loop of all variables to be written.
503!-- Output dimensions chosen
504    if = 1
505    l = MAX( 2, LEN_TRIM( do2d(av,if) ) )
506    do2d_mode = do2d(av,if)(l-1:l)
507
508    DO  WHILE ( do2d(av,if)(1:1) /= ' ' )
509
510       IF ( do2d_mode == mode )  THEN
511!
512!--       Set flag to steer output of radiation, land-surface, or user-defined
513!--       quantities
514          found = .FALSE.
515
516          nzb_do = nzb
517          nzt_do = nzt+1
518!
519!--       Before each output, set array local_pf to fill value
520          local_pf = fill_value
521!
522!--       Set masking flag for topography for not resorted arrays
523          flag_nr = 0
524         
525!
526!--       Store the array chosen on the temporary array.
527          resorted = .FALSE.
528          SELECT CASE ( TRIM( do2d(av,if) ) )
529             CASE ( 'e_xy', 'e_xz', 'e_yz' )
530                IF ( av == 0 )  THEN
531                   to_be_resorted => e
532                ELSE
533                   IF ( .NOT. ALLOCATED( e_av ) ) THEN
534                      ALLOCATE( e_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
535                      e_av = REAL( fill_value, KIND = wp )
536                   ENDIF
537                   to_be_resorted => e_av
538                ENDIF
539                IF ( mode == 'xy' )  level_z = zu
540
541             CASE ( 'lpt_xy', 'lpt_xz', 'lpt_yz' )
542                IF ( av == 0 )  THEN
543                   to_be_resorted => pt
544                ELSE
545                   IF ( .NOT. ALLOCATED( lpt_av ) ) THEN
546                      ALLOCATE( lpt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
547                      lpt_av = REAL( fill_value, KIND = wp )
548                   ENDIF
549                   to_be_resorted => lpt_av
550                ENDIF
551                IF ( mode == 'xy' )  level_z = zu
552
553             CASE ( 'lwp*_xy' )        ! 2d-array
554                IF ( av == 0 )  THEN
555                   DO  i = nxl, nxr
556                      DO  j = nys, nyn
557                         local_pf(i,j,nzb+1) = SUM( ql(nzb:nzt,j,i) *          &
558                                                    dzw(1:nzt+1) )
559                      ENDDO
560                   ENDDO
561                ELSE
562                   IF ( .NOT. ALLOCATED( lwp_av ) ) THEN
563                      ALLOCATE( lwp_av(nysg:nyng,nxlg:nxrg) )
564                      lwp_av = REAL( fill_value, KIND = wp )
565                   ENDIF
566                   DO  i = nxl, nxr
567                      DO  j = nys, nyn
568                         local_pf(i,j,nzb+1) = lwp_av(j,i)
569                      ENDDO
570                   ENDDO
571                ENDIF
572                resorted = .TRUE.
573                two_d = .TRUE.
574                level_z(nzb+1) = zu(nzb+1)
575
576             CASE ( 'nc_xy', 'nc_xz', 'nc_yz' )
577                IF ( av == 0 )  THEN
578                   to_be_resorted => nc
579                ELSE
580                   IF ( .NOT. ALLOCATED( nc_av ) ) THEN
581                      ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
582                      nc_av = REAL( fill_value, KIND = wp )
583                   ENDIF
584                   to_be_resorted => nc_av
585                ENDIF
586                IF ( mode == 'xy' )  level_z = zu
587
588             CASE ( 'nr_xy', 'nr_xz', 'nr_yz' )
589                IF ( av == 0 )  THEN
590                   to_be_resorted => nr
591                ELSE
592                   IF ( .NOT. ALLOCATED( nr_av ) ) THEN
593                      ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
594                      nr_av = REAL( fill_value, KIND = wp )
595                   ENDIF
596                   to_be_resorted => nr_av
597                ENDIF
598                IF ( mode == 'xy' )  level_z = zu
599
600             CASE ( 'ghf*_xy' )        ! 2d-array
601                IF ( av == 0 )  THEN
602                   DO  m = 1, surf_lsm_h%ns
603                      i                   = surf_lsm_h%i(m)           
604                      j                   = surf_lsm_h%j(m)
605                      local_pf(i,j,nzb+1) = surf_lsm_h%ghf(m)
606                   ENDDO
607                   DO  m = 1, surf_usm_h%ns
608                      i                   = surf_usm_h%i(m)           
609                      j                   = surf_usm_h%j(m)
610                      local_pf(i,j,nzb+1) = surf_usm_h%frac(ind_veg_wall,m)  *  &
611                                            surf_usm_h%wghf_eb(m)        +      &
612                                            surf_usm_h%frac(ind_pav_green,m) *  &
613                                            surf_usm_h%wghf_eb_green(m)  +      &
614                                            surf_usm_h%frac(ind_wat_win,m)   *  &
615                                            surf_usm_h%wghf_eb_window(m)
616                   ENDDO
617                ELSE
618                   IF ( .NOT. ALLOCATED( ghf_av ) ) THEN
619                      ALLOCATE( ghf_av(nysg:nyng,nxlg:nxrg) )
620                      ghf_av = REAL( fill_value, KIND = wp )
621                   ENDIF
622                   DO  i = nxl, nxr
623                      DO  j = nys, nyn
624                         local_pf(i,j,nzb+1) = ghf_av(j,i)
625                      ENDDO
626                   ENDDO
627                ENDIF
628
629                resorted = .TRUE.
630                two_d = .TRUE.
631                level_z(nzb+1) = zu(nzb+1)
632
633             CASE ( 'ol*_xy' )        ! 2d-array
634                IF ( av == 0 ) THEN
635                   DO  m = 1, surf_def_h(0)%ns
636                      i = surf_def_h(0)%i(m)
637                      j = surf_def_h(0)%j(m)
638                      local_pf(i,j,nzb+1) = surf_def_h(0)%ol(m)
639                   ENDDO
640                   DO  m = 1, surf_lsm_h%ns
641                      i = surf_lsm_h%i(m)
642                      j = surf_lsm_h%j(m)
643                      local_pf(i,j,nzb+1) = surf_lsm_h%ol(m)
644                   ENDDO
645                   DO  m = 1, surf_usm_h%ns
646                      i = surf_usm_h%i(m)
647                      j = surf_usm_h%j(m)
648                      local_pf(i,j,nzb+1) = surf_usm_h%ol(m)
649                   ENDDO
650                ELSE
651                   IF ( .NOT. ALLOCATED( ol_av ) ) THEN
652                      ALLOCATE( ol_av(nysg:nyng,nxlg:nxrg) )
653                      ol_av = REAL( fill_value, KIND = wp )
654                   ENDIF
655                   DO  i = nxl, nxr
656                      DO  j = nys, nyn
657                         local_pf(i,j,nzb+1) = ol_av(j,i)
658                      ENDDO
659                   ENDDO
660                ENDIF
661                resorted = .TRUE.
662                two_d = .TRUE.
663                level_z(nzb+1) = zu(nzb+1)
664
665             CASE ( 'p_xy', 'p_xz', 'p_yz' )
666                IF ( av == 0 )  THEN
667                   IF ( psolver /= 'sor' )  CALL exchange_horiz( p, nbgp )
668                   to_be_resorted => p
669                ELSE
670                   IF ( .NOT. ALLOCATED( p_av ) ) THEN
671                      ALLOCATE( p_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
672                      p_av = REAL( fill_value, KIND = wp )
673                   ENDIF
674                   IF ( psolver /= 'sor' )  CALL exchange_horiz( p_av, nbgp )
675                   to_be_resorted => p_av
676                ENDIF
677                IF ( mode == 'xy' )  level_z = zu
678
679             CASE ( 'pc_xy', 'pc_xz', 'pc_yz' )  ! particle concentration
680                IF ( av == 0 )  THEN
681                   IF ( simulated_time >= particle_advection_start )  THEN
682                      tend = prt_count
683!                      CALL exchange_horiz( tend, nbgp )
684                   ELSE
685                      tend = 0.0_wp
686                   ENDIF
687                   DO  i = nxl, nxr
688                      DO  j = nys, nyn
689                         DO  k = nzb, nzt+1
690                            local_pf(i,j,k) = tend(k,j,i)
691                         ENDDO
692                      ENDDO
693                   ENDDO
694                   resorted = .TRUE.
695                ELSE
696                   IF ( .NOT. ALLOCATED( pc_av ) ) THEN
697                      ALLOCATE( pc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
698                      pc_av = REAL( fill_value, KIND = wp )
699                   ENDIF
700!                   CALL exchange_horiz( pc_av, nbgp )
701                   to_be_resorted => pc_av
702                ENDIF
703
704             CASE ( 'pr_xy', 'pr_xz', 'pr_yz' )  ! mean particle radius (effective radius)
705                IF ( av == 0 )  THEN
706                   IF ( simulated_time >= particle_advection_start )  THEN
707                      DO  i = nxl, nxr
708                         DO  j = nys, nyn
709                            DO  k = nzb, nzt+1
710                               number_of_particles = prt_count(k,j,i)
711                               IF (number_of_particles <= 0)  CYCLE
712                               particles => grid_particles(k,j,i)%particles(1:number_of_particles)
713                               s_r2 = 0.0_wp
714                               s_r3 = 0.0_wp
715                               DO  n = 1, number_of_particles
716                                  IF ( particles(n)%particle_mask )  THEN
717                                     s_r2 = s_r2 + particles(n)%radius**2 * &
718                                            particles(n)%weight_factor
719                                     s_r3 = s_r3 + particles(n)%radius**3 * &
720                                            particles(n)%weight_factor
721                                  ENDIF
722                               ENDDO
723                               IF ( s_r2 > 0.0_wp )  THEN
724                                  mean_r = s_r3 / s_r2
725                               ELSE
726                                  mean_r = 0.0_wp
727                               ENDIF
728                               tend(k,j,i) = mean_r
729                            ENDDO
730                         ENDDO
731                      ENDDO
732!                      CALL exchange_horiz( tend, nbgp )
733                   ELSE
734                      tend = 0.0_wp
735                   ENDIF
736                   DO  i = nxl, nxr
737                      DO  j = nys, nyn
738                         DO  k = nzb, nzt+1
739                            local_pf(i,j,k) = tend(k,j,i)
740                         ENDDO
741                      ENDDO
742                   ENDDO
743                   resorted = .TRUE.
744                ELSE
745                   IF ( .NOT. ALLOCATED( pr_av ) ) THEN
746                      ALLOCATE( pr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
747                      pr_av = REAL( fill_value, KIND = wp )
748                   ENDIF
749!                   CALL exchange_horiz( pr_av, nbgp )
750                   to_be_resorted => pr_av
751                ENDIF
752
753             CASE ( 'pra*_xy' )        ! 2d-array / integral quantity => no av
754!                CALL exchange_horiz_2d( precipitation_amount )
755                   DO  i = nxl, nxr
756                      DO  j = nys, nyn
757                      local_pf(i,j,nzb+1) =  precipitation_amount(j,i)
758                   ENDDO
759                ENDDO
760                precipitation_amount = 0.0_wp   ! reset for next integ. interval
761                resorted = .TRUE.
762                two_d = .TRUE.
763                level_z(nzb+1) = zu(nzb+1)
764
765             CASE ( 'prr_xy', 'prr_xz', 'prr_yz' )
766                IF ( av == 0 )  THEN
767!                   CALL exchange_horiz( prr, nbgp )
768                   DO  i = nxl, nxr
769                      DO  j = nys, nyn
770                         DO  k = nzb, nzt+1
771                            local_pf(i,j,k) = prr(k,j,i) * hyrho(nzb+1)
772                         ENDDO
773                      ENDDO
774                   ENDDO
775                ELSE
776                   IF ( .NOT. ALLOCATED( prr_av ) ) THEN
777                      ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
778                      prr_av = REAL( fill_value, KIND = wp )
779                   ENDIF
780!                   CALL exchange_horiz( prr_av, nbgp )
781                   DO  i = nxl, nxr
782                      DO  j = nys, nyn
783                         DO  k = nzb, nzt+1
784                            local_pf(i,j,k) = prr_av(k,j,i) * hyrho(nzb+1)
785                         ENDDO
786                      ENDDO
787                   ENDDO
788                ENDIF
789                resorted = .TRUE.
790                IF ( mode == 'xy' )  level_z = zu
791
792             CASE ( 'pt_xy', 'pt_xz', 'pt_yz' )
793                IF ( av == 0 )  THEN
794                   IF ( .NOT. cloud_physics ) THEN
795                      to_be_resorted => pt
796                   ELSE
797                   DO  i = nxl, nxr
798                      DO  j = nys, nyn
799                            DO  k = nzb, nzt+1
800                               local_pf(i,j,k) = pt(k,j,i) + l_d_cp *          &
801                                                             pt_d_t(k) *       &
802                                                             ql(k,j,i)
803                            ENDDO
804                         ENDDO
805                      ENDDO
806                      resorted = .TRUE.
807                   ENDIF
808                ELSE
809                   IF ( .NOT. ALLOCATED( pt_av ) ) THEN
810                      ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
811                      pt_av = REAL( fill_value, KIND = wp )
812                   ENDIF
813                   to_be_resorted => pt_av
814                ENDIF
815                IF ( mode == 'xy' )  level_z = zu
816
817             CASE ( 'q_xy', 'q_xz', 'q_yz' )
818                IF ( av == 0 )  THEN
819                   to_be_resorted => q
820                ELSE
821                   IF ( .NOT. ALLOCATED( q_av ) ) THEN
822                      ALLOCATE( q_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
823                      q_av = REAL( fill_value, KIND = wp )
824                   ENDIF
825                   to_be_resorted => q_av
826                ENDIF
827                IF ( mode == 'xy' )  level_z = zu
828
829             CASE ( 'qc_xy', 'qc_xz', 'qc_yz' )
830                IF ( av == 0 )  THEN
831                   to_be_resorted => qc
832                ELSE
833                   IF ( .NOT. ALLOCATED( qc_av ) ) THEN
834                      ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
835                      qc_av = REAL( fill_value, KIND = wp )
836                   ENDIF
837                   to_be_resorted => qc_av
838                ENDIF
839                IF ( mode == 'xy' )  level_z = zu
840
841             CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
842                IF ( av == 0 )  THEN
843                   to_be_resorted => ql
844                ELSE
845                   IF ( .NOT. ALLOCATED( ql_av ) ) THEN
846                      ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
847                      ql_av = REAL( fill_value, KIND = wp )
848                   ENDIF
849                   to_be_resorted => ql_av
850                ENDIF
851                IF ( mode == 'xy' )  level_z = zu
852
853             CASE ( 'ql_c_xy', 'ql_c_xz', 'ql_c_yz' )
854                IF ( av == 0 )  THEN
855                   to_be_resorted => ql_c
856                ELSE
857                   IF ( .NOT. ALLOCATED( ql_c_av ) ) THEN
858                      ALLOCATE( ql_c_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
859                      ql_c_av = REAL( fill_value, KIND = wp )
860                   ENDIF
861                   to_be_resorted => ql_c_av
862                ENDIF
863                IF ( mode == 'xy' )  level_z = zu
864
865             CASE ( 'ql_v_xy', 'ql_v_xz', 'ql_v_yz' )
866                IF ( av == 0 )  THEN
867                   to_be_resorted => ql_v
868                ELSE
869                   IF ( .NOT. ALLOCATED( ql_v_av ) ) THEN
870                      ALLOCATE( ql_v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
871                      ql_v_av = REAL( fill_value, KIND = wp )
872                   ENDIF
873                   to_be_resorted => ql_v_av
874                ENDIF
875                IF ( mode == 'xy' )  level_z = zu
876
877             CASE ( 'ql_vp_xy', 'ql_vp_xz', 'ql_vp_yz' )
878                IF ( av == 0 )  THEN
879                   IF ( simulated_time >= particle_advection_start )  THEN
880                      DO  i = nxl, nxr
881                         DO  j = nys, nyn
882                            DO  k = nzb, nzt+1
883                               number_of_particles = prt_count(k,j,i)
884                               IF (number_of_particles <= 0)  CYCLE
885                               particles => grid_particles(k,j,i)%particles(1:number_of_particles)
886                               DO  n = 1, number_of_particles
887                                  IF ( particles(n)%particle_mask )  THEN
888                                     tend(k,j,i) =  tend(k,j,i) +                 &
889                                                    particles(n)%weight_factor /  &
890                                                    prt_count(k,j,i)
891                                  ENDIF
892                               ENDDO
893                            ENDDO
894                         ENDDO
895                      ENDDO
896!                      CALL exchange_horiz( tend, nbgp )
897                   ELSE
898                      tend = 0.0_wp
899                   ENDIF
900                   DO  i = nxl, nxr
901                      DO  j = nys, nyn
902                         DO  k = nzb, nzt+1
903                            local_pf(i,j,k) = tend(k,j,i)
904                         ENDDO
905                      ENDDO
906                   ENDDO
907                   resorted = .TRUE.
908                ELSE
909                   IF ( .NOT. ALLOCATED( ql_vp_av ) ) THEN
910                      ALLOCATE( ql_vp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
911                      ql_vp_av = REAL( fill_value, KIND = wp )
912                   ENDIF
913!                   CALL exchange_horiz( ql_vp_av, nbgp )
914                   to_be_resorted => ql_vp_av
915                ENDIF
916                IF ( mode == 'xy' )  level_z = zu
917
918             CASE ( 'qr_xy', 'qr_xz', 'qr_yz' )
919                IF ( av == 0 )  THEN
920                   to_be_resorted => qr
921                ELSE
922                   IF ( .NOT. ALLOCATED( qr_av ) ) THEN
923                      ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
924                      qr_av = REAL( fill_value, KIND = wp )
925                   ENDIF
926                   to_be_resorted => qr_av
927                ENDIF
928                IF ( mode == 'xy' )  level_z = zu
929
930             CASE ( 'qsws*_xy' )        ! 2d-array
931                IF ( av == 0 ) THEN
932!
933!--                In case of default surfaces, clean-up flux by density.
934!--                In case of land- and urban-surfaces, convert fluxes into
935!--                dynamic units
936                   DO  m = 1, surf_def_h(0)%ns
937                      i = surf_def_h(0)%i(m)
938                      j = surf_def_h(0)%j(m)
939                      k = surf_def_h(0)%k(m)
940                      local_pf(i,j,nzb+1) = surf_def_h(0)%qsws(m) *            &
941                                            waterflux_output_conversion(k)
942                   ENDDO
943                   DO  m = 1, surf_lsm_h%ns
944                      i = surf_lsm_h%i(m)
945                      j = surf_lsm_h%j(m)
946                      k = surf_lsm_h%k(m)
947                      local_pf(i,j,nzb+1) = surf_lsm_h%qsws(m) * l_v
948                   ENDDO
949                   DO  m = 1, surf_usm_h%ns
950                      i = surf_usm_h%i(m)
951                      j = surf_usm_h%j(m)
952                      k = surf_usm_h%k(m)
953                      local_pf(i,j,nzb+1) = surf_usm_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, two_d, 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, mode,&
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.