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

Last change on this file since 3619 was 3597, checked in by maronga, 6 years ago

revised calculation of near surface air potential temperature

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