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

Last change on this file since 1329 was 1329, checked in by raasch, 10 years ago

last commit documented

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