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

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

Output of surface temperature tsurf* at urban- and natural-type surfaces

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