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

Last change on this file since 1070 was 1037, checked in by raasch, 12 years ago

last commit documented

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