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

Last change on this file since 3643 was 3637, checked in by knoop, 6 years ago

M Makefile

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