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

Last change on this file since 4017 was 4017, checked in by schwenkel, 5 years ago

Modularization of all lagrangian particle model code components

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