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

Last change on this file since 1054 was 1054, checked in by hoffmann, 11 years ago

last commit documented

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