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

Last change on this file since 1291 was 1116, checked in by hoffmann, 12 years ago

last commit documented

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