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

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

last commit documented

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