source: palm/trunk/SOURCE/data_output_3d.f90 @ 1433

Last change on this file since 1433 was 1360, checked in by hoffmann, 11 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 20.2 KB
RevLine 
[1]1 SUBROUTINE data_output_3d( 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!
[1310]17! Copyright 1997-2014 Leibniz Universitaet Hannover
[1036]18!--------------------------------------------------------------------------------!
19!
[254]20! Current revisions:
[1106]21! ------------------
[1329]22!
[1360]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: data_output_3d.f90 1360 2014-04-11 17:20:32Z keck $
27!
[1360]28! 1359 2014-04-11 17:15:14Z hoffmann
29! New particle structure integrated.
30!
[1354]31! 1353 2014-04-08 15:21:23Z heinze
32! REAL constants provided with KIND-attribute
33!
[1329]34! 1327 2014-03-21 11:00:16Z raasch
35! parts concerning avs output removed,
36! -netcdf output queries
37!
[1321]38! 1320 2014-03-20 08:40:49Z raasch
[1320]39! ONLY-attribute added to USE-statements,
40! kind-parameters added to all INTEGER and REAL declaration statements,
41! kinds are defined in new module kinds,
42! old module precision_kind is removed,
43! revision history before 2012 removed,
44! comment fields (!:) to be used for variable explanations added to
45! all variable declaration statements
[674]46!
[1319]47! 1318 2014-03-17 13:35:16Z raasch
48! barrier argument removed from cpu_log,
49! module interfaces removed
50!
[1309]51! 1308 2014-03-13 14:58:42Z fricke
52! Check, if the limit of the time dimension is exceeded for parallel output
53! To increase the performance for parallel output, the following is done:
54! - Update of time axis is only done by PE0
55!
[1245]56! 1244 2013-10-31 08:16:56Z raasch
57! Bugfix for index bounds in case of 3d-parallel output
58!
[1116]59! 1115 2013-03-26 18:16:16Z hoffmann
60! ql is calculated by calc_liquid_water_content
61!
[1107]62! 1106 2013-03-04 05:31:38Z raasch
63! array_kind renamed precision_kind
64!
[1077]65! 1076 2012-12-05 08:30:18Z hoffmann
66! Bugfix in output of ql
67!
[1054]68! 1053 2012-11-13 17:11:03Z hoffmann
69! +nr, qr, prr, qc and averaged quantities
70!
[1037]71! 1036 2012-10-22 13:43:42Z raasch
72! code put under GPL (PALM 3.9)
73!
[1035]74! 1031 2012-10-19 14:35:30Z raasch
75! netCDF4 without parallel file support implemented
76!
[1008]77! 1007 2012-09-19 14:30:36Z franke
78! Bugfix: missing calculation of ql_vp added
79!
[1]80! Revision 1.1  1997/09/03 06:29:36  raasch
81! Initial revision
82!
83!
84! Description:
85! ------------
[1031]86! Output of the 3D-arrays in netCDF and/or AVS format.
[1]87!------------------------------------------------------------------------------!
88
[1320]89    USE arrays_3d,                                                             &
90        ONLY:  e, nr, p, pt, q, qc, ql, ql_c, ql_v, qr, rho, sa, tend, u, v,   &
91               vpt, w
92       
[1]93    USE averaging
[1320]94       
95    USE cloud_parameters,                                                      &
96        ONLY:  l_d_cp, prr, pt_d_t
97       
98    USE control_parameters,                                                    &
[1327]99        ONLY:  avs_data_file, cloud_physics, do3d, do3d_avs_n,                 &
100               do3d_no, do3d_time_count, io_blocks, io_group,                  &
101               message_string, netcdf_data_format, ntdim_3d,                   &
[1320]102               nz_do3d, plot_3d_precision, psolver, simulated_time,            &
103               simulated_time_chr, skip_do_avs, time_since_reference_point
104       
105    USE cpulog,                                                                &
106        ONLY:  log_point, cpu_log
107       
108    USE indices,                                                               &
109        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzt,  &
110               nzb
111       
112    USE kinds
113   
[1]114    USE netcdf_control
[1320]115       
116    USE particle_attributes,                                                   &
[1359]117        ONLY:  grid_particles, number_of_particles, particles,                 &
118               particle_advection_start, prt_count
[1320]119       
[1]120    USE pegrid
121
122    IMPLICIT NONE
123
[1320]124    CHARACTER (LEN=9) ::  simulated_time_mod  !:
[1]125
[1320]126    INTEGER(iwp) ::  av        !:
127    INTEGER(iwp) ::  i         !:
128    INTEGER(iwp) ::  if        !:
129    INTEGER(iwp) ::  j         !:
130    INTEGER(iwp) ::  k         !:
131    INTEGER(iwp) ::  n         !:
132    INTEGER(iwp) ::  pos       !:
133    INTEGER(iwp) ::  prec      !:
134    INTEGER(iwp) ::  psi       !:
[1]135
[1320]136    LOGICAL      ::  found     !:
137    LOGICAL      ::  resorted  !:
[1]138
[1320]139    REAL(wp)     ::  mean_r    !:
[1359]140    REAL(wp)     ::  s_r2      !:
[1320]141    REAL(wp)     ::  s_r3      !:
[1]142
[1320]143    REAL(sp), DIMENSION(:,:,:), ALLOCATABLE ::  local_pf  !:
[1]144
[1320]145    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !:
[1]146
147!
148!-- Return, if nothing to output
149    IF ( do3d_no(av) == 0 )  RETURN
[1308]150!
151!-- For parallel netcdf output the time axis must be limited. Return, if this
152!-- limit is exceeded. This could be the case, if the simulated time exceeds
153!-- the given end time by the length of the given output interval.
154    IF ( netcdf_data_format > 4 )  THEN
155       IF ( do3d_time_count(av) + 1 > ntdim_3d(av) )  THEN
156          WRITE ( message_string, * ) 'Output of 3d data is not given at t=',  &
157                                      simulated_time, '&because the maximum ', & 
158                                      'number of output time levels is ',      &
159                                      'exceeded.'
160          CALL message( 'data_output_3d', 'PA0387', 0, 1, 0, 6, 0 )         
161          RETURN
162       ENDIF
163    ENDIF
[1]164
165    CALL cpu_log (log_point(14),'data_output_3d','start')
166
167!
168!-- Open output file.
169!-- Also creates coordinate and fld-file for AVS.
[1031]170!-- For classic or 64bit netCDF output or output of other (old) data formats,
[493]171!-- for a run on more than one PE, each PE opens its own file and
[1]172!-- writes the data of its subdomain in binary format (regardless of the format
173!-- the user has requested). After the run, these files are combined to one
174!-- file by combine_plot_fields in the format requested by the user (netcdf
[493]175!-- and/or avs).
[1031]176!-- For netCDF4/HDF5 output, data is written in parallel into one file.
[1327]177    IF ( netcdf_data_format < 5 )  THEN
178       CALL check_open( 30 )
179       IF ( myid == 0 )  CALL check_open( 106+av*10 )
[493]180    ELSE
[1327]181       CALL check_open( 106+av*10 )
[493]182    ENDIF
[1]183
184!
185!-- Allocate a temporary array with the desired output dimensions.
[667]186    ALLOCATE( local_pf(nxlg:nxrg,nysg:nyng,nzb:nz_do3d) )
[1]187
188!
[1031]189!-- Update the netCDF time axis
[1308]190!-- In case of parallel output, this is only done by PE0 to increase the
191!-- performance.
[1]192#if defined( __netcdf )
[1308]193    do3d_time_count(av) = do3d_time_count(av) + 1
194    IF ( myid == 0 )  THEN
[1327]195       nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_time_3d(av),           &
196                               (/ time_since_reference_point /),            &
197                               start = (/ do3d_time_count(av) /),           &
198                               count = (/ 1 /) )
199       CALL handle_netcdf_error( 'data_output_3d', 376 )
[1]200    ENDIF
201#endif
202
203!
204!-- Loop over all variables to be written.
205    if = 1
206
207    DO  WHILE ( do3d(av,if)(1:1) /= ' ' )
208!
209!--    Store the array chosen on the temporary array.
210       resorted = .FALSE.
211       SELECT CASE ( TRIM( do3d(av,if) ) )
212
213          CASE ( 'e' )
214             IF ( av == 0 )  THEN
215                to_be_resorted => e
216             ELSE
217                to_be_resorted => e_av
218             ENDIF
219
[771]220          CASE ( 'lpt' )
221             IF ( av == 0 )  THEN
222                to_be_resorted => pt
223             ELSE
224                to_be_resorted => lpt_av
225             ENDIF
226
[1053]227          CASE ( 'nr' )
228             IF ( av == 0 )  THEN
229                to_be_resorted => nr
230             ELSE
231                to_be_resorted => nr_av
232             ENDIF
233
[1]234          CASE ( 'p' )
235             IF ( av == 0 )  THEN
[727]236                IF ( psolver /= 'sor' )  CALL exchange_horiz( p, nbgp )
[1]237                to_be_resorted => p
238             ELSE
[727]239                IF ( psolver /= 'sor' )  CALL exchange_horiz( p_av, nbgp )
[1]240                to_be_resorted => p_av
241             ENDIF
242
243          CASE ( 'pc' )  ! particle concentration (requires ghostpoint exchange)
244             IF ( av == 0 )  THEN
[1359]245                IF ( simulated_time >= particle_advection_start )  THEN
246                   tend = prt_count
247                   CALL exchange_horiz( tend, nbgp )
248                ELSE
249                   tend = 0.0_wp
250                ENDIF
[667]251                DO  i = nxlg, nxrg
252                   DO  j = nysg, nyng
[1]253                      DO  k = nzb, nz_do3d
254                         local_pf(i,j,k) = tend(k,j,i)
255                      ENDDO
256                   ENDDO
257                ENDDO
258                resorted = .TRUE.
259             ELSE
[667]260                CALL exchange_horiz( pc_av, nbgp )
[1]261                to_be_resorted => pc_av
262             ENDIF
263
[1359]264          CASE ( 'pr' )  ! mean particle radius (effective radius)
[1]265             IF ( av == 0 )  THEN
[1359]266                IF ( simulated_time >= particle_advection_start )  THEN
267                   DO  i = nxl, nxr
268                      DO  j = nys, nyn
269                         DO  k = nzb, nz_do3d
270                            number_of_particles = prt_count(k,j,i)
271                            IF (number_of_particles <= 0)  CYCLE
272                            particles => grid_particles(k,j,i)%particles(1:number_of_particles)
273                            s_r2 = 0.0_wp
274                            s_r3 = 0.0_wp
275                            DO  n = 1, number_of_particles
276                               IF ( particles(n)%particle_mask )  THEN
277                                  s_r2 = s_r2 + particles(n)%radius**2 * &
278                                         particles(n)%weight_factor
279                                  s_r3 = s_r3 + particles(n)%radius**3 * &
280                                         particles(n)%weight_factor
281                               ENDIF
282                            ENDDO
283                            IF ( s_r2 > 0.0_wp )  THEN
284                               mean_r = s_r3 / s_r2
285                            ELSE
286                               mean_r = 0.0_wp
287                            ENDIF
288                            tend(k,j,i) = mean_r
[1]289                         ENDDO
290                      ENDDO
291                   ENDDO
[1359]292                   CALL exchange_horiz( tend, nbgp )
293                ELSE
294                   tend = 0.0_wp
295                ENDIF
[667]296                DO  i = nxlg, nxrg
297                   DO  j = nysg, nyng
[790]298                      DO  k = nzb, nz_do3d
[1]299                         local_pf(i,j,k) = tend(k,j,i)
300                      ENDDO
301                   ENDDO
302                ENDDO
303                resorted = .TRUE.
304             ELSE
[667]305                CALL exchange_horiz( pr_av, nbgp )
[1]306                to_be_resorted => pr_av
307             ENDIF
308
[1053]309          CASE ( 'prr' )
310             IF ( av == 0 )  THEN
311                CALL exchange_horiz( prr, nbgp )
312                DO  i = nxlg, nxrg
313                   DO  j = nysg, nyng
314                      DO  k = nzb, nzt+1
315                         local_pf(i,j,k) = prr(k,j,i)
316                      ENDDO
317                   ENDDO
318                ENDDO
319             ELSE
320                CALL exchange_horiz( prr_av, nbgp )
321                DO  i = nxlg, nxrg
322                   DO  j = nysg, nyng
323                      DO  k = nzb, nzt+1
324                         local_pf(i,j,k) = prr_av(k,j,i)
325                      ENDDO
326                   ENDDO
327                ENDDO
328             ENDIF
329             resorted = .TRUE.
330
[1]331          CASE ( 'pt' )
332             IF ( av == 0 )  THEN
333                IF ( .NOT. cloud_physics ) THEN
334                   to_be_resorted => pt
335                ELSE
[667]336                   DO  i = nxlg, nxrg
337                      DO  j = nysg, nyng
[1]338                         DO  k = nzb, nz_do3d
[1320]339                            local_pf(i,j,k) = pt(k,j,i) + l_d_cp *             &
340                                                          pt_d_t(k) *          &
[1]341                                                          ql(k,j,i)
342                         ENDDO
343                      ENDDO
344                   ENDDO
345                   resorted = .TRUE.
346                ENDIF
347             ELSE
348                to_be_resorted => pt_av
349             ENDIF
350
351          CASE ( 'q' )
352             IF ( av == 0 )  THEN
353                to_be_resorted => q
354             ELSE
355                to_be_resorted => q_av
356             ENDIF
[691]357
[1053]358          CASE ( 'qc' )
[1]359             IF ( av == 0 )  THEN
[1115]360                to_be_resorted => qc
[1]361             ELSE
[1115]362                to_be_resorted => qc_av
[1]363             ENDIF
364
[1053]365          CASE ( 'ql' )
366             IF ( av == 0 )  THEN
[1115]367                to_be_resorted => ql
[1053]368             ELSE
[1115]369                to_be_resorted => ql_av
[1053]370             ENDIF
371
[1]372          CASE ( 'ql_c' )
373             IF ( av == 0 )  THEN
374                to_be_resorted => ql_c
375             ELSE
376                to_be_resorted => ql_c_av
377             ENDIF
378
379          CASE ( 'ql_v' )
380             IF ( av == 0 )  THEN
381                to_be_resorted => ql_v
382             ELSE
383                to_be_resorted => ql_v_av
384             ENDIF
385
386          CASE ( 'ql_vp' )
387             IF ( av == 0 )  THEN
[1359]388                IF ( simulated_time >= particle_advection_start )  THEN
389                   DO  i = nxl, nxr
390                      DO  j = nys, nyn
391                         DO  k = nzb, nz_do3d
392                            number_of_particles = prt_count(k,j,i)
393                            IF (number_of_particles <= 0)  CYCLE
394                            particles => grid_particles(k,j,i)%particles(1:number_of_particles)
395                            DO  n = 1, number_of_particles
396                               IF ( particles(n)%particle_mask )  THEN
397                                  tend(k,j,i) =  tend(k,j,i) +                 &
398                                                 particles(n)%weight_factor /  &
399                                                 prt_count(k,j,i)
400                               ENDIF
401                            ENDDO
[1007]402                         ENDDO
403                      ENDDO
404                   ENDDO
[1359]405                   CALL exchange_horiz( tend, nbgp )
406                ELSE
407                   tend = 0.0_wp
408                ENDIF
[1007]409                DO  i = nxlg, nxrg
410                   DO  j = nysg, nyng
411                      DO  k = nzb, nz_do3d
412                         local_pf(i,j,k) = tend(k,j,i)
413                      ENDDO
414                   ENDDO
415                ENDDO
416                resorted = .TRUE.
[1]417             ELSE
[1007]418                CALL exchange_horiz( ql_vp_av, nbgp )
[1]419                to_be_resorted => ql_vp_av
420             ENDIF
421
[1053]422          CASE ( 'qr' )
423             IF ( av == 0 )  THEN
424                to_be_resorted => qr
425             ELSE
426                to_be_resorted => qr_av
427             ENDIF
428
[1]429          CASE ( 'qv' )
430             IF ( av == 0 )  THEN
[667]431                DO  i = nxlg, nxrg
432                   DO  j = nysg, nyng
[1]433                      DO  k = nzb, nz_do3d
434                         local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
435                      ENDDO
436                   ENDDO
437                ENDDO
438                resorted = .TRUE.
439             ELSE
440                to_be_resorted => qv_av
441             ENDIF
442
[96]443          CASE ( 'rho' )
444             IF ( av == 0 )  THEN
445                to_be_resorted => rho
446             ELSE
447                to_be_resorted => rho_av
448             ENDIF
[691]449
[1]450          CASE ( 's' )
451             IF ( av == 0 )  THEN
452                to_be_resorted => q
453             ELSE
[355]454                to_be_resorted => s_av
[1]455             ENDIF
[691]456
[96]457          CASE ( 'sa' )
458             IF ( av == 0 )  THEN
459                to_be_resorted => sa
460             ELSE
461                to_be_resorted => sa_av
462             ENDIF
[691]463
[1]464          CASE ( 'u' )
465             IF ( av == 0 )  THEN
466                to_be_resorted => u
467             ELSE
468                to_be_resorted => u_av
469             ENDIF
470
471          CASE ( 'v' )
472             IF ( av == 0 )  THEN
473                to_be_resorted => v
474             ELSE
475                to_be_resorted => v_av
476             ENDIF
477
478          CASE ( 'vpt' )
479             IF ( av == 0 )  THEN
480                to_be_resorted => vpt
481             ELSE
482                to_be_resorted => vpt_av
483             ENDIF
484
485          CASE ( 'w' )
486             IF ( av == 0 )  THEN
487                to_be_resorted => w
488             ELSE
489                to_be_resorted => w_av
490             ENDIF
491
492          CASE DEFAULT
493!
494!--          User defined quantity
[1320]495             CALL user_data_output_3d( av, do3d(av,if), found, local_pf,       &
[1]496                                       nz_do3d )
497             resorted = .TRUE.
498
[254]499             IF ( .NOT. found )  THEN
[1320]500                message_string =  'no output available for: ' //               &
[274]501                                  TRIM( do3d(av,if) )
[254]502                CALL message( 'data_output_3d', 'PA0182', 0, 0, 0, 6, 0 )
[1]503             ENDIF
504
505       END SELECT
506
507!
508!--    Resort the array to be output, if not done above
509       IF ( .NOT. resorted )  THEN
[667]510          DO  i = nxlg, nxrg
511             DO  j = nysg, nyng
[1]512                DO  k = nzb, nz_do3d
513                   local_pf(i,j,k) = to_be_resorted(k,j,i)
514                ENDDO
515             ENDDO
516          ENDDO
517       ENDIF
518
519!
[1327]520!--    Output of the 3D-array
521#if defined( __parallel )
522       IF ( netcdf_data_format < 5 )  THEN
[1]523!
[1327]524!--       Non-parallel netCDF output. Data is output in parallel in
525!--       FORTRAN binary format here, and later collected into one file by
526!--       combine_plot_fields
527          IF ( myid == 0 )  THEN
528             WRITE ( 30 )  time_since_reference_point,                   &
529                           do3d_time_count(av), av
[1]530          ENDIF
[1327]531          DO  i = 0, io_blocks-1
532             IF ( i == io_group )  THEN
533                WRITE ( 30 )  nxlg, nxrg, nysg, nyng, nzb, nz_do3d
534                WRITE ( 30 )  local_pf
535             ENDIF
[1]536#if defined( __parallel )
[1327]537             CALL MPI_BARRIER( comm2d, ierr )
[759]538#endif
[1327]539          ENDDO
[559]540
[1327]541       ELSE
[646]542#if defined( __netcdf )
[493]543!
[1327]544!--       Parallel output in netCDF4/HDF5 format.
545!--       Do not output redundant ghost point data except for the
546!--       boundaries of the total domain.
547          IF ( nxr == nx  .AND.  nyn /= ny )  THEN
548             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
549                               local_pf(nxl:nxr+1,nys:nyn,nzb:nz_do3d),  &
550                start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /),  &
551                count = (/ nxr-nxl+2, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
552          ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
553             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
554                               local_pf(nxl:nxr,nys:nyn+1,nzb:nz_do3d),  &
555                start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /),  &
556                count = (/ nxr-nxl+1, nyn-nys+2, nz_do3d-nzb+1, 1 /) )
557          ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
558             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
559                             local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d),  &
560                start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /),  &
561                count = (/ nxr-nxl+2, nyn-nys+2, nz_do3d-nzb+1, 1 /) )
562          ELSE
563             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
564                                 local_pf(nxl:nxr,nys:nyn,nzb:nz_do3d),  &
565                start = (/ nxl+1, nys+1, nzb+1, do3d_time_count(av) /),  &
566                count = (/ nxr-nxl+1, nyn-nys+1, nz_do3d-nzb+1, 1 /) )
567          ENDIF
568          CALL handle_netcdf_error( 'data_output_3d', 386 )
[646]569#endif
[1327]570       ENDIF
[1]571#else
572#if defined( __netcdf )
[1327]573       nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),        &
574                         local_pf(nxl:nxr+1,nys:nyn+1,nzb:nz_do3d),      &
575                         start = (/ 1, 1, 1, do3d_time_count(av) /),     &
576                         count = (/ nx+2, ny+2, nz_do3d-nzb+1, 1 /) )
577       CALL handle_netcdf_error( 'data_output_3d', 446 )
[1]578#endif
579#endif
580
581       if = if + 1
582
583    ENDDO
584
585!
586!-- Deallocate temporary array.
587    DEALLOCATE( local_pf )
588
589
[1318]590    CALL cpu_log( log_point(14), 'data_output_3d', 'stop' )
[1]591
592!
593!-- Formats.
[1320]5943300 FORMAT ('variable ',I4,'  file=',A,'  filetype=unformatted  skip=',I12/   &
[1]595             'label = ',A,A)
596
597 END SUBROUTINE data_output_3d
Note: See TracBrowser for help on using the repository browser.