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

Last change on this file since 3560 was 3554, checked in by gronemeier, 5 years ago

renamed variable if to ivar; add variable description

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