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

Last change on this file since 3401 was 3377, checked in by knoop, 6 years ago

Removed check that prohibits use of bulk_cloud_model with topography

  • 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.7 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 3377 2018-10-19 12:26:54Z eckhard $
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 ( cloud_droplets )  THEN
1130          WRITE( action, '(A)' )  'cloud_droplets = .TRUE.'
1131       ENDIF
1132       IF ( .NOT. constant_flux_layer )  THEN
1133          WRITE( action, '(A)' )  'constant_flux_layer = .FALSE.'
1134       ENDIF
1135       IF ( action /= ' ' )  THEN
1136          message_string = 'a non-flat topography does not allow ' //          &
1137                           TRIM( action )
1138          CALL message( 'check_parameters', 'PA0014', 1, 2, 0, 6, 0 )
1139       ENDIF
1140
1141    ENDIF
1142
1143!
1144!-- Check approximation
1145    IF ( TRIM( approximation ) /= 'boussinesq'   .AND.                         &
1146         TRIM( approximation ) /= 'anelastic' )  THEN
1147       message_string = 'unknown approximation: approximation = "' //          &
1148                        TRIM( approximation ) // '"'
1149       CALL message( 'check_parameters', 'PA0446', 1, 2, 0, 6, 0 )
1150    ENDIF
1151
1152!
1153!-- Check approximation requirements
1154    IF ( TRIM( approximation ) == 'anelastic'   .AND.                          &
1155         TRIM( momentum_advec ) /= 'ws-scheme' )  THEN
1156       message_string = 'Anelastic approximation requires ' //                 &
1157                        'momentum_advec = "ws-scheme"'
1158       CALL message( 'check_parameters', 'PA0447', 1, 2, 0, 6, 0 )
1159    ENDIF
1160    IF ( TRIM( approximation ) == 'anelastic'   .AND.                          &
1161         TRIM( psolver ) == 'multigrid' )  THEN
1162       message_string = 'Anelastic approximation currently only supports ' //  &
1163                        'psolver = "poisfft", ' //                             &
1164                        'psolver = "sor" and ' //                              &
1165                        'psolver = "multigrid_noopt"'
1166       CALL message( 'check_parameters', 'PA0448', 1, 2, 0, 6, 0 )
1167    ENDIF
1168    IF ( TRIM( approximation ) == 'anelastic'   .AND.                          &
1169         conserve_volume_flow )  THEN
1170       message_string = 'Anelastic approximation is not allowed with ' //      &
1171                        'conserve_volume_flow = .TRUE.'
1172       CALL message( 'check_parameters', 'PA0449', 1, 2, 0, 6, 0 )
1173    ENDIF
1174
1175!
1176!-- Check flux input mode
1177    IF ( TRIM( flux_input_mode ) /= 'dynamic'    .AND.                         &
1178         TRIM( flux_input_mode ) /= 'kinematic'  .AND.                         &
1179         TRIM( flux_input_mode ) /= 'approximation-specific' )  THEN
1180       message_string = 'unknown flux input mode: flux_input_mode = "' //      &
1181                        TRIM( flux_input_mode ) // '"'
1182       CALL message( 'check_parameters', 'PA0450', 1, 2, 0, 6, 0 )
1183    ENDIF
1184!
1185!-- Set flux input mode according to approximation if applicable
1186    IF ( TRIM( flux_input_mode ) == 'approximation-specific' )  THEN
1187       IF ( TRIM( approximation ) == 'anelastic' )  THEN
1188          flux_input_mode = 'dynamic'
1189       ELSEIF ( TRIM( approximation ) == 'boussinesq' )  THEN
1190          flux_input_mode = 'kinematic'
1191       ENDIF
1192    ENDIF
1193
1194!
1195!-- Check flux output mode
1196    IF ( TRIM( flux_output_mode ) /= 'dynamic'    .AND.                        &
1197         TRIM( flux_output_mode ) /= 'kinematic'  .AND.                        &
1198         TRIM( flux_output_mode ) /= 'approximation-specific' )  THEN
1199       message_string = 'unknown flux output mode: flux_output_mode = "' //    &
1200                        TRIM( flux_output_mode ) // '"'
1201       CALL message( 'check_parameters', 'PA0451', 1, 2, 0, 6, 0 )
1202    ENDIF
1203!
1204!-- Set flux output mode according to approximation if applicable
1205    IF ( TRIM( flux_output_mode ) == 'approximation-specific' )  THEN
1206       IF ( TRIM( approximation ) == 'anelastic' )  THEN
1207          flux_output_mode = 'dynamic'
1208       ELSEIF ( TRIM( approximation ) == 'boussinesq' )  THEN
1209          flux_output_mode = 'kinematic'
1210       ENDIF
1211    ENDIF
1212
1213
1214!
1215!-- When the land- or urban-surface model is used, the flux output must be
1216!-- dynamic.
1217    IF ( land_surface  .OR.  urban_surface )  THEN
1218       flux_output_mode = 'dynamic'
1219    ENDIF
1220
1221!
1222!-- Set the flux output units according to flux_output_mode
1223    IF ( TRIM( flux_output_mode ) == 'kinematic' ) THEN
1224        heatflux_output_unit              = 'K m/s'
1225        waterflux_output_unit             = 'kg/kg m/s'
1226        momentumflux_output_unit          = 'm2/s2'
1227    ELSEIF ( TRIM( flux_output_mode ) == 'dynamic' ) THEN
1228        heatflux_output_unit              = 'W/m2'
1229        waterflux_output_unit             = 'W/m2'
1230        momentumflux_output_unit          = 'N/m2'
1231    ENDIF
1232
1233!
1234!-- set time series output units for fluxes
1235    dots_unit(14:16) = heatflux_output_unit
1236    dots_unit(21)    = waterflux_output_unit
1237    dots_unit(19:20) = momentumflux_output_unit
1238
1239!
1240!-- Check whether there are any illegal values
1241!-- Pressure solver:
1242    IF ( psolver /= 'poisfft'  .AND.  psolver /= 'sor'  .AND.                  &
1243         psolver /= 'multigrid'  .AND.  psolver /= 'multigrid_noopt' )  THEN
1244       message_string = 'unknown solver for perturbation pressure: psolver' // &
1245                        ' = "' // TRIM( psolver ) // '"'
1246       CALL message( 'check_parameters', 'PA0016', 1, 2, 0, 6, 0 )
1247    ENDIF
1248
1249    IF ( psolver(1:9) == 'multigrid' )  THEN
1250       IF ( cycle_mg == 'w' )  THEN
1251          gamma_mg = 2
1252       ELSEIF ( cycle_mg == 'v' )  THEN
1253          gamma_mg = 1
1254       ELSE
1255          message_string = 'unknown multigrid cycle: cycle_mg = "' //          &
1256                           TRIM( cycle_mg ) // '"'
1257          CALL message( 'check_parameters', 'PA0020', 1, 2, 0, 6, 0 )
1258       ENDIF
1259    ENDIF
1260
1261    IF ( fft_method /= 'singleton-algorithm'  .AND.                            &
1262         fft_method /= 'temperton-algorithm'  .AND.                            &
1263         fft_method /= 'fftw'                 .AND.                            &
1264         fft_method /= 'system-specific' )  THEN
1265       message_string = 'unknown fft-algorithm: fft_method = "' //             &
1266                        TRIM( fft_method ) // '"'
1267       CALL message( 'check_parameters', 'PA0021', 1, 2, 0, 6, 0 )
1268    ENDIF
1269
1270    IF( momentum_advec == 'ws-scheme' .AND.                                    &
1271        .NOT. call_psolver_at_all_substeps  ) THEN
1272        message_string = 'psolver must be called at each RK3 substep when "'// &
1273                      TRIM(momentum_advec) // ' "is used for momentum_advec'
1274        CALL message( 'check_parameters', 'PA0344', 1, 2, 0, 6, 0 )
1275    END IF
1276!
1277!-- Advection schemes:
1278    IF ( momentum_advec /= 'pw-scheme'  .AND.                                  & 
1279         momentum_advec /= 'ws-scheme'  .AND.                                  &
1280         momentum_advec /= 'up-scheme' )                                       &
1281    THEN
1282       message_string = 'unknown advection scheme: momentum_advec = "' //      &
1283                        TRIM( momentum_advec ) // '"'
1284       CALL message( 'check_parameters', 'PA0022', 1, 2, 0, 6, 0 )
1285    ENDIF
1286    IF ( ( momentum_advec == 'ws-scheme' .OR.  scalar_advec == 'ws-scheme' )   &
1287           .AND. ( timestep_scheme == 'euler' .OR.                             &
1288                   timestep_scheme == 'runge-kutta-2' ) )                      &
1289    THEN
1290       message_string = 'momentum_advec or scalar_advec = "'                   &
1291         // TRIM( momentum_advec ) // '" is not allowed with ' //              &
1292         'timestep_scheme = "' // TRIM( timestep_scheme ) // '"'
1293       CALL message( 'check_parameters', 'PA0023', 1, 2, 0, 6, 0 )
1294    ENDIF
1295    IF ( scalar_advec /= 'pw-scheme'  .AND.  scalar_advec /= 'ws-scheme' .AND. &
1296         scalar_advec /= 'bc-scheme' .AND. scalar_advec /= 'up-scheme' )       &
1297    THEN
1298       message_string = 'unknown advection scheme: scalar_advec = "' //        &
1299                        TRIM( scalar_advec ) // '"'
1300       CALL message( 'check_parameters', 'PA0024', 1, 2, 0, 6, 0 )
1301    ENDIF
1302    IF ( scalar_advec == 'bc-scheme'  .AND.  loop_optimization == 'cache' )    &
1303    THEN
1304       message_string = 'advection_scheme scalar_advec = "'                    &
1305         // TRIM( scalar_advec ) // '" not implemented for ' //                &
1306         'loop_optimization = "' // TRIM( loop_optimization ) // '"'
1307       CALL message( 'check_parameters', 'PA0026', 1, 2, 0, 6, 0 )
1308    ENDIF
1309
1310    IF ( use_sgs_for_particles  .AND.  .NOT. cloud_droplets  .AND.             &
1311         .NOT. use_upstream_for_tke  .AND.                                     &
1312         scalar_advec /= 'ws-scheme'                                           &
1313       )  THEN
1314       use_upstream_for_tke = .TRUE.
1315       message_string = 'use_upstream_for_tke is set to .TRUE. because ' //    &
1316                        'use_sgs_for_particles = .TRUE. '          //          &
1317                        'and scalar_advec /= ws-scheme'
1318       CALL message( 'check_parameters', 'PA0025', 0, 1, 0, 6, 0 )
1319    ENDIF
1320
1321!
1322!-- Set LOGICAL switches to enhance performance
1323    IF ( momentum_advec == 'ws-scheme' )  ws_scheme_mom = .TRUE.
1324    IF ( scalar_advec   == 'ws-scheme' )  ws_scheme_sca = .TRUE.
1325
1326
1327!
1328!-- Timestep schemes:
1329    SELECT CASE ( TRIM( timestep_scheme ) )
1330
1331       CASE ( 'euler' )
1332          intermediate_timestep_count_max = 1
1333
1334       CASE ( 'runge-kutta-2' )
1335          intermediate_timestep_count_max = 2
1336
1337       CASE ( 'runge-kutta-3' )
1338          intermediate_timestep_count_max = 3
1339
1340       CASE DEFAULT
1341          message_string = 'unknown timestep scheme: timestep_scheme = "' //   &
1342                           TRIM( timestep_scheme ) // '"'
1343          CALL message( 'check_parameters', 'PA0027', 1, 2, 0, 6, 0 )
1344
1345    END SELECT
1346
1347    IF ( (momentum_advec /= 'pw-scheme' .AND. momentum_advec /= 'ws-scheme')   &
1348         .AND. timestep_scheme(1:5) == 'runge' ) THEN
1349       message_string = 'momentum advection scheme "' // &
1350                        TRIM( momentum_advec ) // '" & does not work with ' // &
1351                        'timestep_scheme "' // TRIM( timestep_scheme ) // '"'
1352       CALL message( 'check_parameters', 'PA0029', 1, 2, 0, 6, 0 )
1353    ENDIF
1354!
1355!-- Check for proper settings for microphysics
1356    IF ( bulk_cloud_model  .AND.  cloud_droplets )  THEN
1357       message_string = 'bulk_cloud_model = .TRUE. is not allowed with ' //    &
1358                        'cloud_droplets = .TRUE.'
1359       CALL message( 'check_parameters', 'PA0442', 1, 2, 0, 6, 0 )
1360    ENDIF
1361
1362!
1363!-- Collision kernels:
1364    SELECT CASE ( TRIM( collision_kernel ) )
1365
1366       CASE ( 'hall', 'hall_fast' )
1367          hall_kernel = .TRUE.
1368
1369       CASE ( 'wang', 'wang_fast' )
1370          wang_kernel = .TRUE.
1371
1372       CASE ( 'none' )
1373
1374
1375       CASE DEFAULT
1376          message_string = 'unknown collision kernel: collision_kernel = "' // &
1377                           TRIM( collision_kernel ) // '"'
1378          CALL message( 'check_parameters', 'PA0350', 1, 2, 0, 6, 0 )
1379
1380    END SELECT
1381    IF ( collision_kernel(6:9) == 'fast' )  use_kernel_tables = .TRUE.
1382
1383!
1384!-- Initializing actions must have been set by the user
1385    IF ( TRIM( initializing_actions ) == '' )  THEN
1386       message_string = 'no value specified for initializing_actions'
1387       CALL message( 'check_parameters', 'PA0149', 1, 2, 0, 6, 0 )
1388    ENDIF
1389
1390    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
1391         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
1392!
1393!--    No restart run: several initialising actions are possible
1394       action = initializing_actions
1395       DO  WHILE ( TRIM( action ) /= '' )
1396          position = INDEX( action, ' ' )
1397          SELECT CASE ( action(1:position-1) )
1398
1399             CASE ( 'set_constant_profiles', 'set_1d-model_profiles',          &
1400                    'by_user', 'initialize_vortex', 'initialize_ptanom',       &
1401                    'initialize_bubble', 'inifor' )
1402                action = action(position+1:)
1403
1404             CASE DEFAULT
1405                message_string = 'initializing_action = "' //                  &
1406                                 TRIM( action ) // '" unknown or not allowed'
1407                CALL message( 'check_parameters', 'PA0030', 1, 2, 0, 6, 0 )
1408
1409          END SELECT
1410       ENDDO
1411    ENDIF
1412
1413    IF ( TRIM( initializing_actions ) == 'initialize_vortex'  .AND.            &
1414         conserve_volume_flow ) THEN
1415         message_string = 'initializing_actions = "initialize_vortex"' //      &
1416                        ' is not allowed with conserve_volume_flow = .T.'
1417       CALL message( 'check_parameters', 'PA0343', 1, 2, 0, 6, 0 )
1418    ENDIF
1419
1420
1421    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.    &
1422         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
1423       message_string = 'initializing_actions = "set_constant_profiles"' //    &
1424                        ' and "set_1d-model_profiles" are not allowed ' //     &
1425                        'simultaneously'
1426       CALL message( 'check_parameters', 'PA0031', 1, 2, 0, 6, 0 )
1427    ENDIF
1428
1429    IF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0  .AND.    &
1430         INDEX( initializing_actions, 'by_user' ) /= 0 )  THEN
1431       message_string = 'initializing_actions = "set_constant_profiles"' //    &
1432                        ' and "by_user" are not allowed simultaneously'
1433       CALL message( 'check_parameters', 'PA0032', 1, 2, 0, 6, 0 )
1434    ENDIF
1435
1436    IF ( INDEX( initializing_actions, 'by_user' ) /= 0  .AND.                  &
1437         INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
1438       message_string = 'initializing_actions = "by_user" and ' //             &
1439                        '"set_1d-model_profiles" are not allowed simultaneously'
1440       CALL message( 'check_parameters', 'PA0033', 1, 2, 0, 6, 0 )
1441    ENDIF
1442!
1443!-- In case of spinup and nested run, spinup end time must be identical
1444!-- in order to have synchronously running simulations.
1445    IF ( nested_run )  THEN
1446#if defined( __parallel )
1447       CALL MPI_ALLREDUCE( spinup_time, spinup_time_max, 1, MPI_REAL,          &
1448                           MPI_MAX, MPI_COMM_WORLD, ierr )
1449       CALL MPI_ALLREDUCE( dt_spinup,   dt_spinup_max,   1, MPI_REAL,          &
1450                           MPI_MAX, MPI_COMM_WORLD, ierr )
1451
1452       IF ( spinup_time /= spinup_time_max  .OR.  dt_spinup /= dt_spinup_max ) &
1453       THEN
1454          message_string = 'In case of nesting, spinup_time and ' //           &
1455                           'dt_spinup must be identical in all parent ' //     &
1456                           'and child domains.'
1457          CALL message( 'check_parameters', 'PA0489', 3, 2, 0, 6, 0 )
1458       ENDIF
1459#endif
1460    ENDIF
1461
1462    IF ( bulk_cloud_model  .AND.  .NOT.  humidity )  THEN
1463       WRITE( message_string, * ) 'bulk_cloud_model = ', bulk_cloud_model,     &
1464              ' is not allowed with humidity = ', humidity
1465       CALL message( 'check_parameters', 'PA0034', 1, 2, 0, 6, 0 )
1466    ENDIF
1467
1468    IF ( humidity  .AND.  sloping_surface )  THEN
1469       message_string = 'humidity = .TRUE. and sloping_surface = .TRUE. ' //   &
1470                        'are not allowed simultaneously'
1471       CALL message( 'check_parameters', 'PA0036', 1, 2, 0, 6, 0 )
1472    ENDIF
1473
1474!   Check for chem_emission_mod parameters setting
1475!   IF ( air_chemistry ) CALL chem_emissions_check_parameters   ! forkel preliminary
1476
1477
1478!   Check for chemitry_model_mod parameters setting
1479    IF ( air_chemistry )  CALL chem_check_parameters
1480
1481
1482!-- Check the module settings
1483    IF ( bulk_cloud_model )     CALL bcm_check_parameters
1484    IF ( gust_module_enabled )  CALL gust_check_parameters
1485    IF ( large_scale_forcing  .OR.  nudging )                                  &
1486                                CALL lsf_nudging_check_parameters
1487    IF ( land_surface )         CALL lsm_check_parameters
1488    IF ( ocean_mode )           CALL ocean_check_parameters
1489    IF ( plant_canopy )         CALL pcm_check_parameters
1490    IF ( radiation )            CALL radiation_check_parameters
1491    IF ( calculate_spectra )    CALL spectra_check_parameters
1492    CALL stg_check_parameters
1493    CALL tcm_check_parameters
1494    IF ( urban_surface )        CALL usm_check_parameters
1495    IF ( wind_turbine )         CALL wtm_check_parameters
1496
1497!
1498!-- In case of no restart run, check initialising parameters and calculate
1499!-- further quantities
1500    IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
1501
1502!
1503!--    Initial profiles for 1D and 3D model, respectively (u,v further below)
1504       pt_init = pt_surface
1505       IF ( humidity       )  q_init  = q_surface
1506       IF ( passive_scalar )  s_init  = s_surface
1507!
1508!--    TODO
1509!--    Russo: Is done in chem_init and would overwrite what is done there
1510!--    --> kanani: Revise
1511!      IF ( air_chemistry )  THEN
1512!        DO lsp = 1, nvar
1513!           chem_species(lsp)%conc_pr_init = cs_surface(lsp)     
1514!        ENDDO
1515!      ENDIF
1516!
1517!--
1518!--    If required, compute initial profile of the geostrophic wind
1519!--    (component ug)
1520       i = 1
1521       gradient = 0.0_wp
1522
1523       IF ( .NOT. ocean_mode )  THEN
1524
1525          ug_vertical_gradient_level_ind(1) = 0
1526          ug(0) = ug_surface
1527          DO  k = 1, nzt+1
1528             IF ( i < 11 )  THEN
1529                IF ( ug_vertical_gradient_level(i) < zu(k)  .AND.              &
1530                     ug_vertical_gradient_level(i) >= 0.0_wp )  THEN
1531                   gradient = ug_vertical_gradient(i) / 100.0_wp
1532                   ug_vertical_gradient_level_ind(i) = k - 1
1533                   i = i + 1
1534                ENDIF
1535             ENDIF
1536             IF ( gradient /= 0.0_wp )  THEN
1537                IF ( k /= 1 )  THEN
1538                   ug(k) = ug(k-1) + dzu(k) * gradient
1539                ELSE
1540                   ug(k) = ug_surface + dzu(k) * gradient
1541                ENDIF
1542             ELSE
1543                ug(k) = ug(k-1)
1544             ENDIF
1545          ENDDO
1546
1547       ELSE
1548
1549          ug_vertical_gradient_level_ind(1) = nzt+1
1550          ug(nzt+1) = ug_surface
1551          DO  k = nzt, nzb, -1
1552             IF ( i < 11 )  THEN
1553                IF ( ug_vertical_gradient_level(i) > zu(k)  .AND.              &
1554                     ug_vertical_gradient_level(i) <= 0.0_wp )  THEN
1555                   gradient = ug_vertical_gradient(i) / 100.0_wp
1556                   ug_vertical_gradient_level_ind(i) = k + 1
1557                   i = i + 1
1558                ENDIF
1559             ENDIF
1560             IF ( gradient /= 0.0_wp )  THEN
1561                IF ( k /= nzt )  THEN
1562                   ug(k) = ug(k+1) - dzu(k+1) * gradient
1563                ELSE
1564                   ug(k)   = ug_surface - 0.5_wp * dzu(k+1) * gradient
1565                   ug(k+1) = ug_surface + 0.5_wp * dzu(k+1) * gradient
1566                ENDIF
1567             ELSE
1568                ug(k) = ug(k+1)
1569             ENDIF
1570          ENDDO
1571
1572       ENDIF
1573
1574!
1575!--    In case of no given gradients for ug, choose a zero gradient
1576       IF ( ug_vertical_gradient_level(1) == -9999999.9_wp )  THEN
1577          ug_vertical_gradient_level(1) = 0.0_wp
1578       ENDIF
1579
1580!
1581!--
1582!--    If required, compute initial profile of the geostrophic wind
1583!--    (component vg)
1584       i = 1
1585       gradient = 0.0_wp
1586
1587       IF ( .NOT. ocean_mode )  THEN
1588
1589          vg_vertical_gradient_level_ind(1) = 0
1590          vg(0) = vg_surface
1591          DO  k = 1, nzt+1
1592             IF ( i < 11 )  THEN
1593                IF ( vg_vertical_gradient_level(i) < zu(k)  .AND.              &
1594                     vg_vertical_gradient_level(i) >= 0.0_wp )  THEN
1595                   gradient = vg_vertical_gradient(i) / 100.0_wp
1596                   vg_vertical_gradient_level_ind(i) = k - 1
1597                   i = i + 1
1598                ENDIF
1599             ENDIF
1600             IF ( gradient /= 0.0_wp )  THEN
1601                IF ( k /= 1 )  THEN
1602                   vg(k) = vg(k-1) + dzu(k) * gradient
1603                ELSE
1604                   vg(k) = vg_surface + dzu(k) * gradient
1605                ENDIF
1606             ELSE
1607                vg(k) = vg(k-1)
1608             ENDIF
1609          ENDDO
1610
1611       ELSE
1612
1613          vg_vertical_gradient_level_ind(1) = nzt+1
1614          vg(nzt+1) = vg_surface
1615          DO  k = nzt, nzb, -1
1616             IF ( i < 11 )  THEN
1617                IF ( vg_vertical_gradient_level(i) > zu(k)  .AND.              &
1618                     vg_vertical_gradient_level(i) <= 0.0_wp )  THEN
1619                   gradient = vg_vertical_gradient(i) / 100.0_wp
1620                   vg_vertical_gradient_level_ind(i) = k + 1
1621                   i = i + 1
1622                ENDIF
1623             ENDIF
1624             IF ( gradient /= 0.0_wp )  THEN
1625                IF ( k /= nzt )  THEN
1626                   vg(k) = vg(k+1) - dzu(k+1) * gradient
1627                ELSE
1628                   vg(k)   = vg_surface - 0.5_wp * dzu(k+1) * gradient
1629                   vg(k+1) = vg_surface + 0.5_wp * dzu(k+1) * gradient
1630                ENDIF
1631             ELSE
1632                vg(k) = vg(k+1)
1633             ENDIF
1634          ENDDO
1635
1636       ENDIF
1637
1638!
1639!--    In case of no given gradients for vg, choose a zero gradient
1640       IF ( vg_vertical_gradient_level(1) == -9999999.9_wp )  THEN
1641          vg_vertical_gradient_level(1) = 0.0_wp
1642       ENDIF
1643
1644!
1645!--    Let the initial wind profiles be the calculated ug/vg profiles or
1646!--    interpolate them from wind profile data (if given)
1647       IF ( u_profile(1) == 9999999.9_wp  .AND.  v_profile(1) == 9999999.9_wp )  THEN
1648
1649          u_init = ug
1650          v_init = vg
1651
1652       ELSEIF ( u_profile(1) == 0.0_wp  .AND.  v_profile(1) == 0.0_wp )  THEN
1653
1654          IF ( uv_heights(1) /= 0.0_wp )  THEN
1655             message_string = 'uv_heights(1) must be 0.0'
1656             CALL message( 'check_parameters', 'PA0345', 1, 2, 0, 6, 0 )
1657          ENDIF
1658
1659          IF ( omega /= 0.0_wp )  THEN
1660             message_string = 'Coriolis force must be switched off (by setting omega=0.0)' //  &
1661                              ' when prescribing the forcing by u_profile and v_profile'
1662             CALL message( 'check_parameters', 'PA0347', 1, 2, 0, 6, 0 )
1663          ENDIF
1664
1665          use_prescribed_profile_data = .TRUE.
1666
1667          kk = 1
1668          u_init(0) = 0.0_wp
1669          v_init(0) = 0.0_wp
1670
1671          DO  k = 1, nz+1
1672
1673             IF ( kk < 200 )  THEN
1674                DO  WHILE ( uv_heights(kk+1) <= zu(k) )
1675                   kk = kk + 1
1676                   IF ( kk == 200 )  EXIT
1677                ENDDO
1678             ENDIF
1679
1680             IF ( kk < 200  .AND.  uv_heights(kk+1) /= 9999999.9_wp )  THEN
1681                u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
1682                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
1683                                       ( u_profile(kk+1) - u_profile(kk) )
1684                v_init(k) = v_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
1685                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
1686                                       ( v_profile(kk+1) - v_profile(kk) )
1687             ELSE
1688                u_init(k) = u_profile(kk)
1689                v_init(k) = v_profile(kk)
1690             ENDIF
1691
1692          ENDDO
1693
1694       ELSE
1695
1696          message_string = 'u_profile(1) and v_profile(1) must be 0.0'
1697          CALL message( 'check_parameters', 'PA0346', 1, 2, 0, 6, 0 )
1698
1699       ENDIF
1700
1701!
1702!--    Compute initial temperature profile using the given temperature gradients
1703       IF (  .NOT.  neutral )  THEN
1704          CALL init_vertical_profiles( pt_vertical_gradient_level_ind,          &
1705                                       pt_vertical_gradient_level,              &
1706                                       pt_vertical_gradient, pt_init,           &
1707                                       pt_surface, bc_pt_t_val )
1708       ENDIF
1709!
1710!--    Compute initial humidity profile using the given humidity gradients
1711       IF ( humidity )  THEN
1712          CALL init_vertical_profiles( q_vertical_gradient_level_ind,          &
1713                                       q_vertical_gradient_level,              &
1714                                       q_vertical_gradient, q_init,            &
1715                                       q_surface, bc_q_t_val )
1716       ENDIF
1717!
1718!--    Compute initial scalar profile using the given scalar gradients
1719       IF ( passive_scalar )  THEN
1720          CALL init_vertical_profiles( s_vertical_gradient_level_ind,          &
1721                                       s_vertical_gradient_level,              &
1722                                       s_vertical_gradient, s_init,            &
1723                                       s_surface, bc_s_t_val )
1724       ENDIF
1725!
1726!--    TODO
1727!--    Compute initial chemistry profile using the given chemical species gradients
1728!--    Russo: Is done in chem_init --> kanani: Revise
1729
1730    ENDIF
1731
1732!
1733!-- Check if the control parameter use_subsidence_tendencies is used correctly
1734    IF ( use_subsidence_tendencies  .AND.  .NOT.  large_scale_subsidence )  THEN
1735       message_string = 'The usage of use_subsidence_tendencies ' //           &
1736                            'requires large_scale_subsidence = .T..'
1737       CALL message( 'check_parameters', 'PA0396', 1, 2, 0, 6, 0 )
1738    ELSEIF ( use_subsidence_tendencies  .AND.  .NOT. large_scale_forcing )  THEN
1739       message_string = 'The usage of use_subsidence_tendencies ' //           &
1740                            'requires large_scale_forcing = .T..'
1741       CALL message( 'check_parameters', 'PA0397', 1, 2, 0, 6, 0 )
1742    ENDIF
1743
1744!
1745!-- Initialize large scale subsidence if required
1746    If ( large_scale_subsidence )  THEN
1747       IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp  .AND.            &
1748                                     .NOT.  large_scale_forcing )  THEN
1749          CALL init_w_subsidence
1750       ENDIF
1751!
1752!--    In case large_scale_forcing is used, profiles for subsidence velocity
1753!--    are read in from file LSF_DATA
1754
1755       IF ( subs_vertical_gradient_level(1) == -9999999.9_wp  .AND.            &
1756            .NOT.  large_scale_forcing )  THEN
1757          message_string = 'There is no default large scale vertical ' //      &
1758                           'velocity profile set. Specify the subsidence ' //  &
1759                           'velocity profile via subs_vertical_gradient ' //   &
1760                           'and subs_vertical_gradient_level.'
1761          CALL message( 'check_parameters', 'PA0380', 1, 2, 0, 6, 0 )
1762       ENDIF
1763    ELSE
1764        IF ( subs_vertical_gradient_level(1) /= -9999999.9_wp )  THEN
1765           message_string = 'Enable usage of large scale subsidence by ' //    &
1766                            'setting large_scale_subsidence = .T..'
1767          CALL message( 'check_parameters', 'PA0381', 1, 2, 0, 6, 0 )
1768        ENDIF
1769    ENDIF
1770
1771!
1772!-- Overwrite latitude if necessary and compute Coriolis parameter.
1773!-- @todo - move initialization of f and fs to coriolis_mod.
1774    IF ( input_pids_static )  THEN
1775       latitude  = init_model%latitude
1776       longitude = init_model%longitude
1777    ENDIF
1778
1779    f  = 2.0_wp * omega * SIN( latitude / 180.0_wp * pi )
1780    fs = 2.0_wp * omega * COS( latitude / 180.0_wp * pi )
1781
1782!
1783!-- Check and set buoyancy related parameters and switches
1784    IF ( reference_state == 'horizontal_average' )  THEN
1785       CONTINUE
1786    ELSEIF ( reference_state == 'initial_profile' )  THEN
1787       use_initial_profile_as_reference = .TRUE.
1788    ELSEIF ( reference_state == 'single_value' )  THEN
1789       use_single_reference_value = .TRUE.
1790       IF ( pt_reference == 9999999.9_wp )  pt_reference = pt_surface
1791       vpt_reference = pt_reference * ( 1.0_wp + 0.61_wp * q_surface )
1792    ELSE
1793       message_string = 'illegal value for reference_state: "' //              &
1794                        TRIM( reference_state ) // '"'
1795       CALL message( 'check_parameters', 'PA0056', 1, 2, 0, 6, 0 )
1796    ENDIF
1797
1798!
1799!-- In case of a given slope, compute the relevant quantities
1800    IF ( alpha_surface /= 0.0_wp )  THEN
1801       IF ( ABS( alpha_surface ) > 90.0_wp )  THEN
1802          WRITE( message_string, * ) 'ABS( alpha_surface = ', alpha_surface,   &
1803                                     ' ) must be < 90.0'
1804          CALL message( 'check_parameters', 'PA0043', 1, 2, 0, 6, 0 )
1805       ENDIF
1806       sloping_surface = .TRUE.
1807       cos_alpha_surface = COS( alpha_surface / 180.0_wp * pi )
1808       sin_alpha_surface = SIN( alpha_surface / 180.0_wp * pi )
1809    ENDIF
1810
1811!
1812!-- Check time step and cfl_factor
1813    IF ( dt /= -1.0_wp )  THEN
1814       IF ( dt <= 0.0_wp )  THEN
1815          WRITE( message_string, * ) 'dt = ', dt , ' <= 0.0'
1816          CALL message( 'check_parameters', 'PA0044', 1, 2, 0, 6, 0 )
1817       ENDIF
1818       dt_3d = dt
1819       dt_fixed = .TRUE.
1820    ENDIF
1821
1822    IF ( cfl_factor <= 0.0_wp  .OR.  cfl_factor > 1.0_wp )  THEN
1823       IF ( cfl_factor == -1.0_wp )  THEN
1824          IF ( timestep_scheme == 'runge-kutta-2' )  THEN
1825             cfl_factor = 0.8_wp
1826          ELSEIF ( timestep_scheme == 'runge-kutta-3' )  THEN
1827             cfl_factor = 0.9_wp
1828          ELSE
1829             cfl_factor = 0.9_wp
1830          ENDIF
1831       ELSE
1832          WRITE( message_string, * ) 'cfl_factor = ', cfl_factor,              &
1833                 ' out of range &0.0 < cfl_factor <= 1.0 is required'
1834          CALL message( 'check_parameters', 'PA0045', 1, 2, 0, 6, 0 )
1835       ENDIF
1836    ENDIF
1837
1838!
1839!-- Store simulated time at begin
1840    simulated_time_at_begin = simulated_time
1841
1842!
1843!-- Store reference time for coupled runs and change the coupling flag,
1844!-- if ...
1845    IF ( simulated_time == 0.0_wp )  THEN
1846       IF ( coupling_start_time == 0.0_wp )  THEN
1847          time_since_reference_point = 0.0_wp
1848       ELSEIF ( time_since_reference_point < 0.0_wp )  THEN
1849          run_coupled = .FALSE.
1850       ENDIF
1851    ENDIF
1852
1853!
1854!-- Set wind speed in the Galilei-transformed system
1855    IF ( galilei_transformation )  THEN
1856       IF ( use_ug_for_galilei_tr                    .AND.                     &
1857            ug_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
1858            ug_vertical_gradient(1) == 0.0_wp        .AND.                     &
1859            vg_vertical_gradient_level(1) == 0.0_wp  .AND.                     &
1860            vg_vertical_gradient(1) == 0.0_wp )  THEN
1861          u_gtrans = ug_surface * 0.6_wp
1862          v_gtrans = vg_surface * 0.6_wp
1863       ELSEIF ( use_ug_for_galilei_tr  .AND.                                   &
1864                ( ug_vertical_gradient_level(1) /= 0.0_wp  .OR.                &
1865                ug_vertical_gradient(1) /= 0.0_wp ) )  THEN
1866          message_string = 'baroclinity (ug) not allowed simultaneously' //    &
1867                           ' with galilei transformation'
1868          CALL message( 'check_parameters', 'PA0046', 1, 2, 0, 6, 0 )
1869       ELSEIF ( use_ug_for_galilei_tr  .AND.                                   &
1870                ( vg_vertical_gradient_level(1) /= 0.0_wp  .OR.                &
1871                vg_vertical_gradient(1) /= 0.0_wp ) )  THEN
1872          message_string = 'baroclinity (vg) not allowed simultaneously' //    &
1873                           ' with galilei transformation'
1874          CALL message( 'check_parameters', 'PA0047', 1, 2, 0, 6, 0 )
1875       ELSE
1876          message_string = 'variable translation speed used for Galilei-' //   &
1877             'transformation, which may cause & instabilities in stably ' //   &
1878             'stratified regions'
1879          CALL message( 'check_parameters', 'PA0048', 0, 1, 0, 6, 0 )
1880       ENDIF
1881    ENDIF
1882
1883!
1884!-- In case of using a prandtl-layer, calculated (or prescribed) surface
1885!-- fluxes have to be used in the diffusion-terms
1886    IF ( constant_flux_layer )  use_surface_fluxes = .TRUE.
1887
1888!
1889!-- Check boundary conditions and set internal variables:
1890!-- Attention: the lateral boundary conditions have been already checked in
1891!-- parin
1892!
1893!-- Non-cyclic lateral boundaries require the multigrid method and Piascek-
1894!-- Willimas or Wicker - Skamarock advection scheme. Several schemes
1895!-- and tools do not work with non-cyclic boundary conditions.
1896    IF ( bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic' )  THEN
1897       IF ( psolver(1:9) /= 'multigrid' )  THEN
1898          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1899                           'psolver = "' // TRIM( psolver ) // '"'
1900          CALL message( 'check_parameters', 'PA0051', 1, 2, 0, 6, 0 )
1901       ENDIF
1902       IF ( momentum_advec /= 'pw-scheme'  .AND.                               &
1903            momentum_advec /= 'ws-scheme' )  THEN
1904
1905          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1906                           'momentum_advec = "' // TRIM( momentum_advec ) // '"'
1907          CALL message( 'check_parameters', 'PA0052', 1, 2, 0, 6, 0 )
1908       ENDIF
1909       IF ( scalar_advec /= 'pw-scheme'  .AND.                                 &
1910            scalar_advec /= 'ws-scheme' )  THEN
1911          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1912                           'scalar_advec = "' // TRIM( scalar_advec ) // '"'
1913          CALL message( 'check_parameters', 'PA0053', 1, 2, 0, 6, 0 )
1914       ENDIF
1915       IF ( galilei_transformation )  THEN
1916          message_string = 'non-cyclic lateral boundaries do not allow ' //    &
1917                           'galilei_transformation = .T.'
1918          CALL message( 'check_parameters', 'PA0054', 1, 2, 0, 6, 0 )
1919       ENDIF
1920    ENDIF
1921
1922!
1923!-- Bottom boundary condition for the turbulent Kinetic energy
1924    IF ( bc_e_b == 'neumann' )  THEN
1925       ibc_e_b = 1
1926    ELSEIF ( bc_e_b == '(u*)**2+neumann' )  THEN
1927       ibc_e_b = 2
1928       IF ( .NOT. constant_flux_layer )  THEN
1929          bc_e_b = 'neumann'
1930          ibc_e_b = 1
1931          message_string = 'boundary condition bc_e_b changed to "' //         &
1932                           TRIM( bc_e_b ) // '"'
1933          CALL message( 'check_parameters', 'PA0057', 0, 1, 0, 6, 0 )
1934       ENDIF
1935    ELSE
1936       message_string = 'unknown boundary condition: bc_e_b = "' //            &
1937                        TRIM( bc_e_b ) // '"'
1938       CALL message( 'check_parameters', 'PA0058', 1, 2, 0, 6, 0 )
1939    ENDIF
1940
1941!
1942!-- Boundary conditions for perturbation pressure
1943    IF ( bc_p_b == 'dirichlet' )  THEN
1944       ibc_p_b = 0
1945    ELSEIF ( bc_p_b == 'neumann' )  THEN
1946       ibc_p_b = 1
1947    ELSE
1948       message_string = 'unknown boundary condition: bc_p_b = "' //            &
1949                        TRIM( bc_p_b ) // '"'
1950       CALL message( 'check_parameters', 'PA0059', 1, 2, 0, 6, 0 )
1951    ENDIF
1952
1953    IF ( bc_p_t == 'dirichlet' )  THEN
1954       ibc_p_t = 0
1955!-- TO_DO: later set bc_p_t to neumann before, in case of nested domain
1956    ELSEIF ( bc_p_t == 'neumann' .OR. bc_p_t == 'nested' )  THEN
1957       ibc_p_t = 1
1958    ELSE
1959       message_string = 'unknown boundary condition: bc_p_t = "' //            &
1960                        TRIM( bc_p_t ) // '"'
1961       CALL message( 'check_parameters', 'PA0061', 1, 2, 0, 6, 0 )
1962    ENDIF
1963
1964!
1965!-- Boundary conditions for potential temperature
1966    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
1967       ibc_pt_b = 2
1968    ELSE
1969       IF ( bc_pt_b == 'dirichlet' )  THEN
1970          ibc_pt_b = 0
1971       ELSEIF ( bc_pt_b == 'neumann' )  THEN
1972          ibc_pt_b = 1
1973       ELSE
1974          message_string = 'unknown boundary condition: bc_pt_b = "' //        &
1975                           TRIM( bc_pt_b ) // '"'
1976          CALL message( 'check_parameters', 'PA0062', 1, 2, 0, 6, 0 )
1977       ENDIF
1978    ENDIF
1979
1980    IF ( bc_pt_t == 'dirichlet' )  THEN
1981       ibc_pt_t = 0
1982    ELSEIF ( bc_pt_t == 'neumann' )  THEN
1983       ibc_pt_t = 1
1984    ELSEIF ( bc_pt_t == 'initial_gradient' )  THEN
1985       ibc_pt_t = 2
1986    ELSEIF ( bc_pt_t == 'nested'  .OR.  bc_pt_t == 'nesting_offline' )  THEN
1987       ibc_pt_t = 3
1988    ELSE
1989       message_string = 'unknown boundary condition: bc_pt_t = "' //           &
1990                        TRIM( bc_pt_t ) // '"'
1991       CALL message( 'check_parameters', 'PA0063', 1, 2, 0, 6, 0 )
1992    ENDIF
1993
1994    IF ( ANY( wall_heatflux /= 0.0_wp )  .AND.                        &
1995         surface_heatflux == 9999999.9_wp )  THEN
1996       message_string = 'wall_heatflux additionally requires ' //     &
1997                        'setting of surface_heatflux'
1998       CALL message( 'check_parameters', 'PA0443', 1, 2, 0, 6, 0 )
1999    ENDIF
2000
2001!
2002!   This IF clause needs revision, got too complex!!
2003    IF ( surface_heatflux == 9999999.9_wp  )  THEN
2004       constant_heatflux = .FALSE.
2005       IF ( large_scale_forcing  .OR.  land_surface  .OR.  urban_surface )  THEN
2006          IF ( ibc_pt_b == 0 )  THEN
2007             constant_heatflux = .FALSE.
2008          ELSEIF ( ibc_pt_b == 1 )  THEN
2009             constant_heatflux = .TRUE.
2010             surface_heatflux = 0.0_wp
2011          ENDIF
2012       ENDIF
2013    ELSE
2014       constant_heatflux = .TRUE.
2015    ENDIF
2016
2017    IF ( top_heatflux     == 9999999.9_wp )  constant_top_heatflux = .FALSE.
2018
2019    IF ( neutral )  THEN
2020
2021       IF ( surface_heatflux /= 0.0_wp  .AND.                                  &
2022            surface_heatflux /= 9999999.9_wp )  THEN
2023          message_string = 'heatflux must not be set for pure neutral flow'
2024          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
2025       ENDIF
2026
2027       IF ( top_heatflux /= 0.0_wp  .AND.  top_heatflux /= 9999999.9_wp )      &
2028       THEN
2029          message_string = 'heatflux must not be set for pure neutral flow'
2030          CALL message( 'check_parameters', 'PA0351', 1, 2, 0, 6, 0 )
2031       ENDIF
2032
2033    ENDIF
2034
2035    IF ( top_momentumflux_u /= 9999999.9_wp  .AND.                             &
2036         top_momentumflux_v /= 9999999.9_wp )  THEN
2037       constant_top_momentumflux = .TRUE.
2038    ELSEIF (  .NOT. ( top_momentumflux_u == 9999999.9_wp  .AND.                &
2039           top_momentumflux_v == 9999999.9_wp ) )  THEN
2040       message_string = 'both, top_momentumflux_u AND top_momentumflux_v ' //  &
2041                        'must be set'
2042       CALL message( 'check_parameters', 'PA0064', 1, 2, 0, 6, 0 )
2043    ENDIF
2044
2045!
2046!-- A given surface temperature implies Dirichlet boundary condition for
2047!-- temperature. In this case specification of a constant heat flux is
2048!-- forbidden.
2049    IF ( ibc_pt_b == 0  .AND.  constant_heatflux  .AND.                        &
2050         surface_heatflux /= 0.0_wp )  THEN
2051       message_string = 'boundary_condition: bc_pt_b = "' // TRIM( bc_pt_b ) //&
2052                        '& is not allowed with constant_heatflux = .TRUE.'
2053       CALL message( 'check_parameters', 'PA0065', 1, 2, 0, 6, 0 )
2054    ENDIF
2055    IF ( constant_heatflux  .AND.  pt_surface_initial_change /= 0.0_wp )  THEN
2056       WRITE ( message_string, * )  'constant_heatflux = .TRUE. is not allo',  &
2057               'wed with pt_surface_initial_change (/=0) = ',                  &
2058               pt_surface_initial_change
2059       CALL message( 'check_parameters', 'PA0066', 1, 2, 0, 6, 0 )
2060    ENDIF
2061
2062!
2063!-- A given temperature at the top implies Dirichlet boundary condition for
2064!-- temperature. In this case specification of a constant heat flux is
2065!-- forbidden.
2066    IF ( ibc_pt_t == 0  .AND.  constant_top_heatflux  .AND.                    &
2067         top_heatflux /= 0.0_wp )  THEN
2068       message_string = 'boundary_condition: bc_pt_t = "' // TRIM( bc_pt_t ) //&
2069                        '" is not allowed with constant_top_heatflux = .TRUE.'
2070       CALL message( 'check_parameters', 'PA0067', 1, 2, 0, 6, 0 )
2071    ENDIF
2072
2073!
2074!-- Set boundary conditions for total water content
2075    IF ( humidity )  THEN
2076
2077       IF ( ANY( wall_humidityflux /= 0.0_wp )  .AND.                        &
2078            surface_waterflux == 9999999.9_wp )  THEN
2079          message_string = 'wall_humidityflux additionally requires ' //     &
2080                           'setting of surface_waterflux'
2081          CALL message( 'check_parameters', 'PA0444', 1, 2, 0, 6, 0 )
2082       ENDIF
2083
2084       CALL set_bc_scalars( 'q', bc_q_b, bc_q_t, ibc_q_b, ibc_q_t,           &
2085                            'PA0071', 'PA0072' )
2086
2087       IF ( surface_waterflux == 9999999.9_wp  )  THEN
2088          constant_waterflux = .FALSE.
2089          IF ( large_scale_forcing .OR. land_surface )  THEN
2090             IF ( ibc_q_b == 0 )  THEN
2091                constant_waterflux = .FALSE.
2092             ELSEIF ( ibc_q_b == 1 )  THEN
2093                constant_waterflux = .TRUE.
2094             ENDIF
2095          ENDIF
2096       ELSE
2097          constant_waterflux = .TRUE.
2098       ENDIF
2099
2100       CALL check_bc_scalars( 'q', bc_q_b, ibc_q_b, 'PA0073', 'PA0074',        &
2101                              constant_waterflux, q_surface_initial_change )
2102
2103    ENDIF
2104
2105    IF ( passive_scalar )  THEN
2106
2107       IF ( ANY( wall_scalarflux /= 0.0_wp )  .AND.                            &
2108            surface_scalarflux == 9999999.9_wp )  THEN
2109          message_string = 'wall_scalarflux additionally requires ' //         &
2110                           'setting of surface_scalarflux'
2111          CALL message( 'check_parameters', 'PA0445', 1, 2, 0, 6, 0 )
2112       ENDIF
2113
2114       IF ( surface_scalarflux == 9999999.9_wp )  constant_scalarflux = .FALSE.
2115
2116       CALL set_bc_scalars( 's', bc_s_b, bc_s_t, ibc_s_b, ibc_s_t,             &
2117                            'PA0071', 'PA0072' )
2118
2119       CALL check_bc_scalars( 's', bc_s_b, ibc_s_b, 'PA0073', 'PA0074',        &
2120                              constant_scalarflux, s_surface_initial_change )
2121
2122       IF ( top_scalarflux == 9999999.9_wp )  constant_top_scalarflux = .FALSE.
2123!
2124!--    A fixed scalar concentration at the top implies Dirichlet boundary
2125!--    condition for scalar. Hence, in this case specification of a constant
2126!--    scalar flux is forbidden.
2127       IF ( ( ibc_s_t == 0 .OR. ibc_s_t == 2 )  .AND.  constant_top_scalarflux &
2128               .AND.  top_scalarflux /= 0.0_wp )  THEN
2129          message_string = 'boundary condition: bc_s_t = "' //                 &
2130                           TRIM( bc_s_t ) // '" is not allowed with ' //       &
2131                           'top_scalarflux /= 0.0'
2132          CALL message( 'check_parameters', 'PA0441', 1, 2, 0, 6, 0 )
2133       ENDIF
2134    ENDIF
2135
2136!
2137!-- Boundary conditions for chemical species
2138    IF ( air_chemistry )  CALL chem_boundary_conds( 'init' )
2139
2140!
2141!-- Boundary conditions for horizontal components of wind speed
2142    IF ( bc_uv_b == 'dirichlet' )  THEN
2143       ibc_uv_b = 0
2144    ELSEIF ( bc_uv_b == 'neumann' )  THEN
2145       ibc_uv_b = 1
2146       IF ( constant_flux_layer )  THEN
2147          message_string = 'boundary condition: bc_uv_b = "' //                &
2148               TRIM( bc_uv_b ) // '" is not allowed with constant_flux_layer'  &
2149               // ' = .TRUE.'
2150          CALL message( 'check_parameters', 'PA0075', 1, 2, 0, 6, 0 )
2151       ENDIF
2152    ELSE
2153       message_string = 'unknown boundary condition: bc_uv_b = "' //           &
2154                        TRIM( bc_uv_b ) // '"'
2155       CALL message( 'check_parameters', 'PA0076', 1, 2, 0, 6, 0 )
2156    ENDIF
2157!
2158!-- In case of coupled simulations u and v at the ground in atmosphere will be
2159!-- assigned with the u and v values of the ocean surface
2160    IF ( coupling_mode == 'atmosphere_to_ocean' )  THEN
2161       ibc_uv_b = 2
2162    ENDIF
2163
2164    IF ( coupling_mode == 'ocean_to_atmosphere' )  THEN
2165       bc_uv_t = 'neumann'
2166       ibc_uv_t = 1
2167    ELSE
2168       IF ( bc_uv_t == 'dirichlet' .OR. bc_uv_t == 'dirichlet_0' )  THEN
2169          ibc_uv_t = 0
2170          IF ( bc_uv_t == 'dirichlet_0' )  THEN
2171!
2172!--          Velocities for the initial u,v-profiles are set zero at the top
2173!--          in case of dirichlet_0 conditions
2174             u_init(nzt+1)    = 0.0_wp
2175             v_init(nzt+1)    = 0.0_wp
2176          ENDIF
2177       ELSEIF ( bc_uv_t == 'neumann' )  THEN
2178          ibc_uv_t = 1
2179       ELSEIF ( bc_uv_t == 'nested'  .OR.  bc_uv_t == 'nesting_offline' )  THEN
2180          ibc_uv_t = 3
2181       ELSE
2182          message_string = 'unknown boundary condition: bc_uv_t = "' //        &
2183                           TRIM( bc_uv_t ) // '"'
2184          CALL message( 'check_parameters', 'PA0077', 1, 2, 0, 6, 0 )
2185       ENDIF
2186    ENDIF
2187
2188!
2189!-- Compute and check, respectively, the Rayleigh Damping parameter
2190    IF ( rayleigh_damping_factor == -1.0_wp )  THEN
2191       rayleigh_damping_factor = 0.0_wp
2192    ELSE
2193       IF ( rayleigh_damping_factor < 0.0_wp  .OR.                             &
2194            rayleigh_damping_factor > 1.0_wp )  THEN
2195          WRITE( message_string, * )  'rayleigh_damping_factor = ',            &
2196                              rayleigh_damping_factor, ' out of range [0.0,1.0]'
2197          CALL message( 'check_parameters', 'PA0078', 1, 2, 0, 6, 0 )
2198       ENDIF
2199    ENDIF
2200
2201    IF ( rayleigh_damping_height == -1.0_wp )  THEN
2202       IF (  .NOT.  ocean_mode )  THEN
2203          rayleigh_damping_height = 0.66666666666_wp * zu(nzt)
2204       ELSE
2205          rayleigh_damping_height = 0.66666666666_wp * zu(nzb)
2206       ENDIF
2207    ELSE
2208       IF (  .NOT.  ocean_mode )  THEN
2209          IF ( rayleigh_damping_height < 0.0_wp  .OR.                          &
2210               rayleigh_damping_height > zu(nzt) )  THEN
2211             WRITE( message_string, * )  'rayleigh_damping_height = ',         &
2212                   rayleigh_damping_height, ' out of range [0.0,', zu(nzt), ']'
2213             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
2214          ENDIF
2215       ELSE
2216          IF ( rayleigh_damping_height > 0.0_wp  .OR.                          &
2217               rayleigh_damping_height < zu(nzb) )  THEN
2218             WRITE( message_string, * )  'rayleigh_damping_height = ',         &
2219                   rayleigh_damping_height, ' out of range [0.0,', zu(nzb), ']'
2220             CALL message( 'check_parameters', 'PA0079', 1, 2, 0, 6, 0 )
2221          ENDIF
2222       ENDIF
2223    ENDIF
2224
2225!
2226!-- Check number of chosen statistic regions
2227    IF ( statistic_regions < 0 )  THEN
2228       WRITE ( message_string, * ) 'number of statistic_regions = ',           &
2229                   statistic_regions+1, ' is not allowed'
2230       CALL message( 'check_parameters', 'PA0082', 1, 2, 0, 6, 0 )
2231    ENDIF
2232    IF ( normalizing_region > statistic_regions  .OR.                          &
2233         normalizing_region < 0)  THEN
2234       WRITE ( message_string, * ) 'normalizing_region = ',                    &
2235                normalizing_region, ' must be >= 0 and <= ',statistic_regions, &
2236                ' (value of statistic_regions)'
2237       CALL message( 'check_parameters', 'PA0083', 1, 2, 0, 6, 0 )
2238    ENDIF
2239
2240!
2241!-- Set the default intervals for data output, if necessary
2242!-- NOTE: dt_dosp has already been set in spectra_parin
2243    IF ( dt_data_output /= 9999999.9_wp )  THEN
2244       IF ( dt_dopr           == 9999999.9_wp )  dt_dopr           = dt_data_output
2245       IF ( dt_dopts          == 9999999.9_wp )  dt_dopts          = dt_data_output
2246       IF ( dt_do2d_xy        == 9999999.9_wp )  dt_do2d_xy        = dt_data_output
2247       IF ( dt_do2d_xz        == 9999999.9_wp )  dt_do2d_xz        = dt_data_output
2248       IF ( dt_do2d_yz        == 9999999.9_wp )  dt_do2d_yz        = dt_data_output
2249       IF ( dt_do3d           == 9999999.9_wp )  dt_do3d           = dt_data_output
2250       IF ( dt_data_output_av == 9999999.9_wp )  dt_data_output_av = dt_data_output
2251       DO  mid = 1, max_masks
2252          IF ( dt_domask(mid) == 9999999.9_wp )  dt_domask(mid)    = dt_data_output
2253       ENDDO
2254    ENDIF
2255
2256!
2257!-- Set the default skip time intervals for data output, if necessary
2258    IF ( skip_time_dopr    == 9999999.9_wp )                                   &
2259                                       skip_time_dopr    = skip_time_data_output
2260    IF ( skip_time_do2d_xy == 9999999.9_wp )                                   &
2261                                       skip_time_do2d_xy = skip_time_data_output
2262    IF ( skip_time_do2d_xz == 9999999.9_wp )                                   &
2263                                       skip_time_do2d_xz = skip_time_data_output
2264    IF ( skip_time_do2d_yz == 9999999.9_wp )                                   &
2265                                       skip_time_do2d_yz = skip_time_data_output
2266    IF ( skip_time_do3d    == 9999999.9_wp )                                   &
2267                                       skip_time_do3d    = skip_time_data_output
2268    IF ( skip_time_data_output_av == 9999999.9_wp )                            &
2269                                skip_time_data_output_av = skip_time_data_output
2270    DO  mid = 1, max_masks
2271       IF ( skip_time_domask(mid) == 9999999.9_wp )                            &
2272                                skip_time_domask(mid)    = skip_time_data_output
2273    ENDDO
2274
2275!
2276!-- Check the average intervals (first for 3d-data, then for profiles)
2277    IF ( averaging_interval > dt_data_output_av )  THEN
2278       WRITE( message_string, * )  'averaging_interval = ',                    &
2279             averaging_interval, ' must be <= dt_data_output_av = ',           &
2280             dt_data_output_av
2281       CALL message( 'check_parameters', 'PA0085', 1, 2, 0, 6, 0 )
2282    ENDIF
2283
2284    IF ( averaging_interval_pr == 9999999.9_wp )  THEN
2285       averaging_interval_pr = averaging_interval
2286    ENDIF
2287
2288    IF ( averaging_interval_pr > dt_dopr )  THEN
2289       WRITE( message_string, * )  'averaging_interval_pr = ',                 &
2290             averaging_interval_pr, ' must be <= dt_dopr = ', dt_dopr
2291       CALL message( 'check_parameters', 'PA0086', 1, 2, 0, 6, 0 )
2292    ENDIF
2293
2294!
2295!-- Set the default interval for profiles entering the temporal average
2296    IF ( dt_averaging_input_pr == 9999999.9_wp )  THEN
2297       dt_averaging_input_pr = dt_averaging_input
2298    ENDIF
2299
2300!
2301!-- Set the default interval for the output of timeseries to a reasonable
2302!-- value (tries to minimize the number of calls of flow_statistics)
2303    IF ( dt_dots == 9999999.9_wp )  THEN
2304       IF ( averaging_interval_pr == 0.0_wp )  THEN
2305          dt_dots = MIN( dt_run_control, dt_dopr )
2306       ELSE
2307          dt_dots = MIN( dt_run_control, dt_averaging_input_pr )
2308       ENDIF
2309    ENDIF
2310
2311!
2312!-- Check the sample rate for averaging (first for 3d-data, then for profiles)
2313    IF ( dt_averaging_input > averaging_interval )  THEN
2314       WRITE( message_string, * )  'dt_averaging_input = ',                    &
2315                dt_averaging_input, ' must be <= averaging_interval = ',       &
2316                averaging_interval
2317       CALL message( 'check_parameters', 'PA0088', 1, 2, 0, 6, 0 )
2318    ENDIF
2319
2320    IF ( dt_averaging_input_pr > averaging_interval_pr )  THEN
2321       WRITE( message_string, * )  'dt_averaging_input_pr = ',                 &
2322                dt_averaging_input_pr, ' must be <= averaging_interval_pr = ', &
2323                averaging_interval_pr
2324       CALL message( 'check_parameters', 'PA0089', 1, 2, 0, 6, 0 )
2325    ENDIF
2326
2327!
2328!-- Determine the number of output profiles and check whether they are
2329!-- permissible
2330    DO  WHILE ( data_output_pr(dopr_n+1) /= '          ' )
2331
2332       dopr_n = dopr_n + 1
2333       i = dopr_n
2334
2335!
2336!--    Determine internal profile number (for hom, homs)
2337!--    and store height levels
2338       SELECT CASE ( TRIM( data_output_pr(i) ) )
2339
2340          CASE ( 'u', '#u' )
2341             dopr_index(i) = 1
2342             dopr_unit(i)  = 'm/s'
2343             hom(:,2,1,:)  = SPREAD( zu, 2, statistic_regions+1 )
2344             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2345                dopr_initial_index(i) = 5
2346                hom(:,2,5,:)          = SPREAD( zu, 2, statistic_regions+1 )
2347                data_output_pr(i)     = data_output_pr(i)(2:)
2348             ENDIF
2349
2350          CASE ( 'v', '#v' )
2351             dopr_index(i) = 2
2352             dopr_unit(i)  = 'm/s'
2353             hom(:,2,2,:)  = SPREAD( zu, 2, statistic_regions+1 )
2354             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2355                dopr_initial_index(i) = 6
2356                hom(:,2,6,:)          = SPREAD( zu, 2, statistic_regions+1 )
2357                data_output_pr(i)     = data_output_pr(i)(2:)
2358             ENDIF
2359
2360          CASE ( 'w' )
2361             dopr_index(i) = 3
2362             dopr_unit(i)  = 'm/s'
2363             hom(:,2,3,:)  = SPREAD( zw, 2, statistic_regions+1 )
2364
2365          CASE ( 'pt', '#pt' )
2366             IF ( .NOT. bulk_cloud_model ) THEN
2367                dopr_index(i) = 4
2368                dopr_unit(i)  = 'K'
2369                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2370                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2371                   dopr_initial_index(i) = 7
2372                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2373                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
2374                   data_output_pr(i)     = data_output_pr(i)(2:)
2375                ENDIF
2376             ELSE
2377                dopr_index(i) = 43
2378                dopr_unit(i)  = 'K'
2379                hom(:,2,43,:)  = SPREAD( zu, 2, statistic_regions+1 )
2380                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2381                   dopr_initial_index(i) = 28
2382                   hom(:,2,28,:)         = SPREAD( zu, 2, statistic_regions+1 )
2383                   hom(nzb,2,28,:)       = 0.0_wp    ! because zu(nzb) is negative
2384                   data_output_pr(i)     = data_output_pr(i)(2:)
2385                ENDIF
2386             ENDIF
2387
2388          CASE ( 'e', '#e' )
2389             dopr_index(i)  = 8
2390             dopr_unit(i)   = 'm2/s2'
2391             hom(:,2,8,:)   = SPREAD( zu, 2, statistic_regions+1 )
2392             hom(nzb,2,8,:) = 0.0_wp
2393             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2394                dopr_initial_index(i) = 8
2395                hom(:,2,8,:)          = SPREAD( zu, 2, statistic_regions+1 )
2396                data_output_pr(i)     = data_output_pr(i)(2:)
2397             ENDIF
2398
2399          CASE ( 'km', '#km' )
2400             dopr_index(i)  = 9
2401             dopr_unit(i)   = 'm2/s'
2402             hom(:,2,9,:)   = SPREAD( zu, 2, statistic_regions+1 )
2403             hom(nzb,2,9,:) = 0.0_wp
2404             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2405                dopr_initial_index(i) = 23
2406                hom(:,2,23,:)         = hom(:,2,9,:)
2407                data_output_pr(i)     = data_output_pr(i)(2:)
2408             ENDIF
2409
2410          CASE ( 'kh', '#kh' )
2411             dopr_index(i)   = 10
2412             dopr_unit(i)    = 'm2/s'
2413             hom(:,2,10,:)   = SPREAD( zu, 2, statistic_regions+1 )
2414             hom(nzb,2,10,:) = 0.0_wp
2415             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2416                dopr_initial_index(i) = 24
2417                hom(:,2,24,:)         = hom(:,2,10,:)
2418                data_output_pr(i)     = data_output_pr(i)(2:)
2419             ENDIF
2420
2421          CASE ( 'l', '#l' )
2422             dopr_index(i)   = 11
2423             dopr_unit(i)    = 'm'
2424             hom(:,2,11,:)   = SPREAD( zu, 2, statistic_regions+1 )
2425             hom(nzb,2,11,:) = 0.0_wp
2426             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2427                dopr_initial_index(i) = 25
2428                hom(:,2,25,:)         = hom(:,2,11,:)
2429                data_output_pr(i)     = data_output_pr(i)(2:)
2430             ENDIF
2431
2432          CASE ( 'w"u"' )
2433             dopr_index(i) = 12
2434             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2435             hom(:,2,12,:) = SPREAD( zw, 2, statistic_regions+1 )
2436             IF ( constant_flux_layer )  hom(nzb,2,12,:) = zu(1)
2437
2438          CASE ( 'w*u*' )
2439             dopr_index(i) = 13
2440             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2441             hom(:,2,13,:) = SPREAD( zw, 2, statistic_regions+1 )
2442
2443          CASE ( 'w"v"' )
2444             dopr_index(i) = 14
2445             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2446             hom(:,2,14,:) = SPREAD( zw, 2, statistic_regions+1 )
2447             IF ( constant_flux_layer )  hom(nzb,2,14,:) = zu(1)
2448
2449          CASE ( 'w*v*' )
2450             dopr_index(i) = 15
2451             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2452             hom(:,2,15,:) = SPREAD( zw, 2, statistic_regions+1 )
2453
2454          CASE ( 'w"pt"' )
2455             dopr_index(i) = 16
2456             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2457             hom(:,2,16,:) = SPREAD( zw, 2, statistic_regions+1 )
2458
2459          CASE ( 'w*pt*' )
2460             dopr_index(i) = 17
2461             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2462             hom(:,2,17,:) = SPREAD( zw, 2, statistic_regions+1 )
2463
2464          CASE ( 'wpt' )
2465             dopr_index(i) = 18
2466             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2467             hom(:,2,18,:) = SPREAD( zw, 2, statistic_regions+1 )
2468
2469          CASE ( 'wu' )
2470             dopr_index(i) = 19
2471             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2472             hom(:,2,19,:) = SPREAD( zw, 2, statistic_regions+1 )
2473             IF ( constant_flux_layer )  hom(nzb,2,19,:) = zu(1)
2474
2475          CASE ( 'wv' )
2476             dopr_index(i) = 20
2477             dopr_unit(i)  = TRIM ( momentumflux_output_unit )
2478             hom(:,2,20,:) = SPREAD( zw, 2, statistic_regions+1 )
2479             IF ( constant_flux_layer )  hom(nzb,2,20,:) = zu(1)
2480
2481          CASE ( 'w*pt*BC' )
2482             dopr_index(i) = 21
2483             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2484             hom(:,2,21,:) = SPREAD( zw, 2, statistic_regions+1 )
2485
2486          CASE ( 'wptBC' )
2487             dopr_index(i) = 22
2488             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2489             hom(:,2,22,:) = SPREAD( zw, 2, statistic_regions+1 )
2490
2491          CASE ( 'u*2' )
2492             dopr_index(i) = 30
2493             dopr_unit(i)  = 'm2/s2'
2494             hom(:,2,30,:) = SPREAD( zu, 2, statistic_regions+1 )
2495
2496          CASE ( 'v*2' )
2497             dopr_index(i) = 31
2498             dopr_unit(i)  = 'm2/s2'
2499             hom(:,2,31,:) = SPREAD( zu, 2, statistic_regions+1 )
2500
2501          CASE ( 'w*2' )
2502             dopr_index(i) = 32
2503             dopr_unit(i)  = 'm2/s2'
2504             hom(:,2,32,:) = SPREAD( zw, 2, statistic_regions+1 )
2505
2506          CASE ( 'pt*2' )
2507             dopr_index(i) = 33
2508             dopr_unit(i)  = 'K2'
2509             hom(:,2,33,:) = SPREAD( zu, 2, statistic_regions+1 )
2510
2511          CASE ( 'e*' )
2512             dopr_index(i) = 34
2513             dopr_unit(i)  = 'm2/s2'
2514             hom(:,2,34,:) = SPREAD( zu, 2, statistic_regions+1 )
2515
2516          CASE ( 'w*2pt*' )
2517             dopr_index(i) = 35
2518             dopr_unit(i)  = 'K m2/s2'
2519             hom(:,2,35,:) = SPREAD( zw, 2, statistic_regions+1 )
2520
2521          CASE ( 'w*pt*2' )
2522             dopr_index(i) = 36
2523             dopr_unit(i)  = 'K2 m/s'
2524             hom(:,2,36,:) = SPREAD( zw, 2, statistic_regions+1 )
2525
2526          CASE ( 'w*e*' )
2527             dopr_index(i) = 37
2528             dopr_unit(i)  = 'm3/s3'
2529             hom(:,2,37,:) = SPREAD( zw, 2, statistic_regions+1 )
2530
2531          CASE ( 'w*3' )
2532             dopr_index(i) = 38
2533             dopr_unit(i)  = 'm3/s3'
2534             hom(:,2,38,:) = SPREAD( zw, 2, statistic_regions+1 )
2535
2536          CASE ( 'Sw' )
2537             dopr_index(i) = 39
2538             dopr_unit(i)  = 'none'
2539             hom(:,2,39,:) = SPREAD( zw, 2, statistic_regions+1 )
2540
2541          CASE ( 'p' )
2542             dopr_index(i) = 40
2543             dopr_unit(i)  = 'Pa'
2544             hom(:,2,40,:) = SPREAD( zu, 2, statistic_regions+1 )
2545
2546          CASE ( 'q', '#q' )
2547             IF ( .NOT. humidity )  THEN
2548                message_string = 'data_output_pr = ' //                        &
2549                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2550                                 'lemented for humidity = .FALSE.'
2551                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2552             ELSE
2553                dopr_index(i) = 41
2554                dopr_unit(i)  = 'kg/kg'
2555                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2556                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2557                   dopr_initial_index(i) = 26
2558                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2559                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2560                   data_output_pr(i)     = data_output_pr(i)(2:)
2561                ENDIF
2562             ENDIF
2563
2564          CASE ( 's', '#s' )
2565             IF ( .NOT. passive_scalar )  THEN
2566                message_string = 'data_output_pr = ' //                        &
2567                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2568                                 'lemented for passive_scalar = .FALSE.'
2569                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2570             ELSE
2571                dopr_index(i) = 115
2572                dopr_unit(i)  = 'kg/m3'
2573                hom(:,2,115,:) = SPREAD( zu, 2, statistic_regions+1 )
2574                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2575                   dopr_initial_index(i) = 121
2576                   hom(:,2,121,:)        = SPREAD( zu, 2, statistic_regions+1 )
2577                   hom(nzb,2,121,:)      = 0.0_wp    ! because zu(nzb) is negative
2578                   data_output_pr(i)     = data_output_pr(i)(2:)
2579                ENDIF
2580             ENDIF
2581
2582          CASE ( 'qv', '#qv' )
2583             IF ( .NOT. bulk_cloud_model ) THEN
2584                dopr_index(i) = 41
2585                dopr_unit(i)  = 'kg/kg'
2586                hom(:,2,41,:) = SPREAD( zu, 2, statistic_regions+1 )
2587                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2588                   dopr_initial_index(i) = 26
2589                   hom(:,2,26,:)         = SPREAD( zu, 2, statistic_regions+1 )
2590                   hom(nzb,2,26,:)       = 0.0_wp    ! because zu(nzb) is negative
2591                   data_output_pr(i)     = data_output_pr(i)(2:)
2592                ENDIF
2593             ELSE
2594                dopr_index(i) = 42
2595                dopr_unit(i)  = 'kg/kg'
2596                hom(:,2,42,:) = SPREAD( zu, 2, statistic_regions+1 )
2597                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2598                   dopr_initial_index(i) = 27
2599                   hom(:,2,27,:)         = SPREAD( zu, 2, statistic_regions+1 )
2600                   hom(nzb,2,27,:)       = 0.0_wp   ! because zu(nzb) is negative
2601                   data_output_pr(i)     = data_output_pr(i)(2:)
2602                ENDIF
2603             ENDIF
2604
2605          CASE ( 'lpt', '#lpt' )
2606             IF ( .NOT. bulk_cloud_model ) THEN
2607                message_string = 'data_output_pr = ' //                        &
2608                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2609                                 'lemented for bulk_cloud_model = .FALSE.'
2610                CALL message( 'check_parameters', 'PA0094', 1, 2, 0, 6, 0 )
2611             ELSE
2612                dopr_index(i) = 4
2613                dopr_unit(i)  = 'K'
2614                hom(:,2,4,:)  = SPREAD( zu, 2, statistic_regions+1 )
2615                IF ( data_output_pr(i)(1:1) == '#' )  THEN
2616                   dopr_initial_index(i) = 7
2617                   hom(:,2,7,:)          = SPREAD( zu, 2, statistic_regions+1 )
2618                   hom(nzb,2,7,:)        = 0.0_wp    ! because zu(nzb) is negative
2619                   data_output_pr(i)     = data_output_pr(i)(2:)
2620                ENDIF
2621             ENDIF
2622
2623          CASE ( 'vpt', '#vpt' )
2624             dopr_index(i) = 44
2625             dopr_unit(i)  = 'K'
2626             hom(:,2,44,:) = SPREAD( zu, 2, statistic_regions+1 )
2627             IF ( data_output_pr(i)(1:1) == '#' )  THEN
2628                dopr_initial_index(i) = 29
2629                hom(:,2,29,:)         = SPREAD( zu, 2, statistic_regions+1 )
2630                hom(nzb,2,29,:)       = 0.0_wp    ! because zu(nzb) is negative
2631                data_output_pr(i)     = data_output_pr(i)(2:)
2632             ENDIF
2633
2634          CASE ( 'w"vpt"' )
2635             dopr_index(i) = 45
2636             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2637             hom(:,2,45,:) = SPREAD( zw, 2, statistic_regions+1 )
2638
2639          CASE ( 'w*vpt*' )
2640             dopr_index(i) = 46
2641             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2642             hom(:,2,46,:) = SPREAD( zw, 2, statistic_regions+1 )
2643
2644          CASE ( 'wvpt' )
2645             dopr_index(i) = 47
2646             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2647             hom(:,2,47,:) = SPREAD( zw, 2, statistic_regions+1 )
2648
2649          CASE ( 'w"q"' )
2650             IF (  .NOT.  humidity )  THEN
2651                message_string = 'data_output_pr = ' //                        &
2652                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2653                                 'lemented for humidity = .FALSE.'
2654                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2655             ELSE
2656                dopr_index(i) = 48
2657                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2658                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2659             ENDIF
2660
2661          CASE ( 'w*q*' )
2662             IF (  .NOT.  humidity )  THEN
2663                message_string = 'data_output_pr = ' //                        &
2664                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2665                                 'lemented for humidity = .FALSE.'
2666                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2667             ELSE
2668                dopr_index(i) = 49
2669                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2670                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2671             ENDIF
2672
2673          CASE ( 'wq' )
2674             IF (  .NOT.  humidity )  THEN
2675                message_string = 'data_output_pr = ' //                        &
2676                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2677                                 'lemented for humidity = .FALSE.'
2678                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2679             ELSE
2680                dopr_index(i) = 50
2681                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2682                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2683             ENDIF
2684
2685          CASE ( 'w"s"' )
2686             IF (  .NOT.  passive_scalar )  THEN
2687                message_string = 'data_output_pr = ' //                        &
2688                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2689                                 'lemented for passive_scalar = .FALSE.'
2690                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2691             ELSE
2692                dopr_index(i) = 117
2693                dopr_unit(i)  = 'kg/m3 m/s'
2694                hom(:,2,117,:) = SPREAD( zw, 2, statistic_regions+1 )
2695             ENDIF
2696
2697          CASE ( 'w*s*' )
2698             IF (  .NOT.  passive_scalar )  THEN
2699                message_string = 'data_output_pr = ' //                        &
2700                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2701                                 'lemented for passive_scalar = .FALSE.'
2702                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2703             ELSE
2704                dopr_index(i) = 114
2705                dopr_unit(i)  = 'kg/m3 m/s'
2706                hom(:,2,114,:) = SPREAD( zw, 2, statistic_regions+1 )
2707             ENDIF
2708
2709          CASE ( 'ws' )
2710             IF (  .NOT.  passive_scalar )  THEN
2711                message_string = 'data_output_pr = ' //                        &
2712                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2713                                 'lemented for passive_scalar = .FALSE.'
2714                CALL message( 'check_parameters', 'PA0093', 1, 2, 0, 6, 0 )
2715             ELSE
2716                dopr_index(i) = 118
2717                dopr_unit(i)  = 'kg/m3 m/s'
2718                hom(:,2,118,:) = SPREAD( zw, 2, statistic_regions+1 )
2719             ENDIF
2720
2721          CASE ( 'w"qv"' )
2722             IF ( humidity  .AND.  .NOT.  bulk_cloud_model )  THEN
2723                dopr_index(i) = 48
2724                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2725                hom(:,2,48,:) = SPREAD( zw, 2, statistic_regions+1 )
2726             ELSEIF ( humidity  .AND.  bulk_cloud_model )  THEN
2727                dopr_index(i) = 51
2728                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2729                hom(:,2,51,:) = SPREAD( zw, 2, statistic_regions+1 )
2730             ELSE
2731                message_string = 'data_output_pr = ' //                        &
2732                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2733                                 'lemented for bulk_cloud_model = .FALSE. ' // &
2734                                 'and humidity = .FALSE.'
2735                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2736             ENDIF
2737
2738          CASE ( 'w*qv*' )
2739             IF ( humidity  .AND.  .NOT. bulk_cloud_model )  THEN
2740                dopr_index(i) = 49
2741                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2742                hom(:,2,49,:) = SPREAD( zw, 2, statistic_regions+1 )
2743             ELSEIF( humidity .AND. bulk_cloud_model ) THEN
2744                dopr_index(i) = 52
2745                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2746                hom(:,2,52,:) = SPREAD( zw, 2, statistic_regions+1 )
2747             ELSE
2748                message_string = 'data_output_pr = ' //                        &
2749                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2750                                 'lemented for bulk_cloud_model = .FALSE. ' // &
2751                                 'and humidity = .FALSE.'
2752                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2753             ENDIF
2754
2755          CASE ( 'wqv' )
2756             IF ( humidity  .AND.  .NOT.  bulk_cloud_model )  THEN
2757                dopr_index(i) = 50
2758                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2759                hom(:,2,50,:) = SPREAD( zw, 2, statistic_regions+1 )
2760             ELSEIF ( humidity  .AND.  bulk_cloud_model )  THEN
2761                dopr_index(i) = 53
2762                dopr_unit(i)  = TRIM ( waterflux_output_unit )
2763                hom(:,2,53,:) = SPREAD( zw, 2, statistic_regions+1 )
2764             ELSE
2765                message_string = 'data_output_pr = ' //                        &
2766                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2767                                 'lemented for bulk_cloud_model = .FALSE. ' // &
2768                                 'and humidity = .FALSE.'
2769                CALL message( 'check_parameters', 'PA0095', 1, 2, 0, 6, 0 )
2770             ENDIF
2771
2772          CASE ( 'ql' )
2773             IF (  .NOT.  bulk_cloud_model  .AND.  .NOT.  cloud_droplets )  THEN
2774                message_string = 'data_output_pr = ' //                        &
2775                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2776                                 'lemented for bulk_cloud_model = .FALSE. ' // &
2777                                 'and cloud_droplets = .FALSE.'
2778                CALL message( 'check_parameters', 'PA0096', 1, 2, 0, 6, 0 )
2779             ELSE
2780                dopr_index(i) = 54
2781                dopr_unit(i)  = 'kg/kg'
2782                hom(:,2,54,:)  = SPREAD( zu, 2, statistic_regions+1 )
2783             ENDIF
2784
2785          CASE ( 'w*u*u*:dz' )
2786             dopr_index(i) = 55
2787             dopr_unit(i)  = 'm2/s3'
2788             hom(:,2,55,:) = SPREAD( zu, 2, statistic_regions+1 )
2789
2790          CASE ( 'w*p*:dz' )
2791             dopr_index(i) = 56
2792             dopr_unit(i)  = 'm2/s3'
2793             hom(:,2,56,:) = SPREAD( zw, 2, statistic_regions+1 )
2794
2795          CASE ( 'w"e:dz' )
2796             dopr_index(i) = 57
2797             dopr_unit(i)  = 'm2/s3'
2798             hom(:,2,57,:) = SPREAD( zu, 2, statistic_regions+1 )
2799
2800          CASE ( 'u"pt"' )
2801             dopr_index(i) = 58
2802             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2803             hom(:,2,58,:) = SPREAD( zu, 2, statistic_regions+1 )
2804
2805          CASE ( 'u*pt*' )
2806             dopr_index(i) = 59
2807             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2808             hom(:,2,59,:) = SPREAD( zu, 2, statistic_regions+1 )
2809
2810          CASE ( 'upt_t' )
2811             dopr_index(i) = 60
2812             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2813             hom(:,2,60,:) = SPREAD( zu, 2, statistic_regions+1 )
2814
2815          CASE ( 'v"pt"' )
2816             dopr_index(i) = 61
2817             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2818             hom(:,2,61,:) = SPREAD( zu, 2, statistic_regions+1 )
2819
2820          CASE ( 'v*pt*' )
2821             dopr_index(i) = 62
2822             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2823             hom(:,2,62,:) = SPREAD( zu, 2, statistic_regions+1 )
2824
2825          CASE ( 'vpt_t' )
2826             dopr_index(i) = 63
2827             dopr_unit(i)  = TRIM ( heatflux_output_unit )
2828             hom(:,2,63,:) = SPREAD( zu, 2, statistic_regions+1 )
2829
2830          CASE ( 'w*p*' )
2831             dopr_index(i) = 68
2832             dopr_unit(i)  = 'm3/s3'
2833             hom(:,2,68,:) = SPREAD( zu, 2, statistic_regions+1 )
2834
2835          CASE ( 'w"e' )
2836             dopr_index(i) = 69
2837             dopr_unit(i)  = 'm3/s3'
2838             hom(:,2,69,:) = SPREAD( zu, 2, statistic_regions+1 )
2839
2840          CASE ( 'q*2' )
2841             IF (  .NOT.  humidity )  THEN
2842                message_string = 'data_output_pr = ' //                        &
2843                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2844                                 'lemented for humidity = .FALSE.'
2845                CALL message( 'check_parameters', 'PA0092', 1, 2, 0, 6, 0 )
2846             ELSE
2847                dopr_index(i) = 70
2848                dopr_unit(i)  = 'kg2/kg2'
2849                hom(:,2,70,:) = SPREAD( zu, 2, statistic_regions+1 )
2850             ENDIF
2851
2852          CASE ( 'hyp' )
2853             dopr_index(i) = 72
2854             dopr_unit(i)  = 'hPa'
2855             hom(:,2,72,:) = SPREAD( zu, 2, statistic_regions+1 )
2856
2857          CASE ( 'rho_air' )
2858             dopr_index(i)  = 119
2859             dopr_unit(i)   = 'kg/m3'
2860             hom(:,2,119,:) = SPREAD( zu, 2, statistic_regions+1 )
2861
2862          CASE ( 'rho_air_zw' )
2863             dopr_index(i)  = 120
2864             dopr_unit(i)   = 'kg/m3'
2865             hom(:,2,120,:) = SPREAD( zw, 2, statistic_regions+1 )
2866
2867          CASE ( 'ug' )
2868             dopr_index(i) = 78
2869             dopr_unit(i)  = 'm/s'
2870             hom(:,2,78,:) = SPREAD( zu, 2, statistic_regions+1 )
2871
2872          CASE ( 'vg' )
2873             dopr_index(i) = 79
2874             dopr_unit(i)  = 'm/s'
2875             hom(:,2,79,:) = SPREAD( zu, 2, statistic_regions+1 )
2876
2877          CASE ( 'w_subs' )
2878             IF (  .NOT.  large_scale_subsidence )  THEN
2879                message_string = 'data_output_pr = ' //                        &
2880                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2881                                 'lemented for large_scale_subsidence = .FALSE.'
2882                CALL message( 'check_parameters', 'PA0382', 1, 2, 0, 6, 0 )
2883             ELSE
2884                dopr_index(i) = 80
2885                dopr_unit(i)  = 'm/s'
2886                hom(:,2,80,:) = SPREAD( zu, 2, statistic_regions+1 )
2887             ENDIF
2888
2889          CASE ( 's*2' )
2890             IF (  .NOT.  passive_scalar )  THEN
2891                message_string = 'data_output_pr = ' //                        &
2892                                 TRIM( data_output_pr(i) ) // ' is not imp' // &
2893                                 'lemented for passive_scalar = .FALSE.'
2894                CALL message( 'check_parameters', 'PA0185', 1, 2, 0, 6, 0 )
2895             ELSE
2896                dopr_index(i) = 116
2897                dopr_unit(i)  = 'kg2/m6'
2898                hom(:,2,116,:) = SPREAD( zu, 2, statistic_regions+1 )
2899             ENDIF
2900
2901          CASE DEFAULT
2902
2903!
2904!--          Check for other modules
2905             IF ( unit == 'illegal'  .AND.  bulk_cloud_model  )  THEN
2906                CALL bcm_check_data_output_pr( data_output_pr(i), i, unit,     &
2907                                               dopr_unit(i) )
2908             ENDIF
2909
2910             IF ( unit == 'illegal' )  THEN
2911                CALL chem_check_data_output_pr( data_output_pr(i), i, unit,    &
2912                                                dopr_unit(i) )
2913             ENDIF
2914
2915             IF ( unit == 'illegal'  .AND.  gust_module_enabled  )  THEN
2916                CALL gust_check_data_output_pr( data_output_pr(i), i, unit,    &
2917                                                dopr_unit(i) )
2918             ENDIF
2919
2920             IF ( unit == 'illegal' )  THEN
2921                CALL lsf_nudging_check_data_output_pr( data_output_pr(i), i,   &
2922                                                       unit, dopr_unit(i) )
2923             ENDIF
2924
2925             CALL lsm_check_data_output_pr( data_output_pr(i), i, unit,        &
2926                                            dopr_unit(i) )
2927
2928             IF ( unit == 'illegal' )  THEN
2929                CALL ocean_check_data_output_pr( data_output_pr(i), i, unit,   &
2930                                                 dopr_unit(i) )
2931             ENDIF
2932
2933             IF ( unit == 'illegal' )  THEN
2934                CALL radiation_check_data_output_pr( data_output_pr(i), i,     &
2935                                                     unit, dopr_unit(i) )
2936             ENDIF
2937
2938!
2939!--          Finally, check for user defined quantities
2940             IF ( unit == 'illegal' )  THEN
2941                unit = ''
2942                CALL user_check_data_output_pr( data_output_pr(i), i, unit )
2943             ENDIF
2944
2945!
2946!--          No valid quantity found
2947             IF ( unit == 'illegal' )  THEN
2948                IF ( data_output_pr_user(1) /= ' ' )  THEN
2949                   message_string = 'illegal value for data_output_pr or ' //  &
2950                                    'data_output_pr_user = "' //               &
2951                                    TRIM( data_output_pr(i) ) // '"'
2952                   CALL message( 'check_parameters', 'PA0097', 1, 2, 0, 6, 0 )
2953                ELSE
2954                   message_string = 'illegal value for data_output_pr = "' //  &
2955                                    TRIM( data_output_pr(i) ) // '"'
2956                   CALL message( 'check_parameters', 'PA0098', 1, 2, 0, 6, 0 )
2957                ENDIF
2958             ENDIF
2959
2960       END SELECT
2961
2962    ENDDO
2963
2964
2965!
2966!-- Append user-defined data output variables to the standard data output
2967    IF ( data_output_user(1) /= ' ' )  THEN
2968       i = 1
2969       DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 500 )
2970          i = i + 1
2971       ENDDO
2972       j = 1
2973       DO  WHILE ( data_output_user(j) /= ' '  .AND.  j <= 500 )
2974          IF ( i > 500 )  THEN
2975             message_string = 'number of output quantitities given by data' // &
2976                '_output and data_output_user exceeds the limit of 500'
2977             CALL message( 'check_parameters', 'PA0102', 1, 2, 0, 6, 0 )
2978          ENDIF
2979          data_output(i) = data_output_user(j)
2980          i = i + 1
2981          j = j + 1
2982       ENDDO
2983    ENDIF
2984
2985!
2986!-- Check and set steering parameters for 2d/3d data output and averaging
2987    i   = 1
2988    DO  WHILE ( data_output(i) /= ' '  .AND.  i <= 500 )
2989!
2990!--    Check for data averaging
2991       ilen = LEN_TRIM( data_output(i) )
2992       j = 0                                                 ! no data averaging
2993       IF ( ilen > 3 )  THEN
2994          IF ( data_output(i)(ilen-2:ilen) == '_av' )  THEN
2995             j = 1                                           ! data averaging
2996             data_output(i) = data_output(i)(1:ilen-3)
2997          ENDIF
2998       ENDIF
2999!
3000!--    Check for cross section or volume data
3001       ilen = LEN_TRIM( data_output(i) )
3002       k = 0                                                   ! 3d data
3003       var = data_output(i)(1:ilen)
3004       IF ( ilen > 3 )  THEN
3005          IF ( data_output(i)(ilen-2:ilen) == '_xy'  .OR.                      &
3006               data_output(i)(ilen-2:ilen) == '_xz'  .OR.                      &
3007               data_output(i)(ilen-2:ilen) == '_yz' )  THEN
3008             k = 1                                             ! 2d data
3009             var = data_output(i)(1:ilen-3)
3010          ENDIF
3011       ENDIF
3012
3013!
3014!--    Check for allowed value and set units
3015       SELECT CASE ( TRIM( var ) )
3016
3017          CASE ( 'e' )
3018             IF ( constant_diffusion )  THEN
3019                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3020                                 'res constant_diffusion = .FALSE.'
3021                CALL message( 'check_parameters', 'PA0103', 1, 2, 0, 6, 0 )
3022             ENDIF
3023             unit = 'm2/s2'
3024
3025          CASE ( 'lpt' )
3026             IF (  .NOT.  bulk_cloud_model )  THEN
3027                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3028                         'res bulk_cloud_model = .TRUE.'
3029                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3030             ENDIF
3031             unit = 'K'
3032
3033          CASE ( 'pc', 'pr' )
3034             IF (  .NOT.  particle_advection )  THEN
3035                message_string = 'output of "' // TRIM( var ) // '" requir' // &
3036                   'es a "particle_parameters"-NAMELIST in the parameter ' //  &
3037                   'file (PARIN)'
3038                CALL message( 'check_parameters', 'PA0104', 1, 2, 0, 6, 0 )
3039             ENDIF
3040             IF ( TRIM( var ) == 'pc' )  unit = 'number'
3041             IF ( TRIM( var ) == 'pr' )  unit = 'm'
3042
3043          CASE ( 'q', 'vpt' )
3044             IF (  .NOT.  humidity )  THEN
3045                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3046                                 'res humidity = .TRUE.'
3047                CALL message( 'check_parameters', 'PA0105', 1, 2, 0, 6, 0 )
3048             ENDIF
3049             IF ( TRIM( var ) == 'q'   )  unit = 'kg/kg'
3050             IF ( TRIM( var ) == 'vpt' )  unit = 'K'
3051
3052          CASE ( 'ql' )
3053             IF ( .NOT.  ( bulk_cloud_model  .OR.  cloud_droplets ) )  THEN
3054                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3055                      'res bulk_cloud_model = .TRUE. or cloud_droplets = .TRUE.'
3056                CALL message( 'check_parameters', 'PA0106', 1, 2, 0, 6, 0 )
3057             ENDIF
3058             unit = 'kg/kg'
3059
3060          CASE ( 'ql_c', 'ql_v', 'ql_vp' )
3061             IF (  .NOT.  cloud_droplets )  THEN
3062                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3063                                 'res cloud_droplets = .TRUE.'
3064                CALL message( 'check_parameters', 'PA0107', 1, 2, 0, 6, 0 )
3065             ENDIF
3066             IF ( TRIM( var ) == 'ql_c'  )  unit = 'kg/kg'
3067             IF ( TRIM( var ) == 'ql_v'  )  unit = 'm3'
3068             IF ( TRIM( var ) == 'ql_vp' )  unit = 'none'
3069
3070          CASE ( 'qv' )
3071             IF (  .NOT.  bulk_cloud_model )  THEN
3072                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3073                                 'res bulk_cloud_model = .TRUE.'
3074                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3075             ENDIF
3076             unit = 'kg/kg'
3077
3078          CASE ( 's' )
3079             IF (  .NOT.  passive_scalar )  THEN
3080                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3081                                 'res passive_scalar = .TRUE.'
3082                CALL message( 'check_parameters', 'PA0110', 1, 2, 0, 6, 0 )
3083             ENDIF
3084             unit = 'kg/m3'
3085
3086          CASE ( 'p', 'pt', 'u', 'v', 'w' )
3087             IF ( TRIM( var ) == 'p'  )  unit = 'Pa'
3088             IF ( TRIM( var ) == 'pt' )  unit = 'K'
3089             IF ( TRIM( var ) == 'u'  )  unit = 'm/s'
3090             IF ( TRIM( var ) == 'v'  )  unit = 'm/s'
3091             IF ( TRIM( var ) == 'w'  )  unit = 'm/s'
3092             CONTINUE
3093
3094          CASE ( 'ghf*', 'lwp*', 'ol*', 'qsws*', 'r_a*',                       &
3095                 'shf*', 'ssws*', 't*', 'tsurf*', 'u*', 'z0*', 'z0h*', 'z0q*' )
3096             IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
3097                message_string = 'illegal value for data_output: "' //         &
3098                                 TRIM( var ) // '" & only 2d-horizontal ' //   &
3099                                 'cross sections are allowed for this value'
3100                CALL message( 'check_parameters', 'PA0111', 1, 2, 0, 6, 0 )
3101             ENDIF
3102
3103             IF ( TRIM( var ) == 'lwp*'  .AND.  .NOT. bulk_cloud_model )  THEN
3104                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3105                                 'res bulk_cloud_model = .TRUE.'
3106                CALL message( 'check_parameters', 'PA0108', 1, 2, 0, 6, 0 )
3107             ENDIF
3108             IF ( TRIM( var ) == 'qsws*'  .AND.  .NOT.  humidity )  THEN
3109                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3110                                 'res humidity = .TRUE.'
3111                CALL message( 'check_parameters', 'PA0322', 1, 2, 0, 6, 0 )
3112             ENDIF
3113
3114             IF ( TRIM( var ) == 'ghf*'  .AND.  .NOT.  land_surface )  THEN
3115                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3116                                 'res land_surface = .TRUE.'
3117                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3118             ENDIF
3119
3120             IF ( ( TRIM( var ) == 'r_a*' .OR.  TRIM( var ) == 'ghf*' )        &
3121                 .AND.  .NOT.  land_surface  .AND.  .NOT.  urban_surface )     &         
3122             THEN
3123                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3124                                 'res land_surface = .TRUE. or ' //            &
3125                                 'urban_surface = .TRUE.'
3126                CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
3127             ENDIF
3128             IF ( TRIM( var ) == 'ssws*'  .AND.  .NOT.  passive_scalar )  THEN
3129                message_string = 'output of "' // TRIM( var ) // '" requi' //  &
3130                                 'res passive_scalar = .TRUE.'
3131                CALL message( 'check_parameters', 'PA0361', 1, 2, 0, 6, 0 )
3132             ENDIF
3133
3134             IF ( TRIM( var ) == 'ghf*'   )  unit = 'W/m2'
3135             IF ( TRIM( var ) == 'lwp*'   )  unit = 'kg/m2'
3136             IF ( TRIM( var ) == 'ol*'    )  unit = 'm'
3137             IF ( TRIM( var ) == 'qsws*'  )  unit = 'kgm/kgs'
3138             IF ( TRIM( var ) == 'r_a*'   )  unit = 's/m'     
3139             IF ( TRIM( var ) == 'shf*'   )  unit = 'K*m/s'
3140             IF ( TRIM( var ) == 'ssws*'  )  unit = 'kg/m2*s'
3141             IF ( TRIM( var ) == 't*'     )  unit = 'K'
3142             IF ( TRIM( var ) == 'tsurf*' )  unit = 'K' 
3143             IF ( TRIM( var ) == 'u*'     )  unit = 'm/s'
3144             IF ( TRIM( var ) == 'z0*'    )  unit = 'm'
3145             IF ( TRIM( var ) == 'z0h*'   )  unit = 'm'
3146!
3147!--          Output of surface latent and sensible heat flux will be in W/m2
3148!--          in case of natural- and urban-type surfaces, even if
3149!--          flux_output_mode is set to kinematic units.
3150             IF ( land_surface  .OR.  urban_surface )  THEN
3151                IF ( TRIM( var ) == 'shf*'   )  unit = 'W/m2'
3152                IF ( TRIM( var ) == 'qsws*'  )  unit = 'W/m2'
3153             ENDIF
3154
3155          CASE DEFAULT
3156
3157             CALL tcm_check_data_output( var, unit )
3158
3159!
3160!--          Check for other modules
3161             IF ( unit == 'illegal'  .AND.  bulk_cloud_model  )  THEN
3162                CALL bcm_check_data_output( var, unit )
3163             ENDIF
3164
3165             IF ( unit == 'illegal'  .AND.  air_chemistry                      &
3166                  .AND.  (var(1:3) == 'kc_' .OR. var(1:3) == 'em_') )  THEN
3167                CALL chem_check_data_output( var, unit, i, ilen, k )
3168             ENDIF
3169
3170             IF ( unit == 'illegal' )  THEN
3171                CALL lsm_check_data_output ( var, unit, i, ilen, k )
3172             ENDIF
3173
3174             IF ( unit == 'illegal' )  THEN
3175                CALL biometeorology_check_data_output( var, unit, i, ilen, k )
3176             ENDIF
3177
3178             IF ( unit == 'illegal'  .AND.  gust_module_enabled  )  THEN
3179                CALL gust_check_data_output ( var, unit )
3180             ENDIF
3181
3182             IF ( unit == 'illegal'  .AND.  ocean_mode )  THEN
3183                CALL ocean_check_data_output( var, unit )
3184             ENDIF
3185
3186             IF ( unit == 'illegal'  .AND.  plant_canopy                       &
3187                                     .AND.  var(1:4) == 'pcm_' )  THEN
3188                CALL pcm_check_data_output( var, unit )
3189             ENDIF
3190
3191             IF ( unit == 'illegal' )  THEN
3192                CALL radiation_check_data_output( var, unit, i, ilen, k )
3193             ENDIF
3194
3195             IF ( unit == 'illegal'  .AND.  urban_surface                      &
3196                                     .AND.  var(1:4) == 'usm_' )  THEN
3197                 CALL usm_check_data_output( var, unit )
3198             ENDIF
3199
3200             IF ( unit == 'illegal'  .AND.  uv_exposure                        &
3201                                     .AND.  var(1:5) == 'uvem_' )  THEN
3202                CALL uvem_check_data_output( var, unit, i, ilen, k )
3203             ENDIF
3204
3205!
3206!--          Finally, check for user-defined quantities
3207             IF ( unit == 'illegal' )  THEN
3208                unit = ''
3209                CALL user_check_data_output( var, unit )
3210             ENDIF
3211
3212             IF ( unit == 'illegal' )  THEN
3213                IF ( data_output_user(1) /= ' ' )  THEN
3214                   message_string = 'illegal value for data_output or ' //     &
3215                         'data_output_user = "' // TRIM( data_output(i) ) // '"'
3216                   CALL message( 'check_parameters', 'PA0114', 1, 2, 0, 6, 0 )
3217                ELSE
3218                   message_string = 'illegal value for data_output = "' //     &
3219                                    TRIM( data_output(i) ) // '"'
3220                   CALL message( 'check_parameters', 'PA0115', 1, 2, 0, 6, 0 )
3221                ENDIF
3222             ENDIF
3223
3224       END SELECT
3225!
3226!--    Set the internal steering parameters appropriately
3227       IF ( k == 0 )  THEN
3228          do3d_no(j)              = do3d_no(j) + 1
3229          do3d(j,do3d_no(j))      = data_output(i)
3230          do3d_unit(j,do3d_no(j)) = unit
3231       ELSE
3232          do2d_no(j)              = do2d_no(j) + 1
3233          do2d(j,do2d_no(j))      = data_output(i)
3234          do2d_unit(j,do2d_no(j)) = unit
3235          IF ( data_output(i)(ilen-2:ilen) == '_xy' )  THEN
3236             data_output_xy(j) = .TRUE.
3237          ENDIF
3238          IF ( data_output(i)(ilen-2:ilen) == '_xz' )  THEN
3239             data_output_xz(j) = .TRUE.
3240          ENDIF
3241          IF ( data_output(i)(ilen-2:ilen) == '_yz' )  THEN
3242             data_output_yz(j) = .TRUE.
3243          ENDIF
3244       ENDIF
3245
3246       IF ( j == 1 )  THEN
3247!
3248!--       Check, if variable is already subject to averaging
3249          found = .FALSE.
3250          DO  k = 1, doav_n
3251             IF ( TRIM( doav(k) ) == TRIM( var ) )  found = .TRUE.
3252          ENDDO
3253
3254          IF ( .NOT. found )  THEN
3255             doav_n = doav_n + 1
3256             doav(doav_n) = var
3257          ENDIF
3258       ENDIF
3259
3260       i = i + 1
3261    ENDDO
3262
3263!
3264!-- Averaged 2d or 3d output requires that an averaging interval has been set
3265    IF ( doav_n > 0  .AND.  averaging_interval == 0.0_wp )  THEN
3266       WRITE( message_string, * )  'output of averaged quantity "',            &
3267                                   TRIM( doav(1) ), '_av" requires to set a ', &
3268                                   'non-zero averaging interval'
3269       CALL message( 'check_parameters', 'PA0323', 1, 2, 0, 6, 0 )
3270    ENDIF
3271
3272!
3273!-- Check sectional planes and store them in one shared array
3274    IF ( ANY( section_xy > nz + 1 ) )  THEN
3275       WRITE( message_string, * )  'section_xy must be <= nz + 1 = ', nz + 1
3276       CALL message( 'check_parameters', 'PA0319', 1, 2, 0, 6, 0 )
3277    ENDIF
3278    IF ( ANY( section_xz > ny + 1 ) )  THEN
3279       WRITE( message_string, * )  'section_xz must be <= ny + 1 = ', ny + 1
3280       CALL message( 'check_parameters', 'PA0320', 1, 2, 0, 6, 0 )
3281    ENDIF
3282    IF ( ANY( section_yz > nx + 1 ) )  THEN
3283       WRITE( message_string, * )  'section_yz must be <= nx + 1 = ', nx + 1
3284       CALL message( 'check_parameters', 'PA0321', 1, 2, 0, 6, 0 )
3285    ENDIF
3286    section(:,1) = section_xy
3287    section(:,2) = section_xz
3288    section(:,3) = section_yz
3289
3290!
3291!-- Upper plot limit for 3D arrays
3292    IF ( nz_do3d == -9999 )  nz_do3d = nzt + 1
3293
3294!
3295!-- Set output format string (used in header)
3296    SELECT CASE ( netcdf_data_format )
3297       CASE ( 1 )
3298          netcdf_data_format_string = 'netCDF classic'
3299       CASE ( 2 )
3300          netcdf_data_format_string = 'netCDF 64bit offset'
3301       CASE ( 3 )
3302          netcdf_data_format_string = 'netCDF4/HDF5'
3303       CASE ( 4 )
3304          netcdf_data_format_string = 'netCDF4/HDF5 classic'
3305       CASE ( 5 )
3306          netcdf_data_format_string = 'parallel netCDF4/HDF5'
3307       CASE ( 6 )
3308          netcdf_data_format_string = 'parallel netCDF4/HDF5 classic'
3309
3310    END SELECT
3311
3312!
3313!-- Check mask conditions
3314    DO mid = 1, max_masks
3315       IF ( data_output_masks(mid,1) /= ' '  .OR.                              &
3316            data_output_masks_user(mid,1) /= ' ' ) THEN
3317          masks = masks + 1
3318       ENDIF
3319    ENDDO
3320
3321    IF ( masks < 0  .OR.  masks > max_masks )  THEN
3322       WRITE( message_string, * )  'illegal value: masks must be >= 0 and ',   &
3323            '<= ', max_masks, ' (=max_masks)'
3324       CALL message( 'check_parameters', 'PA0325', 1, 2, 0, 6, 0 )
3325    ENDIF
3326    IF ( masks > 0 )  THEN
3327       mask_scale(1) = mask_scale_x
3328       mask_scale(2) = mask_scale_y
3329       mask_scale(3) = mask_scale_z
3330       IF ( ANY( mask_scale <= 0.0_wp ) )  THEN
3331          WRITE( message_string, * )                                           &
3332               'illegal value: mask_scale_x, mask_scale_y and mask_scale_z',   &
3333               'must be > 0.0'
3334          CALL message( 'check_parameters', 'PA0326', 1, 2, 0, 6, 0 )
3335       ENDIF
3336!
3337!--    Generate masks for masked data output
3338!--    Parallel netcdf output is not tested so far for masked data, hence
3339!--    netcdf_data_format is switched back to non-parallel output.
3340       netcdf_data_format_save = netcdf_data_format
3341       IF ( netcdf_data_format > 4 )  THEN
3342          IF ( netcdf_data_format == 5 ) netcdf_data_format = 3
3343          IF ( netcdf_data_format == 6 ) netcdf_data_format = 4
3344          message_string = 'netCDF file formats '//                            &
3345                           '5 (parallel netCDF 4) and ' //                     &
3346                           '6 (parallel netCDF 4 Classic model) '//            &
3347                           '& are currently not supported (not yet tested) ' //&
3348                           'for masked data. &Using respective non-parallel' //&
3349                           ' output for masked data.'
3350          CALL message( 'check_parameters', 'PA0383', 0, 0, 0, 6, 0 )
3351       ENDIF
3352       CALL init_masks
3353       netcdf_data_format = netcdf_data_format_save
3354    ENDIF
3355
3356!
3357!-- Check the NetCDF data format
3358    IF ( netcdf_data_format > 2 )  THEN
3359#if defined( __netcdf4 )
3360       CONTINUE
3361#else
3362       message_string = 'netCDF: netCDF4 format requested but no ' //          &
3363                        'cpp-directive __netcdf4 given & switch '  //          &
3364                        'back to 64-bit offset format'
3365       CALL message( 'check_parameters', 'PA0171', 0, 1, 0, 6, 0 )
3366       netcdf_data_format = 2
3367#endif
3368    ENDIF
3369    IF ( netcdf_data_format > 4 )  THEN
3370#if defined( __netcdf4 ) && defined( __netcdf4_parallel )
3371       CONTINUE
3372#else
3373       message_string = 'netCDF: netCDF4 parallel output requested but no ' // &
3374                        'cpp-directive __netcdf4_parallel given & switch '   //&
3375                        'back to netCDF4 non-parallel output'
3376       CALL message( 'check_parameters', 'PA0099', 0, 1, 0, 6, 0 )
3377       netcdf_data_format = netcdf_data_format - 2
3378#endif
3379    ENDIF
3380
3381!
3382!-- Calculate fixed number of output time levels for parallel netcdf output.
3383!-- The time dimension has to be defined as limited for parallel output,
3384!-- because otherwise the I/O performance drops significantly.
3385    IF ( netcdf_data_format > 4 )  THEN
3386
3387!
3388!--    Check if any of the follwoing data output interval is 0.0s, which is
3389!--    not allowed for parallel output.
3390       CALL check_dt_do( dt_do3d,           'dt_do3d'           )
3391       CALL check_dt_do( dt_do2d_xy,        'dt_do2d_xy'        )
3392       CALL check_dt_do( dt_do2d_xz,        'dt_do2d_xz'        )
3393       CALL check_dt_do( dt_do2d_yz,        'dt_do2d_yz'        )
3394       CALL check_dt_do( dt_data_output_av, 'dt_data_output_av' )
3395
3396!--    Set needed time levels (ntdim) to
3397!--    saved time levels + to be saved time levels.
3398       ntdim_3d(0) = do3d_time_count(0) + CEILING(                             &
3399                     ( end_time - MAX( skip_time_do3d,                         &
3400                                       simulated_time_at_begin )               &
3401                     ) / dt_do3d )
3402       IF ( do3d_at_begin ) ntdim_3d(0) = ntdim_3d(0) + 1
3403
3404       ntdim_3d(1) = do3d_time_count(1) + CEILING(                             &
3405                     ( end_time - MAX( skip_time_data_output_av,               &
3406                                       simulated_time_at_begin )               &
3407                     ) / dt_data_output_av )
3408
3409       ntdim_2d_xy(0) = do2d_xy_time_count(0) + CEILING(                       &
3410                        ( end_time - MAX( skip_time_do2d_xy,                   &
3411                                          simulated_time_at_begin )            &
3412                        ) / dt_do2d_xy )
3413
3414       ntdim_2d_xz(0) = do2d_xz_time_count(0) + CEILING(                       &
3415                        ( end_time - MAX( skip_time_do2d_xz,                   &
3416                                          simulated_time_at_begin )            &
3417                        ) / dt_do2d_xz )
3418
3419       ntdim_2d_yz(0) = do2d_yz_time_count(0) + CEILING(                       &
3420                        ( end_time - MAX( skip_time_do2d_yz,                   &
3421                                          simulated_time_at_begin )            &
3422                        ) / dt_do2d_yz )
3423
3424       IF ( do2d_at_begin )  THEN
3425          ntdim_2d_xy(0) = ntdim_2d_xy(0) + 1
3426          ntdim_2d_xz(0) = ntdim_2d_xz(0) + 1
3427          ntdim_2d_yz(0) = ntdim_2d_yz(0) + 1
3428       ENDIF
3429
3430       ntdim_2d_xy(1) = do2d_xy_time_count(1) + CEILING(                       &
3431                        ( end_time - MAX( skip_time_data_output_av,            &
3432                                          simulated_time_at_begin )            &
3433                        ) / dt_data_output_av )
3434
3435       ntdim_2d_xz(1) = do2d_xz_time_count(1) + CEILING(                       &
3436                        ( end_time - MAX( skip_time_data_output_av,            &
3437                                          simulated_time_at_begin )            &
3438                                 ) / dt_data_output_av )
3439
3440       ntdim_2d_yz(1) = do2d_yz_time_count(1) + CEILING(                       &
3441                        ( end_time - MAX( skip_time_data_output_av,            &
3442                                          simulated_time_at_begin )            &
3443                        ) / dt_data_output_av )
3444
3445    ENDIF
3446
3447!
3448!-- Check, whether a constant diffusion coefficient shall be used
3449    IF ( km_constant /= -1.0_wp )  THEN
3450       IF ( km_constant < 0.0_wp )  THEN
3451          WRITE( message_string, * )  'km_constant = ', km_constant, ' < 0.0'
3452          CALL message( 'check_parameters', 'PA0121', 1, 2, 0, 6, 0 )
3453       ELSE
3454          IF ( prandtl_number < 0.0_wp )  THEN
3455             WRITE( message_string, * )  'prandtl_number = ', prandtl_number,  &
3456                                         ' < 0.0'
3457             CALL message( 'check_parameters', 'PA0122', 1, 2, 0, 6, 0 )
3458          ENDIF
3459          constant_diffusion = .TRUE.
3460
3461          IF ( constant_flux_layer )  THEN
3462             message_string = 'constant_flux_layer is not allowed with fixed ' &
3463                              // 'value of km'
3464             CALL message( 'check_parameters', 'PA0123', 1, 2, 0, 6, 0 )
3465          ENDIF
3466       ENDIF
3467    ENDIF
3468
3469!
3470!-- In case of non-cyclic lateral boundaries and a damping layer for the
3471!-- potential temperature, check the width of the damping layer
3472    IF ( bc_lr /= 'cyclic' ) THEN
3473       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
3474            pt_damping_width > REAL( (nx+1) * dx ) )  THEN
3475          message_string = 'pt_damping_width out of range'
3476          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3477       ENDIF
3478    ENDIF
3479
3480    IF ( bc_ns /= 'cyclic' )  THEN
3481       IF ( pt_damping_width < 0.0_wp  .OR.                                    &
3482            pt_damping_width > REAL( (ny+1) * dy ) )  THEN
3483          message_string = 'pt_damping_width out of range'
3484          CALL message( 'check_parameters', 'PA0124', 1, 2, 0, 6, 0 )
3485       ENDIF
3486    ENDIF
3487
3488!
3489!-- Check value range for zeta = z/L
3490    IF ( zeta_min >= zeta_max )  THEN
3491       WRITE( message_string, * )  'zeta_min = ', zeta_min, ' must be less ',  &
3492                                   'than zeta_max = ', zeta_max
3493       CALL message( 'check_parameters', 'PA0125', 1, 2, 0, 6, 0 )
3494    ENDIF
3495
3496!
3497!-- Check random generator
3498    IF ( (random_generator /= 'system-specific'      .AND.                     &
3499          random_generator /= 'random-parallel'   )  .AND.                     &
3500          random_generator /= 'numerical-recipes' )  THEN
3501       message_string = 'unknown random generator: random_generator = "' //    &
3502                        TRIM( random_generator ) // '"'
3503       CALL message( 'check_parameters', 'PA0135', 1, 2, 0, 6, 0 )
3504    ENDIF
3505
3506!
3507!-- Determine upper and lower hight level indices for random perturbations
3508    IF ( disturbance_level_b == -9999999.9_wp )  THEN
3509       IF ( ocean_mode )  THEN
3510          disturbance_level_b     = zu((nzt*2)/3)
3511          disturbance_level_ind_b = ( nzt * 2 ) / 3
3512       ELSE
3513          disturbance_level_b     = zu(nzb+3)
3514          disturbance_level_ind_b = nzb + 3
3515       ENDIF
3516    ELSEIF ( disturbance_level_b < zu(3) )  THEN
3517       WRITE( message_string, * )  'disturbance_level_b = ',                   &
3518                           disturbance_level_b, ' must be >= ', zu(3), '(zu(3))'
3519       CALL message( 'check_parameters', 'PA0126', 1, 2, 0, 6, 0 )
3520    ELSEIF ( disturbance_level_b > zu(nzt-2) )  THEN
3521       WRITE( message_string, * )  'disturbance_level_b = ',                   &
3522                   disturbance_level_b, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3523       CALL message( 'check_parameters', 'PA0127', 1, 2, 0, 6, 0 )
3524    ELSE
3525       DO  k = 3, nzt-2
3526          IF ( disturbance_level_b <= zu(k) )  THEN
3527             disturbance_level_ind_b = k
3528             EXIT
3529          ENDIF
3530       ENDDO
3531    ENDIF
3532
3533    IF ( disturbance_level_t == -9999999.9_wp )  THEN
3534       IF ( ocean_mode )  THEN
3535          disturbance_level_t     = zu(nzt-3)
3536          disturbance_level_ind_t = nzt - 3
3537       ELSE
3538          disturbance_level_t     = zu(nzt/3)
3539          disturbance_level_ind_t = nzt / 3
3540       ENDIF
3541    ELSEIF ( disturbance_level_t > zu(nzt-2) )  THEN
3542       WRITE( message_string, * )  'disturbance_level_t = ',                   &
3543                   disturbance_level_t, ' must be <= ', zu(nzt-2), '(zu(nzt-2))'
3544       CALL message( 'check_parameters', 'PA0128', 1, 2, 0, 6, 0 )
3545    ELSEIF ( disturbance_level_t < disturbance_level_b )  THEN
3546       WRITE( message_string, * )  'disturbance_level_t = ',                   &
3547                   disturbance_level_t, ' must be >= disturbance_level_b = ',  &
3548                   disturbance_level_b
3549       CALL message( 'check_parameters', 'PA0129', 1, 2, 0, 6, 0 )
3550    ELSE
3551       DO  k = 3, nzt-2
3552          IF ( disturbance_level_t <= zu(k) )  THEN
3553             disturbance_level_ind_t = k
3554             EXIT
3555          ENDIF
3556       ENDDO
3557    ENDIF
3558
3559!
3560!-- Check again whether the levels determined this way are ok.
3561!-- Error may occur at automatic determination and too few grid points in
3562!-- z-direction.
3563    IF ( disturbance_level_ind_t < disturbance_level_ind_b )  THEN
3564       WRITE( message_string, * )  'disturbance_level_ind_t = ',               &
3565                disturbance_level_ind_t, ' must be >= ',                       &
3566                'disturbance_level_ind_b = ', disturbance_level_ind_b
3567       CALL message( 'check_parameters', 'PA0130', 1, 2, 0, 6, 0 )
3568    ENDIF
3569
3570!
3571!-- Determine the horizontal index range for random perturbations.
3572!-- In case of non-cyclic horizontal boundaries, no perturbations are imposed
3573!-- near the inflow and the perturbation area is further limited to ...(1)
3574!-- after the initial phase of the flow.
3575
3576    IF ( bc_lr /= 'cyclic' )  THEN
3577       IF ( inflow_disturbance_begin == -1 )  THEN
3578          inflow_disturbance_begin = MIN( 10, nx/2 )
3579       ENDIF
3580       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > nx )&
3581       THEN
3582          message_string = 'inflow_disturbance_begin out of range'
3583          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
3584       ENDIF
3585       IF ( inflow_disturbance_end == -1 )  THEN
3586          inflow_disturbance_end = MIN( 100, 3*nx/4 )
3587       ENDIF
3588       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > nx )    &
3589       THEN
3590          message_string = 'inflow_disturbance_end out of range'
3591          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
3592       ENDIF
3593    ELSEIF ( bc_ns /= 'cyclic' )  THEN
3594       IF ( inflow_disturbance_begin == -1 )  THEN
3595          inflow_disturbance_begin = MIN( 10, ny/2 )
3596       ENDIF
3597       IF ( inflow_disturbance_begin < 0  .OR.  inflow_disturbance_begin > ny )&
3598       THEN
3599          message_string = 'inflow_disturbance_begin out of range'
3600          CALL message( 'check_parameters', 'PA0131', 1, 2, 0, 6, 0 )
3601       ENDIF
3602       IF ( inflow_disturbance_end == -1 )  THEN
3603          inflow_disturbance_end = MIN( 100, 3*ny/4 )
3604       ENDIF
3605       IF ( inflow_disturbance_end < 0  .OR.  inflow_disturbance_end > ny )    &
3606       THEN
3607          message_string = 'inflow_disturbance_end out of range'
3608          CALL message( 'check_parameters', 'PA0132', 1, 2, 0, 6, 0 )
3609       ENDIF
3610    ENDIF
3611
3612    IF ( random_generator == 'random-parallel' )  THEN
3613       dist_nxl = nxl;  dist_nxr = nxr
3614       dist_nys = nys;  dist_nyn = nyn
3615       IF ( bc_lr == 'radiation/dirichlet' )  THEN
3616          dist_nxr    = MIN( nx - inflow_disturbance_begin, nxr )
3617          dist_nxl(1) = MAX( nx - inflow_disturbance_end, nxl )
3618       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
3619          dist_nxl    = MAX( inflow_disturbance_begin, nxl )
3620          dist_nxr(1) = MIN( inflow_disturbance_end, nxr )
3621       ELSEIF ( bc_lr == 'nested'  .OR.  bc_lr == 'nesting_offline' )  THEN
3622          dist_nxl    = MAX( inflow_disturbance_begin, nxl )
3623          dist_nxr    = MIN( nx - inflow_disturbance_begin, nxr )
3624       ENDIF
3625       IF ( bc_ns == 'dirichlet/radiation' )  THEN
3626          dist_nyn    = MIN( ny - inflow_disturbance_begin, nyn )
3627          dist_nys(1) = MAX( ny - inflow_disturbance_end, nys )
3628       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
3629          dist_nys    = MAX( inflow_disturbance_begin, nys )
3630          dist_nyn(1) = MIN( inflow_disturbance_end, nyn )
3631       ELSEIF ( bc_ns == 'nested'  .OR.  bc_ns == 'nesting_offline' )  THEN
3632          dist_nys    = MAX( inflow_disturbance_begin, nys )
3633          dist_nyn    = MIN( ny - inflow_disturbance_begin, nyn )
3634       ENDIF
3635    ELSE
3636       dist_nxl = 0;  dist_nxr = nx
3637       dist_nys = 0;  dist_nyn = ny
3638       IF ( bc_lr == 'radiation/dirichlet' )  THEN
3639          dist_nxr    = nx - inflow_disturbance_begin
3640          dist_nxl(1) = nx - inflow_disturbance_end
3641       ELSEIF ( bc_lr == 'dirichlet/radiation' )  THEN
3642          dist_nxl    = inflow_disturbance_begin
3643          dist_nxr(1) = inflow_disturbance_end
3644       ELSEIF ( bc_lr == 'nested'  .OR.  bc_lr == 'nesting_offline' )  THEN
3645          dist_nxr    = nx - inflow_disturbance_begin
3646          dist_nxl    = inflow_disturbance_begin
3647       ENDIF
3648       IF ( bc_ns == 'dirichlet/radiation' )  THEN
3649          dist_nyn    = ny - inflow_disturbance_begin
3650          dist_nys(1) = ny - inflow_disturbance_end
3651       ELSEIF ( bc_ns == 'radiation/dirichlet' )  THEN
3652          dist_nys    = inflow_disturbance_begin
3653          dist_nyn(1) = inflow_disturbance_end
3654       ELSEIF ( bc_ns == 'nested'  .OR.  bc_ns == 'nesting_offline' )  THEN
3655          dist_nyn    = ny - inflow_disturbance_begin
3656          dist_nys    = inflow_disturbance_begin
3657       ENDIF
3658    ENDIF
3659
3660!
3661!-- A turbulent inflow requires Dirichlet conditions at the respective inflow
3662!-- boundary (so far, a turbulent inflow is realized from the left side only)
3663    IF ( turbulent_inflow  .AND.  bc_lr /= 'dirichlet/radiation' )  THEN
3664       message_string = 'turbulent_inflow = .T. requires a Dirichlet ' //      &
3665                        'condition at the inflow boundary'
3666       CALL message( 'check_parameters', 'PA0133', 1, 2, 0, 6, 0 )
3667    ENDIF
3668
3669!
3670!-- Turbulent inflow requires that 3d arrays have been cyclically filled with
3671!-- data from prerun in the first main run
3672    IF ( turbulent_inflow  .AND.  initializing_actions /= 'cyclic_fill'  .AND. &
3673         initializing_actions /= 'read_restart_data' )  THEN
3674       message_string = 'turbulent_inflow = .T. requires ' //                  &
3675                        'initializing_actions = ''cyclic_fill'' or ' //        &
3676                        'initializing_actions = ''read_restart_data'' '
3677       CALL message( 'check_parameters', 'PA0055', 1, 2, 0, 6, 0 )
3678    ENDIF
3679
3680!
3681!-- In case of turbulent inflow calculate the index of the recycling plane
3682    IF ( turbulent_inflow )  THEN
3683       IF ( recycling_width <= dx  .OR.  recycling_width >= nx * dx )  THEN
3684          WRITE( message_string, * )  'illegal value for recycling_width: ',   &
3685                                      recycling_width
3686          CALL message( 'check_parameters', 'PA0134', 1, 2, 0, 6, 0 )
3687       ENDIF
3688!
3689!--    Calculate the index
3690       recycling_plane = recycling_width / dx
3691!
3692!--    Because the y-shift is done with a distance of INT( npey / 2 ) no shift
3693!--    is possible if there is only one PE in y direction.
3694       IF ( recycling_yshift .AND. pdims(2) < 2 )  THEN
3695          WRITE( message_string, * )  'recycling_yshift = .T. requires more',  &
3696                                      ' than one processor in y direction'
3697          CALL message( 'check_parameters', 'PA0421', 1, 2, 0, 6, 0 )
3698       ENDIF
3699    ENDIF
3700
3701
3702    IF ( turbulent_outflow )  THEN
3703!
3704!--    Turbulent outflow requires Dirichlet conditions at the respective inflow
3705!--    boundary (so far, a turbulent outflow is realized at the right side only)
3706       IF ( bc_lr /= 'dirichlet/radiation' )  THEN
3707          message_string = 'turbulent_outflow = .T. requires ' //              &
3708                           'bc_lr = "dirichlet/radiation"'
3709          CALL message( 'check_parameters', 'PA0038', 1, 2, 0, 6, 0 )
3710       ENDIF
3711!
3712!--    The ouflow-source plane must lay inside the model domain
3713       IF ( outflow_source_plane < dx  .OR.  &
3714            outflow_source_plane > nx * dx )  THEN
3715          WRITE( message_string, * )  'illegal value for outflow_source'//     &
3716                                      '_plane: ', outflow_source_plane
3717          CALL message( 'check_parameters', 'PA0145', 1, 2, 0, 6, 0 )
3718       ENDIF
3719    ENDIF
3720
3721!
3722!-- Determine damping level index for 1D model
3723    IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
3724       IF ( damp_level_1d == -1.0_wp )  THEN
3725          damp_level_1d     = zu(nzt+1)
3726          damp_level_ind_1d = nzt + 1
3727       ELSEIF ( damp_level_1d < 0.0_wp  .OR.  damp_level_1d > zu(nzt+1) )  THEN
3728          WRITE( message_string, * )  'damp_level_1d = ', damp_level_1d,       &
3729                 ' must be >= 0.0 and <= ', zu(nzt+1), '(zu(nzt+1))'
3730          CALL message( 'check_parameters', 'PA0136', 1, 2, 0, 6, 0 )
3731       ELSE
3732          DO  k = 1, nzt+1
3733             IF ( damp_level_1d <= zu(k) )  THEN
3734                damp_level_ind_1d = k
3735                EXIT
3736             ENDIF
3737          ENDDO
3738       ENDIF
3739    ENDIF
3740
3741!
3742!-- Check some other 1d-model parameters
3743    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                   &
3744         TRIM( mixing_length_1d ) /= 'blackadar' )  THEN
3745       message_string = 'mixing_length_1d = "' // TRIM( mixing_length_1d ) //  &
3746                        '" is unknown'
3747       CALL message( 'check_parameters', 'PA0137', 1, 2, 0, 6, 0 )
3748    ENDIF
3749    IF ( TRIM( dissipation_1d ) /= 'as_in_3d_model'  .AND.                     &
3750         TRIM( dissipation_1d ) /= 'detering'  .AND.                           &
3751         TRIM( dissipation_1d ) /= 'prognostic' )  THEN
3752       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) //      &
3753                        '" is unknown'
3754       CALL message( 'check_parameters', 'PA0138', 1, 2, 0, 6, 0 )
3755    ENDIF
3756    IF ( TRIM( mixing_length_1d ) /= 'as_in_3d_model'  .AND.                   &
3757         TRIM( dissipation_1d ) == 'as_in_3d_model' )  THEN
3758       message_string = 'dissipation_1d = "' // TRIM( dissipation_1d ) //      &
3759                        '" requires mixing_length_1d = "as_in_3d_model"'
3760       CALL message( 'check_parameters', 'PA0485', 1, 2, 0, 6, 0 )
3761    ENDIF
3762
3763!
3764!-- Set time for the next user defined restart (time_restart is the
3765!-- internal parameter for steering restart events)
3766    IF ( restart_time /= 9999999.9_wp )  THEN
3767       IF ( restart_time > time_since_reference_point )  THEN
3768          time_restart = restart_time
3769       ENDIF
3770    ELSE
3771!
3772!--    In case of a restart run, set internal parameter to default (no restart)
3773!--    if the NAMELIST-parameter restart_time is omitted
3774       time_restart = 9999999.9_wp
3775    ENDIF
3776
3777!
3778!-- Check pressure gradient conditions
3779    IF ( dp_external  .AND.  conserve_volume_flow )  THEN
3780       WRITE( message_string, * )  'Both dp_external and conserve_volume_flo', &
3781            'w are .TRUE. but one of them must be .FALSE.'
3782       CALL message( 'check_parameters', 'PA0150', 1, 2, 0, 6, 0 )
3783    ENDIF
3784    IF ( dp_external )  THEN
3785       IF ( dp_level_b < zu(nzb)  .OR.  dp_level_b > zu(nzt) )  THEN
3786          WRITE( message_string, * )  'dp_level_b = ', dp_level_b, ' is out ', &
3787               ' of range [zu(nzb), zu(nzt)]'
3788          CALL message( 'check_parameters', 'PA0151', 1, 2, 0, 6, 0 )
3789       ENDIF
3790       IF ( .NOT. ANY( dpdxy /= 0.0_wp ) )  THEN
3791          WRITE( message_string, * )  'dp_external is .TRUE. but dpdxy is ze', &
3792               'ro, i.e. the external pressure gradient will not be applied'
3793          CALL message( 'check_parameters', 'PA0152', 0, 1, 0, 6, 0 )
3794       ENDIF
3795    ENDIF
3796    IF ( ANY( dpdxy /= 0.0_wp )  .AND.  .NOT.  dp_external )  THEN
3797       WRITE( message_string, * )  'dpdxy is nonzero but dp_external is ',     &
3798            '.FALSE., i.e. the external pressure gradient & will not be applied'
3799       CALL message( 'check_parameters', 'PA0153', 0, 1, 0, 6, 0 )
3800    ENDIF
3801    IF ( conserve_volume_flow )  THEN
3802       IF ( TRIM( conserve_volume_flow_mode ) == 'default' )  THEN
3803
3804          conserve_volume_flow_mode = 'initial_profiles'
3805
3806       ELSEIF ( TRIM( conserve_volume_flow_mode ) /= 'initial_profiles' .AND.  &
3807            TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' )  THEN
3808          WRITE( message_string, * )  'unknown conserve_volume_flow_mode: ',   &
3809               conserve_volume_flow_mode
3810          CALL message( 'check_parameters', 'PA0154', 1, 2, 0, 6, 0 )
3811       ENDIF
3812       IF ( (bc_lr /= 'cyclic'  .OR.  bc_ns /= 'cyclic')  .AND.                &
3813          TRIM( conserve_volume_flow_mode ) == 'bulk_velocity' )  THEN
3814          WRITE( message_string, * )  'non-cyclic boundary conditions ',       &
3815               'require  conserve_volume_flow_mode = ''initial_profiles'''
3816          CALL message( 'check_parameters', 'PA0155', 1, 2, 0, 6, 0 )
3817       ENDIF
3818    ENDIF
3819    IF ( ( u_bulk /= 0.0_wp  .OR.  v_bulk /= 0.0_wp )  .AND.                   &
3820         ( .NOT. conserve_volume_flow  .OR.                                    &
3821         TRIM( conserve_volume_flow_mode ) /= 'bulk_velocity' ) )  THEN
3822       WRITE( message_string, * )  'nonzero bulk velocity requires ',          &
3823            'conserve_volume_flow = .T. and ',                                 &
3824            'conserve_volume_flow_mode = ''bulk_velocity'''
3825       CALL message( 'check_parameters', 'PA0157', 1, 2, 0, 6, 0 )
3826    ENDIF
3827
3828!
3829!-- Check particle attributes
3830    IF ( particle_color /= 'none' )  THEN
3831       IF ( particle_color /= 'absuv'  .AND.  particle_color /= 'pt*'  .AND.   &
3832            particle_color /= 'z' )  THEN
3833          message_string = 'illegal value for parameter particle_color: ' //   &
3834                           TRIM( particle_color)
3835          CALL message( 'check_parameters', 'PA0313', 1, 2, 0, 6, 0 )
3836       ELSE
3837          IF ( color_interval(2) <= color_interval(1) )  THEN
3838             message_string = 'color_interval(2) <= color_interval(1)'
3839             CALL message( 'check_parameters', 'PA0315', 1, 2, 0, 6, 0 )
3840          ENDIF
3841       ENDIF
3842    ENDIF
3843
3844    IF ( particle_dvrpsize /= 'none' )  THEN
3845       IF ( particle_dvrpsize /= 'absw' )  THEN
3846          message_string = 'illegal value for parameter particle_dvrpsize:' // &
3847                           ' ' // TRIM( particle_dvrpsize)
3848          CALL message( 'check_parameters', 'PA0314', 1, 2, 0, 6, 0 )
3849       ELSE
3850          IF ( dvrpsize_interval(2) <= dvrpsize_interval(1) )  THEN
3851             message_string = 'dvrpsize_interval(2) <= dvrpsize_interval(1)'
3852             CALL message( 'check_parameters', 'PA0316', 1, 2, 0, 6, 0 )
3853          ENDIF
3854       ENDIF
3855    ENDIF
3856
3857!
3858!-- Prevent empty time records in volume, cross-section and masked data in case
3859!-- of non-parallel netcdf-output in restart runs
3860    IF ( netcdf_data_format < 5 )  THEN
3861       IF ( TRIM( initializing_actions ) == 'read_restart_data' )  THEN
3862          do3d_time_count    = 0
3863          do2d_xy_time_count = 0
3864          do2d_xz_time_count = 0
3865          do2d_yz_time_count = 0
3866          domask_time_count  = 0
3867       ENDIF
3868    ENDIF
3869
3870!
3871!-- Check for valid setting of most_method
3872    IF ( TRIM( most_method ) /= 'circular'  .AND.                              &
3873         TRIM( most_method ) /= 'newton'    .AND.                              &
3874         TRIM( most_method ) /= 'lookup' )  THEN
3875       message_string = 'most_method = "' // TRIM( most_method ) //            &
3876                        '" is unknown'
3877       CALL message( 'check_parameters', 'PA0416', 1, 2, 0, 6, 0 )
3878    ENDIF
3879
3880!
3881!-- Check roughness length, which has to be smaller than dz/2
3882    IF ( ( constant_flux_layer .OR.  &
3883           INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )       &
3884         .AND. roughness_length >= 0.5 * dz(1) )  THEN
3885       message_string = 'roughness_length must be smaller than dz/2'
3886       CALL message( 'check_parameters', 'PA0424', 1, 2, 0, 6, 0 )
3887    ENDIF
3888
3889!
3890!-- Vertical nesting: check fine and coarse grid compatibility for data exchange
3891    IF ( vnested )  CALL vnest_check_parameters
3892
3893!
3894!-- Check if topography is read from file in case of complex terrain simulations
3895    IF ( complex_terrain  .AND.  TRIM( topography ) /= 'read_from_file' )  THEN
3896       message_string = 'complex_terrain requires topography' //               &
3897                        ' = ''read_from_file'''
3898       CALL message( 'check_parameters', 'PA0295', 1, 2, 0, 6, 0 )
3899    ENDIF
3900
3901!
3902!-- Check if vertical grid stretching is switched off in case of complex
3903!-- terrain simulations
3904    IF ( complex_terrain  .AND.                                                &
3905         ANY( dz_stretch_level_start /= -9999999.9_wp ) )  THEN
3906       message_string = 'Vertical grid stretching is not allowed for ' //      &
3907                        'complex_terrain = .T.'
3908       CALL message( 'check_parameters', 'PA0473', 1, 2, 0, 6, 0 )
3909    ENDIF
3910
3911    CALL location_message( 'finished', .TRUE. )
3912!
3913!-- Check &userpar parameters
3914    CALL user_check_parameters
3915
3916 CONTAINS
3917
3918!------------------------------------------------------------------------------!
3919! Description:
3920! ------------
3921!> Check the length of data output intervals. In case of parallel NetCDF output
3922!> the time levels of the output files need to be fixed. Therefore setting the
3923!> output interval to 0.0s (usually used to output each timestep) is not
3924!> possible as long as a non-fixed timestep is used.
3925!------------------------------------------------------------------------------!
3926
3927    SUBROUTINE check_dt_do( dt_do, dt_do_name )
3928
3929       IMPLICIT NONE
3930
3931       CHARACTER (LEN=*), INTENT (IN) :: dt_do_name !< parin variable name
3932
3933       REAL(wp), INTENT (INOUT)       :: dt_do      !< data output interval
3934
3935       IF ( dt_do == 0.0_wp )  THEN
3936          IF ( dt_fixed )  THEN
3937             WRITE( message_string, '(A,F9.4,A)' )  'Output at every '  //     &
3938                    'timestep is wanted (' // dt_do_name // ' = 0.0).&'//      &
3939                    'The output interval is set to the fixed timestep dt '//   &
3940                    '= ', dt, 's.'
3941             CALL message( 'check_parameters', 'PA0060', 0, 0, 0, 6, 0 )
3942             dt_do = dt
3943          ELSE
3944             message_string = dt_do_name // ' = 0.0 while using a ' //         &
3945                              'variable timestep and parallel netCDF4 ' //     &
3946                              'is not allowed.'
3947             CALL message( 'check_parameters', 'PA0081', 1, 2, 0, 6, 0 )
3948          ENDIF
3949       ENDIF
3950
3951    END SUBROUTINE check_dt_do
3952
3953
3954
3955!------------------------------------------------------------------------------!
3956! Description:
3957! ------------
3958!> Set the bottom and top boundary conditions for humidity and scalars.
3959!------------------------------------------------------------------------------!
3960
3961    SUBROUTINE set_bc_scalars( sq, bc_b, bc_t, ibc_b, ibc_t, err_nr_b, err_nr_t )
3962
3963
3964       IMPLICIT NONE
3965
3966       CHARACTER (LEN=1)   ::  sq         !< name of scalar quantity
3967       CHARACTER (LEN=*)   ::  bc_b       !< bottom boundary condition
3968       CHARACTER (LEN=*)   ::  bc_t       !< top boundary condition
3969       CHARACTER (LEN=*)   ::  err_nr_b   !< error number if bottom bc is unknown
3970       CHARACTER (LEN=*)   ::  err_nr_t   !< error number if top bc is unknown
3971
3972       INTEGER(iwp)        ::  ibc_b      !< index for bottom boundary condition
3973       INTEGER(iwp)        ::  ibc_t      !< index for top boundary condition
3974
3975!
3976!--    Set Integer flags and check for possilbe errorneous settings for bottom
3977!--    boundary condition
3978       IF ( bc_b == 'dirichlet' )  THEN
3979          ibc_b = 0
3980       ELSEIF ( bc_b == 'neumann' )  THEN
3981          ibc_b = 1
3982       ELSE
3983          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
3984                           '_b ="' // TRIM( bc_b ) // '"'
3985          CALL message( 'check_parameters', err_nr_b, 1, 2, 0, 6, 0 )
3986       ENDIF
3987!
3988!--    Set Integer flags and check for possilbe errorneous settings for top
3989!--    boundary condition
3990       IF ( bc_t == 'dirichlet' )  THEN
3991          ibc_t = 0
3992       ELSEIF ( bc_t == 'neumann' )  THEN
3993          ibc_t = 1
3994       ELSEIF ( bc_t == 'initial_gradient' )  THEN
3995          ibc_t = 2
3996       ELSEIF ( bc_t == 'nested'  .OR.  bc_t == 'nesting_offline' )  THEN
3997          ibc_t = 3
3998       ELSE
3999          message_string = 'unknown boundary condition: bc_' // TRIM( sq ) //  &
4000                           '_t ="' // TRIM( bc_t ) // '"'
4001          CALL message( 'check_parameters', err_nr_t, 1, 2, 0, 6, 0 )
4002       ENDIF
4003
4004
4005    END SUBROUTINE set_bc_scalars
4006
4007
4008
4009!------------------------------------------------------------------------------!
4010! Description:
4011! ------------
4012!> Check for consistent settings of bottom boundary conditions for humidity
4013!> and scalars.
4014!------------------------------------------------------------------------------!
4015
4016    SUBROUTINE check_bc_scalars( sq, bc_b, ibc_b,                      &
4017                                 err_nr_1, err_nr_2,                   &
4018                                 constant_flux, surface_initial_change )
4019
4020
4021       IMPLICIT NONE
4022
4023       CHARACTER (LEN=1)   ::  sq                       !< name of scalar quantity
4024       CHARACTER (LEN=*)   ::  bc_b                     !< bottom boundary condition
4025       CHARACTER (LEN=*)   ::  err_nr_1                 !< error number of first error
4026       CHARACTER (LEN=*)   ::  err_nr_2                 !< error number of second error
4027
4028       INTEGER(iwp)        ::  ibc_b                    !< index of bottom boundary condition
4029
4030       LOGICAL             ::  constant_flux            !< flag for constant-flux layer
4031
4032       REAL(wp)            ::  surface_initial_change   !< value of initial change at the surface
4033
4034!
4035!--    A given surface value implies Dirichlet boundary condition for
4036!--    the respective quantity. In this case specification of a constant flux is
4037!--    forbidden. However, an exception is made for large-scale forcing as well
4038!--    as land-surface model.
4039       IF ( .NOT. land_surface  .AND.  .NOT. large_scale_forcing )  THEN
4040          IF ( ibc_b == 0  .AND.  constant_flux )  THEN
4041             message_string = 'boundary condition: bc_' // TRIM( sq ) //       &
4042                              '_b ' // '= "' // TRIM( bc_b ) //                &
4043                              '" is not allowed with prescribed surface flux'
4044             CALL message( 'check_parameters', err_nr_1, 1, 2, 0, 6, 0 )
4045          ENDIF
4046       ENDIF
4047       IF ( constant_flux  .AND.  surface_initial_change /= 0.0_wp )  THEN
4048          WRITE( message_string, * )  'a prescribed surface flux is not allo', &
4049                 'wed with ', sq, '_surface_initial_change (/=0) = ',          &
4050                 surface_initial_change
4051          CALL message( 'check_parameters', err_nr_2, 1, 2, 0, 6, 0 )
4052       ENDIF
4053
4054
4055    END SUBROUTINE check_bc_scalars
4056
4057
4058
4059 END SUBROUTINE check_parameters
Note: See TracBrowser for help on using the repository browser.