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

Last change on this file since 3298 was 3294, checked in by raasch, 6 years ago

modularization of the ocean code

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