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

Last change on this file since 4142 was 4142, checked in by suehring, 2 years ago

Consider spinup in number of output timesteps for averaged 2D output

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