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

Last change on this file since 775 was 772, checked in by heinze, 13 years ago

last commit documented

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