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

Last change on this file since 3465 was 3448, checked in by kanani, 5 years ago

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

  • Property svn:keywords set to Id
File size: 31.6 KB
RevLine 
[1682]1!> @file data_output_3d.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[2000]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.
[1036]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!
[2718]17! Copyright 1997-2018 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[254]20! Current revisions:
[1106]21! ------------------
[1961]22!
[3049]23!
[1552]24! Former revisions:
25! -----------------
26! $Id: data_output_3d.f90 3448 2018-10-29 18:14:31Z kanani $
[3448]27! Adjustment of biometeorology calls
28!
29! 3421 2018-10-24 18:39:32Z gronemeier
[3421]30! Renamed output variables
31!
32! 3419 2018-10-24 17:27:31Z gronemeier
[3405]33! bugfix: nx, ny are required in non-parallel case
34!
35! 3337 2018-10-12 15:17:09Z kanani
[3337]36! (from branch resler)
37! Add Biometeorology
38!
39! 3294 2018-10-01 02:37:10Z raasch
[3294]40! changes concerning modularization of ocean option
41!
42! 3274 2018-09-24 15:42:55Z knoop
[3274]43! Modularization of all bulk cloud physics code components
44!
45! 3241 2018-09-12 15:02:00Z raasch
[3241]46! unused variables and format statements removed
47!
48! 3049 2018-05-29 13:52:36Z Giersch
[3049]49! Error messages revised
50!
51! 3045 2018-05-28 07:55:41Z Giersch
[3045]52! Error message revised
53!
54! 3014 2018-05-09 08:42:38Z maronga
[3014]55! Added nzb_do and nzt_do for some modules for 3d output
56!
57! 3004 2018-04-27 12:33:25Z Giersch
[3004]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
[2967]62! bugfix: missing parallel cpp-directives added
63!
64! 2817 2018-02-19 16:32:21Z knoop
[2817]65! Preliminary gust module interface implemented
66!
67! 2766 2018-01-22 17:17:47Z kanani
[2766]68! Removed preprocessor directive __chem
69!
70! 2756 2018-01-16 18:11:14Z suehring
[2756]71! Fill values for 3D output of chemical species introduced.
72!
73! 2746 2018-01-15 12:06:04Z suehring
[2746]74! Move flag plant canopy to modules
75!
76! 2718 2018-01-02 08:49:38Z maronga
[2716]77! Corrected "Former revisions" section
78!
79! 2696 2017-12-14 17:12:51Z kanani
80! Change in file header (GPL part)
[2696]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
[2512]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
[2292]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
[1552]96!
[2233]97! 2232 2017-05-30 17:47:52Z suehring
98! Adjustments to new topography concept
99!
[2210]100! 2209 2017-04-19 09:34:46Z kanani
101! Added plant canopy model output
102!
[2032]103! 2031 2016-10-21 15:11:58Z knoop
104! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
105!
[2012]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!
[2008]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!
[2001]115! 2000 2016-08-20 18:09:15Z knoop
116! Forced header and separation lines into 80 columns
117!
[1981]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!
[1977]122! 1976 2016-07-27 13:28:04Z maronga
123! Output of radiation quantities is now done directly in the respective module
124!
[1973]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!
[1961]129! 1960 2016-07-12 16:34:24Z suehring
130! Scalar surface flux added
131!
[1851]132! 1849 2016-04-08 11:33:18Z hoffmann
133! prr moved to arrays_3d
134!
[1823]135! 1822 2016-04-07 07:49:42Z hoffmann
136! prr vertical dimensions set to nzb_do to nzt_do. Unused variables deleted.
137!
[1809]138! 1808 2016-04-05 19:44:00Z raasch
139! test output removed
140!
[1784]141! 1783 2016-03-06 18:36:17Z raasch
142! name change of netcdf routines and module + related changes
143!
[1746]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!
[1692]147! 1691 2015-10-26 16:17:44Z maronga
148! Added output of radiative heating rates for RRTMG
149!
[1683]150! 1682 2015-10-07 23:56:08Z knoop
151! Code annotations made doxygen readable
152!
[1586]153! 1585 2015-04-30 07:05:52Z maronga
154! Added support for RRTMG
155!
[1552]156! 1551 2015-03-03 14:18:16Z maronga
[1551]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.
[1329]161!
[1360]162! 1359 2014-04-11 17:15:14Z hoffmann
163! New particle structure integrated.
164!
[1354]165! 1353 2014-04-08 15:21:23Z heinze
166! REAL constants provided with KIND-attribute
167!
[1329]168! 1327 2014-03-21 11:00:16Z raasch
169! parts concerning avs output removed,
170! -netcdf output queries
171!
[1321]172! 1320 2014-03-20 08:40:49Z raasch
[1320]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
[674]180!
[1319]181! 1318 2014-03-17 13:35:16Z raasch
182! barrier argument removed from cpu_log,
183! module interfaces removed
184!
[1309]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!
[1245]190! 1244 2013-10-31 08:16:56Z raasch
191! Bugfix for index bounds in case of 3d-parallel output
192!
[1116]193! 1115 2013-03-26 18:16:16Z hoffmann
194! ql is calculated by calc_liquid_water_content
195!
[1107]196! 1106 2013-03-04 05:31:38Z raasch
197! array_kind renamed precision_kind
198!
[1077]199! 1076 2012-12-05 08:30:18Z hoffmann
200! Bugfix in output of ql
201!
[1054]202! 1053 2012-11-13 17:11:03Z hoffmann
203! +nr, qr, prr, qc and averaged quantities
204!
[1037]205! 1036 2012-10-22 13:43:42Z raasch
206! code put under GPL (PALM 3.9)
207!
[1035]208! 1031 2012-10-19 14:35:30Z raasch
209! netCDF4 without parallel file support implemented
210!
[1008]211! 1007 2012-09-19 14:30:36Z franke
212! Bugfix: missing calculation of ql_vp added
213!
[1]214! Revision 1.1  1997/09/03 06:29:36  raasch
215! Initial revision
216!
217!
218! Description:
219! ------------
[1682]220!> Output of the 3D-arrays in netCDF and/or AVS format.
[1]221!------------------------------------------------------------------------------!
[1682]222 SUBROUTINE data_output_3d( av )
223 
[1]224
[1320]225    USE arrays_3d,                                                             &
[3294]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
[3274]228
[1]229    USE averaging
[3274]230
231    USE basic_constants_and_equations_mod,                                     &
232        ONLY:  lv_d_cp
233
[3448]234    USE biometeorology_mod,                                                    &
235        ONLY:  biom_data_output_3d
236
[3274]237    USE bulk_cloud_model_mod,                                                  &
238        ONLY:  bulk_cloud_model, bcm_data_output_3d
239
[2696]240    USE chemistry_model_mod,                                                   &
241        ONLY:  chem_data_output_3d
242
[1320]243    USE control_parameters,                                                    &
[3448]244        ONLY:  air_chemistry, biometeorology, do3d, do3d_no, do3d_time_count,  &
[2696]245               io_blocks, io_group, land_surface, message_string,              &
[3294]246               ntdim_3d, nz_do3d,  ocean_mode, plant_canopy,                   &
[2232]247               psolver, simulated_time, time_since_reference_point,            &
248               urban_surface, varnamelength
[3274]249
[1320]250    USE cpulog,                                                                &
251        ONLY:  log_point, cpu_log
[2817]252
253    USE gust_mod,                                                              &
254        ONLY: gust_data_output_3d, gust_module_enabled
[3274]255
[3405]256#if defined( __parallel )
[1320]257    USE indices,                                                               &
[3241]258        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,     &
259               wall_flags_0
[3405]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
[3274]265
[1320]266    USE kinds
[3274]267
[1551]268    USE land_surface_model_mod,                                                &
[2232]269        ONLY: lsm_data_output_3d, nzb_soil, nzt_soil
[1551]270
[1783]271#if defined( __netcdf )
272    USE NETCDF
273#endif
274
275    USE netcdf_interface,                                                      &
[2696]276        ONLY:  fill_value, id_set_3d, id_var_do3d, id_var_time_3d, nc_stat,    &
[1783]277               netcdf_data_format, netcdf_handle_error
[1320]278       
[3294]279    USE ocean_mod,                                                             &
280        ONLY:  ocean_data_output_3d
281
[1320]282    USE particle_attributes,                                                   &
[1359]283        ONLY:  grid_particles, number_of_particles, particles,                 &
284               particle_advection_start, prt_count
[1320]285       
[1]286    USE pegrid
287
[2209]288    USE plant_canopy_model_mod,                                                &
[2746]289        ONLY:  pcm_data_output_3d
[2209]290       
[1585]291    USE radiation_model_mod,                                                   &
[2696]292        ONLY:  nzub, nzut, radiation, radiation_data_output_3d
[1585]293
[2696]294    USE turbulence_closure_mod,                                                &
295        ONLY:  tcm_data_output_3d
296
[2007]297    USE urban_surface_mod,                                                     &
[2696]298        ONLY:  usm_data_output_3d
[1585]299
[2007]300
[1]301    IMPLICIT NONE
302
[1682]303    INTEGER(iwp) ::  av        !<
[2696]304    INTEGER(iwp) ::  flag_nr   !< number of masking flag
[1682]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
[1]312
[1682]313    LOGICAL      ::  found     !<
314    LOGICAL      ::  resorted  !<
[1]315
[1682]316    REAL(wp)     ::  mean_r    !<
317    REAL(wp)     ::  s_r2      !<
318    REAL(wp)     ::  s_r3      !<
[1]319
[1682]320    REAL(sp), DIMENSION(:,:,:), ALLOCATABLE ::  local_pf  !<
[1]321
[1682]322    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !<
[1]323
[2011]324    CHARACTER (LEN=varnamelength) ::  trimvar  !< TRIM of output-variable string
[2007]325
[1]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.
[2512]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.
[1031]337!-- For netCDF4/HDF5 output, data is written in parallel into one file.
[1327]338    IF ( netcdf_data_format < 5 )  THEN
[2967]339#if defined( __parallel )
[1327]340       CALL check_open( 30 )
[2967]341#endif
[1327]342       IF ( myid == 0 )  CALL check_open( 106+av*10 )
[493]343    ELSE
[1327]344       CALL check_open( 106+av*10 )
[493]345    ENDIF
[1]346
347!
[1745]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=',  &
[3046]354                                      simulated_time, '&because the maximum ', & 
[1745]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!
[1031]364!-- Update the netCDF time axis
[1308]365!-- In case of parallel output, this is only done by PE0 to increase the
366!-- performance.
[1]367#if defined( __netcdf )
[1308]368    do3d_time_count(av) = do3d_time_count(av) + 1
369    IF ( myid == 0 )  THEN
[1327]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 /) )
[1783]374       CALL netcdf_handle_error( 'data_output_3d', 376 )
[1]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) /= ' ' )
[2007]383
[1]384!
[2007]385!--    Temporary solution to account for data output within the new urban
386!--    surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar ).
[1]387!--    Store the array chosen on the temporary array.
[2007]388       trimvar = TRIM( do3d(av,if) )
[2011]389       IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
[2007]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
[1551]399!
[1980]400!--    Set flag to steer output of radiation, land-surface, or user-defined
401!--    quantities
402       found = .FALSE.
403!
[1551]404!--    Allocate a temporary array with the desired output dimensions.
[2512]405       ALLOCATE( local_pf(nxl:nxr,nys:nyn,nzb_do:nzt_do) )
[2696]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
[1551]412
[2007]413       SELECT CASE ( trimvar )
[1]414
415          CASE ( 'e' )
416             IF ( av == 0 )  THEN
417                to_be_resorted => e
418             ELSE
[3004]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
[1]423                to_be_resorted => e_av
424             ENDIF
425
[3421]426          CASE ( 'thetal' )
[771]427             IF ( av == 0 )  THEN
428                to_be_resorted => pt
429             ELSE
[3004]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
[771]434                to_be_resorted => lpt_av
435             ENDIF
436
[1]437          CASE ( 'p' )
438             IF ( av == 0 )  THEN
[727]439                IF ( psolver /= 'sor' )  CALL exchange_horiz( p, nbgp )
[1]440                to_be_resorted => p
441             ELSE
[3004]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
[727]446                IF ( psolver /= 'sor' )  CALL exchange_horiz( p_av, nbgp )
[1]447                to_be_resorted => p_av
448             ENDIF
449
450          CASE ( 'pc' )  ! particle concentration (requires ghostpoint exchange)
451             IF ( av == 0 )  THEN
[1359]452                IF ( simulated_time >= particle_advection_start )  THEN
453                   tend = prt_count
454                ELSE
455                   tend = 0.0_wp
456                ENDIF
[2512]457                DO  i = nxl, nxr
458                   DO  j = nys, nyn
[1551]459                      DO  k = nzb_do, nzt_do
[1]460                         local_pf(i,j,k) = tend(k,j,i)
461                      ENDDO
462                   ENDDO
463                ENDDO
464                resorted = .TRUE.
465             ELSE
[3004]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
[1]470                to_be_resorted => pc_av
471             ENDIF
472
[1359]473          CASE ( 'pr' )  ! mean particle radius (effective radius)
[1]474             IF ( av == 0 )  THEN
[1359]475                IF ( simulated_time >= particle_advection_start )  THEN
476                   DO  i = nxl, nxr
477                      DO  j = nys, nyn
[1551]478                         DO  k = nzb_do, nzt_do
[1359]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
[1]498                         ENDDO
499                      ENDDO
500                   ENDDO
[1359]501                ELSE
502                   tend = 0.0_wp
503                ENDIF
[2512]504                DO  i = nxl, nxr
505                   DO  j = nys, nyn
[1551]506                      DO  k = nzb_do, nzt_do
[1]507                         local_pf(i,j,k) = tend(k,j,i)
508                      ENDDO
509                   ENDDO
510                ENDDO
511                resorted = .TRUE.
512             ELSE
[3004]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
[1]517                to_be_resorted => pr_av
518             ENDIF
519
[3421]520          CASE ( 'theta' )
[1]521             IF ( av == 0 )  THEN
[3274]522                IF ( .NOT. bulk_cloud_model ) THEN
[1]523                   to_be_resorted => pt
524                ELSE
[2512]525                   DO  i = nxl, nxr
526                      DO  j = nys, nyn
[1551]527                         DO  k = nzb_do, nzt_do
[3274]528                            local_pf(i,j,k) = pt(k,j,i) + lv_d_cp *            &
529                                                          d_exner(k) *         &
[1]530                                                          ql(k,j,i)
531                         ENDDO
532                      ENDDO
533                   ENDDO
534                   resorted = .TRUE.
535                ENDIF
536             ELSE
[3004]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
[1]541                to_be_resorted => pt_av
542             ENDIF
543
544          CASE ( 'q' )
545             IF ( av == 0 )  THEN
546                to_be_resorted => q
547             ELSE
[3004]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
[1]552                to_be_resorted => q_av
553             ENDIF
[691]554
[1053]555          CASE ( 'ql' )
556             IF ( av == 0 )  THEN
[1115]557                to_be_resorted => ql
[1053]558             ELSE
[3004]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
[1115]563                to_be_resorted => ql_av
[1053]564             ENDIF
565
[1]566          CASE ( 'ql_c' )
567             IF ( av == 0 )  THEN
568                to_be_resorted => ql_c
569             ELSE
[3004]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
[1]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
[3004]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
[1]585                to_be_resorted => ql_v_av
586             ENDIF
587
588          CASE ( 'ql_vp' )
589             IF ( av == 0 )  THEN
[1359]590                IF ( simulated_time >= particle_advection_start )  THEN
591                   DO  i = nxl, nxr
592                      DO  j = nys, nyn
[1551]593                         DO  k = nzb_do, nzt_do
[1359]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
[1007]604                         ENDDO
605                      ENDDO
606                   ENDDO
[1359]607                ELSE
608                   tend = 0.0_wp
609                ENDIF
[2512]610                DO  i = nxl, nxr
611                   DO  j = nys, nyn
[1551]612                      DO  k = nzb_do, nzt_do
[1007]613                         local_pf(i,j,k) = tend(k,j,i)
614                      ENDDO
615                   ENDDO
616                ENDDO
617                resorted = .TRUE.
[1]618             ELSE
[3004]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
[1]623                to_be_resorted => ql_vp_av
624             ENDIF
625
626          CASE ( 'qv' )
627             IF ( av == 0 )  THEN
[2512]628                DO  i = nxl, nxr
629                   DO  j = nys, nyn
[1551]630                      DO  k = nzb_do, nzt_do
[1]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
[3004]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
[1]641                to_be_resorted => qv_av
642             ENDIF
643
644          CASE ( 's' )
645             IF ( av == 0 )  THEN
[1960]646                to_be_resorted => s
[1]647             ELSE
[3004]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
[355]652                to_be_resorted => s_av
[1]653             ENDIF
[691]654
[1]655          CASE ( 'u' )
[2696]656             flag_nr = 1
[1]657             IF ( av == 0 )  THEN
658                to_be_resorted => u
659             ELSE
[3004]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
[1]664                to_be_resorted => u_av
665             ENDIF
666
667          CASE ( 'v' )
[2696]668             flag_nr = 2
[1]669             IF ( av == 0 )  THEN
670                to_be_resorted => v
671             ELSE
[3004]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
[1]676                to_be_resorted => v_av
677             ENDIF
678
[3421]679          CASE ( 'thetav' )
[1]680             IF ( av == 0 )  THEN
681                to_be_resorted => vpt
682             ELSE
[3004]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
[1]687                to_be_resorted => vpt_av
688             ENDIF
689
690          CASE ( 'w' )
[2696]691             flag_nr = 3
[1]692             IF ( av == 0 )  THEN
693                to_be_resorted => w
694             ELSE
[3004]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
[1]699                to_be_resorted => w_av
700             ENDIF
[2007]701!             
702!--       Block of urban surface model outputs   
703          CASE ( 'usm_output' )
[3274]704             CALL usm_data_output_3d( av, do3d(av,if), found, local_pf,        &
[2007]705                                         nzb_do, nzt_do )
[1]706
707          CASE DEFAULT
[1972]708
[1]709!
[3294]710!--          Quantities of other modules
[3274]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
[3294]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
[3274]729             IF ( .NOT. found  .AND.  land_surface )  THEN
[1972]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 )
[2512]736                ALLOCATE( local_pf(nxl:nxr,nys:nyn,nzb_do:nzt_do) )
[2696]737                local_pf = fill_value
[1972]738
739                CALL lsm_data_output_3d( av, do3d(av,if), found, local_pf )
740                resorted = .TRUE.
[1976]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 )
[2512]749                   ALLOCATE( local_pf(nxl:nxr,nys:nyn,nzb_do:nzt_do) )                 
[1976]750                ENDIF
751
[1972]752             ENDIF
753
[3294]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 )
[1976]757                resorted = .TRUE.
758             ENDIF
759
[3294]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 )
[2817]763                resorted = .TRUE.
764             ENDIF
765
[3294]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 )
[2696]769                resorted = .TRUE.
770             ENDIF
771
[3294]772             IF ( .NOT. found )  THEN
773                CALL tcm_data_output_3d( av, do3d(av,if), found, local_pf,     &
774                                         nzb_do, nzt_do )
[2209]775                resorted = .TRUE.
776             ENDIF
777
[3448]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
[2209]783!
[3294]784!--          User defined quantities
[1972]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
[1]790
[254]791             IF ( .NOT. found )  THEN
[1320]792                message_string =  'no output available for: ' //               &
[274]793                                  TRIM( do3d(av,if) )
[254]794                CALL message( 'data_output_3d', 'PA0182', 0, 0, 0, 6, 0 )
[1]795             ENDIF
796
797       END SELECT
798
799!
800!--    Resort the array to be output, if not done above
801       IF ( .NOT. resorted )  THEN
[2512]802          DO  i = nxl, nxr
803             DO  j = nys, nyn
[1551]804                DO  k = nzb_do, nzt_do
[2696]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 ) )
[1]809                ENDDO
810             ENDDO
811          ENDDO
812       ENDIF
813
814!
[1327]815!--    Output of the 3D-array
816#if defined( __parallel )
817       IF ( netcdf_data_format < 5 )  THEN
[1]818!
[1327]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
[1]825          ENDIF
[1327]826          DO  i = 0, io_blocks-1
827             IF ( i == io_group )  THEN
[2512]828                WRITE ( 30 )  nxl, nxr, nys, nyn, nzb_do, nzt_do
[1551]829                WRITE ( 30 )  local_pf(:,:,nzb_do:nzt_do)
[1327]830             ENDIF
[1972]831
[1327]832             CALL MPI_BARRIER( comm2d, ierr )
[1972]833
[1327]834          ENDDO
[559]835
[1327]836       ELSE
[646]837#if defined( __netcdf )
[493]838!
[1327]839!--       Parallel output in netCDF4/HDF5 format.
[2512]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
[1327]856             nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),  &
[1551]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 /) )
[2512]860!          ENDIF
[1783]861          CALL netcdf_handle_error( 'data_output_3d', 386 )
[646]862#endif
[1327]863       ENDIF
[1]864#else
865#if defined( __netcdf )
[1327]866       nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_do3d(av,if),        &
[1551]867                         local_pf(nxl:nxr+1,nys:nyn+1,nzb_do:nzt_do),        &
[1327]868                         start = (/ 1, 1, 1, do3d_time_count(av) /),     &
[2512]869                         count = (/ nx+1, ny+1, nzt_do-nzb_do+1, 1 /) )
[1783]870       CALL netcdf_handle_error( 'data_output_3d', 446 )
[1]871#endif
872#endif
873
874       if = if + 1
875
876!
[1551]877!--    Deallocate temporary array
878       DEALLOCATE ( local_pf )
[1]879
[1551]880    ENDDO
[1]881
[1318]882    CALL cpu_log( log_point(14), 'data_output_3d', 'stop' )
[1]883
884 END SUBROUTINE data_output_3d
Note: See TracBrowser for help on using the repository browser.