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

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

last commit documented

  • Property svn:keywords set to Id
File size: 16.6 KB
RevLine 
[298]1 SUBROUTINE data_output_mask( av )
2
3!------------------------------------------------------------------------------!
[484]4! Current revisions:
[298]5! -----------------
6!
[1008]7!
[298]8! Former revisions:
9! -----------------
10! $Id: data_output_mask.f90 1008 2012-09-19 14:49:14Z letzel $
11!
[1008]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!
[772]16! 771 2011-10-27 10:56:21Z heinze
17! +lpt
18!
[668]19! 667 2010-12-23 12:06:00Z suehring/gryschka
20! Calls of exchange_horiz are modified.
21!
[565]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!
[494]26! 493 2010-03-01 08:30:24Z raasch
27! netcdf_format_mask* and format_parallel_io replaced by netcdf_data_format
28!
[482]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!
[449]32! 410 2009-12-04 17:05:40Z letzel
33! Initial version
[298]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
[493]72!
[298]73!-- Open output file.
[493]74    IF ( netcdf_output  .AND.  ( myid == 0  .OR.  netcdf_data_format > 2 ) ) &
75    THEN
[564]76       CALL check_open( 200+mid+av*max_masks )
[409]77    ENDIF 
[298]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
[493]92    IF ( netcdf_output  .AND.  ( myid == 0  .OR.  netcdf_data_format > 2 ) ) &
93    THEN
[298]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 /) )
[564]98       CALL handle_netcdf_error( 'data_output_mask', 460 )
[298]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
[493]108       IF ( netcdf_data_format < 3   .AND.  myid == 0  .AND.  if > 1 )  THEN
[298]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
[771]125          CASE ( 'lpt' )
126             IF ( av == 0 )  THEN
127                to_be_resorted => pt
128             ELSE
129                to_be_resorted => lpt_av
130             ENDIF
131
[298]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
[667]142                CALL exchange_horiz( tend, nbgp )
[298]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
[667]153                CALL exchange_horiz( pc_av, nbgp )
[298]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
[1007]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
[298]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
[667]180                CALL exchange_horiz( tend, nbgp )
[298]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
[667]191                CALL exchange_horiz( pr_av, nbgp )
[298]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
[1007]222
[298]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
[1007]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.
[298]268             ELSE
[1007]269                CALL exchange_horiz( ql_vp_av, nbgp )
[298]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
[1007]295
[298]296          CASE ( 's' )
297             IF ( av == 0 )  THEN
298                to_be_resorted => q
299             ELSE
[356]300                to_be_resorted => s_av
[298]301             ENDIF
[1007]302
[298]303          CASE ( 'sa' )
304             IF ( av == 0 )  THEN
305                to_be_resorted => sa
306             ELSE
307                to_be_resorted => sa_av
308             ENDIF
[1007]309
[298]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) )
[564]347                CALL message( 'data_output_mask', 'PA0327', 0, 0, 0, 6, 0 )
[298]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
[493]373!--    netcdf_data_format.
[298]374#if defined( __parallel )
375#if defined( __netcdf4 )
[493]376       IF ( netcdf_data_format > 2 )  THEN
[298]377!
378!--       (1) a. Parallel I/O using NetCDF 4 (not yet tested)
379          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),  &
[340]380               id_var_domask(mid,av,if),  &
[409]381               local_pf,  &
[340]382               start = (/ mask_start_l(mid,1), mask_start_l(mid,2),  &
[409]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 /) )
[564]386          CALL handle_netcdf_error( 'data_output_mask', 461 )
[298]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 /) )
[564]433             CALL handle_netcdf_error( 'data_output_mask', 462 )
[298]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),       &
[475]474            local_pf, &
[298]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 /) )
[564]478       CALL handle_netcdf_error( 'data_output_mask', 463 )
[298]479#endif
480
481       if = if + 1
[667]482
[298]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.