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

Last change on this file since 1784 was 1784, checked in by raasch, 8 years ago

last commit documented

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