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

Last change on this file since 1035 was 1035, checked in by raasch, 11 years ago

revisions r1031 and r1034 documented

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