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

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

last commit documented

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