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

Last change on this file since 3987 was 3987, checked in by kanani, 2 years ago

clean up location, debug and error 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 3987 2019-05-22 09:52:13Z kanani $
27! Introduce alternative switch for debug output during timestepping
28!
29! 3885 2019-04-11 11:29:34Z kanani
30! Changes related to global restructuring of location messages and introduction
31! of additional debug messages
32!
33! 3814 2019-03-26 08:40:31Z pavelkrc
34! unused variables removed
35!
36! 3655 2019-01-07 16:51:22Z knoop
37! Bugfix: use time_since_reference_point instead of simulated_time (relevant
38! when using wall/soil spinup)
39!
40! 3637 2018-12-20 01:51:36Z knoop
41! Implementation of the PALM module interface
42!
43! 3589 2018-11-30 15:09:51Z suehring
44! Move the control parameter "salsa" from salsa_mod to control_parameters
45! (M. Kurppa)
46!
47! 3582 2018-11-29 19:16:36Z suehring
48! add variable description; rename variable 'if' into 'ivar'
49!
50! 3525 2018-11-14 16:06:14Z kanani
51! Changes related to clean-up of biometeorology (dom_dwd_user)
52!
53! 3467 2018-10-30 19:05:21Z suehring
54! Implementation of a new aerosol module salsa.
55!
56! 3448 2018-10-29 18:14:31Z kanani
57! Adjustment of biometeorology calls
58!
59! 3421 2018-10-24 18:39:32Z gronemeier
60! Renamed output variables
61!
62! 3419 2018-10-24 17:27:31Z gronemeier
63! bugfix: nx, ny are required in non-parallel case
64!
65! 3337 2018-10-12 15:17:09Z kanani
66! (from branch resler)
67! Add Biometeorology
68!
69! 3294 2018-10-01 02:37:10Z raasch
70! changes concerning modularization of ocean option
71!
72! 3274 2018-09-24 15:42:55Z knoop
73! Modularization of all bulk cloud physics code components
74!
75! 3241 2018-09-12 15:02:00Z raasch
76! unused variables and format statements removed
77!
78! 3049 2018-05-29 13:52:36Z Giersch
79! Error messages revised
80!
81! 3045 2018-05-28 07:55:41Z Giersch
82! Error message revised
83!
84! 3014 2018-05-09 08:42:38Z maronga
85! Added nzb_do and nzt_do for some modules for 3d output
86!
87! 3004 2018-04-27 12:33:25Z Giersch
88! Allocation checks implemented (averaged data will be assigned to fill values
89! if no allocation happened so far)
90!
91! 2967 2018-04-13 11:22:08Z raasch
92! bugfix: missing parallel cpp-directives added
93!
94! 2817 2018-02-19 16:32:21Z knoop
95! Preliminary gust module interface implemented
96!
97! 2766 2018-01-22 17:17:47Z kanani
98! Removed preprocessor directive __chem
99!
100! 2756 2018-01-16 18:11:14Z suehring
101! Fill values for 3D output of chemical species introduced.
102!
103! 2746 2018-01-15 12:06:04Z suehring
104! Move flag plant canopy to modules
105!
106! 2718 2018-01-02 08:49:38Z maronga
107! Corrected "Former revisions" section
108!
109! 2696 2017-12-14 17:12:51Z kanani
110! Change in file header (GPL part)
111! Implementation of turbulence_closure_mod (TG)
112! Implementation of chemistry module (FK)
113! Set fill values at topography grid points or e.g. non-natural-type surface
114! in case of LSM output (MS)
115!
116! 2512 2017-10-04 08:26:59Z raasch
117! upper bounds of 3d output changed from nx+1,ny+1 to nx,ny
118! no output of ghost layer data
119!
120! 2292 2017-06-20 09:51:42Z schwenkel
121! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
122! includes two more prognostic equations for cloud drop concentration (nc) 
123! and cloud water content (qc).
124!
125! 2233 2017-05-30 18:08:54Z suehring
126!
127! 2232 2017-05-30 17:47:52Z suehring
128! Adjustments to new topography concept
129!
130! 2209 2017-04-19 09:34:46Z kanani
131! Added plant canopy model output
132!
133! 2031 2016-10-21 15:11:58Z knoop
134! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
135!
136! 2011 2016-09-19 17:29:57Z kanani
137! Flag urban_surface is now defined in module control_parameters,
138! changed prefix for urban surface model output to "usm_",
139! introduced control parameter varnamelength for LEN of trimvar.
140!
141! 2007 2016-08-24 15:47:17Z kanani
142! Added support for new urban surface model (temporary modifications of
143! SELECT CASE ( ) necessary, see variable trimvar)
144!
145! 2000 2016-08-20 18:09:15Z knoop
146! Forced header and separation lines into 80 columns
147!
148! 1980 2016-07-29 15:51:57Z suehring
149! Bugfix, in order to steer user-defined output, setting flag found explicitly
150! to .F.
151!
152! 1976 2016-07-27 13:28:04Z maronga
153! Output of radiation quantities is now done directly in the respective module
154!
155! 1972 2016-07-26 07:52:02Z maronga
156! Output of land surface quantities is now done directly in the respective module.
157! Unnecessary directive __parallel removed.
158!
159! 1960 2016-07-12 16:34:24Z suehring
160! Scalar surface flux added
161!
162! 1849 2016-04-08 11:33:18Z hoffmann
163! prr moved to arrays_3d
164!
165! 1822 2016-04-07 07:49:42Z hoffmann
166! prr vertical dimensions set to nzb_do to nzt_do. Unused variables deleted.
167!
168! 1808 2016-04-05 19:44:00Z raasch
169! test output removed
170!
171! 1783 2016-03-06 18:36:17Z raasch
172! name change of netcdf routines and module + related changes
173!
174! 1745 2016-02-05 13:06:51Z gronemeier
175! Bugfix: test if time axis limit exceeds moved to point after call of check_open
176!
177! 1691 2015-10-26 16:17:44Z maronga
178! Added output of radiative heating rates for RRTMG
179!
180! 1682 2015-10-07 23:56:08Z knoop
181! Code annotations made doxygen readable
182!
183! 1585 2015-04-30 07:05:52Z maronga
184! Added support for RRTMG
185!
186! 1551 2015-03-03 14:18:16Z maronga
187! Added suppport for land surface model and radiation model output. In the course
188! of this action, the limits for vertical loops have been changed (from nzb and
189! nzt+1 to nzb_do and nzt_do, respectively in order to allow soil model output).
190! Moreover, a new vertical grid zs was introduced.
191!
192! 1359 2014-04-11 17:15:14Z hoffmann
193! New particle structure integrated.
194!
195! 1353 2014-04-08 15:21:23Z heinze
196! REAL constants provided with KIND-attribute
197!
198! 1327 2014-03-21 11:00:16Z raasch
199! parts concerning avs output removed,
200! -netcdf output queries
201!
202! 1320 2014-03-20 08:40:49Z raasch
203! ONLY-attribute added to USE-statements,
204! kind-parameters added to all INTEGER and REAL declaration statements,
205! kinds are defined in new module kinds,
206! old module precision_kind is removed,
207! revision history before 2012 removed,
208! comment fields (!:) to be used for variable explanations added to
209! all variable declaration statements
210!
211! 1318 2014-03-17 13:35:16Z raasch
212! barrier argument removed from cpu_log,
213! module interfaces removed
214!
215! 1308 2014-03-13 14:58:42Z fricke
216! Check, if the limit of the time dimension is exceeded for parallel output
217! To increase the performance for parallel output, the following is done:
218! - Update of time axis is only done by PE0
219!
220! 1244 2013-10-31 08:16:56Z raasch
221! Bugfix for index bounds in case of 3d-parallel output
222!
223! 1115 2013-03-26 18:16:16Z hoffmann
224! ql is calculated by calc_liquid_water_content
225!
226! 1106 2013-03-04 05:31:38Z raasch
227! array_kind renamed precision_kind
228!
229! 1076 2012-12-05 08:30:18Z hoffmann
230! Bugfix in output of ql
231!
232! 1053 2012-11-13 17:11:03Z hoffmann
233! +nr, qr, prr, qc and averaged quantities
234!
235! 1036 2012-10-22 13:43:42Z raasch
236! code put under GPL (PALM 3.9)
237!
238! 1031 2012-10-19 14:35:30Z raasch
239! netCDF4 without parallel file support implemented
240!
241! 1007 2012-09-19 14:30:36Z franke
242! Bugfix: missing calculation of ql_vp added
243!
244! Revision 1.1  1997/09/03 06:29:36  raasch
245! Initial revision
246!
247!
248! Description:
249! ------------
250!> Output of the 3D-arrays in netCDF and/or AVS format.
251!------------------------------------------------------------------------------!
252 SUBROUTINE data_output_3d( av )
253 
254
255    USE arrays_3d,                                                             &
256        ONLY:  d_exner, e, p, pt, q, ql, ql_c, ql_v, s, tend, u, v, vpt, w
257
258    USE averaging
259
260    USE basic_constants_and_equations_mod,                                     &
261        ONLY:  lv_d_cp
262
263    USE bulk_cloud_model_mod,                                                  &
264        ONLY:  bulk_cloud_model
265
266    USE control_parameters,                                                    &
267        ONLY:  debug_output_timestep,                                          &
268               do3d, do3d_no, do3d_time_count, io_blocks, io_group,            &
269               land_surface, message_string, ntdim_3d, nz_do3d, psolver,       &
270               time_since_reference_point, urban_surface, varnamelength
271
272    USE cpulog,                                                                &
273        ONLY:  log_point, cpu_log
274
275#if defined( __parallel )
276    USE indices,                                                               &
277        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,     &
278               wall_flags_0
279#else
280    USE indices,                                                               &
281        ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb,  &
282               nzt, wall_flags_0
283#endif
284
285    USE kinds
286
287    USE land_surface_model_mod,                                                &
288        ONLY: lsm_data_output_3d, nzb_soil, nzt_soil
289
290    USE module_interface,                                                      &
291        ONLY:  module_interface_data_output_3d
292
293#if defined( __netcdf )
294    USE NETCDF
295#endif
296
297    USE netcdf_interface,                                                      &
298        ONLY:  fill_value, id_set_3d, id_var_do3d, id_var_time_3d, nc_stat,    &
299               netcdf_data_format, netcdf_handle_error
300
301    USE particle_attributes,                                                   &
302        ONLY:  grid_particles, number_of_particles, particles,                 &
303               particle_advection_start, prt_count
304
305    USE pegrid
306
307    USE radiation_model_mod,                                                   &
308        ONLY:  nz_urban_b, nz_urban_t
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    IF ( debug_output_timestep )  CALL debug_message( 'data_output_3d', 'start' )
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_timestep )  CALL debug_message( 'data_output_3d', 'end' )
879
880
881 END SUBROUTINE data_output_3d
Note: See TracBrowser for help on using the repository browser.