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

Last change on this file since 3138 was 3049, checked in by Giersch, 7 years ago

Revision history corrected

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