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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

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