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

Last change on this file since 2191 was 2191, checked in by raasch, 7 years ago

last commit documented

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