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

Last change on this file since 2805 was 2805, checked in by suehring, 4 years ago

Bugfix in re-mapping surface-element data in case of restarts; bugfix in initialization of water-surface roughness; bugfix - uninitialized resistance at urban-type surfaces

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