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

Last change on this file since 3298 was 3298, checked in by kanani, 5 years ago

Merge chemistry branch at r3297 to trunk

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