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

Last change on this file since 3554 was 3554, checked in by gronemeier, 5 years ago

renamed variable if to ivar; add variable description

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