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

Last change on this file since 3014 was 3014, checked in by maronga, 6 years ago

series of bugfixes

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