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

Last change on this file since 3943 was 3943, checked in by maronga, 5 years ago

bugfixes in urban surface model; output of greenz roof transpiration added/corrected; minor formatting improvements

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