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

Last change on this file since 2750 was 2743, checked in by suehring, 7 years ago

In case of natural- and urban-type surfaces output surfaces fluxes in W/m2

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