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

Last change on this file since 328 was 291, checked in by raasch, 16 years ago

changes for coupling with independent precursor runs; z_i calculation with Sullivan criterion

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