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

Last change on this file since 673 was 673, checked in by suehring, 13 years ago

Right computation of the pressure using Runge-Kutta weighting coefficients. Consideration of the pressure gradient during the time integration removed. Removed bugfix concerning velocity variances.

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