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

Last change on this file since 4155 was 4155, checked in by raasch, 2 years ago

bugfix for 3d-output in serial mode (ghost points must not be written)

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