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

Last change on this file since 3579 was 3554, checked in by gronemeier, 6 years ago

renamed variable if to ivar; add variable description

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