source: palm/trunk/SOURCE/data_output_mask.f90 @ 4453

Last change on this file since 4453 was 4444, checked in by raasch, 4 years ago

bugfix: cpp-directives for serial mode added

  • Property svn:keywords set to Id
File size: 32.5 KB
Line 
1!> @file data_output_mask.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2020 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: data_output_mask.f90 4444 2020-03-05 15:59:50Z raasch $
27! bugfix: cpp-directives for serial mode added
28!
29! 4377 2020-01-15 11:10:51Z gronemeier
30! bugfix: set fill value for output according to wall_flags_total_0 for
31!         non-terrain following output
32!
33! 4360 2020-01-07 11:25:50Z suehring
34! Introduction of wall_flags_total_0, which currently sets bits based on static
35! topography information used in wall_flags_static_0
36!
37! 4331 2019-12-10 18:25:02Z suehring
38! Formatting adjustment
39!
40! 4329 2019-12-10 15:46:36Z motisi
41! Renamed wall_flags_0 to wall_flags_static_0
42!
43! 4246 2019-09-30 09:27:52Z pavelkrc
44! Corrected "Former revisions" section
45!
46! 4168 2019-08-16 13:50:17Z suehring
47! Remove variable grid
48!
49! 4167 2019-08-16 11:01:48Z suehring
50! Changed behaviour of masked output over surface to follow terrain and ignore
51! buildings (J.Resler, T.Gronemeier)
52!
53! 4069 2019-07-01 14:05:51Z Giersch
54! Masked output running index mid has been introduced as a local variable to
55! avoid runtime error (Loop variable has been modified) in time_integration
56!
57! 4039 2019-06-18 10:32:41Z suehring
58! Modularize diagnostic output
59!
60! 3994 2019-05-22 18:08:09Z suehring
61! output of turbulence intensity added
62!
63! 3665 2019-01-10 08:28:24Z raasch
64! unused variables removed
65!
66! 3655 2019-01-07 16:51:22Z knoop
67! Fix output time levels (use time_since_reference_point)
68!
69! 410 2009-12-04 17:05:40Z letzel
70! Initial version
71!
72! Description:
73! ------------
74!> Masked data output in netCDF format for current mask (current value of mid).
75!------------------------------------------------------------------------------!
76 SUBROUTINE data_output_mask( av, mid )
77
78
79
80#if defined( __netcdf )
81    USE arrays_3d,                                                             &
82        ONLY:  e, nc, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, sa,  &
83               tend, u, v, vpt, w, d_exner
84
85    USE averaging,                                                             &
86        ONLY:  e_av, lpt_av, nc_av, nr_av, p_av, pc_av, pr_av, pt_av, q_av,    &
87               qc_av, ql_av, ql_c_av, ql_v_av, ql_vp_av, qv_av, qr_av,         &
88               rho_ocean_av, s_av, sa_av, u_av, v_av, vpt_av, w_av
89
90    USE basic_constants_and_equations_mod,                                     &
91        ONLY:  lv_d_cp
92
93    USE chemistry_model_mod,                                                   &
94        ONLY:  chem_data_output_mask
95
96    USE control_parameters,                                                    &
97        ONLY:  air_chemistry, domask, domask_no, domask_time_count, mask_i,    &
98               mask_j, mask_k, mask_size_l, mask_surface,                                                   &
99               max_masks, message_string, nz_do3d, salsa,                      &
100               time_since_reference_point
101
102#if defined( __parallel )
103    USE control_parameters,                                                    &
104        ONLY:  mask_size, mask_start_l
105#endif
106
107    USE diagnostic_output_quantities_mod,                                      &
108        ONLY:  doq_output_mask
109
110    USE cpulog,                                                                &
111        ONLY:  cpu_log, log_point
112
113    USE indices,                                                               &
114        ONLY:  nbgp, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_total_0
115
116    USE kinds
117
118    USE bulk_cloud_model_mod,                                                  &
119        ONLY:  bulk_cloud_model
120
121    USE NETCDF
122
123    USE netcdf_interface,                                                      &
124        ONLY:  fill_value, id_set_mask, id_var_domask, id_var_time_mask,       &
125               nc_stat, netcdf_data_format, netcdf_handle_error
126
127    USE particle_attributes,                                                   &
128        ONLY:  grid_particles, number_of_particles, particles,                 &
129               particle_advection_start, prt_count
130
131    USE pegrid
132
133    USE radiation_model_mod,                                                   &
134        ONLY:  radiation, radiation_data_output_mask
135
136    USE salsa_mod,                                                             &
137        ONLY:  salsa_data_output_mask
138
139
140    IMPLICIT NONE
141
142    INTEGER(iwp) ::  av                      !< flag for (non-)average output
143    INTEGER(iwp) ::  flag_nr                 !< number of masking flag
144    INTEGER(iwp) ::  i                       !< loop index
145    INTEGER(iwp) ::  ivar                    !< variable index
146    INTEGER(iwp) ::  j                       !< loop index
147    INTEGER(iwp) ::  k                       !< loop index
148    INTEGER(iwp) ::  im                      !< loop index for masked variables
149    INTEGER(iwp) ::  jm                      !< loop index for masked variables
150    INTEGER(iwp) ::  kk                      !< vertical index
151    INTEGER(iwp) ::  mid                     !< masked output running index
152    INTEGER(iwp) ::  n                       !< loop index
153    INTEGER(iwp) ::  netcdf_data_format_save !< value of netcdf_data_format
154    INTEGER(iwp) ::  ktt                     !< k index of highest terrain surface
155#if defined( __parallel )
156    INTEGER(iwp) ::  ngp                     !< number of grid points of an output slice
157    INTEGER(iwp) ::  sender                  !< PE id of sending PE
158    INTEGER(iwp) ::  ind(6)                  !< index limits (lower/upper bounds) of array 'local_2d'
159#endif
160
161    LOGICAL ::  found      !< true if output variable was found
162    LOGICAL ::  resorted   !< true if variable is resorted
163
164    REAL(wp) ::  mean_r    !< mean particle radius
165    REAL(wp) ::  s_r2      !< sum( particle-radius**2 )
166    REAL(wp) ::  s_r3      !< sum( particle-radius**3 )
167
168    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  local_pf    !< output array
169#if defined( __parallel )
170    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  total_pf    !< collected output array
171#endif
172    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to array which shall be output
173
174!
175!-- Return, if nothing to output
176    IF ( domask_no(mid,av) == 0 )  RETURN
177
178    CALL cpu_log (log_point(49),'data_output_mask','start')
179
180!
181!-- Parallel netcdf output is not tested so far for masked data, hence
182!-- netcdf_data_format is switched back to non-paralell output.
183    netcdf_data_format_save = netcdf_data_format
184    IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
185    IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
186
187!
188!-- Open output file.
189    IF ( myid == 0  .OR.  netcdf_data_format > 4 )  THEN
190       CALL check_open( 200+mid+av*max_masks )
191    ENDIF
192
193!
194!-- Allocate total and local output arrays.
195#if defined( __parallel )
196    IF ( myid == 0 )  THEN
197       ALLOCATE( total_pf(mask_size(mid,1),mask_size(mid,2),mask_size(mid,3)) )
198    ENDIF
199#endif
200    ALLOCATE( local_pf(mask_size_l(mid,1),mask_size_l(mid,2), &
201                       mask_size_l(mid,3)) )
202
203!
204!-- Update the netCDF time axis.
205    domask_time_count(mid,av) = domask_time_count(mid,av) + 1
206    IF ( myid == 0  .OR.  netcdf_data_format > 4 )  THEN
207       nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_time_mask(mid,av), &
208                               (/ time_since_reference_point /),              &
209                               start = (/ domask_time_count(mid,av) /),       &
210                               count = (/ 1 /) )
211       CALL netcdf_handle_error( 'data_output_mask', 460 )
212    ENDIF
213
214!
215!-- Loop over all variables to be written.
216    ivar = 1
217
218    DO  WHILE ( domask(mid,av,ivar)(1:1) /= ' ' )
219!
220!--    Reallocate local_pf on PE 0 since its shape changes during MPI exchange
221       IF ( netcdf_data_format < 5   .AND.  myid == 0  .AND.  ivar > 1 )  THEN
222          DEALLOCATE( local_pf )
223          ALLOCATE( local_pf(mask_size_l(mid,1),mask_size_l(mid,2), &
224                             mask_size_l(mid,3)) )
225       ENDIF
226!
227!--    Set masking flag for topography for not resorted arrays
228       flag_nr = 0
229!
230!--    Store the variable chosen.
231       resorted = .FALSE.
232       SELECT CASE ( TRIM( domask(mid,av,ivar) ) )
233
234          CASE ( 'e' )
235             IF ( av == 0 )  THEN
236                to_be_resorted => e
237             ELSE
238                to_be_resorted => e_av
239             ENDIF
240
241          CASE ( 'thetal' )
242             IF ( av == 0 )  THEN
243                to_be_resorted => pt
244             ELSE
245                to_be_resorted => lpt_av
246             ENDIF
247
248          CASE ( 'nc' )
249             IF ( av == 0 )  THEN
250                to_be_resorted => nc
251             ELSE
252                to_be_resorted => nc_av
253             ENDIF
254
255          CASE ( 'nr' )
256             IF ( av == 0 )  THEN
257                to_be_resorted => nr
258             ELSE
259                to_be_resorted => nr_av
260             ENDIF
261
262          CASE ( 'p' )
263             IF ( av == 0 )  THEN
264                to_be_resorted => p
265             ELSE
266                to_be_resorted => p_av
267             ENDIF
268
269          CASE ( 'pc' )  ! particle concentration (requires ghostpoint exchange)
270             IF ( av == 0 )  THEN
271                tend = prt_count
272                CALL exchange_horiz( tend, nbgp )
273                IF ( .NOT. mask_surface(mid) )  THEN
274                   DO  i = 1, mask_size_l(mid,1)
275                      DO  j = 1, mask_size_l(mid,2)
276                         DO  k = 1, mask_size_l(mid,3)
277                            local_pf(i,j,k) =  tend(mask_k(mid,k), &
278                                      mask_j(mid,j),mask_i(mid,i))
279                         ENDDO
280                      ENDDO
281                   ENDDO
282                ELSE
283!
284!--                Terrain-following masked output
285                   DO  i = 1, mask_size_l(mid,1)
286                      DO  j = 1, mask_size_l(mid,2)
287!--                      Get k index of the highest terraing surface
288                         im = mask_i(mid,i)
289                         jm = mask_j(mid,j)
290                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )),&
291                                       DIM = 1 ) - 1
292                         DO  k = 1, mask_size_l(mid,3)
293                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
294!--                         Set value if not in building
295                            IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
296                               local_pf(i,j,k) = fill_value
297                            ELSE
298                               local_pf(i,j,k) =  tend(kk,jm,im)
299                            ENDIF
300                         ENDDO
301                      ENDDO
302                   ENDDO
303                ENDIF
304                resorted = .TRUE.
305             ELSE
306                CALL exchange_horiz( pc_av, nbgp )
307                to_be_resorted => pc_av
308             ENDIF
309
310          CASE ( 'pr' )  ! mean particle radius (effective radius)
311             IF ( av == 0 )  THEN
312                IF ( time_since_reference_point >= particle_advection_start )  THEN
313                   DO  i = nxl, nxr
314                      DO  j = nys, nyn
315                         DO  k = nzb, nz_do3d
316                            number_of_particles = prt_count(k,j,i)
317                            IF (number_of_particles <= 0)  CYCLE
318                            particles => grid_particles(k,j,i)%particles(1:number_of_particles)
319                            s_r2 = 0.0_wp
320                            s_r3 = 0.0_wp
321                            DO  n = 1, number_of_particles
322                               IF ( particles(n)%particle_mask )  THEN
323                                  s_r2 = s_r2 + grid_particles(k,j,i)%particles(n)%radius**2 * &
324                                         grid_particles(k,j,i)%particles(n)%weight_factor
325                                  s_r3 = s_r3 + grid_particles(k,j,i)%particles(n)%radius**3 * &
326                                         grid_particles(k,j,i)%particles(n)%weight_factor
327                               ENDIF
328                            ENDDO
329                            IF ( s_r2 > 0.0_wp )  THEN
330                               mean_r = s_r3 / s_r2
331                            ELSE
332                               mean_r = 0.0_wp
333                            ENDIF
334                            tend(k,j,i) = mean_r
335                         ENDDO
336                      ENDDO
337                   ENDDO
338                   CALL exchange_horiz( tend, nbgp )
339                ELSE
340                   tend = 0.0_wp
341                ENDIF
342                IF ( .NOT. mask_surface(mid) )  THEN
343                   DO  i = 1, mask_size_l(mid,1)
344                      DO  j = 1, mask_size_l(mid,2)
345                         DO  k = 1, mask_size_l(mid,3)
346                            local_pf(i,j,k) =  tend(mask_k(mid,k), &
347                                      mask_j(mid,j),mask_i(mid,i))
348                         ENDDO
349                      ENDDO
350                   ENDDO
351                ELSE
352!
353!--                Terrain-following masked output
354                   DO  i = 1, mask_size_l(mid,1)
355                      DO  j = 1, mask_size_l(mid,2)
356!--                      Get k index of the highest terraing surface
357                         im = mask_i(mid,i)
358                         jm = mask_j(mid,j)
359                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
360                                         DIM = 1 ) - 1
361                         DO  k = 1, mask_size_l(mid,3)
362                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
363!--                         Set value if not in building
364                            IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
365                               local_pf(i,j,k) = fill_value
366                            ELSE
367                               local_pf(i,j,k) =  tend(kk,jm,im)
368                            ENDIF
369                         ENDDO
370                      ENDDO
371                   ENDDO
372                ENDIF
373                resorted = .TRUE.
374             ELSE
375                CALL exchange_horiz( pr_av, nbgp )
376                to_be_resorted => pr_av
377             ENDIF
378
379          CASE ( 'theta' )
380             IF ( av == 0 )  THEN
381                IF ( .NOT. bulk_cloud_model ) THEN
382                   to_be_resorted => pt
383                ELSE
384                   IF ( .NOT. mask_surface(mid) )  THEN
385                      DO  i = 1, mask_size_l(mid,1)
386                         DO  j = 1, mask_size_l(mid,2)
387                            DO  k = 1, mask_size_l(mid,3)
388                               local_pf(i,j,k) =  &
389                                  pt(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) &
390                                  + lv_d_cp * d_exner(mask_k(mid,k)) *          &
391                                    ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
392                            ENDDO
393                         ENDDO
394                      ENDDO
395                   ELSE
396!
397!--                   Terrain-following masked output
398                      DO  i = 1, mask_size_l(mid,1)
399                         DO  j = 1, mask_size_l(mid,2)
400!--                         Get k index of the highest terraing surface
401                            im = mask_i(mid,i)
402                            jm = mask_j(mid,j)
403                            ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
404                                             DIM = 1 ) - 1
405                            DO  k = 1, mask_size_l(mid,3)
406                               kk = MIN( ktt+mask_k(mid,k), nzt+1 )
407!--                            Set value if not in building
408                               IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
409                                  local_pf(i,j,k) = fill_value
410                               ELSE
411                                  local_pf(i,j,k) = pt(kk,jm,im) + lv_d_cp * d_exner(kk) * ql(kk,jm,im)
412                               ENDIF
413                            ENDDO
414                         ENDDO
415                      ENDDO
416                   ENDIF
417                   resorted = .TRUE.
418                ENDIF
419             ELSE
420                to_be_resorted => pt_av
421             ENDIF
422
423          CASE ( 'q' )
424             IF ( av == 0 )  THEN
425                to_be_resorted => q
426             ELSE
427                to_be_resorted => q_av
428             ENDIF
429
430          CASE ( 'qc' )
431             IF ( av == 0 )  THEN
432                to_be_resorted => qc
433             ELSE
434                to_be_resorted => qc_av
435             ENDIF
436
437          CASE ( 'ql' )
438             IF ( av == 0 )  THEN
439                to_be_resorted => ql
440             ELSE
441                to_be_resorted => ql_av
442             ENDIF
443
444          CASE ( 'ql_c' )
445             IF ( av == 0 )  THEN
446                to_be_resorted => ql_c
447             ELSE
448                to_be_resorted => ql_c_av
449             ENDIF
450
451          CASE ( 'ql_v' )
452             IF ( av == 0 )  THEN
453                to_be_resorted => ql_v
454             ELSE
455                to_be_resorted => ql_v_av
456             ENDIF
457
458          CASE ( 'ql_vp' )
459             IF ( av == 0 )  THEN
460                IF ( time_since_reference_point >= particle_advection_start )  THEN
461                   DO  i = nxl, nxr
462                      DO  j = nys, nyn
463                         DO  k = nzb, nz_do3d
464                            number_of_particles = prt_count(k,j,i)
465                            IF (number_of_particles <= 0)  CYCLE
466                            particles => grid_particles(k,j,i)%particles(1:number_of_particles)
467                            DO  n = 1, number_of_particles
468                               IF ( particles(n)%particle_mask )  THEN
469                                  tend(k,j,i) = tend(k,j,i) + &
470                                          particles(n)%weight_factor / &
471                                          prt_count(k,j,i)
472                               ENDIF
473                            ENDDO
474                         ENDDO
475                      ENDDO
476                   ENDDO
477                   CALL exchange_horiz( tend, nbgp )
478                ELSE
479                   tend = 0.0_wp
480                ENDIF
481                IF ( .NOT. mask_surface(mid) )  THEN
482                   DO  i = 1, mask_size_l(mid,1)
483                      DO  j = 1, mask_size_l(mid,2)
484                         DO  k = 1, mask_size_l(mid,3)
485                            local_pf(i,j,k) =  tend(mask_k(mid,k), &
486                                      mask_j(mid,j),mask_i(mid,i))
487                         ENDDO
488                      ENDDO
489                   ENDDO
490                ELSE
491!
492!--                Terrain-following masked output
493                   DO  i = 1, mask_size_l(mid,1)
494                      DO  j = 1, mask_size_l(mid,2)
495!--                      Get k index of the highest terraing surface
496                         im = mask_i(mid,i)
497                         jm = mask_j(mid,j)
498                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
499                                          DIM = 1 ) - 1
500                         DO  k = 1, mask_size_l(mid,3)
501                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
502!--                         Set value if not in building
503                            IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
504                               local_pf(i,j,k) = fill_value
505                            ELSE
506                               local_pf(i,j,k) = tend(kk,jm,im)
507                            ENDIF
508                         ENDDO
509                      ENDDO
510                   ENDDO
511                ENDIF
512                resorted = .TRUE.
513             ELSE
514                CALL exchange_horiz( ql_vp_av, nbgp )
515                to_be_resorted => ql_vp_av
516             ENDIF
517
518          CASE ( 'qv' )
519             IF ( av == 0 )  THEN
520                IF ( .NOT. mask_surface(mid) )  THEN
521                   DO  i = 1, mask_size_l(mid,1)
522                      DO  j = 1, mask_size_l(mid,2)
523                         DO  k = 1, mask_size_l(mid,3)
524                            local_pf(i,j,k) =  &
525                                 q(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) -  &
526                                 ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
527                         ENDDO
528                      ENDDO
529                   ENDDO
530                ELSE
531!
532!--                Terrain-following masked output
533                   DO  i = 1, mask_size_l(mid,1)
534                      DO  j = 1, mask_size_l(mid,2)
535!--                      Get k index of the highest terraing surface
536                         im = mask_i(mid,i)
537                         jm = mask_j(mid,j)
538                         ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
539                                          DIM = 1 ) - 1
540                         DO  k = 1, mask_size_l(mid,3)
541                            kk = MIN( ktt+mask_k(mid,k), nzt+1 )
542!--                         Set value if not in building
543                            IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
544                               local_pf(i,j,k) = fill_value
545                            ELSE
546                               local_pf(i,j,k) = q(kk,jm,im) - ql(kk,jm,im)
547                            ENDIF
548                         ENDDO
549                      ENDDO
550                   ENDDO
551                ENDIF
552                resorted = .TRUE.
553             ELSE
554                to_be_resorted => qv_av
555             ENDIF
556
557          CASE ( 'qr' )
558             IF ( av == 0 )  THEN
559                to_be_resorted => qr
560             ELSE
561                to_be_resorted => qr_av
562             ENDIF
563
564          CASE ( 'rho_sea_water' )
565             IF ( av == 0 )  THEN
566                to_be_resorted => rho_ocean
567             ELSE
568                to_be_resorted => rho_ocean_av
569             ENDIF
570
571          CASE ( 's' )
572             IF ( av == 0 )  THEN
573                to_be_resorted => s
574             ELSE
575                to_be_resorted => s_av
576             ENDIF
577
578          CASE ( 'sa' )
579             IF ( av == 0 )  THEN
580                to_be_resorted => sa
581             ELSE
582                to_be_resorted => sa_av
583             ENDIF
584
585          CASE ( 'u' )
586             flag_nr = 1
587             IF ( av == 0 )  THEN
588                to_be_resorted => u
589             ELSE
590                to_be_resorted => u_av
591             ENDIF
592
593          CASE ( 'v' )
594             flag_nr = 2
595             IF ( av == 0 )  THEN
596                to_be_resorted => v
597             ELSE
598                to_be_resorted => v_av
599             ENDIF
600
601          CASE ( 'thetav' )
602             IF ( av == 0 )  THEN
603                to_be_resorted => vpt
604             ELSE
605                to_be_resorted => vpt_av
606             ENDIF
607
608          CASE ( 'w' )
609             flag_nr = 3
610             IF ( av == 0 )  THEN
611                to_be_resorted => w
612             ELSE
613                to_be_resorted => w_av
614             ENDIF
615
616          CASE DEFAULT
617!
618!--          Set flag to steer output of radiation, land-surface, or user-defined
619!--          quantities
620             found = .FALSE.
621!
622!--          Radiation quantity
623             IF ( .NOT. found  .AND. radiation )  THEN
624                CALL radiation_data_output_mask(av, domask(mid,av,ivar), found,&
625                                                local_pf, mid )
626             ENDIF
627
628             IF ( .NOT. found  .AND. air_chemistry )  THEN
629                CALL chem_data_output_mask(av, domask(mid,av,ivar), found,     &
630                                           local_pf, mid )
631             ENDIF
632!
633!--          Check for diagnostic quantities
634             IF ( .NOT. found )  THEN
635                CALL doq_output_mask( av, domask(mid,av,ivar), found, local_pf,   &
636                                      mid)
637             ENDIF
638!
639!--          SALSA quantities
640             IF ( .NOT. found .AND. salsa )  THEN
641                CALL salsa_data_output_mask( av, domask(mid,av,ivar), found,   &
642                                             local_pf, mid )
643             ENDIF
644!
645!--          User defined quantity
646             IF ( .NOT. found )  THEN
647                CALL user_data_output_mask(av, domask(mid,av,ivar), found,     &
648                                           local_pf, mid )
649             ENDIF
650
651             resorted = .TRUE.
652
653             IF ( .NOT. found )  THEN
654                WRITE ( message_string, * ) 'no masked output available for: ',&
655                                            TRIM( domask(mid,av,ivar) )
656                CALL message( 'data_output_mask', 'PA0327', 0, 0, 0, 6, 0 )
657             ENDIF
658
659       END SELECT
660
661!
662!--    Resort the array to be output, if not done above
663       IF ( .NOT. resorted )  THEN
664          IF ( .NOT. mask_surface(mid) )  THEN
665!
666!--          Default masked output
667             DO  i = 1, mask_size_l(mid,1)
668                DO  j = 1, mask_size_l(mid,2)
669                   DO  k = 1, mask_size_l(mid,3)
670                      local_pf(i,j,k) = MERGE( to_be_resorted(mask_k(mid,k),  &
671                                                              mask_j(mid,j),  &
672                                                              mask_i(mid,i)), &
673                                               REAL( fill_value, KIND = wp ), &
674                                               BTEST( wall_flags_total_0(     &
675                                                              mask_k(mid,k),  &
676                                                              mask_j(mid,j),  &
677                                                              mask_i(mid,i)), &
678                                                      flag_nr ) )
679                   ENDDO
680                ENDDO
681             ENDDO
682
683          ELSE
684!
685!--          Terrain-following masked output
686             DO  i = 1, mask_size_l(mid,1)
687                DO  j = 1, mask_size_l(mid,2)
688!--                Get k index of the highest terraing surface
689                   im = mask_i(mid,i)
690                   jm = mask_j(mid,j)
691                   ktt = MINLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,jm,im), 5 )), &
692                                    DIM = 1 ) - 1
693                   DO  k = 1, mask_size_l(mid,3)
694                      kk = MIN( ktt+mask_k(mid,k), nzt+1 )
695!--                   Set value if not in building
696                      IF ( BTEST( wall_flags_total_0(kk,jm,im), 6 ) )  THEN
697                         local_pf(i,j,k) = fill_value
698                      ELSE
699                         local_pf(i,j,k) = to_be_resorted(kk,jm,im)
700                      ENDIF
701                   ENDDO
702                ENDDO
703             ENDDO
704
705          ENDIF
706       ENDIF
707
708!
709!--    I/O block. I/O methods are implemented
710!--    (1) for parallel execution
711!--     a. with netCDF 4 parallel I/O-enabled library
712!--     b. with netCDF 3 library
713!--    (2) for serial execution.
714!--    The choice of method depends on the correct setting of preprocessor
715!--    directives __parallel and __netcdf4_parallel as well as on the parameter
716!--    netcdf_data_format.
717#if defined( __parallel )
718#if defined( __netcdf4_parallel )
719       IF ( netcdf_data_format > 4 )  THEN
720!
721!--       (1) a. Parallel I/O using netCDF 4 (not yet tested)
722          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),                         &
723               id_var_domask(mid,av,ivar), local_pf,                           &
724               start = (/ mask_start_l(mid,1), mask_start_l(mid,2),            &
725                          mask_start_l(mid,3), domask_time_count(mid,av) /),   &
726               count = (/ mask_size_l(mid,1), mask_size_l(mid,2),              &
727                          mask_size_l(mid,3), 1 /) )
728          CALL netcdf_handle_error( 'data_output_mask', 461 )
729       ELSE
730#endif
731!
732!--       (1) b. Conventional I/O only through PE0
733!--       PE0 receives partial arrays from all processors of the respective mask
734!--       and outputs them. Here a barrier has to be set, because otherwise
735!--       "-MPI- FATAL: Remote protocol queue full" may occur.
736          CALL MPI_BARRIER( comm2d, ierr )
737
738          ngp = mask_size_l(mid,1) * mask_size_l(mid,2) * mask_size_l(mid,3)
739          IF ( myid == 0 )  THEN
740!
741!--          Local array can be relocated directly.
742             total_pf( &
743               mask_start_l(mid,1):mask_start_l(mid,1)+mask_size_l(mid,1)-1, &
744               mask_start_l(mid,2):mask_start_l(mid,2)+mask_size_l(mid,2)-1, &
745               mask_start_l(mid,3):mask_start_l(mid,3)+mask_size_l(mid,3)-1 ) &
746               = local_pf
747!
748!--          Receive data from all other PEs.
749             DO  n = 1, numprocs-1
750!
751!--             Receive index limits first, then array.
752!--             Index limits are received in arbitrary order from the PEs.
753                CALL MPI_RECV( ind(1), 6, MPI_INTEGER, MPI_ANY_SOURCE, 0,  &
754                     comm2d, status, ierr )
755!
756!--             Not all PEs have data for the mask
757                IF ( ind(1) /= -9999 )  THEN
758                   ngp = ( ind(2)-ind(1)+1 ) * (ind(4)-ind(3)+1 ) *  &
759                         ( ind(6)-ind(5)+1 )
760                   sender = status(MPI_SOURCE)
761                   DEALLOCATE( local_pf )
762                   ALLOCATE(local_pf(ind(1):ind(2),ind(3):ind(4),ind(5):ind(6)))
763                   CALL MPI_RECV( local_pf(ind(1),ind(3),ind(5)), ngp,  &
764                        MPI_REAL, sender, 1, comm2d, status, ierr )
765                   total_pf(ind(1):ind(2),ind(3):ind(4),ind(5):ind(6)) &
766                        = local_pf
767                ENDIF
768             ENDDO
769
770             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),                      &
771                  id_var_domask(mid,av,ivar), total_pf,                        &
772                  start = (/ 1, 1, 1, domask_time_count(mid,av) /),            &
773                  count = (/ mask_size(mid,1), mask_size(mid,2),               &
774                             mask_size(mid,3), 1 /) )
775             CALL netcdf_handle_error( 'data_output_mask', 462 )
776
777          ELSE
778!
779!--          If at least part of the mask resides on the PE, send the index
780!--          limits for the target array, otherwise send -9999 to PE0.
781             IF ( mask_size_l(mid,1) > 0 .AND.  mask_size_l(mid,2) > 0 .AND. &
782                  mask_size_l(mid,3) > 0  ) &
783                  THEN
784                ind(1) = mask_start_l(mid,1)
785                ind(2) = mask_start_l(mid,1) + mask_size_l(mid,1) - 1
786                ind(3) = mask_start_l(mid,2)
787                ind(4) = mask_start_l(mid,2) + mask_size_l(mid,2) - 1
788                ind(5) = mask_start_l(mid,3)
789                ind(6) = mask_start_l(mid,3) + mask_size_l(mid,3) - 1
790             ELSE
791                ind(1) = -9999; ind(2) = -9999
792                ind(3) = -9999; ind(4) = -9999
793                ind(5) = -9999; ind(6) = -9999
794             ENDIF
795             CALL MPI_SEND( ind(1), 6, MPI_INTEGER, 0, 0, comm2d, ierr )
796!
797!--          If applicable, send data to PE0.
798             IF ( ind(1) /= -9999 )  THEN
799                CALL MPI_SEND( local_pf(1,1,1), ngp, MPI_REAL, 0, 1, comm2d, &
800                     ierr )
801             ENDIF
802          ENDIF
803!
804!--       A barrier has to be set, because otherwise some PEs may proceed too
805!--       fast so that PE0 may receive wrong data on tag 0.
806          CALL MPI_BARRIER( comm2d, ierr )
807#if defined( __netcdf4_parallel )
808       ENDIF
809#endif
810#else
811!
812!--    (2) For serial execution of PALM, the single processor (PE0) holds all
813!--    data and writes them directly to file.
814       nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),                            &
815                               id_var_domask(mid,av,ivar), local_pf,           &
816                             start = (/ 1, 1, 1, domask_time_count(mid,av) /), &
817                             count = (/ mask_size_l(mid,1), mask_size_l(mid,2),&
818                               mask_size_l(mid,3), 1 /) )
819       CALL netcdf_handle_error( 'data_output_mask', 463 )
820#endif
821
822       ivar = ivar + 1
823
824    ENDDO
825
826!
827!-- Deallocate temporary arrays.
828    DEALLOCATE( local_pf )
829#if defined( __parallel )
830    IF ( myid == 0 )  THEN
831       DEALLOCATE( total_pf )
832    ENDIF
833#endif
834
835!
836!-- Switch back to original format given by user (see beginning of this routine)
837    netcdf_data_format = netcdf_data_format_save
838
839    CALL cpu_log( log_point(49), 'data_output_mask', 'stop' )
840#endif
841
842
843 END SUBROUTINE data_output_mask
Note: See TracBrowser for help on using the repository browser.