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

Last change on this file since 3534 was 3525, checked in by kanani, 5 years ago

Changes related to clean-up of biometeorology (by dom_dwd_user)

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