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

Last change on this file since 1327 was 1327, checked in by raasch, 8 years ago

Changed:


-s real64 removed (.mrun.config.hlrnIII)
-r8 removed (.mrun.config.imuk)
deleted: .mrun.config.imuk_ice2_netcdf4 .mrun.config.imuk_hlrn

REAL constants defined as wp-kind in modules

"baroclinicity" renamed "baroclinity", "ocean version" replaced by
"ocean mode"

code parts concerning old output formats "iso2d" and "avs" removed.
netCDF is the only remaining output format.

Errors:


bugfix: duplicate error message 56 removed

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