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

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