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

Last change on this file since 791 was 791, checked in by raasch, 13 years ago

last commit documented

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