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

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