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

Last change on this file since 4111 was 4048, checked in by knoop, 5 years ago

Moved turbulence_closure_mod calls into module_interface

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