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

Last change on this file since 2797 was 2797, checked in by suehring, 4 years ago

Output of ground-heat flux at natural- and urban-type surfaces in one output variable; enable restart data of _av variables that belong to both land- and urban-surface model

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