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

Last change on this file since 3994 was 3994, checked in by suehring, 2 years ago

new module for diagnostic output quantities added + output of turbulence intensity

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