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

Last change on this file since 3475 was 3467, checked in by suehring, 6 years ago

Branch salsa @3446 re-integrated into trunk

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