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

Last change on this file since 3924 was 3885, checked in by kanani, 6 years ago

restructure/add location/debug messages

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