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

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

last commit documented

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