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

Last change on this file since 3031 was 3014, checked in by maronga, 6 years ago

series of bugfixes

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