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

Last change on this file since 2756 was 2756, checked in by suehring, 6 years ago

Fill values for 3D data output of chemical species introduced

  • Property svn:keywords set to Id
File size: 27.1 KB
Line 
1!> @file data_output_3d.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! 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-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! ------------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: data_output_3d.f90 2756 2018-01-16 18:11:14Z suehring $
27! Fill values for 3D output of chemical species introduced.
28!
29! 2746 2018-01-15 12:06:04Z suehring
30! Move flag plant canopy to modules
31!
32! 2718 2018-01-02 08:49:38Z maronga
33! Corrected "Former revisions" section
34!
35! 2696 2017-12-14 17:12:51Z kanani
36! Change in file header (GPL part)
37! Implementation of turbulence_closure_mod (TG)
38! Implementation of chemistry module (FK)
39! Set fill values at topography grid points or e.g. non-natural-type surface
40! in case of LSM output (MS)
41!
42! 2512 2017-10-04 08:26:59Z raasch
43! upper bounds of 3d output changed from nx+1,ny+1 to nx,ny
44! no output of ghost layer data
45!
46! 2292 2017-06-20 09:51:42Z schwenkel
47! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
48! includes two more prognostic equations for cloud drop concentration (nc) 
49! and cloud water content (qc).
50!
51! 2233 2017-05-30 18:08:54Z suehring
52!
53! 2232 2017-05-30 17:47:52Z suehring
54! Adjustments to new topography concept
55!
56! 2209 2017-04-19 09:34:46Z kanani
57! Added plant canopy model output
58!
59! 2031 2016-10-21 15:11:58Z knoop
60! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
61!
62! 2011 2016-09-19 17:29:57Z kanani
63! Flag urban_surface is now defined in module control_parameters,
64! changed prefix for urban surface model output to "usm_",
65! introduced control parameter varnamelength for LEN of trimvar.
66!
67! 2007 2016-08-24 15:47:17Z kanani
68! Added support for new urban surface model (temporary modifications of
69! SELECT CASE ( ) necessary, see variable trimvar)
70!
71! 2000 2016-08-20 18:09:15Z knoop
72! Forced header and separation lines into 80 columns
73!
74! 1980 2016-07-29 15:51:57Z suehring
75! Bugfix, in order to steer user-defined output, setting flag found explicitly
76! to .F.
77!
78! 1976 2016-07-27 13:28:04Z maronga
79! Output of radiation quantities is now done directly in the respective module
80!
81! 1972 2016-07-26 07:52:02Z maronga
82! Output of land surface quantities is now done directly in the respective module.
83! Unnecessary directive __parallel removed.
84!
85! 1960 2016-07-12 16:34:24Z suehring
86! Scalar surface flux added
87!
88! 1849 2016-04-08 11:33:18Z hoffmann
89! prr moved to arrays_3d
90!
91! 1822 2016-04-07 07:49:42Z hoffmann
92! prr vertical dimensions set to nzb_do to nzt_do. Unused variables deleted.
93!
94! 1808 2016-04-05 19:44:00Z raasch
95! test output removed
96!
97! 1783 2016-03-06 18:36:17Z raasch
98! name change of netcdf routines and module + related changes
99!
100! 1745 2016-02-05 13:06:51Z gronemeier
101! Bugfix: test if time axis limit exceeds moved to point after call of check_open
102!
103! 1691 2015-10-26 16:17:44Z maronga
104! Added output of radiative heating rates for RRTMG
105!
106! 1682 2015-10-07 23:56:08Z knoop
107! Code annotations made doxygen readable
108!
109! 1585 2015-04-30 07:05:52Z maronga
110! Added support for RRTMG
111!
112! 1551 2015-03-03 14:18:16Z maronga
113! Added suppport for land surface model and radiation model output. In the course
114! of this action, the limits for vertical loops have been changed (from nzb and
115! nzt+1 to nzb_do and nzt_do, respectively in order to allow soil model output).
116! Moreover, a new vertical grid zs was introduced.
117!
118! 1359 2014-04-11 17:15:14Z hoffmann
119! New particle structure integrated.
120!
121! 1353 2014-04-08 15:21:23Z heinze
122! REAL constants provided with KIND-attribute
123!
124! 1327 2014-03-21 11:00:16Z raasch
125! parts concerning avs output removed,
126! -netcdf output queries
127!
128! 1320 2014-03-20 08:40:49Z raasch
129! ONLY-attribute added to USE-statements,
130! kind-parameters added to all INTEGER and REAL declaration statements,
131! kinds are defined in new module kinds,
132! old module precision_kind is removed,
133! revision history before 2012 removed,
134! comment fields (!:) to be used for variable explanations added to
135! all variable declaration statements
136!
137! 1318 2014-03-17 13:35:16Z raasch
138! barrier argument removed from cpu_log,
139! module interfaces removed
140!
141! 1308 2014-03-13 14:58:42Z fricke
142! Check, if the limit of the time dimension is exceeded for parallel output
143! To increase the performance for parallel output, the following is done:
144! - Update of time axis is only done by PE0
145!
146! 1244 2013-10-31 08:16:56Z raasch
147! Bugfix for index bounds in case of 3d-parallel output
148!
149! 1115 2013-03-26 18:16:16Z hoffmann
150! ql is calculated by calc_liquid_water_content
151!
152! 1106 2013-03-04 05:31:38Z raasch
153! array_kind renamed precision_kind
154!
155! 1076 2012-12-05 08:30:18Z hoffmann
156! Bugfix in output of ql
157!
158! 1053 2012-11-13 17:11:03Z hoffmann
159! +nr, qr, prr, qc and averaged quantities
160!
161! 1036 2012-10-22 13:43:42Z raasch
162! code put under GPL (PALM 3.9)
163!
164! 1031 2012-10-19 14:35:30Z raasch
165! netCDF4 without parallel file support implemented
166!
167! 1007 2012-09-19 14:30:36Z franke
168! Bugfix: missing calculation of ql_vp added
169!
170! Revision 1.1  1997/09/03 06:29:36  raasch
171! Initial revision
172!
173!
174! Description:
175! ------------
176!> Output of the 3D-arrays in netCDF and/or AVS format.
177!------------------------------------------------------------------------------!
178 SUBROUTINE data_output_3d( av )
179 
180
181    USE arrays_3d,                                                             &
182        ONLY:  e, nc, nr, p, pt, prr, q, qc, ql, ql_c, ql_v, qr, rho_ocean, s, &
183               sa, tend, u, v, vpt, w
184       
185    USE averaging
186       
187#if defined( __chem )
188    USE chemistry_model_mod,                                                   &
189        ONLY:  chem_data_output_3d
190#endif       
191
192    USE cloud_parameters,                                                      &
193        ONLY:  l_d_cp, pt_d_t
194       
195    USE control_parameters,                                                    &
196        ONLY:  air_chemistry, cloud_physics, do3d, do3d_no, do3d_time_count,   &
197               io_blocks, io_group, land_surface, message_string,              &
198               ntdim_3d, nz_do3d,  plant_canopy,                               &
199               psolver, simulated_time, time_since_reference_point,            &
200               urban_surface, varnamelength
201       
202    USE cpulog,                                                                &
203        ONLY:  log_point, cpu_log
204       
205    USE indices,                                                               &
206        ONLY:  nbgp, nx, nxl, nxr, ny, nyn, nys, nzb, wall_flags_0
207       
208    USE kinds
209   
210    USE land_surface_model_mod,                                                &
211        ONLY: lsm_data_output_3d, nzb_soil, nzt_soil
212
213#if defined( __netcdf )
214    USE NETCDF
215#endif
216
217    USE netcdf_interface,                                                      &
218        ONLY:  fill_value, id_set_3d, id_var_do3d, id_var_time_3d, nc_stat,    &
219               netcdf_data_format, netcdf_handle_error
220       
221    USE particle_attributes,                                                   &
222        ONLY:  grid_particles, number_of_particles, particles,                 &
223               particle_advection_start, prt_count
224       
225    USE pegrid
226
227    USE plant_canopy_model_mod,                                                &
228        ONLY:  pcm_data_output_3d
229       
230    USE radiation_model_mod,                                                   &
231        ONLY:  nzub, nzut, radiation, radiation_data_output_3d
232
233    USE turbulence_closure_mod,                                                &
234        ONLY:  tcm_data_output_3d
235
236    USE urban_surface_mod,                                                     &
237        ONLY:  usm_data_output_3d
238
239
240    IMPLICIT NONE
241
242    INTEGER(iwp) ::  av        !<
243    INTEGER(iwp) ::  flag_nr   !< number of masking flag
244    INTEGER(iwp) ::  i         !<
245    INTEGER(iwp) ::  if        !<
246    INTEGER(iwp) ::  j         !<
247    INTEGER(iwp) ::  k         !<
248    INTEGER(iwp) ::  n         !<
249    INTEGER(iwp) ::  nzb_do    !< vertical lower limit for data output
250    INTEGER(iwp) ::  nzt_do    !< vertical upper limit for data output
251
252    LOGICAL      ::  found     !<
253    LOGICAL      ::  resorted  !<
254
255    REAL(wp)     ::  mean_r    !<
256    REAL(wp)     ::  s_r2      !<
257    REAL(wp)     ::  s_r3      !<
258
259    REAL(sp), DIMENSION(:,:,:), ALLOCATABLE ::  local_pf  !<
260
261    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !<
262
263    CHARACTER (LEN=varnamelength) ::  trimvar  !< TRIM of output-variable string
264
265!
266!-- Return, if nothing to output
267    IF ( do3d_no(av) == 0 )  RETURN
268
269    CALL cpu_log (log_point(14),'data_output_3d','start')
270
271!
272!-- Open output file.
273!-- For classic or 64bit netCDF output on more than one PE, each PE opens its
274!-- own file and writes the data of its subdomain in binary format. After the
275!-- run, these files are combined to one NetCDF file by combine_plot_fields.
276!-- For netCDF4/HDF5 output, data is written in parallel into one file.
277    IF ( netcdf_data_format < 5 )  THEN
278       CALL check_open( 30 )
279       IF ( myid == 0 )  CALL check_open( 106+av*10 )
280    ELSE
281       CALL check_open( 106+av*10 )
282    ENDIF
283
284!
285!-- For parallel netcdf output the time axis must be limited. Return, if this
286!-- limit is exceeded. This could be the case, if the simulated time exceeds
287!-- the given end time by the length of the given output interval.
288    IF ( netcdf_data_format > 4 )  THEN
289       IF ( do3d_time_count(av) + 1 > ntdim_3d(av) )  THEN
290          WRITE ( message_string, * ) 'Output of 3d data is not given at t=',  &
291                                      simulated_time, '&because the maximum ', & 
292                                      'number of output time levels is ',      &
293                                      'exceeded.'
294          CALL message( 'data_output_3d', 'PA0387', 0, 1, 0, 6, 0 )
295          CALL cpu_log( log_point(14), 'data_output_3d', 'stop' )
296          RETURN
297       ENDIF
298    ENDIF
299
300!
301!-- Update the netCDF time axis
302!-- In case of parallel output, this is only done by PE0 to increase the
303!-- performance.
304#if defined( __netcdf )
305    do3d_time_count(av) = do3d_time_count(av) + 1
306    IF ( myid == 0 )  THEN
307       nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_time_3d(av),           &
308                               (/ time_since_reference_point /),            &
309                               start = (/ do3d_time_count(av) /),           &
310                               count = (/ 1 /) )
311       CALL netcdf_handle_error( 'data_output_3d', 376 )
312    ENDIF
313#endif
314
315!
316!-- Loop over all variables to be written.
317    if = 1
318
319    DO  WHILE ( do3d(av,if)(1:1) /= ' ' )
320
321!
322!--    Temporary solution to account for data output within the new urban
323!--    surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar ).
324!--    Store the array chosen on the temporary array.
325       trimvar = TRIM( do3d(av,if) )
326       IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
327          trimvar = 'usm_output'
328          resorted = .TRUE.
329          nzb_do   = nzub
330          nzt_do   = nzut
331       ELSE
332          resorted = .FALSE.
333          nzb_do   = nzb
334          nzt_do   = nz_do3d
335       ENDIF
336!
337!--    Set flag to steer output of radiation, land-surface, or user-defined
338!--    quantities
339       found = .FALSE.
340!
341!--    Allocate a temporary array with the desired output dimensions.
342       ALLOCATE( local_pf(nxl:nxr,nys:nyn,nzb_do:nzt_do) )
343!
344!--    Before each output, set array local_pf to fill value
345       local_pf = fill_value
346!
347!--    Set masking flag for topography for not resorted arrays
348       flag_nr = 0
349
350       SELECT CASE ( trimvar )
351
352          CASE ( 'e' )
353             IF ( av == 0 )  THEN
354                to_be_resorted => e
355             ELSE
356                to_be_resorted => e_av
357             ENDIF
358
359          CASE ( 'lpt' )
360             IF ( av == 0 )  THEN
361                to_be_resorted => pt
362             ELSE
363                to_be_resorted => lpt_av
364             ENDIF
365
366          CASE ( 'nc' )
367             IF ( av == 0 )  THEN
368                to_be_resorted => nc
369             ELSE
370                to_be_resorted => nc_av
371             ENDIF
372
373          CASE ( 'nr' )
374             IF ( av == 0 )  THEN
375                to_be_resorted => nr
376             ELSE
377                to_be_resorted => nr_av
378             ENDIF
379
380          CASE ( 'p' )
381             IF ( av == 0 )  THEN
382                IF ( psolver /= 'sor' )  CALL exchange_horiz( p, nbgp )
383                to_be_resorted => p
384             ELSE
385                IF ( psolver /= 'sor' )  CALL exchange_horiz( p_av, nbgp )
386                to_be_resorted => p_av
387             ENDIF
388
389          CASE ( 'pc' )  ! particle concentration (requires ghostpoint exchange)
390             IF ( av == 0 )  THEN
391                IF ( simulated_time >= particle_advection_start )  THEN
392                   tend = prt_count
393                ELSE
394                   tend = 0.0_wp
395                ENDIF
396                DO  i = nxl, nxr
397                   DO  j = nys, nyn
398                      DO  k = nzb_do, nzt_do
399                         local_pf(i,j,k) = tend(k,j,i)
400                      ENDDO
401                   ENDDO
402                ENDDO
403                resorted = .TRUE.
404             ELSE
405                to_be_resorted => pc_av
406             ENDIF
407
408          CASE ( 'pr' )  ! mean particle radius (effective radius)
409             IF ( av == 0 )  THEN
410                IF ( simulated_time >= particle_advection_start )  THEN
411                   DO  i = nxl, nxr
412                      DO  j = nys, nyn
413                         DO  k = nzb_do, nzt_do
414                            number_of_particles = prt_count(k,j,i)
415                            IF (number_of_particles <= 0)  CYCLE
416                            particles => grid_particles(k,j,i)%particles(1:number_of_particles)
417                            s_r2 = 0.0_wp
418                            s_r3 = 0.0_wp
419                            DO  n = 1, number_of_particles
420                               IF ( particles(n)%particle_mask )  THEN
421                                  s_r2 = s_r2 + particles(n)%radius**2 * &
422                                         particles(n)%weight_factor
423                                  s_r3 = s_r3 + particles(n)%radius**3 * &
424                                         particles(n)%weight_factor
425                               ENDIF
426                            ENDDO
427                            IF ( s_r2 > 0.0_wp )  THEN
428                               mean_r = s_r3 / s_r2
429                            ELSE
430                               mean_r = 0.0_wp
431                            ENDIF
432                            tend(k,j,i) = mean_r
433                         ENDDO
434                      ENDDO
435                   ENDDO
436                ELSE
437                   tend = 0.0_wp
438                ENDIF
439                DO  i = nxl, nxr
440                   DO  j = nys, nyn
441                      DO  k = nzb_do, nzt_do
442                         local_pf(i,j,k) = tend(k,j,i)
443                      ENDDO
444                   ENDDO
445                ENDDO
446                resorted = .TRUE.
447             ELSE
448                to_be_resorted => pr_av
449             ENDIF
450
451          CASE ( 'prr' )
452             IF ( av == 0 )  THEN
453                DO  i = nxl, nxr
454                   DO  j = nys, nyn
455                      DO  k = nzb_do, nzt_do
456                         local_pf(i,j,k) = prr(k,j,i)
457                      ENDDO
458                   ENDDO
459                ENDDO
460             ELSE
461                DO  i = nxl, nxr
462                   DO  j = nys, nyn
463                      DO  k = nzb_do, nzt_do
464                         local_pf(i,j,k) = prr_av(k,j,i)
465                      ENDDO
466                   ENDDO
467                ENDDO
468             ENDIF
469             resorted = .TRUE.
470
471          CASE ( 'pt' )
472             IF ( av == 0 )  THEN
473                IF ( .NOT. cloud_physics ) THEN
474                   to_be_resorted => pt
475                ELSE
476                   DO  i = nxl, nxr
477                      DO  j = nys, nyn
478                         DO  k = nzb_do, nzt_do
479                            local_pf(i,j,k) = pt(k,j,i) + l_d_cp *             &
480                                                          pt_d_t(k) *          &
481                                                          ql(k,j,i)
482                         ENDDO
483                      ENDDO
484                   ENDDO
485                   resorted = .TRUE.
486                ENDIF
487             ELSE
488                to_be_resorted => pt_av
489             ENDIF
490
491          CASE ( 'q' )
492             IF ( av == 0 )  THEN
493                to_be_resorted => q
494             ELSE
495                to_be_resorted => q_av
496             ENDIF
497
498          CASE ( 'qc' )
499             IF ( av == 0 )  THEN
500                to_be_resorted => qc
501             ELSE
502                to_be_resorted => qc_av
503             ENDIF
504
505          CASE ( 'ql' )
506             IF ( av == 0 )  THEN
507                to_be_resorted => ql
508             ELSE
509                to_be_resorted => ql_av
510             ENDIF
511
512          CASE ( 'ql_c' )
513             IF ( av == 0 )  THEN
514                to_be_resorted => ql_c
515             ELSE
516                to_be_resorted => ql_c_av
517             ENDIF
518
519          CASE ( 'ql_v' )
520             IF ( av == 0 )  THEN
521                to_be_resorted => ql_v
522             ELSE
523                to_be_resorted => ql_v_av
524             ENDIF
525
526          CASE ( 'ql_vp' )
527             IF ( av == 0 )  THEN
528                IF ( simulated_time >= particle_advection_start )  THEN
529                   DO  i = nxl, nxr
530                      DO  j = nys, nyn
531                         DO  k = nzb_do, nzt_do
532                            number_of_particles = prt_count(k,j,i)
533                            IF (number_of_particles <= 0)  CYCLE
534                            particles => grid_particles(k,j,i)%particles(1:number_of_particles)
535                            DO  n = 1, number_of_particles
536                               IF ( particles(n)%particle_mask )  THEN
537                                  tend(k,j,i) =  tend(k,j,i) +                 &
538                                                 particles(n)%weight_factor /  &
539                                                 prt_count(k,j,i)
540                               ENDIF
541                            ENDDO
542                         ENDDO
543                      ENDDO
544                   ENDDO
545                ELSE
546                   tend = 0.0_wp
547                ENDIF
548                DO  i = nxl, nxr
549                   DO  j = nys, nyn
550                      DO  k = nzb_do, nzt_do
551                         local_pf(i,j,k) = tend(k,j,i)
552                      ENDDO
553                   ENDDO
554                ENDDO
555                resorted = .TRUE.
556             ELSE
557                to_be_resorted => ql_vp_av
558             ENDIF
559
560          CASE ( 'qr' )
561             IF ( av == 0 )  THEN
562                to_be_resorted => qr
563             ELSE
564                to_be_resorted => qr_av
565             ENDIF
566
567          CASE ( 'qv' )
568             IF ( av == 0 )  THEN
569                DO  i = nxl, nxr
570                   DO  j = nys, nyn
571                      DO  k = nzb_do, nzt_do
572                         local_pf(i,j,k) = q(k,j,i) - ql(k,j,i)
573                      ENDDO
574                   ENDDO
575                ENDDO
576                resorted = .TRUE.
577             ELSE
578                to_be_resorted => qv_av
579             ENDIF
580
581          CASE ( 'rho_ocean' )
582             IF ( av == 0 )  THEN
583                to_be_resorted => rho_ocean
584             ELSE
585                to_be_resorted => rho_ocean_av
586             ENDIF
587
588          CASE ( 's' )
589             IF ( av == 0 )  THEN
590                to_be_resorted => s
591             ELSE
592                to_be_resorted => s_av
593             ENDIF
594
595          CASE ( 'sa' )
596             IF ( av == 0 )  THEN
597                to_be_resorted => sa
598             ELSE
599                to_be_resorted => sa_av
600             ENDIF
601
602          CASE ( 'u' )
603             flag_nr = 1
604             IF ( av == 0 )  THEN
605                to_be_resorted => u
606             ELSE
607                to_be_resorted => u_av
608             ENDIF
609
610          CASE ( 'v' )
611             flag_nr = 2
612             IF ( av == 0 )  THEN
613                to_be_resorted => v
614             ELSE
615                to_be_resorted => v_av
616             ENDIF
617
618          CASE ( 'vpt' )
619             IF ( av == 0 )  THEN
620                to_be_resorted => vpt
621             ELSE
622                to_be_resorted => vpt_av
623             ENDIF
624
625          CASE ( 'w' )
626             flag_nr = 3
627             IF ( av == 0 )  THEN
628                to_be_resorted => w
629             ELSE
630                to_be_resorted => w_av
631             ENDIF
632!             
633!--       Block of urban surface model outputs   
634          CASE ( 'usm_output' )
635             CALL usm_data_output_3d( av, do3d(av,if), found, local_pf,     &
636                                         nzb_do, nzt_do )
637
638          CASE DEFAULT
639
640!
641!--          Land surface quantity
642             IF ( land_surface )  THEN
643!
644!--             For soil model quantities, it is required to re-allocate local_pf
645                nzb_do = nzb_soil
646                nzt_do = nzt_soil
647
648                DEALLOCATE ( local_pf )
649                ALLOCATE( local_pf(nxl:nxr,nys:nyn,nzb_do:nzt_do) )
650                local_pf = fill_value
651
652                CALL lsm_data_output_3d( av, do3d(av,if), found, local_pf )
653                resorted = .TRUE.
654
655!
656!--             If no soil model variable was found, re-allocate local_pf
657                IF ( .NOT. found )  THEN
658                   nzb_do = nzb
659                   nzt_do = nz_do3d
660
661                   DEALLOCATE ( local_pf )
662                   ALLOCATE( local_pf(nxl:nxr,nys:nyn,nzb_do:nzt_do) )                 
663                ENDIF
664
665             ENDIF
666
667             IF ( .NOT. found )  THEN
668                CALL tcm_data_output_3d( av, do3d(av,if), found, local_pf )
669                resorted = .TRUE.
670             ENDIF
671
672!
673!--          Radiation quantity
674             IF ( .NOT. found  .AND.  radiation )  THEN
675                CALL radiation_data_output_3d( av, do3d(av,if), found,         &
676                                               local_pf )
677                resorted = .TRUE.
678             ENDIF
679
680!
681!--          Chemistry model output
682#if defined( __chem )
683             IF ( .NOT. found  .AND.  air_chemistry )  THEN
684                CALL chem_data_output_3d( av, do3d(av,if), found,              &
685                                          local_pf, fill_value )
686                resorted = .TRUE.
687             ENDIF
688#endif
689
690!
691!--          Plant canopy model output
692             IF ( .NOT. found  .AND.  plant_canopy )  THEN
693                CALL pcm_data_output_3d( av, do3d(av,if), found, local_pf,     &
694                                         fill_value )
695                resorted = .TRUE.
696             ENDIF
697
698!
699!--          User defined quantity
700             IF ( .NOT. found )  THEN
701                CALL user_data_output_3d( av, do3d(av,if), found, local_pf,    &
702                                          nzb_do, nzt_do )
703                resorted = .TRUE.
704             ENDIF
705
706             IF ( .NOT. found )  THEN
707                message_string =  'no output available for: ' //               &
708                                  TRIM( do3d(av,if) )
709                CALL message( 'data_output_3d', 'PA0182', 0, 0, 0, 6, 0 )
710             ENDIF
711
712       END SELECT
713
714!
715!--    Resort the array to be output, if not done above
716       IF ( .NOT. resorted )  THEN
717          DO  i = nxl, nxr
718             DO  j = nys, nyn
719                DO  k = nzb_do, nzt_do
720                   local_pf(i,j,k) = MERGE(                                    &
721                                      to_be_resorted(k,j,i),                   &
722                                      REAL( fill_value, KIND = wp ),           &
723                                      BTEST( wall_flags_0(k,j,i), flag_nr ) )
724                ENDDO
725             ENDDO
726          ENDDO
727       ENDIF
728
729!
730!--    Output of the 3D-array
731#if defined( __parallel )
732       IF ( netcdf_data_format < 5 )  THEN
733!
734!--       Non-parallel netCDF output. Data is output in parallel in
735!--       FORTRAN binary format here, and later collected into one file by
736!--       combine_plot_fields
737          IF ( myid == 0 )  THEN
738             WRITE ( 30 )  time_since_reference_point,                   &
739                           do3d_time_count(av), av
740          ENDIF
741          DO  i = 0, io_blocks-1
742             IF ( i == io_group )  THEN
743                WRITE ( 30 )  nxl, nxr, nys, nyn, nzb_do, nzt_do
744                WRITE ( 30 )  local_pf(:,:,nzb_do:nzt_do)
745             ENDIF
746
747             CALL MPI_BARRIER( comm2d, ierr )
748
749          ENDDO
750
751       ELSE
752#if defined( __netcdf )
753!
754!--       Parallel output in netCDF4/HDF5 format.
755!          IF ( nxr == nx  .AND.  nyn /= ny )  THEN
756!             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
757!                               local_pf(nxl:nxr+1,nys:nyn,nzb_do:nzt_do),    &
758!                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
759!                count = (/ nxr-nxl+2, nyn-nys+1, nzt_do-nzb_do+1, 1 /) )
760!          ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
761!             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
762!                               local_pf(nxl:nxr,nys:nyn+1,nzb_do:nzt_do),    &
763!                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
764!                count = (/ nxr-nxl+1, nyn-nys+2, nzt_do-nzb_do+1, 1 /) )
765!          ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
766!             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
767!                             local_pf(nxl:nxr+1,nys:nyn+1,nzb_do:nzt_do  ),  &
768!                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
769!                count = (/ nxr-nxl+2, nyn-nys+2, nzt_do-nzb_do+1, 1 /) )
770!          ELSE
771             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
772                                 local_pf(nxl:nxr,nys:nyn,nzb_do:nzt_do),    &
773                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
774                count = (/ nxr-nxl+1, nyn-nys+1, nzt_do-nzb_do+1, 1 /) )
775!          ENDIF
776          CALL netcdf_handle_error( 'data_output_3d', 386 )
777#endif
778       ENDIF
779#else
780#if defined( __netcdf )
781       nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),        &
782                         local_pf(nxl:nxr+1,nys:nyn+1,nzb_do:nzt_do),        &
783                         start = (/ 1, 1, 1, do3d_time_count(av) /),     &
784                         count = (/ nx+1, ny+1, nzt_do-nzb_do+1, 1 /) )
785       CALL netcdf_handle_error( 'data_output_3d', 446 )
786#endif
787#endif
788
789       if = if + 1
790
791!
792!--    Deallocate temporary array
793       DEALLOCATE ( local_pf )
794
795    ENDDO
796
797    CALL cpu_log( log_point(14), 'data_output_3d', 'stop' )
798
799!
800!-- Formats.
8013300 FORMAT ('variable ',I4,'  file=',A,'  filetype=unformatted  skip=',I12/   &
802             'label = ',A,A)
803
804 END SUBROUTINE data_output_3d
Note: See TracBrowser for help on using the repository browser.