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

Last change on this file since 1319 was 1319, checked in by raasch, 10 years ago

last commit documented

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