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

Last change on this file since 2963 was 2963, checked in by suehring, 6 years ago

Minor revision of static input file checks, bugfix in initialization of surface-fractions in LSM; minor bugfix in initialization of albedo at window-surfaces; for clearer access of albedo and emissivity introduce index for vegetation/wall, pavement/green-wall and water/window surfaces

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