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

Last change on this file since 3302 was 3302, checked in by raasch, 4 years ago

Craik-Leibovich force and wave breaking effects added to ocean mode, header output for ocean mode

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