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

Last change on this file since 3761 was 3761, checked in by raasch, 3 years ago

unused variables removed, OpenACC directives re-formatted, statements added to avoid compiler warnings

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