source: palm/trunk/SOURCE/data_output_2d.f90

Last change on this file was 4828, checked in by Giersch, 3 years ago

Copyright updated to year 2021, interface pmc_sort removed to accelarate the nesting code

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