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

Last change on this file since 771 was 771, checked in by heinze, 12 years ago

Output of liquid water potential temperature in case of cloud_physics=.T. enabled

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