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

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

further modularization of land surface model (2D/3D output and restart data). Bugfix for restart runs without land surface model

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