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

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

last commit documented

  • Property svn:keywords set to Id
File size: 17.6 KB
RevLine 
[298]1 SUBROUTINE data_output_mask( av )
2
[1036]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!
[484]20! Current revisions:
[298]21! -----------------
22!
[1035]23!
[298]24! Former revisions:
25! -----------------
26! $Id: data_output_mask.f90 1037 2012-10-22 14:10:22Z maronga $
27!
[1037]28! 1036 2012-10-22 13:43:42Z raasch
29! code put under GPL (PALM 3.9)
30!
[1035]31! 1031 2012-10-19 14:35:30Z raasch
32! netCDF4 without parallel file support implemented
33!
[1008]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!
[772]38! 771 2011-10-27 10:56:21Z heinze
39! +lpt
40!
[668]41! 667 2010-12-23 12:06:00Z suehring/gryschka
42! Calls of exchange_horiz are modified.
43!
[565]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!
[494]48! 493 2010-03-01 08:30:24Z raasch
49! netcdf_format_mask* and format_parallel_io replaced by netcdf_data_format
50!
[482]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!
[449]54! 410 2009-12-04 17:05:40Z letzel
55! Initial version
[298]56!
57! Description:
58! ------------
[1031]59! Masked data output in netCDF format for current mask (current value of mid).
[298]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
[493]94!
[298]95!-- Open output file.
[1031]96    IF ( netcdf_output  .AND.  ( myid == 0  .OR.  netcdf_data_format > 4 ) ) &
[493]97    THEN
[564]98       CALL check_open( 200+mid+av*max_masks )
[409]99    ENDIF 
[298]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!
[1031]112!-- Update the netCDF time axis.
[298]113    domask_time_count(mid,av) = domask_time_count(mid,av) + 1
[1031]114    IF ( netcdf_output  .AND.  ( myid == 0  .OR.  netcdf_data_format > 4 ) ) &
[493]115    THEN
[298]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 /) )
[564]120       CALL handle_netcdf_error( 'data_output_mask', 460 )
[298]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
[1031]130       IF ( netcdf_data_format < 5   .AND.  myid == 0  .AND.  if > 1 )  THEN
[298]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
[771]147          CASE ( 'lpt' )
148             IF ( av == 0 )  THEN
149                to_be_resorted => pt
150             ELSE
151                to_be_resorted => lpt_av
152             ENDIF
153
[298]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
[667]164                CALL exchange_horiz( tend, nbgp )
[298]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
[667]175                CALL exchange_horiz( pc_av, nbgp )
[298]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
[1007]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
[298]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
[667]202                CALL exchange_horiz( tend, nbgp )
[298]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
[667]213                CALL exchange_horiz( pr_av, nbgp )
[298]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
[1007]244
[298]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
[1007]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.
[298]290             ELSE
[1007]291                CALL exchange_horiz( ql_vp_av, nbgp )
[298]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
[1007]317
[298]318          CASE ( 's' )
319             IF ( av == 0 )  THEN
320                to_be_resorted => q
321             ELSE
[356]322                to_be_resorted => s_av
[298]323             ENDIF
[1007]324
[298]325          CASE ( 'sa' )
326             IF ( av == 0 )  THEN
327                to_be_resorted => sa
328             ELSE
329                to_be_resorted => sa_av
330             ENDIF
[1007]331
[298]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) )
[564]369                CALL message( 'data_output_mask', 'PA0327', 0, 0, 0, 6, 0 )
[298]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
[1031]390!--     a. with netCDF 4 parallel I/O-enabled library
391!--     b. with netCDF 3 library
[298]392!--    (2) for serial execution.
393!--    The choice of method depends on the correct setting of preprocessor
[1031]394!--    directives __parallel and __netcdf4_parallel as well as on the parameter
[493]395!--    netcdf_data_format.
[298]396#if defined( __parallel )
[1031]397#if defined( __netcdf4_parallel )
398       IF ( netcdf_data_format > 4 )  THEN
[298]399!
[1031]400!--       (1) a. Parallel I/O using netCDF 4 (not yet tested)
[298]401          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),  &
[340]402               id_var_domask(mid,av,if),  &
[409]403               local_pf,  &
[340]404               start = (/ mask_start_l(mid,1), mask_start_l(mid,2),  &
[409]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 /) )
[564]408          CALL handle_netcdf_error( 'data_output_mask', 461 )
[298]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 /) )
[564]455             CALL handle_netcdf_error( 'data_output_mask', 462 )
[298]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 )
[1031]487#if defined( __netcdf4_parallel )
[298]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),       &
[475]496            local_pf, &
[298]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 /) )
[564]500       CALL handle_netcdf_error( 'data_output_mask', 463 )
[298]501#endif
502
503       if = if + 1
[667]504
[298]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.