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

Last change on this file since 3962 was 3933, checked in by kanani, 5 years ago

Bugfixes and clean-up for output quantity theta_2m*

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