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

Last change on this file since 3582 was 3582, checked in by suehring, 5 years ago

Merge branch salsa with trunk

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