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

Last change on this file since 3918 was 3885, checked in by kanani, 6 years ago

restructure/add location/debug messages

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