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

Last change on this file since 3574 was 3569, checked in by kanani, 6 years ago

Fix for biomet output (ticket:757), merge of uv_exposure into biometeorology_mod

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