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

Last change on this file since 3045 was 3045, checked in by Giersch, 3 years ago

Code adjusted according to coding standards, renamed namelists, error messages revised until PA0347, output CASE 108 disabled

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