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

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

Read information from statitic driver for resolved vegetation independently from land- or urban-surface model

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