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

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

Merge chemistry branch at r3297 to trunk

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