source: palm/tags/release-3.3/SOURCE/data_output_2d.f90 @ 3890

Last change on this file since 3890 was 98, checked in by raasch, 17 years ago

updating comments and rc-file

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