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

Last change on this file since 3337 was 3337, checked in by kanani, 3 years ago

reintegrate branch resler to trunk

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