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

Last change on this file since 1691 was 1691, checked in by maronga, 6 years ago

various bugfixes and modifications of the atmosphere-land-surface-radiation interaction. Completely re-written routine to calculate surface fluxes (surface_layer_fluxes.f90) that replaces prandtl_fluxes. Minor formatting corrections and renamings

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