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

Last change on this file since 96 was 96, checked in by raasch, 14 years ago

more preliminary uncomplete changes for ocean version

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