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

Last change on this file since 3677 was 3655, checked in by knoop, 6 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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