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

Last change on this file since 3322 was 3312, checked in by knoop, 6 years ago

Added proper exit code to PALM and fixed exit code handling by palmbuild and palmrun

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