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

Last change on this file since 255 was 254, checked in by heinze, 16 years ago

Output of messages replaced by message handling routine.

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