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

Last change on this file since 3010 was 3004, checked in by Giersch, 7 years ago

precipitation_rate removed, further allocation checks for data output of averaged quantities implemented, double CALL of flow_statistics at the beginning of time_integration removed, further minor bugfixes, comments added

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