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

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