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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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