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

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

minor bugfixes for radiation model. bugfix in subjob

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