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

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

Bugfix initialization of %pt_surface array; Output of surface temperature also for default-type surfaces

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