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

Last change on this file since 1586 was 1586, checked in by maronga, 9 years ago

last commit documented

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