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

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

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