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

Last change on this file since 3448 was 3448, checked in by kanani, 5 years ago

Implementation of human thermal indices (from branch biomet_p2 at r3444)

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