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

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

code has been put under the GNU General Public License (v3)

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