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

Last change on this file since 2735 was 2735, checked in by suehring, 6 years ago

Output of resistance also urban-type surfaces

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