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

Last change on this file since 263 was 263, checked in by heinze, 15 years ago

Output of NetCDF messages with aid of message handling routine.

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