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

Last change on this file since 3456 was 3448, checked in by kanani, 6 years ago

Implementation of human thermal indices (from branch biomet_p2 at r3444)

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