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

Last change on this file since 3045 was 3045, checked in by Giersch, 6 years ago

Code adjusted according to coding standards, renamed namelists, error messages revised until PA0347, output CASE 108 disabled

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