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

Last change on this file since 1524 was 1439, checked in by heinze, 10 years ago

last commit documented

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