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

Last change on this file since 75 was 75, checked in by raasch, 14 years ago

preliminary update for changes concerning non-cyclic boundary conditions

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