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

Last change on this file since 1007 was 1007, checked in by franke, 12 years ago

Bugfixes


Missing calculation of mean particle weighting factor for output added. (data_output_2d, data_output_3d, data_output_mask, sum_up_3d_data)
Calculation of mean particle radius for output now considers the weighting factor. (data_output_mask)
Calculation of sugrid-scale buoyancy flux for humidity and cloud droplets corrected. (flow_statistics)
Factor in calculation of enhancement factor for collision efficencies corrected. (lpm_collision_kernels)
Calculation of buoyancy production now considers the liquid water mixing ratio in case of cloud droplets. (production_e)

Changes


Calculation of buoyancy flux for humidity in case of WS-scheme is now using turbulent fluxes of WS-scheme. (flow_statistics)
Calculation of the collision kernels now in SI units. (lpm_collision_kernels)

  • Property svn:keywords set to Id
File size: 16.6 KB
Line 
1 SUBROUTINE data_output_mask( av )
2
3!------------------------------------------------------------------------------!
4! Current revisions:
5! -----------------
6! Bugfix: calculation of pr must depend on the particle weighting factor,
7! missing calculation of ql_vp added
8!
9! Former revisions:
10! -----------------
11! $Id: data_output_mask.f90 1007 2012-09-19 14:30:36Z franke $
12!
13! 771 2011-10-27 10:56:21Z heinze
14! +lpt
15!
16! 667 2010-12-23 12:06:00Z suehring/gryschka
17! Calls of exchange_horiz are modified.
18!
19! 564 2010-09-30 13:18:59Z helmke
20! start number of mask output files changed to 201, netcdf message identifiers
21! of masked output changed, palm message identifiers of masked output changed
22!
23! 493 2010-03-01 08:30:24Z raasch
24! netcdf_format_mask* and format_parallel_io replaced by netcdf_data_format
25!
26! 475 2010-02-04 02:26:16Z raasch
27! Bugfix in serial branch: arguments from array local_pf removed in N90_PUT_VAR
28!
29! 410 2009-12-04 17:05:40Z letzel
30! Initial version
31!
32! Description:
33! ------------
34! Masked data output in NetCDF format for current mask (current value of mid).
35!------------------------------------------------------------------------------!
36
37#if defined( __netcdf )
38    USE arrays_3d
39    USE averaging
40    USE cloud_parameters
41    USE control_parameters
42    USE cpulog
43    USE grid_variables
44    USE indices
45    USE interfaces
46    USE netcdf
47    USE netcdf_control
48    USE particle_attributes
49    USE pegrid
50
51    IMPLICIT NONE
52
53    INTEGER ::  av, ngp, file_id, i, if, is, j, k, l, n, psi, s, sender, &
54                ind(6)
55    LOGICAL ::  found, resorted
56    REAL    ::  mean_r, s_r3, s_r4
57    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  local_pf
58#if defined( __parallel )
59    REAL, DIMENSION(:,:,:), ALLOCATABLE ::  total_pf
60#endif
61    REAL, DIMENSION(:,:,:), POINTER ::  to_be_resorted
62
63!
64!-- Return, if nothing to output
65    IF ( domask_no(mid,av) == 0 )  RETURN
66
67    CALL cpu_log (log_point(49),'data_output_mask','start')
68
69!
70!-- Open output file.
71    IF ( netcdf_output  .AND.  ( myid == 0  .OR.  netcdf_data_format > 2 ) ) &
72    THEN
73       CALL check_open( 200+mid+av*max_masks )
74    ENDIF 
75
76!
77!-- Allocate total and local output arrays.
78#if defined( __parallel )
79    IF ( myid == 0 )  THEN
80       ALLOCATE( total_pf(mask_size(mid,1),mask_size(mid,2),mask_size(mid,3)) )
81    ENDIF
82#endif
83    ALLOCATE( local_pf(mask_size_l(mid,1),mask_size_l(mid,2), &
84                       mask_size_l(mid,3)) )
85
86!
87!-- Update the NetCDF time axis.
88    domask_time_count(mid,av) = domask_time_count(mid,av) + 1
89    IF ( netcdf_output  .AND.  ( myid == 0  .OR.  netcdf_data_format > 2 ) ) &
90    THEN
91       nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_time_mask(mid,av), &
92                               (/ simulated_time /),                          &
93                               start = (/ domask_time_count(mid,av) /),       &
94                               count = (/ 1 /) )
95       CALL handle_netcdf_error( 'data_output_mask', 460 )
96    ENDIF
97
98!
99!-- Loop over all variables to be written.
100    if = 1
101
102    DO  WHILE ( domask(mid,av,if)(1:1) /= ' ' )
103!
104!--    Reallocate local_pf on PE 0 since its shape changes during MPI exchange
105       IF ( netcdf_data_format < 3   .AND.  myid == 0  .AND.  if > 1 )  THEN
106          DEALLOCATE( local_pf )
107          ALLOCATE( local_pf(mask_size_l(mid,1),mask_size_l(mid,2), &
108                             mask_size_l(mid,3)) )
109       ENDIF
110!
111!--    Store the variable chosen.
112       resorted = .FALSE.
113       SELECT CASE ( TRIM( domask(mid,av,if) ) )
114
115          CASE ( 'e' )
116             IF ( av == 0 )  THEN
117                to_be_resorted => e
118             ELSE
119                to_be_resorted => e_av
120             ENDIF
121
122          CASE ( 'lpt' )
123             IF ( av == 0 )  THEN
124                to_be_resorted => pt
125             ELSE
126                to_be_resorted => lpt_av
127             ENDIF
128
129          CASE ( 'p' )
130             IF ( av == 0 )  THEN
131                to_be_resorted => p
132             ELSE
133                to_be_resorted => p_av
134             ENDIF
135
136          CASE ( 'pc' )  ! particle concentration (requires ghostpoint exchange)
137             IF ( av == 0 )  THEN
138                tend = prt_count
139                CALL exchange_horiz( tend, nbgp )
140                DO  i = 1, mask_size_l(mid,1)
141                   DO  j = 1, mask_size_l(mid,2)
142                      DO  k = 1, mask_size_l(mid,3)
143                         local_pf(i,j,k) =  tend(mask_k(mid,k), &
144                                   mask_j(mid,j),mask_i(mid,i))
145                      ENDDO
146                   ENDDO
147                ENDDO
148                resorted = .TRUE.
149             ELSE
150                CALL exchange_horiz( pc_av, nbgp )
151                to_be_resorted => pc_av
152             ENDIF
153
154          CASE ( 'pr' )  ! mean particle radius
155             IF ( av == 0 )  THEN
156                DO  i = nxl, nxr
157                   DO  j = nys, nyn
158                      DO  k = nzb, nzt+1
159                         psi = prt_start_index(k,j,i)
160                         s_r3 = 0.0
161                         s_r4 = 0.0
162                         DO  n = psi, psi+prt_count(k,j,i)-1
163                            s_r3 = s_r3 + particles(n)%radius**3 * &
164                                          particles(n)%weight_factor
165                            s_r4 = s_r4 + particles(n)%radius**4 * &
166                                          particles(n)%weight_factor
167                         ENDDO
168                         IF ( s_r3 /= 0.0 )  THEN
169                            mean_r = s_r4 / s_r3
170                         ELSE
171                            mean_r = 0.0
172                         ENDIF
173                         tend(k,j,i) = mean_r
174                      ENDDO
175                   ENDDO
176                ENDDO
177                CALL exchange_horiz( tend, nbgp )
178                DO  i = 1, mask_size_l(mid,1)
179                   DO  j = 1, mask_size_l(mid,2)
180                      DO  k = 1, mask_size_l(mid,3)
181                         local_pf(i,j,k) =  tend(mask_k(mid,k), &
182                                   mask_j(mid,j),mask_i(mid,i))
183                      ENDDO
184                   ENDDO
185                ENDDO
186                resorted = .TRUE.
187             ELSE
188                CALL exchange_horiz( pr_av, nbgp )
189                to_be_resorted => pr_av
190             ENDIF
191
192          CASE ( 'pt' )
193             IF ( av == 0 )  THEN
194                IF ( .NOT. cloud_physics ) THEN
195                   to_be_resorted => pt
196                ELSE
197                   DO  i = 1, mask_size_l(mid,1)
198                      DO  j = 1, mask_size_l(mid,2)
199                         DO  k = 1, mask_size_l(mid,3)
200                            local_pf(i,j,k) =  &
201                                 pt(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) &
202                                 + l_d_cp * pt_d_t(mask_k(mid,k)) * &
203                                   ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
204                         ENDDO
205                      ENDDO
206                   ENDDO
207                   resorted = .TRUE.
208                ENDIF
209             ELSE
210                to_be_resorted => pt_av
211             ENDIF
212
213          CASE ( 'q' )
214             IF ( av == 0 )  THEN
215                to_be_resorted => q
216             ELSE
217                to_be_resorted => q_av
218             ENDIF
219
220          CASE ( 'ql' )
221             IF ( av == 0 )  THEN
222                to_be_resorted => ql
223             ELSE
224                to_be_resorted => ql_av
225             ENDIF
226
227          CASE ( 'ql_c' )
228             IF ( av == 0 )  THEN
229                to_be_resorted => ql_c
230             ELSE
231                to_be_resorted => ql_c_av
232             ENDIF
233
234          CASE ( 'ql_v' )
235             IF ( av == 0 )  THEN
236                to_be_resorted => ql_v
237             ELSE
238                to_be_resorted => ql_v_av
239             ENDIF
240
241          CASE ( 'ql_vp' )
242             IF ( av == 0 )  THEN
243                DO  i = nxl, nxr
244                   DO  j = nys, nyn
245                      DO  k = nzb, nz_do3d
246                         psi = prt_start_index(k,j,i)
247                         DO  n = psi, psi+prt_count(k,j,i)-1
248                            tend(k,j,i) = tend(k,j,i) + &
249                                          particles(n)%weight_factor / &
250                                          prt_count(k,j,i)
251                         ENDDO
252                      ENDDO
253                   ENDDO
254                ENDDO
255                CALL exchange_horiz( tend, nbgp )
256                DO  i = 1, mask_size_l(mid,1)
257                   DO  j = 1, mask_size_l(mid,2)
258                      DO  k = 1, mask_size_l(mid,3)
259                         local_pf(i,j,k) =  tend(mask_k(mid,k), &
260                                   mask_j(mid,j),mask_i(mid,i))
261                      ENDDO
262                   ENDDO
263                ENDDO
264                resorted = .TRUE.
265             ELSE
266                CALL exchange_horiz( ql_vp_av, nbgp )
267                to_be_resorted => ql_vp_av
268             ENDIF
269
270          CASE ( 'qv' )
271             IF ( av == 0 )  THEN
272                DO  i = 1, mask_size_l(mid,1)
273                   DO  j = 1, mask_size_l(mid,2)
274                      DO  k = 1, mask_size_l(mid,3)
275                         local_pf(i,j,k) =  &
276                              q(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i)) -  &
277                              ql(mask_k(mid,k),mask_j(mid,j),mask_i(mid,i))
278                      ENDDO
279                   ENDDO
280                ENDDO
281                resorted = .TRUE.
282             ELSE
283                to_be_resorted => qv_av
284             ENDIF
285
286          CASE ( 'rho' )
287             IF ( av == 0 )  THEN
288                to_be_resorted => rho
289             ELSE
290                to_be_resorted => rho_av
291             ENDIF
292
293          CASE ( 's' )
294             IF ( av == 0 )  THEN
295                to_be_resorted => q
296             ELSE
297                to_be_resorted => s_av
298             ENDIF
299
300          CASE ( 'sa' )
301             IF ( av == 0 )  THEN
302                to_be_resorted => sa
303             ELSE
304                to_be_resorted => sa_av
305             ENDIF
306
307          CASE ( 'u' )
308             IF ( av == 0 )  THEN
309                to_be_resorted => u
310             ELSE
311                to_be_resorted => u_av
312             ENDIF
313
314          CASE ( 'v' )
315             IF ( av == 0 )  THEN
316                to_be_resorted => v
317             ELSE
318                to_be_resorted => v_av
319             ENDIF
320
321          CASE ( 'vpt' )
322             IF ( av == 0 )  THEN
323                to_be_resorted => vpt
324             ELSE
325                to_be_resorted => vpt_av
326             ENDIF
327
328          CASE ( 'w' )
329             IF ( av == 0 )  THEN
330                to_be_resorted => w
331             ELSE
332                to_be_resorted => w_av
333             ENDIF
334
335          CASE DEFAULT
336!
337!--          User defined quantity
338             CALL user_data_output_mask(av, domask(mid,av,if), found, local_pf )
339             resorted = .TRUE.
340
341             IF ( .NOT. found )  THEN
342                WRITE ( message_string, * ) 'no output available for: ', &
343                                            TRIM( domask(mid,av,if) )
344                CALL message( 'data_output_mask', 'PA0327', 0, 0, 0, 6, 0 )
345             ENDIF
346
347       END SELECT
348
349!
350!--    Resort the array to be output, if not done above
351       IF ( .NOT. resorted )  THEN
352          DO  i = 1, mask_size_l(mid,1)
353             DO  j = 1, mask_size_l(mid,2)
354                DO  k = 1, mask_size_l(mid,3)
355                   local_pf(i,j,k) =  to_be_resorted(mask_k(mid,k), &
356                                      mask_j(mid,j),mask_i(mid,i))
357                ENDDO
358             ENDDO
359          ENDDO
360       ENDIF
361
362!
363!--    I/O block. I/O methods are implemented
364!--    (1) for parallel execution
365!--     a. with NetCDF 4 parallel I/O-enabled library
366!--     b. with NetCDF 3 library
367!--    (2) for serial execution.
368!--    The choice of method depends on the correct setting of preprocessor
369!--    directives __parallel and __netcdf4 as well as on the parameter
370!--    netcdf_data_format.
371#if defined( __parallel )
372#if defined( __netcdf4 )
373       IF ( netcdf_data_format > 2 )  THEN
374!
375!--       (1) a. Parallel I/O using NetCDF 4 (not yet tested)
376          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),  &
377               id_var_domask(mid,av,if),  &
378               local_pf,  &
379               start = (/ mask_start_l(mid,1), mask_start_l(mid,2),  &
380                          mask_start_l(mid,3), domask_time_count(mid,av) /),  &
381               count = (/ mask_size_l(mid,1), mask_size_l(mid,2),  &
382                          mask_size_l(mid,3), 1 /) )
383          CALL handle_netcdf_error( 'data_output_mask', 461 )
384       ELSE
385#endif
386!
387!--       (1) b. Conventional I/O only through PE0
388!--       PE0 receives partial arrays from all processors of the respective mask
389!--       and outputs them. Here a barrier has to be set, because otherwise
390!--       "-MPI- FATAL: Remote protocol queue full" may occur.
391          CALL MPI_BARRIER( comm2d, ierr )
392
393          ngp = mask_size_l(mid,1) * mask_size_l(mid,2) * mask_size_l(mid,3)
394          IF ( myid == 0 )  THEN
395!
396!--          Local array can be relocated directly.
397             total_pf( &
398               mask_start_l(mid,1):mask_start_l(mid,1)+mask_size_l(mid,1)-1, &
399               mask_start_l(mid,2):mask_start_l(mid,2)+mask_size_l(mid,2)-1, &
400               mask_start_l(mid,3):mask_start_l(mid,3)+mask_size_l(mid,3)-1 ) &
401               = local_pf
402!
403!--          Receive data from all other PEs.
404             DO  n = 1, numprocs-1
405!
406!--             Receive index limits first, then array.
407!--             Index limits are received in arbitrary order from the PEs.
408                CALL MPI_RECV( ind(1), 6, MPI_INTEGER, MPI_ANY_SOURCE, 0,  &
409                     comm2d, status, ierr )
410!
411!--             Not all PEs have data for the mask
412                IF ( ind(1) /= -9999 )  THEN
413                   ngp = ( ind(2)-ind(1)+1 ) * (ind(4)-ind(3)+1 ) *  &
414                         ( ind(6)-ind(5)+1 )
415                   sender = status(MPI_SOURCE)
416                   DEALLOCATE( local_pf )
417                   ALLOCATE(local_pf(ind(1):ind(2),ind(3):ind(4),ind(5):ind(6)))
418                   CALL MPI_RECV( local_pf(ind(1),ind(3),ind(5)), ngp,  &
419                        MPI_REAL, sender, 1, comm2d, status, ierr )
420                   total_pf(ind(1):ind(2),ind(3):ind(4),ind(5):ind(6)) &
421                        = local_pf
422                ENDIF
423             ENDDO
424
425             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),  &
426                  id_var_domask(mid,av,if), total_pf, &
427                  start = (/ 1, 1, 1, domask_time_count(mid,av) /), &
428                  count = (/ mask_size(mid,1), mask_size(mid,2), &
429                             mask_size(mid,3), 1 /) )
430             CALL handle_netcdf_error( 'data_output_mask', 462 )
431
432          ELSE
433!
434!--          If at least part of the mask resides on the PE, send the index
435!--          limits for the target array, otherwise send -9999 to PE0.
436             IF ( mask_size_l(mid,1) > 0 .AND.  mask_size_l(mid,2) > 0 .AND. &
437                  mask_size_l(mid,3) > 0  ) &
438                  THEN
439                ind(1) = mask_start_l(mid,1)
440                ind(2) = mask_start_l(mid,1) + mask_size_l(mid,1) - 1
441                ind(3) = mask_start_l(mid,2)
442                ind(4) = mask_start_l(mid,2) + mask_size_l(mid,2) - 1
443                ind(5) = mask_start_l(mid,3)
444                ind(6) = mask_start_l(mid,3) + mask_size_l(mid,3) - 1
445             ELSE
446                ind(1) = -9999; ind(2) = -9999
447                ind(3) = -9999; ind(4) = -9999
448                ind(5) = -9999; ind(6) = -9999
449             ENDIF
450             CALL MPI_SEND( ind(1), 6, MPI_INTEGER, 0, 0, comm2d, ierr )
451!
452!--          If applicable, send data to PE0.
453             IF ( ind(1) /= -9999 )  THEN
454                CALL MPI_SEND( local_pf(1,1,1), ngp, MPI_REAL, 0, 1, comm2d, &
455                     ierr )
456             ENDIF
457          ENDIF
458!
459!--       A barrier has to be set, because otherwise some PEs may proceed too
460!--       fast so that PE0 may receive wrong data on tag 0.
461          CALL MPI_BARRIER( comm2d, ierr )
462#if defined( __netcdf4 )
463       ENDIF
464#endif
465#else
466!
467!--    (2) For serial execution of PALM, the single processor (PE0) holds all
468!--    data and writes them directly to file.
469       nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),  &
470            id_var_domask(mid,av,if),       &
471            local_pf, &
472            start = (/ 1, 1, 1, domask_time_count(mid,av) /), &
473            count = (/ mask_size_l(mid,1), mask_size_l(mid,2), &
474                       mask_size_l(mid,3), 1 /) )
475       CALL handle_netcdf_error( 'data_output_mask', 463 )
476#endif
477
478       if = if + 1
479
480    ENDDO
481
482!
483!-- Deallocate temporary arrays.
484    DEALLOCATE( local_pf )
485#if defined( __parallel )
486    IF ( myid == 0 )  THEN
487       DEALLOCATE( total_pf )
488    ENDIF
489#endif
490
491
492    CALL cpu_log (log_point(49),'data_output_mask','stop','nobarrier')
493#endif
494
495 END SUBROUTINE data_output_mask
Note: See TracBrowser for help on using the repository browser.