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

Last change on this file since 1353 was 1353, checked in by heinze, 7 years ago

REAL constants provided with KIND-attribute

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