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

Last change on this file since 3545 was 3545, checked in by gronemeier, 3 years ago

removed debug output; removed rans_mode from namelist; altered order of check_parameter-calls of modules

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