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

Last change on this file since 3994 was 3994, checked in by suehring, 5 years ago

new module for diagnostic output quantities added + output of turbulence intensity

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