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

Last change on this file since 876 was 772, checked in by heinze, 13 years ago

last commit documented

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