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

Last change on this file since 2742 was 2742, checked in by suehring, 6 years ago

Output of surface temperature tsurf* at urban- and natural-type surfaces

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