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

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

netCDF4 without parallel file support implemented
additional define string netcdf4_parallel is required to switch on parallel file support

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