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

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

Diagnostic output: Define arrays over ghost points in order to allow for standard mpi-io treatment. By this modularization of restart-data input is possible with the module interface. Move input of restart data to doq_rrd_local. Enable mpi-io for restart data. Bugfix: add missing restart input of wtheta_av, wq_av, wu_av, and wv_av.

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