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

Last change on this file since 2798 was 2766, checked in by kanani, 6 years ago

Removal of chem directive, plus minor changes

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