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

Last change on this file since 3262 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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