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

Last change on this file since 3582 was 3582, checked in by suehring, 3 years ago

Merge branch salsa with trunk

  • Property svn:keywords set to Id
File size: 96.8 KB
Line 
1!> @file data_output_2d.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22! Move the control parameter "salsa" from salsa_mod to control_parameters
23! (M. Kurppa)
24!
25! Former revisions:
26! -----------------
27! $Id: data_output_2d.f90 3582 2018-11-29 19:16:36Z suehring $
28! Remove fill_value from bio_data_output_2d call
29! dom_dwd_user, Schrempf:
30! Clean up of biometeorology calls,
31! remove uv exposure model code, this is now part of biometeorology_mod.
32!
33! 3554 2018-11-22 11:24:52Z gronemeier
34! - add variable description
35! - rename variable 'if' into 'ivar'
36! - removed namelist LOCAL
37! - removed variable rtext
38!
39! 3525 2018-11-14 16:06:14Z kanani
40! Changes related to clean-up of biometeorology (dom_dwd_user)
41!
42! 3467 2018-10-30 19:05:21Z suehring
43! Implementation of a new aerosol module salsa.
44!
45! 3448 2018-10-29 18:14:31Z kanani
46! Add biometeorology
47!
48! 3421 2018-10-24 18:39:32Z gronemeier
49! Renamed output variables
50!
51! 3419 2018-10-24 17:27:31Z gronemeier
52! minor formatting (kanani)
53! chem_data_output_2d subroutine added (basit)
54!
55! 3294 2018-10-01 02:37:10Z raasch
56! changes concerning modularization of ocean option
57!
58! 3274 2018-09-24 15:42:55Z knoop
59! Modularization of all bulk cloud physics code components
60!
61! 3241 2018-09-12 15:02:00Z raasch
62! unused variables removed
63!
64! 3176 2018-07-26 17:12:48Z suehring
65! Remove output of latent heat flux at urban-surfaces and set fill values
66! instead
67!
68! 3052 2018-05-31 06:11:20Z raasch
69! Do not open FORTRAN binary files in case of parallel netCDF I/O
70!
71! 3049 2018-05-29 13:52:36Z Giersch
72! Error messages revised
73!
74! 3045 2018-05-28 07:55:41Z Giersch
75! Error messages revised
76!
77! 3014 2018-05-09 08:42:38Z maronga
78! Added nzb_do and nzt_do for some modules for 2d output
79!
80! 3004 2018-04-27 12:33:25Z Giersch
81! precipitation_rate removed, case prr*_xy removed, to_be_resorted have to point
82! to ql_vp_av and not to ql_vp, allocation checks implemented (averaged data
83! will be assigned to fill values if no allocation happened so far)   
84!
85! 2963 2018-04-12 14:47:44Z suehring
86! Introduce index for vegetation/wall, pavement/green-wall and water/window
87! surfaces, for clearer access of surface fraction, albedo, emissivity, etc. .
88!
89! 2817 2018-02-19 16:32:21Z knoop
90! Preliminary gust module interface implemented
91!
92! 2805 2018-02-14 17:00:09Z suehring
93! Consider also default-type surfaces for surface temperature output.
94!
95! 2797 2018-02-08 13:24:35Z suehring
96! Enable output of ground-heat flux also at urban surfaces.
97!
98! 2743 2018-01-12 16:03:39Z suehring
99! In case of natural- and urban-type surfaces output surfaces fluxes in W/m2.
100!
101! 2742 2018-01-12 14:59:47Z suehring
102! Enable output of surface temperature
103!
104! 2735 2018-01-11 12:01:27Z suehring
105! output of r_a moved from land-surface to consider also urban-type surfaces
106!
107! 2718 2018-01-02 08:49:38Z maronga
108! Corrected "Former revisions" section
109!
110! 2696 2017-12-14 17:12:51Z kanani
111! Change in file header (GPL part)
112! Implementation of uv exposure model (FK)
113! Implementation of turbulence_closure_mod (TG)
114! Set fill values at topography grid points or e.g. non-natural-type surface
115! in case of LSM output (MS)
116!
117! 2512 2017-10-04 08:26:59Z raasch
118! upper bounds of cross section output changed from nx+1,ny+1 to nx,ny
119! no output of ghost layer data
120!
121! 2292 2017-06-20 09:51:42Z schwenkel
122! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
123! includes two more prognostic equations for cloud drop concentration (nc) 
124! and cloud water content (qc).
125!
126! 2277 2017-06-12 10:47:51Z kanani
127! Removed unused variables do2d_xy_n, do2d_xz_n, do2d_yz_n
128!
129! 2233 2017-05-30 18:08:54Z suehring
130!
131! 2232 2017-05-30 17:47:52Z suehring
132! Adjustments to new surface concept
133!
134!
135! 2190 2017-03-21 12:16:43Z raasch
136! bugfix for r2031: string rho replaced by rho_ocean
137!
138! 2031 2016-10-21 15:11:58Z knoop
139! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
140!
141! 2000 2016-08-20 18:09:15Z knoop
142! Forced header and separation lines into 80 columns
143!
144! 1980 2016-07-29 15:51:57Z suehring
145! Bugfix, in order to steer user-defined output, setting flag found explicitly
146! to .F.
147!
148! 1976 2016-07-27 13:28:04Z maronga
149! Output of radiation quantities is now done directly in the respective module
150!
151! 1972 2016-07-26 07:52:02Z maronga
152! Output of land surface quantities is now done directly in the respective
153! module
154!
155! 1960 2016-07-12 16:34:24Z suehring
156! Scalar surface flux added
157! Rename INTEGER variable s into s_ind, as s is already assigned to scalar
158!
159! 1849 2016-04-08 11:33:18Z hoffmann
160! precipitation_amount, precipitation_rate, prr moved to arrays_3d
161!
162! 1822 2016-04-07 07:49:42Z hoffmann
163! Output of bulk cloud physics simplified.
164!
165! 1788 2016-03-10 11:01:04Z maronga
166! Added output of z0q
167!
168! 1783 2016-03-06 18:36:17Z raasch
169! name change of netcdf routines and module + related changes
170!
171! 1745 2016-02-05 13:06:51Z gronemeier
172! Bugfix: test if time axis limit exceeds moved to point after call of check_open
173!
174! 1703 2015-11-02 12:38:44Z raasch
175! bugfix for output of single (*) xy-sections in case of parallel netcdf I/O
176!
177! 1701 2015-11-02 07:43:04Z maronga
178! Bugfix in output of RRTGM data
179!
180! 1691 2015-10-26 16:17:44Z maronga
181! Added output of Obukhov length (ol) and radiative heating rates  for RRTMG.
182! Formatting corrections.
183!
184! 1682 2015-10-07 23:56:08Z knoop
185! Code annotations made doxygen readable
186!
187! 1585 2015-04-30 07:05:52Z maronga
188! Added support for RRTMG
189!
190! 1555 2015-03-04 17:44:27Z maronga
191! Added output of r_a and r_s
192!
193! 1551 2015-03-03 14:18:16Z maronga
194! Added suppport for land surface model and radiation model output. In the course
195! of this action, the limits for vertical loops have been changed (from nzb and
196! nzt+1 to nzb_do and nzt_do, respectively in order to allow soil model output).
197! Moreover, a new vertical grid zs was introduced.
198!
199! 1359 2014-04-11 17:15:14Z hoffmann
200! New particle structure integrated.
201!
202! 1353 2014-04-08 15:21:23Z heinze
203! REAL constants provided with KIND-attribute
204!
205! 1327 2014-03-21 11:00:16Z raasch
206! parts concerning iso2d output removed,
207! -netcdf output queries
208!
209! 1320 2014-03-20 08:40:49Z raasch
210! ONLY-attribute added to USE-statements,
211! kind-parameters added to all INTEGER and REAL declaration statements,
212! kinds are defined in new module kinds,
213! revision history before 2012 removed,
214! comment fields (!:) to be used for variable explanations added to
215! all variable declaration statements
216!
217! 1318 2014-03-17 13:35:16Z raasch
218! barrier argument removed from cpu_log.
219! module interfaces removed
220!
221! 1311 2014-03-14 12:13:39Z heinze
222! bugfix: close #if defined( __netcdf )
223!
224! 1308 2014-03-13 14:58:42Z fricke
225! +local_2d_sections, local_2d_sections_l, ns
226! Check, if the limit of the time dimension is exceeded for parallel output
227! To increase the performance for parallel output, the following is done:
228! - Update of time axis is only done by PE0
229! - Cross sections are first stored on a local array and are written
230!   collectively to the output file by all PEs.
231!
232! 1115 2013-03-26 18:16:16Z hoffmann
233! ql is calculated by calc_liquid_water_content
234!
235! 1076 2012-12-05 08:30:18Z hoffmann
236! Bugfix in output of ql
237!
238! 1065 2012-11-22 17:42:36Z hoffmann
239! Bugfix: Output of cross sections of ql
240!
241! 1053 2012-11-13 17:11:03Z hoffmann
242! +qr, nr, qc and cross sections
243!
244! 1036 2012-10-22 13:43:42Z raasch
245! code put under GPL (PALM 3.9)
246!
247! 1031 2012-10-19 14:35:30Z raasch
248! netCDF4 without parallel file support implemented
249!
250! 1007 2012-09-19 14:30:36Z franke
251! Bugfix: missing calculation of ql_vp added
252!
253! 978 2012-08-09 08:28:32Z fricke
254! +z0h
255!
256! Revision 1.1  1997/08/11 06:24:09  raasch
257! Initial revision
258!
259!
260! Description:
261! ------------
262!> Data output of cross-sections in netCDF format or binary format
263!> to be later converted to NetCDF by helper routine combine_plot_fields.
264!> Attention: The position of the sectional planes is still not always computed
265!> ---------  correctly. (zu is used always)!
266!------------------------------------------------------------------------------!
267 SUBROUTINE data_output_2d( mode, av )
268 
269
270    USE arrays_3d,                                                             &
271        ONLY:  dzw, e, heatflux_output_conversion, nc, nr, p, pt,              &
272               precipitation_amount, prr, q, qc, ql, ql_c, ql_v, qr, s, tend,  &
273               u, v, vpt, w, zu, zw, waterflux_output_conversion, hyrho, d_exner
274
275    USE averaging
276
277    USE basic_constants_and_equations_mod,                                     &
278        ONLY:  c_p, lv_d_cp, l_v
279
280    USE biometeorology_mod,                                                    &
281        ONLY:  bio_data_output_2d
282
283    USE bulk_cloud_model_mod,                                                  &
284        ONLY:  bulk_cloud_model, bcm_data_output_2d
285
286    USE chemistry_model_mod,                                                   &
287        ONLY:  chem_data_output_2d
288
289    USE control_parameters,                                                    &
290        ONLY:  air_chemistry, biometeorology, data_output_2d_on_each_pe,       &
291               data_output_xy, data_output_xz, data_output_yz, do2d,           &
292               do2d_xy_last_time, do2d_xy_time_count,                          &
293               do2d_xz_last_time, do2d_xz_time_count,                          &
294               do2d_yz_last_time, do2d_yz_time_count,                          &
295               ibc_uv_b, io_blocks, io_group, land_surface, message_string,    &
296               ntdim_2d_xy, ntdim_2d_xz, ntdim_2d_yz,                          &
297               ocean_mode, psolver, salsa, section, simulated_time,            &
298               time_since_reference_point
299
300    USE cpulog,                                                                &
301        ONLY:  cpu_log, log_point
302       
303    USE gust_mod,                                                              &
304        ONLY:  gust_data_output_2d, gust_module_enabled
305
306    USE indices,                                                               &
307        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,       &
308               nzb, nzt, wall_flags_0
309
310    USE kinds
311
312    USE land_surface_model_mod,                                                &
313        ONLY:  lsm_data_output_2d, zs
314
315#if defined( __netcdf )
316    USE NETCDF
317#endif
318
319    USE netcdf_interface,                                                      &
320        ONLY:  fill_value, id_set_xy, id_set_xz, id_set_yz, id_var_do2d,       &
321               id_var_time_xy, id_var_time_xz, id_var_time_yz, nc_stat,        &
322               netcdf_data_format, netcdf_handle_error
323
324    USE ocean_mod,                                                             &
325        ONLY:  ocean_data_output_2d
326
327    USE particle_attributes,                                                   &
328        ONLY:  grid_particles, number_of_particles, particle_advection_start,  &
329               particles, prt_count
330   
331    USE pegrid
332
333    USE radiation_model_mod,                                                   &
334        ONLY:  radiation, radiation_data_output_2d
335   
336    USE salsa_mod,                                                             &
337        ONLY:  salsa_data_output_2d   
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=', simulated_time, '&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=', simulated_time, '&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=', simulated_time, '&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 ( simulated_time >= 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 ( simulated_time >= 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 ( simulated_time >= 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 ( 'w_xy', 'w_xz', 'w_yz' )
1212                flag_nr = 3
1213                IF ( av == 0 )  THEN
1214                   to_be_resorted => w
1215                ELSE
1216                   IF ( .NOT. ALLOCATED( w_av ) ) THEN
1217                      ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
1218                      w_av = REAL( fill_value, KIND = wp )
1219                   ENDIF
1220                   to_be_resorted => w_av
1221                ENDIF
1222                IF ( mode == 'xy' )  level_z = zw
1223
1224             CASE ( 'z0*_xy' )        ! 2d-array
1225                IF ( av == 0 ) THEN
1226                   DO  m = 1, surf_def_h(0)%ns
1227                      i = surf_def_h(0)%i(m)
1228                      j = surf_def_h(0)%j(m)
1229                      local_pf(i,j,nzb+1) = surf_def_h(0)%z0(m)
1230                   ENDDO
1231                   DO  m = 1, surf_lsm_h%ns
1232                      i = surf_lsm_h%i(m)
1233                      j = surf_lsm_h%j(m)
1234                      local_pf(i,j,nzb+1) = surf_lsm_h%z0(m)
1235                   ENDDO
1236                   DO  m = 1, surf_usm_h%ns
1237                      i = surf_usm_h%i(m)
1238                      j = surf_usm_h%j(m)
1239                      local_pf(i,j,nzb+1) = surf_usm_h%z0(m)
1240                   ENDDO
1241                ELSE
1242                   IF ( .NOT. ALLOCATED( z0_av ) ) THEN
1243                      ALLOCATE( z0_av(nysg:nyng,nxlg:nxrg) )
1244                      z0_av = REAL( fill_value, KIND = wp )
1245                   ENDIF
1246                   DO  i = nxl, nxr
1247                      DO  j = nys, nyn
1248                         local_pf(i,j,nzb+1) =  z0_av(j,i)
1249                      ENDDO
1250                   ENDDO
1251                ENDIF
1252                resorted = .TRUE.
1253                two_d = .TRUE.
1254                level_z(nzb+1) = zu(nzb+1)
1255
1256             CASE ( 'z0h*_xy' )        ! 2d-array
1257                IF ( av == 0 ) THEN
1258                   DO  m = 1, surf_def_h(0)%ns
1259                      i = surf_def_h(0)%i(m)
1260                      j = surf_def_h(0)%j(m)
1261                      local_pf(i,j,nzb+1) = surf_def_h(0)%z0h(m)
1262                   ENDDO
1263                   DO  m = 1, surf_lsm_h%ns
1264                      i = surf_lsm_h%i(m)
1265                      j = surf_lsm_h%j(m)
1266                      local_pf(i,j,nzb+1) = surf_lsm_h%z0h(m)
1267                   ENDDO
1268                   DO  m = 1, surf_usm_h%ns
1269                      i = surf_usm_h%i(m)
1270                      j = surf_usm_h%j(m)
1271                      local_pf(i,j,nzb+1) = surf_usm_h%z0h(m)
1272                   ENDDO
1273                ELSE
1274                   IF ( .NOT. ALLOCATED( z0h_av ) ) THEN
1275                      ALLOCATE( z0h_av(nysg:nyng,nxlg:nxrg) )
1276                      z0h_av = REAL( fill_value, KIND = wp )
1277                   ENDIF
1278                   DO  i = nxl, nxr
1279                      DO  j = nys, nyn
1280                         local_pf(i,j,nzb+1) =  z0h_av(j,i)
1281                      ENDDO
1282                   ENDDO
1283                ENDIF
1284                resorted = .TRUE.
1285                two_d = .TRUE.
1286                level_z(nzb+1) = zu(nzb+1)
1287
1288             CASE ( 'z0q*_xy' )        ! 2d-array
1289                IF ( av == 0 ) THEN
1290                   DO  m = 1, surf_def_h(0)%ns
1291                      i = surf_def_h(0)%i(m)
1292                      j = surf_def_h(0)%j(m)
1293                      local_pf(i,j,nzb+1) = surf_def_h(0)%z0q(m)
1294                   ENDDO
1295                   DO  m = 1, surf_lsm_h%ns
1296                      i = surf_lsm_h%i(m)
1297                      j = surf_lsm_h%j(m)
1298                      local_pf(i,j,nzb+1) = surf_lsm_h%z0q(m)
1299                   ENDDO
1300                   DO  m = 1, surf_usm_h%ns
1301                      i = surf_usm_h%i(m)
1302                      j = surf_usm_h%j(m)
1303                      local_pf(i,j,nzb+1) = surf_usm_h%z0q(m)
1304                   ENDDO
1305                ELSE
1306                   IF ( .NOT. ALLOCATED( z0q_av ) ) THEN
1307                      ALLOCATE( z0q_av(nysg:nyng,nxlg:nxrg) )
1308                      z0q_av = REAL( fill_value, KIND = wp )
1309                   ENDIF
1310                   DO  i = nxl, nxr
1311                      DO  j = nys, nyn
1312                         local_pf(i,j,nzb+1) =  z0q_av(j,i)
1313                      ENDDO
1314                   ENDDO
1315                ENDIF
1316                resorted = .TRUE.
1317                two_d = .TRUE.
1318                level_z(nzb+1) = zu(nzb+1)
1319
1320             CASE DEFAULT
1321
1322!
1323!--             Quantities of other modules
1324                IF ( .NOT. found  .AND.  bulk_cloud_model )  THEN
1325                   CALL bcm_data_output_2d( av, do2d(av,ivar), found, grid, mode,&
1326                                            local_pf, two_d, nzb_do, nzt_do )
1327                ENDIF
1328
1329                IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
1330                   CALL gust_data_output_2d( av, do2d(av,ivar), found, grid,     &
1331                                             local_pf, two_d, nzb_do, nzt_do )
1332                ENDIF
1333
1334                IF ( .NOT. found  .AND.  biometeorology )  THEN
1335                   CALL bio_data_output_2d( av, do2d(av,ivar), found, grid,      &
1336                                            local_pf, two_d, nzb_do, nzt_do )
1337                ENDIF
1338
1339                IF ( .NOT. found  .AND.  land_surface )  THEN
1340                   CALL lsm_data_output_2d( av, do2d(av,ivar), found, grid, mode,&
1341                                            local_pf, two_d, nzb_do, nzt_do )
1342                ENDIF
1343
1344                IF ( .NOT. found  .AND.  ocean_mode )  THEN
1345                   CALL ocean_data_output_2d( av, do2d(av,ivar), found, grid,    &
1346                                              mode, local_pf, nzb_do, nzt_do )
1347                ENDIF
1348
1349                IF ( .NOT. found  .AND.  radiation )  THEN
1350                   CALL radiation_data_output_2d( av, do2d(av,ivar), found, grid,&
1351                                                  mode, local_pf, two_d,       &
1352                                                  nzb_do, nzt_do  )
1353                ENDIF
1354               
1355                IF ( .NOT. found  .AND.  salsa )  THEN
1356                   CALL salsa_data_output_2d( av, do2d(av,ivar), found, grid,  &
1357                                              mode, local_pf, two_d, nzb_do,   &
1358                                              nzt_do)
1359                ENDIF                 
1360
1361                IF ( .NOT. found  .AND.  air_chemistry )  THEN
1362                   CALL chem_data_output_2d( av, do2d(av,ivar), found, grid, mode, &
1363                                             local_pf, two_d, nzb_do, nzt_do, fill_value )
1364                ENDIF
1365!
1366!--             User defined quantities
1367                IF ( .NOT. found )  THEN
1368                   CALL user_data_output_2d( av, do2d(av,ivar), found, grid,     &
1369                                             local_pf, two_d, nzb_do, nzt_do )
1370                ENDIF
1371
1372                resorted = .TRUE.
1373
1374                IF ( grid == 'zu' )  THEN
1375                   IF ( mode == 'xy' )  level_z = zu
1376                ELSEIF ( grid == 'zw' )  THEN
1377                   IF ( mode == 'xy' )  level_z = zw
1378                ELSEIF ( grid == 'zu1' ) THEN
1379                   IF ( mode == 'xy' )  level_z(nzb+1) = zu(nzb+1)
1380                ELSEIF ( grid == 'zs' ) THEN
1381                   IF ( mode == 'xy' )  level_z = zs
1382                ENDIF
1383
1384                IF ( .NOT. found )  THEN
1385                   message_string = 'no output provided for: ' //              &
1386                                    TRIM( do2d(av,ivar) )
1387                   CALL message( 'data_output_2d', 'PA0181', 0, 0, 0, 6, 0 )
1388                ENDIF
1389
1390          END SELECT
1391
1392!
1393!--       Resort the array to be output, if not done above. Flag topography
1394!--       grid points with fill values, using the corresponding maksing flag.
1395          IF ( .NOT. resorted )  THEN
1396             DO  i = nxl, nxr
1397                DO  j = nys, nyn
1398                   DO  k = nzb_do, nzt_do
1399                      local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),          &
1400                                               REAL( fill_value, KIND = wp ),  &
1401                                               BTEST( wall_flags_0(k,j,i),     &
1402                                                      flag_nr ) ) 
1403                   ENDDO
1404                ENDDO
1405             ENDDO
1406          ENDIF
1407
1408!
1409!--       Output of the individual cross-sections, depending on the cross-
1410!--       section mode chosen.
1411          is = 1
1412   loop1: DO WHILE ( section(is,s_ind) /= -9999  .OR.  two_d )
1413
1414             SELECT CASE ( mode )
1415
1416                CASE ( 'xy' )
1417!
1418!--                Determine the cross section index
1419                   IF ( two_d )  THEN
1420                      layer_xy = nzb+1
1421                   ELSE
1422                      layer_xy = section(is,s_ind)
1423                   ENDIF
1424
1425!
1426!--                Exit the loop for layers beyond the data output domain
1427!--                (used for soil model)
1428                   IF ( layer_xy > nzt_do )  THEN
1429                      EXIT loop1
1430                   ENDIF
1431
1432!
1433!--                Update the netCDF xy cross section time axis.
1434!--                In case of parallel output, this is only done by PE0
1435!--                to increase the performance.
1436                   IF ( simulated_time /= do2d_xy_last_time(av) )  THEN
1437                      do2d_xy_time_count(av) = do2d_xy_time_count(av) + 1
1438                      do2d_xy_last_time(av)  = simulated_time
1439                      IF ( myid == 0 )  THEN
1440                         IF ( .NOT. data_output_2d_on_each_pe  &
1441                              .OR.  netcdf_data_format > 4 )   &
1442                         THEN
1443#if defined( __netcdf )
1444                            nc_stat = NF90_PUT_VAR( id_set_xy(av),             &
1445                                                    id_var_time_xy(av),        &
1446                                             (/ time_since_reference_point /), &
1447                                         start = (/ do2d_xy_time_count(av) /), &
1448                                                    count = (/ 1 /) )
1449                            CALL netcdf_handle_error( 'data_output_2d', 53 )
1450#endif
1451                         ENDIF
1452                      ENDIF
1453                   ENDIF
1454!
1455!--                If required, carry out averaging along z
1456                   IF ( section(is,s_ind) == -1  .AND.  .NOT. two_d )  THEN
1457
1458                      local_2d = 0.0_wp
1459!
1460!--                   Carry out the averaging (all data are on the PE)
1461                      DO  k = nzb_do, nzt_do
1462                         DO  j = nys, nyn
1463                            DO  i = nxl, nxr
1464                               local_2d(i,j) = local_2d(i,j) + local_pf(i,j,k)
1465                            ENDDO
1466                         ENDDO
1467                      ENDDO
1468
1469                      local_2d = local_2d / ( nzt_do - nzb_do + 1.0_wp)
1470
1471                   ELSE
1472!
1473!--                   Just store the respective section on the local array
1474                      local_2d = local_pf(:,:,layer_xy)
1475
1476                   ENDIF
1477
1478#if defined( __parallel )
1479                   IF ( netcdf_data_format > 4 )  THEN
1480!
1481!--                   Parallel output in netCDF4/HDF5 format.
1482                      IF ( two_d ) THEN
1483                         iis = 1
1484                      ELSE
1485                         iis = is
1486                      ENDIF
1487
1488#if defined( __netcdf )
1489!
1490!--                   For parallel output, all cross sections are first stored
1491!--                   here on a local array and will be written to the output
1492!--                   file afterwards to increase the performance.
1493                      DO  i = nxl, nxr
1494                         DO  j = nys, nyn
1495                            local_2d_sections(i,j,iis) = local_2d(i,j)
1496                         ENDDO
1497                      ENDDO
1498#endif
1499                   ELSE
1500
1501                      IF ( data_output_2d_on_each_pe )  THEN
1502!
1503!--                      Output of partial arrays on each PE
1504#if defined( __netcdf )
1505                         IF ( myid == 0 )  THEN
1506                            WRITE ( 21 )  time_since_reference_point,          &
1507                                          do2d_xy_time_count(av), av
1508                         ENDIF
1509#endif
1510                         DO  i = 0, io_blocks-1
1511                            IF ( i == io_group )  THEN
1512                               WRITE ( 21 )  nxl, nxr, nys, nyn, nys, nyn
1513                               WRITE ( 21 )  local_2d
1514                            ENDIF
1515#if defined( __parallel )
1516                            CALL MPI_BARRIER( comm2d, ierr )
1517#endif
1518                         ENDDO
1519
1520                      ELSE
1521!
1522!--                      PE0 receives partial arrays from all processors and
1523!--                      then outputs them. Here a barrier has to be set,
1524!--                      because otherwise "-MPI- FATAL: Remote protocol queue
1525!--                      full" may occur.
1526                         CALL MPI_BARRIER( comm2d, ierr )
1527
1528                         ngp = ( nxr-nxl+1 ) * ( nyn-nys+1 )
1529                         IF ( myid == 0 )  THEN
1530!
1531!--                         Local array can be relocated directly.
1532                            total_2d(nxl:nxr,nys:nyn) = local_2d
1533!
1534!--                         Receive data from all other PEs.
1535                            DO  n = 1, numprocs-1
1536!
1537!--                            Receive index limits first, then array.
1538!--                            Index limits are received in arbitrary order from
1539!--                            the PEs.
1540                               CALL MPI_RECV( ind(1), 4, MPI_INTEGER,          &
1541                                              MPI_ANY_SOURCE, 0, comm2d,       &
1542                                              status, ierr )
1543                               sender = status(MPI_SOURCE)
1544                               DEALLOCATE( local_2d )
1545                               ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
1546                               CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp,    &
1547                                              MPI_REAL, sender, 1, comm2d,     &
1548                                              status, ierr )
1549                               total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
1550                            ENDDO
1551!
1552!--                         Relocate the local array for the next loop increment
1553                            DEALLOCATE( local_2d )
1554                            ALLOCATE( local_2d(nxl:nxr,nys:nyn) )
1555
1556#if defined( __netcdf )
1557                            IF ( two_d ) THEN
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, 1, do2d_xy_time_count(av) /), &
1562                                             count = (/ nx+1, ny+1, 1, 1 /) )
1563                            ELSE
1564                               nc_stat = NF90_PUT_VAR( id_set_xy(av),       &
1565                                                       id_var_do2d(av,ivar),  &
1566                                                       total_2d(0:nx,0:ny), &
1567                            start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
1568                                             count = (/ nx+1, ny+1, 1, 1 /) )
1569                            ENDIF
1570                            CALL netcdf_handle_error( 'data_output_2d', 54 )
1571#endif
1572
1573                         ELSE
1574!
1575!--                         First send the local index limits to PE0
1576                            ind(1) = nxl; ind(2) = nxr
1577                            ind(3) = nys; ind(4) = nyn
1578                            CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0,       &
1579                                           comm2d, ierr )
1580!
1581!--                         Send data to PE0
1582                            CALL MPI_SEND( local_2d(nxl,nys), ngp,             &
1583                                           MPI_REAL, 0, 1, comm2d, ierr )
1584                         ENDIF
1585!
1586!--                      A barrier has to be set, because otherwise some PEs may
1587!--                      proceed too fast so that PE0 may receive wrong data on
1588!--                      tag 0
1589                         CALL MPI_BARRIER( comm2d, ierr )
1590                      ENDIF
1591
1592                   ENDIF
1593#else
1594#if defined( __netcdf )
1595                   IF ( two_d ) THEN
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, 1, do2d_xy_time_count(av) /), &
1600                                           count = (/ nx+1, ny+1, 1, 1 /) )
1601                   ELSE
1602                      nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
1603                                              id_var_do2d(av,ivar),           &
1604                                              local_2d(nxl:nxr,nys:nyn),    &
1605                            start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
1606                                           count = (/ nx+1, ny+1, 1, 1 /) )
1607                   ENDIF
1608                   CALL netcdf_handle_error( 'data_output_2d', 447 )
1609#endif
1610#endif
1611
1612!
1613!--                For 2D-arrays (e.g. u*) only one cross-section is available.
1614!--                Hence exit loop of output levels.
1615                   IF ( two_d )  THEN
1616                      IF ( netcdf_data_format < 5 )  two_d = .FALSE.
1617                      EXIT loop1
1618                   ENDIF
1619
1620                CASE ( 'xz' )
1621!
1622!--                Update the netCDF xz cross section time axis.
1623!--                In case of parallel output, this is only done by PE0
1624!--                to increase the performance.
1625                   IF ( simulated_time /= do2d_xz_last_time(av) )  THEN
1626                      do2d_xz_time_count(av) = do2d_xz_time_count(av) + 1
1627                      do2d_xz_last_time(av)  = simulated_time
1628                      IF ( myid == 0 )  THEN
1629                         IF ( .NOT. data_output_2d_on_each_pe  &
1630                              .OR.  netcdf_data_format > 4 )   &
1631                         THEN
1632#if defined( __netcdf )
1633                            nc_stat = NF90_PUT_VAR( id_set_xz(av),             &
1634                                                    id_var_time_xz(av),        &
1635                                             (/ time_since_reference_point /), &
1636                                         start = (/ do2d_xz_time_count(av) /), &
1637                                                    count = (/ 1 /) )
1638                            CALL netcdf_handle_error( 'data_output_2d', 56 )
1639#endif
1640                         ENDIF
1641                      ENDIF
1642                   ENDIF
1643
1644!
1645!--                If required, carry out averaging along y
1646                   IF ( section(is,s_ind) == -1 )  THEN
1647
1648                      ALLOCATE( local_2d_l(nxl:nxr,nzb_do:nzt_do) )
1649                      local_2d_l = 0.0_wp
1650                      ngp = ( nxr-nxl + 1 ) * ( nzt_do-nzb_do + 1 )
1651!
1652!--                   First local averaging on the PE
1653                      DO  k = nzb_do, nzt_do
1654                         DO  j = nys, nyn
1655                            DO  i = nxl, nxr
1656                               local_2d_l(i,k) = local_2d_l(i,k) +             &
1657                                                 local_pf(i,j,k)
1658                            ENDDO
1659                         ENDDO
1660                      ENDDO
1661#if defined( __parallel )
1662!
1663!--                   Now do the averaging over all PEs along y
1664                      IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1665                      CALL MPI_ALLREDUCE( local_2d_l(nxl,nzb_do),                &
1666                                          local_2d(nxl,nzb_do), ngp, MPI_REAL,   &
1667                                          MPI_SUM, comm1dy, ierr )
1668#else
1669                      local_2d = local_2d_l
1670#endif
1671                      local_2d = local_2d / ( ny + 1.0_wp )
1672
1673                      DEALLOCATE( local_2d_l )
1674
1675                   ELSE
1676!
1677!--                   Just store the respective section on the local array
1678!--                   (but only if it is available on this PE!)
1679                      IF ( section(is,s_ind) >= nys  .AND.  section(is,s_ind) <= nyn ) &
1680                      THEN
1681                         local_2d = local_pf(:,section(is,s_ind),nzb_do:nzt_do)
1682                      ENDIF
1683
1684                   ENDIF
1685
1686#if defined( __parallel )
1687                   IF ( netcdf_data_format > 4 )  THEN
1688!
1689!--                   Output in netCDF4/HDF5 format.
1690!--                   Output only on those PEs where the respective cross
1691!--                   sections reside. Cross sections averaged along y are
1692!--                   output on the respective first PE along y (myidy=0).
1693                      IF ( ( section(is,s_ind) >= nys  .AND.                   &
1694                             section(is,s_ind) <= nyn )  .OR.                  &
1695                           ( section(is,s_ind) == -1  .AND.  myidy == 0 ) )  THEN
1696#if defined( __netcdf )
1697!
1698!--                      For parallel output, all cross sections are first
1699!--                      stored here on a local array and will be written to the
1700!--                      output file afterwards to increase the performance.
1701                         DO  i = nxl, nxr
1702                            DO  k = nzb_do, nzt_do
1703                               local_2d_sections_l(i,is,k) = local_2d(i,k)
1704                            ENDDO
1705                         ENDDO
1706#endif
1707                      ENDIF
1708
1709                   ELSE
1710
1711                      IF ( data_output_2d_on_each_pe )  THEN
1712!
1713!--                      Output of partial arrays on each PE. If the cross
1714!--                      section does not reside on the PE, output special
1715!--                      index values.
1716#if defined( __netcdf )
1717                         IF ( myid == 0 )  THEN
1718                            WRITE ( 22 )  time_since_reference_point,          &
1719                                          do2d_xz_time_count(av), av
1720                         ENDIF
1721#endif
1722                         DO  i = 0, io_blocks-1
1723                            IF ( i == io_group )  THEN
1724                               IF ( ( section(is,s_ind) >= nys  .AND.          &
1725                                      section(is,s_ind) <= nyn )  .OR.         &
1726                                    ( section(is,s_ind) == -1  .AND.           &
1727                                      nys-1 == -1 ) )                          &
1728                               THEN
1729                                  WRITE (22)  nxl, nxr, nzb_do, nzt_do, nzb, nzt+1
1730                                  WRITE (22)  local_2d
1731                               ELSE
1732                                  WRITE (22)  -1, -1, -1, -1, -1, -1
1733                               ENDIF
1734                            ENDIF
1735#if defined( __parallel )
1736                            CALL MPI_BARRIER( comm2d, ierr )
1737#endif
1738                         ENDDO
1739
1740                      ELSE
1741!
1742!--                      PE0 receives partial arrays from all processors of the
1743!--                      respective cross section and outputs them. Here a
1744!--                      barrier has to be set, because otherwise
1745!--                      "-MPI- FATAL: Remote protocol queue full" may occur.
1746                         CALL MPI_BARRIER( comm2d, ierr )
1747
1748                         ngp = ( nxr-nxl + 1 ) * ( nzt_do-nzb_do + 1 )
1749                         IF ( myid == 0 )  THEN
1750!
1751!--                         Local array can be relocated directly.
1752                            IF ( ( section(is,s_ind) >= nys  .AND.              &
1753                                   section(is,s_ind) <= nyn )  .OR.             &
1754                                 ( section(is,s_ind) == -1  .AND.               &
1755                                   nys-1 == -1 ) )  THEN
1756                               total_2d(nxl:nxr,nzb_do:nzt_do) = local_2d
1757                            ENDIF
1758!
1759!--                         Receive data from all other PEs.
1760                            DO  n = 1, numprocs-1
1761!
1762!--                            Receive index limits first, then array.
1763!--                            Index limits are received in arbitrary order from
1764!--                            the PEs.
1765                               CALL MPI_RECV( ind(1), 4, MPI_INTEGER,          &
1766                                              MPI_ANY_SOURCE, 0, comm2d,       &
1767                                              status, ierr )
1768!
1769!--                            Not all PEs have data for XZ-cross-section.
1770                               IF ( ind(1) /= -9999 )  THEN
1771                                  sender = status(MPI_SOURCE)
1772                                  DEALLOCATE( local_2d )
1773                                  ALLOCATE( local_2d(ind(1):ind(2),            &
1774                                                     ind(3):ind(4)) )
1775                                  CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1776                                                 MPI_REAL, sender, 1, comm2d,  &
1777                                                 status, ierr )
1778                                  total_2d(ind(1):ind(2),ind(3):ind(4)) =      &
1779                                                                        local_2d
1780                               ENDIF
1781                            ENDDO
1782!
1783!--                         Relocate the local array for the next loop increment
1784                            DEALLOCATE( local_2d )
1785                            ALLOCATE( local_2d(nxl:nxr,nzb_do:nzt_do) )
1786
1787#if defined( __netcdf )
1788                            nc_stat = NF90_PUT_VAR( id_set_xz(av),             &
1789                                                 id_var_do2d(av,ivar),           &
1790                                                 total_2d(0:nx,nzb_do:nzt_do), &
1791                               start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
1792                                          count = (/ nx+1, 1, nzt_do-nzb_do+1, 1 /) )
1793                            CALL netcdf_handle_error( 'data_output_2d', 58 )
1794#endif
1795
1796                         ELSE
1797!
1798!--                         If the cross section resides on the PE, send the
1799!--                         local index limits, otherwise send -9999 to PE0.
1800                            IF ( ( section(is,s_ind) >= nys  .AND.              &
1801                                   section(is,s_ind) <= nyn )  .OR.             &
1802                                 ( section(is,s_ind) == -1  .AND.  nys-1 == -1 ) ) &
1803                            THEN
1804                               ind(1) = nxl; ind(2) = nxr
1805                               ind(3) = nzb_do;   ind(4) = nzt_do
1806                            ELSE
1807                               ind(1) = -9999; ind(2) = -9999
1808                               ind(3) = -9999; ind(4) = -9999
1809                            ENDIF
1810                            CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0,       &
1811                                           comm2d, ierr )
1812!
1813!--                         If applicable, send data to PE0.
1814                            IF ( ind(1) /= -9999 )  THEN
1815                               CALL MPI_SEND( local_2d(nxl,nzb_do), ngp,         &
1816                                              MPI_REAL, 0, 1, comm2d, ierr )
1817                            ENDIF
1818                         ENDIF
1819!
1820!--                      A barrier has to be set, because otherwise some PEs may
1821!--                      proceed too fast so that PE0 may receive wrong data on
1822!--                      tag 0
1823                         CALL MPI_BARRIER( comm2d, ierr )
1824                      ENDIF
1825
1826                   ENDIF
1827#else
1828#if defined( __netcdf )
1829                   nc_stat = NF90_PUT_VAR( id_set_xz(av),                   &
1830                                           id_var_do2d(av,ivar),              &
1831                                           local_2d(nxl:nxr,nzb_do:nzt_do), &
1832                            start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
1833                                       count = (/ nx+1, 1, nzt_do-nzb_do+1, 1 /) )
1834                   CALL netcdf_handle_error( 'data_output_2d', 451 )
1835#endif
1836#endif
1837
1838                CASE ( 'yz' )
1839!
1840!--                Update the netCDF yz cross section time axis.
1841!--                In case of parallel output, this is only done by PE0
1842!--                to increase the performance.
1843                   IF ( simulated_time /= do2d_yz_last_time(av) )  THEN
1844                      do2d_yz_time_count(av) = do2d_yz_time_count(av) + 1
1845                      do2d_yz_last_time(av)  = simulated_time
1846                      IF ( myid == 0 )  THEN
1847                         IF ( .NOT. data_output_2d_on_each_pe  &
1848                              .OR.  netcdf_data_format > 4 )   &
1849                         THEN
1850#if defined( __netcdf )
1851                            nc_stat = NF90_PUT_VAR( id_set_yz(av),             &
1852                                                    id_var_time_yz(av),        &
1853                                             (/ time_since_reference_point /), &
1854                                         start = (/ do2d_yz_time_count(av) /), &
1855                                                    count = (/ 1 /) )
1856                            CALL netcdf_handle_error( 'data_output_2d', 59 )
1857#endif
1858                         ENDIF
1859                      ENDIF
1860                   ENDIF
1861
1862!
1863!--                If required, carry out averaging along x
1864                   IF ( section(is,s_ind) == -1 )  THEN
1865
1866                      ALLOCATE( local_2d_l(nys:nyn,nzb_do:nzt_do) )
1867                      local_2d_l = 0.0_wp
1868                      ngp = ( nyn-nys+1 ) * ( nzt_do-nzb_do+1 )
1869!
1870!--                   First local averaging on the PE
1871                      DO  k = nzb_do, nzt_do
1872                         DO  j = nys, nyn
1873                            DO  i = nxl, nxr
1874                               local_2d_l(j,k) = local_2d_l(j,k) +             &
1875                                                 local_pf(i,j,k)
1876                            ENDDO
1877                         ENDDO
1878                      ENDDO
1879#if defined( __parallel )
1880!
1881!--                   Now do the averaging over all PEs along x
1882                      IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
1883                      CALL MPI_ALLREDUCE( local_2d_l(nys,nzb_do),                &
1884                                          local_2d(nys,nzb_do), ngp, MPI_REAL,   &
1885                                          MPI_SUM, comm1dx, ierr )
1886#else
1887                      local_2d = local_2d_l
1888#endif
1889                      local_2d = local_2d / ( nx + 1.0_wp )
1890
1891                      DEALLOCATE( local_2d_l )
1892
1893                   ELSE
1894!
1895!--                   Just store the respective section on the local array
1896!--                   (but only if it is available on this PE!)
1897                      IF ( section(is,s_ind) >= nxl  .AND.  section(is,s_ind) <= nxr ) &
1898                      THEN
1899                         local_2d = local_pf(section(is,s_ind),:,nzb_do:nzt_do)
1900                      ENDIF
1901
1902                   ENDIF
1903
1904#if defined( __parallel )
1905                   IF ( netcdf_data_format > 4 )  THEN
1906!
1907!--                   Output in netCDF4/HDF5 format.
1908!--                   Output only on those PEs where the respective cross
1909!--                   sections reside. Cross sections averaged along x are
1910!--                   output on the respective first PE along x (myidx=0).
1911                      IF ( ( section(is,s_ind) >= nxl  .AND.                       &
1912                             section(is,s_ind) <= nxr )  .OR.                      &
1913                           ( section(is,s_ind) == -1  .AND.  myidx == 0 ) )  THEN
1914#if defined( __netcdf )
1915!
1916!--                      For parallel output, all cross sections are first
1917!--                      stored here on a local array and will be written to the
1918!--                      output file afterwards to increase the performance.
1919                         DO  j = nys, nyn
1920                            DO  k = nzb_do, nzt_do
1921                               local_2d_sections_l(is,j,k) = local_2d(j,k)
1922                            ENDDO
1923                         ENDDO
1924#endif
1925                      ENDIF
1926
1927                   ELSE
1928
1929                      IF ( data_output_2d_on_each_pe )  THEN
1930!
1931!--                      Output of partial arrays on each PE. If the cross
1932!--                      section does not reside on the PE, output special
1933!--                      index values.
1934#if defined( __netcdf )
1935                         IF ( myid == 0 )  THEN
1936                            WRITE ( 23 )  time_since_reference_point,          &
1937                                          do2d_yz_time_count(av), av
1938                         ENDIF
1939#endif
1940                         DO  i = 0, io_blocks-1
1941                            IF ( i == io_group )  THEN
1942                               IF ( ( section(is,s_ind) >= nxl  .AND.          &
1943                                      section(is,s_ind) <= nxr )  .OR.         &
1944                                    ( section(is,s_ind) == -1  .AND.           &
1945                                      nxl-1 == -1 ) )                          &
1946                               THEN
1947                                  WRITE (23)  nys, nyn, nzb_do, nzt_do, nzb, nzt+1
1948                                  WRITE (23)  local_2d
1949                               ELSE
1950                                  WRITE (23)  -1, -1, -1, -1, -1, -1
1951                               ENDIF
1952                            ENDIF
1953#if defined( __parallel )
1954                            CALL MPI_BARRIER( comm2d, ierr )
1955#endif
1956                         ENDDO
1957
1958                      ELSE
1959!
1960!--                      PE0 receives partial arrays from all processors of the
1961!--                      respective cross section and outputs them. Here a
1962!--                      barrier has to be set, because otherwise
1963!--                      "-MPI- FATAL: Remote protocol queue full" may occur.
1964                         CALL MPI_BARRIER( comm2d, ierr )
1965
1966                         ngp = ( nyn-nys+1 ) * ( nzt_do-nzb_do+1 )
1967                         IF ( myid == 0 )  THEN
1968!
1969!--                         Local array can be relocated directly.
1970                            IF ( ( section(is,s_ind) >= nxl  .AND.             &
1971                                   section(is,s_ind) <= nxr )   .OR.           &
1972                                 ( section(is,s_ind) == -1  .AND.  nxl-1 == -1 ) ) &
1973                            THEN
1974                               total_2d(nys:nyn,nzb_do:nzt_do) = local_2d
1975                            ENDIF
1976!
1977!--                         Receive data from all other PEs.
1978                            DO  n = 1, numprocs-1
1979!
1980!--                            Receive index limits first, then array.
1981!--                            Index limits are received in arbitrary order from
1982!--                            the PEs.
1983                               CALL MPI_RECV( ind(1), 4, MPI_INTEGER,          &
1984                                              MPI_ANY_SOURCE, 0, comm2d,       &
1985                                              status, ierr )
1986!
1987!--                            Not all PEs have data for YZ-cross-section.
1988                               IF ( ind(1) /= -9999 )  THEN
1989                                  sender = status(MPI_SOURCE)
1990                                  DEALLOCATE( local_2d )
1991                                  ALLOCATE( local_2d(ind(1):ind(2),            &
1992                                                     ind(3):ind(4)) )
1993                                  CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1994                                                 MPI_REAL, sender, 1, comm2d,  &
1995                                                 status, ierr )
1996                                  total_2d(ind(1):ind(2),ind(3):ind(4)) =      &
1997                                                                        local_2d
1998                               ENDIF
1999                            ENDDO
2000!
2001!--                         Relocate the local array for the next loop increment
2002                            DEALLOCATE( local_2d )
2003                            ALLOCATE( local_2d(nys:nyn,nzb_do:nzt_do) )
2004
2005#if defined( __netcdf )
2006                            nc_stat = NF90_PUT_VAR( id_set_yz(av),             &
2007                                                 id_var_do2d(av,ivar),           &
2008                                                 total_2d(0:ny,nzb_do:nzt_do), &
2009                            start = (/ is, 1, 1, do2d_yz_time_count(av) /),    &
2010                                       count = (/ 1, ny+1, nzt_do-nzb_do+1, 1 /) )
2011                            CALL netcdf_handle_error( 'data_output_2d', 61 )
2012#endif
2013
2014                         ELSE
2015!
2016!--                         If the cross section resides on the PE, send the
2017!--                         local index limits, otherwise send -9999 to PE0.
2018                            IF ( ( section(is,s_ind) >= nxl  .AND.              &
2019                                   section(is,s_ind) <= nxr )  .OR.             &
2020                                 ( section(is,s_ind) == -1  .AND.  nxl-1 == -1 ) ) &
2021                            THEN
2022                               ind(1) = nys; ind(2) = nyn
2023                               ind(3) = nzb_do;   ind(4) = nzt_do
2024                            ELSE
2025                               ind(1) = -9999; ind(2) = -9999
2026                               ind(3) = -9999; ind(4) = -9999
2027                            ENDIF
2028                            CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0,       &
2029                                           comm2d, ierr )
2030!
2031!--                         If applicable, send data to PE0.
2032                            IF ( ind(1) /= -9999 )  THEN
2033                               CALL MPI_SEND( local_2d(nys,nzb_do), ngp,         &
2034                                              MPI_REAL, 0, 1, comm2d, ierr )
2035                            ENDIF
2036                         ENDIF
2037!
2038!--                      A barrier has to be set, because otherwise some PEs may
2039!--                      proceed too fast so that PE0 may receive wrong data on
2040!--                      tag 0
2041                         CALL MPI_BARRIER( comm2d, ierr )
2042                      ENDIF
2043
2044                   ENDIF
2045#else
2046#if defined( __netcdf )
2047                   nc_stat = NF90_PUT_VAR( id_set_yz(av),                   &
2048                                           id_var_do2d(av,ivar),              &
2049                                           local_2d(nys:nyn,nzb_do:nzt_do), &
2050                            start = (/ is, 1, 1, do2d_xz_time_count(av) /), &
2051                                           count = (/ 1, ny+1, nzt_do-nzb_do+1, 1 /) )
2052                   CALL netcdf_handle_error( 'data_output_2d', 452 )
2053#endif
2054#endif
2055
2056             END SELECT
2057
2058             is = is + 1
2059          ENDDO loop1
2060
2061!
2062!--       For parallel output, all data were collected before on a local array
2063!--       and are written now to the netcdf file. This must be done to increase
2064!--       the performance of the parallel output.
2065#if defined( __netcdf )
2066          IF ( netcdf_data_format > 4 )  THEN
2067
2068                SELECT CASE ( mode )
2069
2070                   CASE ( 'xy' )
2071                      IF ( two_d ) THEN
2072                         nis = 1
2073                         two_d = .FALSE.
2074                      ELSE
2075                         nis = ns
2076                      ENDIF
2077!
2078!--                   Do not output redundant ghost point data except for the
2079!--                   boundaries of the total domain.
2080!                      IF ( nxr == nx  .AND.  nyn /= ny )  THEN
2081!                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
2082!                                                 id_var_do2d(av,ivar),           &
2083!                                                 local_2d_sections(nxl:nxr+1,  &
2084!                                                    nys:nyn,1:nis),            &
2085!                                                 start = (/ nxl+1, nys+1, 1,   &
2086!                                                    do2d_xy_time_count(av) /), &
2087!                                                 count = (/ nxr-nxl+2,         &
2088!                                                            nyn-nys+1, nis, 1  &
2089!                                                          /) )
2090!                      ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
2091!                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
2092!                                                 id_var_do2d(av,ivar),           &
2093!                                                 local_2d_sections(nxl:nxr,    &
2094!                                                    nys:nyn+1,1:nis),          &
2095!                                                 start = (/ nxl+1, nys+1, 1,   &
2096!                                                    do2d_xy_time_count(av) /), &
2097!                                                 count = (/ nxr-nxl+1,         &
2098!                                                            nyn-nys+2, nis, 1  &
2099!                                                          /) )
2100!                      ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
2101!                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
2102!                                                 id_var_do2d(av,ivar),           &
2103!                                                 local_2d_sections(nxl:nxr+1,  &
2104!                                                    nys:nyn+1,1:nis),          &
2105!                                                 start = (/ nxl+1, nys+1, 1,   &
2106!                                                    do2d_xy_time_count(av) /), &
2107!                                                 count = (/ nxr-nxl+2,         &
2108!                                                            nyn-nys+2, nis, 1  &
2109!                                                          /) )
2110!                      ELSE
2111                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
2112                                                 id_var_do2d(av,ivar),           &
2113                                                 local_2d_sections(nxl:nxr,    &
2114                                                    nys:nyn,1:nis),            &
2115                                                 start = (/ nxl+1, nys+1, 1,   &
2116                                                    do2d_xy_time_count(av) /), &
2117                                                 count = (/ nxr-nxl+1,         &
2118                                                            nyn-nys+1, nis, 1  &
2119                                                          /) )
2120!                      ENDIF   
2121
2122                      CALL netcdf_handle_error( 'data_output_2d', 55 )
2123
2124                   CASE ( 'xz' )
2125!
2126!--                   First, all PEs get the information of all cross-sections.
2127!--                   Then the data are written to the output file by all PEs
2128!--                   while NF90_COLLECTIVE is set in subroutine
2129!--                   define_netcdf_header. Although redundant information are
2130!--                   written to the output file in that case, the performance
2131!--                   is significantly better compared to the case where only
2132!--                   the first row of PEs in x-direction (myidx = 0) is given
2133!--                   the output while NF90_INDEPENDENT is set.
2134                      IF ( npey /= 1 )  THEN
2135                         
2136#if defined( __parallel )
2137!
2138!--                      Distribute data over all PEs along y
2139                         ngp = ( nxr-nxl+1 ) * ( nzt_do-nzb_do+1 ) * ns
2140                         IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2141                         CALL MPI_ALLREDUCE( local_2d_sections_l(nxl,1,nzb_do),  &
2142                                             local_2d_sections(nxl,1,nzb_do),    &
2143                                             ngp, MPI_REAL, MPI_SUM, comm1dy,  &
2144                                             ierr )
2145#else
2146                         local_2d_sections = local_2d_sections_l
2147#endif
2148                      ENDIF
2149!
2150!--                   Do not output redundant ghost point data except for the
2151!--                   boundaries of the total domain.
2152!                      IF ( nxr == nx )  THEN
2153!                         nc_stat = NF90_PUT_VAR( id_set_xz(av),                &
2154!                                             id_var_do2d(av,ivar),               &
2155!                                             local_2d_sections(nxl:nxr+1,1:ns, &
2156!                                                nzb_do:nzt_do),                &
2157!                                             start = (/ nxl+1, 1, 1,           &
2158!                                                do2d_xz_time_count(av) /),     &
2159!                                             count = (/ nxr-nxl+2, ns, nzt_do-nzb_do+1,  &
2160!                                                        1 /) )
2161!                      ELSE
2162                         nc_stat = NF90_PUT_VAR( id_set_xz(av),                &
2163                                             id_var_do2d(av,ivar),               &
2164                                             local_2d_sections(nxl:nxr,1:ns,   &
2165                                                nzb_do:nzt_do),                &
2166                                             start = (/ nxl+1, 1, 1,           &
2167                                                do2d_xz_time_count(av) /),     &
2168                                             count = (/ nxr-nxl+1, ns, nzt_do-nzb_do+1,  &
2169                                                1 /) )
2170!                      ENDIF
2171
2172                      CALL netcdf_handle_error( 'data_output_2d', 57 )
2173
2174                   CASE ( 'yz' )
2175!
2176!--                   First, all PEs get the information of all cross-sections.
2177!--                   Then the data are written to the output file by all PEs
2178!--                   while NF90_COLLECTIVE is set in subroutine
2179!--                   define_netcdf_header. Although redundant information are
2180!--                   written to the output file in that case, the performance
2181!--                   is significantly better compared to the case where only
2182!--                   the first row of PEs in y-direction (myidy = 0) is given
2183!--                   the output while NF90_INDEPENDENT is set.
2184                      IF ( npex /= 1 )  THEN
2185
2186#if defined( __parallel )
2187!
2188!--                      Distribute data over all PEs along x
2189                         ngp = ( nyn-nys+1 ) * ( nzt-nzb + 2 ) * ns
2190                         IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
2191                         CALL MPI_ALLREDUCE( local_2d_sections_l(1,nys,nzb_do),  &
2192                                             local_2d_sections(1,nys,nzb_do),    &
2193                                             ngp, MPI_REAL, MPI_SUM, comm1dx,  &
2194                                             ierr )
2195#else
2196                         local_2d_sections = local_2d_sections_l
2197#endif
2198                      ENDIF
2199!
2200!--                   Do not output redundant ghost point data except for the
2201!--                   boundaries of the total domain.
2202!                      IF ( nyn == ny )  THEN
2203!                         nc_stat = NF90_PUT_VAR( id_set_yz(av),                &
2204!                                             id_var_do2d(av,ivar),               &
2205!                                             local_2d_sections(1:ns,           &
2206!                                                nys:nyn+1,nzb_do:nzt_do),      &
2207!                                             start = (/ 1, nys+1, 1,           &
2208!                                                do2d_yz_time_count(av) /),     &
2209!                                             count = (/ ns, nyn-nys+2,         &
2210!                                                        nzt_do-nzb_do+1, 1 /) )
2211!                      ELSE
2212                         nc_stat = NF90_PUT_VAR( id_set_yz(av),                &
2213                                             id_var_do2d(av,ivar),               &
2214                                             local_2d_sections(1:ns,nys:nyn,   &
2215                                                nzb_do:nzt_do),                &
2216                                             start = (/ 1, nys+1, 1,           &
2217                                                do2d_yz_time_count(av) /),     &
2218                                             count = (/ ns, nyn-nys+1,         &
2219                                                        nzt_do-nzb_do+1, 1 /) )
2220!                      ENDIF
2221
2222                      CALL netcdf_handle_error( 'data_output_2d', 60 )
2223
2224                   CASE DEFAULT
2225                      message_string = 'unknown cross-section: ' // TRIM( mode )
2226                      CALL message( 'data_output_2d', 'PA0180', 1, 2, 0, 6, 0 )
2227
2228                END SELECT                     
2229
2230          ENDIF
2231#endif
2232       ENDIF
2233
2234       ivar = ivar + 1
2235       l = MAX( 2, LEN_TRIM( do2d(av,ivar) ) )
2236       do2d_mode = do2d(av,ivar)(l-1:l)
2237
2238    ENDDO
2239
2240!
2241!-- Deallocate temporary arrays.
2242    IF ( ALLOCATED( level_z ) )  DEALLOCATE( level_z )
2243    IF ( netcdf_data_format > 4 )  THEN
2244       DEALLOCATE( local_pf, local_2d, local_2d_sections )
2245       IF( mode == 'xz' .OR. mode == 'yz' ) DEALLOCATE( local_2d_sections_l )
2246    ENDIF
2247#if defined( __parallel )
2248    IF ( .NOT.  data_output_2d_on_each_pe  .AND.  myid == 0 )  THEN
2249       DEALLOCATE( total_2d )
2250    ENDIF
2251#endif
2252
2253!
2254!-- Close plot output file.
2255    file_id = 20 + s_ind
2256
2257    IF ( data_output_2d_on_each_pe )  THEN
2258       DO  i = 0, io_blocks-1
2259          IF ( i == io_group )  THEN
2260             CALL close_file( file_id )
2261          ENDIF
2262#if defined( __parallel )
2263          CALL MPI_BARRIER( comm2d, ierr )
2264#endif
2265       ENDDO
2266    ELSE
2267       IF ( myid == 0 )  CALL close_file( file_id )
2268    ENDIF
2269
2270    CALL cpu_log( log_point(3), 'data_output_2d', 'stop' )
2271
2272 END SUBROUTINE data_output_2d
Note: See TracBrowser for help on using the repository browser.