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

Last change on this file since 3418 was 3405, checked in by raasch, 6 years ago

bugfix: BIND attribute added to derived type particle_type, bugfix: nx, ny are required in non-parallel case

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