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

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

Merge with branch resler: biomet- output of bio_mrt added; plant_canopy - separate vertical dimension for 3D output (to save disk space); radiation - remove unused plant canopy variables; urban-surface model - do not add anthropogenic heat during wall spin-up

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