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

Last change on this file since 2817 was 2817, checked in by knoop, 4 years ago

Preliminary gust module interface implemented

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