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

Last change on this file since 1115 was 1115, checked in by hoffmann, 8 years ago

optimization of two-moments cloud physics

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