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

Last change on this file since 3383 was 3337, checked in by kanani, 6 years ago

reintegrate branch resler to trunk

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