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

Last change on this file since 4088 was 4048, checked in by knoop, 5 years ago

Moved turbulence_closure_mod calls into module_interface

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