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

Last change on this file since 758 was 730, checked in by heinze, 14 years ago

last commit documented

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