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

Last change on this file since 3676 was 3668, checked in by maronga, 6 years ago

removed most_methods circular and lookup. added improved version of palm_csd

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