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

Last change on this file since 3882 was 3766, checked in by raasch, 5 years ago

unused_variables removed, bugfix in im_define_netcdf_grid argument list, trim added to avoid truncation compiler warnings, save attribute added to local targets to avoid outlive pointer target warning, first argument removed from module_interface_rrd_*, file module_interface reformatted with respect to coding standards, bugfix in surface_data_output_rrd_local (variable k removed)

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