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

Last change on this file since 3646 was 3646, checked in by kanani, 5 years ago

Bugfix: replace simulated_time by time_since_reference_point where required

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