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

Last change on this file since 3582 was 3582, checked in by suehring, 5 years ago

Merge branch salsa with trunk

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