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

Last change on this file since 16 was 4, checked in by raasch, 18 years ago

Id keyword set as property for all *.f90 files

  • Property svn:keywords set to Id
File size: 44.4 KB
Line 
1 SUBROUTINE data_output_2d( mode, av )
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: data_output_2d.f90 4 2007-02-13 11:33:16Z 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 ( 'pt_xy', 'pt_xz', 'pt_yz' )
260                IF ( av == 0 )  THEN
261                   IF ( .NOT. cloud_physics ) THEN
262                      to_be_resorted => pt
263                   ELSE
264                      DO  i = nxl-1, nxr+1
265                         DO  j = nys-1, nyn+1
266                            DO  k = nzb, nzt+1
267                               local_pf(i,j,k) = pt(k,j,i) + l_d_cp *    &
268                                                             pt_d_t(k) * &
269                                                             ql(k,j,i)
270                            ENDDO
271                         ENDDO
272                      ENDDO
273                      resorted = .TRUE.
274                   ENDIF
275                ELSE
276                   to_be_resorted => pt_av
277                ENDIF
278                IF ( mode == 'xy' )  level_z = zu
279
280             CASE ( 'q_xy', 'q_xz', 'q_yz' )
281                IF ( av == 0 )  THEN
282                   to_be_resorted => q
283                ELSE
284                   to_be_resorted => q_av
285                ENDIF
286                IF ( mode == 'xy' )  level_z = zu
287
288             CASE ( 'ql_xy', 'ql_xz', 'ql_yz' )
289                IF ( av == 0 )  THEN
290                   to_be_resorted => ql
291                ELSE
292                   to_be_resorted => ql_av
293                ENDIF
294                IF ( mode == 'xy' )  level_z = zu
295
296             CASE ( 'ql_c_xy', 'ql_c_xz', 'ql_c_yz' )
297                IF ( av == 0 )  THEN
298                   to_be_resorted => ql_c
299                ELSE
300                   to_be_resorted => ql_c_av
301                ENDIF
302                IF ( mode == 'xy' )  level_z = zu
303
304             CASE ( 'ql_v_xy', 'ql_v_xz', 'ql_v_yz' )
305                IF ( av == 0 )  THEN
306                   to_be_resorted => ql_v
307                ELSE
308                   to_be_resorted => ql_v_av
309                ENDIF
310                IF ( mode == 'xy' )  level_z = zu
311
312             CASE ( 'ql_vp_xy', 'ql_vp_xz', 'ql_vp_yz' )
313                IF ( av == 0 )  THEN
314                   to_be_resorted => ql_vp
315                ELSE
316                   to_be_resorted => ql_vp_av
317                ENDIF
318                IF ( mode == 'xy' )  level_z = zu
319
320             CASE ( 'qv_xy', 'qv_xz', 'qv_yz' )
321                IF ( av == 0 )  THEN
322                   DO  i = nxl-1, nxr+1
323                      DO  j = nys-1, nyn+1
324                         DO  k = nzb, nzt+1
325                            local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
326                         ENDDO
327                      ENDDO
328                   ENDDO
329                   resorted = .TRUE.
330                ELSE
331                   to_be_resorted => qv_av
332                ENDIF
333                IF ( mode == 'xy' )  level_z = zu
334
335             CASE ( 's_xy', 's_xz', 's_yz' )
336                IF ( av == 0 )  THEN
337                   to_be_resorted => q
338                ELSE
339                   to_be_resorted => q_av
340                ENDIF
341
342             CASE ( 't*_xy' )        ! 2d-array
343                IF ( av == 0 )  THEN
344                   DO  i = nxl-1, nxr+1
345                      DO  j = nys-1, nyn+1
346                         local_pf(i,j,nzb+1) = ts(j,i)
347                      ENDDO
348                   ENDDO
349                ELSE
350                   DO  i = nxl-1, nxr+1
351                      DO  j = nys-1, nyn+1
352                         local_pf(i,j,nzb+1) = ts_av(j,i)
353                      ENDDO
354                   ENDDO
355                ENDIF
356                resorted = .TRUE.
357                two_d = .TRUE.
358                level_z(nzb+1) = zu(nzb+1)
359
360             CASE ( 'u_xy', 'u_xz', 'u_yz' )
361                IF ( av == 0 )  THEN
362                   to_be_resorted => u
363                ELSE
364                   to_be_resorted => u_av
365                ENDIF
366                IF ( mode == 'xy' )  level_z = zu
367!
368!--             Substitute the values generated by "mirror" boundary condition
369!--             at the bottom boundary by the real surface values.
370                IF ( do2d(av,if) == 'u_xz'  .OR.  do2d(av,if) == 'u_yz' )  THEN
371                   IF ( ibc_uv_b == 0 )  local_pf(:,:,nzb) = 0.0
372                ENDIF
373
374             CASE ( 'u*_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) = us(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) = us_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 ( 'v_xy', 'v_xz', 'v_yz' )
393                IF ( av == 0 )  THEN
394                   to_be_resorted => v
395                ELSE
396                   to_be_resorted => v_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) == 'v_xz'  .OR.  do2d(av,if) == 'v_yz' )  THEN
403                   IF ( ibc_uv_b == 0 )  local_pf(:,:,nzb) = 0.0
404                ENDIF
405
406             CASE ( 'vpt_xy', 'vpt_xz', 'vpt_yz' )
407                IF ( av == 0 )  THEN
408                   to_be_resorted => vpt
409                ELSE
410                   to_be_resorted => vpt_av
411                ENDIF
412                IF ( mode == 'xy' )  level_z = zu
413
414             CASE ( 'w_xy', 'w_xz', 'w_yz' )
415                IF ( av == 0 )  THEN
416                   to_be_resorted => w
417                ELSE
418                   to_be_resorted => w_av
419                ENDIF
420                IF ( mode == 'xy' )  level_z = zw
421
422             CASE DEFAULT
423!
424!--             User defined quantity
425                CALL user_data_output_2d( av, do2d(av,if), found, grid, &
426                                          local_pf )
427                resorted = .TRUE.
428
429                IF ( grid == 'zu' )  THEN
430                   IF ( mode == 'xy' )  level_z = zu
431                ELSEIF ( grid == 'zw' )  THEN
432                   IF ( mode == 'xy' )  level_z = zw
433                ENDIF
434
435                IF ( .NOT. found )  THEN
436                   PRINT*, '+++ data_output_2d: no output provided for: ', &
437                                do2d(av,if)
438                ENDIF
439
440          END SELECT
441
442!
443!--       Resort the array to be output, if not done above
444          IF ( .NOT. resorted )  THEN
445             DO  i = nxl-1, nxr+1
446                DO  j = nys-1, nyn+1
447                   DO  k = nzb, nzt+1
448                      local_pf(i,j,k) = to_be_resorted(k,j,i)
449                   ENDDO
450                ENDDO
451             ENDDO
452          ENDIF
453
454!
455!--       Output of the individual cross-sections, depending on the cross-
456!--       section mode chosen.
457          is = 1
458   loop1: DO  WHILE ( section(is,s) /= -9999  .OR.  two_d )
459
460             SELECT CASE ( mode )
461
462                CASE ( 'xy' )
463!
464!--                Determine the cross section index
465                   IF ( two_d )  THEN
466                      layer_xy = nzb+1
467                   ELSE
468                      layer_xy = section(is,s)
469                   ENDIF
470
471!
472!--                Update the NetCDF xy cross section time axis
473                   IF ( myid == 0 )  THEN
474                      IF ( simulated_time /= do2d_xy_last_time(av) )  THEN
475                         do2d_xy_time_count(av) = do2d_xy_time_count(av) + 1
476                         do2d_xy_last_time(av)  = simulated_time
477                         IF ( .NOT. data_output_2d_on_each_pe  .AND. &
478                              netcdf_output )  THEN
479#if defined( __netcdf )
480                            nc_stat = NF90_PUT_VAR( id_set_xy(av),             &
481                                                    id_var_time_xy(av),        &
482                                                    (/ simulated_time /),      &
483                                         start = (/ do2d_xy_time_count(av) /), &
484                                                    count = (/ 1 /) )
485                            IF ( nc_stat /= NF90_NOERR )  THEN
486                               CALL handle_netcdf_error( 53 )
487                            ENDIF
488#endif
489                         ENDIF
490                      ENDIF
491                   ENDIF
492!
493!--                If required, carry out averaging along z
494                   IF ( section(is,s) == -1 )  THEN
495
496                      local_2d = 0.0
497!
498!--                   Carry out the averaging (all data are on the PE)
499                      DO  k = nzb, nzt+1
500                         DO  j = nys-1, nyn+1
501                            DO  i = nxl-1, nxr+1
502                               local_2d(i,j) = local_2d(i,j) + local_pf(i,j,k)
503                            ENDDO
504                         ENDDO
505                      ENDDO
506
507                      local_2d = local_2d / ( nzt -nzb + 2.0 )
508
509                   ELSE
510!
511!--                   Just store the respective section on the local array
512                      local_2d = local_pf(:,:,layer_xy)
513
514                   ENDIF
515
516#if defined( __parallel )
517                   IF ( data_output_2d_on_each_pe )  THEN
518!
519!--                   Output of partial arrays on each PE
520#if defined( __netcdf )
521                      IF ( netcdf_output  .AND.  myid == 0 )  THEN
522                         WRITE ( 21 )  simulated_time, do2d_xy_time_count(av), &
523                                       av
524                      ENDIF
525#endif
526                      WRITE ( 21 )  nxl-1, nxr+1, nys-1, nyn+1
527                      WRITE ( 21 )  local_2d
528
529                   ELSE
530!
531!--                   PE0 receives partial arrays from all processors and then
532!--                   outputs them. Here a barrier has to be set, because
533!--                   otherwise "-MPI- FATAL: Remote protocol queue full" may
534!--                   occur.
535                      CALL MPI_BARRIER( comm2d, ierr )
536
537                      ngp = ( nxr-nxl+3 ) * ( nyn-nys+3 )
538                      IF ( myid == 0 )  THEN
539!
540!--                      Local array can be relocated directly.
541                         total_2d(nxl-1:nxr+1,nys-1:nyn+1) = local_2d
542!
543!--                      Receive data from all other PEs.
544                         DO  n = 1, numprocs-1
545!
546!--                         Receive index limits first, then array.
547!--                         Index limits are received in arbitrary order from
548!--                         the PEs.
549                            CALL MPI_RECV( ind(1), 4, MPI_INTEGER,            &
550                                           MPI_ANY_SOURCE, 0, comm2d, status, &
551                                           ierr )
552                            sender = status(MPI_SOURCE)
553                            DEALLOCATE( local_2d )
554                            ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
555                            CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp,      &
556                                           MPI_REAL, sender, 1, comm2d,       &
557                                           status, ierr )
558                            total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
559                         ENDDO
560!
561!--                      Output of the total cross-section.
562                         IF ( iso2d_output ) WRITE (21)  total_2d(0:nx+1,0:ny+1)
563!
564!--                      Relocate the local array for the next loop increment
565                         DEALLOCATE( local_2d )
566                         ALLOCATE( local_2d(nxl-1:nxr+1,nys-1:nyn+1) )
567
568#if defined( __netcdf )
569                         IF ( netcdf_output )  THEN
570                            IF ( two_d ) THEN
571                               nc_stat = NF90_PUT_VAR( id_set_xy(av),          &
572                                                       id_var_do2d(av,if),     &
573                                                      total_2d(0:nx+1,0:ny+1), &
574                                start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
575                                                count = (/ nx+2, ny+2, 1, 1 /) )
576                            ELSE
577                               nc_stat = NF90_PUT_VAR( id_set_xy(av),          &
578                                                       id_var_do2d(av,if),     &
579                                                      total_2d(0:nx+1,0:ny+1), &
580                               start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
581                                                count = (/ nx+2, ny+2, 1, 1 /) )
582                            ENDIF
583                            IF ( nc_stat /= NF90_NOERR )  &
584                                                  CALL handle_netcdf_error( 54 )
585                         ENDIF
586#endif
587
588                      ELSE
589!
590!--                      First send the local index limits to PE0
591                         ind(1) = nxl-1; ind(2) = nxr+1
592                         ind(3) = nys-1; ind(4) = nyn+1
593                         CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, comm2d, &
594                                        ierr )
595!
596!--                      Send data to PE0
597                         CALL MPI_SEND( local_2d(nxl-1,nys-1), ngp, MPI_REAL, &
598                                        0, 1, comm2d, ierr )
599                      ENDIF
600!
601!--                   A barrier has to be set, because otherwise some PEs may
602!--                   proceed too fast so that PE0 may receive wrong data on
603!--                   tag 0
604                      CALL MPI_BARRIER( comm2d, ierr )
605                   ENDIF
606#else
607                   IF ( iso2d_output )  THEN
608                      WRITE (21)  local_2d(nxl:nxr+1,nys:nyn+1)
609                   ENDIF
610#if defined( __netcdf )
611                   IF ( netcdf_output )  THEN
612                      IF ( two_d ) THEN
613                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
614                                                 id_var_do2d(av,if),           &
615                                                local_2d(nxl:nxr+1,nys:nyn+1), &
616                                start = (/ 1, 1, 1, do2d_xy_time_count(av) /), &
617                                              count = (/ nx+2, ny+2, 1, 1 /) )
618                      ELSE
619                         nc_stat = NF90_PUT_VAR( id_set_xy(av),                &
620                                                 id_var_do2d(av,if),           &
621                                                local_2d(nxl:nxr+1,nys:nyn+1), &
622                               start = (/ 1, 1, is, do2d_xy_time_count(av) /), &
623                                              count = (/ nx+2, ny+2, 1, 1 /) )
624                      ENDIF
625                      IF ( nc_stat /= NF90_NOERR )  &
626                                                  CALL handle_netcdf_error( 55 )
627                   ENDIF
628#endif
629#endif
630                   do2d_xy_n = do2d_xy_n + 1
631!
632!--                Write LOCAL parameter set for ISO2D.
633                   IF ( myid == 0  .AND.  iso2d_output )  THEN
634                      IF ( section(is,s) /= -1 )  THEN
635                         WRITE ( section_chr, '(''z = '',F7.2,'' m  (GP '',I3, &
636                                               &'')'')'                        &
637                               )  level_z(layer_xy), layer_xy
638                      ELSE
639                         section_chr = 'averaged along z'
640                      ENDIF
641                      IF ( av == 0 )  THEN
642                         rtext = TRIM( do2d(av,if) ) // '  t = ' //    &
643                                 TRIM( simulated_time_chr ) // '  ' // &
644                                 TRIM( section_chr )
645                      ELSE
646                         rtext = TRIM( do2d(av,if) ) // '  averaged t = ' // &
647                                 TRIM( simulated_time_chr ) // '  ' //       &
648                                 TRIM( section_chr )
649                      ENDIF
650                      WRITE (27,LOCAL)
651                   ENDIF
652!
653!--                For 2D-arrays (e.g. u*) only one cross-section is available.
654!--                Hence exit loop of output levels.
655                   IF ( two_d )  THEN
656                      two_d = .FALSE.
657                      EXIT loop1
658                   ENDIF
659
660                CASE ( 'xz' )
661!
662!--                Update the NetCDF xy cross section time axis
663                   IF ( myid == 0 )  THEN
664                      IF ( simulated_time /= do2d_xz_last_time(av) )  THEN
665                         do2d_xz_time_count(av) = do2d_xz_time_count(av) + 1
666                         do2d_xz_last_time(av)  = simulated_time
667                         IF ( .NOT. data_output_2d_on_each_pe  .AND. &
668                              netcdf_output )  THEN
669#if defined( __netcdf )
670                            nc_stat = NF90_PUT_VAR( id_set_xz(av),             &
671                                                    id_var_time_xz(av),        &
672                                                    (/ simulated_time /),      &
673                                         start = (/ do2d_xz_time_count(av) /), &
674                                                    count = (/ 1 /) )
675                            IF ( nc_stat /= NF90_NOERR )  THEN
676                               CALL handle_netcdf_error( 56 )
677                            ENDIF
678#endif
679                         ENDIF
680                      ENDIF
681                   ENDIF
682!
683!--                If required, carry out averaging along y
684                   IF ( section(is,s) == -1 )  THEN
685
686                      ALLOCATE( local_2d_l(nxl-1:nxr+1,nzb:nzt+1) )
687                      local_2d_l = 0.0
688                      ngp = ( nxr-nxl+3 ) * ( nzt-nzb+2 )
689!
690!--                   First local averaging on the PE
691                      DO  k = nzb, nzt+1
692                         DO  j = nys, nyn
693                            DO  i = nxl-1, nxr+1
694                               local_2d_l(i,k) = local_2d_l(i,k) + &
695                                                 local_pf(i,j,k)
696                            ENDDO
697                         ENDDO
698                      ENDDO
699#if defined( __parallel )
700!
701!--                   Now do the averaging over all PEs along y
702                      CALL MPI_ALLREDUCE( local_2d_l(nxl-1,nzb),              &
703                                          local_2d(nxl-1,nzb), ngp, MPI_REAL, &
704                                          MPI_SUM, comm1dy, ierr )
705#else
706                      local_2d = local_2d_l
707#endif
708                      local_2d = local_2d / ( ny + 1.0 )
709
710                      DEALLOCATE( local_2d_l )
711
712                   ELSE
713!
714!--                   Just store the respective section on the local array
715!--                   (but only if it is available on this PE!)
716                      IF ( section(is,s) >= nys  .AND.  section(is,s) <= nyn ) &
717                      THEN
718                         local_2d = local_pf(:,section(is,s),nzb:nzt+1)
719                      ENDIF
720
721                   ENDIF
722
723#if defined( __parallel )
724                   IF ( data_output_2d_on_each_pe )  THEN
725!
726!--                   Output of partial arrays on each PE. If the cross section
727!--                   does not reside on the PE, output special index values.
728#if defined( __netcdf )
729                      IF ( netcdf_output  .AND.  myid == 0 )  THEN
730                         WRITE ( 22 )  simulated_time, do2d_xz_time_count(av), &
731                                       av
732                      ENDIF
733#endif
734                      IF ( ( section(is,s)>=nys .AND. section(is,s)<=nyn ) .OR.&
735                           ( section(is,s) == -1  .AND.  nys-1 == -1 ) )       &
736                      THEN
737                         WRITE (22)  nxl-1, nxr+1, nzb, nzt+1
738                         WRITE (22)  local_2d
739                      ELSE
740                         WRITE (22)  -1, -1, -1, -1
741                      ENDIF
742
743                   ELSE
744!
745!--                   PE0 receives partial arrays from all processors of the
746!--                   respective cross section and outputs them. Here a
747!--                   barrier has to be set, because otherwise
748!--                   "-MPI- FATAL: Remote protocol queue full" may occur.
749                      CALL MPI_BARRIER( comm2d, ierr )
750
751                      ngp = ( nxr-nxl+3 ) * ( nzt-nzb+2 )
752                      IF ( myid == 0 )  THEN
753!
754!--                      Local array can be relocated directly.
755                         IF ( ( section(is,s)>=nys .AND. section(is,s)<=nyn )  &
756                            .OR. ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
757                         THEN
758                            total_2d(nxl-1:nxr+1,nzb:nzt+1) = local_2d
759                         ENDIF
760!
761!--                      Receive data from all other PEs.
762                         DO  n = 1, numprocs-1
763!
764!--                         Receive index limits first, then array.
765!--                         Index limits are received in arbitrary order from
766!--                         the PEs.
767                            CALL MPI_RECV( ind(1), 4, MPI_INTEGER,            &
768                                           MPI_ANY_SOURCE, 0, comm2d, status, &
769                                           ierr )
770!
771!--                         Not all PEs have data for XZ-cross-section.
772                            IF ( ind(1) /= -9999 )  THEN
773                               sender = status(MPI_SOURCE)
774                               DEALLOCATE( local_2d )
775                               ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
776                               CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
777                                              MPI_REAL, sender, 1, comm2d,  &
778                                              status, ierr )
779                               total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
780                            ENDIF
781                         ENDDO
782!
783!--                      Output of the total cross-section.
784                         IF ( iso2d_output )  THEN
785                            WRITE (22)  total_2d(0:nx+1,nzb:nzt+1)
786                         ENDIF
787!
788!--                      Relocate the local array for the next loop increment
789                         DEALLOCATE( local_2d )
790                         ALLOCATE( local_2d(nxl-1:nxr+1,nzb:nzt+1) )
791
792#if defined( __netcdf )
793                         IF ( netcdf_output )  THEN
794                            nc_stat = NF90_PUT_VAR( id_set_xz(av),             &
795                                                    id_var_do2d(av,if),        &
796                                                    total_2d(0:nx+1,nzb:nzt+1),&
797                               start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
798                                                count = (/ nx+2, 1, nz+2, 1 /) )
799                            IF ( nc_stat /= NF90_NOERR ) &
800                                                  CALL handle_netcdf_error( 57 )
801                         ENDIF
802#endif
803
804                      ELSE
805!
806!--                      If the cross section resides on the PE, send the
807!--                      local index limits, otherwise send -9999 to PE0.
808                         IF ( ( section(is,s)>=nys .AND. section(is,s)<=nyn )  &
809                            .OR. ( section(is,s) == -1  .AND.  nys-1 == -1 ) ) &
810                         THEN
811                            ind(1) = nxl-1; ind(2) = nxr+1
812                            ind(3) = nzb;   ind(4) = nzt+1
813                         ELSE
814                            ind(1) = -9999; ind(2) = -9999
815                            ind(3) = -9999; ind(4) = -9999
816                         ENDIF
817                         CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, comm2d, &
818                                        ierr )
819!
820!--                      If applicable, send data to PE0.
821                         IF ( ind(1) /= -9999 )  THEN
822                            CALL MPI_SEND( local_2d(nxl-1,nzb), ngp, MPI_REAL, &
823                                           0, 1, comm2d, ierr )
824                         ENDIF
825                      ENDIF
826!
827!--                   A barrier has to be set, because otherwise some PEs may
828!--                   proceed too fast so that PE0 may receive wrong data on
829!--                   tag 0
830                      CALL MPI_BARRIER( comm2d, ierr )
831                   ENDIF
832#else
833                   IF ( iso2d_output )  THEN
834                      WRITE (22)  local_2d(nxl:nxr+1,nzb:nzt+1)
835                   ENDIF
836#if defined( __netcdf )
837                   IF ( netcdf_output )  THEN
838                      nc_stat = NF90_PUT_VAR( id_set_xz(av),                   &
839                                              id_var_do2d(av,if),              &
840                                              local_2d(nxl:nxr+1,nzb:nzt+1),   &
841                               start = (/ 1, is, 1, do2d_xz_time_count(av) /), &
842                                              count = (/ nx+2, 1, nz+2, 1 /) )
843                      IF ( nc_stat /= NF90_NOERR )  &
844                                                  CALL handle_netcdf_error( 58 )
845                   ENDIF
846#endif
847#endif
848                   do2d_xz_n = do2d_xz_n + 1
849!
850!--                Write LOCAL-parameter set for ISO2D.
851                   IF ( myid == 0  .AND.  iso2d_output )  THEN
852                      IF ( section(is,s) /= -1 )  THEN
853                         WRITE ( section_chr, '(''y = '',F8.2,'' m  (GP '',I3, &
854                                               &'')'')'                        &
855                               )  section(is,s)*dy, section(is,s)
856                      ELSE
857                         section_chr = 'averaged along y'
858                      ENDIF
859                      IF ( av == 0 )  THEN
860                         rtext = TRIM( do2d(av,if) ) // '  t = ' //    &
861                                 TRIM( simulated_time_chr ) // '  ' // &
862                                 TRIM( section_chr )
863                      ELSE
864                         rtext = TRIM( do2d(av,if) ) // '  averaged t = ' // &
865                                 TRIM( simulated_time_chr ) // '  ' //       &
866                                 TRIM( section_chr )
867                      ENDIF
868                      WRITE (28,LOCAL)
869                   ENDIF
870
871                CASE ( 'yz' )
872!
873!--                Update the NetCDF xy cross section time axis
874                   IF ( myid == 0 )  THEN
875                      IF ( simulated_time /= do2d_yz_last_time(av) )  THEN
876                         do2d_yz_time_count(av) = do2d_yz_time_count(av) + 1
877                         do2d_yz_last_time(av)  = simulated_time
878                         IF ( .NOT. data_output_2d_on_each_pe  .AND. &
879                              netcdf_output )  THEN
880#if defined( __netcdf )
881                            nc_stat = NF90_PUT_VAR( id_set_yz(av),             &
882                                                    id_var_time_yz(av),        &
883                                                    (/ simulated_time /),      &
884                                         start = (/ do2d_yz_time_count(av) /), &
885                                                    count = (/ 1 /) )
886                            IF ( nc_stat /= NF90_NOERR )  THEN
887                               CALL handle_netcdf_error( 59 )
888                            ENDIF
889#endif
890                         ENDIF
891                      ENDIF
892                   ENDIF
893!
894!--                If required, carry out averaging along x
895                   IF ( section(is,s) == -1 )  THEN
896
897                      ALLOCATE( local_2d_l(nys-1:nyn+1,nzb:nzt+1) )
898                      local_2d_l = 0.0
899                      ngp = ( nyn-nys+3 ) * ( nzt-nzb+2 )
900!
901!--                   First local averaging on the PE
902                      DO  k = nzb, nzt+1
903                         DO  j = nys-1, nyn+1
904                            DO  i = nxl, nxr
905                               local_2d_l(j,k) = local_2d_l(j,k) + &
906                                                 local_pf(i,j,k)
907                            ENDDO
908                         ENDDO
909                      ENDDO
910#if defined( __parallel )
911!
912!--                   Now do the averaging over all PEs along x
913                      CALL MPI_ALLREDUCE( local_2d_l(nys-1,nzb),              &
914                                          local_2d(nys-1,nzb), ngp, MPI_REAL, &
915                                          MPI_SUM, comm1dx, ierr )
916#else
917                      local_2d = local_2d_l
918#endif
919                      local_2d = local_2d / ( nx + 1.0 )
920
921                      DEALLOCATE( local_2d_l )
922
923                   ELSE
924!
925!--                   Just store the respective section on the local array
926!--                   (but only if it is available on this PE!)
927                      IF ( section(is,s) >= nxl  .AND.  section(is,s) <= nxr ) &
928                      THEN
929                         local_2d = local_pf(section(is,s),:,nzb:nzt+1)
930                      ENDIF
931
932                   ENDIF
933
934#if defined( __parallel )
935                   IF ( data_output_2d_on_each_pe )  THEN
936!
937!--                   Output of partial arrays on each PE. If the cross section
938!--                   does not reside on the PE, output special index values.
939#if defined( __netcdf )
940                      IF ( netcdf_output  .AND.  myid == 0 )  THEN
941                         WRITE ( 23 )  simulated_time, do2d_yz_time_count(av), &
942                                       av
943                      ENDIF
944#endif
945                      IF ( ( section(is,s)>=nxl .AND. section(is,s)<=nxr ) .OR.&
946                           ( section(is,s) ==  -1  .AND.  nxl-1 == -1 ) )      &
947                      THEN
948                         WRITE (23)  nys-1, nyn+1, nzb, nzt+1
949                         WRITE (23)  local_2d
950                      ELSE
951                         WRITE (23)  -1, -1, -1, -1
952                      ENDIF
953
954                   ELSE
955!
956!--                   PE0 receives partial arrays from all processors of the
957!--                   respective cross section and outputs them. Here a
958!--                   barrier has to be set, because otherwise
959!--                   "-MPI- FATAL: Remote protocol queue full" may occur.
960                      CALL MPI_BARRIER( comm2d, ierr )
961
962                      ngp = ( nyn-nys+3 ) * ( nzt-nzb+2 )
963                      IF ( myid == 0 )  THEN
964!
965!--                      Local array can be relocated directly.
966                         IF ( ( section(is,s)>=nxl .AND. section(is,s)<=nxr )  &
967                           .OR. ( section(is,s) ==  -1  .AND.  nxl-1 == -1 ) ) &
968                         THEN
969                            total_2d(nys-1:nyn+1,nzb:nzt+1) = local_2d
970                         ENDIF
971!
972!--                      Receive data from all other PEs.
973                         DO  n = 1, numprocs-1
974!
975!--                         Receive index limits first, then array.
976!--                         Index limits are received in arbitrary order from
977!--                         the PEs.
978                            CALL MPI_RECV( ind(1), 4, MPI_INTEGER,            &
979                                           MPI_ANY_SOURCE, 0, comm2d, status, &
980                                           ierr )
981!
982!--                         Not all PEs have data for YZ-cross-section.
983                            IF ( ind(1) /= -9999 )  THEN
984                               sender = status(MPI_SOURCE)
985                               DEALLOCATE( local_2d )
986                               ALLOCATE( local_2d(ind(1):ind(2),ind(3):ind(4)) )
987                               CALL MPI_RECV( local_2d(ind(1),ind(3)), ngp, &
988                                              MPI_REAL, sender, 1, comm2d,  &
989                                              status, ierr )
990                               total_2d(ind(1):ind(2),ind(3):ind(4)) = local_2d
991                            ENDIF
992                         ENDDO
993!
994!--                      Output of the total cross-section.
995                         IF ( iso2d_output )  THEN
996                            WRITE (23)  total_2d(0:ny+1,nzb:nzt+1)
997                         ENDIF
998!
999!--                      Relocate the local array for the next loop increment
1000                         DEALLOCATE( local_2d )
1001                         ALLOCATE( local_2d(nys-1:nyn+1,nzb:nzt+1) )
1002
1003#if defined( __netcdf )
1004                         IF ( netcdf_output )  THEN
1005                            nc_stat = NF90_PUT_VAR( id_set_yz(av),             &
1006                                                    id_var_do2d(av,if),        &
1007                                                    total_2d(0:ny+1,nzb:nzt+1),&
1008                               start = (/ is, 1, 1, do2d_yz_time_count(av) /), &
1009                                                count = (/ 1, ny+2, nz+2, 1 /) )
1010                            IF ( nc_stat /= NF90_NOERR ) &
1011                                                  CALL handle_netcdf_error( 60 )
1012                         ENDIF
1013#endif
1014
1015                      ELSE
1016!
1017!--                      If the cross section resides on the PE, send the
1018!--                      local index limits, otherwise send -9999 to PE0.
1019                         IF ( ( section(is,s)>=nxl .AND. section(is,s)<=nxr )  &
1020                           .OR. ( section(is,s) ==  -1  .AND.  nxl-1 == -1 ) ) &
1021                         THEN
1022                            ind(1) = nys-1; ind(2) = nyn+1
1023                            ind(3) = nzb;   ind(4) = nzt+1
1024                         ELSE
1025                            ind(1) = -9999; ind(2) = -9999
1026                            ind(3) = -9999; ind(4) = -9999
1027                         ENDIF
1028                         CALL MPI_SEND( ind(1), 4, MPI_INTEGER, 0, 0, comm2d, &
1029                                        ierr )
1030!
1031!--                      If applicable, send data to PE0.
1032                         IF ( ind(1) /= -9999 )  THEN
1033                            CALL MPI_SEND( local_2d(nys-1,nzb), ngp, MPI_REAL, &
1034                                           0, 1, comm2d, ierr )
1035                         ENDIF
1036                      ENDIF
1037!
1038!--                   A barrier has to be set, because otherwise some PEs may
1039!--                   proceed too fast so that PE0 may receive wrong data on
1040!--                   tag 0
1041                      CALL MPI_BARRIER( comm2d, ierr )
1042                   ENDIF
1043#else
1044                   IF ( iso2d_output )  THEN
1045                      WRITE (23)  local_2d(nys:nyn+1,nzb:nzt+1)
1046                   ENDIF
1047#if defined( __netcdf )
1048                   IF ( netcdf_output )  THEN
1049                      nc_stat = NF90_PUT_VAR( id_set_yz(av),                   &
1050                                              id_var_do2d(av,if),              &
1051                                              local_2d(nys:nyn+1,nzb:nzt+1),   &
1052                               start = (/ is, 1, 1, do2d_xz_time_count(av) /), &
1053                                              count = (/ 1, ny+2, nz+2, 1 /) )
1054                      IF ( nc_stat /= NF90_NOERR )  &
1055                                                  CALL handle_netcdf_error( 61 )
1056                   ENDIF
1057#endif
1058#endif
1059                   do2d_yz_n = do2d_yz_n + 1
1060!
1061!--                Write LOCAL-parameter set for ISO2D.
1062                   IF ( myid == 0  .AND.  iso2d_output )  THEN
1063                      IF ( section(is,s) /= -1 )  THEN
1064                         WRITE ( section_chr, '(''x = '',F8.2,'' m  (GP '',I3, &
1065                                               &'')'')'                        &
1066                               )  section(is,s)*dx, section(is,s)
1067                      ELSE
1068                         section_chr = 'averaged along x'
1069                      ENDIF
1070                      IF ( av == 0 )  THEN
1071                         rtext = TRIM( do2d(av,if) ) // '  t = ' //    &
1072                                 TRIM( simulated_time_chr ) // '  ' // &
1073                                 TRIM( section_chr )
1074                      ELSE
1075                         rtext = TRIM( do2d(av,if) ) // '  averaged t = ' // &
1076                                 TRIM( simulated_time_chr ) // '  ' //       &
1077                                 TRIM( section_chr )
1078                      ENDIF
1079                      WRITE (29,LOCAL)
1080                   ENDIF
1081
1082             END SELECT
1083
1084             is = is + 1
1085          ENDDO loop1
1086
1087       ENDIF
1088
1089       if = if + 1
1090       l = MAX( 2, LEN_TRIM( do2d(av,if) ) )
1091       do2d_mode = do2d(av,if)(l-1:l)
1092
1093    ENDDO
1094
1095!
1096!-- Deallocate temporary arrays.
1097    IF ( ALLOCATED( level_z ) )  DEALLOCATE( level_z )
1098    DEALLOCATE( local_pf, local_2d )
1099#if defined( __parallel )
1100    IF ( .NOT.  data_output_2d_on_each_pe  .AND.  myid == 0 )  THEN
1101       DEALLOCATE( total_2d )
1102    ENDIF
1103#endif
1104
1105!
1106!-- Close plot output file.
1107    file_id = 20 + s
1108
1109    IF ( data_output_2d_on_each_pe )  THEN
1110       CALL close_file( file_id )
1111    ELSE
1112       IF ( myid == 0 )  CALL close_file( file_id )
1113    ENDIF
1114
1115
1116    CALL cpu_log (log_point(3),'data_output_2d','stop','nobarrier')
1117
1118 END SUBROUTINE data_output_2d
Note: See TracBrowser for help on using the repository browser.