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

Last change on this file since 771 was 771, checked in by heinze, 12 years ago

Output of liquid water potential temperature in case of cloud_physics=.T. enabled

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