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

Last change on this file since 4125 was 4069, checked in by Giersch, 5 years ago

Bugfix for masked output, compiler warning removed, test case for wind turbine model revised

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