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

Last change on this file since 3058 was 3049, checked in by Giersch, 7 years ago

Revision history corrected

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