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

Last change on this file since 790 was 790, checked in by raasch, 12 years ago

Bugfix for output of mean particle radius + preliminary works for implementing the Wang collision kernel

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