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

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

Major bugfix in calculation of ol and ts at building roofs; bugfix in restart data for surface elements; in 2D data output, mask latent heat flux at urban-type surfaces

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