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

Last change on this file since 771 was 771, checked in by heinze, 12 years ago

Output of liquid water potential temperature in case of cloud_physics=.T. enabled

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