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

Last change on this file since 3641 was 3637, checked in by knoop, 5 years ago

M Makefile

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