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

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

Bugfix, in order to steer user-defined output, setting flag found explicitly

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