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

Last change on this file since 3676 was 3655, checked in by knoop, 6 years ago

Bugfix: made "unit" and "found" intend INOUT in module interface subroutines + automatic copyright update

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