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

Last change on this file since 4039 was 4039, checked in by suehring, 5 years ago

diagnostic output: Modularize diagnostic output, rename subroutines; formatting adjustments; allocate arrays only when required; add output of uu, vv, ww to enable variance calculation via temporal EC method; radiation: bugfix in masked data output; flow_statistics: Correct conversion to kinematic vertical scalar fluxes in case of pw-scheme and statistic regions

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