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

Last change on this file since 219 was 215, checked in by raasch, 16 years ago

precompilation mechanism completely revised: now one depository per configuration block, further change of output messages

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