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

Last change on this file since 4178 was 4176, checked in by oliver.maas, 2 years ago

added recycle_absolute_quantities to initialization_parameters namelist in parin.f90, bugfix: replace PA184 by PA0184 in check_parameters.f90

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