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

Last change on this file since 3896 was 3885, checked in by kanani, 6 years ago

restructure/add location/debug messages

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