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

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

last commit documented

  • Property svn:keywords set to Id
File size: 82.6 KB
RevLine 
[1]1 SUBROUTINE data_output_2d( mode, av )
2
[1036]3!--------------------------------------------------------------------------------!
4! This file is part of PALM.
5!
6! PALM is free software: you can redistribute it and/or modify it under the terms
7! of the GNU General Public License as published by the Free Software Foundation,
8! either version 3 of the License, or (at your option) any later version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
[1310]17! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]18!--------------------------------------------------------------------------------!
19!
[254]20! Current revisions:
[1]21! -----------------
[1552]22!
23!
24! Former revisions:
25! -----------------
26! $Id: data_output_2d.f90 1552 2015-03-03 14:27:15Z 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,     &
141               qsws_veg_eb, qsws_veg_eb_av, shf_eb, shf_eb_av, t_soil,         &
142               t_soil_av, zs
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
[1]937             CASE ( 's_xy', 's_xz', 's_yz' )
938                IF ( av == 0 )  THEN
939                   to_be_resorted => q
940                ELSE
[355]941                   to_be_resorted => s_av
[1]942                ENDIF
943
[96]944             CASE ( 'sa_xy', 'sa_xz', 'sa_yz' )
945                IF ( av == 0 )  THEN
946                   to_be_resorted => sa
947                ELSE
948                   to_be_resorted => sa_av
949                ENDIF
950
[354]951             CASE ( 'shf*_xy' )        ! 2d-array
952                IF ( av == 0 ) THEN
[667]953                   DO  i = nxlg, nxrg
954                      DO  j = nysg, nyng
[354]955                         local_pf(i,j,nzb+1) =  shf(j,i)
956                      ENDDO
957                   ENDDO
958                ELSE
[667]959                   DO  i = nxlg, nxrg
960                      DO  j = nysg, nyng
[354]961                         local_pf(i,j,nzb+1) =  shf_av(j,i)
962                      ENDDO
963                   ENDDO
964                ENDIF
965                resorted = .TRUE.
966                two_d = .TRUE.
967                level_z(nzb+1) = zu(nzb+1)
968
[1551]969             CASE ( 'shf_eb*_xy' )        ! 2d-array
970                IF ( av == 0 ) THEN
971                   DO  i = nxlg, nxrg
972                      DO  j = nysg, nyng
973                         local_pf(i,j,nzb+1) =  shf_eb(j,i)
974                      ENDDO
975                   ENDDO
976                ELSE
977                   DO  i = nxlg, nxrg
978                      DO  j = nysg, nyng
979                         local_pf(i,j,nzb+1) =  shf_eb_av(j,i)
980                      ENDDO
981                   ENDDO
982                ENDIF
983                resorted = .TRUE.
984                two_d = .TRUE.
985                level_z(nzb+1) = zu(nzb+1)
986
[1]987             CASE ( 't*_xy' )        ! 2d-array
988                IF ( av == 0 )  THEN
[667]989                   DO  i = nxlg, nxrg
990                      DO  j = nysg, nyng
[1]991                         local_pf(i,j,nzb+1) = ts(j,i)
992                      ENDDO
993                   ENDDO
994                ELSE
[667]995                   DO  i = nxlg, nxrg
996                      DO  j = nysg, nyng
[1]997                         local_pf(i,j,nzb+1) = ts_av(j,i)
998                      ENDDO
999                   ENDDO
1000                ENDIF
1001                resorted = .TRUE.
1002                two_d = .TRUE.
1003                level_z(nzb+1) = zu(nzb+1)
1004
[1551]1005             CASE ( 't_soil_xy', 't_soil_xz', 't_soil_yz' )
1006                nzb_do = nzb_soil
1007                nzt_do = nzt_soil
1008                IF ( av == 0 )  THEN
1009                   to_be_resorted => t_soil
1010                ELSE
1011                   to_be_resorted => t_soil_av
1012                ENDIF
1013                IF ( mode == 'xy' )  level_z = zs
1014
[1]1015             CASE ( 'u_xy', 'u_xz', 'u_yz' )
1016                IF ( av == 0 )  THEN
1017                   to_be_resorted => u
1018                ELSE
1019                   to_be_resorted => u_av
1020                ENDIF
1021                IF ( mode == 'xy' )  level_z = zu
1022!
1023!--             Substitute the values generated by "mirror" boundary condition
1024!--             at the bottom boundary by the real surface values.
1025                IF ( do2d(av,if) == 'u_xz'  .OR.  do2d(av,if) == 'u_yz' )  THEN
[1353]1026                   IF ( ibc_uv_b == 0 )  local_pf(:,:,nzb) = 0.0_wp
[1]1027                ENDIF
1028
1029             CASE ( 'u*_xy' )        ! 2d-array
1030                IF ( av == 0 )  THEN
[667]1031                   DO  i = nxlg, nxrg
1032                      DO  j = nysg, nyng
[1]1033                         local_pf(i,j,nzb+1) = us(j,i)
1034                      ENDDO
1035                   ENDDO
1036                ELSE
[667]1037                   DO  i = nxlg, nxrg
1038                      DO  j = nysg, nyng
[1]1039                         local_pf(i,j,nzb+1) = us_av(j,i)
1040                      ENDDO
1041                   ENDDO
1042                ENDIF
1043                resorted = .TRUE.
1044                two_d = .TRUE.
1045                level_z(nzb+1) = zu(nzb+1)
1046
1047             CASE ( 'v_xy', 'v_xz', 'v_yz' )
1048                IF ( av == 0 )  THEN
1049                   to_be_resorted => v
1050                ELSE
1051                   to_be_resorted => v_av
1052                ENDIF
1053                IF ( mode == 'xy' )  level_z = zu
1054!
1055!--             Substitute the values generated by "mirror" boundary condition
1056!--             at the bottom boundary by the real surface values.
1057                IF ( do2d(av,if) == 'v_xz'  .OR.  do2d(av,if) == 'v_yz' )  THEN
[1353]1058                   IF ( ibc_uv_b == 0 )  local_pf(:,:,nzb) = 0.0_wp
[1]1059                ENDIF
1060
1061             CASE ( 'vpt_xy', 'vpt_xz', 'vpt_yz' )
1062                IF ( av == 0 )  THEN
1063                   to_be_resorted => vpt
1064                ELSE
1065                   to_be_resorted => vpt_av
1066                ENDIF
1067                IF ( mode == 'xy' )  level_z = zu
1068
1069             CASE ( 'w_xy', 'w_xz', 'w_yz' )
1070                IF ( av == 0 )  THEN
1071                   to_be_resorted => w
1072                ELSE
1073                   to_be_resorted => w_av
1074                ENDIF
1075                IF ( mode == 'xy' )  level_z = zw
1076
[72]1077             CASE ( 'z0*_xy' )        ! 2d-array
1078                IF ( av == 0 ) THEN
[667]1079                   DO  i = nxlg, nxrg
1080                      DO  j = nysg, nyng
[72]1081                         local_pf(i,j,nzb+1) =  z0(j,i)
1082                      ENDDO
1083                   ENDDO
1084                ELSE
[667]1085                   DO  i = nxlg, nxrg
1086                      DO  j = nysg, nyng
[72]1087                         local_pf(i,j,nzb+1) =  z0_av(j,i)
1088                      ENDDO
1089                   ENDDO
1090                ENDIF
1091                resorted = .TRUE.
1092                two_d = .TRUE.
1093                level_z(nzb+1) = zu(nzb+1)
1094
[978]1095             CASE ( 'z0h*_xy' )        ! 2d-array
1096                IF ( av == 0 ) THEN
1097                   DO  i = nxlg, nxrg
1098                      DO  j = nysg, nyng
1099                         local_pf(i,j,nzb+1) =  z0h(j,i)
1100                      ENDDO
1101                   ENDDO
1102                ELSE
1103                   DO  i = nxlg, nxrg
1104                      DO  j = nysg, nyng
1105                         local_pf(i,j,nzb+1) =  z0h_av(j,i)
1106                      ENDDO
1107                   ENDDO
1108                ENDIF
1109                resorted = .TRUE.
1110                two_d = .TRUE.
1111                level_z(nzb+1) = zu(nzb+1)
1112
[1]1113             CASE DEFAULT
1114!
1115!--             User defined quantity
[1320]1116                CALL user_data_output_2d( av, do2d(av,if), found, grid,        &
[1551]1117                                          local_pf, two_d, nzb_do, nzt_do )
[1]1118                resorted = .TRUE.
1119
1120                IF ( grid == 'zu' )  THEN
1121                   IF ( mode == 'xy' )  level_z = zu
1122                ELSEIF ( grid == 'zw' )  THEN
1123                   IF ( mode == 'xy' )  level_z = zw
[343]1124                ELSEIF ( grid == 'zu1' ) THEN
1125                   IF ( mode == 'xy' )  level_z(nzb+1) = zu(nzb+1)
[1551]1126                ELSEIF ( grid == 'zs' ) THEN
1127                   IF ( mode == 'xy' )  level_z = zs
[1]1128                ENDIF
1129
1130                IF ( .NOT. found )  THEN
[1320]1131                   message_string = 'no output provided for: ' //              &
[274]1132                                    TRIM( do2d(av,if) )
[254]1133                   CALL message( 'data_output_2d', 'PA0181', 0, 0, 0, 6, 0 )
[1]1134                ENDIF
1135
1136          END SELECT
1137
1138!
1139!--       Resort the array to be output, if not done above
1140          IF ( .NOT. resorted )  THEN
[667]1141             DO  i = nxlg, nxrg
1142                DO  j = nysg, nyng
[1551]1143                   DO  k = nzb_do, nzt_do
[1]1144                      local_pf(i,j,k) = to_be_resorted(k,j,i)
1145                   ENDDO
1146                ENDDO
1147             ENDDO
1148          ENDIF
1149
1150!
1151!--       Output of the individual cross-sections, depending on the cross-
1152!--       section mode chosen.
1153          is = 1
[1551]1154   loop1: DO WHILE ( section(is,s) /= -9999  .OR.  two_d )
[1]1155
1156             SELECT CASE ( mode )
1157
1158                CASE ( 'xy' )
1159!
1160!--                Determine the cross section index
1161                   IF ( two_d )  THEN
1162                      layer_xy = nzb+1
1163                   ELSE
1164                      layer_xy = section(is,s)
1165                   ENDIF
1166
1167!
[1551]1168!--                Exit the loop for layers beyond the data output domain
1169!--                (used for soil model)
1170                   IF ( layer_xy .GT. nzt_do )  THEN
1171                      EXIT loop1
1172                   ENDIF
1173
1174!
[1308]1175!--                Update the netCDF xy cross section time axis.
1176!--                In case of parallel output, this is only done by PE0
1177!--                to increase the performance.
1178                   IF ( simulated_time /= do2d_xy_last_time(av) )  THEN
1179                      do2d_xy_time_count(av) = do2d_xy_time_count(av) + 1
1180                      do2d_xy_last_time(av)  = simulated_time
1181                      IF ( myid == 0 )  THEN
[1327]1182                         IF ( .NOT. data_output_2d_on_each_pe  &
1183                              .OR.  netcdf_data_format > 4 )   &
[493]1184                         THEN
[1]1185#if defined( __netcdf )
1186                            nc_stat = NF90_PUT_VAR( id_set_xy(av),             &
1187                                                    id_var_time_xy(av),        &
[291]1188                                             (/ time_since_reference_point /), &
[1]1189                                         start = (/ do2d_xy_time_count(av) /), &
1190                                                    count = (/ 1 /) )
[493]1191                            CALL handle_netcdf_error( 'data_output_2d', 53 )
[1]1192#endif
1193                         ENDIF
1194                      ENDIF
1195                   ENDIF
1196!
1197!--                If required, carry out averaging along z
[336]1198                   IF ( section(is,s) == -1  .AND.  .NOT. two_d )  THEN
[1]1199
[1353]1200                      local_2d = 0.0_wp
[1]1201!
1202!--                   Carry out the averaging (all data are on the PE)
[1551]1203                      DO  k = nzb_do, nzt_do
[667]1204                         DO  j = nysg, nyng
1205                            DO  i = nxlg, nxrg
[1]1206                               local_2d(i,j) = local_2d(i,j) + local_pf(i,j,k)
1207                            ENDDO
1208                         ENDDO
1209                      ENDDO
1210
[1551]1211                      local_2d = local_2d / ( nzt_do - nzb_do + 1.0_wp)
[1]1212
1213                   ELSE
1214!
1215!--                   Just store the respective section on the local array
1216                      local_2d = local_pf(:,:,layer_xy)
1217
1218                   ENDIF
1219
1220#if defined( __parallel )
[1327]1221                   IF ( netcdf_data_format > 4 )  THEN
[1]1222!
[1031]1223!--                   Parallel output in netCDF4/HDF5 format.
[493]1224                      IF ( two_d ) THEN
1225                         iis = 1
1226                      ELSE
1227                         iis = is
1228                      ENDIF
1229
[1]1230#if defined( __netcdf )
[1308]1231!
1232!--                   For parallel output, all cross sections are first stored
1233!--                   here on a local array and will be written to the output
1234!--                   file afterwards to increase the performance.
1235                      DO  i = nxlg, nxrg
1236                         DO  j = nysg, nyng
1237                            local_2d_sections(i,j,iis) = local_2d(i,j)
1238                         ENDDO
1239                      ENDDO
[1]1240#endif
[493]1241                   ELSE
[1]1242
[493]1243                      IF ( data_output_2d_on_each_pe )  THEN
[1]1244!
[493]1245!--                      Output of partial arrays on each PE
1246#if defined( __netcdf )
[1327]1247                         IF ( myid == 0 )  THEN
[1320]1248                            WRITE ( 21 )  time_since_reference_point,          &
[493]1249                                          do2d_xy_time_count(av), av
1250                         ENDIF
1251#endif
[759]1252                         DO  i = 0, io_blocks-1
1253                            IF ( i == io_group )  THEN
[1551]1254                               WRITE ( 21 )  nxlg, nxrg, nysg, nyng, nysg, nyng
[759]1255                               WRITE ( 21 )  local_2d
1256                            ENDIF
1257#if defined( __parallel )
1258                            CALL MPI_BARRIER( comm2d, ierr )
1259#endif
1260                         ENDDO
[559]1261
[493]1262                      ELSE
[1]1263!
[493]1264!--                      PE0 receives partial arrays from all processors and
1265!--                      then outputs them. Here a barrier has to be set,
1266!--                      because otherwise "-MPI- FATAL: Remote protocol queue
1267!--                      full" may occur.
1268                         CALL MPI_BARRIER( comm2d, ierr )
1269
[667]1270                         ngp = ( nxrg-nxlg+1 ) * ( nyng-nysg+1 )
[493]1271                         IF ( myid == 0 )  THEN
[1]1272!
[493]1273!--                         Local array can be relocated directly.
[667]1274                            total_2d(nxlg:nxrg,nysg:nyng) = local_2d
[1]1275!
[493]1276!--                         Receive data from all other PEs.
1277                            DO  n = 1, numprocs-1
[1]1278!
[493]1279!--                            Receive index limits first, then array.
1280!--                            Index limits are received in arbitrary order from
1281!--                            the PEs.
[1320]1282                               CALL MPI_RECV( ind(1), 4, MPI_INTEGER,          &
1283                                              MPI_ANY_SOURCE, 0, comm2d,       &
[493]1284                                              status, ierr )
1285                               sender = status(MPI_SOURCE)
1286                               DEALLOCATE( local_2d )
1287                               ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
[1320]1288                               CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp,    &
1289                                              MPI_REAL, sender, 1, comm2d,     &
[493]1290                                              status, ierr )
1291                               total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
1292                            ENDDO
[1]1293!
[493]1294!--                         Relocate the local array for the next loop increment
1295                            DEALLOCATE( local_2d )
[667]1296                            ALLOCATE( local_2d(nxlg:nxrg,nysg:nyng) )
[1]1297
1298#if defined( __netcdf )
[1327]1299                            IF ( two_d ) THEN
1300                               nc_stat = NF90_PUT_VAR( id_set_xy(av),       &
1301                                                       id_var_do2d(av,if),  &
1302                                                   total_2d(0:nx+1,0:ny+1), &
1303                             start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
1304                                             count = (/ nx+2, ny+2, 1, 1 /) )
1305                            ELSE
1306                               nc_stat = NF90_PUT_VAR( id_set_xy(av),       &
1307                                                       id_var_do2d(av,if),  &
1308                                                   total_2d(0:nx+1,0:ny+1), &
1309                            start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
1310                                             count = (/ nx+2, ny+2, 1, 1 /) )
[1]1311                            ENDIF
[1327]1312                            CALL handle_netcdf_error( 'data_output_2d', 54 )
[1]1313#endif
1314
[493]1315                         ELSE
[1]1316!
[493]1317!--                         First send the local index limits to PE0
[667]1318                            ind(1) = nxlg; ind(2) = nxrg
1319                            ind(3) = nysg; ind(4) = nyng
[1320]1320                            CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0,       &
[493]1321                                           comm2d, ierr )
[1]1322!
[493]1323!--                         Send data to PE0
[1320]1324                            CALL MPI_SEND( local_2d(nxlg,nysg), ngp,           &
[493]1325                                           MPI_REAL, 0, 1, comm2d, ierr )
1326                         ENDIF
1327!
1328!--                      A barrier has to be set, because otherwise some PEs may
1329!--                      proceed too fast so that PE0 may receive wrong data on
1330!--                      tag 0
1331                         CALL MPI_BARRIER( comm2d, ierr )
[1]1332                      ENDIF
[493]1333
[1]1334                   ENDIF
1335#else
1336#if defined( __netcdf )
[1327]1337                   IF ( two_d ) THEN
1338                      nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
1339                                              id_var_do2d(av,if),           &
1340                                             local_2d(nxl:nxr+1,nys:nyn+1), &
1341                             start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
1342                                           count = (/ nx+2, ny+2, 1, 1 /) )
1343                   ELSE
1344                      nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
1345                                              id_var_do2d(av,if),           &
1346                                             local_2d(nxl:nxr+1,nys:nyn+1), &
1347                            start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
1348                                           count = (/ nx+2, ny+2, 1, 1 /) )
[1]1349                   ENDIF
[1327]1350                   CALL handle_netcdf_error( 'data_output_2d', 447 )
[1]1351#endif
1352#endif
1353                   do2d_xy_n = do2d_xy_n + 1
1354!
1355!--                For 2D-arrays (e.g. u*) only one cross-section is available.
1356!--                Hence exit loop of output levels.
1357                   IF ( two_d )  THEN
1358                      two_d = .FALSE.
1359                      EXIT loop1
1360                   ENDIF
1361
1362                CASE ( 'xz' )
1363!
[1308]1364!--                Update the netCDF xz cross section time axis.
1365!--                In case of parallel output, this is only done by PE0
1366!--                to increase the performance.
1367                   IF ( simulated_time /= do2d_xz_last_time(av) )  THEN
1368                      do2d_xz_time_count(av) = do2d_xz_time_count(av) + 1
1369                      do2d_xz_last_time(av)  = simulated_time
1370                      IF ( myid == 0 )  THEN
[1327]1371                         IF ( .NOT. data_output_2d_on_each_pe  &
1372                              .OR.  netcdf_data_format > 4 )   &
[493]1373                         THEN
[1]1374#if defined( __netcdf )
1375                            nc_stat = NF90_PUT_VAR( id_set_xz(av),             &
1376                                                    id_var_time_xz(av),        &
[291]1377                                             (/ time_since_reference_point /), &
[1]1378                                         start = (/ do2d_xz_time_count(av) /), &
1379                                                    count = (/ 1 /) )
[493]1380                            CALL handle_netcdf_error( 'data_output_2d', 56 )
[1]1381#endif
1382                         ENDIF
1383                      ENDIF
1384                   ENDIF
[667]1385
[1]1386!
1387!--                If required, carry out averaging along y
1388                   IF ( section(is,s) == -1 )  THEN
1389
[1551]1390                      ALLOCATE( local_2d_l(nxlg:nxrg,nzb_do:nzt_do) )
[1353]1391                      local_2d_l = 0.0_wp
[1551]1392                      ngp = ( nxrg-nxlg + 1 ) * ( nzt_do-nzb_do + 1 )
[1]1393!
1394!--                   First local averaging on the PE
[1551]1395                      DO  k = nzb_do, nzt_do
[1]1396                         DO  j = nys, nyn
[667]1397                            DO  i = nxlg, nxrg
[1320]1398                               local_2d_l(i,k) = local_2d_l(i,k) +             &
[1]1399                                                 local_pf(i,j,k)
1400                            ENDDO
1401                         ENDDO
1402                      ENDDO
1403#if defined( __parallel )
1404!
1405!--                   Now do the averaging over all PEs along y
[622]1406                      IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1551]1407                      CALL MPI_ALLREDUCE( local_2d_l(nxlg,nzb_do),                &
1408                                          local_2d(nxlg,nzb_do), ngp, MPI_REAL,   &
[1]1409                                          MPI_SUM, comm1dy, ierr )
1410#else
1411                      local_2d = local_2d_l
1412#endif
[1353]1413                      local_2d = local_2d / ( ny + 1.0_wp )
[1]1414
1415                      DEALLOCATE( local_2d_l )
1416
1417                   ELSE
1418!
1419!--                   Just store the respective section on the local array
1420!--                   (but only if it is available on this PE!)
1421                      IF ( section(is,s) >= nys  .AND.  section(is,s) <= nyn ) &
1422                      THEN
[1551]1423                         local_2d = local_pf(:,section(is,s),nzb_do:nzt_do)
[1]1424                      ENDIF
1425
1426                   ENDIF
1427
1428#if defined( __parallel )
[1327]1429                   IF ( netcdf_data_format > 4 )  THEN
[1]1430!
[1031]1431!--                   Output in netCDF4/HDF5 format.
[493]1432!--                   Output only on those PEs where the respective cross
1433!--                   sections reside. Cross sections averaged along y are
1434!--                   output on the respective first PE along y (myidy=0).
[1320]1435                      IF ( ( section(is,s) >= nys  .AND.                       &
1436                             section(is,s) <= nyn )  .OR.                      &
[493]1437                           ( section(is,s) == -1  .AND.  myidy == 0 ) )  THEN
[1]1438#if defined( __netcdf )
[493]1439!
[1308]1440!--                      For parallel output, all cross sections are first
1441!--                      stored here on a local array and will be written to the
1442!--                      output file afterwards to increase the performance.
1443                         DO  i = nxlg, nxrg
[1551]1444                            DO  k = nzb_do, nzt_do
[1308]1445                               local_2d_sections_l(i,is,k) = local_2d(i,k)
1446                            ENDDO
1447                         ENDDO
[1]1448#endif
1449                      ENDIF
1450
1451                   ELSE
1452
[493]1453                      IF ( data_output_2d_on_each_pe )  THEN
[1]1454!
[493]1455!--                      Output of partial arrays on each PE. If the cross
1456!--                      section does not reside on the PE, output special
1457!--                      index values.
1458#if defined( __netcdf )
[1327]1459                         IF ( myid == 0 )  THEN
[1320]1460                            WRITE ( 22 )  time_since_reference_point,          &
[493]1461                                          do2d_xz_time_count(av), av
1462                         ENDIF
1463#endif
[759]1464                         DO  i = 0, io_blocks-1
1465                            IF ( i == io_group )  THEN
[1320]1466                               IF ( ( section(is,s) >= nys  .AND.              &
1467                                      section(is,s) <= nyn )  .OR.             &
1468                                    ( section(is,s) == -1  .AND.               &
1469                                      nys-1 == -1 ) )                          &
[759]1470                               THEN
[1551]1471                                  WRITE (22)  nxlg, nxrg, nzb_do, nzt_do, nzb, nzt+1
[759]1472                                  WRITE (22)  local_2d
1473                               ELSE
[1551]1474                                  WRITE (22)  -1, -1, -1, -1, -1, -1
[759]1475                               ENDIF
1476                            ENDIF
1477#if defined( __parallel )
1478                            CALL MPI_BARRIER( comm2d, ierr )
1479#endif
1480                         ENDDO
[493]1481
1482                      ELSE
[1]1483!
[493]1484!--                      PE0 receives partial arrays from all processors of the
1485!--                      respective cross section and outputs them. Here a
1486!--                      barrier has to be set, because otherwise
1487!--                      "-MPI- FATAL: Remote protocol queue full" may occur.
1488                         CALL MPI_BARRIER( comm2d, ierr )
1489
[1551]1490                         ngp = ( nxrg-nxlg + 1 ) * ( nzt_do-nzb_do + 1 )
[493]1491                         IF ( myid == 0 )  THEN
[1]1492!
[493]1493!--                         Local array can be relocated directly.
1494                            IF ( ( section(is,s) >= nys  .AND.                 &
1495                                   section(is,s) <= nyn )  .OR.                &
1496                                 ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
1497                            THEN
[1551]1498                               total_2d(nxlg:nxrg,nzb_do:nzt_do) = local_2d
[493]1499                            ENDIF
[1]1500!
[493]1501!--                         Receive data from all other PEs.
1502                            DO  n = 1, numprocs-1
1503!
1504!--                            Receive index limits first, then array.
1505!--                            Index limits are received in arbitrary order from
1506!--                            the PEs.
[1320]1507                               CALL MPI_RECV( ind(1), 4, MPI_INTEGER,          &
1508                                              MPI_ANY_SOURCE, 0, comm2d,       &
[1]1509                                              status, ierr )
[493]1510!
1511!--                            Not all PEs have data for XZ-cross-section.
1512                               IF ( ind(1) /= -9999 )  THEN
1513                                  sender = status(MPI_SOURCE)
1514                                  DEALLOCATE( local_2d )
[1320]1515                                  ALLOCATE( local_2d(ind(1):ind(2),            &
[493]1516                                                     ind(3):ind(4)) )
1517                                  CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1518                                                 MPI_REAL, sender, 1, comm2d,  &
1519                                                 status, ierr )
[1320]1520                                  total_2d(ind(1):ind(2),ind(3):ind(4)) =      &
[493]1521                                                                        local_2d
1522                               ENDIF
1523                            ENDDO
1524!
1525!--                         Relocate the local array for the next loop increment
1526                            DEALLOCATE( local_2d )
[1551]1527                            ALLOCATE( local_2d(nxlg:nxrg,nzb_do:nzt_do) )
[1]1528
1529#if defined( __netcdf )
[1327]1530                            nc_stat = NF90_PUT_VAR( id_set_xz(av),          &
1531                                                 id_var_do2d(av,if),        &
[1551]1532                                                 total_2d(0:nx+1,nzb_do:nzt_do),&
[1327]1533                            start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
[1551]1534                                             count = (/ nx+2, 1, nzt_do-nzb_do+1, 1 /) )
[1327]1535                            CALL handle_netcdf_error( 'data_output_2d', 58 )
[1]1536#endif
1537
[493]1538                         ELSE
[1]1539!
[493]1540!--                         If the cross section resides on the PE, send the
1541!--                         local index limits, otherwise send -9999 to PE0.
1542                            IF ( ( section(is,s) >= nys  .AND.                 &
1543                                   section(is,s) <= nyn )  .OR.                &
1544                                 ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
1545                            THEN
[667]1546                               ind(1) = nxlg; ind(2) = nxrg
[1551]1547                               ind(3) = nzb_do;   ind(4) = nzt_do
[493]1548                            ELSE
1549                               ind(1) = -9999; ind(2) = -9999
1550                               ind(3) = -9999; ind(4) = -9999
1551                            ENDIF
[1320]1552                            CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0,       &
[493]1553                                           comm2d, ierr )
1554!
1555!--                         If applicable, send data to PE0.
1556                            IF ( ind(1) /= -9999 )  THEN
[1551]1557                               CALL MPI_SEND( local_2d(nxlg,nzb_do), ngp,         &
[493]1558                                              MPI_REAL, 0, 1, comm2d, ierr )
1559                            ENDIF
[1]1560                         ENDIF
1561!
[493]1562!--                      A barrier has to be set, because otherwise some PEs may
1563!--                      proceed too fast so that PE0 may receive wrong data on
1564!--                      tag 0
1565                         CALL MPI_BARRIER( comm2d, ierr )
[1]1566                      ENDIF
[493]1567
[1]1568                   ENDIF
1569#else
1570#if defined( __netcdf )
[1327]1571                   nc_stat = NF90_PUT_VAR( id_set_xz(av),                   &
1572                                           id_var_do2d(av,if),              &
[1551]1573                                           local_2d(nxl:nxr+1,nzb_do:nzt_do),   &
[1327]1574                            start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
[1551]1575                                           count = (/ nx+2, 1, nzt_do-nzb_do+1, 1 /) )
[1327]1576                   CALL handle_netcdf_error( 'data_output_2d', 451 )
[1]1577#endif
1578#endif
1579                   do2d_xz_n = do2d_xz_n + 1
1580
1581                CASE ( 'yz' )
1582!
[1308]1583!--                Update the netCDF yz cross section time axis.
1584!--                In case of parallel output, this is only done by PE0
1585!--                to increase the performance.
1586                   IF ( simulated_time /= do2d_yz_last_time(av) )  THEN
1587                      do2d_yz_time_count(av) = do2d_yz_time_count(av) + 1
1588                      do2d_yz_last_time(av)  = simulated_time
1589                      IF ( myid == 0 )  THEN
[1327]1590                         IF ( .NOT. data_output_2d_on_each_pe  &
1591                              .OR.  netcdf_data_format > 4 )   &
[493]1592                         THEN
[1]1593#if defined( __netcdf )
1594                            nc_stat = NF90_PUT_VAR( id_set_yz(av),             &
1595                                                    id_var_time_yz(av),        &
[291]1596                                             (/ time_since_reference_point /), &
[1]1597                                         start = (/ do2d_yz_time_count(av) /), &
1598                                                    count = (/ 1 /) )
[263]1599                            CALL handle_netcdf_error( 'data_output_2d', 59 )
[1]1600#endif
1601                         ENDIF
1602                      ENDIF
[1308]1603                   ENDIF
[493]1604
[1]1605!
1606!--                If required, carry out averaging along x
1607                   IF ( section(is,s) == -1 )  THEN
1608
[1551]1609                      ALLOCATE( local_2d_l(nysg:nyng,nzb_do:nzt_do) )
[1353]1610                      local_2d_l = 0.0_wp
[1551]1611                      ngp = ( nyng-nysg+1 ) * ( nzt_do-nzb_do+1 )
[1]1612!
1613!--                   First local averaging on the PE
[1551]1614                      DO  k = nzb_do, nzt_do
[667]1615                         DO  j = nysg, nyng
[1]1616                            DO  i = nxl, nxr
[1320]1617                               local_2d_l(j,k) = local_2d_l(j,k) +             &
[1]1618                                                 local_pf(i,j,k)
1619                            ENDDO
1620                         ENDDO
1621                      ENDDO
1622#if defined( __parallel )
1623!
1624!--                   Now do the averaging over all PEs along x
[622]1625                      IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1551]1626                      CALL MPI_ALLREDUCE( local_2d_l(nysg,nzb_do),                &
1627                                          local_2d(nysg,nzb_do), ngp, MPI_REAL,   &
[1]1628                                          MPI_SUM, comm1dx, ierr )
1629#else
1630                      local_2d = local_2d_l
1631#endif
[1353]1632                      local_2d = local_2d / ( nx + 1.0_wp )
[1]1633
1634                      DEALLOCATE( local_2d_l )
1635
1636                   ELSE
1637!
1638!--                   Just store the respective section on the local array
1639!--                   (but only if it is available on this PE!)
1640                      IF ( section(is,s) >= nxl  .AND.  section(is,s) <= nxr ) &
1641                      THEN
[1551]1642                         local_2d = local_pf(section(is,s),:,nzb_do:nzt_do)
[1]1643                      ENDIF
1644
1645                   ENDIF
1646
1647#if defined( __parallel )
[1327]1648                   IF ( netcdf_data_format > 4 )  THEN
[1]1649!
[1031]1650!--                   Output in netCDF4/HDF5 format.
[493]1651!--                   Output only on those PEs where the respective cross
1652!--                   sections reside. Cross sections averaged along x are
1653!--                   output on the respective first PE along x (myidx=0).
[1320]1654                      IF ( ( section(is,s) >= nxl  .AND.                       &
1655                             section(is,s) <= nxr )  .OR.                      &
[493]1656                           ( section(is,s) == -1  .AND.  myidx == 0 ) )  THEN
[1]1657#if defined( __netcdf )
[493]1658!
[1308]1659!--                      For parallel output, all cross sections are first
1660!--                      stored here on a local array and will be written to the
1661!--                      output file afterwards to increase the performance.
1662                         DO  j = nysg, nyng
[1551]1663                            DO  k = nzb_do, nzt_do
[1308]1664                               local_2d_sections_l(is,j,k) = local_2d(j,k)
1665                            ENDDO
1666                         ENDDO
[1]1667#endif
1668                      ENDIF
1669
1670                   ELSE
1671
[493]1672                      IF ( data_output_2d_on_each_pe )  THEN
[1]1673!
[493]1674!--                      Output of partial arrays on each PE. If the cross
1675!--                      section does not reside on the PE, output special
1676!--                      index values.
1677#if defined( __netcdf )
[1327]1678                         IF ( myid == 0 )  THEN
[1320]1679                            WRITE ( 23 )  time_since_reference_point,          &
[493]1680                                          do2d_yz_time_count(av), av
1681                         ENDIF
1682#endif
[759]1683                         DO  i = 0, io_blocks-1
1684                            IF ( i == io_group )  THEN
[1320]1685                               IF ( ( section(is,s) >= nxl  .AND.              &
1686                                      section(is,s) <= nxr )  .OR.             &
1687                                    ( section(is,s) == -1  .AND.               &
1688                                      nxl-1 == -1 ) )                          &
[759]1689                               THEN
[1551]1690                                  WRITE (23)  nysg, nyng, nzb_do, nzt_do, nzb, nzt+1
[759]1691                                  WRITE (23)  local_2d
1692                               ELSE
[1551]1693                                  WRITE (23)  -1, -1, -1, -1, -1, -1
[759]1694                               ENDIF
1695                            ENDIF
1696#if defined( __parallel )
1697                            CALL MPI_BARRIER( comm2d, ierr )
1698#endif
1699                         ENDDO
[493]1700
1701                      ELSE
[1]1702!
[493]1703!--                      PE0 receives partial arrays from all processors of the
1704!--                      respective cross section and outputs them. Here a
1705!--                      barrier has to be set, because otherwise
1706!--                      "-MPI- FATAL: Remote protocol queue full" may occur.
1707                         CALL MPI_BARRIER( comm2d, ierr )
1708
[1551]1709                         ngp = ( nyng-nysg+1 ) * ( nzt_do-nzb_do+1 )
[493]1710                         IF ( myid == 0 )  THEN
[1]1711!
[493]1712!--                         Local array can be relocated directly.
1713                            IF ( ( section(is,s) >= nxl  .AND.                 &
1714                                   section(is,s) <= nxr )   .OR.               &
1715                                 ( section(is,s) == -1  .AND.  nxl-1 == -1 ) ) &
1716                            THEN
[1551]1717                               total_2d(nysg:nyng,nzb_do:nzt_do) = local_2d
[493]1718                            ENDIF
[1]1719!
[493]1720!--                         Receive data from all other PEs.
1721                            DO  n = 1, numprocs-1
1722!
1723!--                            Receive index limits first, then array.
1724!--                            Index limits are received in arbitrary order from
1725!--                            the PEs.
[1320]1726                               CALL MPI_RECV( ind(1), 4, MPI_INTEGER,          &
1727                                              MPI_ANY_SOURCE, 0, comm2d,       &
[1]1728                                              status, ierr )
[493]1729!
1730!--                            Not all PEs have data for YZ-cross-section.
1731                               IF ( ind(1) /= -9999 )  THEN
1732                                  sender = status(MPI_SOURCE)
1733                                  DEALLOCATE( local_2d )
[1320]1734                                  ALLOCATE( local_2d(ind(1):ind(2),            &
[493]1735                                                     ind(3):ind(4)) )
1736                                  CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
1737                                                 MPI_REAL, sender, 1, comm2d,  &
1738                                                 status, ierr )
[1320]1739                                  total_2d(ind(1):ind(2),ind(3):ind(4)) =      &
[493]1740                                                                        local_2d
1741                               ENDIF
1742                            ENDDO
1743!
1744!--                         Relocate the local array for the next loop increment
1745                            DEALLOCATE( local_2d )
[1551]1746                            ALLOCATE( local_2d(nysg:nyng,nzb_do:nzt_do) )
[1]1747
1748#if defined( __netcdf )
[1327]1749                            nc_stat = NF90_PUT_VAR( id_set_yz(av),          &
1750                                                 id_var_do2d(av,if),        &
[1551]1751                                                 total_2d(0:ny+1,nzb_do:nzt_do),&
[1327]1752                            start = (/ is, 1, 1, do2d_yz_time_count(av) /), &
[1551]1753                                             count = (/ 1, ny+2, nzt_do-nzb_do+1, 1 /) )
[1327]1754                            CALL handle_netcdf_error( 'data_output_2d', 61 )
[1]1755#endif
1756
[493]1757                         ELSE
[1]1758!
[493]1759!--                         If the cross section resides on the PE, send the
1760!--                         local index limits, otherwise send -9999 to PE0.
1761                            IF ( ( section(is,s) >= nxl  .AND.                 &
1762                                   section(is,s) <= nxr )  .OR.                &
1763                                 ( section(is,s) == -1  .AND.  nxl-1 == -1 ) ) &
1764                            THEN
[667]1765                               ind(1) = nysg; ind(2) = nyng
[1551]1766                               ind(3) = nzb_do;   ind(4) = nzt_do
[493]1767                            ELSE
1768                               ind(1) = -9999; ind(2) = -9999
1769                               ind(3) = -9999; ind(4) = -9999
1770                            ENDIF
[1320]1771                            CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0,       &
[493]1772                                           comm2d, ierr )
1773!
1774!--                         If applicable, send data to PE0.
1775                            IF ( ind(1) /= -9999 )  THEN
[1551]1776                               CALL MPI_SEND( local_2d(nysg,nzb_do), ngp,         &
[493]1777                                              MPI_REAL, 0, 1, comm2d, ierr )
1778                            ENDIF
[1]1779                         ENDIF
1780!
[493]1781!--                      A barrier has to be set, because otherwise some PEs may
1782!--                      proceed too fast so that PE0 may receive wrong data on
1783!--                      tag 0
1784                         CALL MPI_BARRIER( comm2d, ierr )
[1]1785                      ENDIF
[493]1786
[1]1787                   ENDIF
1788#else
1789#if defined( __netcdf )
[1327]1790                   nc_stat = NF90_PUT_VAR( id_set_yz(av),                   &
1791                                           id_var_do2d(av,if),              &
[1551]1792                                           local_2d(nys:nyn+1,nzb_do:nzt_do),   &
[1327]1793                            start = (/ is, 1, 1, do2d_xz_time_count(av) /), &
[1551]1794                                           count = (/ 1, ny+2, nzt_do-nzb_do+1, 1 /) )
[1327]1795                   CALL handle_netcdf_error( 'data_output_2d', 452 )
[1]1796#endif
1797#endif
1798                   do2d_yz_n = do2d_yz_n + 1
1799
1800             END SELECT
1801
1802             is = is + 1
1803          ENDDO loop1
1804
[1308]1805!
1806!--       For parallel output, all data were collected before on a local array
1807!--       and are written now to the netcdf file. This must be done to increase
1808!--       the performance of the parallel output.
1809#if defined( __netcdf )
[1327]1810          IF ( netcdf_data_format > 4 )  THEN
[1308]1811
1812                SELECT CASE ( mode )
1813
1814                   CASE ( 'xy' )
1815                      IF ( two_d ) THEN
1816                         iis = 1
1817                      ELSE
1818                         iis = is-1
1819                      ENDIF
1820!
1821!--                   Do not output redundant ghost point data except for the
1822!--                   boundaries of the total domain.
1823                      IF ( nxr == nx  .AND.  nyn /= ny )  THEN
1824                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
1825                                                 id_var_do2d(av,if),           &
1826                                                 local_2d_sections(nxl:nxr+1,  &
1827                                                    nys:nyn,1:ns),             &
1828                                                 start = (/ nxl+1, nys+1, 1,   &
1829                                                    do2d_xy_time_count(av) /), &
1830                                                 count = (/ nxr-nxl+2,         &
1831                                                            nyn-nys+1, ns, 1   &
1832                                                          /) )
1833                      ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
1834                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
1835                                                 id_var_do2d(av,if),           &
1836                                                 local_2d_sections(nxl:nxr,    &
1837                                                    nys:nyn+1,1:ns),           &
1838                                                 start = (/ nxl+1, nys+1, 1,   &
1839                                                    do2d_xy_time_count(av) /), &
1840                                                 count = (/ nxr-nxl+1,         &
1841                                                            nyn-nys+2, ns, 1   &
1842                                                          /) )
1843                      ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
1844                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
1845                                                 id_var_do2d(av,if),           &
1846                                                 local_2d_sections(nxl:nxr+1,  &
1847                                                    nys:nyn+1,1:ns),           &
1848                                                 start = (/ nxl+1, nys+1, 1,   &
1849                                                    do2d_xy_time_count(av) /), &
1850                                                 count = (/ nxr-nxl+2,         &
1851                                                            nyn-nys+2, ns, 1   &
1852                                                          /) )
1853                      ELSE
1854                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
1855                                                 id_var_do2d(av,if),           &
1856                                                 local_2d_sections(nxl:nxr,    &
1857                                                    nys:nyn,1:ns),             &
1858                                                 start = (/ nxl+1, nys+1, 1,   &
1859                                                    do2d_xy_time_count(av) /), &
1860                                                 count = (/ nxr-nxl+1,         &
1861                                                            nyn-nys+1, ns, 1   &
1862                                                          /) )
1863                      ENDIF   
1864
1865                      CALL handle_netcdf_error( 'data_output_2d', 55 ) 
1866
1867                   CASE ( 'xz' )
1868!
1869!--                   First, all PEs get the information of all cross-sections.
1870!--                   Then the data are written to the output file by all PEs
1871!--                   while NF90_COLLECTIVE is set in subroutine
1872!--                   define_netcdf_header. Although redundant information are
1873!--                   written to the output file in that case, the performance
1874!--                   is significantly better compared to the case where only
1875!--                   the first row of PEs in x-direction (myidx = 0) is given
1876!--                   the output while NF90_INDEPENDENT is set.
1877                      IF ( npey /= 1 )  THEN
1878                         
1879#if defined( __parallel )
1880!
1881!--                      Distribute data over all PEs along y
[1551]1882                         ngp = ( nxrg-nxlg+1 ) * ( nzt_do-nzb_do+1 ) * ns
[1308]1883                         IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1551]1884                         CALL MPI_ALLREDUCE( local_2d_sections_l(nxlg,1,nzb_do),  &
1885                                             local_2d_sections(nxlg,1,nzb_do),    &
[1308]1886                                             ngp, MPI_REAL, MPI_SUM, comm1dy,  &
1887                                             ierr )
1888#else
1889                         local_2d_sections = local_2d_sections_l
1890#endif
1891                      ENDIF
1892!
1893!--                   Do not output redundant ghost point data except for the
1894!--                   boundaries of the total domain.
1895                      IF ( nxr == nx )  THEN
1896                         nc_stat = NF90_PUT_VAR( id_set_xz(av),                &
1897                                             id_var_do2d(av,if),               & 
1898                                             local_2d_sections(nxl:nxr+1,1:ns, &
[1551]1899                                                nzb_do:nzt_do),                &
[1308]1900                                             start = (/ nxl+1, 1, 1,           &
1901                                                do2d_xz_time_count(av) /),     &
[1551]1902                                             count = (/ nxr-nxl+2, ns, nzt_do-nzb_do+1,  &
[1308]1903                                                        1 /) )
1904                      ELSE
1905                         nc_stat = NF90_PUT_VAR( id_set_xz(av),                &
1906                                             id_var_do2d(av,if),               &
1907                                             local_2d_sections(nxl:nxr,1:ns,   &
[1551]1908                                                nzb_do:nzt_do),                &
[1308]1909                                             start = (/ nxl+1, 1, 1,           &
1910                                                do2d_xz_time_count(av) /),     &
[1551]1911                                             count = (/ nxr-nxl+1, ns, nzt_do-nzb_do+1,  &
[1308]1912                                                1 /) )
1913                      ENDIF
1914
1915                      CALL handle_netcdf_error( 'data_output_2d', 57 )
1916
1917                   CASE ( 'yz' )
1918!
1919!--                   First, all PEs get the information of all cross-sections.
1920!--                   Then the data are written to the output file by all PEs
1921!--                   while NF90_COLLECTIVE is set in subroutine
1922!--                   define_netcdf_header. Although redundant information are
1923!--                   written to the output file in that case, the performance
1924!--                   is significantly better compared to the case where only
1925!--                   the first row of PEs in y-direction (myidy = 0) is given
1926!--                   the output while NF90_INDEPENDENT is set.
1927                      IF ( npex /= 1 )  THEN
1928
1929#if defined( __parallel )
1930!
1931!--                      Distribute data over all PEs along x
1932                         ngp = ( nyng-nysg+1 ) * ( nzt-nzb + 2 ) * ns
1933                         IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr )
[1551]1934                         CALL MPI_ALLREDUCE( local_2d_sections_l(1,nysg,nzb_do),  &
1935                                             local_2d_sections(1,nysg,nzb_do),    &
[1308]1936                                             ngp, MPI_REAL, MPI_SUM, comm1dx,  &
1937                                             ierr )
1938#else
1939                         local_2d_sections = local_2d_sections_l
1940#endif
1941                      ENDIF
1942!
1943!--                   Do not output redundant ghost point data except for the
1944!--                   boundaries of the total domain.
1945                      IF ( nyn == ny )  THEN
1946                         nc_stat = NF90_PUT_VAR( id_set_yz(av),                &
1947                                             id_var_do2d(av,if),               &
1948                                             local_2d_sections(1:ns,           &
[1551]1949                                                nys:nyn+1,nzb_do:nzt_do),      &
[1308]1950                                             start = (/ 1, nys+1, 1,           &
1951                                                do2d_yz_time_count(av) /),     &
1952                                             count = (/ ns, nyn-nys+2,         &
[1551]1953                                                        nzt_do-nzb_do+1, 1 /) )
[1308]1954                      ELSE
1955                         nc_stat = NF90_PUT_VAR( id_set_yz(av),                &
1956                                             id_var_do2d(av,if),               &
1957                                             local_2d_sections(1:ns,nys:nyn,   &
[1551]1958                                                nzb_do:nzt_do),                &
[1308]1959                                             start = (/ 1, nys+1, 1,           &
1960                                                do2d_yz_time_count(av) /),     &
1961                                             count = (/ ns, nyn-nys+1,         &
[1551]1962                                                        nzt_do-nzb_do+1, 1 /) )
[1308]1963                      ENDIF
1964
1965                      CALL handle_netcdf_error( 'data_output_2d', 60 )
1966
1967                   CASE DEFAULT
1968                      message_string = 'unknown cross-section: ' // TRIM( mode )
1969                      CALL message( 'data_output_2d', 'PA0180', 1, 2, 0, 6, 0 )
1970
1971                END SELECT                     
1972
1973          ENDIF
[1311]1974#endif
[1]1975       ENDIF
1976
1977       if = if + 1
1978       l = MAX( 2, LEN_TRIM( do2d(av,if) ) )
1979       do2d_mode = do2d(av,if)(l-1:l)
1980
1981    ENDDO
1982
1983!
1984!-- Deallocate temporary arrays.
1985    IF ( ALLOCATED( level_z ) )  DEALLOCATE( level_z )
[1308]1986    IF ( netcdf_data_format > 4 )  THEN
1987       DEALLOCATE( local_pf, local_2d, local_2d_sections )
1988       IF( mode == 'xz' .OR. mode == 'yz' ) DEALLOCATE( local_2d_sections_l )
1989    ENDIF
[1]1990#if defined( __parallel )
1991    IF ( .NOT.  data_output_2d_on_each_pe  .AND.  myid == 0 )  THEN
1992       DEALLOCATE( total_2d )
1993    ENDIF
1994#endif
1995
1996!
1997!-- Close plot output file.
1998    file_id = 20 + s
1999
2000    IF ( data_output_2d_on_each_pe )  THEN
[759]2001       DO  i = 0, io_blocks-1
2002          IF ( i == io_group )  THEN
2003             CALL close_file( file_id )
2004          ENDIF
2005#if defined( __parallel )
2006          CALL MPI_BARRIER( comm2d, ierr )
2007#endif
2008       ENDDO
[1]2009    ELSE
2010       IF ( myid == 0 )  CALL close_file( file_id )
2011    ENDIF
2012
[1318]2013    CALL cpu_log( log_point(3), 'data_output_2d', 'stop' )
[1]2014
2015 END SUBROUTINE data_output_2d
Note: See TracBrowser for help on using the repository browser.