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

Last change on this file since 3004 was 3004, checked in by Giersch, 6 years ago

precipitation_rate removed, further allocation checks for data output of averaged quantities implemented, double CALL of flow_statistics at the beginning of time_integration removed, further minor bugfixes, comments added

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