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

Last change on this file since 2718 was 2718, checked in by maronga, 4 years ago

deleting of deprecated files; headers updated where needed

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