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

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

several small bugfixes; some more dvrp changes

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