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

Last change on this file since 3529 was 3529, checked in by gronemeier, 6 years ago

change date format in output files; add global attributes; change fill_value; move definition of UTM and lon/lat into subroutine; change attributes of time variable; read optional attributes from input netcdf file; update test cases

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