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

Last change on this file since 484 was 484, checked in by raasch, 14 years ago

typo in file headers removed

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