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

Last change on this file since 3606 was 3597, checked in by maronga, 5 years ago

revised calculation of near surface air potential temperature

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