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

Last change on this file since 2993 was 2967, checked in by raasch, 7 years ago

bugfix: missing parallel cpp-directives added

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