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

Last change on this file since 631 was 623, checked in by raasch, 14 years ago

last commit documented

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