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

Last change on this file since 3700 was 3700, checked in by knoop, 4 years ago

Moved user_define_netdf_grid into user_module.f90
Added module interface for the definition of additional timeseries

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