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

Last change on this file since 3679 was 3655, checked in by knoop, 6 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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