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

Last change on this file since 1065 was 1065, checked in by hoffmann, 9 years ago

cloud physics: rain sedimentation and turbulence effects

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