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

Last change on this file since 3347 was 3347, checked in by suehring, 3 years ago

Offline nesting revised and separated from large_scale_forcing_mod; Time-dependent synthetic turbulence generator; bugfixes in USM/LSM radiation model and init_pegrid

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