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

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

bugfix for r4155

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