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

Last change on this file since 3525 was 3525, checked in by kanani, 4 years ago

Changes related to clean-up of biometeorology (by dom_dwd_user)

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