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

Last change on this file since 1972 was 1972, checked in by maronga, 8 years ago

further modularization of land surface model (2D/3D output and restart data)

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