source: palm/tags/release-3.9/SOURCE/data_output_2d.f90 @ 3968

Last change on this file since 3968 was 1037, checked in by raasch, 11 years ago

last commit documented

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