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

Last change on this file since 1272 was 1093, 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! -----------------
[1093]22!
[298]23!
24! Former revisions:
25! -----------------
26! $Id: data_output_mask.f90 1093 2013-02-02 12:58:49Z witha $
27!
[1093]28! 1092 2013-02-02 11:24:22Z raasch
29! unused variables removed
30!
[1037]31! 1036 2012-10-22 13:43:42Z raasch
32! code put under GPL (PALM 3.9)
33!
[1035]34! 1031 2012-10-19 14:35:30Z raasch
35! netCDF4 without parallel file support implemented
36!
[1008]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!
[772]41! 771 2011-10-27 10:56:21Z heinze
42! +lpt
43!
[668]44! 667 2010-12-23 12:06:00Z suehring/gryschka
45! Calls of exchange_horiz are modified.
46!
[565]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!
[494]51! 493 2010-03-01 08:30:24Z raasch
52! netcdf_format_mask* and format_parallel_io replaced by netcdf_data_format
53!
[482]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!
[449]57! 410 2009-12-04 17:05:40Z letzel
58! Initial version
[298]59!
60! Description:
61! ------------
[1031]62! Masked data output in netCDF format for current mask (current value of mid).
[298]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
[1092]81    INTEGER ::  av, ngp, i, if, j, k, n, psi, sender, &
[298]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
[493]97!
[298]98!-- Open output file.
[1031]99    IF ( netcdf_output  .AND.  ( myid == 0  .OR.  netcdf_data_format > 4 ) ) &
[493]100    THEN
[564]101       CALL check_open( 200+mid+av*max_masks )
[409]102    ENDIF 
[298]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!
[1031]115!-- Update the netCDF time axis.
[298]116    domask_time_count(mid,av) = domask_time_count(mid,av) + 1
[1031]117    IF ( netcdf_output  .AND.  ( myid == 0  .OR.  netcdf_data_format > 4 ) ) &
[493]118    THEN
[298]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 /) )
[564]123       CALL handle_netcdf_error( 'data_output_mask', 460 )
[298]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
[1031]133       IF ( netcdf_data_format < 5   .AND.  myid == 0  .AND.  if > 1 )  THEN
[298]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
[771]150          CASE ( 'lpt' )
151             IF ( av == 0 )  THEN
152                to_be_resorted => pt
153             ELSE
154                to_be_resorted => lpt_av
155             ENDIF
156
[298]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
[667]167                CALL exchange_horiz( tend, nbgp )
[298]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
[667]178                CALL exchange_horiz( pc_av, nbgp )
[298]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
[1007]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
[298]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
[667]205                CALL exchange_horiz( tend, nbgp )
[298]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
[667]216                CALL exchange_horiz( pr_av, nbgp )
[298]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
[1007]247
[298]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
[1007]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.
[298]293             ELSE
[1007]294                CALL exchange_horiz( ql_vp_av, nbgp )
[298]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
[1007]320
[298]321          CASE ( 's' )
322             IF ( av == 0 )  THEN
323                to_be_resorted => q
324             ELSE
[356]325                to_be_resorted => s_av
[298]326             ENDIF
[1007]327
[298]328          CASE ( 'sa' )
329             IF ( av == 0 )  THEN
330                to_be_resorted => sa
331             ELSE
332                to_be_resorted => sa_av
333             ENDIF
[1007]334
[298]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) )
[564]372                CALL message( 'data_output_mask', 'PA0327', 0, 0, 0, 6, 0 )
[298]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
[1031]393!--     a. with netCDF 4 parallel I/O-enabled library
394!--     b. with netCDF 3 library
[298]395!--    (2) for serial execution.
396!--    The choice of method depends on the correct setting of preprocessor
[1031]397!--    directives __parallel and __netcdf4_parallel as well as on the parameter
[493]398!--    netcdf_data_format.
[298]399#if defined( __parallel )
[1031]400#if defined( __netcdf4_parallel )
401       IF ( netcdf_data_format > 4 )  THEN
[298]402!
[1031]403!--       (1) a. Parallel I/O using netCDF 4 (not yet tested)
[298]404          nc_stat = NF90_PUT_VAR( id_set_mask(mid,av),  &
[340]405               id_var_domask(mid,av,if),  &
[409]406               local_pf,  &
[340]407               start = (/ mask_start_l(mid,1), mask_start_l(mid,2),  &
[409]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 /) )
[564]411          CALL handle_netcdf_error( 'data_output_mask', 461 )
[298]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 /) )
[564]458             CALL handle_netcdf_error( 'data_output_mask', 462 )
[298]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 )
[1031]490#if defined( __netcdf4_parallel )
[298]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),       &
[475]499            local_pf, &
[298]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 /) )
[564]503       CALL handle_netcdf_error( 'data_output_mask', 463 )
[298]504#endif
505
506       if = if + 1
[667]507
[298]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.