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

Last change on this file since 285 was 274, checked in by heinze, 15 years ago

Indentation of the message calls corrected

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