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

Last change on this file since 3250 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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