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

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

preliminary changes for precipitation output

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