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

Last change on this file since 587 was 559, checked in by weinreis, 14 years ago


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