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

Last change on this file since 3004 was 3004, checked in by Giersch, 6 years ago

precipitation_rate removed, further allocation checks for data output of averaged quantities implemented, double CALL of flow_statistics at the beginning of time_integration removed, further minor bugfixes, comments added

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