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

Last change on this file since 3294 was 3294, checked in by raasch, 6 years ago

modularization of the ocean code

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