source: palm/trunk/SOURCE/check_parameters.f90 @ 3049

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

Revision history corrected

  • Property svn:keywords set to Id
  • Property svn:mergeinfo set to False
    /palm/branches/forwind/SOURCE/check_parameters.f901564-1913
    /palm/branches/palm4u/SOURCE/check_parameters.f902540-2692
    /palm/branches/rans/SOURCE/check_parameters.f902078-2917
File size: 181.3 KB
Line 
1!> @file check_parameters.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: check_parameters.f90 3049 2018-05-29 13:52:36Z Giersch $
27! add variable description
28!
29! 3046 2018-05-29 08:02:15Z Giersch
30! Error messages revised
31!
32! 3045 2018-05-28 07:55:41Z Giersch
33! Error messages revised
34!
35! 3035 2018-05-24 09:35:20Z schwenkel
36! Add option to initialize warm air bubble close to surface
37!
38! 3034 2018-05-24 08:41:20Z raasch
39! bugfix: check that initializing_actions has been set
40!
41! 2980 2018-04-17 15:19:27Z suehring
42! Further improvement for spinup checks.
43!
44! 2974 2018-04-16 12:59:52Z gronemeier
45! Bugfix: check if dt_data_output_av is zero in case of parallel NetCDF output
46!
47! 2970 2018-04-13 15:09:23Z suehring
48! Bugfix in old large-scale forcing mode
49!
50! 2964 2018-04-12 16:04:03Z Giersch
51! Calculation of fixed number of output time levels for parallel netcdf output
52! has been revised (based on calculations in netcdf_interface_mod)
53!
54! 2938 2018-03-27 15:52:42Z suehring
55! - Revise start and end indices for imposing random disturbances in case of
56!   nesting or large-scale forcing.
57! - Remove check for inifor initialization in case of nested runs.
58! - Adapt call to check for synthetic turbulence geneartor settings.
59!
60! 2936 2018-03-27 14:49:27Z suehring
61! Check spinup in case of nested runs, spinup time and timestep must be
62! identical to assure synchronuous simulations.
63!
64! 2932 2018-03-26 09:39:22Z maronga
65! renamed particles_par to particle_parameters
66!
67! 2918 2018-03-21 15:52:14Z gronemeier
68! Add check for 1D model
69!
70! 2883 2018-03-14 08:29:10Z Giersch
71! dt_dopr_listing is not set to the default value zero anymore
72!
73! 2851 2018-03-05 14:39:31Z maronga
74! Bugfix: calculation of output time levels in case of restart runs (parallel
75! NetCDF)
76!
77! 2836 2018-02-26 13:40:05Z Giersch
78! dt_dopr_listing is set to the default value zero
79!
80! 2817 2018-02-19 16:32:21Z knoop
81! Preliminary gust module interface implemented
82!
83! 2798 2018-02-09 17:16:39Z suehring
84! Consider also default-type surfaces for surface temperature output.
85!
86! 2797 2018-02-08 13:24:35Z suehring
87! Enable output of ground-heat flux also at urban surfaces.
88!
89! 2776 2018-01-31 10:44:42Z Giersch
90! Variable synthetic_turbulence_generator has been abbreviated
91!
92! 2773 2018-01-30 14:12:54Z suehring
93! Check for consistent initialization in nesting mode added.
94!
95! 2766 2018-01-22 17:17:47Z kanani
96! Removed preprocessor directive __chem
97!
98! 2765 2018-01-22 11:34:58Z maronga
99! Renamed simulation_time_since_reference to
100! time_to_be_simulated_from_reference_point
101!
102! 2746 2018-01-15 12:06:04Z suehring
103! Move flag plant canopy to modules
104!
105! 2743 2018-01-12 16:03:39Z suehring
106! In case of natural- and urban-type surfaces output surfaces fluxes in W/m2.
107!
108! 2742 2018-01-12 14:59:47Z suehring
109! Enable output of surface temperature
110!
111! 2735 2018-01-11 12:01:27Z suehring
112! output of r_a moved from land-surface to consider also urban-type surfaces
113!
114! 2718 2018-01-02 08:49:38Z maronga
115! Corrected "Former revisions" section
116!
117! 2696 2017-12-14 17:12:51Z kanani
118! Change in file header (GPL part)
119! Implementation of uv exposure model (FK)
120! + new possible value for dissipation_1d
121! Added checks for turbulence_closure_mod (TG)
122! Implementation of chemistry module (FK)
123!
124! 2689 2017-12-12 17:46:55Z Giersch
125! Bugfix in if query
126!
127! 2688 2017-12-12 17:27:04Z Giersch
128! Check if humidity is set to TRUE in the _p3d file for coupled runs
129!
130! 2669 2017-12-06 16:03:27Z raasch
131! mrun-string replaced by palmrun
132!
133! 2628 2017-11-20 12:40:38Z schwenkel
134! Enabled particle advection with grid stretching -> Removed parameter check
135!
136! 2575 2017-10-24 09:57:58Z maronga
137! Renamed phi --> latitude
138!
139! 2564 2017-10-19 15:56:56Z Giersch
140! Variable wind_turbine was added to control_parameters.
141!
142! 2550 2017-10-16 17:12:01Z boeske
143! Added checks for complex terrain simulations
144!
145! 2513 2017-10-04 09:24:39Z kanani
146! Bugfix for some dopr(_initial)_index values and units connected to
147! passive-scalar output
148!
149! 2508 2017-10-02 08:57:09Z suehring
150! Bugfix, change default value of vertical_gradient level in order to consider
151! also ocean runs
152!
153! 2422 2017-09-08 08:25:41Z raasch
154! error message in case of missing "restart" file activation string
155!
156! 2375 2017-08-29 14:10:28Z schwenkel
157! Added aerosol for bulk microphysics
158!
159! 2365 2017-08-21 14:59:59Z kanani
160! Vertical grid nesting implemented: Check coupling mode. Generate file header
161! (SadiqHuq)
162!
163! 2354 2017-08-17 10:49:36Z schwenkel
164! Bugfix correlated to lsm_check_data_output_pr.
165! If-statement for following checks is essential, otherwise units for lsm output
166! are set to 'illegal' and palm will be aborted.
167!
168! 2348 2017-08-10 10:40:10Z kanani
169! New: Check for simultaneous use of geostrophic wind and u_profile/v_profile
170!
171! 2345 2017-08-09 11:50:30Z Giersch
172! Remove error message PA0156 and the conserve_volume_flow_mode option
173! inflow_profile
174!
175! 2339 2017-08-07 13:55:26Z raasch
176! corrected timestamp in header
177!
178! 2338 2017-08-07 12:15:38Z gronemeier
179! Modularize 1D model
180!
181! 2329 2017-08-03 14:24:56Z knoop
182! Bugfix: index corrected for rho_air and rho_air_zw output
183!
184! 2320 2017-07-21 12:47:43Z suehring
185! Modularize large-scale forcing and nudging
186!
187! 2312 2017-07-14 20:26:51Z hoffmann
188! PA0349 and PA0420 removed.
189!
190! 2300 2017-06-29 13:31:14Z raasch
191! host-specific settings and checks removed
192!
193! 2292 2017-06-20 09:51:42Z schwenkel
194! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
195! includes two more prognostic equations for cloud drop concentration (nc)
196! and cloud water content (qc).
197!
198! 2274 2017-06-09 13:27:48Z Giersch
199! Changed error messages
200!
201! 2271 2017-06-09 12:34:55Z sward
202! roughness-length check altered
203! Error messages fixed
204!
205! 2270 2017-06-09 12:18:47Z maronga
206! Revised numbering (removed 2 timeseries)
207!
208! 2259 2017-06-08 09:09:11Z gronemeier
209! Implemented synthetic turbulence generator
210!
211! 2251 2017-06-06 15:10:46Z Giersch
212!
213! 2250 2017-06-06 15:07:12Z Giersch
214! Doxygen comment added
215!
216! 2248 2017-06-06 13:52:54Z sward
217! Error message fixed
218!
219! 2209 2017-04-19 09:34:46Z kanani
220! Check for plant canopy model output
221!
222! 2200 2017-04-11 11:37:51Z suehring
223! monotonic_adjustment removed
224!
225! 2178 2017-03-17 11:07:39Z hellstea
226! Index limits for perturbations are now set also in case of nested
227! boundary conditions
228!
229! 2169 2017-03-06 18:16:35Z suehring
230! Bugfix, move setting for topography grid convention to init_grid
231!
232! 2118 2017-01-17 16:38:49Z raasch
233! OpenACC related parts of code removed
234!
235! 2088 2016-12-19 16:30:25Z suehring
236! Bugfix in initial salinity profile
237!
238! 2084 2016-12-09 15:59:42Z knoop
239! Removed anelastic + multigrid error message
240!
241! 2052 2016-11-08 15:14:59Z gronemeier
242! Bugfix: remove setting of default value for recycling_width
243!
244! 2050 2016-11-08 15:00:55Z gronemeier
245! Implement turbulent outflow condition
246!
247! 2044 2016-11-02 16:44:25Z knoop
248! Added error code for anelastic approximation
249!
250! 2042 2016-11-02 13:47:31Z suehring
251! Additional checks for wall_heatflux, wall_humidityflux and wall_scalarflux.
252! Bugfix, check for constant_scalarflux.
253!
254! 2040 2016-10-26 16:58:09Z gronemeier
255! Removed check for statistic_regions > 9.
256!
257! 2037 2016-10-26 11:15:40Z knoop
258! Anelastic approximation implemented
259!
260! 2031 2016-10-21 15:11:58Z knoop
261! renamed variable rho to rho_ocean
262!
263! 2026 2016-10-18 10:27:02Z suehring
264! Bugfix, enable output of s*2
265!
266! 2024 2016-10-12 16:42:37Z kanani
267! Added missing CASE, error message and unit for ssws*,
268! increased number of possible output quantities to 500.
269!
270! 2011 2016-09-19 17:29:57Z kanani
271! Flag urban_surface is now defined in module control_parameters,
272! changed prefix for urban surface model output to "usm_",
273! added flag lsf_exception (inipar-Namelist parameter) to allow explicit
274! enabling of large scale forcing together with buildings on flat terrain,
275! introduced control parameter varnamelength for LEN of var.
276!
277! 2007 2016-08-24 15:47:17Z kanani
278! Added checks for the urban surface model,
279! increased counter in DO WHILE loop over data_output (for urban surface output)
280!
281! 2000 2016-08-20 18:09:15Z knoop
282! Forced header and separation lines into 80 columns
283!
284! 1994 2016-08-15 09:52:21Z suehring
285! Add missing check for cloud_physics and cloud_droplets
286!
287! 1992 2016-08-12 15:14:59Z suehring
288! New checks for top_scalarflux
289! Bugfixes concerning data output of profiles in case of passive_scalar
290!
291! 1984 2016-08-01 14:48:05Z suehring
292! Bugfix: checking for bottom and top boundary condition for humidity and scalars
293!
294! 1972 2016-07-26 07:52:02Z maronga
295! Removed check of lai* (done in land surface module)
296!
297! 1970 2016-07-18 14:27:12Z suehring
298! Bugfix, check for ibc_q_b and constant_waterflux in case of land-surface scheme
299!
300! 1962 2016-07-13 09:16:44Z suehring
301! typo removed
302!
303! 1960 2016-07-12 16:34:24Z suehring
304! Separate checks and calculations for humidity and passive scalar. Therefore,
305! a subroutine to calculate vertical gradients is introduced, in order  to reduce
306! redundance.
307! Additional check - large-scale forcing is not implemented for passive scalars
308! so far.
309!
310! 1931 2016-06-10 12:06:59Z suehring
311! Rename multigrid into multigrid_noopt and multigrid_fast into multigrid
312!
313! 1929 2016-06-09 16:25:25Z suehring
314! Bugfix in check for use_upstream_for_tke
315!
316! 1918 2016-05-27 14:35:57Z raasch
317! setting of a fixed reference state ('single_value') for ocean runs removed
318!
319! 1914 2016-05-26 14:44:07Z witha
320! Added checks for the wind turbine model
321!
322! 1841 2016-04-07 19:14:06Z raasch
323! redundant location message removed
324!
325! 1833 2016-04-07 14:23:03Z raasch
326! check of spectra quantities moved to spectra_mod
327!
328! 1829 2016-04-07 12:16:29Z maronga
329! Bugfix: output of user defined data required reset of the string 'unit'
330!
331! 1826 2016-04-07 12:01:39Z maronga
332! Moved checks for radiation model to the respective module.
333! Moved checks for plant canopy model to the respective module.
334! Bugfix for check of too large roughness_length
335!
336! 1824 2016-04-07 09:55:29Z gronemeier
337! Check if roughness_length < dz/2. Moved location_message(finished) to the end.
338!
339! 1822 2016-04-07 07:49:42Z hoffmann
340! PALM collision kernel deleted. Collision algorithms are checked for correct
341! spelling.
342!
343! Tails removed.
344! !
345! Checks for use_sgs_for_particles adopted for the use of droplets with
346! use_sgs_for_particles adopted.
347!
348! Unused variables removed.
349!
350! 1817 2016-04-06 15:44:20Z maronga
351! Moved checks for land_surface model to the respective module
352!
353! 1806 2016-04-05 18:55:35Z gronemeier
354! Check for recycling_yshift
355!
356! 1804 2016-04-05 16:30:18Z maronga
357! Removed code for parameter file check (__check)
358!
359! 1795 2016-03-18 15:00:53Z raasch
360! Bugfix: setting of initial scalar profile in ocean
361!
362! 1788 2016-03-10 11:01:04Z maronga
363! Added check for use of most_method = 'lookup' in combination with water
364! surface presribed in the land surface model. Added output of z0q.
365! Syntax layout improved.
366!
367! 1786 2016-03-08 05:49:27Z raasch
368! cpp-direktives for spectra removed, check of spectra level removed
369!
370! 1783 2016-03-06 18:36:17Z raasch
371! netcdf variables and module name changed,
372! check of netcdf precision removed (is done in the netcdf module)
373!
374! 1764 2016-02-28 12:45:19Z raasch
375! output of nest id in run description header,
376! bugfix: check of validity of lateral boundary conditions moved to parin
377!
378! 1762 2016-02-25 12:31:13Z hellstea
379! Introduction of nested domain feature
380!
381! 1745 2016-02-05 13:06:51Z gronemeier
382! Bugfix: check data output intervals to be /= 0.0 in case of parallel NetCDF4
383!
384! 1701 2015-11-02 07:43:04Z maronga
385! Bugfix: definition of rad_net timeseries was missing
386!
387! 1691 2015-10-26 16:17:44Z maronga
388! Added output of Obukhov length (ol) and radiative heating rates for RRTMG.
389! Added checks for use of radiation / lsm with topography.
390!
391! 1682 2015-10-07 23:56:08Z knoop
392! Code annotations made doxygen readable
393!
394! 1606 2015-06-29 10:43:37Z maronga
395! Added check for use of RRTMG without netCDF.
396!
397! 1587 2015-05-04 14:19:01Z maronga
398! Added check for set of albedos when using RRTMG
399!
400! 1585 2015-04-30 07:05:52Z maronga
401! Added support for RRTMG
402!
403! 1575 2015-03-27 09:56:27Z raasch
404! multigrid_fast added as allowed pressure solver
405!
406! 1573 2015-03-23 17:53:37Z suehring
407! Bugfix: check for advection schemes in case of non-cyclic boundary conditions
408!
409! 1557 2015-03-05 16:43:04Z suehring
410! Added checks for monotonic limiter
411!
412! 1555 2015-03-04 17:44:27Z maronga
413! Added output of r_a and r_s. Renumbering of LSM PA-messages.
414!
415! 1553 2015-03-03 17:33:54Z maronga
416! Removed check for missing soil temperature profile as default values were added.
417!
418! 1551 2015-03-03 14:18:16Z maronga
419! Added various checks for land surface and radiation model. In the course of this
420! action, the length of the variable var has to be increased
421!
422! 1504 2014-12-04 09:23:49Z maronga
423! Bugfix: check for large_scale forcing got lost.
424!
425! 1500 2014-12-03 17:42:41Z maronga
426! Boundary conditions changed to dirichlet for land surface model
427!
428! 1496 2014-12-02 17:25:50Z maronga
429! Added checks for the land surface model
430!
431! 1484 2014-10-21 10:53:05Z kanani
432! Changes due to new module structure of the plant canopy model:
433!   module plant_canopy_model_mod added,
434!   checks regarding usage of new method for leaf area density profile
435!   construction added,
436!   lad-profile construction moved to new subroutine init_plant_canopy within
437!   the module plant_canopy_model_mod,
438!   drag_coefficient renamed to canopy_drag_coeff.
439! Missing KIND-attribute for REAL constant added
440!
441! 1455 2014-08-29 10:47:47Z heinze
442! empty time records in volume, cross-section and masked data output prevented
443! in case of non-parallel netcdf-output in restart runs
444!
445! 1429 2014-07-15 12:53:45Z knoop
446! run_description_header exended to provide ensemble_member_nr if specified
447!
448! 1425 2014-07-05 10:57:53Z knoop
449! bugfix: perturbation domain modified for parallel random number generator
450!
451! 1402 2014-05-09 14:25:13Z raasch
452! location messages modified
453!
454! 1400 2014-05-09 14:03:54Z knoop
455! Check random generator extended by option random-parallel
456!
457! 1384 2014-05-02 14:31:06Z raasch
458! location messages added
459!
460! 1365 2014-04-22 15:03:56Z boeske
461! Usage of large scale forcing for pt and q enabled:
462! output for profiles of large scale advection (td_lsa_lpt, td_lsa_q),
463! large scale subsidence (td_sub_lpt, td_sub_q)
464! and nudging tendencies (td_nud_lpt, td_nud_q, td_nud_u and td_nud_v) added,
465! check if use_subsidence_tendencies is used correctly
466!
467! 1361 2014-04-16 15:17:48Z hoffmann
468! PA0363 removed
469! PA0362 changed
470!
471! 1359 2014-04-11 17:15:14Z hoffmann
472! Do not allow the execution of PALM with use_particle_tails, since particle
473! tails are currently not supported by our new particle structure.
474!
475! PA0084 not necessary for new particle structure
476!
477! 1353 2014-04-08 15:21:23Z heinze
478! REAL constants provided with KIND-attribute
479!
480! 1330 2014-03-24 17:29:32Z suehring
481! In case of SGS-particle velocity advection of TKE is also allowed with
482! dissipative 5th-order scheme.
483!
484! 1327 2014-03-21 11:00:16Z raasch
485! "baroclinicity" renamed "baroclinity", "ocean version" replaced by "ocean mode"
486! bugfix: duplicate error message 56 removed,
487! check of data_output_format and do3d_compress removed
488!
489! 1322 2014-03-20 16:38:49Z raasch
490! some REAL constants defined as wp-kind
491!
492! 1320 2014-03-20 08:40:49Z raasch
493! Kind-parameters added to all INTEGER and REAL declaration statements,
494! kinds are defined in new module kinds,
495! revision history before 2012 removed,
496! comment fields (!:) to be used for variable explanations added to
497! all variable declaration statements
498!
499! 1308 2014-03-13 14:58:42Z fricke
500! +netcdf_data_format_save
501! Calculate fixed number of output time levels for parallel netcdf output.
502! For masked data, parallel netcdf output is not tested so far, hence
503! netcdf_data_format is switched back to non-paralell output.
504!
505! 1299 2014-03-06 13:15:21Z heinze
506! enable usage of large_scale subsidence in combination with large_scale_forcing
507! output for profile of large scale vertical velocity w_subs added
508!
509! 1276 2014-01-15 13:40:41Z heinze
510! Use LSF_DATA also in case of Dirichlet bottom boundary condition for scalars
511!
512! 1241 2013-10-30 11:36:58Z heinze
513! output for profiles of ug and vg added
514! set surface_heatflux and surface_waterflux also in dependence on
515! large_scale_forcing
516! checks for nudging and large scale forcing from external file
517!
518! 1236 2013-09-27 07:21:13Z raasch
519! check number of spectra levels
520!
521! 1216 2013-08-26 09:31:42Z raasch
522! check for transpose_compute_overlap (temporary)
523!
524! 1214 2013-08-21 12:29:17Z kanani
525! additional check for simultaneous use of vertical grid stretching
526! and particle advection
527!
528! 1212 2013-08-15 08:46:27Z raasch
529! checks for poisfft_hybrid removed
530!
531! 1210 2013-08-14 10:58:20Z raasch
532! check for fftw
533!
534! 1179 2013-06-14 05:57:58Z raasch
535! checks and settings of buoyancy parameters and switches revised,
536! initial profile for rho_ocean added to hom (id=77)
537!
538! 1174 2013-05-31 10:28:08Z gryschka
539! Bugfix in computing initial profiles for ug, vg, lad, q in case of Atmosphere
540!
541! 1159 2013-05-21 11:58:22Z fricke
542! bc_lr/ns_dirneu/neudir removed
543!
544! 1115 2013-03-26 18:16:16Z hoffmann
545! unused variables removed
546! drizzle can be used without precipitation
547!
548! 1111 2013-03-08 23:54:10Z raasch
549! ibc_p_b = 2 removed
550!
551! 1103 2013-02-20 02:15:53Z raasch
552! Bugfix: turbulent inflow must not require cyclic fill in restart runs
553!
554! 1092 2013-02-02 11:24:22Z raasch
555! unused variables removed
556!
557! 1069 2012-11-28 16:18:43Z maronga
558! allow usage of topography in combination with cloud physics
559!
560! 1065 2012-11-22 17:42:36Z hoffmann
561! Bugfix: It is not allowed to use cloud_scheme = seifert_beheng without
562!         precipitation in order to save computational resources.
563!
564! 1060 2012-11-21 07:19:51Z raasch
565! additional check for parameter turbulent_inflow
566!
567! 1053 2012-11-13 17:11:03Z hoffmann
568! necessary changes for the new two-moment cloud physics scheme added:
569! - check cloud physics scheme (Kessler or Seifert and Beheng)
570! - plant_canopy is not allowed
571! - currently, only cache loop_optimization is allowed
572! - initial profiles of nr, qr
573! - boundary condition of nr, qr
574! - check output quantities (qr, nr, prr)
575!
576! 1036 2012-10-22 13:43:42Z raasch
577! code put under GPL (PALM 3.9)
578!
579! 1031/1034 2012-10-22 11:32:49Z raasch
580! check of netcdf4 parallel file support
581!
582! 1019 2012-09-28 06:46:45Z raasch
583! non-optimized version of prognostic_equations not allowed any more
584!
585! 1015 2012-09-27 09:23:24Z raasch
586! acc allowed for loop optimization,
587! checks for adjustment of mixing length to the Prandtl mixing length removed
588!
589! 1003 2012-09-14 14:35:53Z raasch
590! checks for cases with unequal subdomain sizes removed
591!
592! 1001 2012-09-13 14:08:46Z raasch
593! all actions concerning leapfrog- and upstream-spline-scheme removed
594!
595! 996 2012-09-07 10:41:47Z raasch
596! little reformatting
597
598! 978 2012-08-09 08:28:32Z fricke
599! setting of bc_lr/ns_dirneu/neudir
600! outflow damping layer removed
601! check for z0h*
602! check for pt_damping_width
603!
604! 964 2012-07-26 09:14:24Z raasch
605! check of old profil-parameters removed
606!
607! 940 2012-07-09 14:31:00Z raasch
608! checks for parameter neutral
609!
610! 924 2012-06-06 07:44:41Z maronga
611! Bugfix: preprocessor directives caused error during compilation
612!
613! 892 2012-05-02 13:51:44Z maronga
614! Bugfix for parameter file check ( excluding __netcdf4 )
615!
616! 866 2012-03-28 06:44:41Z raasch
617! use only 60% of the geostrophic wind as translation speed in case of Galilean
618! transformation and use_ug_for_galilei_tr = .T. in order to mimimize the
619! timestep
620!
621! 861 2012-03-26 14:18:34Z suehring
622! Check for topography and ws-scheme removed.
623! Check for loop_optimization = 'vector' and ws-scheme removed.
624!
625! 845 2012-03-07 10:23:05Z maronga
626! Bugfix: exclude __netcdf4 directive part from namelist file check compilation
627!
628! 828 2012-02-21 12:00:36Z raasch
629! check of collision_kernel extended
630!
631! 825 2012-02-19 03:03:44Z raasch
632! check for collision_kernel and curvature_solution_effects
633!
634! 809 2012-01-30 13:32:58Z maronga
635! Bugfix: replaced .AND. and .NOT. with && and ! in the preprocessor directives
636!
637! 807 2012-01-25 11:53:51Z maronga
638! New cpp directive "__check" implemented which is used by check_namelist_files
639!
640! Revision 1.1  1997/08/26 06:29:23  raasch
641! Initial revision
642!
643!
644! Description:
645! ------------
646!> Check control parameters and deduce further quantities.
647!------------------------------------------------------------------------------!
648 SUBROUTINE check_parameters
649
650
651    USE arrays_3d
652    USE chemistry_model_mod,                                                   &
653        ONLY:  chem_boundary_conds, chem_check_data_output,                    &
654               chem_check_data_output_pr, chem_species                 
655    USE chem_modules
656    USE cloud_parameters
657    USE constants
658    USE control_parameters
659    USE dvrp_variables
660    USE grid_variables
661    USE gust_mod,                                                              &
662        ONLY: gust_check_data_output, gust_check_data_output_pr,               &
663              gust_check_parameters, gust_module_enabled
664    USE indices
665    USE land_surface_model_mod,                                                &
666        ONLY: lsm_check_data_output, lsm_check_data_output_pr,                 &
667              lsm_check_parameters
668
669    USE lsf_nudging_mod,                                                       &
670        ONLY:  lsf_nudging_check_parameters, lsf_nudging_check_data_output_pr
671
672    USE kinds
673
674    USE model_1d_mod,                                                          &
675        ONLY:  damp_level_1d, damp_level_ind_1d
676
677    USE netcdf_data_input_mod,                                                 &
678        ONLY:  init_model, input_pids_static, netcdf_data_input_check_dynamic, &
679               netcdf_data_input_check_static
680
681    USE netcdf_interface,                                                      &
682        ONLY:  dopr_unit, do2d_unit, do3d_unit, netcdf_data_format,            &
683               netcdf_data_format_string, dots_unit, heatflux_output_unit,     &
684               waterflux_output_unit, momentumflux_output_unit
685
686    USE particle_attributes
687    USE pegrid
688    USE plant_canopy_model_mod,                                                &
689        ONLY:  pcm_check_data_output, pcm_check_parameters
690
691    USE pmc_interface,                                                         &
692        ONLY:  cpl_id, nested_run
693
694    USE profil_parameter
695    USE radiation_model_mod,                                                   &
696        ONLY:  radiation, radiation_check_data_output,                         &
697               radiation_check_data_output_pr, radiation_check_parameters
698    USE spectra_mod,                                                           &
699        ONLY:  calculate_spectra, spectra_check_parameters
700    USE statistics
701    USE subsidence_mod
702    USE statistics
703    USE synthetic_turbulence_generator_mod,                                    &
704        ONLY:  stg_check_parameters
705    USE transpose_indices
706    USE turbulence_closure_mod,                                                &
707        ONLY:  tcm_check_data_output, tcm_check_parameters
708    USE urban_surface_mod,                                                     &
709        ONLY:  usm_check_data_output, usm_check_parameters
710    USE uv_exposure_model_mod,                                                 &
711        ONLY:  uvem_check_data_output
712    USE wind_turbine_model_mod,                                                &
713        ONLY:  wtm_check_parameters
714    USE vertical_nesting_mod,                                                  &
715        ONLY:  vnested, vnest_check_parameters
716
717
718    IMPLICIT NONE
719
720    CHARACTER (LEN=varnamelength)  ::  var           !< variable name
721    CHARACTER (LEN=7)   ::  unit                     !< unit of variable
722    CHARACTER (LEN=8)   ::  date                     !< current date string
723    CHARACTER (LEN=10)  ::  time                     !< current time string
724    CHARACTER (LEN=20)  ::  ensemble_string          !< string containing number of ensemble member
725    CHARACTER (LEN=15)  ::  nest_string              !< string containing id of nested domain
726    CHARACTER (LEN=40)  ::  coupling_string          !< string containing type of coupling
727    CHARACTER (LEN=100) ::  action                   !< flag string
728
729    INTEGER(iwp) ::  i                               !< loop index
730    INTEGER(iwp) ::  ilen                            !< string length
731    INTEGER(iwp) ::  j                               !< loop index
732    INTEGER(iwp) ::  k                               !< loop index
733    INTEGER(iwp) ::  kk                              !< loop index
734    INTEGER(iwp) ::  netcdf_data_format_save         !< initial value of netcdf_data_format
735    INTEGER(iwp) ::  position                        !< index position of string
736    INTEGER(iwp) ::  lsp                             !< running index for chem spcs.
737
738    LOGICAL     ::  found                            !< flag, true if output variable is already marked for averaging
739
740    REAL(wp)    ::  dt_spinup_max                    !< maximum spinup timestep in nested domains
741    REAL(wp)    ::  dum                              !< dummy variable
742    REAL(wp)    ::  gradient                         !< local gradient
743    REAL(wp)    ::  remote = 0.0_wp                  !< MPI id of remote processor
744    REAL(wp)    ::  spinup_time_max                  !< maximum spinup time in nested domains
745    REAL(wp)    ::  time_to_be_simulated_from_reference_point  !< time to be simulated from reference point
746
747
748    CALL location_message( 'checking parameters', .FALSE. )
749!
750!-- At first, check static and dynamic input for consistency
751    CALL netcdf_data_input_check_dynamic
752    CALL netcdf_data_input_check_static
753!
754!-- Check for overlap combinations, which are not realized yet
755    IF ( transpose_compute_overlap )  THEN
756       IF ( numprocs == 1 )  STOP '+++ transpose-compute-overlap not implemented for single PE runs'
757    ENDIF
758
759!
760!-- Check the coupling mode
761    IF ( coupling_mode /= 'uncoupled'            .AND.                         &
762         coupling_mode /= 'precursor_atmos'      .AND.                         &
763         coupling_mode /= 'precursor_ocean'      .AND.                         &
764         coupling_mode /= 'vnested_crse'         .AND.                         &
765         coupling_mode /= 'vnested_fine'         .AND.                         &
766         coupling_mode /= 'atmosphere_to_ocean'  .AND.                         &
767         coupling_mode /= 'ocean_to_atmosphere' )  THEN
768       message_string = 'illegal coupling mode: ' // TRIM( coupling_mode )
769       CALL message( 'check_parameters', 'PA0002', 1, 2, 0, 6, 0 )
770    ENDIF
771
772!
773!-- Check if humidity is set to TRUE in case of the atmospheric run (for coupled runs)
774    IF ( coupling_mode == 'atmosphere_to_ocean' .AND. .NOT. humidity) THEN
775       message_string = ' Humidity has to be set to .T. in the _p3d file ' //  &
776                        'for coupled runs between ocean and atmosphere.'
777       CALL message( 'check_parameters', 'PA0476', 1, 2, 0, 6, 0 )
778    ENDIF
779   
780!
781!-- Check dt_coupling, restart_time, dt_restart, end_time, dx, dy, nx and ny
782    IF ( coupling_mode /= 'uncoupled'       .AND.                              &
783         coupling_mode(1:8) /= 'vnested_'   .AND.                              &
784         coupling_mode /= 'precursor_atmos' .AND.                              &
785         coupling_mode /= 'precursor_ocean' )  THEN
786
787       IF ( dt_coupling == 9999999.9_wp )  THEN
788          message_string = 'dt_coupling is not set but required for coup' //   &
789                           'ling mode "' //  TRIM( coupling_mode ) // '"'
790          CALL message( 'check_parameters', 'PA0003', 1, 2, 0, 6, 0 )
791       ENDIF
792
793#if defined( __parallel )
794
795
796       IF ( myid == 0 ) THEN
797          CALL MPI_SEND( dt_coupling, 1, MPI_REAL, target_id, 11, comm_inter,  &
798                         ierr )
799          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 11, comm_inter,       &
800                         status, ierr )
801       ENDIF
802       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
803
804       IF ( dt_coupling /= remote )  THEN
805          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
806                 '": dt_coupling = ', dt_coupling, '& is not equal to ',       &
807                 'dt_coupling_remote = ', remote
808          CALL message( 'check_parameters', 'PA0004', 1, 2, 0, 6, 0 )
809       ENDIF
810       IF ( dt_coupling <= 0.0_wp )  THEN
811
812          IF ( myid == 0  ) THEN
813             CALL MPI_SEND( dt_max, 1, MPI_REAL, target_id, 19, comm_inter, ierr )
814             CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 19, comm_inter,    &
815                            status, ierr )
816          ENDIF
817          CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
818
819          dt_coupling = MAX( dt_max, remote )
820          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
821                 '": dt_coupling <= 0.0 & is not allowed and is reset to ',    &
822                 'MAX(dt_max(A,O)) = ', dt_coupling
823          CALL message( 'check_parameters', 'PA0005', 0, 1, 0, 6, 0 )
824       ENDIF
825
826       IF ( myid == 0 ) THEN
827          CALL MPI_SEND( restart_time, 1, MPI_REAL, target_id, 12, comm_inter, &
828                         ierr )
829          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 12, comm_inter,       &
830                         status, ierr )
831       ENDIF
832       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
833
834       IF ( restart_time /= remote )  THEN
835          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
836                 '": restart_time = ', restart_time, '& is not equal to ',     &
837                 'restart_time_remote = ', remote
838          CALL message( 'check_parameters', 'PA0006', 1, 2, 0, 6, 0 )
839       ENDIF
840
841       IF ( myid == 0 ) THEN
842          CALL MPI_SEND( dt_restart, 1, MPI_REAL, target_id, 13, comm_inter,   &
843                         ierr )
844          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 13, comm_inter,       &
845                         status, ierr )
846       ENDIF
847       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
848
849       IF ( dt_restart /= remote )  THEN
850          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
851                 '": dt_restart = ', dt_restart, '& is not equal to ',         &
852                 'dt_restart_remote = ', remote
853          CALL message( 'check_parameters', 'PA0007', 1, 2, 0, 6, 0 )
854       ENDIF
855
856       time_to_be_simulated_from_reference_point = end_time-coupling_start_time
857
858       IF  ( myid == 0 ) THEN
859          CALL MPI_SEND( time_to_be_simulated_from_reference_point, 1,         &
860                         MPI_REAL, target_id, 14, comm_inter, ierr )
861          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 14, comm_inter,       &
862                         status, ierr )
863       ENDIF
864       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
865
866       IF ( time_to_be_simulated_from_reference_point /= remote )  THEN
867          WRITE( message_string, * ) 'coupling mode "', TRIM( coupling_mode ), &
868                 '": time_to_be_simulated_from_reference_point = ',            &
869                 time_to_be_simulated_from_reference_point, '& is not equal ', &
870                 'to time_to_be_simulated_from_reference_point_remote = ',     &
871                 remote
872          CALL message( 'check_parameters', 'PA0008', 1, 2, 0, 6, 0 )
873       ENDIF
874
875       IF ( myid == 0 ) THEN
876          CALL MPI_SEND( dx, 1, MPI_REAL, target_id, 15, comm_inter, ierr )
877          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 15, comm_inter,       &
878                                                             status, ierr )
879       ENDIF
880       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
881
882
883       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
884
885          IF ( dx < remote ) THEN
886             WRITE( message_string, * ) 'coupling mode "',                     &
887                   TRIM( coupling_mode ),                                      &
888           '": dx in Atmosphere is not equal to or not larger than dx in ocean'
889             CALL message( 'check_parameters', 'PA0009', 1, 2, 0, 6, 0 )
890          ENDIF
891
892          IF ( (nx_a+1)*dx /= (nx_o+1)*remote )  THEN
893             WRITE( message_string, * ) 'coupling mode "',                     &
894                    TRIM( coupling_mode ),                                     &
895             '": Domain size in x-direction is not equal in ocean and atmosphere'
896             CALL message( 'check_parameters', 'PA0010', 1, 2, 0, 6, 0 )
897          ENDIF
898
899       ENDIF
900
901       IF ( myid == 0) THEN
902          CALL MPI_SEND( dy, 1, MPI_REAL, target_id, 16, comm_inter, ierr )
903          CALL MPI_RECV( remote, 1, MPI_REAL, target_id, 16, comm_inter,       &
904                         status, ierr )
905       ENDIF
906       CALL MPI_BCAST( remote, 1, MPI_REAL, 0, comm2d, ierr)
907
908       IF ( coupling_mode == 'atmosphere_to_ocean') THEN
909
910          IF ( dy < remote )  THEN
911             WRITE( message_string, * ) 'coupling mode "',                     &
912                    TRIM( coupling_mode ),                                     &
913                 '": dy in Atmosphere is not equal to or not larger than dy in ocean'
914             CALL message( 'check_parameters', 'PA0011', 1, 2, 0, 6, 0 )
915          ENDIF
916
917          IF ( (ny_a+1)*dy /= (ny_o+1)*remote )  THEN
918             WRITE( message_string, * ) 'coupling mode "',                     &
919                   TRIM( coupling_mode ),                                      &
920             '": Domain size in y-direction is not equal in ocean and atmosphere'
921             CALL message( 'check_parameters', 'PA0012', 1, 2, 0, 6, 0 )
922          ENDIF
923
924          IF ( MOD(nx_o+1,nx_a+1) /= 0 )  THEN
925             WRITE( message_string, * ) 'coupling mode "',                     &
926                   TRIM( coupling_mode ),                                      &
927             '": nx+1 in ocean is not divisible by nx+1 in',                   &
928             ' atmosphere without remainder'
929             CALL message( 'check_parameters', 'PA0339', 1, 2, 0, 6, 0 )
930          ENDIF
931
932          IF ( MOD(ny_o+1,ny_a+1) /= 0 )  THEN
933             WRITE( message_string, * ) 'coupling mode "',                     &
934                   TRIM( coupling_mode ),                                      &
935             '": ny+1 in ocean is not divisible by ny+1 in', &
936             ' atmosphere without remainder'
937
938             CALL message( 'check_parameters', 'PA0340', 1, 2, 0, 6, 0 )
939          ENDIF
940
941       ENDIF
942#else
943       WRITE( message_string, * ) 'coupling requires PALM to be compiled with',&
944            ' cpp-option "-D__parallel"'
945       CALL message( 'check_parameters', 'PA0141', 1, 2, 0, 6, 0 )
946#endif
947    ENDIF
948
949#if defined( __parallel )
950!
951!-- Exchange via intercommunicator
952    IF ( coupling_mode == 'atmosphere_to_ocean' .AND. myid == 0 )  THEN
953       CALL MPI_SEND( humidity, 1, MPI_LOGICAL, target_id, 19, comm_inter,     &
954                      ierr )
955    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' .AND. myid == 0)  THEN
956       CALL MPI_RECV( humidity_remote, 1, MPI_LOGICAL, target_id, 19,          &
957                      comm_inter, status, ierr )
958    ENDIF
959    CALL MPI_BCAST( humidity_remote, 1, MPI_LOGICAL, 0, comm2d, ierr)
960
961#endif
962
963!
964!-- User settings for restart times requires that "restart" has been given as
965!-- file activation string. Otherwise, binary output would not be saved by
966!-- palmrun.
967    IF (  ( restart_time /= 9999999.9_wp  .OR.  dt_restart /= 9999999.9_wp )   &
968         .AND.  .NOT. write_binary )  THEN
969       WRITE( message_string, * ) 'manual restart settings requires file ',    &
970                                  'activation string "restart"'
971       CALL message( 'check_parameters', 'PA0001', 1, 2, 0, 6, 0 )
972    ENDIF
973
974
975!
976!-- Generate the file header which is used as a header for most of PALM's
977!-- output files
978    CALL DATE_AND_TIME( date, time )
979    run_date = date(7:8)//'-'//date(5:6)//'-'//date(3:4)
980    run_time = time(1:2)//':'//time(3:4)//':'//time(5:6)
981    IF ( coupling_mode == 'uncoupled' )  THEN
982       coupling_string = ''
983    ELSEIF ( coupling_mode == 'vnested_crse' )  THEN
984       coupling_string = ' nested (coarse)'
985    ELSEIF ( coupling_mode == 'vnested_fine' )  THEN
986       coupling_string = ' nested (fine)'
987    ELSEIF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
988       coupling_string = ' coupled (atmosphere)'
989    ELSEIF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
990       coupling_string = ' coupled (ocean)'
991    ENDIF
992    IF ( ensemble_member_nr /= 0 )  THEN
993       WRITE( ensemble_string, '(2X,A,I2.2)' )  'en-no: ', ensemble_member_nr
994    ELSE
995       ensemble_string = ''
996    ENDIF
997    IF ( nested_run )  THEN
998       WRITE( nest_string, '(2X,A,I2.2)' )  'nest-id: ', cpl_id
999    ELSE
1000       nest_string = ''
1001    ENDIF
1002
1003    WRITE ( run_description_header,                                            &
1004            '(A,2X,A,2X,A,A,A,I2.2,A,A,A,2X,A,A,2X,A,1X,A)' )                  &
1005          TRIM( version ), TRIM( revision ), 'run: ',                          &
1006          TRIM( run_identifier ), '.', runnr, TRIM( coupling_string ),         &
1007          TRIM( nest_string ), TRIM( ensemble_string), 'host: ', TRIM( host ), &
1008          run_date, run_time
1009
1010!
1011!-- Check the general loop optimization method
1012    SELECT CASE ( TRIM( loop_optimization ) )
1013
1014       CASE ( 'cache', 'vector' )
1015          CONTINUE
1016
1017       CASE DEFAULT
1018          message_string = 'illegal value given for loop_optimization: "' //   &
1019                           TRIM( loop_optimization ) // '"'
1020          CALL message( 'check_parameters', 'PA0013', 1, 2, 0, 6, 0 )
1021
1022    END SELECT
1023
1024!
1025!-- Check topography setting (check for illegal parameter combinations)
1026    IF ( topography /= 'flat' )  THEN
1027       action = ' '
1028       IF ( scalar_advec /= 'pw-scheme' .AND. scalar_advec /= 'ws-scheme'      &
1029          )  THEN
1030          WRITE( action, '(A,A)' )  'scalar_advec = ', scalar_advec
1031       ENDIF
1032       IF ( momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme' )&
1033       THEN
1034          WRITE( action, '(A,A)' )  'momentum_advec = ', momentum_advec
1035       ENDIF
1036       IF ( psolver == 'sor' )  THEN
1037          WRITE( action, '(A,A)' )  'psolver = ', psolver
1038       ENDIF
1039       IF ( sloping_surface )  THEN
1040          WRITE( action, '(A)' )  'sloping surface = .TRUE.'
1041       ENDIF
1042       IF ( galilei_transformation )  THEN
1043          WRITE( action, '(A)' )  'galilei_transformation = .TRUE.'
1044       ENDIF
1045       IF ( cloud_physics )  THEN
1046          WRITE( action, '(A)' )  'cloud_physics = .TRUE.'
1047       ENDIF
1048       IF ( cloud_droplets )  THEN
1049          WRITE( action, '(A)' )  'cloud_droplets = .TRUE.'
1050       ENDIF
1051       IF ( .NOT. constant_flux_layer )  THEN
1052          WRITE( action, '(A)' )  'constant_flux_layer = .FALSE.'
1053       ENDIF
1054       IF ( action /= ' ' )  THEN
1055          message_string = 'a non-flat topography does not allow ' //          &
1056                           TRIM( action )
1057          CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 )
1058       ENDIF
1059
1060    ENDIF
1061
1062!
1063!-- Check turbulence closure setup
1064    CALL tcm_check_parameters
1065!
1066!-- Check approximation
1067    IF ( TRIM( approximation ) /= 'boussinesq'   .AND.                         &
1068         TRIM( approximation ) /= 'anelastic' )  THEN
1069       message_string = 'unknown approximation: approximation = "' //          &
1070                        TRIM( approximation ) // '"'
1071       CALL message( 'check_parameters', 'PA0446', 1, 2, 0, 6, 0 )
1072    ENDIF
1073
1074!
1075!-- Check approximation requirements
1076    IF ( TRIM( approximation ) == 'anelastic'   .AND.                          &
1077         TRIM( momentum_advec ) /= 'ws-scheme' )  THEN
1078       message_string = 'Anelastic approximation requires ' //                 &
1079                        'momentum_advec = "ws-scheme"'
1080       CALL message( 'check_parameters', 'PA0447', 1, 2, 0, 6, 0 )
1081    ENDIF
1082    IF ( TRIM( approximation ) == 'anelastic'   .AND.                          &
1083         TRIM( psolver ) == 'multigrid' )  THEN
1084       message_string = 'Anelastic approximation currently only supports ' //  &
1085                        'psolver = "poisfft", ' //                             &
1086                        'psolver = "sor" and ' //                              &
1087                        'psolver = "multigrid_noopt"'
1088       CALL message( 'check_parameters', 'PA0448', 1, 2, 0, 6, 0 )
1089    ENDIF
1090    IF ( TRIM( approximation ) == 'anelastic'   .AND.                          &
1091         conserve_volume_flow )  THEN
1092       message_string = 'Anelastic approximation is not allowed with ' //      &
1093                        'conserve_volume_flow = .TRUE.'
1094       CALL message( 'check_parameters', 'PA0449', 1, 2, 0, 6, 0 )
1095    ENDIF
1096
1097!
1098!-- Check flux input mode
1099    IF ( TRIM( flux_input_mode ) /= 'dynamic'    .AND.                         &
1100         TRIM( flux_input_mode ) /= 'kinematic'  .AND.                         &
1101         TRIM( flux_input_mode ) /= 'approximation-specific' )  THEN
1102       message_string = 'unknown flux input mode: flux_input_mode = "' //      &
1103                        TRIM( flux_input_mode ) // '"'
1104       CALL message( 'check_parameters', 'PA0450', 1, 2, 0, 6, 0 )
1105    ENDIF
1106!-- Set flux input mode according to approximation if applicable
1107    IF ( TRIM( flux_input_mode ) == 'approximation-specific' )  THEN
1108       IF ( TRIM( approximation ) == 'anelastic' )  THEN
1109          flux_input_mode = 'dynamic'
1110       ELSEIF ( TRIM( approximation ) == 'boussinesq' )  THEN
1111          flux_input_mode = 'kinematic'
1112       ENDIF
1113    ENDIF
1114
1115!
1116!-- Check flux output mode
1117    IF ( TRIM( flux_output_mode ) /= 'dynamic'    .AND.                        &
1118         TRIM( flux_output_mode ) /= 'kinematic'  .AND.                        &
1119         TRIM( flux_output_mode ) /= 'approximation-specific' )  THEN
1120       message_string = 'unknown flux output mode: flux_output_mode = "' //    &
1121                        TRIM( flux_output_mode ) // '"'
1122       CALL message( 'check_parameters', 'PA0451', 1, 2, 0, 6, 0 )
1123    ENDIF
1124!-- Set flux output mode according to approximation if applicable
1125    IF ( TRIM( flux_output_mode ) == 'approximation-specific' )  THEN
1126       IF ( TRIM( approximation ) == 'anelastic' )  THEN
1127          flux_output_mode = 'dynamic'
1128       ELSEIF ( TRIM( approximation ) == 'boussinesq' )  THEN
1129          flux_output_mode = 'kinematic'
1130       ENDIF
1131    ENDIF
1132
1133
1134!
1135!-- When the land- or urban-surface model is used, the flux output must be
1136!-- dynamic.
1137    IF ( land_surface  .OR.  urban_surface )  THEN
1138       flux_output_mode = 'dynamic'
1139    ENDIF
1140
1141!
1142!-- set the flux output units according to flux_output_mode
1143    IF ( TRIM( flux_output_mode ) == 'kinematic' ) THEN
1144        heatflux_output_unit              = 'K m/s'
1145        waterflux_output_unit             = 'kg/kg m/s'
1146        momentumflux_output_unit          = 'm2/s2'
1147    ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN
1148        heatflux_output_unit              = 'W/m2'
1149        waterflux_output_unit             = 'W/m2'
1150        momentumflux_output_unit          = 'N/m2'
1151    ENDIF
1152
1153!-- set time series output units for fluxes
1154    dots_unit(14:16) = heatflux_output_unit
1155    dots_unit(21)    = waterflux_output_unit
1156    dots_unit(19:20) = momentumflux_output_unit
1157
1158!
1159!-- Check ocean setting
1160    IF ( TRIM( coupling_mode ) == 'uncoupled'  .AND.                           &
1161         TRIM( coupling_char ) == '_O' .AND.                                   &
1162         .NOT. ocean )  THEN
1163
1164!
1165!--    Check whether an (uncoupled) atmospheric run has been declared as an
1166!--    ocean run (this setting is done via mrun-option -y)
1167       message_string = 'ocean = .F. does not allow coupling_char = "' //      &
1168                        TRIM( coupling_char ) // '" set by palmrun-option "-y"'
1169       CALL message( 'check_parameters', 'PA0317', 1, 2, 0, 6, 0 )
1170
1171    ENDIF
1172!
1173!-- Check cloud scheme
1174    IF ( cloud_scheme == 'saturation_adjust' )  THEN
1175       microphysics_sat_adjust = .TRUE.
1176       microphysics_seifert    = .FALSE.
1177       microphysics_kessler    = .FALSE.
1178       precipitation           = .FALSE.
1179    ELSEIF ( cloud_scheme == 'seifert_beheng' )  THEN
1180       microphysics_sat_adjust = .FALSE.
1181       microphysics_seifert    = .TRUE.
1182       microphysics_kessler    = .FALSE.
1183       microphysics_morrison  = .FALSE.
1184       precipitation           = .TRUE.
1185    ELSEIF ( cloud_scheme == 'kessler' )  THEN
1186       microphysics_sat_adjust = .FALSE.
1187       microphysics_seifert    = .FALSE.
1188       microphysics_kessler    = .TRUE.
1189       microphysics_morrison   = .FALSE.
1190       precipitation           = .TRUE.
1191    ELSEIF ( cloud_scheme == 'morrison' )  THEN
1192       microphysics_sat_adjust = .FALSE.
1193       microphysics_seifert    = .TRUE.
1194       microphysics_kessler    = .FALSE.
1195       microphysics_morrison   = .TRUE.
1196       precipitation           = .TRUE.
1197    ELSE
1198       message_string = 'unknown cloud microphysics scheme cloud_scheme ="' // &
1199                        TRIM( cloud_scheme ) // '"'
1200       CALL message( 'check_parameters', 'PA0357', 1, 2, 0, 6, 0 )
1201    ENDIF
1202!
1203!-- Check aerosol
1204    IF ( aerosol_bulk == 'nacl' )  THEN
1205       aerosol_nacl   = .TRUE.
1206       aerosol_c3h4o4 = .FALSE.
1207       aerosol_nh4no3 = .FALSE.
1208    ELSEIF ( aerosol_bulk == 'c3h4o4' )  THEN
1209       aerosol_nacl   = .FALSE.
1210       aerosol_c3h4o4 = .TRUE.
1211       aerosol_nh4no3 = .FALSE.
1212    ELSEIF ( aerosol_bulk == 'nh4no3' )  THEN
1213       aerosol_nacl   = .FALSE.
1214       aerosol_c3h4o4 = .FALSE.
1215       aerosol_nh4no3 = .TRUE.
1216    ELSE
1217       message_string = 'unknown aerosol = "' // TRIM( aerosol_bulk ) // '"'
1218       CALL message( 'check_parameters', 'PA0469', 1, 2, 0, 6, 0 )
1219    ENDIF
1220
1221!
1222!-- Check whether there are any illegal values
1223!-- Pressure solver:
1224    IF ( psolver /= 'poisfft'  .AND.  psolver /= 'sor'  .AND.                  &
1225         psolver /= 'multigrid'  .AND.  psolver /= 'multigrid_noopt' )  THEN
1226       message_string = 'unknown solver for perturbation pressure: psolver' // &
1227                        ' = "' // TRIM( psolver ) // '"'
1228       CALL message( 'check_parameters', 'PA0016', 1, 2, 0, 6, 0 )
1229    ENDIF
1230
1231    IF ( psolver(1:9) == 'multigrid' )  THEN
1232       IF ( cycle_mg == 'w' )  THEN
1233          gamma_mg = 2
1234       ELSEIF ( cycle_mg == 'v' )  THEN
1235          gamma_mg = 1
1236       ELSE
1237          message_string = 'unknown multigrid cycle: cycle_mg = "' //          &
1238                           TRIM( cycle_mg ) // '"'
1239          CALL message( 'check_parameters', 'PA0020', 1, 2, 0, 6, 0 )
1240       ENDIF
1241    ENDIF
1242
1243    IF ( fft_method /= 'singleton-algorithm'  .AND.                            &
1244         fft_method /= 'temperton-algorithm'  .AND.                            &
1245         fft_method /= 'fftw'                 .AND.                            &
1246         fft_method /= 'system-specific' )  THEN
1247       message_string = 'unknown fft-algorithm: fft_method = "' //             &
1248                        TRIM( fft_method ) // '"'
1249       CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 )
1250    ENDIF
1251
1252    IF( momentum_advec == 'ws-scheme' .AND.                                    &
1253        .NOT. call_psolver_at_all_substeps  ) THEN
1254        message_string = 'psolver must be called at each RK3 substep when "'// &
1255                      TRIM(momentum_advec) // ' "is used for momentum_advec'
1256        CALL message( 'check_parameters', 'PA0344', 1, 2, 0, 6, 0 )
1257    END IF
1258!
1259!-- Advection schemes:
1260    IF ( momentum_advec /= 'pw-scheme'  .AND.                                  & 
1261         momentum_advec /= 'ws-scheme'  .AND.                                  &
1262         momentum_advec /= 'up-scheme' )                                       &
1263    THEN
1264       message_string = 'unknown advection scheme: momentum_advec = "' //      &
1265                        TRIM( momentum_advec ) // '"'
1266       CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 )
1267    ENDIF
1268    IF ( ( momentum_advec == 'ws-scheme' .OR.  scalar_advec == 'ws-scheme' )   &
1269           .AND. ( timestep_scheme == 'euler' .OR.                             &
1270                   timestep_scheme == 'runge-kutta-2' ) )                      &
1271    THEN
1272       message_string = 'momentum_advec or scalar_advec = "'                   &
1273         // TRIM( momentum_advec ) // '" is not allowed with ' //              &
1274         'timestep_scheme = "' // TRIM( timestep_scheme ) // '"'
1275       CALL message( 'check_parameters', 'PA0023', 1, 2, 0, 6, 0 )
1276    ENDIF
1277    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme' .AND. &
1278         scalar_advec /= 'bc-scheme' .AND. scalar_advec /= 'up-scheme' )       &
1279    THEN
1280       message_string = 'unknown advection scheme: scalar_advec = "' //        &
1281                        TRIM( scalar_advec ) // '"'
1282       CALL message( 'check_parameters', 'PA0024', 1, 2, 0, 6, 0 )
1283    ENDIF
1284    IF ( scalar_advec == 'bc-scheme'  .AND.  loop_optimization == 'cache' )    &
1285    THEN
1286       message_string = 'advection_scheme scalar_advec = "'                    &
1287         // TRIM( scalar_advec ) // '" not implemented for ' //                &
1288         'loop_optimization = "' // TRIM( loop_optimization ) // '"'
1289       CALL message( 'check_parameters', 'PA0026', 1, 2, 0, 6, 0 )
1290    ENDIF
1291
1292    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets  .AND.             &
1293         .NOT. use_upstream_for_tke  .AND.                                     &
1294         scalar_advec /= 'ws-scheme'                                           &
1295       )  THEN
1296       use_upstream_for_tke = .TRUE.
1297       message_string = 'use_upstream_for_tke is set to .TRUE. because ' //    &
1298                        'use_sgs_for_particles = .TRUE. '          //          &
1299                        'and scalar_advec /= ws-scheme'
1300       CALL message( 'check_parameters', 'PA0025', 0, 1, 0, 6, 0 )
1301    ENDIF
1302
1303!
1304!-- Set LOGICAL switches to enhance performance
1305    IF ( momentum_advec == 'ws-scheme' )  ws_scheme_mom = .TRUE.
1306    IF ( scalar_advec   == 'ws-scheme' )  ws_scheme_sca = .TRUE.
1307
1308
1309!
1310!-- Timestep schemes:
1311    SELECT CASE ( TRIM( timestep_scheme ) )
1312
1313       CASE ( 'euler' )
1314          intermediate_timestep_count_max = 1
1315
1316       CASE ( 'runge-kutta-2' )
1317          intermediate_timestep_count_max = 2
1318
1319       CASE ( 'runge-kutta-3' )
1320          intermediate_timestep_count_max = 3
1321
1322       CASE DEFAULT
1323          message_string = 'unknown timestep scheme: timestep_scheme = "' //   &
1324                           TRIM( timestep_scheme ) // '"'
1325          CALL message( 'check_parameters', 'PA0027', 1, 2, 0, 6, 0 )
1326
1327    END SELECT
1328
1329    IF ( (momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme')   &
1330         .AND. timestep_scheme(1:5) == 'runge' ) THEN
1331       message_string = 'momentum advection scheme "' // &
1332                        TRIM( momentum_advec ) // '" & does not work with ' // &
1333                        'timestep_scheme "' // TRIM( timestep_scheme ) // '"'
1334       CALL message( 'check_parameters', 'PA0029', 1, 2, 0, 6, 0 )
1335    ENDIF
1336!
1337!-- Check for proper settings for microphysics
1338    IF ( cloud_physics  .AND.  cloud_droplets )  THEN
1339       message_string = 'cloud_physics = .TRUE. is not allowed with ' //       &
1340                        'cloud_droplets = .TRUE.'
1341       CALL message( 'check_parameters', 'PA0442', 1, 2, 0, 6, 0 )
1342    ENDIF
1343
1344!
1345!-- Collision kernels:
1346    SELECT CASE ( TRIM( collision_kernel ) )
1347
1348       CASE ( 'hall', 'hall_fast' )
1349          hall_kernel = .TRUE.
1350
1351       CASE ( 'wang', 'wang_fast' )
1352          wang_kernel = .TRUE.
1353
1354       CASE ( 'none' )
1355
1356
1357       CASE DEFAULT
1358          message_string = 'unknown collision kernel: collision_kernel = "' // &
1359                           TRIM( collision_kernel ) // '"'
1360          CALL message( 'check_parameters', 'PA0350', 1, 2, 0, 6, 0 )
1361
1362    END SELECT
1363    IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
1364
1365!
1366!-- Initializing actions must have been set by the user
1367    IF ( TRIM( initializing_actions ) == '' )  THEN
1368       message_string = 'no value specified for initializing_actions'
1369       CALL message( 'check_parameters', 'PA0149', 1, 2, 0, 6, 0 )
1370    ENDIF
1371
1372    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
1373         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1374!
1375!--    No restart run: several initialising actions are possible
1376       action = initializing_actions
1377       DO  WHILE ( TRIM( action ) /= '' )
1378          position = INDEX( action, ' ' )
1379          SELECT CASE ( action(1:position-1) )
1380
1381             CASE ( 'set_constant_profiles', 'set_1d-model_profiles',          &
1382                    'by_user', 'initialize_vortex', 'initialize_ptanom',       &
1383                    'initialize_bubble', 'inifor' )
1384                action = action(position+1:)
1385
1386             CASE DEFAULT
1387                message_string = 'initializing_action = "' //                  &
1388                                 TRIM( action ) // '" unknown or not allowed'
1389                CALL message( 'check_parameters', 'PA0030', 1, 2, 0, 6, 0 )
1390
1391          END SELECT
1392       ENDDO
1393    ENDIF
1394
1395    IF ( TRIM( initializing_actions ) == 'initialize_vortex'  .AND.            &
1396         conserve_volume_flow ) THEN
1397         message_string = 'initializing_actions = "initialize_vortex"' //      &
1398                        ' is not allowed with conserve_volume_flow = .T.'
1399       CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )
1400    ENDIF
1401
1402
1403    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.    &
1404         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
1405       message_string = 'initializing_actions = "set_constant_profiles"' //    &
1406                        ' and "set_1d-model_profiles" are not allowed ' //     &
1407                        'simultaneously'
1408       CALL message( 'check_parameters', 'PA0031', 1, 2, 0, 6, 0 )
1409    ENDIF
1410
1411    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.    &
1412         INDEX( initializing_actions, 'by_user' ) /= 0 )  THEN
1413       message_string = 'initializing_actions = "set_constant_profiles"' //    &
1414                        ' and "by_user" are not allowed simultaneously'
1415       CALL message( 'check_parameters', 'PA0032', 1, 2, 0, 6, 0 )
1416    ENDIF
1417
1418    IF ( INDEX( initializing_actions, 'by_user' ) /= 0  .AND.                  &
1419         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
1420       message_string = 'initializing_actions = "by_user" and ' //             &
1421                        '"set_1d-model_profiles" are not allowed simultaneously'
1422       CALL message( 'check_parameters', 'PA0033', 1, 2, 0, 6, 0 )
1423    ENDIF
1424!
1425!-- In case of spinup and nested run, spinup end time must be identical
1426!-- in order to have synchronously running simulations.
1427    IF ( nested_run )  THEN
1428#if defined( __parallel )
1429       CALL MPI_ALLREDUCE( spinup_time, spinup_time_max, 1, MPI_REAL,          &
1430                           MPI_MAX, MPI_COMM_WORLD, ierr )
1431       CALL MPI_ALLREDUCE( dt_spinup,   dt_spinup_max,   1, MPI_REAL,          &
1432                           MPI_MAX, MPI_COMM_WORLD, ierr )
1433       IF ( spinup_time /= spinup_time_max  .OR.  dt_spinup /= dt_spinup_max ) &
1434       THEN
1435          message_string = 'In case of nesting, spinup_time and ' //           &
1436                           'dt_spinup must be identical in all parent ' //     &
1437                           'and child domains.'
1438          CALL message( 'check_parameters', 'PA0489', 3, 2, 0, 6, 0 )
1439       ENDIF
1440#endif
1441    ENDIF
1442
1443    IF ( cloud_physics  .AND.  .NOT.  humidity )  THEN
1444       WRITE( message_string, * ) 'cloud_physics = ', cloud_physics, ' is ',   &
1445              'not allowed with humidity = ', humidity
1446       CALL message( 'check_parameters', 'PA0034', 1, 2, 0, 6, 0 )
1447    ENDIF
1448
1449    IF ( humidity  .AND.  sloping_surface )  THEN
1450       message_string = 'humidity = .TRUE. and sloping_surface = .TRUE. ' //   &
1451                        'are not allowed simultaneously'
1452       CALL message( 'check_parameters', 'PA0036', 1, 2, 0, 6, 0 )
1453    ENDIF
1454
1455!
1456!-- Check for synthetic turbulence generator settings
1457    CALL stg_check_parameters
1458!
1459!-- When plant canopy model is used, peform addtional checks
1460    IF ( plant_canopy )  CALL pcm_check_parameters
1461
1462!
1463!-- Additional checks for spectra
1464    IF ( calculate_spectra )  CALL spectra_check_parameters
1465!
1466!-- When land surface model is used, perform additional checks
1467    IF ( land_surface )  CALL lsm_check_parameters
1468
1469!
1470!-- When urban surface model is used, perform additional checks
1471    IF ( urban_surface )  CALL usm_check_parameters
1472
1473!
1474!-- If wind turbine model is used, perform additional checks
1475    IF ( wind_turbine )  CALL wtm_check_parameters
1476!
1477!-- When radiation model is used, peform addtional checks
1478    IF ( radiation )  CALL radiation_check_parameters
1479!
1480!-- When gust module is used, perform additional checks
1481    IF ( gust_module_enabled )  CALL gust_check_parameters
1482!
1483!-- When large-scale forcing or nudging is used, peform addtional checks
1484    IF ( large_scale_forcing  .OR.  nudging )  CALL lsf_nudging_check_parameters
1485
1486!
1487!-- In case of no model continuation run, check initialising parameters and
1488!-- deduce further quantities
1489    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1490
1491!
1492!--    Initial profiles for 1D and 3D model, respectively (u,v further below)
1493       pt_init = pt_surface
1494       IF ( humidity       )  q_init  = q_surface
1495       IF ( ocean          )  sa_init = sa_surface
1496       IF ( passive_scalar )  s_init  = s_surface
1497       IF ( air_chemistry )  THEN
1498         DO lsp = 1, nvar
1499            chem_species(lsp)%conc_pr_init = cs_surface(lsp)     
1500         ENDDO
1501       ENDIF
1502!
1503!--
1504!--    If required, compute initial profile of the geostrophic wind
1505!--    (component ug)
1506       i = 1
1507       gradient = 0.0_wp
1508
1509       IF (  .NOT.  ocean )  THEN
1510
1511          ug_vertical_gradient_level_ind(1) = 0
1512          ug(0) = ug_surface
1513          DO  k = 1, nzt+1
1514             IF ( i < 11 )  THEN
1515                IF ( ug_vertical_gradient_level(i) < zu(k)  .AND.              &
1516                     ug_vertical_gradient_level(i) >= 0.0_wp )  THEN
1517                   gradient = ug_vertical_gradient(i) / 100.0_wp
1518                   ug_vertical_gradient_level_ind(i) = k - 1
1519                   i = i + 1
1520                ENDIF
1521             ENDIF
1522             IF ( gradient /= 0.0_wp )  THEN
1523                IF ( k /= 1 )  THEN
1524                   ug(k) = ug(k-1) + dzu(k) * gradient
1525                ELSE
1526                   ug(k) = ug_surface + dzu(k) * gradient
1527                ENDIF
1528             ELSE
1529                ug(k) = ug(k-1)
1530             ENDIF
1531          ENDDO
1532
1533       ELSE
1534
1535          ug_vertical_gradient_level_ind(1) = nzt+1
1536          ug(nzt+1) = ug_surface
1537          DO  k = nzt, nzb, -1
1538             IF ( i < 11 )  THEN
1539                IF ( ug_vertical_gradient_level(i) > zu(k)  .AND.              &
1540                     ug_vertical_gradient_level(i) <= 0.0_wp )  THEN
1541                   gradient = ug_vertical_gradient(i) / 100.0_wp
1542                   ug_vertical_gradient_level_ind(i) = k + 1
1543                   i = i + 1
1544                ENDIF
1545             ENDIF
1546             IF ( gradient /= 0.0_wp )  THEN
1547                IF ( k /= nzt )  THEN
1548                   ug(k) = ug(k+1) - dzu(k+1) * gradient
1549                ELSE
1550                   ug(k)   = ug_surface - 0.5_wp * dzu(k+1) * gradient
1551                   ug(k+1) = ug_surface + 0.5_wp * dzu(k+1) * gradient
1552                ENDIF
1553             ELSE
1554                ug(k) = ug(k+1)
1555             ENDIF
1556          ENDDO
1557
1558       ENDIF
1559
1560!
1561!--    In case of no given gradients for ug, choose a zero gradient
1562       IF ( ug_vertical_gradient_level(1) == -9999999.9_wp )  THEN
1563          ug_vertical_gradient_level(1) = 0.0_wp
1564       ENDIF
1565
1566!
1567!--
1568!--    If required, compute initial profile of the geostrophic wind
1569!--    (component vg)
1570       i = 1
1571       gradient = 0.0_wp
1572
1573       IF (  .NOT.  ocean )  THEN
1574
1575          vg_vertical_gradient_level_ind(1) = 0
1576          vg(0) = vg_surface
1577          DO  k = 1, nzt+1
1578             IF ( i < 11 )  THEN
1579                IF ( vg_vertical_gradient_level(i) < zu(k)  .AND.              &
1580                     vg_vertical_gradient_level(i) >= 0.0_wp )  THEN
1581                   gradient = vg_vertical_gradient(i) / 100.0_wp
1582                   vg_vertical_gradient_level_ind(i) = k - 1
1583                   i = i + 1
1584                ENDIF
1585             ENDIF
1586             IF ( gradient /= 0.0_wp )  THEN
1587                IF ( k /= 1 )  THEN
1588                   vg(k) = vg(k-1) + dzu(k) * gradient
1589                ELSE
1590                   vg(k) = vg_surface + dzu(k) * gradient
1591                ENDIF
1592             ELSE
1593                vg(k) = vg(k-1)
1594             ENDIF
1595          ENDDO
1596
1597       ELSE
1598
1599          vg_vertical_gradient_level_ind(1) = nzt+1
1600          vg(nzt+1) = vg_surface
1601          DO  k = nzt, nzb, -1
1602             IF ( i < 11 )  THEN
1603                IF ( vg_vertical_gradient_level(i) > zu(k)  .AND.              &
1604                     vg_vertical_gradient_level(i) <= 0.0_wp )  THEN
1605                   gradient = vg_vertical_gradient(i) / 100.0_wp
1606                   vg_vertical_gradient_level_ind(i) = k + 1
1607                   i = i + 1
1608                ENDIF
1609             ENDIF
1610             IF ( gradient /= 0.0_wp )  THEN
1611                IF ( k /= nzt )  THEN
1612                   vg(k) = vg(k+1) - dzu(k+1) * gradient
1613                ELSE
1614                   vg(k)   = vg_surface - 0.5_wp * dzu(k+1) * gradient
1615                   vg(k+1) = vg_surface + 0.5_wp * dzu(k+1) * gradient
1616                ENDIF
1617             ELSE
1618                vg(k) = vg(k+1)
1619             ENDIF
1620          ENDDO
1621
1622       ENDIF
1623
1624!
1625!--    In case of no given gradients for vg, choose a zero gradient
1626       IF ( vg_vertical_gradient_level(1) == -9999999.9_wp )  THEN
1627          vg_vertical_gradient_level(1) = 0.0_wp
1628       ENDIF
1629
1630!
1631!--    Let the initial wind profiles be the calculated ug/vg profiles or
1632!--    interpolate them from wind profile data (if given)
1633       IF ( u_profile(1) == 9999999.9_wp  .AND.  v_profile(1) == 9999999.9_wp )  THEN
1634
1635          u_init = ug
1636          v_init = vg
1637
1638       ELSEIF ( u_profile(1) == 0.0_wp  .AND.  v_profile(1) == 0.0_wp )  THEN
1639
1640          IF ( uv_heights(1) /= 0.0_wp )  THEN
1641             message_string = 'uv_heights(1) must be 0.0'
1642             CALL message( 'check_parameters', 'PA0345', 1, 2, 0, 6, 0 )
1643          ENDIF
1644
1645          IF ( omega /= 0.0_wp )  THEN
1646             message_string = 'Coriolis force must be switched off (by setting omega=0.0)' //  &
1647                              ' when prescribing the forcing by u_profile and v_profile'
1648             CALL message( 'check_parameters', 'PA0347', 1, 2, 0, 6, 0 )
1649          ENDIF
1650
1651          use_prescribed_profile_data = .TRUE.
1652
1653          kk = 1
1654          u_init(0) = 0.0_wp
1655          v_init(0) = 0.0_wp
1656
1657          DO  k = 1, nz+1
1658
1659             IF ( kk < 100 )  THEN
1660                DO  WHILE ( uv_heights(kk+1) <= zu(k) )
1661                   kk = kk + 1
1662                   IF ( kk == 100 )  EXIT
1663                ENDDO
1664             ENDIF
1665
1666             IF ( kk < 100  .AND.  uv_heights(kk+1) /= 9999999.9_wp )  THEN
1667                u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
1668                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
1669                                       ( u_profile(kk+1) - u_profile(kk) )
1670                v_init(k) = v_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
1671                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
1672                                       ( v_profile(kk+1) - v_profile(kk) )
1673             ELSE
1674                u_init(k) = u_profile(kk)
1675                v_init(k) = v_profile(kk)
1676             ENDIF
1677
1678          ENDDO
1679
1680       ELSE
1681
1682          message_string = 'u_profile(1) and v_profile(1) must be 0.0'
1683          CALL message( 'check_parameters', 'PA0346', 1, 2, 0, 6, 0 )
1684
1685       ENDIF
1686
1687!
1688!--    Compute initial temperature profile using the given temperature gradients
1689       IF (  .NOT.  neutral )  THEN
1690          CALL init_vertical_profiles( pt_vertical_gradient_level_ind,          &
1691                                       pt_vertical_gradient_level,              &
1692                                       pt_vertical_gradient, pt_init,           &
1693                                       pt_surface, bc_pt_t_val )
1694       ENDIF
1695!
1696!--    Compute initial humidity profile using the given humidity gradients
1697       IF ( humidity )  THEN
1698          CALL init_vertical_profiles( q_vertical_gradient_level_ind,          &
1699                                       q_vertical_gradient_level,              &
1700                                       q_vertical_gradient, q_init,            &
1701                                       q_surface, bc_q_t_val )
1702       ENDIF
1703!
1704!--    Compute initial scalar profile using the given scalar gradients
1705       IF ( passive_scalar )  THEN
1706          CALL init_vertical_profiles( s_vertical_gradient_level_ind,          &
1707                                       s_vertical_gradient_level,              &
1708                                       s_vertical_gradient, s_init,            &
1709                                       s_surface, bc_s_t_val )
1710       ENDIF
1711!
1712!--    Compute initial chemistry profile using the given chemical species gradients
1713       IF ( air_chemistry )  THEN                                                         
1714         DO lsp = 1, nvar
1715          CALL init_vertical_profiles( cs_vertical_gradient_level_ind(lsp,:),  &
1716                                       cs_vertical_gradient_level(lsp,:),      &
1717                                       cs_vertical_gradient(lsp,:),            &
1718                                       chem_species(lsp)%conc_pr_init,         &
1719                                       cs_surface(lsp), bc_cs_t_val(lsp) )
1720         ENDDO
1721       ENDIF
1722!
1723!
1724!--    If required, compute initial salinity profile using the given salinity
1725!--    gradients
1726       IF ( ocean )  THEN
1727          CALL init_vertical_profiles( sa_vertical_gradient_level_ind,          &
1728                                       sa_vertical_gradient_level,              &
1729                                       sa_vertical_gradient, sa_init,           &
1730                                       sa_surface, dum )
1731       ENDIF
1732
1733
1734    ENDIF
1735
1736!
1737!-- Check if the control parameter use_subsidence_tendencies is used correctly
1738    IF ( use_subsidence_tendencies  .AND.  .NOT.  large_scale_subsidence )  THEN
1739       message_string = 'The usage of use_subsidence_tendencies ' //           &
1740                            'requires large_scale_subsidence = .T..'
1741       CALL message( 'check_parameters', 'PA0396', 1, 2, 0, 6, 0 )
1742    ELSEIF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_forcing )  THEN
1743       message_string = 'The usage of use_subsidence_tendencies ' //           &
1744                            'requires large_scale_forcing = .T..'
1745       CALL message( 'check_parameters', 'PA0397', 1, 2, 0, 6, 0 )
1746    ENDIF
1747
1748!
1749!-- Initialize large scale subsidence if required
1750    If ( large_scale_subsidence )  THEN
1751       IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp  .AND.            &
1752                                     .NOT.  large_scale_forcing )  THEN
1753          CALL init_w_subsidence
1754       ENDIF
1755!
1756!--    In case large_scale_forcing is used, profiles for subsidence velocity
1757!--    are read in from file LSF_DATA
1758
1759       IF ( subs_vertical_gradient_level(1) == -9999999.9_wp  .AND.            &
1760            .NOT.  large_scale_forcing )  THEN
1761          message_string = 'There is no default large scale vertical ' //      &
1762                           'velocity profile set. Specify the subsidence ' //  &
1763                           'velocity profile via subs_vertical_gradient ' //   &
1764                           'and subs_vertical_gradient_level.'
1765          CALL message( 'check_parameters', 'PA0380', 1, 2, 0, 6, 0 )
1766       ENDIF
1767    ELSE
1768        IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp )  THEN
1769           message_string = 'Enable usage of large scale subsidence by ' //    &
1770                            'setting large_scale_subsidence = .T..'
1771          CALL message( 'check_parameters', 'PA0381', 1, 2, 0, 6, 0 )
1772        ENDIF
1773    ENDIF
1774
1775!
1776!-- Overwrite latitude if necessary and compute Coriolis parameter.
1777!-- To do - move initialization of f and fs to coriolis_mod.
1778    IF ( input_pids_static )  THEN
1779       latitude  = init_model%latitude
1780       longitude = init_model%longitude
1781    ENDIF
1782
1783    f  = 2.0_wp * omega * SIN( latitude / 180.0_wp * pi )
1784    fs = 2.0_wp * omega * COS( latitude / 180.0_wp * pi )
1785
1786!
1787!-- Check and set buoyancy related parameters and switches
1788    IF ( reference_state == 'horizontal_average' )  THEN
1789       CONTINUE
1790    ELSEIF ( reference_state == 'initial_profile' )  THEN
1791       use_initial_profile_as_reference = .TRUE.
1792    ELSEIF ( reference_state == 'single_value' )  THEN
1793       use_single_reference_value = .TRUE.
1794       IF ( pt_reference == 9999999.9_wp )  pt_reference = pt_surface
1795       vpt_reference = pt_reference * ( 1.0_wp + 0.61_wp * q_surface )
1796    ELSE
1797       message_string = 'illegal value for reference_state: "' //              &
1798                        TRIM( reference_state ) // '"'
1799       CALL message( 'check_parameters', 'PA0056', 1, 2, 0, 6, 0 )
1800    ENDIF
1801
1802!
1803!-- Sign of buoyancy/stability terms
1804    IF ( ocean )  atmos_ocean_sign = -1.0_wp
1805
1806!
1807!-- Ocean version must use flux boundary conditions at the top
1808    IF ( ocean .AND. .NOT. use_top_fluxes )  THEN
1809       message_string = 'use_top_fluxes must be .TRUE. in ocean mode'
1810       CALL message( 'check_parameters', 'PA0042', 1, 2, 0, 6, 0 )
1811    ENDIF
1812
1813!
1814!-- In case of a given slope, compute the relevant quantities
1815    IF ( alpha_surface /= 0.0_wp )  THEN
1816       IF ( ABS( alpha_surface ) > 90.0_wp )  THEN
1817          WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface,   &
1818                                     ' ) must be < 90.0'
1819          CALL message( 'check_parameters', 'PA0043', 1, 2, 0, 6, 0 )
1820       ENDIF
1821       sloping_surface = .TRUE.
1822       cos_alpha_surface = COS( alpha_surface / 180.0_wp * pi )
1823       sin_alpha_surface = SIN( alpha_surface / 180.0_wp * pi )
1824    ENDIF
1825
1826!
1827!-- Check time step and cfl_factor
1828    IF ( dt /= -1.0_wp )  THEN
1829       IF ( dt <= 0.0_wp )  THEN
1830          WRITE( message_string, * ) 'dt = ', dt , ' <= 0.0'
1831          CALL message( 'check_parameters', 'PA0044', 1, 2, 0, 6, 0 )
1832       ENDIF
1833       dt_3d = dt
1834       dt_fixed = .TRUE.
1835    ENDIF
1836
1837    IF ( cfl_factor <= 0.0_wp  .OR.  cfl_factor > 1.0_wp )  THEN
1838       IF ( cfl_factor == -1.0_wp )  THEN
1839          IF ( timestep_scheme == 'runge-kutta-2' )  THEN
1840             cfl_factor = 0.8_wp
1841          ELSEIF ( timestep_scheme == 'runge-kutta-3' )  THEN
1842             cfl_factor = 0.9_wp
1843          ELSE
1844             cfl_factor = 0.9_wp
1845          ENDIF
1846       ELSE
1847          WRITE( message_string, * ) 'cfl_factor = ', cfl_factor,              &
1848                 ' out of range &0.0 < cfl_factor <= 1.0 is required'
1849          CALL message( 'check_parameters', 'PA0045', 1, 2, 0, 6, 0 )
1850       ENDIF
1851    ENDIF
1852
1853!
1854!-- Store simulated time at begin
1855    simulated_time_at_begin = simulated_time
1856
1857!
1858!-- Store reference time for coupled runs and change the coupling flag,
1859!-- if ...
1860    IF ( simulated_time == 0.0_wp )  THEN
1861       IF ( coupling_start_time == 0.0_wp )  THEN
1862          time_since_reference_point = 0.0_wp
1863       ELSEIF ( time_since_reference_point < 0.0_wp )  THEN
1864          run_coupled = .FALSE.
1865       ENDIF
1866    ENDIF
1867
1868!
1869!-- Set wind speed in the Galilei-transformed system
1870    IF ( galilei_transformation )  THEN
1871       IF ( use_ug_for_galilei_tr                    .AND.                     &
1872            ug_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
1873            ug_vertical_gradient(1) == 0.0_wp        .AND.                     &
1874            vg_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
1875            vg_vertical_gradient(1) == 0.0_wp )  THEN
1876          u_gtrans = ug_surface * 0.6_wp
1877          v_gtrans = vg_surface * 0.6_wp
1878       ELSEIF ( use_ug_for_galilei_tr  .AND.                                   &
1879                ( ug_vertical_gradient_level(1) /= 0.0_wp  .OR.                &
1880                ug_vertical_gradient(1) /= 0.0_wp ) )  THEN
1881          message_string = 'baroclinity (ug) not allowed simultaneously' //    &
1882                           ' with galilei transformation'
1883          CALL message( 'check_parameters', 'PA0046', 1, 2, 0, 6, 0 )
1884       ELSEIF ( use_ug_for_galilei_tr  .AND.                                   &
1885                ( vg_vertical_gradient_level(1) /= 0.0_wp  .OR.                &
1886                vg_vertical_gradient(1) /= 0.0_wp ) )  THEN
1887          message_string = 'baroclinity (vg) not allowed simultaneously' //    &
1888                           ' with galilei transformation'
1889          CALL message( 'check_parameters', 'PA0047', 1, 2, 0, 6, 0 )
1890       ELSE
1891          message_string = 'variable translation speed used for Galilei-' //   &
1892             'transformation, which may cause & instabilities in stably ' //   &
1893             'stratified regions'
1894          CALL message( 'check_parameters', 'PA0048', 0, 1, 0, 6, 0 )
1895       ENDIF
1896    ENDIF
1897
1898!
1899!-- In case of using a prandtl-layer, calculated (or prescribed) surface
1900!-- fluxes have to be used in the diffusion-terms
1901    IF ( constant_flux_layer )  use_surface_fluxes = .TRUE.
1902
1903!
1904!-- Check boundary conditions and set internal variables:
1905!-- Attention: the lateral boundary conditions have been already checked in
1906!-- parin
1907!
1908!-- Non-cyclic lateral boundaries require the multigrid method and Piascek-
1909!-- Willimas or Wicker - Skamarock advection scheme. Several schemes
1910!-- and tools do not work with non-cyclic boundary conditions.
1911    IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
1912       IF ( psolver(1:9) /= 'multigrid' )  THEN
1913          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1914                           'psolver = "' // TRIM( psolver ) // '"'
1915          CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 )
1916       ENDIF
1917       IF ( momentum_advec /= 'pw-scheme'  .AND.                               &
1918            momentum_advec /= 'ws-scheme' )  THEN
1919
1920          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1921                           'momentum_advec = "' // TRIM( momentum_advec ) // '"'
1922          CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 )
1923       ENDIF
1924       IF ( scalar_advec /= 'pw-scheme'  .AND.                                 &
1925            scalar_advec /= 'ws-scheme' )  THEN
1926          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1927                           'scalar_advec = "' // TRIM( scalar_advec ) // '"'
1928          CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 )
1929       ENDIF
1930       IF ( galilei_transformation )  THEN
1931          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1932                           'galilei_transformation = .T.'
1933          CALL message( 'check_parameters', 'PA0054', 1, 2, 0, 6, 0 )
1934       ENDIF
1935    ENDIF
1936
1937!
1938!-- Bottom boundary condition for the turbulent Kinetic energy
1939    IF ( bc_e_b == 'neumann' )  THEN
1940       ibc_e_b = 1
1941    ELSEIF ( bc_e_b == '(u*)**2+neumann' )  THEN
1942       ibc_e_b = 2
1943       IF ( .NOT. constant_flux_layer )  THEN
1944          bc_e_b = 'neumann'
1945          ibc_e_b = 1
1946          message_string = 'boundary condition bc_e_b changed to "' //         &
1947                           TRIM( bc_e_b ) // '"'
1948          CALL message( 'check_parameters', 'PA0057', 0, 1, 0, 6, 0 )
1949       ENDIF
1950    ELSE
1951       message_string = 'unknown boundary condition: bc_e_b = "' //            &
1952                        TRIM( bc_e_b ) // '"'
1953       CALL message( 'check_parameters', 'PA0058', 1, 2, 0, 6, 0 )
1954    ENDIF
1955
1956!
1957!-- Boundary conditions for perturbation pressure
1958    IF ( bc_p_b == 'dirichlet' )  THEN
1959       ibc_p_b = 0
1960    ELSEIF ( bc_p_b == 'neumann' )  THEN
1961       ibc_p_b = 1
1962    ELSE
1963       message_string = 'unknown boundary condition: bc_p_b = "' //            &
1964                        TRIM( bc_p_b ) // '"'
1965       CALL message( 'check_parameters', 'PA0059', 1, 2, 0, 6, 0 )
1966    ENDIF
1967
1968    IF ( bc_p_t == 'dirichlet' )  THEN
1969       ibc_p_t = 0
1970!-- TO_DO: later set bc_p_t to neumann before, in case of nested domain
1971    ELSEIF ( bc_p_t == 'neumann' .OR. bc_p_t == 'nested'  .OR.                 &
1972             bc_p_t == 'forcing'  )  THEN
1973       ibc_p_t = 1
1974    ELSE
1975       message_string = 'unknown boundary condition: bc_p_t = "' //            &
1976                        TRIM( bc_p_t ) // '"'
1977       CALL message( 'check_parameters', 'PA0061', 1, 2, 0, 6, 0 )
1978    ENDIF
1979
1980!
1981!-- Boundary conditions for potential temperature
1982    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
1983       ibc_pt_b = 2
1984    ELSE
1985       IF ( bc_pt_b == 'dirichlet' )  THEN
1986          ibc_pt_b = 0
1987       ELSEIF ( bc_pt_b == 'neumann' )  THEN
1988          ibc_pt_b = 1
1989       ELSE
1990          message_string = 'unknown boundary condition: bc_pt_b = "' //        &
1991                           TRIM( bc_pt_b ) // '"'
1992          CALL message( 'check_parameters', 'PA0062', 1, 2, 0, 6, 0 )
1993       ENDIF
1994    ENDIF
1995
1996    IF ( bc_pt_t == 'dirichlet' )  THEN
1997       ibc_pt_t = 0
1998    ELSEIF ( bc_pt_t == 'neumann' )  THEN
1999       ibc_pt_t = 1
2000    ELSEIF ( bc_pt_t == 'initial_gradient' )  THEN
2001       ibc_pt_t = 2
2002    ELSEIF ( bc_pt_t == 'nested'  .OR.  bc_pt_t == 'forcing' )  THEN
2003       ibc_pt_t = 3
2004    ELSE
2005       message_string = 'unknown boundary condition: bc_pt_t = "' //           &
2006                        TRIM( bc_pt_t ) // '"'
2007       CALL message( 'check_parameters', 'PA0063', 1, 2, 0, 6, 0 )
2008    ENDIF
2009
2010    IF ( ANY( wall_heatflux /= 0.0_wp )  .AND.                        &
2011         surface_heatflux == 9999999.9_wp )  THEN
2012       message_string = 'wall_heatflux additionally requires ' //     &
2013                        'setting of surface_heatflux'
2014       CALL message( 'check_parameters', 'PA0443', 1, 2, 0, 6, 0 )
2015    ENDIF
2016
2017!
2018!   This IF clause needs revision, got too complex!!
2019    IF ( surface_heatflux == 9999999.9_wp  )  THEN
2020       constant_heatflux = .FALSE.
2021       IF ( large_scale_forcing  .OR.  land_surface  .OR.  urban_surface )  THEN
2022          IF ( ibc_pt_b == 0 )  THEN
2023             constant_heatflux = .FALSE.
2024          ELSEIF ( ibc_pt_b == 1 )  THEN
2025             constant_heatflux = .TRUE.
2026             surface_heatflux = 0.0_wp
2027          ENDIF
2028       ENDIF
2029    ELSE
2030       constant_heatflux = .TRUE.
2031    ENDIF
2032
2033    IF ( top_heatflux     == 9999999.9_wp )  constant_top_heatflux = .FALSE.
2034
2035    IF ( neutral )  THEN
2036
2037       IF ( surface_heatflux /= 0.0_wp  .AND.                                  &
2038            surface_heatflux /= 9999999.9_wp )  THEN
2039          message_string = 'heatflux must not be set for pure neutral flow'
2040          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
2041       ENDIF
2042
2043       IF ( top_heatflux /= 0.0_wp  .AND.  top_heatflux /= 9999999.9_wp )      &
2044       THEN
2045          message_string = 'heatflux must not be set for pure neutral flow'
2046          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
2047       ENDIF
2048
2049    ENDIF
2050
2051    IF ( top_momentumflux_u /= 9999999.9_wp  .AND.                             &
2052         top_momentumflux_v /= 9999999.9_wp )  THEN
2053       constant_top_momentumflux = .TRUE.
2054    ELSEIF (  .NOT. ( top_momentumflux_u == 9999999.9_wp  .AND.                &
2055           top_momentumflux_v == 9999999.9_wp ) )  THEN
2056       message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' //  &
2057                        'must be set'
2058       CALL message( 'check_parameters', 'PA0064', 1, 2, 0, 6, 0 )
2059    ENDIF
2060
2061!
2062!-- A given surface temperature implies Dirichlet boundary condition for
2063!-- temperature. In this case specification of a constant heat flux is
2064!-- forbidden.
2065    IF ( ibc_pt_b == 0  .AND.  constant_heatflux  .AND.                        &
2066         surface_heatflux /= 0.0_wp )  THEN
2067       message_string = 'boundary_condition: bc_pt_b = "' // TRIM( bc_pt_b ) //&
2068                        '& is not allowed with constant_heatflux = .TRUE.'
2069       CALL message( 'check_parameters', 'PA0065', 1, 2, 0, 6, 0 )
2070    ENDIF
2071    IF ( constant_heatflux  .AND.  pt_surface_initial_change /= 0.0_wp )  THEN
2072       WRITE ( message_string, * )  'constant_heatflux = .TRUE. is not allo',  &
2073               'wed with pt_surface_initial_change (/=0) = ',                  &
2074               pt_surface_initial_change
2075       CALL message( 'check_parameters', 'PA0066', 1, 2, 0, 6, 0 )
2076    ENDIF
2077
2078!
2079!-- A given temperature at the top implies Dirichlet boundary condition for
2080!-- temperature. In this case specification of a constant heat flux is
2081!-- forbidden.
2082    IF ( ibc_pt_t == 0  .AND.  constant_top_heatflux  .AND.                    &
2083         top_heatflux /= 0.0_wp )  THEN
2084       message_string = 'boundary_condition: bc_pt_t = "' // TRIM( bc_pt_t ) //&
2085                        '" is not allowed with constant_top_heatflux = .TRUE.'
2086       CALL message( 'check_parameters', 'PA0067', 1, 2, 0, 6, 0 )
2087    ENDIF
2088
2089!
2090!-- Boundary conditions for salinity
2091    IF ( ocean )  THEN
2092       IF ( bc_sa_t == 'dirichlet' )  THEN
2093          ibc_sa_t = 0
2094       ELSEIF ( bc_sa_t == 'neumann' )  THEN
2095          ibc_sa_t = 1
2096       ELSE
2097          message_string = 'unknown boundary condition: bc_sa_t = "' //        &
2098                           TRIM( bc_sa_t ) // '"'
2099          CALL message( 'check_parameters', 'PA0068', 1, 2, 0, 6, 0 )
2100       ENDIF
2101
2102       IF ( top_salinityflux == 9999999.9_wp )  constant_top_salinityflux = .FALSE.
2103       IF ( ibc_sa_t == 1  .AND.  top_salinityflux == 9999999.9_wp )  THEN
2104          message_string = 'boundary condition: bc_sa_t = "' //                &
2105                           TRIM( bc_sa_t ) // '" requires to set ' //          &
2106                           'top_salinityflux'
2107          CALL message( 'check_parameters', 'PA0069', 1, 2, 0, 6, 0 )
2108       ENDIF
2109
2110!
2111!--    A fixed salinity at the top implies Dirichlet boundary condition for
2112!--    salinity. In this case specification of a constant salinity flux is
2113!--    forbidden.
2114       IF ( ibc_sa_t == 0  .AND.  constant_top_salinityflux  .AND.             &
2115            top_salinityflux /= 0.0_wp )  THEN
2116          message_string = 'boundary condition: bc_sa_t = "' //                &
2117                           TRIM( bc_sa_t ) // '" is not allowed with ' //      &
2118                           'top_salinityflux /= 0.0'
2119          CALL message( 'check_parameters', 'PA0070', 1, 2, 0, 6, 0 )
2120       ENDIF
2121
2122    ENDIF
2123
2124!
2125!-- Set boundary conditions for total water content
2126    IF ( humidity )  THEN
2127
2128       IF ( ANY( wall_humidityflux /= 0.0_wp )  .AND.                        &
2129            surface_waterflux == 9999999.9_wp )  THEN
2130          message_string = 'wall_humidityflux additionally requires ' //     &
2131                           'setting of surface_waterflux'
2132          CALL message( 'check_parameters', 'PA0444', 1, 2, 0, 6, 0 )
2133       ENDIF
2134
2135       CALL set_bc_scalars( 'q', bc_q_b, bc_q_t, ibc_q_b, ibc_q_t,           &
2136                            'PA0071', 'PA0072' )
2137
2138       IF ( surface_waterflux == 9999999.9_wp  )  THEN
2139          constant_waterflux = .FALSE.
2140          IF ( large_scale_forcing .OR. land_surface )  THEN
2141             IF ( ibc_q_b == 0 )  THEN
2142                constant_waterflux = .FALSE.
2143             ELSEIF ( ibc_q_b == 1 )  THEN
2144                constant_waterflux = .TRUE.
2145             ENDIF
2146          ENDIF
2147       ELSE
2148          constant_waterflux = .TRUE.
2149       ENDIF
2150
2151       CALL check_bc_scalars( 'q', bc_q_b, ibc_q_b, 'PA0073', 'PA0074',        &
2152                              constant_waterflux, q_surface_initial_change )
2153
2154    ENDIF
2155
2156    IF ( passive_scalar )  THEN
2157
2158       IF ( ANY( wall_scalarflux /= 0.0_wp )  .AND.                            &
2159            surface_scalarflux == 9999999.9_wp )  THEN
2160          message_string = 'wall_scalarflux additionally requires ' //         &
2161                           'setting of surface_scalarflux'
2162          CALL message( 'check_parameters', 'PA0445', 1, 2, 0, 6, 0 )
2163       ENDIF
2164
2165       IF ( surface_scalarflux == 9999999.9_wp )  constant_scalarflux = .FALSE.
2166
2167       CALL set_bc_scalars( 's', bc_s_b, bc_s_t, ibc_s_b, ibc_s_t,             &
2168                            'PA0071', 'PA0072' )
2169
2170       CALL check_bc_scalars( 's', bc_s_b, ibc_s_b, 'PA0073', 'PA0074',        &
2171                              constant_scalarflux, s_surface_initial_change )
2172
2173       IF ( top_scalarflux == 9999999.9_wp )  constant_top_scalarflux = .FALSE.
2174!
2175!--    A fixed scalar concentration at the top implies Dirichlet boundary
2176!--    condition for scalar. Hence, in this case specification of a constant
2177!--    scalar flux is forbidden.
2178       IF ( ( ibc_s_t == 0 .OR. ibc_s_t == 2 )  .AND.  constant_top_scalarflux &
2179               .AND.  top_scalarflux /= 0.0_wp )  THEN
2180          message_string = 'boundary condition: bc_s_t = "' //                 &
2181                           TRIM( bc_sa_t ) // '" is not allowed with ' //      &
2182                           'top_scalarflux /= 0.0'
2183          CALL message( 'check_parameters', 'PA0441', 1, 2, 0, 6, 0 )
2184       ENDIF
2185    ENDIF
2186
2187!
2188!-- Boundary conditions for chemical species
2189    IF ( air_chemistry )  CALL chem_boundary_conds( 'init' )
2190
2191!
2192!-- Boundary conditions for horizontal components of wind speed
2193    IF ( bc_uv_b == 'dirichlet' )  THEN
2194       ibc_uv_b = 0
2195    ELSEIF ( bc_uv_b == 'neumann' )  THEN
2196       ibc_uv_b = 1
2197       IF ( constant_flux_layer )  THEN
2198          message_string = 'boundary condition: bc_uv_b = "' //                &
2199               TRIM( bc_uv_b ) // '" is not allowed with constant_flux_layer'  &
2200               // ' = .TRUE.'
2201          CALL message( 'check_parameters', 'PA0075', 1, 2, 0, 6, 0 )
2202       ENDIF
2203    ELSE
2204       message_string = 'unknown boundary condition: bc_uv_b = "' //           &
2205                        TRIM( bc_uv_b ) // '"'
2206       CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 )
2207    ENDIF
2208!
2209!-- In case of coupled simulations u and v at the ground in atmosphere will be
2210!-- assigned with the u and v values of the ocean surface
2211    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
2212       ibc_uv_b = 2
2213    ENDIF
2214
2215    IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
2216       bc_uv_t = 'neumann'
2217       ibc_uv_t = 1
2218    ELSE
2219       IF ( bc_uv_t == 'dirichlet' .OR. bc_uv_t == 'dirichlet_0' )  THEN
2220          ibc_uv_t = 0
2221          IF ( bc_uv_t == 'dirichlet_0' )  THEN
2222!
2223!--          Velocities for the initial u,v-profiles are set zero at the top
2224!--          in case of dirichlet_0 conditions
2225             u_init(nzt+1)    = 0.0_wp
2226             v_init(nzt+1)    = 0.0_wp
2227          ENDIF
2228       ELSEIF ( bc_uv_t == 'neumann' )  THEN
2229          ibc_uv_t = 1
2230       ELSEIF ( bc_uv_t == 'nested'  .OR.  bc_uv_t == 'forcing' )  THEN
2231          ibc_uv_t = 3
2232       ELSE
2233          message_string = 'unknown boundary condition: bc_uv_t = "' //        &
2234                           TRIM( bc_uv_t ) // '"'
2235          CALL message( 'check_parameters', 'PA0077', 1, 2, 0, 6, 0 )
2236       ENDIF
2237    ENDIF
2238
2239!
2240!-- Compute and check, respectively, the Rayleigh Damping parameter
2241    IF ( rayleigh_damping_factor == -1.0_wp )  THEN
2242       rayleigh_damping_factor = 0.0_wp
2243    ELSE
2244       IF ( rayleigh_damping_factor < 0.0_wp  .OR.                             &
2245            rayleigh_damping_factor > 1.0_wp )  THEN
2246          WRITE( message_string, * )  'rayleigh_damping_factor = ',            &
2247                              rayleigh_damping_factor, ' out of range [0.0,1.0]'
2248          CALL message( 'check_parameters', 'PA0078', 1, 2, 0, 6, 0 )
2249       ENDIF
2250    ENDIF
2251
2252    IF ( rayleigh_damping_height == -1.0_wp )  THEN
2253       IF (  .NOT.  ocean )  THEN
2254          rayleigh_damping_height = 0.66666666666_wp * zu(nzt)
2255       ELSE
2256          rayleigh_damping_height = 0.66666666666_wp * zu(nzb)
2257       ENDIF
2258    ELSE
2259       IF (  .NOT.  ocean )  THEN
2260          IF ( rayleigh_damping_height < 0.0_wp  .OR.                          &
2261               rayleigh_damping_height > zu(nzt) )  THEN
2262             WRITE( message_string, * )  'rayleigh_damping_height = ',         &
2263                   rayleigh_damping_height, ' out of range [0.0,', zu(nzt), ']'
2264             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
2265          ENDIF
2266       ELSE
2267          IF ( rayleigh_damping_height > 0.0_wp  .OR.                          &
2268               rayleigh_damping_height < zu(nzb) )  THEN
2269             WRITE( message_string, * )  'rayleigh_damping_height = ',         &
2270                   rayleigh_damping_height, ' out of range [0.0,', zu(nzb), ']'
2271             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
2272          ENDIF
2273       ENDIF
2274    ENDIF
2275
2276!
2277!-- Check number of chosen statistic regions
2278    IF ( statistic_regions < 0 )  THEN
2279       WRITE ( message_string, * ) 'number of statistic_regions = ',           &
2280                   statistic_regions+1, ' is not allowed'
2281       CALL message( 'check_parameters', 'PA0082', 1, 2, 0, 6, 0 )
2282    ENDIF
2283    IF ( normalizing_region > statistic_regions  .OR.                          &
2284         normalizing_region < 0)  THEN
2285       WRITE ( message_string, * ) 'normalizing_region = ',                    &
2286                normalizing_region, ' must be >= 0 and <= ',statistic_regions, &
2287                ' (value of statistic_regions)'
2288       CALL message( 'check_parameters', 'PA0083', 1, 2, 0, 6, 0 )
2289    ENDIF
2290
2291!
2292!-- Set the default intervals for data output, if necessary
2293!-- NOTE: dt_dosp has already been set in spectra_parin
2294    IF ( dt_data_output /= 9999999.9_wp )  THEN
2295       IF ( dt_dopr           == 9999999.9_wp )  dt_dopr           = dt_data_output
2296       IF ( dt_dopts          == 9999999.9_wp )  dt_dopts          = dt_data_output
2297       IF ( dt_do2d_xy        == 9999999.9_wp )  dt_do2d_xy        = dt_data_output
2298       IF ( dt_do2d_xz        == 9999999.9_wp )  dt_do2d_xz        = dt_data_output
2299       IF ( dt_do2d_yz        == 9999999.9_wp )  dt_do2d_yz        = dt_data_output
2300       IF ( dt_do3d           == 9999999.9_wp )  dt_do3d           = dt_data_output
2301       IF ( dt_data_output_av == 9999999.9_wp )  dt_data_output_av = dt_data_output
2302       DO  mid = 1, max_masks
2303          IF ( dt_domask(mid) == 9999999.9_wp )  dt_domask(mid)    = dt_data_output
2304       ENDDO
2305    ENDIF
2306
2307!
2308!-- Set the default skip time intervals for data output, if necessary
2309    IF ( skip_time_dopr    == 9999999.9_wp )                                   &
2310                                       skip_time_dopr    = skip_time_data_output
2311    IF ( skip_time_do2d_xy == 9999999.9_wp )                                   &
2312                                       skip_time_do2d_xy = skip_time_data_output
2313    IF ( skip_time_do2d_xz == 9999999.9_wp )                                   &
2314                                       skip_time_do2d_xz = skip_time_data_output
2315    IF ( skip_time_do2d_yz == 9999999.9_wp )                                   &
2316                                       skip_time_do2d_yz = skip_time_data_output
2317    IF ( skip_time_do3d    == 9999999.9_wp )                                   &
2318                                       skip_time_do3d    = skip_time_data_output
2319    IF ( skip_time_data_output_av == 9999999.9_wp )                            &
2320                                skip_time_data_output_av = skip_time_data_output
2321    DO  mid = 1, max_masks
2322       IF ( skip_time_domask(mid) == 9999999.9_wp )                            &
2323                                skip_time_domask(mid)    = skip_time_data_output
2324    ENDDO
2325
2326!
2327!-- Check the average intervals (first for 3d-data, then for profiles)
2328    IF ( averaging_interval > dt_data_output_av )  THEN
2329       WRITE( message_string, * )  'averaging_interval = ',                    &
2330             averaging_interval, ' must be <= dt_data_output_av = ',           &
2331             dt_data_output_av
2332       CALL message( 'check_parameters', 'PA0085', 1, 2, 0, 6, 0 )
2333    ENDIF
2334
2335    IF ( averaging_interval_pr == 9999999.9_wp )  THEN
2336       averaging_interval_pr = averaging_interval
2337    ENDIF
2338
2339    IF ( averaging_interval_pr > dt_dopr )  THEN
2340       WRITE( message_string, * )  'averaging_interval_pr = ',                 &
2341             averaging_interval_pr, ' must be <= dt_dopr = ', dt_dopr
2342       CALL message( 'check_parameters', 'PA0086', 1, 2, 0, 6, 0 )
2343    ENDIF
2344
2345!
2346!-- Set the default interval for profiles entering the temporal average
2347    IF ( dt_averaging_input_pr == 9999999.9_wp )  THEN
2348       dt_averaging_input_pr = dt_averaging_input
2349    ENDIF
2350
2351!
2352!-- Set the default interval for the output of timeseries to a reasonable
2353!-- value (tries to minimize the number of calls of flow_statistics)
2354    IF ( dt_dots == 9999999.9_wp )  THEN
2355       IF ( averaging_interval_pr == 0.0_wp )  THEN
2356          dt_dots = MIN( dt_run_control, dt_dopr )
2357       ELSE
2358          dt_dots = MIN( dt_run_control, dt_averaging_input_pr )
2359       ENDIF
2360    ENDIF
2361
2362!
2363!-- Check the sample rate for averaging (first for 3d-data, then for profiles)
2364    IF ( dt_averaging_input > averaging_interval )  THEN
2365       WRITE( message_string, * )  'dt_averaging_input = ',                    &
2366                dt_averaging_input, ' must be <= averaging_interval = ',       &
2367                averaging_interval
2368       CALL message( 'check_parameters', 'PA0088', 1, 2, 0, 6, 0 )
2369    ENDIF
2370
2371    IF ( dt_averaging_input_pr > averaging_interval_pr )  THEN
2372       WRITE( message_string, * )  'dt_averaging_input_pr = ',                 &
2373                dt_averaging_input_pr, ' must be <= averaging_interval_pr = ', &
2374                averaging_interval_pr
2375       CALL message( 'check_parameters', 'PA0089', 1, 2, 0, 6, 0 )
2376    ENDIF
2377
2378!
2379!-- Set the default value for the integration interval of precipitation amount
2380    IF ( microphysics_seifert  .OR.  microphysics_kessler )  THEN
2381       IF ( precipitation_amount_interval == 9999999.9_wp )  THEN
2382          precipitation_amount_interval = dt_do2d_xy
2383       ELSE
2384          IF ( precipitation_amount_interval > dt_do2d_xy )  THEN
2385             WRITE( message_string, * )  'precipitation_amount_interval = ',   &
2386                 precipitation_amount_interval, ' must not be larger than ',   &
2387                 'dt_do2d_xy = ', dt_do2d_xy
2388             CALL message( 'check_parameters', 'PA0090', 1, 2, 0, 6, 0 )
2389          ENDIF
2390       ENDIF
2391    ENDIF
2392
2393!
2394!-- Determine the number of output profiles and check whether they are
2395!-- permissible
2396    DO  WHILE ( data_output_pr(dopr_n+1) /= '          ' )
2397
2398       dopr_n = dopr_n + 1
2399       i = dopr_n
2400
2401!
2402!--    Determine internal profile number (for hom, homs)
2403!--    and store height levels
2404       SELECT CASE ( TRIM( data_output_pr(i) ) )
2405
2406          CASE ( 'u', '#u' )
2407             dopr_index(i) = 1
2408             dopr_unit(i)  = 'm/s'
2409             hom(:,2,1,:)  = SPREAD( zu, 2, statistic_regions+1 )
2410             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2411                dopr_initial_index(i) = 5
2412                hom(:,2,5,:)          = SPREAD( zu, 2, statistic_regions+1 )
2413                data_output_pr(i)     = data_output_pr(i)(2:)
2414             ENDIF
2415
2416          CASE ( 'v', '#v' )
2417             dopr_index(i) = 2
2418             dopr_unit(i)  = 'm/s'
2419             hom(:,2,2,:)  = SPREAD( zu, 2, statistic_regions+1 )
2420             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2421                dopr_initial_index(i) = 6
2422                hom(:,2,6,:)          = SPREAD( zu, 2, statistic_regions+1 )
2423                data_output_pr(i)     = data_output_pr(i)(2:)
2424             ENDIF
2425
2426          CASE ( 'w' )
2427             dopr_index(i) = 3
2428             dopr_unit(i)  = 'm/s'
2429             hom(:,2,3,:)  = SPREAD( zw, 2, statistic_regions+1 )
2430
2431          CASE ( 'pt', '#pt' )
2432             IF ( .NOT. cloud_physics ) THEN
2433                dopr_index(i) = 4
2434                dopr_unit(i)  = 'K'
2435                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2436                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2437                   dopr_initial_index(i) = 7
2438                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2439                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
2440                   data_output_pr(i)     = data_output_pr(i)(2:)
2441                ENDIF
2442             ELSE
2443                dopr_index(i) = 43
2444                dopr_unit(i)  = 'K'
2445                hom(:,2,43,:)  = SPREAD( zu, 2, statistic_regions+1 )
2446                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2447                   dopr_initial_index(i) = 28
2448                   hom(:,2,28,:)         = SPREAD( zu, 2, statistic_regions+1 )
2449                   hom(nzb,2,28,:)       = 0.0_wp    ! because zu(nzb) is negative
2450                   data_output_pr(i)     = data_output_pr(i)(2:)
2451                ENDIF
2452             ENDIF
2453
2454          CASE ( 'e' )
2455             dopr_index(i)  = 8
2456             dopr_unit(i)   = 'm2/s2'
2457             hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
2458             hom(nzb,2,8,:) = 0.0_wp
2459
2460          CASE ( 'km', '#km' )
2461             dopr_index(i)  = 9
2462             dopr_unit(i)   = 'm2/s'
2463             hom(:,2,9,:)   = SPREAD( zu, 2, statistic_regions+1 )
2464             hom(nzb,2,9,:) = 0.0_wp
2465             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2466                dopr_initial_index(i) = 23
2467                hom(:,2,23,:)         = hom(:,2,9,:)
2468                data_output_pr(i)     = data_output_pr(i)(2:)
2469             ENDIF
2470
2471          CASE ( 'kh', '#kh' )
2472             dopr_index(i)   = 10
2473             dopr_unit(i)    = 'm2/s'
2474             hom(:,2,10,:)   = SPREAD( zu, 2, statistic_regions+1 )
2475             hom(nzb,2,10,:) = 0.0_wp
2476             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2477                dopr_initial_index(i) = 24
2478                hom(:,2,24,:)         = hom(:,2,10,:)
2479                data_output_pr(i)     = data_output_pr(i)(2:)
2480             ENDIF
2481
2482          CASE ( 'l', '#l' )
2483             dopr_index(i)   = 11
2484             dopr_unit(i)    = 'm'
2485             hom(:,2,11,:)   = SPREAD( zu, 2, statistic_regions+1 )
2486             hom(nzb,2,11,:) = 0.0_wp
2487             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2488                dopr_initial_index(i) = 25
2489                hom(:,2,25,:)         = hom(:,2,11,:)
2490                data_output_pr(i)     = data_output_pr(i)(2:)
2491             ENDIF
2492
2493          CASE ( 'w"u"' )
2494             dopr_index(i) = 12
2495             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2496             hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
2497             IF ( constant_flux_layer )  hom(nzb,2,12,:) = zu(1)
2498
2499          CASE ( 'w*u*' )
2500             dopr_index(i) = 13
2501             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2502             hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
2503
2504          CASE ( 'w"v"' )
2505             dopr_index(i) = 14
2506             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2507             hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
2508             IF ( constant_flux_layer )  hom(nzb,2,14,:) = zu(1)
2509
2510          CASE ( 'w*v*' )
2511             dopr_index(i) = 15
2512             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2513             hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
2514
2515          CASE ( 'w"pt"' )
2516             dopr_index(i) = 16
2517             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2518             hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
2519
2520          CASE ( 'w*pt*' )
2521             dopr_index(i) = 17
2522             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2523             hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
2524
2525          CASE ( 'wpt' )
2526             dopr_index(i) = 18
2527             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2528             hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
2529
2530          CASE ( 'wu' )
2531             dopr_index(i) = 19
2532             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2533             hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
2534             IF ( constant_flux_layer )  hom(nzb,2,19,:) = zu(1)
2535
2536          CASE ( 'wv' )
2537             dopr_index(i) = 20
2538             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2539             hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
2540             IF ( constant_flux_layer )  hom(nzb,2,20,:) = zu(1)
2541
2542          CASE ( 'w*pt*BC' )
2543             dopr_index(i) = 21
2544             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2545             hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
2546
2547          CASE ( 'wptBC' )
2548             dopr_index(i) = 22
2549             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2550             hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
2551
2552          CASE ( 'sa', '#sa' )
2553             IF ( .NOT. ocean )  THEN
2554                message_string = 'data_output_pr = ' //                        &
2555                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2556                                 'lemented for ocean = .FALSE.'
2557                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2558             ELSE
2559                dopr_index(i) = 23
2560                dopr_unit(i)  = 'psu'
2561                hom(:,2,23,:) = SPREAD( zu, 2, statistic_regions+1 )
2562                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2563                   dopr_initial_index(i) = 26
2564                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2565                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2566                   data_output_pr(i)     = data_output_pr(i)(2:)
2567                ENDIF
2568             ENDIF
2569
2570          CASE ( 'u*2' )
2571             dopr_index(i) = 30
2572             dopr_unit(i)  = 'm2/s2'
2573             hom(:,2,30,:) = SPREAD( zu, 2, statistic_regions+1 )
2574
2575          CASE ( 'v*2' )
2576             dopr_index(i) = 31
2577             dopr_unit(i)  = 'm2/s2'
2578             hom(:,2,31,:) = SPREAD( zu, 2, statistic_regions+1 )
2579
2580          CASE ( 'w*2' )
2581             dopr_index(i) = 32
2582             dopr_unit(i)  = 'm2/s2'
2583             hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
2584
2585          CASE ( 'pt*2' )
2586             dopr_index(i) = 33
2587             dopr_unit(i)  = 'K2'
2588             hom(:,2,33,:) = SPREAD( zu, 2, statistic_regions+1 )
2589
2590          CASE ( 'e*' )
2591             dopr_index(i) = 34
2592             dopr_unit(i)  = 'm2/s2'
2593             hom(:,2,34,:) = SPREAD( zu, 2, statistic_regions+1 )
2594
2595          CASE ( 'w*2pt*' )
2596             dopr_index(i) = 35
2597             dopr_unit(i)  = 'K m2/s2'
2598             hom(:,2,35,:) = SPREAD( zw, 2, statistic_regions+1 )
2599
2600          CASE ( 'w*pt*2' )
2601             dopr_index(i) = 36
2602             dopr_unit(i)  = 'K2 m/s'
2603             hom(:,2,36,:) = SPREAD( zw, 2, statistic_regions+1 )
2604
2605          CASE ( 'w*e*' )
2606             dopr_index(i) = 37
2607             dopr_unit(i)  = 'm3/s3'
2608             hom(:,2,37,:) = SPREAD( zw, 2, statistic_regions+1 )
2609
2610          CASE ( 'w*3' )
2611             dopr_index(i) = 38
2612             dopr_unit(i)  = 'm3/s3'
2613             hom(:,2,38,:) = SPREAD( zw, 2, statistic_regions+1 )
2614
2615          CASE ( 'Sw' )
2616             dopr_index(i) = 39
2617             dopr_unit(i)  = 'none'
2618             hom(:,2,39,:) = SPREAD( zw, 2, statistic_regions+1 )
2619
2620          CASE ( 'p' )
2621             dopr_index(i) = 40
2622             dopr_unit(i)  = 'Pa'
2623             hom(:,2,40,:) = SPREAD( zu, 2, statistic_regions+1 )
2624
2625          CASE ( 'q', '#q' )
2626             IF ( .NOT. humidity )  THEN
2627                message_string = 'data_output_pr = ' //                        &
2628                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2629                                 'lemented for humidity = .FALSE.'
2630                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2631             ELSE
2632                dopr_index(i) = 41
2633                dopr_unit(i)  = 'kg/kg'
2634                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2635                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2636                   dopr_initial_index(i) = 26
2637                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2638                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2639                   data_output_pr(i)     = data_output_pr(i)(2:)
2640                ENDIF
2641             ENDIF
2642
2643          CASE ( 's', '#s' )
2644             IF ( .NOT. passive_scalar )  THEN
2645                message_string = 'data_output_pr = ' //                        &
2646                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2647                                 'lemented for passive_scalar = .FALSE.'
2648                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2649             ELSE
2650                dopr_index(i) = 115
2651                dopr_unit(i)  = 'kg/m3'
2652                hom(:,2,115,:) = SPREAD( zu, 2, statistic_regions+1 )
2653                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2654                   dopr_initial_index(i) = 121
2655                   hom(:,2,121,:)        = SPREAD( zu, 2, statistic_regions+1 )
2656                   hom(nzb,2,121,:)      = 0.0_wp    ! because zu(nzb) is negative
2657                   data_output_pr(i)     = data_output_pr(i)(2:)
2658                ENDIF
2659             ENDIF
2660
2661          CASE ( 'qv', '#qv' )
2662             IF ( .NOT. cloud_physics ) THEN
2663                dopr_index(i) = 41
2664                dopr_unit(i)  = 'kg/kg'
2665                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2666                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2667                   dopr_initial_index(i) = 26
2668                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2669                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2670                   data_output_pr(i)     = data_output_pr(i)(2:)
2671                ENDIF
2672             ELSE
2673                dopr_index(i) = 42
2674                dopr_unit(i)  = 'kg/kg'
2675                hom(:,2,42,:) = SPREAD( zu, 2, statistic_regions+1 )
2676                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2677                   dopr_initial_index(i) = 27
2678                   hom(:,2,27,:)         = SPREAD( zu, 2, statistic_regions+1 )
2679                   hom(nzb,2,27,:)       = 0.0_wp   ! because zu(nzb) is negative
2680                   data_output_pr(i)     = data_output_pr(i)(2:)
2681                ENDIF
2682             ENDIF
2683
2684          CASE ( 'lpt', '#lpt' )
2685             IF ( .NOT. cloud_physics ) THEN
2686                message_string = 'data_output_pr = ' //                        &
2687                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2688                                 'lemented for cloud_physics = .FALSE.'
2689                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2690             ELSE
2691                dopr_index(i) = 4
2692                dopr_unit(i)  = 'K'
2693                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2694                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2695                   dopr_initial_index(i) = 7
2696                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2697                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
2698                   data_output_pr(i)     = data_output_pr(i)(2:)
2699                ENDIF
2700             ENDIF
2701
2702          CASE ( 'vpt', '#vpt' )
2703             dopr_index(i) = 44
2704             dopr_unit(i)  = 'K'
2705             hom(:,2,44,:) = SPREAD( zu, 2, statistic_regions+1 )
2706             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2707                dopr_initial_index(i) = 29
2708                hom(:,2,29,:)         = SPREAD( zu, 2, statistic_regions+1 )
2709                hom(nzb,2,29,:)       = 0.0_wp    ! because zu(nzb) is negative
2710                data_output_pr(i)     = data_output_pr(i)(2:)
2711             ENDIF
2712
2713          CASE ( 'w"vpt"' )
2714             dopr_index(i) = 45
2715             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2716             hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
2717
2718          CASE ( 'w*vpt*' )
2719             dopr_index(i) = 46
2720             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2721             hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
2722
2723          CASE ( 'wvpt' )
2724             dopr_index(i) = 47
2725             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2726             hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
2727
2728          CASE ( 'w"q"' )
2729             IF (  .NOT.  humidity )  THEN
2730                message_string = 'data_output_pr = ' //                        &
2731                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2732                                 'lemented for humidity = .FALSE.'
2733                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2734             ELSE
2735                dopr_index(i) = 48
2736                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2737                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2738             ENDIF
2739
2740          CASE ( 'w*q*' )
2741             IF (  .NOT.  humidity )  THEN
2742                message_string = 'data_output_pr = ' //                        &
2743                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2744                                 'lemented for humidity = .FALSE.'
2745                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2746             ELSE
2747                dopr_index(i) = 49
2748                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2749                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2750             ENDIF
2751
2752          CASE ( 'wq' )
2753             IF (  .NOT.  humidity )  THEN
2754                message_string = 'data_output_pr = ' //                        &
2755                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2756                                 'lemented for humidity = .FALSE.'
2757                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2758             ELSE
2759                dopr_index(i) = 50
2760                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2761                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2762             ENDIF
2763
2764          CASE ( 'w"s"' )
2765             IF (  .NOT.  passive_scalar )  THEN
2766                message_string = 'data_output_pr = ' //                        &
2767                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2768                                 'lemented for passive_scalar = .FALSE.'
2769                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2770             ELSE
2771                dopr_index(i) = 117
2772                dopr_unit(i)  = 'kg/m3 m/s'
2773                hom(:,2,117,:) = SPREAD( zw, 2, statistic_regions+1 )
2774             ENDIF
2775
2776          CASE ( 'w*s*' )
2777             IF (  .NOT.  passive_scalar )  THEN
2778                message_string = 'data_output_pr = ' //                        &
2779                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2780                                 'lemented for passive_scalar = .FALSE.'
2781                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2782             ELSE
2783                dopr_index(i) = 114
2784                dopr_unit(i)  = 'kg/m3 m/s'
2785                hom(:,2,114,:) = SPREAD( zw, 2, statistic_regions+1 )
2786             ENDIF
2787
2788          CASE ( 'ws' )
2789             IF (  .NOT.  passive_scalar )  THEN
2790                message_string = 'data_output_pr = ' //                        &
2791                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2792                                 'lemented for passive_scalar = .FALSE.'
2793                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2794             ELSE
2795                dopr_index(i) = 118
2796                dopr_unit(i)  = 'kg/m3 m/s'
2797                hom(:,2,118,:) = SPREAD( zw, 2, statistic_regions+1 )
2798             ENDIF
2799
2800          CASE ( 'w"qv"' )
2801             IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
2802                dopr_index(i) = 48
2803                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2804                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2805             ELSEIF ( humidity  .AND.  cloud_physics )  THEN
2806                dopr_index(i) = 51
2807                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2808                hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
2809             ELSE
2810                message_string = 'data_output_pr = ' //                        &
2811                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2812                                 'lemented for cloud_physics = .FALSE. an'  // &
2813                                 'd humidity = .FALSE.'
2814                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2815             ENDIF
2816
2817          CASE ( 'w*qv*' )
2818             IF ( humidity  .AND.  .NOT. cloud_physics )                       &
2819             THEN
2820                dopr_index(i) = 49
2821                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2822                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2823             ELSEIF( humidity .AND. cloud_physics ) THEN
2824                dopr_index(i) = 52
2825                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2826                hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
2827             ELSE
2828                message_string = 'data_output_pr = ' //                        &
2829                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2830                                 'lemented for cloud_physics = .FALSE. an'  // &
2831                                 'd humidity = .FALSE.'
2832                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2833             ENDIF
2834
2835          CASE ( 'wqv' )
2836             IF ( humidity  .AND.  .NOT.  cloud_physics )  THEN
2837                dopr_index(i) = 50
2838                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2839                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2840             ELSEIF ( humidity  .AND.  cloud_physics )  THEN
2841                dopr_index(i) = 53
2842                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2843                hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
2844             ELSE
2845                message_string = 'data_output_pr = ' //                        &
2846                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2847                                 'lemented for cloud_physics = .FALSE. an'  // &
2848                                 'd humidity = .FALSE.'
2849                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2850             ENDIF
2851
2852          CASE ( 'ql' )
2853             IF (  .NOT.  cloud_physics  .AND.  .NOT.  cloud_droplets )  THEN
2854                message_string = 'data_output_pr = ' //                        &
2855                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2856                                 'lemented for cloud_physics = .FALSE. and' // &
2857                                 'cloud_droplets = .FALSE.'
2858                CALL message( 'check_parameters', 'PA0096', 1, 2, 0, 6, 0 )
2859             ELSE
2860                dopr_index(i) = 54
2861                dopr_unit(i)  = 'kg/kg'
2862                hom(:,2,54,:)  = SPREAD( zu, 2, statistic_regions+1 )
2863             ENDIF
2864
2865          CASE ( 'w*u*u*:dz' )
2866             dopr_index(i) = 55
2867             dopr_unit(i)  = 'm2/s3'
2868             hom(:,2,55,:) = SPREAD( zu, 2, statistic_regions+1 )
2869
2870          CASE ( 'w*p*:dz' )
2871             dopr_index(i) = 56
2872             dopr_unit(i)  = 'm2/s3'
2873             hom(:,2,56,:) = SPREAD( zw, 2, statistic_regions+1 )
2874
2875          CASE ( 'w"e:dz' )
2876             dopr_index(i) = 57
2877             dopr_unit(i)  = 'm2/s3'
2878             hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
2879
2880
2881          CASE ( 'u"pt"' )
2882             dopr_index(i) = 58
2883             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2884             hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
2885
2886          CASE ( 'u*pt*' )
2887             dopr_index(i) = 59
2888             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2889             hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
2890
2891          CASE ( 'upt_t' )
2892             dopr_index(i) = 60
2893             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2894             hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
2895
2896          CASE ( 'v"pt"' )
2897             dopr_index(i) = 61
2898             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2899             hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
2900
2901          CASE ( 'v*pt*' )
2902             dopr_index(i) = 62
2903             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2904             hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
2905
2906          CASE ( 'vpt_t' )
2907             dopr_index(i) = 63
2908             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2909             hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
2910
2911          CASE ( 'rho_ocean' )
2912             IF (  .NOT.  ocean ) THEN
2913                message_string = 'data_output_pr = ' //                        &
2914                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2915                                 'lemented for ocean = .FALSE.'
2916                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2917             ELSE
2918                dopr_index(i) = 64
2919                dopr_unit(i)  = 'kg/m3'
2920                hom(:,2,64,:) = SPREAD( zu, 2, statistic_regions+1 )
2921                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2922                   dopr_initial_index(i) = 77
2923                   hom(:,2,77,:)         = SPREAD( zu, 2, statistic_regions+1 )
2924                   hom(nzb,2,77,:)       = 0.0_wp    ! because zu(nzb) is negative
2925                   data_output_pr(i)     = data_output_pr(i)(2:)
2926                ENDIF
2927             ENDIF
2928
2929          CASE ( 'w"sa"' )
2930             IF (  .NOT.  ocean ) THEN
2931                message_string = 'data_output_pr = ' //                        &
2932                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2933                                 'lemented for ocean = .FALSE.'
2934                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2935             ELSE
2936                dopr_index(i) = 65
2937                dopr_unit(i)  = 'psu m/s'
2938                hom(:,2,65,:) = SPREAD( zw, 2, statistic_regions+1 )
2939             ENDIF
2940
2941          CASE ( 'w*sa*' )
2942             IF (  .NOT. ocean  ) THEN
2943                message_string = 'data_output_pr = ' //                        &
2944                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2945                                 'lemented for ocean = .FALSE.'
2946                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2947             ELSE
2948                dopr_index(i) = 66
2949                dopr_unit(i)  = 'psu m/s'
2950                hom(:,2,66,:) = SPREAD( zw, 2, statistic_regions+1 )
2951             ENDIF
2952
2953          CASE ( 'wsa' )
2954             IF (  .NOT.  ocean ) THEN
2955                message_string = 'data_output_pr = ' //                        &
2956                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2957                                 'lemented for ocean = .FALSE.'
2958                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2959             ELSE
2960                dopr_index(i) = 67
2961                dopr_unit(i)  = 'psu m/s'
2962                hom(:,2,67,:) = SPREAD( zw, 2, statistic_regions+1 )
2963             ENDIF
2964
2965          CASE ( 'w*p*' )
2966             dopr_index(i) = 68
2967             dopr_unit(i)  = 'm3/s3'
2968             hom(:,2,68,:) = SPREAD( zu, 2, statistic_regions+1 )
2969
2970          CASE ( 'w"e' )
2971             dopr_index(i) = 69
2972             dopr_unit(i)  = 'm3/s3'
2973             hom(:,2,69,:) = SPREAD( zu, 2, statistic_regions+1 )
2974
2975          CASE ( 'q*2' )
2976             IF (  .NOT.  humidity )  THEN
2977                message_string = 'data_output_pr = ' //                        &
2978                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2979                                 'lemented for humidity = .FALSE.'
2980                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2981             ELSE
2982                dopr_index(i) = 70
2983                dopr_unit(i)  = 'kg2/kg2'
2984                hom(:,2,70,:) = SPREAD( zu, 2, statistic_regions+1 )
2985             ENDIF
2986
2987          CASE ( 'prho' )
2988             IF (  .NOT.  ocean ) THEN
2989                message_string = 'data_output_pr = ' //                        &
2990                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2991                                 'lemented for ocean = .FALSE.'
2992                CALL message( 'check_parameters', 'PA0091', 1, 2, 0, 6, 0 )
2993             ELSE
2994                dopr_index(i) = 71
2995                dopr_unit(i)  = 'kg/m3'
2996                hom(:,2,71,:) = SPREAD( zu, 2, statistic_regions+1 )
2997             ENDIF
2998
2999          CASE ( 'hyp' )
3000             dopr_index(i) = 72
3001             dopr_unit(i)  = 'hPa'
3002             hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
3003
3004          CASE ( 'rho_air' )
3005             dopr_index(i)  = 119
3006             dopr_unit(i)   = 'kg/m3'
3007             hom(:,2,119,:) = SPREAD( zu, 2, statistic_regions+1 )
3008
3009          CASE ( 'rho_air_zw' )
3010             dopr_index(i)  = 120
3011             dopr_unit(i)   = 'kg/m3'
3012             hom(:,2,120,:) = SPREAD( zw, 2, statistic_regions+1 )
3013
3014          CASE ( 'nc' )
3015             IF (  .NOT.  cloud_physics )  THEN
3016                message_string = 'data_output_pr = ' //                        &
3017                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3018                                 'lemented for cloud_physics = .FALSE.'
3019                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
3020             ELSEIF ( .NOT.  microphysics_morrison )  THEN
3021                message_string = 'data_output_pr = ' //                        &
3022                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3023                                 'lemented for cloud_scheme /= morrison'
3024                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
3025             ELSE
3026                dopr_index(i) = 89
3027                dopr_unit(i)  = '1/m3'
3028                hom(:,2,89,:)  = SPREAD( zu, 2, statistic_regions+1 )
3029             ENDIF
3030
3031          CASE ( 'nr' )
3032             IF (  .NOT.  cloud_physics )  THEN
3033                message_string = 'data_output_pr = ' //                        &
3034                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3035                                 'lemented for cloud_physics = .FALSE.'
3036                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
3037             ELSEIF ( .NOT.  microphysics_seifert )  THEN
3038                message_string = 'data_output_pr = ' //                        &
3039                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3040                                 'lemented for cloud_scheme /= seifert_beheng'
3041                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
3042             ELSE
3043                dopr_index(i) = 73
3044                dopr_unit(i)  = '1/m3'
3045                hom(:,2,73,:)  = SPREAD( zu, 2, statistic_regions+1 )
3046             ENDIF
3047
3048          CASE ( 'qr' )
3049             IF (  .NOT.  cloud_physics )  THEN
3050                message_string = 'data_output_pr = ' //                        &
3051                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3052                                 'lemented for cloud_physics = .FALSE.'
3053                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
3054             ELSEIF ( .NOT.  microphysics_seifert )  THEN
3055                message_string = 'data_output_pr = ' //                        &
3056                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3057                                 'lemented for cloud_scheme /= seifert_beheng'
3058                CALL message( 'check_parameters', 'PA0358', 1, 2, 0, 6, 0 )
3059             ELSE
3060                dopr_index(i) = 74
3061                dopr_unit(i)  = 'kg/kg'
3062                hom(:,2,74,:)  = SPREAD( zu, 2, statistic_regions+1 )
3063             ENDIF
3064
3065          CASE ( 'qc' )
3066             IF (  .NOT.  cloud_physics )  THEN
3067                message_string = 'data_output_pr = ' //                        &
3068                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3069                                 'lemented for cloud_physics = .FALSE.'
3070                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
3071             ELSE
3072                dopr_index(i) = 75
3073                dopr_unit(i)  = 'kg/kg'
3074                hom(:,2,75,:)  = SPREAD( zu, 2, statistic_regions+1 )
3075             ENDIF
3076
3077          CASE ( 'prr' )
3078             IF (  .NOT.  cloud_physics )  THEN
3079                message_string = 'data_output_pr = ' //                        &
3080                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3081                                 'lemented for cloud_physics = .FALSE.'
3082                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
3083             ELSEIF ( microphysics_sat_adjust )  THEN
3084                message_string = 'data_output_pr = ' //                        &
3085                                 TRIM( data_output_pr(i) ) // ' is not ava' // &
3086                                 'ilable for cloud_scheme = saturation_adjust'
3087                CALL message( 'check_parameters', 'PA0422', 1, 2, 0, 6, 0 )
3088             ELSE
3089                dopr_index(i) = 76
3090                dopr_unit(i)  = 'kg/kg m/s'
3091                hom(:,2,76,:)  = SPREAD( zu, 2, statistic_regions+1 )
3092             ENDIF
3093
3094          CASE ( 'ug' )
3095             dopr_index(i) = 78
3096             dopr_unit(i)  = 'm/s'
3097             hom(:,2,78,:) = SPREAD( zu, 2, statistic_regions+1 )
3098
3099          CASE ( 'vg' )
3100             dopr_index(i) = 79
3101             dopr_unit(i)  = 'm/s'
3102             hom(:,2,79,:) = SPREAD( zu, 2, statistic_regions+1 )
3103
3104          CASE ( 'w_subs' )
3105             IF (  .NOT.  large_scale_subsidence )  THEN
3106                message_string = 'data_output_pr = ' //                        &
3107                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3108                                 'lemented for large_scale_subsidence = .FALSE.'
3109                CALL message( 'check_parameters', 'PA0382', 1, 2, 0, 6, 0 )
3110             ELSE
3111                dopr_index(i) = 80
3112                dopr_unit(i)  = 'm/s'
3113                hom(:,2,80,:) = SPREAD( zu, 2, statistic_regions+1 )
3114             ENDIF
3115
3116          CASE ( 's*2' )
3117             IF (  .NOT.  passive_scalar )  THEN
3118                message_string = 'data_output_pr = ' //                        &
3119                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
3120                                 'lemented for passive_scalar = .FALSE.'
3121                CALL message( 'check_parameters', 'PA0185', 1, 2, 0, 6, 0 )
3122             ELSE
3123                dopr_index(i) = 116
3124                dopr_unit(i)  = 'kg2/m6'
3125                hom(:,2,116,:) = SPREAD( zu, 2, statistic_regions+1 )
3126             ENDIF
3127
3128
3129
3130          CASE DEFAULT
3131
3132             CALL lsm_check_data_output_pr( data_output_pr(i), i, unit,        &
3133                                            dopr_unit(i) )
3134
3135             IF ( unit == 'illegal' )  THEN
3136                CALL lsf_nudging_check_data_output_pr( data_output_pr(i), i,   &
3137                                                    unit, dopr_unit(i) )
3138             ENDIF
3139
3140             IF ( unit == 'illegal' )  THEN
3141                CALL radiation_check_data_output_pr( data_output_pr(i), i,     &
3142                                                     unit, dopr_unit(i) )
3143             ENDIF
3144!
3145!--          Block of gust module profile outputs
3146             IF ( unit == 'illegal'  .AND.  gust_module_enabled  )  THEN
3147                   CALL gust_check_data_output_pr( data_output_pr(i), i, unit, &
3148                                                   dopr_unit(i) )
3149             ENDIF
3150
3151             IF ( unit == 'illegal' )  THEN                                       
3152                CALL chem_check_data_output_pr( data_output_pr(i), i,          &
3153                                                unit, dopr_unit(i) )
3154             ENDIF
3155
3156             IF ( unit == 'illegal' )  THEN
3157                unit = ''
3158                CALL user_check_data_output_pr( data_output_pr(i), i, unit )
3159             ENDIF
3160
3161             IF ( unit == 'illegal' )  THEN
3162                IF ( data_output_pr_user(1) /= ' ' )  THEN
3163                   message_string = 'illegal value for data_output_pr or ' //  &
3164                                    'data_output_pr_user = "' //               &
3165                                    TRIM( data_output_pr(i) ) // '"'
3166                   CALL message( 'check_parameters', 'PA0097', 1, 2, 0, 6, 0 )
3167                ELSE
3168                   message_string = 'illegal value for data_output_pr = "' //  &
3169                                    TRIM( data_output_pr(i) ) // '"'
3170                   CALL message( 'check_parameters', 'PA0098', 1, 2, 0, 6, 0 )
3171                ENDIF
3172             ENDIF
3173
3174       END SELECT
3175
3176    ENDDO
3177
3178
3179!
3180!-- Append user-defined data output variables to the standard data output
3181    IF ( data_output_user(1) /= ' ' )  THEN
3182       i = 1
3183       DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 500 )
3184          i = i + 1
3185       ENDDO
3186       j = 1
3187       DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 500 )
3188          IF ( i > 500 )  THEN
3189             message_string = 'number of output quantitities given by data' // &
3190                '_output and data_output_user exceeds the limit of 500'
3191             CALL message( 'check_parameters', 'PA0102', 1, 2, 0, 6, 0 )
3192          ENDIF
3193          data_output(i) = data_output_user(j)
3194          i = i + 1
3195          j = j + 1
3196       ENDDO
3197    ENDIF
3198
3199!
3200!-- Check and set steering parameters for 2d/3d data output and averaging
3201    i   = 1
3202    DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 500 )
3203!
3204!--    Check for data averaging
3205       ilen = LEN_TRIM( data_output(i) )
3206       j = 0                                                 ! no data averaging
3207       IF ( ilen > 3 )  THEN
3208          IF ( data_output(i)(ilen-2:ilen) == '_av' )  THEN
3209             j = 1                                           ! data averaging
3210             data_output(i) = data_output(i)(1:ilen-3)
3211          ENDIF
3212       ENDIF
3213!
3214!--    Check for cross section or volume data
3215       ilen = LEN_TRIM( data_output(i) )
3216       k = 0                                                   ! 3d data
3217       var = data_output(i)(1:ilen)
3218       IF ( ilen > 3 )  THEN
3219          IF ( data_output(i)(ilen-2:ilen) == '_xy'  .OR.                      &
3220               data_output(i)(ilen-2:ilen) == '_xz'  .OR.                      &
3221               data_output(i)(ilen-2:ilen) == '_yz' )  THEN
3222             k = 1                                             ! 2d data
3223             var = data_output(i)(1:ilen-3)
3224          ENDIF
3225       ENDIF
3226
3227!
3228!--    Check for allowed value and set units
3229       SELECT CASE ( TRIM( var ) )
3230
3231          CASE ( 'e' )
3232             IF ( constant_diffusion )  THEN
3233                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3234                                 'res constant_diffusion = .FALSE.'
3235                CALL message( 'check_parameters', 'PA0103', 1, 2, 0, 6, 0 )
3236             ENDIF
3237             unit = 'm2/s2'
3238
3239          CASE ( 'lpt' )
3240             IF (  .NOT.  cloud_physics )  THEN
3241                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3242                         'res cloud_physics = .TRUE.'
3243                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3244             ENDIF
3245             unit = 'K'
3246
3247          CASE ( 'nc' )
3248             IF (  .NOT.  cloud_physics )  THEN
3249                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3250                         'res cloud_physics = .TRUE.'
3251                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3252             ELSEIF ( .NOT.  microphysics_morrison )  THEN
3253                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3254                         'res = morrison '
3255                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3256             ENDIF
3257             unit = '1/m3'
3258
3259          CASE ( 'nr' )
3260             IF (  .NOT.  cloud_physics )  THEN
3261                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3262                         'res cloud_physics = .TRUE.'
3263                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3264             ELSEIF ( .NOT.  microphysics_seifert )  THEN
3265                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3266                         'res cloud_scheme = seifert_beheng'
3267                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3268             ENDIF
3269             unit = '1/m3'
3270
3271          CASE ( 'pc', 'pr' )
3272             IF (  .NOT.  particle_advection )  THEN
3273                message_string = 'output of "' // TRIM( var ) // '" requir' // &
3274                   'es a "particle_parameters"-NAMELIST in the parameter ' //  &
3275                   'file (PARIN)'
3276                CALL message( 'check_parameters', 'PA0104', 1, 2, 0, 6, 0 )
3277             ENDIF
3278             IF ( TRIM( var ) == 'pc' )  unit = 'number'
3279             IF ( TRIM( var ) == 'pr' )  unit = 'm'
3280
3281          CASE ( 'prr' )
3282             IF (  .NOT.  cloud_physics )  THEN
3283                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3284                         'res cloud_physics = .TRUE.'
3285                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3286             ELSEIF ( microphysics_sat_adjust )  THEN
3287                message_string = 'output of "' // TRIM( var ) // '" is ' //    &
3288                         'not available for cloud_scheme = saturation_adjust'
3289                CALL message( 'check_parameters', 'PA0423', 1, 2, 0, 6, 0 )
3290             ENDIF
3291             unit = 'kg/kg m/s'
3292
3293          CASE ( 'q', 'vpt' )
3294             IF (  .NOT.  humidity )  THEN
3295                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3296                                 'res humidity = .TRUE.'
3297                CALL message( 'check_parameters', 'PA0105', 1, 2, 0, 6, 0 )
3298             ENDIF
3299             IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
3300             IF ( TRIM( var ) == 'vpt' )  unit = 'K'
3301
3302          CASE ( 'qc' )
3303             IF (  .NOT.  cloud_physics )  THEN
3304                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3305                         'res cloud_physics = .TRUE.'
3306                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3307             ENDIF
3308             unit = 'kg/kg'
3309
3310          CASE ( 'ql' )
3311             IF ( .NOT.  ( cloud_physics  .OR.  cloud_droplets ) )  THEN
3312                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3313                         'res cloud_physics = .TRUE. or cloud_droplets = .TRUE.'
3314                CALL message( 'check_parameters', 'PA0106', 1, 2, 0, 6, 0 )
3315             ENDIF
3316             unit = 'kg/kg'
3317
3318          CASE ( 'ql_c', 'ql_v', 'ql_vp' )
3319             IF (  .NOT.  cloud_droplets )  THEN
3320                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3321                                 'res cloud_droplets = .TRUE.'
3322                CALL message( 'check_parameters', 'PA0107', 1, 2, 0, 6, 0 )
3323             ENDIF
3324             IF ( TRIM( var ) == 'ql_c'  )  unit = 'kg/kg'
3325             IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
3326             IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
3327
3328          CASE ( 'qr' )
3329             IF (  .NOT.  cloud_physics )  THEN
3330                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3331                         'res cloud_physics = .TRUE.'
3332                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3333             ELSEIF ( .NOT.  microphysics_seifert ) THEN
3334                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3335                         'res cloud_scheme = seifert_beheng'
3336                CALL message( 'check_parameters', 'PA0359', 1, 2, 0, 6, 0 )
3337             ENDIF
3338             unit = 'kg/kg'
3339
3340          CASE ( 'qv' )
3341             IF (  .NOT.  cloud_physics )  THEN
3342                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3343                                 'res cloud_physics = .TRUE.'
3344                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3345             ENDIF
3346             unit = 'kg/kg'
3347
3348          CASE ( 'rho_ocean' )
3349             IF (  .NOT.  ocean )  THEN
3350                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3351                                 'res ocean = .TRUE.'
3352                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
3353             ENDIF
3354             unit = 'kg/m3'
3355
3356          CASE ( 's' )
3357             IF (  .NOT.  passive_scalar )  THEN
3358                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3359                                 'res passive_scalar = .TRUE.'
3360                CALL message( 'check_parameters', 'PA0110', 1, 2, 0, 6, 0 )
3361             ENDIF
3362             unit = 'kg/m3'
3363
3364          CASE ( 'sa' )
3365             IF (  .NOT.  ocean )  THEN
3366                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3367                                 'res ocean = .TRUE.'
3368                CALL message( 'check_parameters', 'PA0109', 1, 2, 0, 6, 0 )
3369             ENDIF
3370             unit = 'psu'
3371
3372          CASE ( 'p', 'pt', 'u', 'v', 'w' )
3373             IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
3374             IF ( TRIM( var ) == 'pt' )  unit = 'K'
3375             IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
3376             IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
3377             IF ( TRIM( var ) == 'w'  )  unit = 'm/s'
3378             CONTINUE
3379
3380          CASE ( 'ghf*', 'lwp*', 'ol*', 'pra*', 'prr*', 'qsws*', 'r_a*',       &
3381                 'shf*', 'ssws*', 't*', 'tsurf*', 'u*', 'z0*', 'z0h*', 'z0q*' )
3382             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
3383                message_string = 'illegal value for data_output: "' //         &
3384                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
3385                                 'cross sections are allowed for this value'
3386                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
3387             ENDIF
3388
3389             IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. cloud_physics )  THEN
3390                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3391                                 'res cloud_physics = .TRUE.'
3392                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3393             ENDIF
3394             IF ( TRIM( var ) == 'pra*'  .AND.                                 &
3395                  .NOT. ( microphysics_kessler .OR. microphysics_seifert ) )  THEN
3396                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3397                                 'res cloud_scheme = kessler or seifert_beheng'
3398                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3399             ENDIF
3400             IF ( TRIM( var ) == 'pra*'  .AND.  j == 1 )  THEN
3401                message_string = 'temporal averaging of precipitation ' //     &
3402                          'amount "' // TRIM( var ) // '" is not possible'
3403                CALL message( 'check_parameters', 'PA0113', 1, 2, 0, 6, 0 )
3404             ENDIF
3405             IF ( TRIM( var ) == 'prr*'  .AND.                                 &
3406                  .NOT. ( microphysics_kessler .OR. microphysics_seifert ) )  THEN
3407                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3408                                 'res cloud_scheme = kessler or seifert_beheng'
3409                CALL message( 'check_parameters', 'PA0112', 1, 2, 0, 6, 0 )
3410             ENDIF
3411             IF ( TRIM( var ) == 'qsws*'  .AND.  .NOT.  humidity )  THEN
3412                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3413                                 'res humidity = .TRUE.'
3414                CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
3415             ENDIF
3416
3417             IF ( TRIM( var ) == 'ghf*'  .AND.  .NOT.  land_surface )  THEN
3418                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3419                                 'res land_surface = .TRUE.'
3420                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3421             ENDIF
3422
3423             IF ( ( TRIM( var ) == 'r_a*' .OR.  TRIM( var ) == 'ghf*' )        &
3424                 .AND.  .NOT.  land_surface  .AND.  .NOT.  urban_surface )     &         
3425             THEN
3426                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3427                                 'res land_surface = .TRUE. or ' //            &
3428                                 'urban_surface = .TRUE.'
3429                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3430             ENDIF
3431             IF ( TRIM( var ) == 'ssws*'  .AND.  .NOT.  passive_scalar )  THEN
3432                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3433                                 'res passive_scalar = .TRUE.'
3434                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
3435             ENDIF
3436
3437             IF ( TRIM( var ) == 'ghf*'   )  unit = 'W/m2'
3438             IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/m2'
3439             IF ( TRIM( var ) == 'ol*'    )  unit = 'm'
3440             IF ( TRIM( var ) == 'pra*'   )  unit = 'mm'
3441             IF ( TRIM( var ) == 'prr*'   )  unit = 'mm/s'
3442             IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
3443             IF ( TRIM( var ) == 'r_a*'   )  unit = 's/m'     
3444             IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
3445             IF ( TRIM( var ) == 'ssws*'  )  unit = 'kg/m2*s'
3446             IF ( TRIM( var ) == 't*'     )  unit = 'K'
3447             IF ( TRIM( var ) == 'tsurf*' )  unit = 'K' 
3448             IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
3449             IF ( TRIM( var ) == 'z0*'    )  unit = 'm'
3450             IF ( TRIM( var ) == 'z0h*'   )  unit = 'm'
3451!
3452!--          Output of surface latent and sensible heat flux will be in W/m2
3453!--          in case of natural- and urban-type surfaces, even if
3454!--          flux_output_mode is set to kinematic units.
3455             IF ( land_surface  .OR.  urban_surface )  THEN
3456                IF ( TRIM( var ) == 'shf*'   )  unit = 'W/m2'
3457                IF ( TRIM( var ) == 'qsws*'  )  unit = 'W/m2'
3458             ENDIF
3459
3460          CASE DEFAULT
3461
3462             CALL tcm_check_data_output ( var, unit, i, ilen, k )
3463
3464             IF ( unit == 'illegal' )  THEN
3465                CALL lsm_check_data_output ( var, unit, i, ilen, k )
3466             ENDIF
3467
3468             IF ( unit == 'illegal' )  THEN
3469                CALL radiation_check_data_output( var, unit, i, ilen, k )
3470             ENDIF
3471
3472!
3473!--          Block of gust module outputs
3474             IF ( unit == 'illegal'  .AND.  gust_module_enabled  )  THEN
3475                CALL gust_check_data_output ( var, unit )
3476             ENDIF
3477
3478!
3479!--          Block of chemistry model outputs
3480             IF ( unit == 'illegal'  .AND.  air_chemistry                      &
3481                                     .AND.  var(1:3) == 'kc_' )  THEN
3482                CALL chem_check_data_output( var, unit, i, ilen, k )
3483             ENDIF
3484
3485!
3486!--          Block of urban surface model outputs
3487             IF ( unit == 'illegal' .AND. urban_surface .AND. var(1:4) == 'usm_' ) THEN
3488                 CALL usm_check_data_output( var, unit )
3489             ENDIF
3490
3491!
3492!--          Block of plant canopy model outputs
3493             IF ( unit == 'illegal' .AND. plant_canopy .AND. var(1:4) == 'pcm_' ) THEN
3494                 CALL pcm_check_data_output( var, unit )
3495             ENDIF
3496
3497!
3498!--          Block of uv exposure model outputs
3499             IF ( unit == 'illegal' .AND. uv_exposure .AND. var(1:5) == 'uvem_' )  THEN
3500                CALL uvem_check_data_output( var, unit, i, ilen, k )
3501             ENDIF
3502
3503             IF ( unit == 'illegal' )  THEN
3504                unit = ''
3505                CALL user_check_data_output( var, unit )
3506             ENDIF
3507
3508             IF ( unit == 'illegal' )  THEN
3509                IF ( data_output_user(1) /= ' ' )  THEN
3510                   message_string = 'illegal value for data_output or ' //     &
3511                         'data_output_user = "' // TRIM( data_output(i) ) // '"'
3512                   CALL message( 'check_parameters', 'PA0114', 1, 2, 0, 6, 0 )
3513                ELSE
3514                   message_string = 'illegal value for data_output = "' //     &
3515                                    TRIM( data_output(i) ) // '"'
3516                   CALL message( 'check_parameters', 'PA0115', 1, 2, 0, 6, 0 )
3517                ENDIF
3518             ENDIF
3519
3520       END SELECT
3521!
3522!--    Set the internal steering parameters appropriately
3523       IF ( k == 0 )  THEN
3524          do3d_no(j)              = do3d_no(j) + 1
3525          do3d(j,do3d_no(j))      = data_output(i)
3526          do3d_unit(j,do3d_no(j)) = unit
3527       ELSE
3528          do2d_no(j)              = do2d_no(j) + 1
3529          do2d(j,do2d_no(j))      = data_output(i)
3530          do2d_unit(j,do2d_no(j)) = unit
3531          IF ( data_output(i)(ilen-2:ilen) == '_xy' )  THEN
3532             data_output_xy(j) = .TRUE.
3533          ENDIF
3534          IF ( data_output(i)(ilen-2:ilen) == '_xz' )  THEN
3535             data_output_xz(j) = .TRUE.
3536          ENDIF
3537          IF ( data_output(i)(ilen-2:ilen) == '_yz' )  THEN
3538             data_output_yz(j) = .TRUE.
3539          ENDIF
3540       ENDIF
3541
3542       IF ( j == 1 )  THEN
3543!
3544!--       Check, if variable is already subject to averaging
3545          found = .FALSE.
3546          DO  k = 1, doav_n
3547             IF ( TRIM( doav(k) ) == TRIM( var ) )  found = .TRUE.
3548          ENDDO
3549
3550          IF ( .NOT. found )  THEN
3551             doav_n = doav_n + 1
3552             doav(doav_n) = var
3553          ENDIF
3554       ENDIF
3555
3556       i = i + 1
3557    ENDDO
3558
3559!
3560!-- Averaged 2d or 3d output requires that an averaging interval has been set
3561    IF ( doav_n > 0  .AND.  averaging_interval == 0.0_wp )  THEN
3562       WRITE( message_string, * )  'output of averaged quantity "',            &
3563                                   TRIM( doav(1) ), '_av" requires to set a ', &
3564                                   'non-zero averaging interval'
3565       CALL message( 'check_parameters', 'PA0323', 1, 2, 0, 6, 0 )
3566    ENDIF
3567
3568!
3569!-- Check sectional planes and store them in one shared array
3570    IF ( ANY( section_xy > nz + 1 ) )  THEN
3571       WRITE( message_string, * )  'section_xy must be <= nz + 1 = ', nz + 1
3572       CALL message( 'check_parameters', 'PA0319', 1, 2, 0, 6, 0 )
3573    ENDIF
3574    IF ( ANY( section_xz > ny + 1 ) )  THEN
3575       WRITE( message_string, * )  'section_xz must be <= ny + 1 = ', ny + 1
3576       CALL message( 'check_parameters', 'PA0320', 1, 2, 0, 6, 0 )
3577    ENDIF
3578    IF ( ANY( section_yz > nx + 1 ) )  THEN
3579       WRITE( message_string, * )  'section_yz must be <= nx + 1 = ', nx + 1
3580       CALL message( 'check_parameters', 'PA0321', 1, 2, 0, 6, 0 )
3581    ENDIF
3582    section(:,1) = section_xy
3583    section(:,2) = section_xz
3584    section(:,3) = section_yz
3585
3586!
3587!-- Upper plot limit for 3D arrays
3588    IF ( nz_do3d == -9999 )  nz_do3d = nzt + 1
3589
3590!
3591!-- Set output format string (used in header)
3592    SELECT CASE ( netcdf_data_format )
3593       CASE ( 1 )
3594          netcdf_data_format_string = 'netCDF classic'
3595       CASE ( 2 )
3596          netcdf_data_format_string = 'netCDF 64bit offset'
3597       CASE ( 3 )
3598          netcdf_data_format_string = 'netCDF4/HDF5'
3599       CASE ( 4 )
3600          netcdf_data_format_string = 'netCDF4/HDF5 classic'
3601       CASE ( 5 )
3602          netcdf_data_format_string = 'parallel netCDF4/HDF5'
3603       CASE ( 6 )
3604          netcdf_data_format_string = 'parallel netCDF4/HDF5 classic'
3605
3606    END SELECT
3607
3608!
3609!-- Check mask conditions
3610    DO mid = 1, max_masks
3611       IF ( data_output_masks(mid,1) /= ' '  .OR.                              &
3612            data_output_masks_user(mid,1) /= ' ' ) THEN
3613          masks = masks + 1
3614       ENDIF
3615    ENDDO
3616
3617    IF ( masks < 0  .OR.  masks > max_masks )  THEN
3618       WRITE( message_string, * )  'illegal value: masks must be >= 0 and ',   &
3619            '<= ', max_masks, ' (=max_masks)'
3620       CALL message( 'check_parameters', 'PA0325', 1, 2, 0, 6, 0 )
3621    ENDIF
3622    IF ( masks > 0 )  THEN
3623       mask_scale(1) = mask_scale_x
3624       mask_scale(2) = mask_scale_y
3625       mask_scale(3) = mask_scale_z
3626       IF ( ANY( mask_scale <= 0.0_wp ) )  THEN
3627          WRITE( message_string, * )                                           &
3628               'illegal value: mask_scale_x, mask_scale_y and mask_scale_z',   &
3629               'must be > 0.0'
3630          CALL message( 'check_parameters', 'PA0326', 1, 2, 0, 6, 0 )
3631       ENDIF
3632!
3633!--    Generate masks for masked data output
3634!--    Parallel netcdf output is not tested so far for masked data, hence
3635!--    netcdf_data_format is switched back to non-parallel output.
3636       netcdf_data_format_save = netcdf_data_format
3637       IF ( netcdf_data_format > 4 )  THEN
3638          IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
3639          IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
3640          message_string = 'netCDF file formats '//                            &
3641                           '5 (parallel netCDF 4) and ' //                     &
3642                           '6 (parallel netCDF 4 Classic model) '//            &
3643                           '& are currently not supported (not yet tested) ' //&
3644                           'for masked data. &Using respective non-parallel' //&
3645                           ' output for masked data.'
3646          CALL message( 'check_parameters', 'PA0383', 0, 0, 0, 6, 0 )
3647       ENDIF
3648       CALL init_masks
3649       netcdf_data_format = netcdf_data_format_save
3650    ENDIF
3651
3652!
3653!-- Check the NetCDF data format
3654    IF ( netcdf_data_format > 2 )  THEN
3655#if defined( __netcdf4 )
3656       CONTINUE
3657#else
3658       message_string = 'netCDF: netCDF4 format requested but no ' //          &
3659                        'cpp-directive __netcdf4 given & switch '  //          &
3660                        'back to 64-bit offset format'
3661       CALL message( 'check_parameters', 'PA0171', 0, 1, 0, 6, 0 )
3662       netcdf_data_format = 2
3663#endif
3664    ENDIF
3665    IF ( netcdf_data_format > 4 )  THEN
3666#if defined( __netcdf4 ) && defined( __netcdf4_parallel )
3667       CONTINUE
3668#else
3669       message_string = 'netCDF: netCDF4 parallel output requested but no ' // &
3670                        'cpp-directive __netcdf4_parallel given & switch '   //&
3671                        'back to netCDF4 non-parallel output'
3672       CALL message( 'check_parameters', 'PA0099', 0, 1, 0, 6, 0 )
3673       netcdf_data_format = netcdf_data_format - 2
3674#endif
3675    ENDIF
3676
3677!
3678!-- Calculate fixed number of output time levels for parallel netcdf output.
3679!-- The time dimension has to be defined as limited for parallel output,
3680!-- because otherwise the I/O performance drops significantly.
3681    IF ( netcdf_data_format > 4 )  THEN
3682
3683!
3684!--    Check if any of the follwoing data output interval is 0.0s, which is
3685!--    not allowed for parallel output.
3686       CALL check_dt_do( dt_do3d,           'dt_do3d'           )
3687       CALL check_dt_do( dt_do2d_xy,        'dt_do2d_xy'        )
3688       CALL check_dt_do( dt_do2d_xz,        'dt_do2d_xz'        )
3689       CALL check_dt_do( dt_do2d_yz,        'dt_do2d_yz'        )
3690       CALL check_dt_do( dt_data_output_av, 'dt_data_output_av' )
3691
3692!--    Set needed time levels (ntdim) to
3693!--    saved time levels + to be saved time levels.
3694       ntdim_3d(0) = do3d_time_count(0) + CEILING(                             &
3695                     ( end_time - MAX( skip_time_do3d,                         &
3696                                       simulated_time_at_begin )               &
3697                     ) / dt_do3d )
3698       IF ( do3d_at_begin ) ntdim_3d(0) = ntdim_3d(0) + 1
3699
3700       ntdim_3d(1) = do3d_time_count(1) + CEILING(                             &
3701                     ( end_time - MAX( skip_time_data_output_av,               &
3702                                       simulated_time_at_begin )               &
3703                     ) / dt_data_output_av )
3704
3705       ntdim_2d_xy(0) = do2d_xy_time_count(0) + CEILING(                       &
3706                        ( end_time - MAX( skip_time_do2d_xy,                   &
3707                                          simulated_time_at_begin )            &
3708                        ) / dt_do2d_xy )
3709
3710       ntdim_2d_xz(0) = do2d_xz_time_count(0) + CEILING(                       &
3711                        ( end_time - MAX( skip_time_do2d_xz,                   &
3712                                          simulated_time_at_begin )            &
3713                        ) / dt_do2d_xz )
3714
3715       ntdim_2d_yz(0) = do2d_yz_time_count(0) + CEILING(                       &
3716                        ( end_time - MAX( skip_time_do2d_yz,                   &
3717                                          simulated_time_at_begin )            &
3718                        ) / dt_do2d_yz )
3719
3720       IF ( do2d_at_begin )  THEN
3721          ntdim_2d_xy(0) = ntdim_2d_xy(0) + 1
3722          ntdim_2d_xz(0) = ntdim_2d_xz(0) + 1
3723          ntdim_2d_yz(0) = ntdim_2d_yz(0) + 1
3724       ENDIF
3725
3726       ntdim_2d_xy(1) = do2d_xy_time_count(1) + CEILING(                       &
3727                        ( end_time - MAX( skip_time_data_output_av,            &
3728                                          simulated_time_at_begin )            &
3729                        ) / dt_data_output_av )
3730
3731       ntdim_2d_xz(1) = do2d_xz_time_count(1) + CEILING(                       &
3732                        ( end_time - MAX( skip_time_data_output_av,            &
3733                                          simulated_time_at_begin )            &
3734                                 ) / dt_data_output_av )
3735
3736       ntdim_2d_yz(1) = do2d_yz_time_count(1) + CEILING(                       &
3737                        ( end_time - MAX( skip_time_data_output_av,            &
3738                                          simulated_time_at_begin )            &
3739                        ) / dt_data_output_av )
3740
3741    ENDIF
3742
3743!
3744!-- Check, whether a constant diffusion coefficient shall be used
3745    IF ( km_constant /= -1.0_wp )  THEN
3746       IF ( km_constant < 0.0_wp )  THEN
3747          WRITE( message_string, * )  'km_constant = ', km_constant, ' < 0.0'
3748          CALL message( 'check_parameters', 'PA0121', 1, 2, 0, 6, 0 )
3749       ELSE
3750          IF ( prandtl_number < 0.0_wp )  THEN
3751             WRITE( message_string, * )  'prandtl_number = ', prandtl_number,  &
3752                                         ' < 0.0'
3753             CALL message( 'check_parameters', 'PA0122', 1, 2, 0, 6, 0 )
3754          ENDIF
3755          constant_diffusion = .TRUE.
3756
3757          IF ( constant_flux_layer )  THEN
3758             message_string = 'constant_flux_layer is not allowed with fixed ' &
3759                              // 'value of km'
3760             CALL message( 'check_parameters', 'PA0123', 1, 2, 0, 6, 0 )
3761          ENDIF
3762       ENDIF
3763    ENDIF
3764
3765!
3766!-- In case of non-cyclic lateral boundaries and a damping layer for the
3767!-- potential temperature, check the width of the damping layer
3768    IF ( bc_lr /= 'cyclic' ) THEN
3769       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
3770            pt_damping_width > REAL( (nx+1) * dx ) )  THEN
3771          message_string = 'pt_damping_width out of range'
3772          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3773       ENDIF
3774    ENDIF
3775
3776    IF ( bc_ns /= 'cyclic' )  THEN
3777       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
3778            pt_damping_width > REAL( (ny+1) * dy ) )  THEN
3779          message_string = 'pt_damping_width out of range'
3780          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3781       ENDIF
3782    ENDIF
3783
3784!
3785!-- Check value range for zeta = z/L
3786    IF ( zeta_min >= zeta_max )  THEN
3787       WRITE( message_string, * )  'zeta_min = ', zeta_min, ' must be less ',  &
3788                                   'than zeta_max = ', zeta_max
3789       CALL message( 'check_parameters', 'PA0125', 1, 2, 0, 6, 0 )
3790    ENDIF
3791
3792!
3793!-- Check random generator
3794    IF ( (random_generator /= 'system-specific'      .AND.                     &
3795          random_generator /= 'random-parallel'   )  .AND.                     &
3796          random_generator /= 'numerical-recipes' )  THEN
3797       message_string = 'unknown random generator: random_generator = "' //    &
3798                        TRIM( random_generator ) // '"'
3799       CALL message( 'check_parameters', 'PA0135', 1, 2, 0, 6, 0 )
3800    ENDIF
3801
3802!
3803!-- Determine upper and lower hight level indices for random perturbations
3804    IF ( disturbance_level_b == -9999999.9_wp )  THEN
3805       IF ( ocean )  THEN
3806          disturbance_level_b     = zu((nzt*2)/3)
3807          disturbance_level_ind_b = ( nzt * 2 ) / 3
3808       ELSE
3809          disturbance_level_b     = zu(nzb+3)
3810          disturbance_level_ind_b = nzb + 3
3811       ENDIF
3812    ELSEIF ( disturbance_level_b < zu(3) )  THEN
3813       WRITE( message_string, * )  'disturbance_level_b = ',                   &
3814                           disturbance_level_b, ' must be >= ', zu(3), '(zu(3))'
3815       CALL message( 'check_parameters', 'PA0126', 1, 2, 0, 6, 0 )
3816    ELSEIF ( disturbance_level_b > zu(nzt-2) )  THEN
3817       WRITE( message_string, * )  'disturbance_level_b = ',                   &
3818                   disturbance_level_b, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3819       CALL message( 'check_parameters', 'PA0127', 1, 2, 0, 6, 0 )
3820    ELSE
3821       DO  k = 3, nzt-2
3822          IF ( disturbance_level_b <= zu(k) )  THEN
3823             disturbance_level_ind_b = k
3824             EXIT
3825          ENDIF
3826       ENDDO
3827    ENDIF
3828
3829    IF ( disturbance_level_t == -9999999.9_wp )  THEN
3830       IF ( ocean )  THEN
3831          disturbance_level_t     = zu(nzt-3)
3832          disturbance_level_ind_t = nzt - 3
3833       ELSE
3834          disturbance_level_t     = zu(nzt/3)
3835          disturbance_level_ind_t = nzt / 3
3836       ENDIF
3837    ELSEIF ( disturbance_level_t > zu(nzt-2) )  THEN
3838       WRITE( message_string, * )  'disturbance_level_t = ',                   &
3839                   disturbance_level_t, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3840       CALL message( 'check_parameters', 'PA0128', 1, 2, 0, 6, 0 )
3841    ELSEIF ( disturbance_level_t < disturbance_level_b )  THEN
3842       WRITE( message_string, * )  'disturbance_level_t = ',                   &
3843                   disturbance_level_t, ' must be >= disturbance_level_b = ',  &
3844                   disturbance_level_b
3845       CALL message( 'check_parameters', 'PA0129', 1, 2, 0, 6, 0 )
3846    ELSE
3847       DO  k = 3, nzt-2
3848          IF ( disturbance_level_t <= zu(k) )  THEN
3849             disturbance_level_ind_t = k
3850             EXIT
3851          ENDIF
3852       ENDDO
3853    ENDIF
3854
3855!
3856!-- Check again whether the levels determined this way are ok.
3857!-- Error may occur at automatic determination and too few grid points in
3858!-- z-direction.
3859    IF ( disturbance_level_ind_t < disturbance_level_ind_b )  THEN
3860       WRITE( message_string, * )  'disturbance_level_ind_t = ',               &
3861                disturbance_level_ind_t, ' must be >= ',                       &
3862                'disturbance_level_ind_b = ', disturbance_level_ind_b
3863       CALL message( 'check_parameters', 'PA0130', 1, 2, 0, 6, 0 )
3864    ENDIF
3865
3866!
3867!-- Determine the horizontal index range for random perturbations.
3868!-- In case of non-cyclic horizontal boundaries, no perturbations are imposed
3869!-- near the inflow and the perturbation area is further limited to ...(1)
3870!-- after the initial phase of the flow.
3871
3872    IF ( bc_lr /= 'cyclic' )  THEN
3873       IF ( inflow_disturbance_begin == -1 )  THEN
3874          inflow_disturbance_begin = MIN( 10, nx/2 )
3875       ENDIF
3876       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > nx )&
3877       THEN
3878          message_string = 'inflow_disturbance_begin out of range'
3879          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
3880       ENDIF
3881       IF ( inflow_disturbance_end == -1 )  THEN
3882          inflow_disturbance_end = MIN( 100, 3*nx/4 )
3883       ENDIF
3884       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > nx )    &
3885       THEN
3886          message_string = 'inflow_disturbance_end out of range'
3887          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
3888       ENDIF
3889    ELSEIF ( bc_ns /= 'cyclic' )  THEN
3890       IF ( inflow_disturbance_begin == -1 )  THEN
3891          inflow_disturbance_begin = MIN( 10, ny/2 )
3892       ENDIF
3893       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > ny )&
3894       THEN
3895          message_string = 'inflow_disturbance_begin out of range'
3896          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
3897       ENDIF
3898       IF ( inflow_disturbance_end == -1 )  THEN
3899          inflow_disturbance_end = MIN( 100, 3*ny/4 )
3900       ENDIF
3901       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > ny )    &
3902       THEN
3903          message_string = 'inflow_disturbance_end out of range'
3904          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
3905       ENDIF
3906    ENDIF
3907
3908    IF ( random_generator == 'random-parallel' )  THEN
3909       dist_nxl = nxl;  dist_nxr = nxr
3910       dist_nys = nys;  dist_nyn = nyn
3911       IF ( bc_lr == 'radiation/dirichlet' )  THEN
3912          dist_nxr    = MIN( nx - inflow_disturbance_begin, nxr )
3913          dist_nxl(1) = MAX( nx - inflow_disturbance_end, nxl )
3914       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
3915          dist_nxl    = MAX( inflow_disturbance_begin, nxl )
3916          dist_nxr(1) = MIN( inflow_disturbance_end, nxr )
3917       ELSEIF ( bc_lr == 'nested'  .OR.  bc_lr == 'forcing' )  THEN
3918          dist_nxl    = MAX( inflow_disturbance_begin, nxl )
3919          dist_nxr    = MIN( nx - inflow_disturbance_begin, nxr )
3920       ENDIF
3921       IF ( bc_ns == 'dirichlet/radiation' )  THEN
3922          dist_nyn    = MIN( ny - inflow_disturbance_begin, nyn )
3923          dist_nys(1) = MAX( ny - inflow_disturbance_end, nys )
3924       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
3925          dist_nys    = MAX( inflow_disturbance_begin, nys )
3926          dist_nyn(1) = MIN( inflow_disturbance_end, nyn )
3927       ELSEIF ( bc_ns == 'nested'  .OR.  bc_ns == 'forcing' )  THEN
3928          dist_nys    = MAX( inflow_disturbance_begin, nys )
3929          dist_nyn    = MIN( ny - inflow_disturbance_begin, nyn )
3930       ENDIF
3931    ELSE
3932       dist_nxl = 0;  dist_nxr = nx
3933       dist_nys = 0;  dist_nyn = ny
3934       IF ( bc_lr == 'radiation/dirichlet' )  THEN
3935          dist_nxr    = nx - inflow_disturbance_begin
3936          dist_nxl(1) = nx - inflow_disturbance_end
3937       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
3938          dist_nxl    = inflow_disturbance_begin
3939          dist_nxr(1) = inflow_disturbance_end
3940       ELSEIF ( bc_lr == 'nested'  .OR.  bc_lr == 'forcing' )  THEN
3941          dist_nxr    = nx - inflow_disturbance_begin
3942          dist_nxl    = inflow_disturbance_begin
3943       ENDIF
3944       IF ( bc_ns == 'dirichlet/radiation' )  THEN
3945          dist_nyn    = ny - inflow_disturbance_begin
3946          dist_nys(1) = ny - inflow_disturbance_end
3947       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
3948          dist_nys    = inflow_disturbance_begin
3949          dist_nyn(1) = inflow_disturbance_end
3950       ELSEIF ( bc_ns == 'nested'  .OR.  bc_ns == 'forcing' )  THEN
3951          dist_nyn    = ny - inflow_disturbance_begin
3952          dist_nys    = inflow_disturbance_begin
3953       ENDIF
3954    ENDIF
3955
3956!
3957!-- A turbulent inflow requires Dirichlet conditions at the respective inflow
3958!-- boundary (so far, a turbulent inflow is realized from the left side only)
3959    IF ( turbulent_inflow  .AND.  bc_lr /= 'dirichlet/radiation' )  THEN
3960       message_string = 'turbulent_inflow = .T. requires a Dirichlet ' //      &
3961                        'condition at the inflow boundary'
3962       CALL message( 'check_parameters', 'PA0133', 1, 2, 0, 6, 0 )
3963    ENDIF
3964
3965!
3966!-- Turbulent inflow requires that 3d arrays have been cyclically filled with
3967!-- data from prerun in the first main run
3968    IF ( turbulent_inflow  .AND.  initializing_actions /= 'cyclic_fill'  .AND. &
3969         initializing_actions /= 'read_restart_data' )  THEN
3970       message_string = 'turbulent_inflow = .T. requires ' //                  &
3971                        'initializing_actions = ''cyclic_fill'' or ' //        &
3972                        'initializing_actions = ''read_restart_data'' '
3973       CALL message( 'check_parameters', 'PA0055', 1, 2, 0, 6, 0 )
3974    ENDIF
3975
3976!
3977!-- In case of turbulent inflow calculate the index of the recycling plane
3978    IF ( turbulent_inflow )  THEN
3979       IF ( recycling_width <= dx  .OR.  recycling_width >= nx * dx )  THEN
3980          WRITE( message_string, * )  'illegal value for recycling_width: ',   &
3981                                      recycling_width
3982          CALL message( 'check_parameters', 'PA0134', 1, 2, 0, 6, 0 )
3983       ENDIF
3984!
3985!--    Calculate the index
3986       recycling_plane = recycling_width / dx
3987!
3988!--    Because the y-shift is done with a distance of INT( npey / 2 ) no shift
3989!--    is possible if there is only one PE in y direction.
3990       IF ( recycling_yshift .AND. pdims(2) < 2 )  THEN
3991          WRITE( message_string, * )  'recycling_yshift = .T. requires more',  &
3992                                      ' than one processor in y direction'
3993          CALL message( 'check_parameters', 'PA0421', 1, 2, 0, 6, 0 )
3994       ENDIF
3995    ENDIF
3996
3997
3998    IF ( turbulent_outflow )  THEN
3999!
4000!--    Turbulent outflow requires Dirichlet conditions at the respective inflow
4001!--    boundary (so far, a turbulent outflow is realized at the right side only)
4002       IF ( bc_lr /= 'dirichlet/radiation' )  THEN
4003          message_string = 'turbulent_outflow = .T. requires ' //              &
4004                           'bc_lr = "dirichlet/radiation"'
4005          CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 )
4006       ENDIF
4007!
4008!--    The ouflow-source plane must lay inside the model domain
4009       IF ( outflow_source_plane < dx  .OR.  &
4010            outflow_source_plane > nx * dx )  THEN
4011          WRITE( message_string, * )  'illegal value for outflow_source'//     &
4012                                      '_plane: ', outflow_source_plane
4013          CALL message( 'check_parameters', 'PA0145', 1, 2, 0, 6, 0 )
4014       ENDIF
4015    ENDIF
4016
4017!
4018!-- Determine damping level index for 1D model
4019    IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
4020       IF ( damp_level_1d == -1.0_wp )  THEN
4021          damp_level_1d     = zu(nzt+1)
4022          damp_level_ind_1d = nzt + 1
4023       ELSEIF ( damp_level_1d < 0.0_wp  .OR.  damp_level_1d > zu(nzt+1) )  THEN
4024          WRITE( message_string, * )  'damp_level_1d = ', damp_level_1d,       &
4025                 ' must be >= 0.0 and <= ', zu(nzt+1), '(zu(nzt+1))'
4026          CALL message( 'check_parameters', 'PA0136', 1, 2, 0, 6, 0 )
4027       ELSE
4028          DO  k = 1, nzt+1
4029             IF ( damp_level_1d <= zu(k) )  THEN
4030                damp_level_ind_1d = k
4031                EXIT
4032             ENDIF
4033          ENDDO
4034       ENDIF
4035    ENDIF
4036
4037!
4038!-- Check some other 1d-model parameters
4039    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                   &
4040         TRIM( mixing_length_1d ) /= 'blackadar' )  THEN
4041       message_string = 'mixing_length_1d = "' // TRIM( mixing_length_1d ) //  &
4042                        '" is unknown'
4043       CALL message( 'check_parameters', 'PA0137', 1, 2, 0, 6, 0 )
4044    ENDIF
4045    IF ( TRIM( dissipation_1d ) /= 'as_in_3d_model'  .AND.                     &
4046         TRIM( dissipation_1d ) /= 'detering'  .AND.                           &
4047         TRIM( dissipation_1d ) /= 'prognostic' )  THEN
4048       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) //      &
4049                        '" is unknown'
4050       CALL message( 'check_parameters', 'PA0138', 1, 2, 0, 6, 0 )
4051    ENDIF
4052    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                   &
4053         TRIM( dissipation_1d ) == 'as_in_3d_model' )  THEN
4054       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) //      &
4055                        '" requires mixing_length_1d = "as_in_3d_model"'
4056       CALL message( 'check_parameters', 'PA0485', 1, 2, 0, 6, 0 )
4057    ENDIF
4058
4059!
4060!-- Set time for the next user defined restart (time_restart is the
4061!-- internal parameter for steering restart events)
4062    IF ( restart_time /= 9999999.9_wp )  THEN
4063       IF ( restart_time > time_since_reference_point )  THEN
4064          time_restart = restart_time
4065       ENDIF
4066    ELSE
4067!
4068!--    In case of a restart run, set internal parameter to default (no restart)
4069!--    if the NAMELIST-parameter restart_time is omitted
4070       time_restart = 9999999.9_wp
4071    ENDIF
4072
4073!
4074!-- Check pressure gradient conditions
4075    IF ( dp_external  .AND.  conserve_volume_flow )  THEN
4076       WRITE( message_string, * )  'Both dp_external and conserve_volume_flo', &
4077            'w are .TRUE. but one of them must be .FALSE.'
4078       CALL message( 'check_parameters', 'PA0150', 1, 2, 0, 6, 0 )
4079    ENDIF
4080    IF ( dp_external )  THEN
4081       IF ( dp_level_b < zu(nzb)  .OR.  dp_level_b > zu(nzt) )  THEN
4082          WRITE( message_string, * )  'dp_level_b = ', dp_level_b, ' is out ', &
4083               ' of range [zu(nzb), zu(nzt)]'
4084          CALL message( 'check_parameters', 'PA0151', 1, 2, 0, 6, 0 )
4085       ENDIF
4086       IF ( .NOT. ANY( dpdxy /= 0.0_wp ) )  THEN
4087          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is ze', &
4088               'ro, i.e. the external pressure gradient will not be applied'
4089          CALL message( 'check_parameters', 'PA0152', 0, 1, 0, 6, 0 )
4090       ENDIF
4091    ENDIF
4092    IF ( ANY( dpdxy /= 0.0_wp )  .AND.  .NOT.  dp_external )  THEN
4093       WRITE( message_string, * )  'dpdxy is nonzero but dp_external is ',     &
4094            '.FALSE., i.e. the external pressure gradient & will not be applied'
4095       CALL message( 'check_parameters', 'PA0153', 0, 1, 0, 6, 0 )
4096    ENDIF
4097    IF ( conserve_volume_flow )  THEN
4098       IF ( TRIM( conserve_volume_flow_mode ) == 'default' )  THEN
4099
4100          conserve_volume_flow_mode = 'initial_profiles'
4101
4102       ELSEIF ( TRIM( conserve_volume_flow_mode ) /= 'initial_profiles' .AND.  &
4103            TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' )  THEN
4104          WRITE( message_string, * )  'unknown conserve_volume_flow_mode: ',   &
4105               conserve_volume_flow_mode
4106          CALL message( 'check_parameters', 'PA0154', 1, 2, 0, 6, 0 )
4107       ENDIF
4108       IF ( (bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic')  .AND.                &
4109          TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
4110          WRITE( message_string, * )  'non-cyclic boundary conditions ',       &
4111               'require  conserve_volume_flow_mode = ''initial_profiles'''
4112          CALL message( 'check_parameters', 'PA0155', 1, 2, 0, 6, 0 )
4113       ENDIF
4114    ENDIF
4115    IF ( ( u_bulk /= 0.0_wp  .OR.  v_bulk /= 0.0_wp )  .AND.                   &
4116         ( .NOT. conserve_volume_flow  .OR.                                    &
4117         TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) )  THEN
4118       WRITE( message_string, * )  'nonzero bulk velocity requires ',          &
4119            'conserve_volume_flow = .T. and ',                                 &
4120            'conserve_volume_flow_mode = ''bulk_velocity'''
4121       CALL message( 'check_parameters', 'PA0157', 1, 2, 0, 6, 0 )
4122    ENDIF
4123
4124!
4125!-- Check particle attributes
4126    IF ( particle_color /= 'none' )  THEN
4127       IF ( particle_color /= 'absuv'  .AND.  particle_color /= 'pt*'  .AND.   &
4128            particle_color /= 'z' )  THEN
4129          message_string = 'illegal value for parameter particle_color: ' //   &
4130                           TRIM( particle_color)
4131          CALL message( 'check_parameters', 'PA0313', 1, 2, 0, 6, 0 )
4132       ELSE
4133          IF ( color_interval(2) <= color_interval(1) )  THEN
4134             message_string = 'color_interval(2) <= color_interval(1)'
4135             CALL message( 'check_parameters', 'PA0315', 1, 2, 0, 6, 0 )
4136          ENDIF
4137       ENDIF
4138    ENDIF
4139
4140    IF ( particle_dvrpsize /= 'none' )  THEN
4141       IF ( particle_dvrpsize /= 'absw' )  THEN
4142          message_string = 'illegal value for parameter particle_dvrpsize:' // &
4143                           ' ' // TRIM( particle_dvrpsize)
4144          CALL message( 'check_parameters', 'PA0314', 1, 2, 0, 6, 0 )
4145       ELSE
4146          IF ( dvrpsize_interval(2) <= dvrpsize_interval(1) )  THEN
4147             message_string = 'dvrpsize_interval(2) <= dvrpsize_interval(1)'
4148             CALL message( 'check_parameters', 'PA0316', 1, 2, 0, 6, 0 )
4149          ENDIF
4150       ENDIF
4151    ENDIF
4152
4153!
4154!-- Prevent empty time records in volume, cross-section and masked data in case
4155!-- of non-parallel netcdf-output in restart runs
4156    IF ( netcdf_data_format < 5 )  THEN
4157       IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
4158          do3d_time_count    = 0
4159          do2d_xy_time_count = 0
4160          do2d_xz_time_count = 0
4161          do2d_yz_time_count = 0
4162          domask_time_count  = 0
4163       ENDIF
4164    ENDIF
4165
4166!
4167!-- Check for valid setting of most_method
4168    IF ( TRIM( most_method ) /= 'circular'  .AND.                              &
4169         TRIM( most_method ) /= 'newton'    .AND.                              &
4170         TRIM( most_method ) /= 'lookup' )  THEN
4171       message_string = 'most_method = "' // TRIM( most_method ) //            &
4172                        '" is unknown'
4173       CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 )
4174    ENDIF
4175
4176!
4177!-- Check roughness length, which has to be smaller than dz/2
4178    IF ( ( constant_flux_layer .OR.  &
4179           INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )       &
4180         .AND. roughness_length >= 0.5 * dz )  THEN
4181       message_string = 'roughness_length must be smaller than dz/2'
4182       CALL message( 'check_parameters', 'PA0424', 1, 2, 0, 6, 0 )
4183    ENDIF
4184
4185!
4186!-- Vertical nesting: check fine and coarse grid compatibility for data exchange
4187    IF ( vnested )  CALL vnest_check_parameters
4188
4189!
4190!-- Check if topography is read from file in case of complex terrain simulations
4191    IF ( complex_terrain  .AND.  TRIM( topography ) /= 'read_from_file' )  THEN
4192       message_string = 'complex_terrain requires topography' //               &
4193                        ' = ''read_from_file'''
4194       CALL message( 'check_parameters', 'PA0295', 1, 2, 0, 6, 0 )
4195    ENDIF
4196
4197!
4198!-- Check if vertical grid stretching is switched off in case of complex
4199!-- terrain simulations
4200    IF ( complex_terrain  .AND.  dz_stretch_level < 100000.0_wp )  THEN
4201       message_string = 'Vertical grid stretching is not allowed for ' //      &
4202                        'complex_terrain = .T.'
4203       CALL message( 'check_parameters', 'PA0473', 1, 2, 0, 6, 0 )
4204    ENDIF
4205
4206    CALL location_message( 'finished', .TRUE. )
4207!
4208!-- Check &userpar parameters
4209    CALL user_check_parameters
4210
4211 CONTAINS
4212
4213!------------------------------------------------------------------------------!
4214! Description:
4215! ------------
4216!> Check the length of data output intervals. In case of parallel NetCDF output
4217!> the time levels of the output files need to be fixed. Therefore setting the
4218!> output interval to 0.0s (usually used to output each timestep) is not
4219!> possible as long as a non-fixed timestep is used.
4220!------------------------------------------------------------------------------!
4221
4222    SUBROUTINE check_dt_do( dt_do, dt_do_name )
4223
4224       IMPLICIT NONE
4225
4226       CHARACTER (LEN=*), INTENT (IN) :: dt_do_name !< parin variable name
4227
4228       REAL(wp), INTENT (INOUT)       :: dt_do      !< data output interval
4229
4230       IF ( dt_do == 0.0_wp )  THEN
4231          IF ( dt_fixed )  THEN
4232             WRITE( message_string, '(A,F9.4,A)' )  'Output at every '  //     &
4233                    'timestep is wanted (' // dt_do_name // ' = 0.0).&'//      &
4234                    'The output interval is set to the fixed timestep dt '//   &
4235                    '= ', dt, 's.'
4236             CALL message( 'check_parameters', 'PA0060', 0, 0, 0, 6, 0 )
4237             dt_do = dt
4238          ELSE
4239             message_string = dt_do_name // ' = 0.0 while using a ' //         &
4240                              'variable timestep and parallel netCDF4 ' //     &
4241                              'is not allowed.'
4242             CALL message( 'check_parameters', 'PA0081', 1, 2, 0, 6, 0 )
4243          ENDIF
4244       ENDIF
4245
4246    END SUBROUTINE check_dt_do
4247
4248
4249
4250
4251!------------------------------------------------------------------------------!
4252! Description:
4253! ------------
4254!> Inititalizes the vertical profiles of scalar quantities.
4255!------------------------------------------------------------------------------!
4256
4257    SUBROUTINE init_vertical_profiles( vertical_gradient_level_ind,            &
4258                                       vertical_gradient_level,                &
4259                                       vertical_gradient,                      &
4260                                       pr_init, surface_value, bc_t_val )
4261
4262
4263       IMPLICIT NONE
4264
4265       INTEGER(iwp) ::  i     !< counter
4266       INTEGER(iwp), DIMENSION(1:10) ::  vertical_gradient_level_ind !< vertical grid indices for gradient levels
4267
4268       REAL(wp)     ::  bc_t_val      !< model top gradient
4269       REAL(wp)     ::  gradient      !< vertica gradient of the respective quantity
4270       REAL(wp)     ::  surface_value !< surface value of the respecitve quantity
4271
4272
4273       REAL(wp), DIMENSION(0:nz+1) ::  pr_init                 !< initialisation profile
4274       REAL(wp), DIMENSION(1:10)   ::  vertical_gradient       !< given vertical gradient
4275       REAL(wp), DIMENSION(1:10)   ::  vertical_gradient_level !< given vertical gradient level
4276
4277       i = 1
4278       gradient = 0.0_wp
4279
4280       IF (  .NOT.  ocean )  THEN
4281
4282          vertical_gradient_level_ind(1) = 0
4283          DO  k = 1, nzt+1
4284             IF ( i < 11 )  THEN
4285                IF ( vertical_gradient_level(i) < zu(k)  .AND.            &
4286                     vertical_gradient_level(i) >= 0.0_wp )  THEN
4287                   gradient = vertical_gradient(i) / 100.0_wp
4288                   vertical_gradient_level_ind(i) = k - 1
4289                   i = i + 1
4290                ENDIF
4291             ENDIF
4292             IF ( gradient /= 0.0_wp )  THEN
4293                IF ( k /= 1 )  THEN
4294                   pr_init(k) = pr_init(k-1) + dzu(k) * gradient
4295                ELSE
4296                   pr_init(k) = pr_init(k-1) + dzu(k) * gradient
4297                ENDIF
4298             ELSE
4299                pr_init(k) = pr_init(k-1)
4300             ENDIF
4301   !
4302   !--       Avoid negative values
4303             IF ( pr_init(k) < 0.0_wp )  THEN
4304                pr_init(k) = 0.0_wp
4305             ENDIF
4306          ENDDO
4307
4308       ELSE
4309
4310          vertical_gradient_level_ind(1) = nzt+1
4311          DO  k = nzt, 0, -1
4312             IF ( i < 11 )  THEN
4313                IF ( vertical_gradient_level(i) > zu(k)  .AND.            &
4314                     vertical_gradient_level(i) <= 0.0_wp )  THEN
4315                   gradient = vertical_gradient(i) / 100.0_wp
4316                   vertical_gradient_level_ind(i) = k + 1
4317                   i = i + 1
4318                ENDIF
4319             ENDIF
4320             IF ( gradient /= 0.0_wp )  THEN
4321                IF ( k /= nzt )  THEN
4322                   pr_init(k) = pr_init(k+1) - dzu(k+1) * gradient
4323                ELSE
4324                   pr_init(k)   = surface_value - 0.5_wp * dzu(k+1) * gradient
4325                   pr_init(k+1) = surface_value + 0.5_wp * dzu(k+1) * gradient
4326                ENDIF
4327             ELSE
4328                pr_init(k) = pr_init(k+1)
4329             ENDIF
4330   !
4331   !--       Avoid negative humidities
4332             IF ( pr_init(k) < 0.0_wp )  THEN
4333                pr_init(k) = 0.0_wp
4334             ENDIF
4335          ENDDO
4336
4337       ENDIF
4338
4339!
4340!--    In case of no given gradients, choose zero gradient conditions
4341       IF ( vertical_gradient_level(1) == -999999.9_wp )  THEN
4342          vertical_gradient_level(1) = 0.0_wp
4343       ENDIF
4344!
4345!--    Store gradient at the top boundary for possible Neumann boundary condition
4346       bc_t_val  = ( pr_init(nzt+1) - pr_init(nzt) ) / dzu(nzt+1)
4347
4348    END SUBROUTINE init_vertical_profiles
4349
4350
4351
4352!------------------------------------------------------------------------------!
4353! Description:
4354! ------------
4355!> Set the bottom and top boundary conditions for humidity and scalars.
4356!------------------------------------------------------------------------------!
4357
4358    SUBROUTINE set_bc_scalars( sq, bc_b, bc_t, ibc_b, ibc_t, err_nr_b, err_nr_t )
4359
4360
4361       IMPLICIT NONE
4362
4363       CHARACTER (LEN=1)   ::  sq         !< name of scalar quantity
4364       CHARACTER (LEN=*)   ::  bc_b       !< bottom boundary condition
4365       CHARACTER (LEN=*)   ::  bc_t       !< top boundary condition
4366       CHARACTER (LEN=*)   ::  err_nr_b   !< error number if bottom bc is unknown
4367       CHARACTER (LEN=*)   ::  err_nr_t   !< error number if top bc is unknown
4368
4369       INTEGER(iwp)        ::  ibc_b      !< index for bottom boundary condition
4370       INTEGER(iwp)        ::  ibc_t      !< index for top boundary condition
4371
4372!
4373!--    Set Integer flags and check for possilbe errorneous settings for bottom
4374!--    boundary condition
4375       IF ( bc_b == 'dirichlet' )  THEN
4376          ibc_b = 0
4377       ELSEIF ( bc_b == 'neumann' )  THEN
4378          ibc_b = 1
4379       ELSE
4380          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
4381                           '_b ="' // TRIM( bc_b ) // '"'
4382          CALL message( 'check_parameters', err_nr_b, 1, 2, 0, 6, 0 )
4383       ENDIF
4384!
4385!--    Set Integer flags and check for possilbe errorneous settings for top
4386!--    boundary condition
4387       IF ( bc_t == 'dirichlet' )  THEN
4388          ibc_t = 0
4389       ELSEIF ( bc_t == 'neumann' )  THEN
4390          ibc_t = 1
4391       ELSEIF ( bc_t == 'initial_gradient' )  THEN
4392          ibc_t = 2
4393       ELSEIF ( bc_t == 'nested'  .OR.  bc_t == 'forcing' )  THEN
4394          ibc_t = 3
4395       ELSE
4396          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
4397                           '_t ="' // TRIM( bc_t ) // '"'
4398          CALL message( 'check_parameters', err_nr_t, 1, 2, 0, 6, 0 )
4399       ENDIF
4400
4401
4402    END SUBROUTINE set_bc_scalars
4403
4404
4405
4406!------------------------------------------------------------------------------!
4407! Description:
4408! ------------
4409!> Check for consistent settings of bottom boundary conditions for humidity
4410!> and scalars.
4411!------------------------------------------------------------------------------!
4412
4413    SUBROUTINE check_bc_scalars( sq, bc_b, ibc_b,                      &
4414                                 err_nr_1, err_nr_2,                   &
4415                                 constant_flux, surface_initial_change )
4416
4417
4418       IMPLICIT NONE
4419
4420       CHARACTER (LEN=1)   ::  sq                       !< name of scalar quantity
4421       CHARACTER (LEN=*)   ::  bc_b                     !< bottom boundary condition
4422       CHARACTER (LEN=*)   ::  err_nr_1                 !< error number of first error
4423       CHARACTER (LEN=*)   ::  err_nr_2                 !< error number of second error
4424
4425       INTEGER(iwp)        ::  ibc_b                    !< index of bottom boundary condition
4426
4427       LOGICAL             ::  constant_flux            !< flag for constant-flux layer
4428
4429       REAL(wp)            ::  surface_initial_change   !< value of initial change at the surface
4430
4431!
4432!--    A given surface value implies Dirichlet boundary condition for
4433!--    the respective quantity. In this case specification of a constant flux is
4434!--    forbidden. However, an exception is made for large-scale forcing as well
4435!--    as land-surface model.
4436       IF ( .NOT. land_surface  .AND.  .NOT. large_scale_forcing )  THEN
4437          IF ( ibc_b == 0  .AND.  constant_flux )  THEN
4438             message_string = 'boundary condition: bc_' // TRIM( sq ) //       &
4439                              '_b ' // '= "' // TRIM( bc_b ) //                &
4440                              '" is not allowed with prescribed surface flux'
4441             CALL message( 'check_parameters', err_nr_1, 1, 2, 0, 6, 0 )
4442          ENDIF
4443       ENDIF
4444       IF ( constant_flux  .AND.  surface_initial_change /= 0.0_wp )  THEN
4445          WRITE( message_string, * )  'a prescribed surface flux is not allo', &
4446                 'wed with ', sq, '_surface_initial_change (/=0) = ',          &
4447                 surface_initial_change
4448          CALL message( 'check_parameters', err_nr_2, 1, 2, 0, 6, 0 )
4449       ENDIF
4450
4451
4452    END SUBROUTINE check_bc_scalars
4453
4454
4455
4456 END SUBROUTINE check_parameters
Note: See TracBrowser for help on using the repository browser.