source: palm/tags/release-3.6/SOURCE/data_output_2d.f90 @ 3979

Last change on this file since 3979 was 226, checked in by raasch, 15 years ago

preparations for the next release

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