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

Last change on this file since 1992 was 1981, checked in by suehring, 8 years ago

last commit documented

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