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

Last change on this file since 3435 was 3421, checked in by gronemeier, 6 years ago

new surface-data output; renamed output variables (pt to theta, rho_air to rho, rho_ocean to rho_sea_water)

  • Property svn:keywords set to Id
File size: 31.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 3421 2018-10-24 18:39:32Z gronemeier $
27! Renamed output variables
28!
29! 3419 2018-10-24 17:27:31Z gronemeier
30! bugfix: nx, ny are required in non-parallel case
31!
32! 3337 2018-10-12 15:17:09Z kanani
33! (from branch resler)
34! Add Biometeorology
35!
36! 3294 2018-10-01 02:37:10Z raasch
37! changes concerning modularization of ocean option
38!
39! 3274 2018-09-24 15:42:55Z knoop
40! Modularization of all bulk cloud physics code components
41!
42! 3241 2018-09-12 15:02:00Z raasch
43! unused variables and format statements removed
44!
45! 3049 2018-05-29 13:52:36Z Giersch
46! Error messages revised
47!
48! 3045 2018-05-28 07:55:41Z Giersch
49! Error message revised
50!
51! 3014 2018-05-09 08:42:38Z maronga
52! Added nzb_do and nzt_do for some modules for 3d output
53!
54! 3004 2018-04-27 12:33:25Z Giersch
55! Allocation checks implemented (averaged data will be assigned to fill values
56! if no allocation happened so far)
57!
58! 2967 2018-04-13 11:22:08Z raasch
59! bugfix: missing parallel cpp-directives added
60!
61! 2817 2018-02-19 16:32:21Z knoop
62! Preliminary gust module interface implemented
63!
64! 2766 2018-01-22 17:17:47Z kanani
65! Removed preprocessor directive __chem
66!
67! 2756 2018-01-16 18:11:14Z suehring
68! Fill values for 3D output of chemical species introduced.
69!
70! 2746 2018-01-15 12:06:04Z suehring
71! Move flag plant canopy to modules
72!
73! 2718 2018-01-02 08:49:38Z maronga
74! Corrected "Former revisions" section
75!
76! 2696 2017-12-14 17:12:51Z kanani
77! Change in file header (GPL part)
78! Implementation of turbulence_closure_mod (TG)
79! Implementation of chemistry module (FK)
80! Set fill values at topography grid points or e.g. non-natural-type surface
81! in case of LSM output (MS)
82!
83! 2512 2017-10-04 08:26:59Z raasch
84! upper bounds of 3d output changed from nx+1,ny+1 to nx,ny
85! no output of ghost layer data
86!
87! 2292 2017-06-20 09:51:42Z schwenkel
88! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
89! includes two more prognostic equations for cloud drop concentration (nc) 
90! and cloud water content (qc).
91!
92! 2233 2017-05-30 18:08:54Z suehring
93!
94! 2232 2017-05-30 17:47:52Z suehring
95! Adjustments to new topography concept
96!
97! 2209 2017-04-19 09:34:46Z kanani
98! Added plant canopy model output
99!
100! 2031 2016-10-21 15:11:58Z knoop
101! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
102!
103! 2011 2016-09-19 17:29:57Z kanani
104! Flag urban_surface is now defined in module control_parameters,
105! changed prefix for urban surface model output to "usm_",
106! introduced control parameter varnamelength for LEN of trimvar.
107!
108! 2007 2016-08-24 15:47:17Z kanani
109! Added support for new urban surface model (temporary modifications of
110! SELECT CASE ( ) necessary, see variable trimvar)
111!
112! 2000 2016-08-20 18:09:15Z knoop
113! Forced header and separation lines into 80 columns
114!
115! 1980 2016-07-29 15:51:57Z suehring
116! Bugfix, in order to steer user-defined output, setting flag found explicitly
117! to .F.
118!
119! 1976 2016-07-27 13:28:04Z maronga
120! Output of radiation quantities is now done directly in the respective module
121!
122! 1972 2016-07-26 07:52:02Z maronga
123! Output of land surface quantities is now done directly in the respective module.
124! Unnecessary directive __parallel removed.
125!
126! 1960 2016-07-12 16:34:24Z suehring
127! Scalar surface flux added
128!
129! 1849 2016-04-08 11:33:18Z hoffmann
130! prr moved to arrays_3d
131!
132! 1822 2016-04-07 07:49:42Z hoffmann
133! prr vertical dimensions set to nzb_do to nzt_do. Unused variables deleted.
134!
135! 1808 2016-04-05 19:44:00Z raasch
136! test output removed
137!
138! 1783 2016-03-06 18:36:17Z raasch
139! name change of netcdf routines and module + related changes
140!
141! 1745 2016-02-05 13:06:51Z gronemeier
142! Bugfix: test if time axis limit exceeds moved to point after call of check_open
143!
144! 1691 2015-10-26 16:17:44Z maronga
145! Added output of radiative heating rates for RRTMG
146!
147! 1682 2015-10-07 23:56:08Z knoop
148! Code annotations made doxygen readable
149!
150! 1585 2015-04-30 07:05:52Z maronga
151! Added support for RRTMG
152!
153! 1551 2015-03-03 14:18:16Z maronga
154! Added suppport for land surface model and radiation model output. In the course
155! of this action, the limits for vertical loops have been changed (from nzb and
156! nzt+1 to nzb_do and nzt_do, respectively in order to allow soil model output).
157! Moreover, a new vertical grid zs was introduced.
158!
159! 1359 2014-04-11 17:15:14Z hoffmann
160! New particle structure integrated.
161!
162! 1353 2014-04-08 15:21:23Z heinze
163! REAL constants provided with KIND-attribute
164!
165! 1327 2014-03-21 11:00:16Z raasch
166! parts concerning avs output removed,
167! -netcdf output queries
168!
169! 1320 2014-03-20 08:40:49Z raasch
170! ONLY-attribute added to USE-statements,
171! kind-parameters added to all INTEGER and REAL declaration statements,
172! kinds are defined in new module kinds,
173! old module precision_kind is removed,
174! revision history before 2012 removed,
175! comment fields (!:) to be used for variable explanations added to
176! all variable declaration statements
177!
178! 1318 2014-03-17 13:35:16Z raasch
179! barrier argument removed from cpu_log,
180! module interfaces removed
181!
182! 1308 2014-03-13 14:58:42Z fricke
183! Check, if the limit of the time dimension is exceeded for parallel output
184! To increase the performance for parallel output, the following is done:
185! - Update of time axis is only done by PE0
186!
187! 1244 2013-10-31 08:16:56Z raasch
188! Bugfix for index bounds in case of 3d-parallel output
189!
190! 1115 2013-03-26 18:16:16Z hoffmann
191! ql is calculated by calc_liquid_water_content
192!
193! 1106 2013-03-04 05:31:38Z raasch
194! array_kind renamed precision_kind
195!
196! 1076 2012-12-05 08:30:18Z hoffmann
197! Bugfix in output of ql
198!
199! 1053 2012-11-13 17:11:03Z hoffmann
200! +nr, qr, prr, qc and averaged quantities
201!
202! 1036 2012-10-22 13:43:42Z raasch
203! code put under GPL (PALM 3.9)
204!
205! 1031 2012-10-19 14:35:30Z raasch
206! netCDF4 without parallel file support implemented
207!
208! 1007 2012-09-19 14:30:36Z franke
209! Bugfix: missing calculation of ql_vp added
210!
211! Revision 1.1  1997/09/03 06:29:36  raasch
212! Initial revision
213!
214!
215! Description:
216! ------------
217!> Output of the 3D-arrays in netCDF and/or AVS format.
218!------------------------------------------------------------------------------!
219 SUBROUTINE data_output_3d( av )
220 
221
222    USE arrays_3d,                                                             &
223        ONLY:  d_exner, e, nc, nr, p, pt, prr, q, qc, ql, ql_c, ql_v, qr, s,   &
224               tend, u, v, vpt, w
225
226    USE averaging
227
228    USE basic_constants_and_equations_mod,                                     &
229        ONLY:  lv_d_cp
230
231    USE bulk_cloud_model_mod,                                                  &
232        ONLY:  bulk_cloud_model, bcm_data_output_3d
233
234    USE chemistry_model_mod,                                                   &
235        ONLY:  chem_data_output_3d
236
237    USE control_parameters,                                                    &
238        ONLY:  air_chemistry, do3d, do3d_no, do3d_time_count,                  &
239               io_blocks, io_group, land_surface, message_string,              &
240               ntdim_3d, nz_do3d,  ocean_mode, plant_canopy,                   &
241               psolver, simulated_time, time_since_reference_point,            &
242               urban_surface, varnamelength
243
244    USE cpulog,                                                                &
245        ONLY:  log_point, cpu_log
246
247    USE gust_mod,                                                              &
248        ONLY: gust_data_output_3d, gust_module_enabled
249
250#if defined( __parallel )
251    USE indices,                                                               &
252        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,     &
253               wall_flags_0
254#else
255    USE indices,                                                               &
256        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb,  &
257               nzt, wall_flags_0
258#endif
259
260    USE kinds
261
262    USE land_surface_model_mod,                                                &
263        ONLY: lsm_data_output_3d, nzb_soil, nzt_soil
264
265#if defined( __netcdf )
266    USE NETCDF
267#endif
268
269    USE netcdf_interface,                                                      &
270        ONLY:  fill_value, id_set_3d, id_var_do3d, id_var_time_3d, nc_stat,    &
271               netcdf_data_format, netcdf_handle_error
272       
273    USE ocean_mod,                                                             &
274        ONLY:  ocean_data_output_3d
275
276    USE particle_attributes,                                                   &
277        ONLY:  grid_particles, number_of_particles, particles,                 &
278               particle_advection_start, prt_count
279       
280    USE pegrid
281
282    USE plant_canopy_model_mod,                                                &
283        ONLY:  pcm_data_output_3d
284       
285    USE radiation_model_mod,                                                   &
286        ONLY:  nzub, nzut, radiation, radiation_data_output_3d
287
288    USE biometeorology_mod,                                                    &
289        ONLY: biometeorology_data_output_3d
290
291    USE turbulence_closure_mod,                                                &
292        ONLY:  tcm_data_output_3d
293
294    USE urban_surface_mod,                                                     &
295        ONLY:  usm_data_output_3d
296
297
298    IMPLICIT NONE
299
300    INTEGER(iwp) ::  av        !<
301    INTEGER(iwp) ::  flag_nr   !< number of masking flag
302    INTEGER(iwp) ::  i         !<
303    INTEGER(iwp) ::  if        !<
304    INTEGER(iwp) ::  j         !<
305    INTEGER(iwp) ::  k         !<
306    INTEGER(iwp) ::  n         !<
307    INTEGER(iwp) ::  nzb_do    !< vertical lower limit for data output
308    INTEGER(iwp) ::  nzt_do    !< vertical upper limit for data output
309
310    LOGICAL      ::  found     !<
311    LOGICAL      ::  resorted  !<
312
313    REAL(wp)     ::  mean_r    !<
314    REAL(wp)     ::  s_r2      !<
315    REAL(wp)     ::  s_r3      !<
316
317    REAL(sp), DIMENSION(:,:,:), ALLOCATABLE ::  local_pf  !<
318
319    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !<
320
321    CHARACTER (LEN=varnamelength) ::  trimvar  !< TRIM of output-variable string
322
323!
324!-- Return, if nothing to output
325    IF ( do3d_no(av) == 0 )  RETURN
326
327    CALL cpu_log (log_point(14),'data_output_3d','start')
328
329!
330!-- Open output file.
331!-- For classic or 64bit netCDF output on more than one PE, each PE opens its
332!-- own file and writes the data of its subdomain in binary format. After the
333!-- run, these files are combined to one NetCDF file by combine_plot_fields.
334!-- For netCDF4/HDF5 output, data is written in parallel into one file.
335    IF ( netcdf_data_format < 5 )  THEN
336#if defined( __parallel )
337       CALL check_open( 30 )
338#endif
339       IF ( myid == 0 )  CALL check_open( 106+av*10 )
340    ELSE
341       CALL check_open( 106+av*10 )
342    ENDIF
343
344!
345!-- For parallel netcdf output the time axis must be limited. Return, if this
346!-- limit is exceeded. This could be the case, if the simulated time exceeds
347!-- the given end time by the length of the given output interval.
348    IF ( netcdf_data_format > 4 )  THEN
349       IF ( do3d_time_count(av) + 1 > ntdim_3d(av) )  THEN
350          WRITE ( message_string, * ) 'Output of 3d data is not given at t=',  &
351                                      simulated_time, '&because the maximum ', & 
352                                      'number of output time levels is ',      &
353                                      'exceeded.'
354          CALL message( 'data_output_3d', 'PA0387', 0, 1, 0, 6, 0 )
355          CALL cpu_log( log_point(14), 'data_output_3d', 'stop' )
356          RETURN
357       ENDIF
358    ENDIF
359
360!
361!-- Update the netCDF time axis
362!-- In case of parallel output, this is only done by PE0 to increase the
363!-- performance.
364#if defined( __netcdf )
365    do3d_time_count(av) = do3d_time_count(av) + 1
366    IF ( myid == 0 )  THEN
367       nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_time_3d(av),           &
368                               (/ time_since_reference_point /),            &
369                               start = (/ do3d_time_count(av) /),           &
370                               count = (/ 1 /) )
371       CALL netcdf_handle_error( 'data_output_3d', 376 )
372    ENDIF
373#endif
374
375!
376!-- Loop over all variables to be written.
377    if = 1
378
379    DO  WHILE ( do3d(av,if)(1:1) /= ' ' )
380
381!
382!--    Temporary solution to account for data output within the new urban
383!--    surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar ).
384!--    Store the array chosen on the temporary array.
385       trimvar = TRIM( do3d(av,if) )
386       IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
387          trimvar = 'usm_output'
388          resorted = .TRUE.
389          nzb_do   = nzub
390          nzt_do   = nzut
391       ELSE
392          resorted = .FALSE.
393          nzb_do   = nzb
394          nzt_do   = nz_do3d
395       ENDIF
396!
397!--    Set flag to steer output of radiation, land-surface, or user-defined
398!--    quantities
399       found = .FALSE.
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!--       Block of urban surface model outputs   
700          CASE ( 'usm_output' )
701             CALL usm_data_output_3d( av, do3d(av,if), found, local_pf,        &
702                                         nzb_do, nzt_do )
703
704          CASE DEFAULT
705
706!
707!--          Quantities of other modules
708             IF ( .NOT. found  .AND.  bulk_cloud_model )  THEN
709                CALL bcm_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.  air_chemistry )  THEN
715                CALL chem_data_output_3d( av, do3d(av,if), found,              &
716                                          local_pf, fill_value, nzb_do, nzt_do )
717                resorted = .TRUE.
718             ENDIF
719
720             IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
721                CALL gust_data_output_3d( av, do3d(av,if), found, local_pf,    &
722                                          nzb_do, nzt_do )
723                resorted = .TRUE.
724             ENDIF
725
726             IF ( .NOT. found  .AND.  land_surface )  THEN
727!
728!--             For soil model quantities, it is required to re-allocate local_pf
729                nzb_do = nzb_soil
730                nzt_do = nzt_soil
731
732                DEALLOCATE ( local_pf )
733                ALLOCATE( local_pf(nxl:nxr,nys:nyn,nzb_do:nzt_do) )
734                local_pf = fill_value
735
736                CALL lsm_data_output_3d( av, do3d(av,if), found, local_pf )
737                resorted = .TRUE.
738
739!
740!--             If no soil model variable was found, re-allocate local_pf
741                IF ( .NOT. found )  THEN
742                   nzb_do = nzb
743                   nzt_do = nz_do3d
744
745                   DEALLOCATE ( local_pf )
746                   ALLOCATE( local_pf(nxl:nxr,nys:nyn,nzb_do:nzt_do) )                 
747                ENDIF
748
749             ENDIF
750
751             IF ( .NOT. found  .AND.  ocean_mode )  THEN
752                CALL ocean_data_output_3d( av, do3d(av,if), found, local_pf,   &
753                                           nzb_do, nzt_do )
754                resorted = .TRUE.
755             ENDIF
756
757             IF ( .NOT. found  .AND.  plant_canopy )  THEN
758                CALL pcm_data_output_3d( av, do3d(av,if), found, local_pf,     &
759                                         fill_value, nzb_do, nzt_do )
760                resorted = .TRUE.
761             ENDIF
762
763             IF ( .NOT. found  .AND.  radiation )  THEN
764                CALL radiation_data_output_3d( av, do3d(av,if), found,         &
765                                               local_pf, nzb_do, nzt_do )
766                resorted = .TRUE.
767             ENDIF
768
769!
770!--          Biometeorology quantity
771             IF ( .NOT. found  .AND.  radiation )  THEN
772                CALL biometeorology_data_output_3d( av, do3d(av,if), found,    &
773                                               local_pf, nzb_do, nzt_do )
774                resorted = .TRUE.
775             ENDIF
776
777             IF ( .NOT. found )  THEN
778                CALL tcm_data_output_3d( av, do3d(av,if), found, local_pf,     &
779                                         nzb_do, nzt_do )
780                resorted = .TRUE.
781             ENDIF
782
783!
784!--          User defined quantities
785             IF ( .NOT. found )  THEN
786                CALL user_data_output_3d( av, do3d(av,if), found, local_pf,    &
787                                          nzb_do, nzt_do )
788                resorted = .TRUE.
789             ENDIF
790
791             IF ( .NOT. found )  THEN
792                message_string =  'no output available for: ' //               &
793                                  TRIM( do3d(av,if) )
794                CALL message( 'data_output_3d', 'PA0182', 0, 0, 0, 6, 0 )
795             ENDIF
796
797       END SELECT
798
799!
800!--    Resort the array to be output, if not done above
801       IF ( .NOT. resorted )  THEN
802          DO  i = nxl, nxr
803             DO  j = nys, nyn
804                DO  k = nzb_do, nzt_do
805                   local_pf(i,j,k) = MERGE(                                    &
806                                      to_be_resorted(k,j,i),                   &
807                                      REAL( fill_value, KIND = wp ),           &
808                                      BTEST( wall_flags_0(k,j,i), flag_nr ) )
809                ENDDO
810             ENDDO
811          ENDDO
812       ENDIF
813
814!
815!--    Output of the 3D-array
816#if defined( __parallel )
817       IF ( netcdf_data_format < 5 )  THEN
818!
819!--       Non-parallel netCDF output. Data is output in parallel in
820!--       FORTRAN binary format here, and later collected into one file by
821!--       combine_plot_fields
822          IF ( myid == 0 )  THEN
823             WRITE ( 30 )  time_since_reference_point,                   &
824                           do3d_time_count(av), av
825          ENDIF
826          DO  i = 0, io_blocks-1
827             IF ( i == io_group )  THEN
828                WRITE ( 30 )  nxl, nxr, nys, nyn, nzb_do, nzt_do
829                WRITE ( 30 )  local_pf(:,:,nzb_do:nzt_do)
830             ENDIF
831
832             CALL MPI_BARRIER( comm2d, ierr )
833
834          ENDDO
835
836       ELSE
837#if defined( __netcdf )
838!
839!--       Parallel output in netCDF4/HDF5 format.
840!          IF ( nxr == nx  .AND.  nyn /= ny )  THEN
841!             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
842!                               local_pf(nxl:nxr+1,nys:nyn,nzb_do:nzt_do),    &
843!                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
844!                count = (/ nxr-nxl+2, nyn-nys+1, nzt_do-nzb_do+1, 1 /) )
845!          ELSEIF ( nxr /= nx  .AND.  nyn == ny )  THEN
846!             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
847!                               local_pf(nxl:nxr,nys:nyn+1,nzb_do:nzt_do),    &
848!                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
849!                count = (/ nxr-nxl+1, nyn-nys+2, nzt_do-nzb_do+1, 1 /) )
850!          ELSEIF ( nxr == nx  .AND.  nyn == ny )  THEN
851!             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
852!                             local_pf(nxl:nxr+1,nys:nyn+1,nzb_do:nzt_do  ),  &
853!                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
854!                count = (/ nxr-nxl+2, nyn-nys+2, nzt_do-nzb_do+1, 1 /) )
855!          ELSE
856             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
857                                 local_pf(nxl:nxr,nys:nyn,nzb_do:nzt_do),    &
858                start = (/ nxl+1, nys+1, nzb_do+1, do3d_time_count(av) /),  &
859                count = (/ nxr-nxl+1, nyn-nys+1, nzt_do-nzb_do+1, 1 /) )
860!          ENDIF
861          CALL netcdf_handle_error( 'data_output_3d', 386 )
862#endif
863       ENDIF
864#else
865#if defined( __netcdf )
866       nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),        &
867                         local_pf(nxl:nxr+1,nys:nyn+1,nzb_do:nzt_do),        &
868                         start = (/ 1, 1, 1, do3d_time_count(av) /),     &
869                         count = (/ nx+1, ny+1, nzt_do-nzb_do+1, 1 /) )
870       CALL netcdf_handle_error( 'data_output_3d', 446 )
871#endif
872#endif
873
874       if = if + 1
875
876!
877!--    Deallocate temporary array
878       DEALLOCATE ( local_pf )
879
880    ENDDO
881
882    CALL cpu_log( log_point(14), 'data_output_3d', 'stop' )
883
884 END SUBROUTINE data_output_3d
Note: See TracBrowser for help on using the repository browser.