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

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

LSM output of r_a and r_s added

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