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

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

New:
---
Output in NetCDF4-format. New d3par-parameter netcdf_data_format.

(check_open, check_parameters, close_file, data_output_2d, data_output_3d, header, modules, netcdf, parin)

Modules to be loaded for compilation (mbuild) or job execution (mrun)
can be given in the configuration file using variable modules. Example:

%modules ifort/11.0.069:netcdf lcsgih parallel

This method replaces the (undocumented) mpilib-variable.

WARNING: All fixed settings of modules in the scripts mbuild, mrun, and subjob
have been removed! Please set the modules variable appropriately in your
configuration file. (mbuild, mrun, subjob)

Changed:


Parameters netcdf_64bit and netcdf_64bit_3d have been removed. Use
netcdf_data_format = 2 for choosing the classic 64bit-offset format (this is
the default). The offset-format can not be set independently for the
3d-output-data any more.

Parameters netcdf_format_mask, netcdf_format_mask_av, and variables
nc_format_mask, format_parallel_io removed. They are replaced by the new
parameter netcdf_data_format. (check_open, close_file,
data_output_mask, header, init_masks, modules, parin)

Errors:


bugfix in trunk/UTIL/Makefile: forgot to compile for interpret_config

Bugfix: timeseries data have to be collected by PE0 (user_statistics)

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