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

Last change on this file since 3467 was 3467, checked in by suehring, 5 years ago

Branch salsa @3446 re-integrated into trunk

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