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

Last change on this file since 3589 was 3589, checked in by suehring, 6 years ago

Remove erroneous UTF encoding; last commit documented

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