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

Last change on this file since 3569 was 3569, checked in by kanani, 3 years ago

Fix for biomet output (ticket:757), merge of uv_exposure into biometeorology_mod

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