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

Last change on this file since 2422 was 2292, checked in by schwenkel, 7 years ago

implementation of new bulk microphysics scheme

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