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

Last change on this file since 2174 was 2101, checked in by suehring, 8 years ago

last commit documented

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