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

Last change on this file since 3589 was 3589, checked in by suehring, 3 years ago

Remove erroneous UTF encoding; last commit documented

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